Strong biases in retrieved atmospheric composition caused by strong day-night chemical heterogeneities

Most planets currently amenable to transit spectroscopy are close enough to their host star to exhibit a relatively strong day to night temperature gradient. For hot planets, this leads to cause a chemical composition dichotomy between the two hemispheres. In the extreme case of ultra hot jupiters, some species, such as molecular hydrogen and water, are strongly dissociated on the day-side while others, such as carbon monoxide, are not. However, most current retrieval algorithm rely on 1D forward models that are unable to model this effect. We thus investigate how the 3D structure of the atmosphere biases the abundances retrieved using commonly used algorithms. We study the case of Wasp-121b as a prototypical ultra hot Jupiter. We use the simulations of this planet performed with the SPARC/MIT global climate model (GCM) and generate transmission spectra that fully account for the 3D structure of the atmosphere with Pytmopsh3R. These spectra are then analyzed using the \taurex retrieval code. We find that such ultra hot jupiter's transmission spectra exhibit muted H$_2$O features that originate in the night-side where the temperature, hence the scale-height, is smaller than on the day-side. However, the spectral features of molecules present on the day-side are boosted by both its high temperature and low mean molecular weight. As a result, the retrieved parameters are strongly biased compared to the ground truth. In particular the [CO]/[H$_2$O] is overestimated by one to three orders of magnitude. This must be kept in mind when using such retrieval analysis to infer the C/O ratio of a planet's atmosphere. We also discuss whether indicators can allow us to infer the 3D structure of an observed atmosphere. Finally we show that HST/WFC3 transmission data of Wasp-121b are compatible with the day-night thermal and compositional dichotomy predicted by models.


Introduction
Since the discovery of the first exoplanets, observations have shown a great diversity of objects, from Earth-like planets to Ultra Hot Jupiters (UHJ). Orbiting very close to their star, UHJs can reach an atmospheric temperature high enough to trigger the thermal dissociation of some of the chemical species, such as H 2 O and H 2 (Lodders & Fegley 2002;Visscher et al. 2006Visscher et al. , 2010. They thus offer the opportunity to study the chemistry and physics of planetary atmospheres under extreme conditions, for which we have no equivalent in the Solar System. These interesting objects will be prime targets for new generation space missions like the James-Webb Space Telescope (JWST) (Beichman et al. 2014) and Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL) (Tinetti et al. 2018).
Currently, only a few UHJs have been discovered and studied using both transit and eclipse spectroscopy (Wright et al. 2012;Haynes et al. 2015;Sheppard et al. 2017;Evans et al. 2017;. However, the analysis of UHJs is not simple, due to their complex chemical composition and dynamics. Bayesian retrieval procedures being computationally-intensive, it is necessary to make strong assumptions to speed-up the atmospheric forward model. But a too simplistic set of assumptions can lead to strongly biased interpretations where the retrieved abundances are much higher than the expected ones. Evans et al. (2017), for example, used transit data to suggest presence of supersolar FeH abundance to fit HST-WFC3 data of Wasp-121b and explained the reduced water shape around 1.3 µm. Parmentier et al. (2018) later suggest that this could be due to the presence of partially CaTiO 3 cloudy atmosphere.
In the present study, we investigate whether the variations in composition inside the atmosphere of UHJs may affect transmission spectroscopy as severely as emission spectroscopy. Those tidally locked planets present a strong day/night contrast both in temperature (Sudarsky et al. 2000;Arcangeli et al. 2018) and in chemical heterogeneities due to the thermal dissociation of certain species such as H 2 O and H 2 Marley et al. 2017). As the thermal dissociation of the species is strongly linked to the temperature, the thermal day-night dichotomy entails a chemical dichotomy, with a day side devoid of water above about 100 mbar and a night side where water is present everywhere. However, other species such as CO requires higher temperature to dissociate (Lodders & Fegley 2002) and are expected to remain constant in the atmosphere which would change the apparent [CO]/[H 2 O] ratio. Caldas et al. (2019) developed a fully 3D model that allows to generate a transmission spectrum which considers the 3D structure of the atmosphere. They show that the light that goes through the planetary atmosphere carries the information of a portion of atmospheres that significantly extend around the limb. They highlight that the rays of light go through the day side first and then to the night side implying that strong 3D effects on the transmission spectrum across the limb occur when those effects along the limb are negligible. Hence, 1D transmission models cannot well take into account chemical heterogeneities describe above thus the need to use 3D transmission models. mal effects. To that purpose, we focus on Wasp-121b (Evans et al. 2016(Evans et al. , 2017Parmentier et al. 2018) as a prototype for UHJs. The complexity of the UHJ atmospheres describes in our GCM model may suggest that we would need a more complex framework to analyze the transmission spectra of those particular planets. Thus, we will need fully 3D forward models to simulate realistic transmission spectra to better fit transit observations of UHJs. Hence, 1D retrieval models such as TauREx are probably biased in their analysis.
Hereafter, we first describe the observational chain used to simulate JWST observations in Sect. 2. Then, in Sect. 3, we explain the numerical experiments done to unravel the biases starting from a very simple, parametric modeling of Wasp-121b to a more elaborate GCM model. Moreover, we explain how thermal dissociation induces strong heterogeneities of composition between the day and night sides of the planet. We describe our retrieval results in Sect. 4. Finally, we highlight the main conclusion of our study and we discuss the limitations of our method and our models in Sect. 5.

SPARC/MITgcm global circulation model
We use the SPARC/MITgcm global circulation model (Showman et al. 2009) to model the atmosphere of WASP-121b. The model solves the primitive equations on a cubic-sphere grid. It has been successfully applied to a wide range of hot Jupiters (Showman et al. 2009;Kataria et al. 2015;Lewis et al. 2017;Parmentier et al. 2013;Parmentier et al. 2016) and has recently been used to study a few ultra-hot Jupiters in details Parmentier et al. 2018;Arcangeli et al. 2019). The model used here is exactly the same as described in . Our pressure ranges from 200 bar to 2 µbar over 53 levels, We use a horizontal resolution of C32, equivalent to an approximate resolution of 128 cells in longitude and 64 in latitude. The radiative transfer is handled with a two-stream radiation scheme (Marley & McKay 1999), with the opacities treated using 8 correlated-k coefficients (Goody & Yung 1989) withing each of our 11 wavelength bins (Kataria et al. 2013).
Chemical equilibrium is assumed when calculating the opacities, meaning that thermal dissociation of molecules, including water and hydrogen is naturally taken into account. However, for practical purposes, we assume that the mean molecular weight and the heat capacity are constant throughout the atmosphere. The change in mean molecular weight is explored a-posteriori when projecting the GCM thermal structure in the Pytmosph3R grid (see Sect. 2.2). This model does not consider the impact of H 2 recombination in the atmosphere which can have a nonnegligible effect on the physics and the dynamics (Tan & Komacek 2019; see Sect. 3.1 for further details). Fig 1, 2 show the temperature maps generated by the GCM which are structured in three main parts: 1. A quasi isothermal annulus around 2500 K from the surface pressure to approximately a pressure of 0.1 bar. We are here in the deep atmosphere, at high pressure, thus the redistribution of energy is very efficient due to jet stream. 2. Above this annulus, a cold night side where the temperature gradient decreases with altitude going from 1000 K to 500 K at whole latitudes. 3. Then, a very hot day side where the temperature gradient decreases with altitude going from 3500 K to 3000 K at whole latitudes. Those high temperatures enlarge the scale height implying a strong asymmetry in altitude between the day and the night side.
The transition between the day and the night side is very sharp, with a quasi isothermal terminator around 2200 K. We used in our study this global thermal structure to build idealize case of , as it will be explain in Sect 3. As introduced before, the temperatures are high enough to allow dissociation of species, especially the water. We calculate the abundance of H 2 O in the GCM simulations using the analytical equations provided by Parmentier et al. (2018) and we plot the abundances maps of H 2 O in Fig 1, 2. We can see on the equatorial and polar maps that the water abundance is linked to the temperature since there is a total absence of water deep in the atmosphere on the day side, while the water abundance reaches solar abundance in the night side. The pressure dependence of water dissociation appears in the limb map, which is quasi-isothermal, since the water abundance reaches the solar abundance in the surface pressure then decreases to less than 10 −12 at the top of the atmosphere (see Fig 5).

Generation of transmission spectra with Pytmosph3R
Transmission spectra are computed using Pytmosph3R as described in Caldas et al. (2019). Pytmosph3R is designed to simulate transmission spectra based on any 3D atmospheric structure, including outputs of a Global Circulation Model (GCM ). It produces transmittance maps at any wavelength that can later be spatially integrated to yield the transmission spectrum. The code can account for molecular opacities using the correlatedk method (Fu & Liou 1992) or monochromatic cross-sections calculated by ExoMol (Yurchenko et al. 2011;Tennyson & Yurchenko 2012;Barton et al. 2013;Yurchenko et al. 2014;Barton et al. 2014). In this study, we only use the latter to ensure a complete compatibility with the TauREx retrieval code (Waldmann et al. 2015b,a). Unless stated otherwise, our simulations only include H 2 , He, H 2 O, and CO. We also take into account H 2 -H 2 and He-H 2 continua and the atmospheric scattering. To simplify the analysis, we did not include TiO and VO to focus on the molecules above.
The abundances of the species in the atmosphere are calculated based on the temperature given by the GCM simulation and following Parmentier et al. (2018) when dissociation occurs as described in Sect 2.1. The molecular weight then derives from these abundances. The gravity and the atmospheric scale height vary with the altitude from the planetary surface, i.e. the radius of the planet at a pressure of 10 bar, as described in Caldas et al. (2019), and are defined, respectively, in Eq (1) and (2): where g 0 is the surface gravity, R p is the planetary radius, T is the temperature and µ is the mean molecular mass.
We define the volume mixing ratio (VMR) as:  ) without H 2 dissociation. Temperature (top) and the water abundance (bottom) for equatorial cut (left), limb cut (middle), and pole cut (right). From center outward, the 5 solid lines are respectively the 1, 434.10 7 , 10 3 , 1, 10 −2 , and 10 −4 Pa pressure levels. The colormap for water abundance maps goes from 5.10 −4 to 10 −7 . Note that the radius of the planet and the atmosphere are shown to scale. We clearly see that the dissociation of H 2 mostly affects the day side of the atmosphere which increases even more the dichotomy between the day and the night side.
where N i is the molecular number density in molecules per volume unit. We also assume that all the others species are present in the atmosphere as trace gases, so we can write the number density of atomic Hydrogen: where N 0 H 2 is the molecular number density when no dissociation occur -e.g. deep in the atmosphere (see Tab. 1).
We ran Pytmosph3R at R=10000 resolution, with R = λ ∆λ , in the spectral range from 0.6 µm to 10 µm. The spectra are then binned down to a resolution of R=100 for the retrievals. We use this low spectral resolution for several reasons: 1. We simulate JWST observations from 0.6 to 10 µm using the low resolution prism mod provided by Near Infrared Spectrograph (NIRSpec) and Mid Infra-Red Instrument (MIRI) instruments (Stevenson et al. 2016). This does not entail that using the prism configuration is the best observational strategy for such a bright target, but provides us with a uniformly sampled spectrum that does not arbitrarily put more weight in the retrieval on some spectral regions; 2. We compared TauREx with the fully homogeneous case planet, as in Caldas et al. (2019), and we saw that the optimum resolution for the best accuracy in the retrieval is the resolution R=100. 3. There are only two absorbing species in our atmosphere (H 2 O and CO) and we do not need to resolve specific lines, but large features in large spectral bands; In all our simulations we used the planetary and stellar parameters shown in Tab 2. In the part of the atmosphere where H 2 does not dissociate, we have a fixed He H 2 ratio at 0.25885. Then, we compute the volume mixing ratios of H and He by combining the equations (3) and (4) where H 2 dissociates.
Since the overall flux is dominated by the star, we assume that the spectral noise is dominated by the stellar photon noise, defined as: where λ 1 and λ 2 are the limiting wavelengths of the bin considered, d is the distance of the star (270 pc for WASP-121) and R * , T * are, respectively, the stellar radius and the effective temperature. D, τ and ∆t are, respectively, the telescope diameter, the system throughput and the integration time, whose values has been fixed for JWST according to Cowan et al. (2015). Using Eq 5, the uncertainty varies from 10ppm to 50ppm between 1 and 10µm, assuming a single WASP-121b transit. Since systematics may prevent us from reaching a 10ppm precision with JWST, wherever the shot noise was lower than 30ppm we assumed a floor noise of 30ppms through the whole spectral domain (Greene et al. 2016).
Note that we compute different 3D structures as input for Pytmosph3R with the parameters described in table 1. We will describe in details those structures in section 3.

Retrieval model: TauREx
The spectral retrievals has been computed with the TauREx retrieval code (Waldmann et al. 2015b,a). TauREx consists in a Bayesian retrieval code which uses the ExoMol line lists (Yurchenko et al. 2011;Tennyson & Yurchenko 2012;Barton et al. 2013;Yurchenko et al. 2014;Barton et al. 2014), HITRAN (Rothman et al. 2009;Gordon et al. 2013) and HITEMP (Gordon et al. 2010). TauREx assumes a plane parallel exoplanetary atmosphere with a one-dimensional distribution elements and temperature-pressure profile. We assume the radius at the base of the model of Wasp-121b as the value for which the atmospheric pressure is 10 bar. We simulate the exoplanetary atmosphere assuming a pressure between 10 6 and 10 −4 Pascals. We consider an atmosphere dominated by Hydrogen and Helium, with a mean molecular weight of 2.3 amu. We also suppose the presence of CO and H 2 O as trace gases.
An important assumption is that TauREx does not take into account thermal dissociation of species and assumes the molecular abundances constant with the atmospheric pressure (see Changeat et al. (2019) for a detailed discussion on the effect of this assumption). We set up the prior range for the molecular volume mixing ratios of CO and H 2 O between 10 −12 and 10 −1 which allow us to explore a large range of solutions. Even if we never add clouds in our input simulation, we set up the prior range for clouds between 10 −12 and 10 −1 to allow the retrieval code to have more freedom. We also feel that this better simulates a real situation where clouds would be put in the retrieval without prior knowledge about their occurrence in the real planet. The allowed range for the planetary radius and the atmospheric temperature (assumed constant in the vertical) are, respectively, 1.3R Jup − 4.0R Jup and 400K − 3600K, with R Jup the radius of Jupiter. The range for the temperature is linked to the minimal and maximal temperature in the GCM simulation of Wasp-121b as shown in Tab 1. As we simulate JWST observations, we set up the range of wavelength between 0.6 to 10 µm. For the non-retrieved parameters such as masses and stellar radius, we use those describe in Tab 2.
We note that the TauREx retrieval code is not designed to take into account the thermal dissociation of H 2 . When H 2 thermal dissociation occurs, the He/H 2 ratio changes as a function of pressure and temperature. The mean atmospheric weight can also be lower than 2amu. In order to be consistent across all our sets of simulations, we used the same configuration file for our spectral retrievals and we left the He/H 2 ratio constant. By doing so we could estimate whether an atmosphere where H 2 dissociation does not occur can explain our input spectrum or not. We will come back to this in the Sect 4.2.3.
Again, some of these assumptions may seem too simplistic in view of the physics we envision in the real atmospheres of UHJs, but we think that it is important to understand how such simplistic assumptions -commonly used nowadays -can bias the conclusions from a retrieval analysis.

Input models: from simple case to GCM simulations
In order to understand and to characterize the biases due to the atmospheric 3D structure of Wasp-121b atmosphere, we decided to start from a simple case, and add, as the study progresses, further levels of complexity. Using this approach, we will show how to distinguish the observed biases. Note that we used the same atmospheric composition in all our cases, which is He, H 2 , H 2 O and CO. The validation of Pytmosph3R 's model with TauREx was done by Caldas et al. (2019). So, we started by considering the atmospheric 3D structure of an idealized Wasp-121b . In Figs 3 and 4 we show the 3D temperature maps used by Pytmosph3R to  First, we create an idealized case of Wasp-121b atmospheric structure that is symmetric around the planet/observer line of sight. The main parameters to simulate this simple structure are T d and T n and β, respectively the day-side, the night-side and the opening angle. T d and T n are correlated respectively to the sub-stellar and the anti-stellar point. The opening angle β describes how the transitions between the day and the night side is computed, the smaller the sharper. We constructed our model by going symmetrically from the day side to the night side temperature considering a linear transition between the two. To avoid as much as possible effects due to vertical temperature gradients, every column of our model follows an isothermal profile above a given pressure level (P iso = 0.13 bar). In this part of the atmosphere, the temperature however varies with respect to the local solar elevation angle, α * . In the deeper parts of the model, in accordance with results from a GCM for Wasp-121b (see Sect. 2.1), the temperature is uniform and set to 2500 K. In summary, for a given pressure and local solar elevation angle, the temperature is given by: For Wasp-121b we considered two end-member cases: 1. The high gradient model (∇ + T ), which goes from 500 K to 3500 K using β = 10 • . This one is supposed to be representative of the day-night gradient present in the GCM without the heterogeneities in the vertical direction and along the limb. 2. The low gradient model (∇ − T ), which goes from 1400 K to 2800 K using β = 20 • . This case is more representative of GCM simulations accounting for the thermal effect of H 2 recombination. We anticipate this case to be conservative and to provide a lower bound for the real bias caused by day-night temperature differences. Fig. 3 and Fig. 4 show the temperature structure of the atmosphere of the cases studied. These maps are plotted to scale. Note that we consider a symmetric atmospheric structure, so the equatorial cut and the polar cut are exactly the same. As we can see, the day side of those atmospheres are strongly inflated in comparison to the night side which will create the horizontal biases that we want to characterize.
The ∇ + T and ∇ − T cases have been chosen to encompass the real day-night gradient present in Wasp-121b (see below) while removing effects due to differences along the limb (for example, east vs. west limb or equator vs. poles). This should allow us to identify the specific effect of these differences by comparing the retrieval results with our idealized cases and with the actual GCM outputs.
The ∇ + T case is close to the GCM case.  and Tan & Komacek (2019) showed that the energy redistribution effect of hydrogen recombination cools the day side, heats the night side and reduces the gradient across the limb. The ∇ − T case thus uses a reduced gradient to account for this effect and present a very conservative estimate of the 3D effect.

Composition effects due to the dissociation of some molecules
In UHJs, molecular abundances are not the same everywhere, since the day-side temperature is high enough to allow thermal dissociation of many species (Lodders & Fegley 2002;Visscher et al. 2006Visscher et al. , 2010Marley et al. 2017). In this work, we take into account the dissociation of two molecular species, first only H 2 O and later also H 2 , using the equations provided by Parmentier et al. (2018) as shown in Fig. 5. At a given temperature, the lower the pressure, the stronger the dissociation and at a given pressure, higher the temperature higher the dissociation. In this study, we assume that dissociated H 2 predominantly forms atomic hydrogen. We know that hydrogen anion can absorb part of the electromagnetic radiation between 1µm and 10µm (Lenzuni et al. 1991;Parmentier et al. 2018). However, we did not include in our computations the absorption contribution from ionized hydrogen. The main reason is that we want to be consistent with TauREx, which does not include this opacity source. Furthermore, the dissociation of CO cannot occur at the temperature regimes in our simulations, i.e. 500K to 3500K. In order to dissociate, CO requires higher temperatures, due to its stronger triple bond (Lodders & Fegley 2002). For this reason, we assumed in all the simulations, a constant CO abundance everywhere in the atmosphere. In the first set of sim-  Fig. 3: Atmospheric structure of the symmetric, idealized case ∇ − T , assuming an absence and the presence of H 2 dissociation (left and right respectively). We show the temperature map (top) and the water abundance map (bottom). The temperature gradient goes from 1400 K to 2800 K. The transmission angle is 20°around the terminator line and have a 2500 K ring at P=0.13 bar from the surface pressure. From center outward, the 5 solid lines are, respectively, the 10 6 , 10 3 , 1, 10 −2 , and 10 −4 Pa pressure levels. The colormap for water abundance maps goes from 5.10 −4 to 10 −7 . Note that the radius of the planet and the atmosphere are shown to scale.
ulations, we considered the ∇ + T and ∇ − T models considering water dissociation alone (H 2 is held constant). As shown in the abundance map of H 2 O in Figs. 3 and 4, there is no more water in the day-side of the planet at pressure lower than P iso . The dissociation stops around the terminator and, in the night-side, the water distribution comes back at a constant value.
In the second set of simulations we took into account both H 2 and H 2 O dissociation. As H 2 does not have significant absorption bands in the 0.6 -10µm wavelength range, the main effects of its dissociation is to decrease the molecular weight, causing an increase of the atmospheric scale height. It also affects the thermochemical properties of the atmosphere because atomic hydrogen recombines on the night side. This redistribution of energy will heat the night side and cool the day side Tan & Komacek 2019). An other important effect is that without H 2 in the atmosphere, collision induced effects from both, H 2 -H 2 and He-H 2 collisions, are less significant. As shown in Figs. 3 and 4, H 2 dissociation increases the atmospheric scale height of the day side by a factor of 1.5 and 2 in the ∇ − T and ∇ + T simulations, respectively. We note that, the temperature is not high enough to induce H 2 dissociation in the night side. This characteristics, creates an even stronger difference between the night side and the day side when H 2 dis-sociates. Consequently, CO is the dominant radiatively active species in the day-side because it does not dissociate at these temperatures.

3D models of Wasp-121b atmosphere using SPARC/MITgcm global circulation model
In the last configuration, we add both the vertical and azimuthal (along the limb) effects. Indeed, we model Wasp-121b using the Global Circulation Model (GCM) SPARC/MITgcm Parmentier et al. (2018) described in Sect 2.1. The chemistry is kept unchanged. As in Sect. 3.1, we studied four different cases considering all the possible combinations where water and H 2 dissociate. Figs. 1 and 2 show the temperature and water abundance maps for the GCM cases with and without H 2 dissociation. We see that the thermal structure of the atmosphere is more complex than the one in the symmetric case, even if the salient features remain the same. The most important modification concerns the terminator which is now asymmetric because of the winds created by the strong temperature dichotomy between the day and night side . Those equatorial jets warm the east side of the atmosphere and cools the west side around the equator, which implies less dissociation of water on the west  Fig. 4: Atmospheric structure of the symmetric, idealized case ∇ + T , assuming an absence of H 2 dissociation. We show the temperature map (top) and the water abundance map (bottom). The temperature gradient goes from 500 K to 3500 K. The transmission angle is 10°around the terminator line and have a 2500 K ring at P=0.13 bar from the surface pressure. From center outward, the 5 solid lines are, respectively, the 10 6 , 10 3 , 1, 10 −2 , and 10 −4 Pa pressure levels. The colormap for water abundance maps goes from 5.10 −4 to 10 −7 . Note that the radius of the planet and the atmosphere are shown to scale. The temperature and abundance maps of the ∇ − T case are almost the same with only a smaller scale height of the atmosphere. because of the cooler temperature. Nevertheless, the abundance maps show that there is still water dissociated at the terminator which could allow the light to probe deeper into the atmosphere on the night side of the planet. Note also that when H 2 dissociation is taken into account, the global behavior of the 3D structure remains the same, but since the day side of the atmosphere has a bigger scale height, it is more inflated.
An important caveat here is that while we model H 2 dissociation to compute the transit spectrum, the dissociation is not taken into account in the GCM itself. This is expected to reduce the day night contrast thus the winds fields Tan & Komacek 2019). So this case is not completely self consistent.

Results
In this section, we describe the results of the transmission spectra generated by Pytmosph3R and the spectral retrieval using our three set of simulations, each one being performed on our three temperature structures (∇ − T , ∇ + T , and GCM):

The 3D transmission spectrum of ultra hot Jupiters
The transmission spectra generated by Pytmosph3R are shown in Fig. 6. Here, we tested two cases with H 2 constant considering or not water dissociation in the atmosphere (Fig. 6a), and a third one with both H 2 O and H 2 dissociation (Fig. 6b). When H 2 O is constant, the water features are more evident and the CO's features are barely visible as it is highlighted comparing the transmission spectra with and without CO in the atmosphere (respectively the blue and the dotted black lines). Since there is CO and H 2 O everywhere, then the transmitted light probes the same region for both molecules, which are at a high altitude in the day side of the atmosphere. In this case then the spectral features are similar to the ones of the 1D case.
When H 2 O dissociates (light blue curve in Fig. 6a), the water features in the spectrum become smaller because it almost disappears from the day side of the atmosphere. Then, the light from the star can probe much deeper regions of the atmosphere in the wavelength range where the CO does not absorb and reaches part of the night side. The CO, which is not dissociated on the day side -where the scale height is large due to the high temperatureappears more evidently. So, as the light goes first through the day side and then through the night side, the transmission spectrum 3 0 0 0 K 2 7 5 0 K 2 5 0 0 K 2 2 5 0 K 200 0 K 3 0 0 0 K 2 7 5 0 K 2 5 0 0 K 2 2 5 0 K 20 00 K  Fig 5 and Fig 6b. Indeed, the main effect of the dissociation of H 2 in the atmosphere is the increasing of the scale height. As the dissociation of H 2 mainly occurs in the day side of the atmosphere as well as the dissociation of H 2 O (see Fig 1 and 2) the region probed is the same in both cases. However, the CO features appears even more clearly. As CO abundance remains constant everywhere in the atmosphere, its features are more significant due to the enlarger of the scale height caused by the dissociation of H 2 in the day side.
Note that Fig. 6 represents the transmission spectra for the case with the GCM input, but the whole behaviours explained above remain similar to the symmetric ∇ + T and ∇ − T cases.

H 2 and H 2 O constant in the atmosphere
The posterior distributions of ∇ + T and GCM simulations are shown in Fig 7. The retrieval code TauREx converges, in every simulation, to a non-degenerate solution. Fig. 8 (Fig 8). This ratio allows us to understand if the dissociation of water affects the spectrum. The [CO]/[H 2 O] ratio is constant when water does not dissociate, which it is found here by TauREx retrieval. The temperatures retrieved suggest that the light is probing the day side of the atmosphere. Indeed, as shown in Fig 4 and Fig 1, the atmosphere presents an inflated day side. The temperature retrieved in the ∇ + T case is higher than the one retrieved in the ∇ − T case which is expected since the temperature of the day side are higher in the ∇ + T case. Using as input the spectrum computed with GCM, TauREx finds a slightly lower temperature than the one found in the two symmetric simulations. Considering that the GCM simulations has been made assuming a non-isothermal temperature-pressure profile, the temperature structure of GCM is indeed more complex than the one in the other two simulations. These results are all consistent with Caldas et al. (2019). Note that the calculated χ 2 is much higher for the GCM case compared to the symmetric cases ∇ + T and ∇ − T .

Thermal dissociation of H 2 O with H 2 constant
When water dissociation is taken into account, TauREx returns significantly different solutions compared to the homogeneous simulations. The results are shown in the yellow parts of Fig. 8.
First, all the retrieved temperatures are lower than the ones found in the constant composition simulations. Because the main constraint that allows TauREx to retrieve the temperature is the amplitude of the spectral bands -the higher the temperature the higher the amplitude -the temperatures retrieved are lower because the water features are coming from a colder region of the atmosphere. In the GCM case, the temperature retrieved is lower than the temperatures of the limb, which suggest that we probe the night side of the atmosphere.
Since CO does not dissociate, it absorbs mainly in the dayside of the atmosphere, where the high temperature leads to larger absorption features as explained in Sect 4.1. Looking at the yellow part of Fig 8, there are non-degenerate solution for TauREx to fit the CO features, considering that it also finds a low temperature because of the presence of water in the cold nightside of the planet. This solution suggests a high CO abundance in the atmosphere and converges around the correct water abundance within 3 σ. This water abundance corresponds indeed to the abundance in the night side of the atmosphere where water is not dissociated (see Figs 1, 3 and 4).
As a consequence, the [CO]/[H 2 O] is around 100 times higher than expected. We note that for the conservative case ∇ − T the [CO]/[H 2 O] is less biased, around 10 times higher than the solar abundance, due to the smaller day to night contrast and the smoother transition between them. To quantify the reliability of the solutions given by TauREx , we plotted the absorption contribution plots of the solutions for the GCM case and the differences in ppm between the generated and the retrieved spectrum in Fig 9. We highlight here in which wavelength bands the fit is good or not. Moreover, the calculatedχ 2 which is below 1 for each simulation, demonstrates a high significance level of the solution given by TauREx (see Fig 8).

Thermal dissociation of H 2 and H 2 O
In Fig 8 we Fig. 6: (Left): Transmission spectra of Wasp-121b at resolution of R = 100 for GCM simulations assuming a constant H 2 abundance in the whole atmosphere. When water dissociation is taken into account (light blue line), the water features become shallower compared to when we assume no water dissociation (blue line). (Right): Transmission spectra of Wasp-121b at resolution of R = 100 for GCM simulations taking into account H 2 Odissociation in the atmosphere. When H 2 dissociation is considered (blue line), the CO features appears more clearly compared to when we neglect H 2 dissociation (light blue line). Black line and grey line correspond to the transmission spectra for an atmosphere without CO for the water constant and dissociated case respectively to highlights the features of CO in the other curves.
constant composition in the entire atmosphere and in the case where only water dissociates in the day side, TauREx finds a consistent solution, even though it does not necessarily converge towards the correct input parameters. On the contrary, when H 2 dissociation is taken into account, TauREx cannot find a suitable solution to explain the atmospheric spectrum. As shown in the lower left panel of Fig. 8, the best reduced χ 2 for in ∇ + T and GCM models with H 2 dissociation are respectively around 30 and 8. The χ 2 test can give us a possible signal that in the planet we are taking into account it occurs H 2 dissociation. However, TauREx manages to find statistically significant solution in the ∇ − T case, but we keep in mind that this simulation is a conservative choice and we think that the truth remains in between ∇ − T and ∇ + T simulations.
We note that the TauREx retrieval code is not designed to take into account the thermal dissociation of H 2 . When H 2 dissociation occurs, the He/H 2 ratio changes as a function of pressure and temperature. The mean atmospheric weight can also be lower than 2amu in the atmospheric regions where H 2 dissociates. In order to be consistent with all the set of simulations, we used the same configuration file for our spectral retrievals and we left the He/H 2 ratio constant. By doing so we could estimate whether an atmosphere where H 2 dissociation does not occur can explain our input spectrum or not.
The H 2 dissociation plays a major role on the spectra and the retrievals results. All the results are shown in the light yellow part of Fig 8. In each simulation, the [CO]/[H 2 O] ratio retrieved is higher than the one calculated when H 2 is constant. As CO abundance remains constant everywhere, its features are, then, more significant as explained in Sect. 4.1 and TauREx finds a best fit model with a low water abundance and a high CO abundance. As a consequence, the [CO]/[H 2 O] ratio increases.
As shown in the light yellow part in Fig 8, when we consider H 2 dissociation in the GCM case TauREx finds non-degenerate solution. The temperature retrieved is a bit higher than the tem-perature we find when H 2 does not dissociate. Indeed, as shown in Fig 6, H 2 dissociation does not alter the amplitude of the water features dramatically. This is because the water absorption comes from the cold night side of the planet where molecular H 2 dominates. Thus, the temperature retrieved still corresponds to the temperature of the night side.
For the other retrieved parameters, we are consistent with the solution found when H 2 is constant. We remember that the dissociation of H 2 affects the transmission spectrum in mainly 2 ways: i) it increases the scale height of the day side of the atmosphere implying stronger features for species who remains in the day side such as CO. We only consider this effect in the ∇ + T simulation; ii) it also significantly affects the thermochemical properties of the atmosphere due to the recombination of H 2 in the night side which redistribute the energy, heating the night side and cooling the day side Tan & Komacek 2019). We took into account both of them in the ∇ − T simulation. H 2 dissociation results also in a decrease of H 2 -H 2 and He-H 2 collisions, hence less intense continuum absorptions. However, in the ∇ + T simulations, TauREx does not manage to converge to a reasonable solution. TauREx cannot explain the line distortion using a 1D atmospheric model and tends to increase the temperature far above the equilibrium temperature of the planet, reaching the upper prior range temperature, i.e. 3600K. (see Fig. 8).
When both H 2 and H 2 O dissociate, it is not possible to find any 1D atmospheric model that matches with the input spectrum in any of the cases under study in this paper.

Discussions
Our results show that a 1D model is always able to find a statistically consistent solution -except for the case in which H 2 dissociation occurs -although it does not always converge towards the correct input solution. Fig. 7: Nest posteriors retrieval from TauREx assuming a constant H 2 abundance in the atmosphere. There are the ∇ + T (left) and GCM cases (right) which assume, respectively, no water dissociation (Top) and water dissociation (Bottom). We retrieved five free parameters which are H 2 O, CO abundances in log10(VMR), clouds pressure in bar, the temperature in Kelvin and the planetary radius in Jupiter's radius. The molecular weight is derived from those parameters. Caldas et al. (2019) showed that day-night temperature contrasts lead to a bias in retrieved temperatures (toward the one of the day side) even when the chemical composition is uniform. In this work, we extend this by showing that when we take into account H 2 O dissociation, we have both composition and temperature biases.
Although a retrieval code can use a 1D atmospheric model to detect the presence of a particular chemical species from an atmosphere with a complex 3D structure, measuring its abundance is more challenging. When temperature and chemical composition vary across the limb, the 1D retrieval cannot find the correct molecular abundances: the best fit parameters can be several orders of magnitude different from the correct input ones, as shown in Fig 8. Elemental ratios (such as C/O, C/N, etc.) are a key parameter determining the chemical processes taking place in planetary atmospheres (Line et al. 2013;Oreshenko et al. 2017). At high temperature, Madhusudhan (2012); Espinoza et al. (2017) showed that the C/O ratio plays a major role in the atmospheric chemistry and non-solar values could explain observations for 6 hot Jupiters . However, these elemental Const. H2O Diss. H2O Diss. H2O -Diss. H2 The fact that our retrieval results are consistent for both the GCM model and our more idealized models of WASP-121b where temperature is constant in the vertical and symmetric around the substellar point reveals that, for UHJs at least, the effect of thermal and compositional heterogeneities across the limb dominate over the vertical ones as well as the ones along the limb. Fixing an atmospheric configuration -i.e. in presence or not of H 2 and/or H 2 O dissociation, corresponding to the red, dark yellow, and light yellow areas in Fig 8 -  its minimum value when molecular dissociation does not occur, increases by a factor 10-100 when H 2 O dissociates, and by a factor 25-1000 when H 2 dissociates as well, irrespective of the details of the vertical structure of equator-to-pole thermal contrasts. This large scale in the [CO]/[H 2 O] ratio retrieved are due to our idealized ∇ + T and ∇ − T simulations which represent extreme cases, thus they encompass the biases.
A general characteristics of our results is that, in all the cases under our study, we cannot converge around the ground truth parameters. When we study a transmission spectrum extracted from a 3D atmosphere even a high signal to noise spectrum does not allow you to know the real chemical distribution of elements with a retrieval code which uses a 1D atmospheric approximation.
A second aspect is that, with aχ 2 -test it is possible to understand whether there is some more complex physical phenomenon in the atmosphere under study. In our case, we see that aχ 2 > 1 was due to the presence of H 2 dissociation in our atmosphere. From this point of view, aχ 2 -test can reveal whether there is an important physical phenomenon that we are not considering.
Note that a limitation of our study concerns the chemical components of the atmosphere. We only consider two absorber species, CO and H 2 O, but we know that other species have been detected in Wasp-121b such as TiO or VO. Plus, we also know that we are hot enough to dissociate those species ) which would add more complexity in the retrieval analysis.

Are there hints of 3D structures in real data?
In Fig. 10 we compare the spectrum of Wasp-121b as observed with HST/WFC3 (Evans et al. 2016) with two sets of synthetic transmission spectra: i) the models developed by Parmentier et al. (2018) that extract the columns at the terminator of their GCM simulation to compute the transmission spectrum, and ii) our models based on the same GCM simulation but with our more complete 3D radiative transfer framework. In addition to the radiatively active species discussed above, we account here for absoptions by TiO, VO, Na and H − using the analytical fit for the thermal dissociation of the species in Parmentier et al. (2018). Compared to Parmentier et al. (2018), only FeH is missing, but it is believed to be less important in this part of the spectrum since we assume a solar composition.
In order to fit HST/WFC3 data, Parmentier et al. (2018) suggests to add CaTiO 3 clouds in the atmosphere to increase the absorption at 1.25 µm and fill the water window region. From our analysis, we see that considering a cloud free 3D structure can equivalently explain water amplitude in the WFC3 bandpass (thick blue line) even if our spectrum does not fit very good the data either. 3D effects can, indeed, deform the shape of the spectrum. The reason is however quite different. Instead of fitting the 1.25 µm window region, here the agreement is reached by reducing the strength of the 1.15 and 1.4 µm water bands thanks to water removal on the day side of the atmosphere. The fact that both models can have the same observed transit depth at these two wavelengths in Fig. 10 is just due to the fact that the radius at the base of the model remains a free parameter allowing us to shift the whole spectra vertically. Note that, expect for this parameter, there is no fitting of the abundances or thermal structure involved in our approach. We shown that because of the water dissociation, 3D effects was not negligible for atmospheres as hot as wasp-121b, even if clouds remains a reasonable assumption to fit the data. This findings suggests that both 3D heterogeneities in the transmission spectra and clouds can be combined to explain the reduction of the water feature. However, this is depending on which altitude the clouds remain and where they are in the atmosphere because even if there is clouds, the light rays could not probe them due to the 3D structure of the atmosphere.

Conclusion
The 3D structure of the planetary atmospheres plays a major role in shaping of the transmission spectra, especially for ultra hot Jupiters. The tidally locked planets have an inflated day side due to both the high temperature and the dissociation of species which increase the scale height of the day side atmosphere. In these atmospheres, the transmitted spectrum is affected by an extended region around the terminator line, which includes part of both the day side and the night side of the atmosphere.
The stellar light probes an atmospheric region which extends significantly toward the limb, depending on the chemical composition as a function of wavelength. If the temperature is not high enough to dissociate all the molecules, such as CO, they will remain everywhere in the atmosphere. Thus, the amplitude of their features in the transmission spectrum will be larger because they come from the hot regions if the inflated day side atmosphere. On the other hand, the molecules which dissociate more easily, such as H 2 O, will only remain in the cold night side of the atmosphere. Therefore, the amplitude of the features of those molecules in the transmission spectrum will be less pronounced, since they come from colder regions in the night side of the atmosphere. A 1D retrieval code, which tries to fit the spectrum with an isothermal atmosphere, cannot unravel the information of a more complex temperature distribution and, then the abundances retrieved are biased. In particular, we demonstrate that for UHJs 1D retrieval would be biased towards high [CO]/[H 2 O] ratio, the higher the ratio the stronger the chemical heterogeneities. For instance, the TauREx retrieval code finds the planetary temperature according to the strength of the main spectral features. Since the amplitude of the spectral features is rather large because of the gas present in the hot day side of the planet, it is natural for a retrieval code to converge towards hotter temperatures. In case of H 2 O dissociation, the only relevant gas present in the hot day-side is the CO. The water absorption is due to the water present in the cold night side of the planet. Therefore, it is reasonable to retrieve a lower temperature in the case in which water dissociates in the day side.
We demonstrated that a three-dimensional atmospheric structure induces spectral distortions impossible to explain with a 1D retrieval. We also demonstrate that the water features can be strongly reduced by the presence of molecular heterogeneities due to a large day/night temperature contrast. For instance, we show that a 3D geometry can explain some of the features observed in the HST/WFC3 wavelength range of Wasp-121b , even in absence of hazes or clouds, thanks to the dissociation of water in the day side of the atmosphere. However, our works do not exclude the presence of clouds in UHJ atmospheres, the 3D effects are a complementary contribution to fill the water windows in the data. Thanks to our results, we think that Wasp-121b would be so a good target for future observations. We propose that 3D structural effects should already be consider when we study hot and ultra hot Jupiter to avoid erroneous conclusions. Moreover, future space missions such as JWST (Beichman et al. 2014) and ARIEL (Tinetti et al. 2018) will probe a very large range in wavelength (from 0.6 to 28 µm for JWST) and, then 3D effects will be even more evident in other parts of the exoplanetary spectra which has not been observed yet.
Finally, after the analysis of an exoplanetary atmosphere with a 1D retrieval tool, theχ 2 -test on a 1D retrieval can raise a warning about an important effect that the 1D model cannot consider, such as effects induced by the 3-dimensional structure of the exoplanetary atmosphere. In such cases, a parametrized 2D retrieval approach may be warranted. Whether there is enough information in the spectrum to actually constrain such a 2D approach remains to be demonstrated.