Robin Baeyens

Leen Decin

Ludmila Carone

Olivia Venot

Marcelino Agúndez

Paul Mollière

Grid of Pseudo-2D Chemistry Models for Tidally-Locked Exoplanets. I. The Role of Vertical and Horizontal Mixing

Keywords: hydrodynamics -planets and satellites, atmospheres -planets and satellites, composition -planets and satellites, gaseous planets

The atmospheres of synchronously rotating exoplanets are intrinsically three-dimensional, and fast vertical and horizontal winds are expected to mix the atmosphere, driving the chemical composition out of equilibrium. Due to the longer computation times associated with multidimensional forward models, horizontal mixing has only been investigated for a few case studies. In this paper, we aim to generalize the impact of horizontal and vertical mixing on the chemistry of exoplanet atmospheres over a large parameter space. We do this by applying a sequence of post-processed forward models for a large grid of synchronously rotating gaseous exoplanets, where we vary the effective temperature (between 400 K and 2600 K), surface gravity, and rotation rate. We find that there is a dichotomy in the horizontal homogeneity of the chemical abundances. Planets with effective temperatures below 1400 K tend to have horizontally homogeneous, vertically quenched chemical compositions, while planets hotter than 1400 K exhibit large compositional day-night differences for molecules such as CH 4 . Furthermore, we find that the planet's rotation rate impacts the planetary climate, and thus also the molecular abundances and transmission spectrum. By employing a hierarchical modelling approach, we assess the relative importance of disequilibrium chemistry on the exoplanet transmission spectrum, and conclude that the temperature has the most profound impact. Temperature differences are also the main cause of limb asymmetries, which we estimate could be observable with the James Webb Space Telescope. This work highlights the value of applying a consistent modelling setup to a broad parameter space in exploratory theoretical research.

INTRODUCTION

The atmospheres of exoplanets hold great potential for answering questions about a planet's history and habitability. Indeed, accurately determining their bulk chemical composition can help constrain the conditions of planet formation, and provide insight into the potential migration history (Öberg et al. 2011;Madhusudhan et al. 2016;[START_REF] Notsu | [END_REF]). In addition, determining the atmospheric composition is a crucial step in achieving one of the ultimate goals of exoplanet astronomy: detecting biosignatures on extrasolar planets (Schwieterman et al. 2018). We are entering the era of comparative exoplanetary science, as developments of new observational facilities, such as the James Webb Space Telescope ★ E-mail: robin.baeyens@kuleuven.be (Greene et al. 2016) and Ariel (Tinetti et al. 2018), will result in progressively more accessible atmospheric detections. Therefore, understanding the complex physical and chemical processes occurring in the atmospheres of exoplanets will be essential to properly interpret new observational results.

The atmospheres of tidally locked hot Jupiters, i.e. giant closein planets that are the most suitable for atmospheric observations, are highly dynamic environments due to the strong one-sided irradiation from their host star. Acting as a planetary heat engine, this irradiation is the driver for atmospheric circulation with very fast (∼km/s) wind flows, which redistribute the heat from the day side toward the night side of the planet (Koll & Komacek 2018). These high wind speeds have been detected with high-resolution transmission spectroscopy in the hot Jupiters HD 209458 b (Snellen et al. 2010) and HD 189733 b (Brogi et al. 2016;Flowers et al. 2019). The dynamic climate is tightly coupled to the atmospheric chemistry through the temperature and the wind field that causes mixing. This interplay potentially yields very different chemical environments between the day and night sides (Ehrenreich et al. 2020). All these aspects make exoplanet atmospheres intrinsically three-dimensional and motivate the need for sophisticated modelling tools.

General circulation models (GCMs), i.e. three-dimensional hydrodynamics codes used to simulate the large-scale climate of planets, have been frequently applied to tidally locked gaseous exoplanets (e.g. Showman et al. 2009;Rauscher & Menou 2012;Charnay et al. 2015;Amundsen et al. 2016;Mendonça et al. 2018a;Carone et al. 2020). A prominent result of many hot Jupiter GCMs is the consistent prediction of a fast equatorial jet stream, advecting air, and thus heat, eastward. The mechanism for this jet stream was identified to be the up-gradient pumping of angular momentum through planetary waves (Showman & Polvani 2011;Tsai et al. 2014), and its existence has been observationally inferred from the offsets in the peak flux emission of several hot Jupiter phase curves (e.g. Knutson et al. 2007;Komacek et al. 2017;Zhang et al. 2018, and references therein). The 3D structure of the atmosphere and the heat advection through fast winds have been demonstrated to introduce temperature gradients and chemical heterogeneities between the day-night terminators, having an important effect on transit observations (Caldas et al. 2019;Pluriel et al. 2020).

Studies with chemical kinetics models -initially onedimensional -have demonstrated that disequilibrium chemistry significantly changes the atmospheric mixing ratios compared to local equilibrium chemistry (e.g. Moses et al. 2011;Venot et al. 2012;Tsai et al. 2017). Vertical mixing, if efficient enough, can cause quenching in the molecular abundances of prevalent species like CO and CH 4 . If this happens, chemical reactions are too slow to achieve chemical equilibrium, and the gas gets replenished dynamically, potentially with dramatic changes in the atmospheric composition at pressures lower than a critical level. Furthermore, photochemical reactions in the upper atmosphere day side of these highly irradiated planets can cause dissociation of molecules like CH 4 and NH 3 , especially in cooler planets, where the reaction rates are comparatively slow and these photochemically active molecules tend to be favoured (Moses 2014). As a next step, pseudo-2D chemical codes employ the same chemical kinetics, but in addition introduce a longitude-dependent background temperature and uniform horizontal advection (Agúndez et al. 2012;Agúndez et al. 2014a). It has been demonstrated that species can also be efficiently mixed horizontally, in addition to vertically, as they are advected by the equatorial jet stream (Agúndez et al. 2014a;Venot et al. 2020b;Moses et al. 2021). This can result in contamination of the night side by photochemically produced species, such as HCN, and can cause a horizontal homogenization of the chemical composition in general. As observational efforts are increasingly directed towards cooler objects, accurately taking into account disequilibrium chemistry will become more important when deriving bulk elemental abundances (Line & Yung 2013;Baudino et al. 2017;Venot et al. 2020b).

Due to the high computational cost associated with both GCMs and chemical kinetics models, previous efforts of coupling dynamics and chemistry in 1D, 2D or 3D models have often relied on simplifying assumptions. An ubiquitous, parametrized approach of implementing dynamic disequilibrium chemistry is to approximate vertical mixing as a one-dimensional diffusion process, tuned by the eddy diffusion coefficient 𝐾 𝑧𝑧 . Previous studies have either considered 𝐾 𝑧𝑧 to be a free parameter (e.g. Miguel & Kaltenegger 2014;Drummond et al. 2016;Tsai et al. 2017), or derived a numer-ical value for 𝐾 𝑧𝑧 from the Eulerian mean overturning circulation (Heng et al. 2011b) or from mixing length theory (e.g. Cooper & Showman 2006;Moses et al. 2011;Agúndez et al. 2012;Venot et al. 2012;Lee et al. 2015;Rimmer & Helling 2016;Kitzmann et al. 2018;Molaverdikhani et al. 2020;Shulyak et al. 2020). In the latter case, 𝐾 𝑧𝑧 = 𝑤(𝑧)𝐿(𝑧) is a commonly used assumption, where 𝑤(𝑧) is the root mean square of the horizontally averaged vertical velocity, and 𝐿 (𝑧) is a characteristic mixing length scale, often taken to be a (fraction of) the atmospheric pressure scale height. It has, however, been argued that the mixing length approximation is only applicable in convective regimes, and is thus too crude to describe the complex motions in the radiative dynamic atmospheres of hot Jupiters. Parmentier et al. (2013) have equipped their GCM with passive tracers, which are advected by the flow, in order to better constrain the amount of (vertical) mixing in the canonical hot Jupiter HD 209458 b. Likewise, Charnay et al. (2015) have used the method of passive tracers to investigate the mixing in GJ 1214 b. Both studies found 𝐾 𝑧𝑧 values up to two orders of magnitude lower than those derived from mixing length theory, and have derived 1D parametrized, pressure-dependent expressions for the vertical mixing strength in these planets. Some theoretical studies have aimed at refining and generalizing the theory of vertical mixing, by deriving analytic expressions for the quench level in brown dwarfs (Bordwell et al. 2018), and for the 𝐾 𝑧𝑧 of both fast-rotating planets (Zhang & Showman 2018a), and tidally-locked, irradiated planets (Zhang & Showman 2018b;Komacek et al. 2019). These studies have shown that 𝐾 𝑧𝑧 generally depends on properties of the chemical species, such as the mean chemical lifetime 𝜏 chem , in addition to the (local) velocity field. Furthermore, Komacek et al. (2019) have demonstrated a positive trend between 𝐾 𝑧𝑧 and equilibrium temperature, as well as different regimes of vertical transport, depending on whether the atmospheric climate is governed by a fast equatorial jet stream or day-to-night-type winds. Recently, statistical evidence has been found for variable mixing efficiencies across the range of hot Jupiter atmospheres (Baxter et al. 2021).

Another, more direct method of coupling chemistry to the dynamic atmosphere has been to implement the interconversion of CO, CH 4 and H 2 O, dominant opacity sources in hot Jupiter atmospheres, in GCMs by prescribing a fixed chemical relaxation time-scale towards equilibrium (Cooper & Showman 2006;Drummond et al. 2018a,b). This approach was expanded upon by Mendonça et al. (2018b) by using different and more accurate chemical relaxation times for each molecule, as determined by the rate-limiting reaction (Tsai et al. 2018). Finally, the most recent development has been the full coupling of a GCM with a reduced chemical kinetics scheme (Drummond et al. 2020). These studies demonstrate the potential of disequilibrium chemistry to homogenize the atmospheric composition both vertically and horizontally. A caveat to chemistry-coupled GCMs is that, due to their limited resolution, only the large-scale advection of chemical species is taken into account. This could underestimate the mixing efficiency if microphysical or sub-grid scale processes, such as turbulence, are important (Menou 2019).

As forward models of exoplanet atmospheres are becoming increasingly more sophisticated and complex, many studies have focused on a small number of exoplanets to apply these models to. Most prominently, these include HD 209458 b (Showman et al. 2008(Showman et al. , 2009;;Rauscher & Menou 2012;Mayne et al. 2014;Agúndez et al. 2014a;Helling et al. 2016;Drummond et al. 2018aDrummond et al. , 2020;;Carone et al. 2020), HD 189733 b (Showman et al. 2008(Showman et al. , 2009;;Agúndez et al. 2014a;Helling et al. 2016;Drummond et al. 2018b;Steinrueck et al. 2019), and WASP-43 b (Kataria et al. 2015;Mendonça et al. 2018a,b;Venot et al. 2020b;Carone et al. 2020;Helling et al. 2020). One major advantage of this practice is the ability to immediately measure the model against past simulation results, which serve as a benchmark. Such model hierarchy is crucial to investigate the complex physical and chemical processes in exoplanet atmospheres.

In order to explore the diversity in exoplanet atmospheres, ensemble studies using GCMs with varying levels of complexity have been carried out. This is especially the case for atmospheric dynamics. Showman et al. (2015) have conducted a parameter study by varying the rotation period for nine 3D climate models of (non-)synchronously rotating giant exoplanets, and Kataria et al. (2016) have generated consistent climate models for nine known hot Jupiters. Furthermore, some studies have endeavoured to assimilate climate model ensembles into analytic theories (Perez-Becker & Showman 2013; Komacek & Showman 2016;Zhang & Showman 2017;Komacek et al. 2019). These studies have established a variety of climate types and highlighted the potential variations in chemistry and cloud composition. Regarding exoplanet atmospheric chemistry, however, due to the computation times involved in multi-dimensional chemical kinetics, most conducted ensemble or parameter studies are one-dimensional (e.g. Moses et al. 2013;Miguel & Kaltenegger 2014;Heng & Lyons 2016). As such, they do not include the potentially important effects of zonal (i.e. longitudinal) quenching or the zonal dependency of photochemistry. To date, one chemical kinetics ensemble study has incorporated zonal mixing and photochemistry, in order to assess the impact of atmospheric temperatures and chemistry on the phase curves of sub-Neptunes (Moses et al. 2021).

In this work, we aim to explore the chemical diversity of exoplanet atmospheres over a wide range of temperatures, gravities and climate regimes, with a specific focus on disequilibrium chemistry due to atmospheric dynamics. In order to find a trade-off between simulation complexity and computation time, we employ a sequence of post-processed models: a 3D GCM with simplified radiative transfer, and a pseudo-2D chemical kinetics code with a reduced chemical network. With this methodology, we present a large grid of 144 climate and chemistry simulations, incorporating zonal and vertical mixing. In this work, we focus exclusively on the effects of zonal and vertical mixing, and we aim to investigate the additional effect of photochemistry in detail in a follow-up paper.

This paper is organized as follows. In Section 2 we present the sequence of modelling tools that is used for our atmospheric simulations. In Section 3, the grid and the selected parameter space are discussed. The simulation results are presented in Section 4 and further discussed in Section 5. In our discussion we focus on the impact of zonal advection (5.1), limb asymmetries (5.2), case studies of selected exoplanets (5.3), observing windows for some molecules (5.4), deep atmospheric quenching (5.5), and brown dwarf-white dwarf binaries (5.6). Furthermore, some of the shortcomings of the models used in this work are discussed in Section 6. Finally, our conclusions are presented in Section 7.

MODELLING TOOLS

For the forward modelling of our grid of planetary atmospheres, we use a range of modelling tools, further discussed below. Each tool relates to a different aspect of the atmosphere, namely radiativeconvective equilibrium (petitCODE, Mollière et al. 2015Mollière et al. , 2017)), dynamical climate (MITgcm, [START_REF] Adcroft | [END_REF]), chemical composition (pseudo-2D chemistry code, Agúndez et al. 2014a), and synthetic spectrum (petitRADTRANS, Mollière et al. 2019). These codes are used in a post-processing sequence, meaning that there is no iteration or feedback between them. This is done in order to conserve the computational efficiency. We discuss the limitations of the post-processing approach in Section 6.1. To achieve cohesion in the modelling, all parameters are passed down the modelling chain consistently, as shown schematically in Fig. 1. Hereby, the number of free parameters is restricted where possible. In Section 3.2, the parameter dependencies are discussed in more detail.

Radiative-Convective Equilibrium

The first step is the construction of several 1D pressure-temperature profiles for each planetary atmosphere. This is achieved with pe-titCODE, a 1D plane-parallel radiative transfer code for exoplanet atmospheres (Mollière et al. 2015(Mollière et al. , 2017)). Starting from a gas mixture of solar metallicity and a carbon-to-oxygen (C/O) ratio of 0.55, the code converges the pressure-temperature profile to a valid solution, under the assumption of radiative-convective equilibrium and chemical equilibrium. We only consider opacities due to absorption, and neglect scattering processes.

The molecular opacities of H 2 , C 2 H 2 , H 2 S (Rothman et al. 2013), H 2 O, CO, CO 2 , OH (Rothman et al. 2010), CH 4 (Yurchenko & Tennyson 2014), NH 3 (Yurchenko et al. 2011), HCN (Harris et al. 2006;Barber et al. 2014), andPH 3 (Sousa-Silva et al. 2015) have been included in the computation, as well as opacities for Na and K (Piskunov et al. 1995), and collision-induced absorption of H 2 -H 2 and H 2 -He. We note that the opacities of TiO, VO, H -and several (ionized) metals have been excluded, despite their possible importance for the thermal structure of strongly irradiated planetary atmospheres (Arcangeli et al. 2018;Lothringer et al. 2018;Parmentier et al. 2018). The reason for this is twofold. First, the prevalence of these absorbers, especially TiO and VO, and their role in creating thermal inversions in the atmospheres of hot irradiated planets is still not fully understood (see e.g. Piette et al. (2020) for a recent discussion). Since our goal is primarily to analyse the chemistry for a broad range of planets in a consistent manner, we opt to leave out these absorbers, rather than complicate the picture. Second, the inclusion of these opacities tends to give rise to numerical difficulties when using petitCODE, especially when the angle of incident radiation is very shallow. The omission of TiO, VO, H -, and several metal-opacities could lead to an underestimation of the day-side temperature and radiative forcing in our models for planets with 𝑇 eff ≥ 2200 K. We return to the implications of excluding these opacities in Section 6.2.

For each planetary atmosphere in the grid, we follow a similar methodology in constructing a full radiative-equilibrium temperature view of the atmosphere as Showman et al. (2008) and Carone et al. (2020). More specifically, for a fixed orbital distance multiple petitCODE models are computed with varying values of 𝜇 = cos 𝜃, where 𝜃 is the angle of incident stellar irradiation. We have adopted values of 𝜇 = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 and 1.0. For the night side of the planet (𝜇 = 0.0) we computed a petitCODE model without any irradiation from the host star. Finally, these models have been gathered so that their temperatures reach the same common deep adiabat, using the entropy-matching method (see Appendix B). This approach results in a 3D view of the radiative-equilibrium temperature of the planet.

Furthermore, we again use the same method as Showman et al. (2008) and Carone et al. (2020) for determining the radiative timescales 𝜏 rad . We use petitCODE to perturb each converged, pressuredependent temperature profile. Next, the radiative time-scale at that particular temperature and pressure is computed with

𝜏 rad ( 𝑝, 𝑇) = Δ𝑇 𝜌𝑐 𝑃 Δ𝐹/Δ𝑧 , (1)
where the Δ𝑇 = 10 K is the temperature perturbation, 𝜌 is the density, 𝑐 𝑃 is the heat capacity and Δ𝐹/Δ𝑧 is the discretized vertical gradient of the flux perturbation at that atmospheric layer. Together with the radiative-equilibrium temperature profiles, the radiative time-scales are then used as input for the radiative forcing of the GCM.

Dynamical Climate

The dynamical climate of the atmosphere is computed with a General Circulation Model (GCM), more specifically the MITgcm 1 [START_REF] Adcroft | [END_REF]. While this GCM was originally developed for the Earth climate, it has been successfully applied to both rocky exoplanets (Carone et al. 2014) and hot Jupiters (Showman et al. 2009;Carone et al. 2020). With MITgcm, we solve the primitive hydrostatic equations on a cubed-sphere grid with a horizontal resolution of 192 × 32, which corresponds to a resolution of 128 × 64 on a latitude-longitude grid. This yields grid cells of about 2.8 • wide. Furthermore, we use 45 vertical layers: forty layers are equally spaced in log-space between 0.1 mbar and 200 bar, and five additional layers with a thickness of 100 bar are added at the bottom of the atmosphere until 700 bar. This setup is the same as that of Carone et al. (2020) and yields a resolution of about three layers per atmospheric scale height. Our climate model is forced using Newtonian relaxation of the 1 mitgcm.org local temperature towards a prescribed radiative-convective equilibrium temperature during the radiative time-scale:

𝑞 𝑐 𝑝 = 𝑇 rad-eq (𝜙, 𝜆, 𝑝) -𝑇 (𝜙, 𝜆, 𝑝) 𝜏 rad ( 𝑝, 𝑇) , ( 2
)
where 𝑞 is the heating rate per unit volume, and 𝑇 and 𝑇 rad-eq are the local model and radiative-equilibrium temperatures respectively, with a dependence on the latitude 𝜙, longitude 𝜆 and pressure 𝑝.

The radiative-equilibrium temperatures and radiative time-scales are derived from the precomputed petitCODE profiles. A sample of petitCODE radiative-equilibrium temperature profiles, used as input for different climate models, is displayed in Appendix D.

Rather than using the actual temperature, in the code the model is stepped forward using the potential temperature 𝜃, an adiabatically corrected temperature for a fluid parcel that is lifted from a reference pressure 𝑝 0 (in our model 𝑝 0 = 700 bar) to a new pressure 𝑝:

𝜃 = 𝑇 𝑝 0 𝑝 𝜅 . (3)
Here, 𝜅 = 𝑅/𝑐 𝑝 is the ratio of the specific gas constant divided by the heat capacity at constant pressure. In the literature (e.g. Holton & Hakim 2013a;Mayne et al. 2014), the potential temperature can be found expressed using the Exner function Π = 𝑝 𝑝 0

𝜅

as 𝜃 = 𝑇/Π. The code is then stepped forward with

𝐹 𝜃 = 𝜃 eq (𝜇, 𝑝) -𝜃 (𝜙, 𝜆, 𝑝) 𝜏 rad ( 𝑝, 𝑇) , (4)
where 𝐹 𝜃 is the radiative forcing term in the energy equation 𝑐 𝑝 Π d𝜃 dt = 𝐹 𝜃 . For each time step, the irradiation angle parameter 𝜇 = cos 𝜙 cos 𝜆 is determined based on the horizontal location of the fluid parcel, assuming that the substellar point is located at (𝜙, 𝜆) = (0 • , 0 • ) and the day side at -90 • < 𝜆 < 90 • . The associated radiative-equilibrium temperature profile is the precomputed petitCODE profile with that 𝜇-value, rounded up. For the whole night side, when |𝜆| > 90 • , the (𝜇 = 0.0)-profile is used. The radiative time-scales 𝜏 rad ( 𝑝, 𝑇) are computed with an interpolation between temperatures, for a fixed pressure, on the precomputed petitCODE grid using a third-degree polynomial. Finally, in order to prevent numerical problems such as unphysical temperatures, we apply a measure to set the temperature equal to the night-side equilibrium profile, whenever it would become negative during the temperature forcing time step. This numerical problem can sometimes arise on the night sides of highly irradiated and fast rotating planets, during the spin-up of the model. Overall, this setup is very similar to the one used by Showman et al. (2008) and especially Carone et al. (2020). We note that the Newtonian relaxation scheme is a simplified method of radiative forcing, but one that has been extensively used to investigate the climates of hot Jupiters, because of its transparency and computational efficiency. This point is further discussed in Section 6.1.

Additionally, we applied friction to our climate simulations. Following the reasoning of Liu & Showman (2013); Komacek & Showman (2016); Carone et al. (2020), we have used a basal drag that serves to anchor the bottom boundary of the simulation domain to the interior of the planet. This practise was demonstrated to reduce the dependence on initial conditions (Liu & Showman 2013) and numerically stabilize the simulation (Carone et al. 2020). We adopt the same description for basal friction as Carone et al. (2020), using a Rayleigh-type friction proportional to the horizontal wind velocity v h :

F v = - 1 𝜏 fric v h . (5)
Here, the friction time-scale 𝜏 fric decreases linearly with pressure, from its maximal value of 𝜏 fric = 1 day at 𝑝 = 700 bar to a value of zero at 𝑝 < 490 bar. The qualitative dynamical state is not sensitive to the exact value of 𝜏 fric (Liu & Showman 2013;Komacek & Showman 2016) and we adopt the value used by Carone et al. (2020). Furthermore, we applied a fourth-order Shapiro filter with a time-scale of 25 seconds, which acts as an explicit numerical diffusion mechanism (hyperviscosity) to smooth grid-scale noise, as is customary with (exoplanet) GCMs (e.g. [START_REF] Jablonowski | Numerical Techniques for Global Atmospheric Models[END_REF]Showman et al. 2009;Heng et al. 2011a;Carone et al. 2020).

All our simulations were run with fixed values for the specific gas constant and heat capacity, namely 𝑅 = 3589 J kg -1 K -1 and 𝑐 𝑝 = 1.3 •10 4 J kg -1 K -1 respectively. This is appropriate for a hydrogen-helium composition atmosphere with mean molecular weight 2.3. The chemical composition of the air does not come into play directly in the Newtonian-forced GCM, but only enters through the radiative temperatures and time-scales computed with petitCODE. All simulations are performed for a fixed radius at the bottom of the model of 9.7•10 4 km (1.35 R Jup ), the radius of a typical inflated hot Jupiter, so that the only differences between the grid simulations are the gravity, rotation rate and radiative forcing.

Following the recommendation of Sainsbury-Martinez et al. ( 2019), we initialized our model with a hot adiabatic temperature profile. This allows for a comparatively fast convergence of the model to a steady state, even for the deep atmospheric layers, which are known to exhibit prohibitively long convergence times (additionally, see the discussions in Mayne et al. 2017;Carone et al. 2020;Wang & Wordsworth 2020). From there, all models were run for 1500 (Earth) days with time steps of 25 seconds, and the output has been averaged over the final 100 days of the simulation. The model convergence is further discussed in Appendix A.

Chemical Composition

To compute the chemical composition of the atmosphere, we use the pseudo-two-dimensional chemical kinetics model developed by Agúndez et al. (2014a). It involves an effectively one-dimensional vertical column of gas rotating along the equator of the planet. In this column, the coupled vertical continuity equations for each chemical species are integrated, incorporating vertical transport and production and loss rates from thermochemical kinetics (Venot et al. 2012;Agúndez et al. 2014a). The chemical reaction rates are computed with respect to a fixed background temperature which is discretized in longitude (zonally) and pressure (vertically). After a certain integration time, which is computed based on the zonal wind speed, the whole vertical column rotates along the equator as a solid body, thereby mimicking horizontal advection by a uniform equatorial jet stream. The integration of the coupled system of continuity equations (eq. (1) in Agúndez et al. 2014a) is then resumed for the new physical situation. This process is continued for a number of rotations around the planet, until the chemical mixing ratios reach a periodic steady state. For the most abundant atmospheric constituents this typically occurs after several tens of rotations. The result is a pseudo-2D, longitude-and pressure-dependent, chemical composition of the exoplanetary atmosphere. The pseudo-2D framework for chemical kinetics (Agúndez et al. 2014a;Venot et al. 2020b) represents an intermediate step between 1D photochemical kinetics models (Moses et al. 2011;Venot et al. 2012) and fully coupled 3D chemical kinetics (Drummond et al. 2020).

For the thermochemical reactions used to compute the production and loss rates, we have used the reduced chemical network by Venot et al. (2020a), which is available from the Kinetic Database for Astrochemistry (KIDA, Wakelam et al. 2012). It is an updated version of the original network used by Venot et al. (2012); Agúndez et al. (2012); Agúndez et al. (2014a), in order to resolve the existing discrepancies between different chemical kinetics schemes concerning the quenching level of CH 4 (see also Moses 2014;Wang et al. 2016). Furthermore, the number of chemical species and reactions has been reduced, with the aim of a numerically efficient coupling to dynamical models. The network includes 44 chemical species, consisting of C, H, N and O, and contains 582 reactions (288 reversible and 6 irreversible reactions). It is applicable to both cold and hot, hydrogen-dominated exoplanet atmospheres and brown dwarfs, and it has been shown to reproduce the dominant neutral species of those atmospheres (H 2 O, CH 4 , CO, CO 2 , NH 3 , HCN and C 2 H 2 ) accurately (Venot et al. 2020a). The main advantage of the reduced network in this study is the numerical speed-up with respect to a full chemical network: we have measured that the model execution is a few tens of times faster, which allows us to employ it in a wide grid of chemical models. However, as a drawback, discrepancies with the full chemical network have been observed in the upper atmospheres of hot exoplanets (Venot et al. 2020a). Furthermore, the network contains no photolysis reactions and has not been developed with photochemistry in mind, focusing instead on the accurate reproduction of thermochemical equilibrium and the dynamical quenching points. Therefore, it is not suitable for simulating the upper atmospheric chemistry. In this study, we do not use photochemistry, instead aiming to investigate the impact of dynamical disequilibrium chemistry over a wide range of parameters. We aim to address the impact of photochemistry in detail in a follow-up study, using the full version of the chemical network.

For each chemical kinetics simulation, the physical input data, such as temperatures and wind speeds, have been extracted from the corresponding GCM simulation. In the pseudo-2D formalism, the zonal wind speed determines the rotation rate of the vertical gas column, and hence the horizontal advection rate. Additionally, information from the horizontal and vertical wind components is used to compute the degree of vertical mixing within this column, as detailed below. To obtain the longitude-and pressure-dependent temperature background, we extracted an equatorial band from the GCM simulation, averaged and weighted by cos 𝜙 in latitude, for latitudes 𝜙 between +20 • and -20 • . The zonal wind advection speeds have been derived from the GCM simulations as well, namely by averaging the local zonal wind speeds zonally, meridionally between +20 • and -20 • , and vertically for pressures 𝑝 < 10 bar, following the approach of Agúndez et al. (2014a). Averaging zonally over all longitudes is justified for climate regimes dominated by a superrotating equatorial jet stream, the validity of which is further discussed in Section 4.2.1. Most derived zonal wind speeds used in the chemical modelling are of the order of 1 -2 km/s. For some of the coldest planets, the derived zonal wind speeds approach zero, or are very slightly negative, i.e. the zonal mean flow is retrograde. In these cases, we have used a minimum wind speed of 10 m/s in the chemical models, in order to avoid numerical problems. The temperature and wind data resulting from the 3D GCM simulations are presented and discussed below, in Section 4.2.

Vertical mixing is incorporated in the chemical kinetics code as a diffusion process, the efficiency of which is given by a pressuredependent eddy diffusion coefficient 𝐾 𝑧𝑧 ( 𝑝). Although the vertical mixing in irradiated exoplanets is governed by the general atmospheric circulation, and thus it is not diffusive (Zhang & Showman 2018b;Komacek et al. 2019), it is still widespread practice to parametrize vertical mixing using 𝐾 𝑧𝑧 . In this work, we compute 𝐾 𝑧𝑧 using the scaling relation derived by Komacek et al. (2019):

𝐾 𝑧𝑧 ∼ 𝑤 2 1 𝜏 chem + 1 𝜏 adv , ( 6
)
where 𝑤 is the vertical wind speed, 𝜏 chem the chemical time-scale, and 𝜏 adv the horizontal advective time-scale of the wind flow. In this equation, we computed the vertical wind speed value by globally averaging the vertical wind from the corresponding GCM simulation on isobars. While it is conceivable that the vertical mixing efficiency varies with location, e.g. between the day and night sides of a planet, it is difficult to assign a local 𝐾 𝑧𝑧 -value to it, as it has been argued that the approximation of mixing as a diffusive process is mostly valid on the global scale (Parmentier et al. 2013;Showman et al. 2020). For this reason, the wind speeds involved in the 𝐾 𝑧𝑧 computation are globally averaged on isobars.

The horizontal advective time-scale 𝜏 adv represents the efficiency of chemical advection by the wind flow on isobars, and thus the extent to which horizontal chemical perturbations can be maintained. Here it is computed as a function of pressure 𝑝:

𝜏 adv ( 𝑝) = 𝐿 𝑣 ℎ ( 𝑝) , ( 7
)
where 𝐿 is a typical horizontal length scale, taken to be the planetary radius (𝑅 𝑝 = 1.35 𝑅 Jup ), and 𝑣 ℎ is the horizontal wind speed, extracted from the GCM simulation and globally averaged on isobars.

Next, the appearance of the chemical time-scale 𝜏 chem in equation (6) implies that a different 𝐾 𝑧𝑧 -profile is associated with each individual chemical species (Zhang & Showman 2018a,b;Komacek et al. 2019). Although it is possible to derive chemical time-scales through careful analysis of the chemical network by determining the rate-limiting step as a function of pressure and temperature (Tsai et al. 2018), this is outside the scope of this work. Hence, we adopt here a constant chemical time-scale of 𝜏 chem = 1.6 • 10 5 s, which is the intermediate value used by Komacek et al. (2019), and we make an a posteriori estimate of the chemical time-scale based on our chemical network, as a sanity check. Tsai et al. (2018) have performed a comparison of the chemical reaction pathway by which CH 4 is transformed into CO, for three different chemical networks (Tsai et al. 2017;Moses et al. 2011;Venot et al. 2012). At temperatures between 1000K and 1500 K and pressures higher than 1 bar, the rate-limiting reaction in the former two networks is the methanol (CH 3 OH) production reaction

OH + CH 3 + M -→ CH 3 OH + M, ( 8
)
where M is any third species. In the Venot et al. (2012)-network this reaction is not the rate-limiting step, because a different reaction in this network accounts for efficient methanol production. However, the more efficient reaction was removed in the updated chemical network used here (Venot et al. 2020a), making it very likely that the rate-limiting reaction in our network is indeed (8). We estimated the chemical time-scale by evaluating the reaction rate of (8) for different converged chemistry models in our grid. We find values between 10 4 s at 10 bar and 10 10 s at 1 bar for warm to hot planets, while the cold (𝑇 eff < 1000 K) planets tend to have longer chemical time-scales. Thus, our adopted value of 𝜏 chem = 1.6 • 10 5 s seems suitable, although it must be remarked that the chemical time-scale can be expected to change by orders of magnitude depending on the planetary parameters or location inside the atmosphere. Finally, we calibrated this scaling relation to the eddy diffusion coefficient profiles that have been determined using 3D simulations with passive tracers, for the planets HD 209458 b (𝐾 𝑧𝑧 = 5 • 10 4 • 𝑝 -0.5 bar m2 /s, Parmentier et al. (2013)) and HD 189733 b (𝐾 𝑧𝑧 = 10 3 • 𝑝 -0.65 bar m 2 /s, Agúndez et al. (2014a)). To achieve this simple calibration, the 𝐾 𝑧𝑧 -profile obtained from eq. ( 6) was multiplied by a scale factor for the two models in our grid closest to the planets HD 209458 b and HD 189733 b. We found that a scale factor of 0.1 results in a good order-of-magnitude agreement for both cases. Therefore, we use

𝐾 𝑧𝑧 ( 𝑝) = 0.1 𝑤 2 ( 𝑝) 1 𝜏 chem + 𝑣 ℎ ( 𝑝) 𝑅 𝑝 (9)
to compute vertical mixing profiles for our chemical kinetics simulations. The 𝐾 𝑧𝑧 -profiles, and how they relate to the different planetary climate models in our grid, are presented and discussed below, in Section 4.2.3. Finally, we remark that 𝐾 𝑧𝑧 can be expected to rise again in the deep atmosphere (𝑝 > 1 bar), more specifically when the atmosphere becomes convective and hence well-mixed. Our 𝐾 𝑧𝑧 -profiles estimated via (9) are based on atmospheric circulation theory, and do not scale well to the deep convective interior, because the low resolution of most exoplanet GCMs cannot capture convective motions. The uncertainty associated with 𝐾 𝑧𝑧 in the deep atmosphere is further discussed in Section 5.5.

The initial state of the chemical kinetics integration is a converged, one-dimensional chemical kinetics simulation, assuming solar metallicity and C/O ratio, and a pressure-temperature profile corresponding to the substellar point of the planet. This choice of initial state results in an efficient convergence to a periodic steady state (Agúndez et al. 2012). Each chemical model was subsequently run for 100 rotations about the equator. Depending on the different zonal wind speeds used in each model, all chemical models have total simulation times between 10 7 s and 10 10 s. This is sufficient for the evolution from a vertically quenched initial state to a zonally and vertically quenched steady-state in the pseudo-2D framework. The model convergence to a steady state is demonstrated and discussed in Appendix A.

Synthetic Spectrum

As a final step, the results from the GCM (temperatures) and the chemical kinetics simulations (molecular abundances) are assimilated to compute synthetic transmission spectra, using the radiative transfer package petitRADTRANS (Mollière et al. 2019). In the spectral synthesis, we include line opacities for H 2 , C 2 H 2 (Rothman et al. 2013), H 2 O, CO, CO 2 , OH (Rothman et al. 2010), CH 4 (Yurchenko & Tennyson 2014), NH 3 (Yurchenko et al. 2011) and HCN (Harris et al. 2006;Barber et al. 2014), as well as Rayleigh scattering opacities for H 2 (Dalgarno & Williams 1962) and He (Chan & Dalgarno 1965), and continuum opacities in the form of collision-induced absorption from H 2 -H 2 and H 2 -He interactions. For all species, the chemical abundance profiles are first computed with the pseudo-2D chemical kinetics code.

Regarding species not included in the chemical network, such as Na and K, we take an agnostic approach. They are not incorporated in the spectral synthesis, despite their ubiquity in transmission spectroscopy (e.g. Sing et al. 2016). Although it is possible to include these species with chemical equilibrium concentrations 2 , this is against the spirit of this work, namely a consistent study of dynamical disequilibrium chemistry. Hence, we omit these opacities and indicate their absence wherever it is relevant.

Using petitRADTRANS, a synthetic transmission spectrum is then computed with a wavelength range between 0.3 𝜇m and 30 𝜇m, and with a spectral resolution of 𝜆/Δ𝜆 = 1000. To this end, a pressure-grid of 100 layers is set up, evenly spaced in log-space between 100 bar and 10 -4 bar. An atmospheric structure is interpolated on to this pressure grid, consisting of a vertical temperature profile, extracted from the GCM simulation, and molecular abundances and mean molecular weight, derived from the chemical kinetics code. Finally, the effective radius is computed with a ray-tracing method (see Mollière et al. (2019) for more details), where a pressure of 10 mbar is used as a reference pressure for the planetary radius 𝑅 𝑝 = 1.35 𝑅 Jup .

In the computation of a synthetic transmission spectrum we take into account contributions from both the evening and morning limb of the exoplanet disc. We calculate the averaged effective transit radius of the planet as

𝑅 ave = √ 𝑅 2 m + 𝑅 2 e 2 , ( 10
)
where 𝑅 m and 𝑅 e are the effective transit radii computed with peti-tRADTRANS for a one-dimensional, vertical atmospheric structure extracted from the morning or evening limb respectively. We note that this approach in constructing a transmission spectrum is rather crude, since we do not incorporate contributions from the high latitudes of the projected disc, or variations in the atmospheric structure along the line of sight. Instead, the transiting area is here approximated as two half discs (morning and evening), which are assumed to each have a spherically uniform atmospheric structure. Thus, the resulting transmission spectrum will be biased towards the physical and chemical conditions of the equatorial region (between ±20 • latitude). The consequences of including a more realistic 3D approach to computing transmission spectra, including line-of-sight variations in the atmospheric structure, have been examined by Caldas et al. (2019).

GRID OF MODELS

Free Parameters

In our grid of exoplanet atmosphere models, we have three varying free parameters: the effective temperature 𝑇 eff of the planet, the surface gravity 𝑔 of the planet, and the stellar type of the host star.

Other necessary model parameters remain fixed (the atmospheric metallicity, C/O ratio), or they are changed self-consistently (rotation rate, Section 3.2). These parameters and their adopted values are described in what follows.

Effective temperature

The effective temperature 𝑇 eff of the planet is defined as the blackbody temperature of an object with the same radius that would emit a flux equivalent to the real planet. Thus, it can be expressed as combination of the equilibrium temperature 𝑇 eq , resulting from energy balance with the stellar radiation, and the intrinsic temperature 𝑇 int , resulting from radiation originating in the planet's interior (Fortney 2018):

𝑇 4 eff = 𝑇 4 eq + 𝑇 4 int . ( 11
)
The effective temperature is arguably the most important free parameter, affecting the radiative, dynamical and especially chemical time-scales. Therefore, we probe a wide range of temperatures, from fairly cold (400 K) to ultra-hot (2600 K) gaseous exoplanets. We choose a uniform and quite fine temperature spacing of 200 K, to properly diagnose climate and chemical regime changes in the atmosphere. This leads to 12 sampled temperatures: 𝑇 eff = 400 K, 600 K, 800 K, 1000 K, 1200 K, 1400 K, 1600 K, 1800 K, 2000 K, 2200 K, 2400 K, and 2600 K.

Surface gravity

The surface gravity depends on the planetary mass and radius, which display large variations, especially in the Neptune-and Jupitermass regime (e.g. Chen & Kipping 2017). Indeed, hot Jupiters are often inflated and typically have surface gravities of ∼10 m/s 2 (or log 𝑔 = 3 in cgs units). However, some exoplanets have a very high density and thus surface gravity, such as WASP-18b (Hellier et al. 2009), which despite its high temperature of 2400 K has a log 𝑔 of 4.3 (Southworth 2010;Stassun et al. 2017). On the other side of the spectrum, extremely low density planets, so-called super-puffs, can have surface gravities as low as log 𝑔 ≈ 2 (Lee & Chiang 2016;Libby-Roberts et al. 2020).

The effect of surface gravity is mainly twofold. First, through the density it affects the radiative forcing, see eq. ( 1). Second, by changing the atmospheric scale height, the opacity of the atmosphere changes as well, which will have a noticeable effect on the planet's spectrum. In particular, the photospheric pressure scales with gravity, and the molecular features probed in transmission spectroscopy are small if the surface gravity is high. In this work, we opt for a coarse sampling of the surface gravity, choosing three values to represent the diversity exhibited by gaseous exoplanets: log 𝑔 = 2, 3 and 4.

Stellar type

The effect of the host star's stellar type on the planet atmosphere is twofold. First, as noted by Mollière et al. (2015), the radiation coming from hotter stars is absorbed at shorter wavelengths than radiation originating from cooler stars. This causes a higher temperature gradient in the former, and a more isothermal atmosphere in the latter case. Second, by changing the stellar temperature in our grid, while keeping the effective temperature of the planet fixed, we effectively vary the distance to the host star, and thus also the rotation period of the planet. This is further detailed in Section 3.2.

Similar to Mollière et al. (2015), we choose four values for the stellar type, namely M5, K5, G5 and F5 main sequence stars. The corresponding stellar effective temperatures are 3100 K, 4250 K, 5650 K and 6500 K respectively. For the same stellar type, the stellar properties -temperature, luminosity, mass and radius -have been correlated based on table 2.4 of [START_REF] Hubeny | Theory of Stellar Atmospheres[END_REF].

Parameter Dependence

Aside from the three free parameters, our grid of tidally locked atmospheric models is designed to be self-consistent in the adopted input parameters. For instance, the incident stellar flux, the equilibrium temperature and the rotation period of the planet are related under the assumption of synchronous rotation: two planets with the same equilibrium temperatures orbiting stars of a different stellar type will have different semi-major axes and orbital periods, and thus also rotation periods. More specifically, the orbital distance 𝑎 of a planet with an effective temperature 𝑇 eff is given by

𝑎 = √ 𝑓 (1 -𝐴 𝐵 )𝐿 * 16𝜋𝜎(𝑇 4 eff -𝑇 4 int ) , (12)
where 𝐴 𝐵 is the Bond albedo of the planet, 𝐿 * is the stellar luminosity, 𝜎 the Stefan-Boltzmann constant and 𝑇 int the intrinsic temperature of the planet. The heat redistribution factor 𝑓 equals 2 when the incoming radiation is distributed over the day-side, and 1 when it is distributed over the whole surface of the planet. In the calculations for our grid, we have assumed a global heat redistribution, and a Bond albedo of 0, i.e. all incident radiation is absorbed by the planet. Furthermore, we adopted 𝑇 int = 200 K (see Appendix B for a more detailed discussion of the intrinsic temperature).

The rotation and orbital period is then computed via Kepler's third law:

𝑃 rot = 𝑃 orb = √ 𝑎 3 4𝜋 2 𝐺 (𝑀 * + 𝑀 𝑝 ) , (13)
with 𝑀 * and 𝑀 𝑝 the stellar and planetary mass respectively, and 𝐺 the gravitational constant. In the calculation of the orbital/rotational periods in our grid, we have kept the planetary mass fixed to 1 𝑀 Jup , so the rotation period does not vary as a function of gravity in our grid. Nevertheless, the planetary contribution to the total mass of the star-planet system is negligible. As a consequence, under the assumption of synchronous rotation, the rotation period of the planet scales with the planetary temperature as 𝑃 rot ∼ 𝑇 -3 eff . An overview of the computed orbital distance and or-bital/rotation period for each model in the grid, as a function of planet temperature and stellar type, is given in Table 1.

SIMULATION RESULTS

Radiative-Convective Equilibrium

The outcome of the petitCODE simulations are temperaturepressure profiles, corresponding to different angles of incidence, for each simulation in the grid. For each combination of surface gravity and stellar type, the night-side temperature profiles are similar, while the profiles corresponding to different incident angles on the day side increase with the simulation's effective temperature. Each temperature profile is put onto a common adiabat in the deep atmosphere, equal to that of the planetary average temperature, using the method described in Appendix B. This results in the typical structure of irradiated exoplanets consisting of an irradiated zone, in which the temperature depends strongly on the incidence angle, connecting to a common, deep adiabatic temperature gradient (Showman et al. 2008;Guillot 2010;Carone et al. 2020). A sample of the radiative-convective equilibrium temperature profiles, used as input for the GCM modelling, is displayed in Appendix D.

Dynamical Climate

Overview of climate regimes

The majority of the MITgcm climate models in our grid of tidallylocked gaseous exoplanets have a climate regime dominated by equatorial superrotation, i.e. a fast (∼km/s) prograde jet stream around the equator; see Figures 2 and3 for 𝑔 = 10 m/s 2 . (The equivalent figures for 𝑔 = 1 m/s 2 and 𝑔 = 100 m/s 2 can be found in Appendix E.) 3 This equatorial superrotation is a well-known and robust outcome of hot Jupiter GCMs (Showman & Polvani 2011;Debras et al. 2020) and it has been established over a broad range of temperatures (Komacek & Showman 2016;Komacek et al. 2017). We thus reconfirm it as the most prominent climate type in our grid.

Some notable deviations to the superrotating climate regime start to appear for very fast-rotating simulations in the grid. These simulations correspond to highly irradiated planets around cool stars, which reside in the lower left corner of Fig. 2. Here, an equatorial retrograde flow starts to appear on the day side. A detailed discussion on this type of wind flow can be found in Carone et al. (2020). Although the zonal mean wind still exhibits equatorial superrotation globally, as indicated by the prograde wind jets at the equator of all models (Fig. 3), locally this retrograde flow is disrupting the superrotating jet stream. In doing so, it converges at the western terminator with the night-side prograde flow, resulting in locally high vertical wind speeds. The models in Fig. 2 that show equatorial retrograde winds at this pressure level all have rotation periods shorter than 8 hours. For the low-gravity models that critical rotation period is even shorter, for the high-gravity models, it is longer. The influence of gravity as additional parameter affecting retrograde flow has been established by Carone et al. (2020) as well.

Additionally, the equatorial superrotation is not always present in the coldest planets of our grid. More specifically, out of all 𝑇 eff = 400 K models in Fig. 2, only the K5-orbiting planet shows signs of equatorial superrotation. The other models exhibit a more direct day-to-night-side flow (G5 and F5) or an equatorial retrograde flow on the day side together with a prograde flow on the night side (M5). The cold models corresponding to low (Fig. E1) and high surface gravity (Fig. E2) show analogous deviations from equatorial superrotation. This may indicate that the radiative forcing, caused by the permanent one-sided heating, for these cold models is insufficient to trigger the standing wave response that is necessary to excite and sustain the equatorial jet stream. The excitation of superrotation in this thermal regime seems to be less robust, and can be expected to depend on many other parameters such as the rotation rate, gravity, and numerical assumptions, like in the cases of Jupiter and Saturn (Showman et al. 2020). It can be questioned, finally, if these very cool, synchronously rotating models provide a proper representation of cool exoplanets, as their synchronization time-scales approach 1 Gyr and higher (see Appendix C).

Another prominent trend is the qualitative change in the zonal wind flow with respect to changing effective temperature and rotation rate. Fig. 3 indicates that the strength of the equatorial jet stream is correlated with the effective temperature of the planet, as the zonal wind speed increases from almost zero for the coldest planets to about 6 km/s for the hottest. Interestingly, the fast-rotating planets (orbiting M5-stars) do not follow this pattern. Instead, the jet stream speed stays constant with effective temperatures increasing above 1200 K. An additional effect of the rotation rate is the variation in the width of the equatorial jet stream for groups of models with the same effective temperature. Indeed, as the rotation period becomes shorter, the latitudinal extent of the jet decreases as well, as was previously found by Showman & Polvani (2011); Showman et al. (2015), who showed that the jet width is proportional to the Rossby deformation radius. The latter quantity depends inversely on the 𝛽-parameter, and thus on the rotation rate. Equivalently, the jet width is reduced due to the exponential decay of the latitudinal extent of the standing equatorial Kelvin waves, needed in the superrotation excitation. The e-folding decay length varies similarly with the 𝛽-parameter (Holton & Hakim 2013b). Finally, we find that the gravity does not have a big impact on the zonal mean wind flow, as indicated by the strong similarities between Figures 3,E3 and E4. Finally, we can discuss the a priori validity of adopting a single, zonally averaged wind speed to represent the zonal advection in the pseudo-2D chemistry code, considering the climate regimes displayed in Fig. 2 and3. As discussed above, the equatorial zonal wind speed is not always uniform. Indeed, longitudinal variations in the horizontal wind speed are noticeable. Most noteworthy are the cases with high rotation rates, where fast retrograde flow can be seen on the day-side. As was hypothesized by Carone et al. (2020), this situation could result in a morning terminator composition that is affected by the day-side state, rather than the night-side state. Moreover, the fast vertical motions associated with the convergence of horizontal wind flow could cause locally high vertical mixing efficiencies, or a dynamical pile-up of clouds. None the less, in the zonal wind plots it can be seen that -on average -the superrotating jet stream dominates the atmospheric circulation, even in the cases with fast retrograde day-side flow. Hence, a uniform zonal wind speed should still be a good, first-order assumption for all cases, except the coolest ones (𝑇 eff = 400 K).

Heat redistribution

The temperature plots of the GCM models, as isobaric map (Fig- ures 2, E1 andE2) or as a function of pressure (Figures 4, E5 andE6), already indicate trends in the difference between the dayand night-side temperatures. To provide a good overview of the heat redistribution in all 3D GCM models, the temperatures of all models have been averaged, and the relative day-night temperature differences are plotted with respect to the planet's effective temperature for different surface gravities and host star types (Fig. 5). In all cases, a clear increasing trend with effective temperature is noticeable, which is in line with previous climate models (Perna et al. 2012;Perez-Becker & Showman 2013;Komacek & Showman 2016) and observations (Cowan & Agol 2011;Komacek et al. 2017). The reason often provided for this trend is that the radiative time-scale and heat advection time-scale both decrease with increasing irradiation, but the radiative time-scale more so. Thus, in these idealized GCM simulations, hotter planets can efficiently radiate away heat advected by the wind flow (see also Showman et al. (2020) for a more comprehensive discussion).

Another parameter influencing the atmospheric heat redistribution in irradiated exoplanets, is the rotation rate. Indeed, planets of the same effective temperature orbiting different stars, and hence having different rotation rates, show clear differences in their day-night temperature contrast (Figures 2 and5). More specifically, higher rotation rates inhibit the heat redistribution, resulting in larger day-night contrasts. This can be understood in terms of the force balance in the hydrodynamic momentum equation between the Coriolis force and the pressure gradient caused by the temperature difference between the day side and night side. In the case of fast rotation, the Coriolis force can support a bigger pressure gradient, resulting in larger day-night temperature contrasts (Komacek et al. 2017). This is clearly demonstrated in Fig. 5, where planets orbiting cooler host stars, with faster planetary rotation rates, exhibit a consistently larger day-night temperature contrast.

Vertical wind speed and 𝐾 𝑧𝑧

The vertical wind speeds in exoplanet atmospheres are a valuable diagnostic to constrain the amount of vertical mixing. As a consequence, understanding the effects of different physical parameters on the vertical wind speeds is key to predict disequilibrium chemistry. Analogous to the temperature maps of Fig. 2, horizontal, isobaric maps of the vertical wind speeds are shown in Fig. 6 for the medium-gravity subsection of the grid. It can be seen that an increase in the effective temperature of the planet results in faster vertical wind speeds. This is to be expected, as the unilateral irradiation is the driver for the atmospheric circulation. The vertical motions at 70 mbar are characterized in general by upwelling flow at the day side, and downward flow occurring in more localized fronts.

Among the most prominent fronts are those forming a 'chevron' pattern at the equator, between 90 • and 180 • on the night side (see also Rauscher & Menou 2010;Zhang & Showman 2018b). They are characterized by convergences in the horizontal wind field, where the wind from high latitudes flows equatorward and meets with the eastward advected air from the day side, resulting in very fast vertical flow in a 'chimney'-like structure (Parmentier et al. 2013;Komacek et al. 2019). Both the location and speed of this vertical chimney seem to be influenced by the Coriolis force, as evidenced by their dependence on the planet's rotation rate (e.g. Fig. 6, 𝑇 eff = 1200 K, along row). Slowly rotating planets display a comparatively fast vertical wind chimney, which is located near the anti-stellar point. On the other hand, planetary atmospheres that rotate faster see these columns shifting westward to the evening terminator. The vertical wind speeds generally also decrease with The relative mean temperature difference between the day-and night side is systematically larger for fast-rotating models (orbiting M5-stars) than slow-rotating models (orbiting F5-stars). Furthermore, it increases with planetary effective temperature. The relative mean temperature difference was computed by subtracting the average night-side temperature from the average day-side temperature, and dividing by the average day-side temperature. To compute the average day-or night-side temperature, we have taken the arithmetic mean of the temperature 𝑇 ( 𝜙, 𝜆, 𝑝) in our GCM model for all pressures 𝑝, latitudes (weighted with cos 𝜙) and longitudes 𝜆 of the day or night side respectively. The three panels show the low-, medium-and high-gravity models.

faster rotation. For the fastest rotators, with 𝑃 rot < 1 day typically, the large-scale flow breaks up into small-scale substructures and the prominence of the chevron front gradually disappears.

Beyond the presence of localized fronts, the vertical wind speeds in general show a clear dependence on the rotation rate, with fast rotators displaying slow vertical wind speeds and vice versa (Fig. 6 along rows). Indeed, the upward wind on the planet's day side reaches roughly 50 m/s in general (not including localized fronts) when the planet is rotating comparatively slow (G5 and F5). For faster rotators (K5), this is already reduced to 10 m/s, with the exception of the dynamically active equatorial band. Finally, in the case of the ultra-fast rotators (M5), vertical wind speeds become negligible in the regions outside the narrow equatorial band. Hence, the rotation rate is critical for a good estimation of the vertical motion in irradiated, tidally-locked exoplanets.

Upon comparing the vertical wind maps of models with different surface gravities, all sub-grids show similar trends with effective temperature and rotation rate/host star type. However, in absolute terms the values are very different (see Fig. 6 for medium, and Figures E7 and E10 for low and high gravities). It is clear that a high surface gravity serves to reduce the vertical wind speeds overall, as the (𝑔 = 100 m/s 2 )-simulations generally do not have vertical wind speeds above 10 m/s. Regarding the low gravity (𝑔 = 1 m/s 2 )simulations, on the other hand, the hot and slowly rotating planets display wind speeds that surpass ±100 m/s on large parts of the isobaric surface plotted in Fig. E7. The inverse dependence between vertical wind speed and the surface gravity is evident, since the vertical wind speed scales with the atmospheric scale height 𝐻 = 𝑅𝑇 𝑔 (Komacek et al. 2019). Hence, high-gravity planets can be expected to have lower eddy diffusion coefficients than their lowgravity counterparts. The vertical advection time-scale can then be estimated through 𝜏 vertical ∼ 𝐻 2 /𝐾 𝑧𝑧 . Given that 𝑤 ∼ 𝐻 and 𝐾 𝑧𝑧 ∼ 𝑤 2 , first-order analytic theory would suggest that the scale height and vertical wind terms cancel out, and the vertical transport efficiency is not dependent on the gravity. However, in our GCM simulations, which span two orders of magnitude in 𝑔, we find that the models with higher gravity generally have lower vertical transport time-scales than models with low gravity (see also Section 4.3.2).

The eddy diffusion coefficients, computed using eq. ( 9) and based on the GCM wind maps, show the same trends as the vertical wind speeds, namely correlated with effective temperature and inversely dependent on the surface gravity (Fig. 7). Indeed, the value of the eddy diffusion coefficient changes by up to four orders of magnitude, depending on the effective temperature of the planet. This is roughly in both qualitative and quantitative agreement with the 𝐾 𝑧𝑧 computations by Komacek et al. (2019) (see their Fig. 7, for 𝜏 chem = 1.6•10 5 s and an infinite drag time-scale). Furthermore, our 𝐾 𝑧𝑧 -profiles show a similar pressure dependence as the power-law parametrization made by Parmentier et al. (2013) (dashed line in Fig. 7). Although the method employed in this work to derive 𝐾 𝑧𝑧 via ( 9) is not as accurate as passive tracers in describing the vertical mixing through atmospheric circulation, our agreement with studies that do involve passive tracers (Parmentier et al. 2013;Komacek et al. 2019) allows us to conclude that this description is sufficiently accurate for the purposes of this study.

While there are clear trends between the order of magnitude 𝐾 𝑧𝑧 value, and atmospheric parameters such as effective temperature and gravity, some of the derived 𝐾 𝑧𝑧 -profiles have anomalous shapes, especially if the effective temperature is low (400-600 K). In these cases, 𝐾 𝑧𝑧 shows relatively little variation with altitude, rather than the steep monotonic increase that characterizes the hotter models. Moreover, in many low-temperature 𝐾 𝑧𝑧 profiles, a localized sharp drop can be seen (e.g. K5, 𝑔 = 10 m/s 2 , Fig. 7). Since the horizontal advection time-scale shows only small variations with pressure, the general morphology of the 𝐾 𝑧𝑧 -profiles is attributed to the mean vertical wind speed. Thus, this sharp drop in 𝐾 𝑧𝑧 is associated with a sign-change in the vertical wind. It appears that the irradiation in the cold models is insufficient to drive a vertically coherent atmospheric circulation. The transition in climate regime is further exemplified in the wind maps (Fig. 2 and3), and the absence of equatorial superrotation as we discussed in Section 4.2.1. It has indeed been demonstrated that a low stellar irradiation results in climate regimes that strongly differ from typical hot Jupiters, with possible transitions between sub-regimes depending on rotation rate (Showman et al. 2015). We conclude that care should be taken when 𝐾 𝑧𝑧 is applied to cool (𝑇 eff < 600 K) exoplanets, as extrapolations from hotter exoplanet atmospheres break down. As the focus in exoplanet science shifts from hot Jupiters to cooler planets, further investigations in the climate regimes of mildly irradiated gaseous planets would be beneficial.

Finally, we note that in the deepest parts of the atmosphere, any clear trends in the 𝐾 𝑧𝑧 -profiles break down. Regardless, the vertical mixing strength in these layers is uncertain, as 𝐾 𝑧𝑧 is expected to increase again when the convective regime is reached (see Section 5.5).

Chemical Composition

Overview

The chemical abundances for some prominent molecules in exoplanetary atmospheres, computed with the pseudo-2D chemical kinetics code, are displayed in Fig. 8 for the simulations with medium surface gravity. The figure gives an overview of the dominant atmospheric constituents, besides H 2 and He. In all cases, it can be seen that H 2 O is one of the most abundant molecules. This is complemented by CO in the hot planets (𝑇 eff ≥ 1400 K), and by CH 4 in the cool planets (𝑇 eff ≤ 1000 K). In between these temperatures, coexistence of H 2 O, CO and CH 4 is possible in equal abundance, but the rotation rate (captured in the host star type) can be seen to play a role in the balance between these three molecules. The transition between CO and CH 4 as the main carbon-bearing species in exoplanet atmospheres is well understood in chemical equilibrium (Lodders & Fegley 2002). Additionally, it is crucial to map this transition in disequilibrium chemistry, since it serves as the basis for a number of simplifying assumptions in exoplanet chemistry studies (e.g. Cooper & Showman 2006;Mendonça et al. 2018b;Drummond et al. 2018a,b;Steinrueck et al. 2019).

Upon viewing the grid of pseudo-2D chemical kinetics results in Fig. 8, a striking dichotomy can be seen. Cool planets with 𝑇 eff < 1400 K do not show longitudinal variations in their chemical composition. Hot planets with 𝑇 eff > 1400 K do show longitudinal variations. This is especially the case for temperaturesensitive molecules, such as CH 4 , NH 3 and HCN, but even H 2 O shows longitudinal variations for those models where the temperature differences become very high, namely the ultra-hot, fast rotating planets. Thus, zonal quenching by the horizontal wind appears to be sufficiently strong to smear all chemical variations on isobars in the low-temperature regime. However, as the effective temperature increases, both the zonal wind speed (Fig. 3) and day-night temperature differences (Fig. 5) increase. The former would keep the chemical composition horizontally homogeneous, whereas the latter would cause chemical variations, as the warm day side tends more towards a chemical equilibrium composition. Above the threshold of 𝑇 eff ≈ 1400 K, it appears that the zonal temperature variation becomes the main cause for zonal chemistry variations. Finally, note that we did not incorporate photochemistry, which could be an additional catalyst for longitudinal chemistry variations, especially in the upper atmosphere (Agúndez et al. 2014a;Venot et al. 2020b).

Vertical mixing ensures an almost constant chemical composition in the majority of the mid to high altitudes, especially in the atmospheric models of cool planets. We observe that, with our description of 𝐾 𝑧𝑧 , the typical quenching pressure lies between 1 bar and 10 bar for the planets with 𝑇 eff < 1600 K. Here, the coldest planets have the deepest quenching points. Thus, although the eddy diffusion coefficient generally increases with effective temperature (Fig. 7), models with higher effective temperatures still have shallower quenching points, because of the faster reaction rates associated with higher temperatures. Furthermore, planets hotter than 1600 K show different turnoff points from their chemical equilibrium composition depending on the longitude probed. Species on the night side are typically still quenched at the 1 bar pressure level, but the quenching level increases as hotter longitudes are viewed. Consequently, for the hottest planetary atmospheres, species at the day-side longitudes remain nearly unquenched and in equilibrium.

In the hottest atmospheres paired with G-and F-stars, the abundances of H 2 O and CO decrease sharply in the deep atmosphere, around 100 bar. The temperatures at this depth reach about 5000 K, which is too hot for most molecules to exist. Hence, H 2 O, CO, but also H 2 revert to their atomic components. A similar phenomenon is seen in the chemical abundances of hot, low-gravity planets (Fig. E9). We note that most NASA coefficients, used to calculate the chemical equilibrium, are only verified up to 5000 K (Mcbride et al. 1993). Moreover, the temperature in the deepest parts of the atmosphere is rather uncertain (see Section 5.5), and transitions in the equation of state can be expected at these conditions (Saumon et al. 1995).

The influence of the rotation/host star type on the chemical composition models manifests itself in three ways. First, a higher rotation rate is related to a slower and narrower equatorial jet (see Section 4.2.1), and thus a lower degree of horizontal quenching in the chemistry models. Second, a higher rotation rate is associated with an inefficient heat redistribution (see Section 4.2.2), resulting in chemical equilibrium compositions that are highly longitudinally dependent. Third, a higher rotation rate correlates to slow vertical wind speeds and consequently lower eddy diffusion coefficients (see Section 4.2.3), leading to a reduced vertical mixing of the chemical species and a lower quenching pressure. In summary, the first and seconds effects complement each other, so that the fast rotating planets show more pronounced horizontal chemical variations, which is confirmed in Fig. 8 (see for instance along rows of 𝑇 eff ≥ 1400 K). The third effect, of rotation leading to reduced quenching levels, can be discerned as well (Fig. 8), in particular by comparing the departure of chemical species from their equilibrium track in the models with 𝑇 eff = 800 K. Here, the quenching level strongly affects the CO-content: a fast rotating planet (host star M5) has a low CO content due to the relatively shallow quenching level, but a slowly rotating planet (host star F5) can be seen to have a deeper quench level, leading to a coexistence of CO, CH 4 and H 2 O. Furthermore, the temperatures of the G5 and F5 models appear to be slightly higher near the quenching pressure, resulting in more CO in equilibrium (Fig. 4). By analogy, for planetary temperatures with 𝑇 eff = 1200 K, CO, CH 4 and H 2 O can coexist if they are fast rotating, but CO and H 2 O dominate if they are slowly rotating. We note that the set of 400 K-models also shows varying levels of COquenching for different host star types. In this case, however, the behaviour is not easily generalized into a trend with rotation rate, due to the more erratic 𝐾 𝑧𝑧 -profiles of these models, highlighted earlier (Section 4.2.3). For hotter planets, in general, the effect of rotation on the quenching level is limited.

Upon comparing the chemical models in Fig. 8 with their low (Fig. E9) and high (Fig. E10) surface gravity counterparts in the grid, agreement with the above trends is found. The dichotomy between horizontally homogeneous models (𝑇 eff < 1400 K) and models showing horizontal chemical gradients (𝑇 eff > 1400 K) is still present. However, the distinction is less clear for the low-gravity case, where a longitudinal dependence is never very pronounced. The high-gravity models, on the other hand, display very strong longitudinal chemistry gradients in the molar fractions of CH 4 , NH 3 and HCN. The increasing chemical inhomogeneity for hot model atmospheres with increasing gravity, is most likely linked to a combination of two effects. First, the heat redistribution is relatively inefficient in high-gravity atmospheres (see Fig. 5), resulting in strong horizontal temperature gradients. Second, in this part of the parameter space, vertical mixing is more efficient than zonal wind advection. This point is discussed in more detail in Section 4.3.2. Thus, the compositional differences set by a longitudinally dependent chemical equilibrium, are maintained throughout the vertical extent of the atmosphere. Finally, the dichotomy threshold temperature of ∼1400 K appears to be independent of the surface gravity.

Time-scale comparison

In order to disentangle and compare the advection effects of vertical mixing and zonal mixing, the typical time-scales associated with both processes can be computed. The vertical mixing time-scale is computed as 𝜏 vertical = 𝐻 2 /𝐾 𝑧𝑧 , where 𝐻 is the atmospheric scale height and 𝐾 𝑧𝑧 the eddy diffusion coefficient computed with (9). The horizontal mixing time-scale is computed as 𝜏 zonal = 2𝜋𝑅 𝑝 /𝑢, where the planet circumference 2𝜋𝑅 𝑝 is taken to be a typical length scale for the zonal wind advection, and 𝑢 is the zonal wind speed. Even though these time-scales are only zeroth order estimates of the mixing efficiencies, they can still provide insight in the relative importance of both mechanisms (see e.g. Drummond et al. 2020).

In the relatively deep part of the atmosphere, up until ∼1 bar, the dominant advection type is zonal mixing (𝜏 zonal < 𝜏 vertical , see Fig. 9). While the zonal advection time-scale stays almost constant with pressure, the time-scale associated with vertical mixing decreases strongly with altitude (see also Fig. 7), so that it becomes shorter than the zonal advection time-scale at lower pressures. The quenching pressures, at which chemical reactions no longer reach an equilibrium, are located in the zonally dominated regime. Consequently, chemical species are first horizontally quenched, and later vertically. This is in agreement with previous studies researching the quenching behaviour of HD 209458 b (Agúndez et al. 2014a;Drummond et al. 2018aDrummond et al. , 2020)).

Generalizing the time-scale analysis of Fig. 9, we have computed the ratio of the zonal and vertical time-scales at 1 bar, for each model in the grid, and plotted the result as a function of the effective temperature, for different model parameter combinations (Fig. 10). It becomes clear that, at 1 bar, the zonal advection time-scale is in most cases shorter than the vertical mixing time-scale. This is a consequence of the fast equatorial jet stream that is excited in almost all models. Thus, if the quenching pressures are located around 1 bar, as is the case for atmospheres with an effective temperature of ∼1400 K or lower, most atmospheres are primarily horizontally quenched, and thus do not exhibit strong chemical gradients on isobars. This is in agreement with our chemical models (Fig. 8). On the other hand, if the quenching pressure is raised to lower pressures, as can be expected for hotter atmospheres, it can be expected that vertical mixing, rather than zonal, will gain importance, given the uniform zonal wind speed and the increase in 𝐾 𝑧𝑧 with pressure.

The zonal-to-vertical mixing ratio seems to be sensitive to the host-star type/rotation rate, although it is difficult to find a clear trend spanning the parameter space. It could be envisaged that slowly rotating atmospheres have zonal-to-vertical mixing ratio's that are relatively high, since they tend to have a faster equatorial jet stream (see Fig. 3). However, they also tend to have higher vertical wind speeds (see Fig. 6), so the relative importance of zonal and vertical mixing is not trivially interpreted. Moreover, in our computation of 𝐾 𝑧𝑧 via eq. ( 6), both effects are already taken into account.

In some regions of the parameter space, vertical mixing is the more dominant process at 1 bar. This occurs in coldest models, for which the equatorial jet stream is insufficiently excited, as well as the slowest rotating (M5-typed), ultra-hot atmospheres. Finally, the models with high surface gravities (100 m/s 2 ) all exhibit a relatively zonal mixing dominates Figure 10. If the log ratio of the zonal advection time-scale 𝜏 zonal and the vertical mixing time-scale 𝜏 vertical is positive, the atmospheric mixing is more efficient in the vertical direction (𝜏 vertical < 𝜏 zonal ). If it is negative, the atmospheric mixing is more efficient in the zonal direction (𝜏 zonal < 𝜏 vertical ). In the ratio computation, the vertical mixing efficiency is sampled at 1 bar, and the zonal mixing efficiency is computed using the uniform zonal wind of the pseudo-2D framework. The three panels show the low-, medium-and high-gravity models.

high vertical mixing efficiency, with many of them having zonalto-vertical time-scale ratio's above 1. The result is a lack of zonal advection for atmospheres with high surface gravity, which is either associated with chemical equilibrium for pressures higher than the quenching pressure, or associated with a zonal-to-vertical timescale ratio above 1 for pressures below the quenching pressure. The corresponding chemical abundance plots likewise show this lack of zonal mixing as strong horizontal gradients in the molar fractions of most species (Fig. E10).

Synthetic Spectrum

A suite of cloud-free transmission spectra, modelled with petitRAD-TRANS, is shown in Fig. 11. The spectra are displayed according to effective temperature, with spectra corresponding to different rotation rates plotted together. We have opted here to plot the effective transit radius (in 𝑅 Jup ), computed from the morning and evening limb contributions (see eq. 10), in order to be able to compare simulations with different rotation rates/host star types directly. The transit depth, on the other hand, which is an observational parameter and is denoted by the ratio of the planetary and stellar radii 𝑅 2 𝑝 /𝑅 2 * , will vary much more among spectra of different rotation rates, because the stellar radius varies self-consistently with the planetary rotation rate in the framework of synchronous rotation.

Upon inspecting the spectra in Fig. 11, from the coldest to the hottest atmospheres, molecular absorption bands can be distinguished for different species, which are in line with the most abundant species in the chemistry models. Water absorption is om-

1.355 1.360 K K K K M5 K5 G5 F5 1.35 1.36 K K K K 1.35 1.36 1.37 K K K K 1.35 1.36 1.37 K K K K 1.350 1.375 K K K K 1.350 1.375
Transit radius (R Jup )

K K K K 1.350 1.375 K K K K 1.350 1.375 K K K K 1.350 1.375 K K K K 1.350 1.375 K K K K 1.350 1.375 K K K K 10 0 10 1
Wavelength (micron)

CO

Figure 11. The synthetic transmission spectra computed with petitRADTRANS for simulations with medium surface gravity (𝑔 = 10 m/s 2 ), ordered by increasing effective temperature from top to bottom. Transmission spectra corresponding to different rotation rates/host star types are colour coded: (red, M5), (orange, K5), (green, G5) and (blue, F5). These simulations, differing in rotation period, have largely overlapping spectra, but in some molecular bands differences can be discerned. A selection of molecular bands and features has been annotated, either above or below the region of their highest opacity. Note that the vertical axes have differing scales. Moreover, the units are 𝑅 Jup for an easy comparison between planets orbiting different stars. As such, the actual transit depth, 𝑅 2 𝑝 /𝑅 2 * , will differ, depending on the stellar radius of the host star 𝑅 * .

nipresent. In the coldest atmospheres (𝑇 eff = 400 K), it is competing with Rayleigh-scattering at the short end of the spectrum (700 -900 nm). As the temperature increases and the contribution of methane diminishes, water absorption bands in the nearinfrared become more pronounced, and at effective temperatures above ∼1400 K, water is the dominant opacity source for the majority of the spectral range. Furthermore, in the mid-infrared, the broad water absorption centred at 6.3 𝜇m is visible in all spectra.

Methane on the other hand, is only visible in the coldest planets.

At 𝑇 eff = 400 K, it is the dominant opacity source in the nearinfrared, with very pronounced molecular bands. As the methane content decreases, however, it becomes more difficult to attribute spectral features to methane or water absorption, due to the large wavelength overlap of their absorption bands. This is illustrated for 𝑇 eff = 1000 K, where both molecules contribute equally to the transmission spectrum. Similarly to methane, ammonia is only visible in transmission spectra corresponding to the coldest atmospheres.

It has a prominent absorption band at 10 𝜇m, which is very pronounced for the 𝑇 eff = 400 K models, but quickly diminishes as the temperature increases. Nevertheless, at this limited wavelength window, it seems to persist in atmospheres with effective temperatures of up to 1400 K. We discuss the outlook of observing ammonia in Section 5.4.3. Finally, molecular signatures of carbon monoxide and carbon dioxide are quite rare, as they both have a relatively poor absorption spectrum in the considered wavelength range. Still, at 4.3 𝜇m, the opacities of both CO and CO 2 are at their maximum, and the resulting opacity bump can be seen in many spectra in Fig. 11. The potential of observing CO and CO 2 is further discussed in Section 5.4.2. We note that the spectra corresponding to atmospheres with a different rotation rate/host star type, are generally very consistent. Nevertheless, a few noteworthy differences in the transmission spectra arise when the rotation rate -and all associated climate and dynamical mixing properties -are changed. Most notably, the coldest model with 𝑇 eff = 400 K shows distinct absorption features at 4.3 𝜇m and 4.6 𝜇m, coming from CO 2 and CO respectively, for stellar types M5 and G5, but these are not present for K5 and F5. This is a strong signature of disequilibrium chemistry for the two former cases. Indeed, the abundances of both CO and CO 2 depend very strongly on the temperature in chemical equilibrium (Fig. 8). Small changes in the quenching level can, consequently, induce order-ofmagnitude-sized variations in the abundances of these species. This results in a spectrum which is very sensitive to the vertical mixing strength.

A similar discrepancy between fast (M5, K5) and slow (G5, F5) rotators can be seen in the 1200 K-spectra. Here, the former show much more pronounced absorption bands of methane than the latter. From Fig. 8 it becomes clear that the methane content in the atmosphere depends critically on the chemical equilibrium composition at the quenching pressure. Indeed, the CH 4 equilibrium abundance slightly decreases as the host star gets hotter. This is caused by the gradual temperature increase with host star type in the deeper parts of the atmosphere (Fig. 4). Looking at the 1200 Kmodels in this figure, the atmosphere goes from a CH 4 -dominated equilibrium (M-star) to a CO-dominated equilibrium (G-, F-star) near the quenching pressure (a few bar). The eddy diffusion coefficient also shows a small gradual increase with host star type for these models (Fig. 7), so the quenching pressure is similar in all four cases. This results in a higher methane fraction throughout the atmosphere of fast rotating planets, and a stronger absorption signature in the spectrum.

Finally, for transmission spectra corresponding to atmospheres of ultra-hot exoplanets (𝑇 eff > 2000 K), the overall transmission depth can be seen to decrease with the rotation rate. Upon closer inspection, we have established that this decrease is mainly caused by the changing temperature structure with rotation rate, rather than the comparatively small influence of water dissociation, which can be seen in Fig. 8. Overall, we show that -in certain cases -the rotation rate of the planet can have a very strong impact on its transmission spectrum, even though it is negligible in most cases. In the future, it could be interesting to see how taking into account the effect that the host star and rotation rate have on the climate and the chemical mixing efficiencies, can potentially bias exoplanet retrieval results, such as the inferred C/O ratio.

DISCUSSION

Zonal Mixing In Context

Previous research into the coupling of chemical kinetics and dynamics, in pseudo-2D (Agúndez et al. 2014a;Venot et al. 2020b) or in 3D (Mendonça et al. 2018b;Drummond et al. 2018aDrummond et al. ,b, 2020)), has highlighted the potential of horizontal (zonal and/or meridional) advection to produce horizontal homogeneous chemical compositions. Indeed, upon comparison with 1D chemical kinetics models, many main molecules in the atmospheres of the exoplanets in consideration no longer show horizontal compositional gradients when horizontal mixing is taken into account. This is especially important for a species like CH 4 , that can show order-of-magnitude variations with longitude when horizontal mixing is not taken into account. Pseudo-2D chemical kinetics models have shown that in such case, contamination from the day side will reduce the overall methane content (Agúndez et al. 2014a;Venot et al. 2020b). For the planet HD 209458 b, the same trend has been confirmed with a 3D GCM coupled to chemical kinetics (Drummond et al. 2020). For the planet HD 189733 b, this study also found a chemically homogeneous atmosphere, although in this cooler planet the methane abundance is expected to be enhanced due to meridional quenching from high latitudes.

To assess the impact of zonal mixing in exoplanetary atmospheres, and disentangle it from other modelling assumptions, we analyse a hierarchy of simulations with increasing complexity. What follows is a short description of each modelling level.

(i) Agnostic 1D average: a one-dimensional radiativeequilibrium temperature profile, computed for planetary averaged irradiation (see red line in Fig. B1). With the denomination 'agnostic', we stress that no preconceived notions about three-dimensional properties of the planetary atmosphere are used in the calculation (as opposed to averaging the results of a 3D GCM simulation). Chemical equilibrium is assumed. The synthetic spectrum is computed for this single vertical structure.

(ii) 3D GCM + equilibrium chemistry: two vertical temperature profiles, corresponding to the morning and evening limbs, extracted from the GCM climate simulation. The profiles are weighted and meridionally averaged about the equator, in the same way as it was done for the pseudo-2D case (see Section 2.3). For each of these, the chemical equilibrium is determined. The synthetic spectrum is computed by taking into account equal contributions from the morning and evening limbs (see Section 2.4).

(iii) 3D GCM + 1D chemical kinetics: same temperature profiles In this model hierarchy, the chemical composition of the morning and evening limbs (left) and the corresponding synthetic transmission spectrum (right) are shown for each modelling assumption. Four modelling assumptions are presented, from top to bottom: First, an 'agnostic' one-dimensional model with equilibrium chemistry is displayed, corresponding to the planetary averaged flux distribution (see also Fig. B1). Second, the pressure-temperature profiles are drawn from the 3D GCM, still assuming equilibrium chemistry. Third, the model is augmented with vertical mixing, resulting in potential departures (solid lines) from equilibrium chemistry (dashed lines). Finally, the introduction of the pseudo-2D framework incorporates horizontal advection of the chemical species, which represents the most sophisticated modelling step employed in this work. All hierarchical steps are based on a planet with 𝑇 eff = 1400 K and 𝑔 = 10 m/s 2 orbiting a G5-star.

as above. For each of these, a 1D chemical kinetics model with vertical mixing is computed. The synthetic spectrum is computed as above.

(iv) 3D GCM + pseudo-2D chemical kinetics: same temperature profiles as above. For the equatorial band, the pseudo-2D chemical kinetics framework is used to compute the chemical composition, incorporating both vertical and zonal mixing. The synthetic spec-trum is computed as above. This is the nominal modelling approach for the grid models in this work.

For a planetary atmosphere of 1400 K, the chemical composition and synthetic transmission spectrum for each of these four modelling steps are shown in Fig. 12. The overlapping transmission spectra demonstrate that the incorporation of the morning and evening limb temperatures, drawn from a 3D climate model, reduces the extent of the molecular features, compared to the simple agnostic model. It can be reduced with respect to a one-dimensional mean temperature profile (Caldas et al. 2019). The reason is that the 3D approach results in colder temperatures for the morning terminator, due to horizontal advection from the night side. Furthermore, retrievals of exoplanet emission spectra, taking into account the 3D information of the day side temperature, result in retrieved temperature profiles that are close to the arithmetic mean temperature distribution (Blecic et al. 2017). Since the limbs will have lower temperatures than the dominating day side, it is reasonable to assume that the temperature probed with transmission spectroscopy will be lower than the arithmetic mean and also lower than the agnostic temperature profile used here.

Additionally, the introduction of disequilibrium chemistry, through vertical and also zonal mixing, manifests itself in the appearance of methane absorption. This is most prominently seen at 3.3 𝜇m. The case with only vertical mixing shows the strongest absorption by methane, and relatively high volume mixing ratios of about 10 -5 for both limbs. When horizontal quenching is taken into account, the mixing ratio and the prominence of the absorption band both drop, as methane is horizontally quenched to the day-side composition, similar to what was found by Agúndez et al. (2014a).

When we apply the same hierarchical methodology to a planet with a cooler (𝑇 eff = 800 K, Fig. 13) and a hotter (𝑇 eff = 2000 K, Fig. 14) atmosphere, a trend starts to emerge. In the case of a cooler atmosphere, which shows methane absorption in the transmission spectrum, the switch from a 1D temperature profile to temperatures derived from a 3D GCM again has the biggest impact on the spectrum. However, subsequent levels of complexity very minorly or negligibly influence the transmission radius. mixing ratios, because of deep quenching, and the increased absorption between 4 and 5 𝜇m accompanying it. Zonal mixing has no additional effect, since the abundances are all quenched quite deep, so that no longitudinal gradients are present.

In the hierarchical case study of a hotter planetary atmosphere, each modelling level results in a distinct composition and spectrum (Fig. 14). In line with the changes seen in the planets of 800 K and 1400 K, the use of 3D GCM temperatures changes the molecular absorption features strongly. Next, vertical mixing quenches the abundances at pressures of 100 mbar to 1 bar, which is higher up in the atmosphere than where it occurs in the planets with lower temperatures, despite more vigorous mixing. In this case, unlike the cool planetary atmosphere of 800 K, horizontal quenching has a large impact on the chemical composition. This is especially the case for CH 4 , reducing its fraction in the atmosphere from more than 10 -5 in the 1D kinetics case to below 10 -6 in the pseudo-2D case, for all pressures at the morning terminator. This change is reflected in the spectrum, where it can be seen that the transit radius reduces in the most prominent methane absorption bands.

In summary, from this discussion it appears that the consideration of more realistic temperature profiles, either drawn from a 3D GCM simulation (e.g. Venot et al. 2020b;Lewis et al. 2020) or self-consistently parametrized (Gandhi & Jermyn 2020), has a profound impact on both the exoplanet chemistry and the transmission signature. This illustrates the importance of taking into account horizontal heat redistribution and terminator asymmetries. Furthermore, the relevance of these effects applies to a wide range of atmospheric temperatures (Fig. 12,13,14). Where the temperatures modify the extent of the transmission features, disequilibrium processes act to increase or reduce absorption in specific wavelength ranges, since they have the potential to alter the concentration of some species by orders of magnitude. In relatively cool planets, vertical mixing can be seen to quench molecular abundances efficiently, as is expected for planets below 1200 K (Line & Yung 2013). Because of this deep quenching, longitudinal gradients in the chemistry are almost non-existent. For this reason, horizontal mixing is quasi-irrelevant in relatively cool planets. In relatively hot planets, on the other hand, horizontal quenching is needed to avoid overpredicting the abundances of species like CH 4 and NH 3 , that tend to favour lower temperatures. The atmospheres of these planets are quenched at lower pressures, which -in combination with their large day-night temperature differences -results in strong horizontal chemical gradients in a 1D framework. Horizontal diffusion acts as a homogenization factor, which in general sets the chemical composition at the limbs closer to the day-side composition.

Limb Asymmetries

The modelling approach used in this work is well suited to explore the presence of morning and evening limb asymmetries in the diversity of exoplanet atmospheres, be it because the limbs differ in temperature, or in chemical composition. To this end, the transmission spectra corresponding to the morning and evening terminators are considered separately, where the atmospheric structure at -90 • and 90 • longitude, respectively, is assumed to constitute the whole projected ring that is probed in transmission spectroscopy (see also Section 2.4).

Detectability

The maximal difference in the transit depths between the morning and evening terminators shows an increase with increasing temperature, for a subset of grid models with medium surface gravity (Fig. 15). Indeed, considering a sample of models of increasing temperature -800 K, 1400 K, and 2000 K -the maximal differences between both limbs are 0.004 R Jup , 0.011 R Jup , and 0.019 R Jup , respectively, for the G5-host star case. These values correspond approximately to 0.3, 0.8, and 1.5 per cent of the radius of a typical inflated hot Jupiter (1.35 R Jup ).

Evidently, the potential for observing such a small difference will depend on the exact physical characteristics of the planet-star system. However, assuming a G5-type host star with a radius of 0.9 R (table 2.4, Hubeny & Mihalas 2014), the estimated transit depth 𝑅 2 𝑝 /𝑅 2 * differences for the three cases shown would amount to 1, 2 and 4 ppm, respectively: too low to detect. Models corresponding to high (𝑔 = 100 m/s 2 ) or low (𝑔 = 1 m/s 2 ) surface gravities generally show a similar correlation with the effective temperature. The amplitude of the limb differences, however, changes, because of the change in scale height. In the low-gravity case, this results in limb asymmetries of 38, 235 and 1951 ppm, which would be comfortably detectable with the James Webb Space Telescope (JWST), which would have a noise floor of 20-50 ppm (Greene et al. 2016).

We note that this represents a very optimistic scenario, since the gravity of most transiting hot Jupiters will be higher than 1 m/s 2 . However, by doing this exercise, we illustrate that the detection of limb asymmetries with JWST should be possible, and hinges on a combination of a comparatively low gravity, a medium to hot atmosphere, and a favourable host star.

It is important to remark that neither clouds nor hazes are taken into account in our models. It has been shown that clouds and hazes can exhibit inhomogeneous distributions (Parmentier et Figure 15. The maximal difference in transit depth between the morning limb (-90 • longitude) and the evening limb (90 • longitude) generally increases with the exoplanet's effective temperature (top panel). The data displayed corresponds to grid models with 𝑔 = 10 m/s 2 , and host star types/rotation rates as indicated with coloured lines. Next, the transmission spectra (grey), along with those corresponding to the morning (blue) and evening (pink) limbs, illustrate in which wavelength ranges limb asymmetries result in the biggest differences (bottom three panels). Spectra for three different grid models are shown: 𝑇 eff = 800 K, 1400 K, and 2000 K. All of them have 𝑔 = 10 m/s 2 and a G5-star host. Powell et al. 2019;Steinrueck et al. 2020), which add an additional layer of complexity in the comparison between the morning and evening limbs. On the upside, the observation of limb asymmetries could provide important information to break retrieval degeneracies (Line & Parmentier 2016;Kempton et al. 2017). Finally, we note that, although our results suggest that ultra-hot (𝑇 eff > 2200 K) Jupiters are the optimal candidates for limb asymmetries, we are careful with this conclusion. Ultra-hot exoplanets are expected to have additional opacity sources (Lothringer et al. 2018;Arcangeli et al. 2018;Parmentier et al. 2018), not included here (see Section 2.1). Moreover, hydrogen gas dissociation on the day side and recombination on the night side has been shown to reduce the daynight temperature contrast in these planets (Tan & Komacek 2019).

Hence, it can be expected that the limb asymmetries associated with these planets are not as large as presented in Fig. 15. This point is discussed further in Section 6.

Thermal asymmetry, but chemical similarity

Upon examining the causes of limb transit depth differences in our models, we find that temperature asymmetries between the morning and evening limbs are much more important than differences in the chemical composition. The morning and evening limbs are at different temperatures because of eastward heat advection. Hot, day-side air and cold, night-side air are both transported eastward, so that the morning limb (-90 • longitude) is generally colder than the evening limb (90 • longitude). This can be readily seen in Fig. 2. Likewise, the separated transmission spectra show that the morning limb has a lower transmission signal in the infrared, being the colder of the two (Fig. 15). We conclude that limb asymmetries are maximized when the heat transport is sufficiently efficient for heat redistribution to take place, but not too efficient, so as to remove all horizontal temperature gradients. Chemical differences between the two limbs, on the other hand, are either not present, or appear to be of secondary importance in the interpretation of limb transit depth asymmetry. As we established in Section 4.3, planets in our grid with 𝑇 eff < 1400 K tend to have a chemical composition that is horizontally homogeneous. Thus, for these cases, any differences in the transit depth of the limb are caused by the temperature structure. Also for hotter cases, which can have horizontal chemical gradients, the chemical composition is not the deciding factor in setting the transit depth. One scenario is that chemical differences between both terminators remain small due to efficient horizontal quenching (see Fig. 14). In another scenario, the atmosphere is too hot for efficient horizontal mixing, and the chemical composition remains close to equilibrium, as would be the case for the hottest planets. In both cases, any difference in the limb compositions is directly attributed to temperature differences. Moreover, the dominant absorbing species, H 2 O and CO, remain approximately homogeneously distributed, even in hot atmospheres. It is not expected that three-dimensional advection will introduce an additional asymmetry in the equatorial regions of the morning and evening limbs (Mendonça et al. 2018b;Drummond et al. 2020).

The conclusion that thermal and not chemical limb differences are causing the transit depth asymmetry, has been formulated earlier by Agúndez et al. (2014a) for the hot Jupiters HD 209458 b and HD 189733 b. Here, we confirm and widen it to a large parameter space. That temperature variations drive the limb asymmetry is visible in the individual transmission spectra (Fig. 15). The amplitude of all molecular features is affected, rather than a select number of absorption bands related to a single species, which would be the case when chemical variations were the main cause of the spectral variation. Interestingly, the limb variation appears to be the smallest around 1 𝜇m, roughly the wavelength range of HST/WFC3. This additionally highlights the need for JWST in order to be able to study the exoplanet atmosphere limbs separately with ingressegress spectroscopy.

Implications for Benchmark JWST Targets

For a selection of planets, we sample the nearest atmospheric model from our grid of climate and chemistry models, and compare it with previous literature, theoretical models and observations. In the comparison with observations, mostly Hubble Space Telescope (HST) transmission spectroscopy, we apply manual vertical shifts in the transit depth, to correct for the fixed planetary radius in our models. We choose to focus on planets that have been heavily studied with 3D GCMs before, and/or will be observed with the James Webb Space Telescope (JWST). Note that we do not aim to provide quantitatively accurate atmospheric models for these specific exoplanets; this would be impossible given the wide grid spacing and model assumptions (see also Section 6). Rather, we aim to put our multi-dimensional modelling approach in a broader context by connecting the generic grid models back to the known diversity of exoplanets.

HD 209458 b

The hot Jupiter HD 209458 b, with its mass of 0.73 𝑀 Jup and radius of 1.39 𝑅 Jup (Stassun et al. 2017), has a gravity of 9.4 m/s 2 . It further has an equilibrium temperature of ∼1450 K (Barstow et al. 2017) and a rotation period of 3.5 days. Going by these measurements, the grid model it corresponds most to, has 𝑇 eff = 1400 K, 𝑔 = 10 m/s 2 and rotation period 2.32 days/host star type G5 (see Table 1). 4 Diagnostics for the temperature, chemical composition and spectrum for this model are shown in Fig. 16. The synthetic transmission spectrum in this figure has been computed assuming a stellar radius of 1.19 M , and a manual vertical shift has been applied to match the data of Sing et al. (2016). This is necessary because all grid models have been computed with the same radius (1.35 R Jup ), so the synthetic transmission spectrum differs from the data by a constant value.

The transmission spectrum of our HD 209458 b-like model agrees reasonably well with the observational data (Fig. 16,panel c). The Spitzer data points at 3.6 𝜇m and 4.5 𝜇m, and the water molecular feature at 1.4 𝜇m are reproduced, although the amplitude of the latter is a bit overestimated, probably because we did not include clouds (Deming et al. 2013). There is no good match between our model spectrum and the HST/STIS data in the optical wavelength range. The optical scattering slopes are determined by a variety of haze particles with non-grey opacities, as well as broad-base sodium and potassium absorption lines (Sing et al. 2016), both of which are not included in our opacities (see Section 2.4).

Regarding the chemical composition (Fig. 16, panel b), our results for HD 209458 b show some quantitative differences with the pseudo-2D model obtained with the same chemical kinetics code by Agúndez et al. (2014a). The upper atmosphere (𝑝 < 10 -3 bar) modelled by Agúndez et al. (2014a) shows distinct longitudinal gradients in many species, because of photochemistry. This effect is not taken into account in this work, and we plan to include it in a future study. Furthermore, we do not find the same strong vertical gradients in the molar fractions of CH 4 , NH 3 and HCN. Instead, we observe an efficient vertical mixing with almost constant molecular abundances. The discrepancy can be explained by the different temperatures employed. The atmospheric temperatures used by Agúndez et al. (2014a) have been derived from a 3D GCM with a temperature inversion (Showman et al. 2009), resulting in a much hotter day side, closer to equilibrium. This results in a strong decrease in CH 4 , NH 3 and HCN. It is now generally accepted that 2012)-network used by Agúndez et al. (2014a). We note that, for HD 209458 b, the variations due to the chemical network are expected to be small (Venot et al. 2020a). Despite the methodology differences, we obtain similar abundances for the dominant molecules, H 2 O and CO, as Agúndez et al. (2014a), and we arrive at the same conclusion that horizontal advection efficiently mixes the atmosphere toward the day-side composition.

The 3D chemistry-coupled GCM simulation of HD 209458 b by Drummond et al. (2020) predicts horizontally homogeneous chemistry, similar to what we find, although there is a difference in the vertical quenching pressures. Molar fractions for CH 4 , NH 3 and HCN become vertically constant around 100-200 mbar in the model of Drummond et al. (2020), whereas this happens deeper in the atmosphere in our work (a few bar). Again, the different chemical network employed by these authors, namely that of Venot et al. (2019), and without C 2 H 2 , could result in slightly different quenching pressures compared to this work, although this effect is expected to be minor for HD 209458 b (Venot et al. 2020a).

Rather, when comparing the vertical advection time-scales (Fig. 9 2020) (1600 K-1800 K), resulting in a shorter scale height and thus potentially more efficient vertical advection of chemical species. Note that the 1 bar-temperatures do not depend on the state of the deep atmosphere (Drummond et al. 2020), and that uncertainties of ∼100 K are typical in hot Jupiter GCMs.

Another possible reason for the comparatively high quenching pressure in our work compared to Drummond et al. (2020) could be that the vertical mixing efficiencies are in fact similar, but zonal mixing is more efficient in the work of Drummond et al. (2020), resulting in an efficient homogenisation of the chemistry to the day-side composition between 0.1 bar and 1 bar. Indeed, the zonal wind speed in our model (2.6 km/s) is somewhat lower than in Drummond et al. (2020) (> 4 km/s), and the vertical and zonal advection time-scales are comparable (Fig. 9).

It appears that temperature and wind speed variations between models -even if they are small -can influence the quenching pressures so as to cause order-of-magnitude quantitative changes in the atmospheric chemistry. Other than the quantitative differences associated with the vertical quenching, we find a good qualitative agreement with the more sophisticated 3D chemical kinetics model of Drummond et al. (2020). A planet such as HD 209458 b seems to be rather sensitive to variations between models, because the chemical and advection time-scales (Fig. 8), as well as the vertical and horizontal advection time-scales (Fig. 9 and 10) are similar. A hierarchy of modelling approaches, in tandem with benchmarking of GCMs and chemistry codes, should be used to alleviate this issue.

HD 189733 b

The irradiated exoplanet HD 189733 b has a mass of 1.13 𝑀 Jup and a radius of 1.13 𝑅 Jup (Stassun et al. 2017), resulting in a surface gravity of 21.9 m/s 2 . Furthermore, it is slightly cooler than HD 209458 b, with an equilibrium temperature of ∼1200 K (Barstow et al. 2017), and a rotation period of 2.2 days. Our best matching grid model has 𝑇 eff = 1200 K, 𝑔 = 10 m/s 2 and rotation period 1.26 days/host star type K5.

In analogy with HD 209458 b, after applying a vertical offset to our synthetic transmission spectrum, the observational data from the HST and Spitzer (Sing et al. 2016) are matched, except wherever clouds or hazes are concerned (Fig. 17,panel c). This includes the very prominent slope at optical wavelengths, which has been attributed to Rayleigh scattering on high altitude haze (Pont et al. 2013;Sing et al. 2016;Steinrueck et al. 2020). Furthermore, since our best-match model has a lower gravity than HD 189733 b, the amplitudes of the synthetic transmission features are too high for the observational data, even without taking into consideration that HD 189733 b is very cloudy.

We expect HD 189733 b to be chemically homogeneous (Fig. 17, panel b), as we have established for planets with 𝑇 eff < 1400 K. This is in agreement with Agúndez et al. (2014a), except for the rising HCN they find with altitude, which is a product of photochemistry on the day side. A notable distinction with their work, are their homogenized values for methane (10 -6 ) and ammonia (5 • 10 -6 ), which are higher in our model (10 -4 and 5 • 10 -5 respectively). We attribute this difference to the climate model: whereas the day side and evening limb in Agúndez et al. (2014a) are still CO-dominated in equilibrium, our model appears to be cooler at the quenching pressure (∼ 1 bar), allowing for a considerable equilibrium fraction of methane at all wavelengths. This strong dependence of the overall chemical composition on the deeper atmospheric temperature is indicative of a potential source of uncertainty in atmospheric models. This issue is further discussed in Section 5.4.1 for CH 4 and Section 5.5 in general.

We further find a good agreement with the work by Drummond et al. ( 2020) for most species. They report a deep meridional quenching of methane, which increases the methane content before it is vertically mixed throughout. This mechanism is a 3D-effect that is not present in our work. Nevertheless, we obtain a similar methane fraction, since the chemical equilibrium in our deep atmosphere already has a sizeable methane fraction before quenching. We conclude that the temperatures between 1 and 10 bar play an important role in the chemical composition. Although the temperature discrepancies between our work, and those of Agúndez et al. (2014a) and Drummond et al. (2020), can likely be attributed to physical model parameters, such as the low gravity of 10 m/s 2 , care must be taken so that the temperatures in these layers are well converged (Sainsbury-Martinez et al. 2019;Drummond et al. 2020), and unaffected by other physical processes that raise or lower the internal temperature (Fortney et al. 2020). We elaborate on this point further in Section 5.4.1.

WASP-43 b

WASP-43 b is a dense and fast rotating hot Jupiter with a mass of 1.998 𝑀 Jup and radius 1.006 𝑅 Jup (Esposito et al. 2017). This results in a gravity of 48.9 m/s 2 . It has a very similar temperature to HD 209458 b, namely ∼1400 K (Esposito et al. 2017). Finally, it has a short rotation period of only 0.81 days. We match it to a grid model with 𝑇 eff = 1400 K, 𝑔 = 100 m/s 2 and rotation period 0.79 days/host star type K5.

The synthetic spectrum of our WASP-43 b-like simulation (Fig. 18,panel c) shows a reasonable match with the observational transmission spectrum (Kreidberg et al. 2014), although the high surface gravity in our best-match model (100 m/s 2 ) results in a lower scale height, and thus, low-amplitude molecular features. Nevertheless, the correspondence between the data and our forward model is acceptable within the error bars, which is aided by the fact that this planet is expected to have a solar composition and is cloud-free within the observable pressure range (Kreidberg et al. 2014).

Because of its comparatively fast rotation, WASP-43 b has a limited heat redistribution, resulting in a relatively cold night side and morning terminator (Fig. 18, panel a, and Carone et al. (2020)). This large day-night contrast results in some chemical gradients around 100 mbar, especially of methane. In general, we find a very good agreement between the chemical composition for a WASP-43 b-like planet (Fig. 18) and another pseudo-2D chemical kinetics model for this planet (Venot et al. 2020b). In the pressure range where the two models overlap, the only big difference is that the horizontal transport seems to be more efficient in the work by Venot et al. (2020b). This is explained by the horizontal wind speed they adopt: 4.6 km/s, compared to 2.9 km/s in this work. This substantial difference stems from different modelling assumptions in the 3D GCM simulation (see Carone et al. 2020, for a detailed comparison). An additional reason for discrepancies could be the different chemical networks used, namely the reduced Venot et al. (2020a)network used in this work and the Moses et al. (2013)-network used by Venot et al. (2020b).

In comparison with a 3D chemical advection model with simplified chemistry (Mendonça et al. 2018b), we find good agreement for H 2 O, CO and CO 2 , and somewhat less for CH 4 . A possible reason could be that the night-side temperatures are higher in the study of Mendonça et al. (2018b), resulting in a methane abundance that is already horizontally uniform, and lower than seen in Fig. 18, panel b. No strong meridional advection is expected for this planet (Mendonça et al. 2018b).

WASP-121 b

The ultra-hot Jupiter WASP-121 b has a mass of 1.183 𝑀 Jup and a radius of 1.865 𝑅 Jup (Delrez et al. 2016), and hence a surface gravity of 8.4 m/s 2 . It is orbiting an F-type star on a comparatively short orbit of 1.3 days, giving rise to its high equilibrium temperature of ∼2400 K. In our grid, it is best matched by a model with 𝑇 eff = 2400 K, 𝑔 = 10 m/s 2 and rotation period 1.13 days/host star type F5.

When the spectrum of our WASP-121 b-like atmosphere (Fig. 19,panel c) is compared to observations (Evans et al. 2016), a discrepancy is seen in both the shape of the water band and the optical wavelengths, despite our good match to the physical temperature and gravity of WASP-121 b. The excess absorption has been interpreted as TiO-and VO-opacity, with potential contributions from FeH (Evans et al. 2016) or a partial cloud cover at the morning terminator (Parmentier et al. 2018). These opacities -typical for ultra-hot planets -are not included in our homogeneous modelling approach, and we discuss the applicability of our grid to ultra-hot atmospheres in Section 6. 2016). The transmission spectrum is cloud-free, and does not contain Na/K opacities. In the chemistry plot, the equilibrium compositions of the morning and evening limb are plotted as thin dashed lines. Kreidberg et al. (2014). The transmission spectrum is cloud-free, and does not contain Na/K opacities. In the chemistry plot, the equilibrium compositions of the morning and evening limb are plotted as thin dashed lines. 2016). The transmission spectrum is cloud-free, and does not contain Na/K opacities. In the chemistry plot, the equilibrium compositions of the morning and evening limb are plotted as thin dashed lines.

Our atmospheric models for a WASP-121 b-like planet display a hot atmosphere, with a day-side temperature of about 2000 K, a poor heat redistribution, and cold night sides (Fig. 19). This results in very low abundances of CH 4 , NH 3 and HCN on the day side, and somewhat increased fractions at the night side. Horizontal quenching is, however, sufficiently quick, and these molecules never attain number fractions over 10 -6 . Since we did not include TiO, VO, or ionized-metal opacities in our models, we do not retrieve a strong day-side thermal inversion, and thus, the day-side temperatures in our model are too low to cause water dissociation, which occurs typically at temperatures above 2500 K (Lothringer et al. 2018;Parmentier et al. 2018). Hence, we find a constant water abundance throughout the atmosphere. Previous models mapping the 2D or 3D chemical composition of ultra-hot Jupiters have either invoked chemical equilibrium (Parmentier et al. 2018) or only vertical mixing (Molaverdikhani et al. 2020), so the question whether water can exhibit horizontal gradients, or limb asymmetries, in ultra-hot atmospheres, remains to be answered. We note that in our hottest models with 𝑇 eff = 2600 K and inefficient heat recirculation (M5-star host), the day-side does reach temperatures high enough for H 2 O dissociation. This is also visible in the chemistry models (Fig. 8). In these cases, it appears that horizontal circulation is not efficient enough to completely mix the atmosphere to the water-depleted, day-side state. It is important to understand the modifying role of dynamical disequilibrium chemistry in molecular dissociation, as it has been shown that 3D effects near the limbs can bias retrievals of ultra-hot Jupiters (Pluriel et al. 2020).

WASP-107 b

The warm sub-Saturn mass planet WASP-107 b has a mass of 0.096 𝑀 Jup and radius of 0.96 𝑅 Jup (Anderson et al. 2017;Piaulet et al. 2020). This results in a very low gravity of 2.6 m/s 2 . Furthermore, its equilibrium temperature is about ∼750 K, with an orbital period of 5.7 days. Hence, we opt to match it to a grid model with 𝑇 eff = 800 K, 𝑔 = 1 m/s 2 and rotation period 4.25 days/host star type K5. Its low gravity and effective temperature make it a suitable case study for our low-gravity models with disequilibrium chemistry.

The synthetic transmission spectrum of our WASP-107 b-like model (Fig. 20,panel c) does not correspond well to the HST transmission spectrum presented by Spake et al. (2018). The main reason for the discrepancy is a thick, opaque cloud or haze layer at altitudes of at most a few mbar (Kreidberg et al. 2018). This severely mutes the water absorption, which would normally be very extended because of the low gravity and high scale height of this planet.

The chemical composition is homogeneous, in accordance with the dichotomy established in Section 4.3. Additionally, the climate does not exhibit strong temperature variations (Fig. 20). Because of vertical mixing, the atmospheric chemistry is CO-dominated under our assumption of solar metallicity and C/O ratio, although a considerable methane fraction is still present. The transmission spectrum still shows very prominent methane absorption.

Interestingly, retrieval modelling by Kreidberg et al. ( 2018) concluded a depletion of methane for this planet. This could indicate a composition that is far removed from solar values, and/or a combination of a high internal temperature with deep vertical mixing (Kreidberg et al. 2018;Fortney et al. 2020). Our calculation of 𝐾 𝑧𝑧 yields 10 8 -10 9 cm 2 /s for this WASP-107 b-like model, which 2018). The transmission spectrum is cloud-free, and does not contain Na/K opacities. In the chemistry plot, the equilibrium compositions of the morning and evening limb are plotted as thin dashed lines. should be high enough to quench methane at low abundances, as demonstrated for the comparable planet WASP-117 b (Carone et al. 2021). However, a low intrinsic temperature of 200 K has been adopted for all models in this work. Hence, we do not observe the reducing behaviour of deep quenching very strongly here.

Notable Detection Windows

Observing methane

Methane is an important molecule to observe in exoplanetary atmospheres, since its abundance -together with that of water -can provide a handle on the C/O ratio, which can be linked to planet formation processes (Öberg et al. 2011;[START_REF] Notsu | [END_REF]. For this reason, several attempts have been made to constrain the methane abundance in exoplanet atmospheres, using Hubble Space Telescope WFC3 transmission spectroscopy. However, many of these studies have resulted in an absence of methane absorption, even in the atmospheres of cool planets with equilibrium temperatures below 1000 K (Kreidberg et al. 2018;Benneke et al. 2019;Carone et al. 2021). In one instance, transmission spectroscopy has provided significant evidence for both water and methane absorption (Chachan et al. 2019). Additionally, in a different planet, with an estimated equilibrium temperature of at most 850 K, a methane detection has been reported using high-resolution spectroscopy (Guilluy et al. 2019), although a re-analysis of the data with a different line list has not resulted in a detection (Gandhi et al. 2020).

The difficulty of observing methane with transmission spectroscopy could be in part due to the largely overlapping absorption bands of water and methane in the HST/WFC3 wavelength range (see Bézard et al. 2020, for a discussion). Moreover, since methane formation is more favourable in relatively cold environments, future methane observations with transmission spectroscopy will likely be plagued by cloud obscuration (Crossfield & Kreidberg 2017;Gao et al. 2020;Taylor et al. 2020). Nevertheless, it is instructive to discuss the detectable presence of methane in the cloudless grid presented in this work, taking into account disequilibrium chemistry through horizontal and vertical mixing.

In the medium-gravity atmospheric models of our grid, methane absorption bands become prominent at effective temperatures of 1200 K or lower. Methane is not necessarily the dominant carbon-bearing molecule in these models, but the CH 4 volume mixing ratio approaches or exceeds 10 -4 (Fig. 8), which should indeed be above the detection threshold (e.g. Fortney et al. 2020). All rotation rates at 1200 K show a high methane abundance with prominent absorption bands (Fig. 11), but the fast rotating models more so, because the quench point is situated at slightly lower pressures. In this temperature regime, meridional quenching, not considered in this work, could result in an additional increase in the methane abundance (Drummond et al. 2020). Even at 1400 K, methane has an opacity comparable to that of water in some distinct methane absorption bands at 2.3 𝜇m, 3.3 𝜇m and 7.8 𝜇m. The former of those, however, is located between two water absorption bands, and may thus be difficult to distinguish from clouds. Hence, the latter appear to be more favourable for methane detections. Furthermore, we note that the distinct methane absorption bands diminish at effective temperatures of 600 K or lower (Fig. 11). At these low temperatures, the spectrum consists of numerous overlapping absorption bands of H 2 O, CH 4 and NH 3 .

When we compare the transmission opacities of methane with our low-and high-gravity atmospheric models, we find that the low-gravity models display the same general trend as above, with methane becoming the dominant opacity source at temperatures equal or below 1200 K. At 1400 K, however, no methane absorption peaks at 2.3 𝜇m, 3.3 𝜇m or 7.8 𝜇m are present. The high-gravity models, on the other hand, can show methane opacities comparable to those of water for relatively hot atmospheres even. The main cause is the comparably low heat redistribution (Fig. 5), and thus cool temperatures near the limbs (Fig. E2), which can give rise to a strong longitudinal dependence of the chemical composition, and methane formation on the night side and near the limbs (Fig. E10). Thus, for these very high-gravity cases, methane could be detectable at effective temperatures of 1400 K or hotter, depending on the rotation rate and horizontal mixing efficiency of the planet. However, it is unlikely that the effect of an increased temperature on the atmospheric scale height is enough to offset the adverse effect of gravity. Thus, our pseudo-2D, cloud-free chemical kinetics models indicate that it may be worthwhile to investigate planets with effective temperatures between 600 K and 1200 K for the presence or absence of methane absorption. The 3.3 𝜇m absorption band, which will become accessible through JWST, appears to be the most suitable. For some hot Jupiters with high densities (𝑔 > 10 m/s 2 ), it may even be valuable to employ transmission spectroscopy to hunt for methane originating from the night side. If night-side clouds and photochemical hazes are not too dense, methane could exhibit a relatively strong signal. We find that the rotation rate has some influence on the detectability of methane in the atmosphere, by raising or lowering the vertical mixing quench pressure. In this context, rapidly rotating planets are more prone to contain appreciable methane concentrations.

We caution, however, that the CH 4 abundance is strongly impacted by the vertical mixing efficiency, which is not straightforward to determine accurately. Moreover, the internal temperature of the planet plays an important role in this analysis. In this work, an intrinsic temperature of 200 K is adopted for all planets, which is likely too low for inflated exoplanets, by several 100 K (Thorngren et al. 2019;Sarkis et al. 2020;Fortney et al. 2020). The deep-atmosphere temperatures determine the vertical chemical equilibrium gradient. If a species' equilibrium abundance increases with altitude, a lower eddy diffusion will result in a higher disequilibrium abundance, and vice-versa. Furthermore, a steeper gradient will result in a higher sensitivity to the eddy diffusion coefficient. Indeed, it has been shown that atmospheric chemistry can be altered by the internal temperature (Agúndez et al. 2014b;Fortney et al. 2020). On top of that, the convergence of deep atmospheric layers is a well-known problem in GCM simulations (Amundsen et al. 2016;Mayne et al. 2017;Carone et al. 2020;Wang & Wordsworth 2020;Showman et al. 2020). As input temperatures for chemical disequilibrium studies are often derived from GCM simulations (e.g. Agúndez et al. 2014a;Venot et al. 2020b;Helling et al. 2019, and this work), care must be taken that the deep (𝑝 > 1 bar) atmospheric temperatures are well-converged, or do not feed back into the chemistry through deep vertical mixing. In this work, we have applied the recommendation laid out by Sainsbury-Martinez et al. (2019) to ensure that a suitable temperature profile is attained in the deep atmosphere (see also the discussion in Drummond et al. 2020). We conclude that a detailed study of the planetary interior, and its impact on the chemistry in a 3D context, would be beneficial. Furthermore, JWST observations of transiting exoplanets with effective temperatures below ∼1200 K will be instrumental in constraining the processes of methane quenching and interior heating, through vertical mixing of the deep atmosphere, or horizontal transport from the night side.

Bump in CO 2 opacity

Carbon dioxide is an important molecule in exoplanet atmospheres, because of its impact on habitability of rocky planets, as a greenhouse gas (e.g. Schwieterman et al. 2018). In irradiated gaseous exoplanets, it can be used as a probe for the atmospheric metallicity and C/O ratio (Zahnle et al. 2009;Moses et al. 2013;Fleury et al. 2020). Furthermore, unlike methane, it is fairly unaffected by chemical disequilibrium (Fig. 8).

Although carbon dioxide only has a few absorption bands, the opacity in these bands can be relatively high, so that CO 2 dominates the opacity at 4.3 𝜇m, even at low mixing ratios of ∼10 -7 . This is illustrated for a medium-gravity model with 𝑇 eff = 1400 K in Fig. 21. The opacity bump at 4.3 𝜇m by CO 2 , and to a lesser extent CO, is also very distinctive in hotter atmospheres (Fig. 11). When assuming a G5-host star radius, we estimate the transit depth difference associated with this CO 2 transmission feature to be about 65 ppm or 1 ppm, for our low (𝑔 = 1 m/s 2 ) and medium (𝑔 = 10 m/s 2 ) gravity models respectively, making it an intriguing, albeit challenging, observing window for this molecule.

Outlook of observing ammonia

Ammonia (NH 3 ) has an interconversion relation with nitrogen gas (N 2 ), so that, in chemical equilibrium, the former is expected to occur in colder, and the latter in hotter exoplanet atmospheres (e.g. Venot et al. 2012;Moses 2014). Together with the photochemically produced HCN, these molecules are expected to be among the main nitrogen-bearing molecules in exoplanet atmospheres. Up to now, detections of nitrogen-bearing molecules in transiting exoplanet atmospheres have been scarce and tentative (Tsiaras et al. 2016;MacDonald & Madhusudhan 2017). Still, there is much potential in nitrogen-detections, since these molecules should be good indicators of disequilibrium chemistry (both through vertical mixing and through photochemistry). Furthermore, information on the atmospheric C/N or N/O ratio would provide valuable constraints for planet formation theories (e.g. Mordasini et al. 2016;Piso et al. 2016;Cridland et al. 2020). Some of the spectra in the cool parameter space of our grid feature an NH 3 opacity peak in the mid-infrared, at 10𝜇m (Fig. 22), making it a potential target for spectroscopy with JWST/MIRI. Because the ammonia peak at 10 𝜇m is quite distinctive, and located in a region with only pseudo-continuous water opacity, it is visible in a large range of synthetic spectra, from 1400 K to 400 K (Fig. 11). In a behaviour analogous to methane (Section 5.4.1), the relative amplitude of the 10 𝜇m ammonia feature appears to increase with gravity. However, the absolute amplitude of the feature, and thus its detectability, decreases with gravity, so medium-to low-gravity planets would still be more suitable to detect ammonia. For the model shown in Fig. 22, we estimate a transit depth difference of 1 ppm between the nominal spectrum and one without NH 3 . For the low-gravity case, the feature amplitude is estimated at 106 ppm. Finally, we did not find a strong dependence of ammonia absorption on the rotation period.

Impact of Deep Atmospheric Quenching

In the case of deep quenching (𝑝 > 1 bar), both the thermal structure of the atmosphere at depth and the vertical mixing strength at depth are critical in setting the disequilibrium composition higher up in the atmosphere. This is highlighted through our comparison of HD 209458 b (Fig. 16) and HD 189733 b (Fig. 17) with earlier simulations using a different thermal structure and 𝐾 𝑧𝑧 -profile (Agúndez et al. 2014a). As pointed out in Section 5.3, discrepancies in these cases are often caused by a different temperature structure at depth, resulting in a different chemical equilibrium composition at the vertical quenching pressure. As such, any uncertainties in the deep thermal profile could be propagated to the chemical composition in general.

If the thermal evolution in the deep part of the atmosphere is adiabatic, the associated vertical temperature gradient tends to be steep, and a small shift in the quenching pressure could result in a big change to the chemistry. Moreover, if the quenching happens near a chemical transition, like in our HD 189733 b simulation (Fig. 17), the impact on the chemical composition is exacerbated. Therefore, it appears that great care should be taken in the selection of both the temperature profile, and the vertical mixing strength at depth, if the quenching pressure is high. We now discuss the determination of these two factors.

Vertical Mixing in the Deep Atmosphere

Giant irradiated exoplanets are expected to consist of a stably stratified radiative atmosphere, overlying a deep convective atmosphere (e.g. Marley & Robinson 2015), with a radiative-convective boundary situated roughly between 1 bar and a few 100 bar, depending on the planet's gravity and incident irradiation (Thorngren et al. 2019). In the radiative part, which includes the photosphere, vertical mixing is expected to be dominated by atmospheric circulation or atmospheric waves, which can be approximated as a 1D diffusion process with coefficient 𝐾 𝑧𝑧 (Parmentier et al. 2013;Charnay et al. 2015;Komacek et al. 2019). The eddy diffusion coefficient generally increases with altitude (Fig. 7), and hence reaches its lowest values when approaching the deep atmosphere (𝑝 1 bar). Convection in the deep atmosphere, however, ensures efficient mixing, and the value of 𝐾 𝑧𝑧 can be expected to increase again as the vertical mixing becomes dominated by small-scale turbulent convection, which is not captured by GCMs.

In the deep convective region, mixing length theory can be used to estimate the eddy diffusivity (𝐾 𝑧𝑧 ∼ 𝑤𝐿) (e.g. Wang et al. 2015;Zhang 2020, and references therein), resulting in values around 10 9 cm 2 /s for Jovian planets. However, knowledge of the vertical velocities and/or the thermal structure, in addition to a good estimate of the mixing length (Smith 1998;Bordwell et al. 2018) are still required to determine 𝐾 𝑧𝑧 in the deep regions of giant exoplanets.

It should thus be clear that is not trivial to determine a vertical mixing profile, especially in the 3D context of an irradiated exoplanet. The pressure of the radiative-convective boundary is location-dependent (Rauscher & Showman 2014, and our Fig. B1), which makes it more likely that on the cold night side quenching pressures fall within the convective zone. Compressible hydrodynamic simulations suggest that convective overshooting is limited (Freytag et al. 2010), so 𝐾 𝑧𝑧 could change quickly close to the radiative-convective boundary, as some previous chemical kinetics studies have assumed (e.g. Moses et al. 2011;Venot et al. 2020b). Furthermore, convectively driven, vertically propagating gravity waves could provide an additional source of vertical mixing above the radiative-convective boundary (Medvedev & Yiğit 2019;Bordwell 2020), but resolutions of typical exoplanet GCMs are too low to capture these waves. The potentially rapid increase in vertical mixing efficiency in the deep interior is not taken into account in our 𝐾 𝑧𝑧 -profiles, but could contribute to the prevalence of deep quenching, and the associated chemical uncertainties. Future hydrodynamic studies focusing on the transition between the deep convective zone and the stratified radiative layers above, would provide valuable constraints for exoplanet disequilibrium chemistry.

Temperatures in the Deep Atmosphere

The assumed temperature profile in the deep part of the atmosphere is crucial to the chemistry, not only locally, but throughout the vertical extent of the atmosphere if quenched. However, unlike the 3D temperature structure at lower pressures, which can be relatively well constrained with a GCM, the deep thermal structure is often ill-constrained by GCMs, because it evolves very slowly (Carone et al. 2020;Mendonça 2020;Wang & Wordsworth 2020;Showman et al. 2020). This can result in temperature-profiles that are only partially converged, featuring a small thermal inversion in the deep (𝑝 ≈ 10 bar) layers (e.g. Agúndez et al. 2014a;Amundsen et al. 2016;Venot et al. 2020b;Drummond et al. 2020). Although it can be argued in some cases that the observable atmosphere is not affected by the deep thermal evolution, this does not appear to be true in general (Carone et al. 2020). Additionally, we stress that care should be taken when unconverged temperature profiles are used as input for disequilibrium chemistry studies, because of the potential for deep quenching.

Simulations by Sainsbury-Martinez et al. (2019) suggest that GCMs tend to converge to a hot, deep adiabatic temperature profile when evolved up to a full steady-state. Our temperature plots indeed show that an adiabatic gradient is reached at depth (Fig. 4). Although this fact is reassuring, there is no guarantee that this adiabat is physically suitable. Fortney et al. (2020) have demonstrated the sensitivity of the chemical disequilibrium abundances to the interior temperature. Irradiation (Thorngren et al. 2019), tidal heating (Agúndez et al. 2014b;Fortney et al. 2020), radius inflation (Tremblin et al. 2017), Ohmic heating (Batygin et al. 2011), or the planet age (Fortney et al. 2020) all impact the interior temperature, and thus potentially also the chemical abundances if quenching occurs deeply.

It this context, a self-consistent model coupling the atmospheric circulation with the deep interior evolution of giant exoplanets in a 3D framework would be very beneficial. Such a coupling is needed in order to study the interface near the radiative-convective boundary, and gauge the impact of the interior temperature on the (observable) atmospheric layers above. Finally, we note that the deep interior is not totally inaccessible to observations. Indeed, the detection of certain refractory species in the gas-phase at high altitudes could place valuable constraints on the physical conditions below (Sing et al. 2019).

Implications for Brown Dwarf -White Dwarf Binaries

By systematically applying our grid parameter combinations of effective temperature and stellar type, under the assumption of synchronous rotation, we have produced numerous atmospheric models with rotation rates below 10 hours (residing in the 'lower left' corner of the grid space, Table 1, Fig. 2 andfollowing). These increasingly short rotation periods are not typical for giant exoplanets, of which the fastest rotators have orbital periods of ∼1 day, whereas these models have rotation periods of the order of hours. However, because of their high equilibrium temperatures and short rotation periods, these ultra-fast rotating models have interesting real-life analogues, namely brown dwarf -white dwarf (BD-WD) binaries, a rare class of objects with similarities to both brown dwarfs and ultra-hot Jupiters (Lothringer & Casewell 2020). This is especially the case for our high gravity set of models, with 𝑔 = 100 m/s 2 . In this section, we highlight the applicability of our models to BD-WD binaries.

Two important distinctions between our modelling approach, stemming from a self-consistent extrapolation of hot Jupiter modelling, and BD-WD binaries are:

• the high intrinsic temperatures in irradiated brown dwarfs;

• the ionizing (UV-)radiation from the white dwarf host.

Indeed, whereas we have adopted an intrinsic temperature of 200 K in our models, values of up to 2000 K have been retrieved for irradiated brown dwarfs (Lothringer & Casewell 2020). This could impact both the circulation (Carone et al. 2020), as well as the chemistry (Fortney et al. 2020) when the vertical quench level is sufficiently deep. Furthermore, as explained in Section 6, the incorporation of longitude-dependent photochemistry could be important for these objects in the future.

Regarding the atmospheric circulation, we compare our results with recent previous studies of BD-WD binary climates (Lee et al. 2020;Tan & Showman 2020). In general, we find qualitatively agreeing circulation regimes with a large day-night temperature contrast, a very thin, eastward zonal-mean equatorial jet stream, and multiple smaller jet structures at high latitudes. Similarly to Lee et al. (2020), we obtain retrograde flow at the equatorial day-side, resulting in high local wind speeds near the morning terminator, which is typical for fast rotation and high gravity (Carone et al. 2020).

From our chemistry modelling, it appears that BD-WD binaries should feature a chemical regime that is only partially zonally quenched. A certain contamination from the day side can be expected, but our models suggest that advection by the zonal mean circulation is not efficient enough to remove all horizontal chemical gradients. Depending on local thermal conditions, these objects could exhibit zonally changing water absorption, if horizontal quenching of thermally dissociated H 2 O is inefficient. Thus, besides the temperature itself, it may be needed to take into account the local chemical composition to interpret phase curve observations (Casewell et al. 2015). We stress that our models do not include a stratosphere or photochemistry, both of which could be present in BD-WD binaries (Lothringer & Casewell 2020), and could impact the circulation and chemical composition, as well as their zonal dependence.

METHODOLOGY LIMITATIONS AND FUTURE WORK

Post-Processing Models

The connection between the radiative transfer calculations with pe-titCODE and the 3D climate model MITgcm is made with the socalled Newtonian cooling (also called Newtonian relaxation) framework. It comprises a parametrized atmospheric forcing, lacking selfconsistent feedback of dynamical processes to the radiative heating and cooling terms. This simplified radiation treatment is ideal for dynamical studies, but less suitable for a direct comparison with observational fluxes (Showman et al. 2020). Comparisons between GCMs coupled with Newtonian cooling or more robust radiative schemes yield qualitatively similar outcomes, but with quantitative differences. The day-night temperature contrast was found to be reduced (compare Showman et al. (2009); Amundsen et al. (2016) with Showman et al. (2008); Mayne et al. (2014)), especially at high altitudes, where radiative time-scales are short. Furthermore, the lack of a self-consistent radiation can sometimes result in the appearance of a weak, dynamically induced temperature inversion. This can be attributed to an insufficient coupling between radiation and dynamics at the hot day side of irradiated exoplanets (Carone et al. 2020). Indeed, such feature disappears when the output of a GCM with Newtonian relaxation (Showman et al. 2008) and with self-consistently coupled radiation (Showman et al. 2009) are compared. Nevertheless, Newtonian cooling has been used widely for parameter and ensemble studies of exoplanet climates, as well as detailed dynamical studies that benefit from a parametrized radiative time-scale, since it is faster and more readily interpreted than a full radiation-coupled climate model (see e.g. Showman et al. 2008;Rauscher & Menou 2010;Mayne et al. 2014;Zhang & Showman 2018a,b;Sainsbury-Martinez et al. 2019;Carone et al. 2020). Moreover, in a comparison made between a suite of models with Newtonian relaxation (Komacek & Showman 2016;Zhang & Showman 2017) and with a more sophisticated double-grey radiation scheme (Komacek et al. 2019), it was noted that the results yielded by those simulations agree qualitatively, as well as quantitatively by up to a factor of ∼2 over a large parameter space. This motivates our choice for Newtonian relaxation as a mechanism to consistently investigate a large grid of planetary atmospheres.

Furthermore, employing a fixed background temperature for the chemical kinetics calculations neglects feedback which could have taken place in a coupled dynamical-chemical model. Two of such feedback mechanisms that could be important, come to mind. First, chemical heat sources or sinks could impact the temperature because of the endo-/exothermic nature of chemical reactions. Second, the presence of chemical elements can influence the temperature via the radiation field through radiative heating and cooling. The impact of this feedback mechanism has been studied in 1D for two hot Jupiters (Drummond et al. 2016) and a hot Neptune (Agúndez et al. 2014b). Both of these studies concluded that temperature differences of up to 100 K are possible for cases of strong chemical disequilibrium, in comparison with post-processed chemistry models. A similar magnitude in temperature variation was found by Steinrueck et al. (2019), who used parametrized disequilibrium chemistry to investigate the feedback to the radiation and temperature field in a 3D GCM. The uncertainties associated with the Newtonian-cooling GCM simulations of this work are about ∼100 K as well. Therefore, we do not expect the post-processing methodology that was used in this work to undermine our conclusions. Still, it would be beneficial if the impact of consistent chemistry-, dynamics-and radiation-coupling could be investigated in a broader parameter space in future studies.

Ultra-Hot Jupiters

There are many indications that ultra-hot Jupiters, highly irradiated exoplanets with day-side temperatures above ∼2200 K, are emerging as a distinct class of planets (Baxter et al. 2020), featuring a range of physical processes that are not included in this study. Indeed, ultrahot Jupiters are expected to have thermal inversions (stratospheres) on the day side, caused by atomic and ionized metal absorption, and/or TiO and VO absorption (Lothringer et al. 2018;Arcangeli et al. 2018;Parmentier et al. 2018). These opacities can give rise to temperature excesses of 1000 K. Thermal inversions, along with the detection of species that cause them, is still a highly active topic of research (e.g. Piette et al. 2020). Furthermore, molecular dissociation on the day side and recombination on the night side of these planets is expected to constitute a substantial part of the heat transport budget (Bell & Cowan 2018), thereby heating the night side of these planets and reducing the day-night temperature (Tan & Komacek 2019;Keating et al. 2019).

These physical effects, which are now in the process of being understood, make it difficult to interpret atmospheric models that span a large range of equilibrium temperatures. More concretely, by excluding thermal inversions and hydrogen recombination, our 3D climate models are expected to have both cooler day sides and cooler night sides than what would be considered realistic, based on the current known ultra-hot Jupiter population. As of yet, it is unclear whether the inclusion of strongly absorbing opacities and hydrogen recombination would increase or decrease the day-night temperature contrast compared to the models presented here. However, recent phase curve observations of ultra-hot Jupiters suggest that the overall effect would be a day-night contrast that is lower than obtained with the simplifications assumed in this work (Mansfield et al. 2020). Nevertheless, these simplifying assumptions were made in favour of having a more continuously varying suite of atmospheric models to investigate, and to build up a hierarchical framework from which more complex models can draw.

Pseudo-2D Framework

The pseudo-2D framework used in this work has the obvious disadvantage that dynamics is only taken into account in a parametrized or simplified way. First, in the vertical direction, the widely-used but ill-determined 𝐾 𝑧𝑧 parameter is used. By using scaling relations (Komacek et al. 2019), we have applied informed estimates of the eddy diffusion coefficient, based on the global vertical wind speed and zonal wind advection, over a large parameter space. However, the vertical mixing could have local efficiency extrema, such as vertical advection 'chimneys' (Parmentier et al. 2013;Zhang & Showman 2018b;Komacek et al. 2019;Carone et al. 2020;Steinrueck et al. 2020), which are not captured in this approach. Since we have established that the chemical composition, in some cases, depends strongly on an accurate determination of the quenching pressure, the ability to make precise predictions of the eddy diffusion coefficient over a wide range of exoplanets will remain important in the future.

Next, in the zonal direction, a uniform, unidirectional wind advection is assumed. This assumption is valid for most of the climates in this grid, as discussed in Section 4.2.1. However, the assumption of a uniform zonal wind restricts the area under scrutiny to the equatorial region, since it is expected to break down in the high-latitude regimes of synchronously rotating exoplanets. As observational data becomes more accurate, taking into account the full 3D geometry of exoplanet transits may become necessary (Caldas et al. 2019;Pluriel et al. 2020).

Last, the meridional direction is absent in our framework. Meridional chemical advection has been studied with a GCM coupled to simplified chemical relaxation (Mendonça et al. 2018b;Drummond et al. 2018a,b), as well as a fully coupled 3D chemical kinetics models (Drummond et al. 2020). We have discussed these research results in our case studies of HD 209458 b, HD 189733 b and WASP-43b (Section 5. 3.1, 5.3.2 and 5.3.3 respectively). In only one case -HD 189733 b -meridional advection has been found to significantly impact the chemical composition. Incidentally, this is the planet with the lowest equilibrium temperature. Thus, it remains an open question if this finding may be extrapolated, and if meridional quenching plays an important role in the chemistry of all cool exoplanet atmospheres.

Despite these simplifying assumptions, solving chemical kinetics in a pseudo-2D framework has the obvious advantage of being reasonably fast compared to 3D GCMs with fully coupled chemistry. This advantage, combined with its more parametrized nature, allows for broad parameter studies, such as this work and Moses et al. (2021). Potentially, detailed microphysical cloud modelling (e.g. Helling et al. 2019) could also benefit from a pseudo-2D framework to incorporate simplified atmospheric dynamics, since cloud formation and gas-phase chemistry are intrinsically coupled. Considering that more observations of phase curves and potential limb asymmetries will become available in the JWST and Ariel era, pseudo-2D chemical kinetics will remain a useful tool to interpret and investigate these spatially resolved data.

Photochemistry

Photochemistry is a potentially important chemical disequilibrium effect that is currently not taken into account in this work. Photochemistry can cause photolysis, molecular dissociation in the upper atmospheres of exoplanets and brown dwarf-white dwarf binaries, leading to the reduction of some molecules, such as CH 4 , and the increase of others, such as HCN (Moses et al. 2011;Venot et al. 2012). Chemical kinetics in combination with the pseudo-2D framework applied in this work is well equipped to take into account the changing irradiation angle with longitude (Agúndez et al. 2014a;Venot et al. 2020b). This is especially the case for the self-consistent variation of the host star type employed in the grid setup of this work. Recently, the need for such photochemical modelling has been demonstrated for the hot Jupiter HAT-P-41 b (Lewis et al. 2020). We thus aim to include photochemistry in a detailed followup study.

CONCLUSIONS

In this work, we have synthesized a large grid of atmospheric models for synchronously rotating, irradiated giant exoplanets. Our cloudfree models include the three-dimensional climate of the planet, together with two-dimensional information about the chemical composition of the atmosphere. From these models, we have synthesized transmission spectra to assess the potential impact on observational data. These simulations were conducted for temperatures ranging from 400 K to 2600 K, gravities between 1 m/s 2 and 100 m/s 2 , and for different rotation rates. We have included chemical disequilibrium effects in the form of vertical mixing, parametrized by the eddy diffusion coefficient 𝐾 𝑧𝑧 , and horizontal mixing, approximated as uniform eastward advection. Our aim was to explore the chemical diversity over a wide range of effective temperatures, gravities and rotation rates, with a specific focus on the role of vertical and horizontal mixing.

Our main finding is that a dichotomy is expected in the chemical compositions of exoplanet atmospheres. Cool planets with effective temperatures below 1400 K are expected to be chemically homogeneous, whereas hot planets with effective temperatures above 1400 K show longitudinal changes in their composition. The established homogeneity in cool planets is caused by deep vertical mixing, where species are quenched before chemical equilibrium leads to strong longitudinal gradients. Furthermore, through a timescale analysis we found that for most cases, the zonal advection is more efficient than vertical mixing at a pressure of 1 bar. In hotter planets with effective temperatures above 1400 K, zonal advection is not efficient enough to completely remove longitudinal chemical gradients, and colder regions, like the night side and morning terminator, tend to be chemically contaminated by air advected from the day side. These longitudinal gradients, mainly of CH 4 , NH 3 and HCN, are exaggerated with increasing gravity.

The rotation rate was found to impact the chemical composition in multiple ways. A high rotation rate leads to a reduced day-tonight heat redistribution, a slower horizontal wind advection, and lower vertical mixing efficiencies. This results in bigger horizontal variations of the chemical composition, and in some cases, a shift in the CO-CH 4 balance. In the latter case, the impact of the rotation rate on the transmission spectrum could be considerable, but in general, it is minor.

By adopting a hierarchical approach in composing the synthetic transmission spectrum, we established that the inclusion of a three-dimensional temperature, informed by a GCM, has the biggest impact on the chemical composition at the exoplanet limbs, and on its transmission signature. Chemical disequilibrium effects are of secondary importance. For cool planets (𝑇 eff < 1400 K), horizontal advection appears to be unimportant if vertical quenching occurs deep enough. For hotter planets, on the other hands, horizontal advection plays an important homogenizing role, as the chemical composition is set close to that expected on the day side of the planet.

Differences in the atmospheric structure at the morning and the evening terminator lead to different contributions of both limbs to the exoplanet transit spectrum. For a wide parameter space, we confirm the conclusion of Agúndez et al. (2014a) that the temperature, rather than the chemical composition, is the main cause for transit depth asymmetries between the morning and evening limbs. Furthermore, the maximal transit depth difference between the two limbs was found to increase with the effective temperature of the atmosphere. We estimate that, for a cloud-free planet, featuring a low gravity, a hot atmosphere and a slow rotation rate, limb differences could be detectable with James Webb Space Telescope ingress and egress transit observations. Finally, we note some of the limitations of our modelling setup, which could be improved upon in future studies. It is likely that the day-night temperature contrast in our hottest planets (𝑇 eff > 2000 K) is overestimated, since we did not include hydrogen dissociation and recombination feedback in our climate models. GCMs with detailed radiative transfer, including some high-opacity species omitted in this work, would be needed to more accurately constrain the thermal structure for these ultra-hot planets. Two other potentially important effects, that are not included in our setup, are meridional advection and photochemistry. Regarding the former, it remains to be seen to what degree, and in which planets, meridional quenching can impact the chemical composition. Regarding the latter, we aim to investigate the impact of photochemistry over a broad parameter range in a detailed follow-up study. Despite these modelling limitations, we generally find reasonable agreement between the atmospheric models presented in this work, and previous studies (Agúndez et al. 2014a;Mendonça et al. 2018b;Venot et al. 2020b;Drummond et al. 2020). We find that the existing discrepancies are often caused by differences in the temperature structure and 𝐾 𝑧𝑧 . This sensitivity highlights that care should be taken in the derivation of the atmospheric thermal structure, especially when vertical mixing is efficient, and the quenching pressures are high (Fortney et al. 2020). In the future, it would be beneficial to investigate the interior temperature impact on the observable atmosphere in a 3D context.

APPENDIX A: ASSESSMENT OF MODEL CONVERGENCE

In order to ensure that the models presented in this work are converged to a steady-state, in which they do not exhibit substantial, non-cyclic evolution, diagnostic parameters have been plotted as a function of simulation time for a typical general circulation model (Fig. A1) and chemical kinetics simulation (Fig. A2).

A1 GCM Convergence

The convergence rate of irradiated exoplanet GCMs is a longstanding subject of scrutiny. Especially the deep atmospheric layers (𝑝 > 10 bar), which have a comparatively long radiative timescale, have been noted to evolve very slowly (e.g. Amundsen et al. 2016), leading to impractically long computation time requirements (Wang & Wordsworth 2020;Mendonça 2020). Moreover, it has been shown that, in some cases, the deep atmospheric evolution can have a dynamical feedback on the observable atmospheric layers above (Carone et al. 2020). In order to combat the slow spin-up times of the deep atmospheric layers, Carone et al. (2020) have prescribed an accelerated forcing for these layers on to a pre-computed adiabatic temperature profile. A different method, that has been proposed by Sainsbury-Martinez et al. (2019), is to initialize the temperature structure with a hot adiabat. We have employed the latter method in this work. In addition to the basal drag applied in our models, this serves to stabilize and spin up the atmosphere on a time-scale of ∼1000 (Earth) days, as can be seen in Fig. A1.

A2 Chemical Kinetics Convergence

The pseudo-2D chemical kinetics model used in this work also relies on a convergence of the initial profile to a steady-state. The model is initialized with a vertically quenched chemistry profile corresponding to the substellar point, following the recommendation of Agúndez et al. (2012), after which it is run for 100 rotations. We demonstrate that 100 rotations are sufficient to ensure the convergence for the models in this grid. The total simulation time of the models in our grid varies depending on the zonal wind speed 𝑢 (see Fig. 3, E3, E4) as 𝑡 = 100 × 2 𝜋𝑅 𝑝 𝑢 , with 𝑅 𝑝 the fixed planet radius of 1.35 R Jup . Hence, models with fast equatorial jet streams are run for shorter times than models with slow equatorial flow. As a consequence, generally, models with high effective temperatures and slow rotation rates are run for ∼10 7 s, while the coolest models are run for ∼10 10 s. A formal requirement for convergence is that the total simulation time is higher than the vertical transport time-scale of the model. Since 𝐾 𝑧𝑧 scales with the effective temperature (Section 4.2.3), the simulation times roughly match the maximal vertical advection times needed for convergence, namely 10 6 -10 8 s for the high-, 10 6 -10 10 s for the medium-, and 10 7 -10 10 s for the low-gravity models. However, since we start the pseudo-2D model from a converged 1D model at substellar point, the initial composition is already close to the expected steady-state, and thus, the actual convergence should be faster.

To empirically check for convergence after 100 rotations, we plotted the number fractions at the morning terminator (-90 • ) for some dominant molecules as a function of the total simulation time after 100 rotations (Fig. A2). The model corresponds to a planet with 𝑇 eff = 1400 K and g=10 m/s 2 . We purposely plot a model with moderate temperature, since it would take the longest time to converge to a periodic steady-state. The required convergence times for cold and hot models would be lower, because of the small longitudinal temperature differences in the former, and short chemical time-scales in the latter case.

The different molecules shown in the evolution plot (Fig. A2) do not all evolve to a steady-state equally fast. Species such as H 2 O and CO do not show deviations from their vertically quenched initial state. Other species, such as CH 4 and NH 3 do show evolution, because their initial concentration is further from the steady-state and they are zonally quenched. The nitrogen-bearing molecules NH 3 and HCN exhibit the slowest evolution. After about 10 7 s, most of the species have attained a final steady-state.

We conclude that an integration time of 100 rotations, as we have used in this work for each chemical model, is sufficient to reach a steady-state when starting from a converged 1D thermochemical kinetics model. time-scales. Therefore, these layers are likely to be well-mixed and have a common interior entropy (Barman et al. 2004). For fully convective models, it is a physical requirement that the irradiated and non-irradiated hemispheres of the same object reach a constant interior entropy at depth, as was argued for irradiated, convective stars in binary systems by Vaz & Nordlund (1985); Brett & Smith (1993). Furthermore, 3D radiative transfer computations of irradiated planets also indicate a common radiative equilibrium in the deep layers of the atmosphere (Hauschildt et al. 2008).

Mixing towards a constant entropy value in the deep interior of convective atmospheres involves horizontal transport of energy from the irradiated side towards the non-irradiated side. This is, however, impossible to model rigorously in the 1D framework of petitCODE. The total emerging flux of a converged petitCODE model is equal to the sum of the intrinsic and irradiation fluxes of the planet (see Mollière et al. 2015):

𝐹 emerging = 𝜎 𝑇 4 int + 𝜇 𝑇 4 * 4 𝑅 * 𝑑 2 , ( B1
)
where 𝑅 * and 𝑇 * are the stellar radius and temperature, 𝑑 is the distance from the star, and 𝜎 is the Stefan-Boltzmann constant. Increasing the angle of incidence 𝜇 is equivalent to increasing the irradiation on the planet. This additional flux can reach the top layers of the convective zone and, since the temperature gradient is set by the adiabatic gradient, the extra heating causes an increase in temperature and entropy throughout the entire deep convective spheres, due to molecular opacities rather than a continuum opacity source (Fortney 2018). In general, radiation originates from pressures between 0.01 and 1 bar.

region. To counteract this and obtain a constant interior entropy, the intrinsic temperature can be lowered for larger angles of incidence and raised for smaller angles of incidence, thereby mimicking the horizontal transport of energy in these deep convective layers (Barman et al. 2004). Thus, the intrinsic temperatures associated with each irradiation angle are reduced or increased to match a common adiabat. In this study, the common adiabatic temperature gradient is taken to be the one associated with the model of globally-averaged energy redistribution (indicated as the red line in Fig. B1). We adopted a constant intrinsic temperature 𝑇 int = 200 K for all planetary-averaged models, to facilitate both the intercomparison between the models in our grid, and the comparison with earlier works using petitCODE (Mollière et al. 2015;Molaverdikhani et al. 2019). However, it is important to note that hot Jupiters may exhibit intrinsic temperatures that are much higher, ranging from 100 K to 700 K, depending on their irradiation and radius inflation (Thorngren et al. 2019). Furthermore, the intrinsic temperature and the extent of radius inflation have been shown to impact the atmospheric circulation of hot Jupiters (Carone et al. 2020). A dedicated study of the impact of the intrinsic temperature, however, is outside the scope of this study. The pressure-temperature profiles of a locus corresponding to a certain irradiation angle then have their intrinsic temperature adjusted, so that they end up on this planetary-averaged adiabat. . The radiative-equilibrium temperature for different angles of incidence, plotted per model as separate lines, increases with effective temperature, except for the coldest profile, which corresponds to the planet night side. Gravity acts to shift the temperature structure vertically to higher/lower pressures. This selection of petitCODE (Mollière et al. 2015) models corresponds to models with effective temperatures of 400 K, 1400 K, and 2400 K, and surface gravities of 1 m/s 2 , 10 m/s 2 , and 100 m/s 2 .

dence. These models have been constructed using the methodology laid out in Appendix B. Increasing the effective temperature of the exoplanet atmosphere leads to an increase in radiative-equilibrium temperature for all angles of incidence, except that of the night side (Fig. D1).

APPENDIX E: ADDITIONAL FIGURES FOR LOW AND HIGH SURFACE GRAVITIES

Here, we present additional figures of the low-(𝑔 = 1 m/s 2 ) and high-(𝑔 = 100 m/s 2 ) gravity models in our grid, which are analogous to the medium gravity models shown in Section 4. These consist of isobaric temperature maps (Fig. E1 and E2), zonal-mean zonal wind plots (Fig. E3 and E4), pressure-temperature profiles (Fig. E5 andE6), vertical wind maps (Fig. E7 and E8), and pseudo-2D chemistry plots (Fig. E9 andE10). For the low-gravity chemistry models, models with effective temperatures above 2000 K are omitted. For these cases, computations with the pseudo-2D chemical kinetics code did not result in physically acceptable solutions. The reason for the numerical problems has been identified to be linked to the atmospheric column height, which becomes disproportionately high in comparison with the planet radius, resulting in unphysically low column densities. Furthermore, planets with a low gravity of 𝑔 = 1 m/s 2 and temperatures above 2000 K are not expected to retain their atmosphere, due to thermal evaporation. Therefore, we choose to omit these ultra-hot models for the low-gravity case.

Figure 1 .

 1 Figure 1. Schematic overview of the different models used in this work, and the post-processing sequence in which they are linked.

Figure 2 .

 2 Figure 2. Temperature maps for the 𝑔 = 10 m/s 2 models in our grid, plotted as isobaric slices at a pressure of 70 mbar. Grey arrows show the direction of the local horizontal winds, and the arrow length is scaled by the wind speed. The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically). The substellar point is located at (0 • , 0 • ).

Figure 3 .

 3 Figure 3. Zonally averaged zonal wind plots for the 𝑔 = 10 m/s 2 models in our grid. Positive values of the wind speed (blue) correspond to prograde wind flow. The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically).

Figure 4 .

 4 Figure 4. Temperature-pressure profiles from the 𝑔 = 10 m/s 2 GCM models in our grid. Temperatures are meridionally averaged between ±20 • , shown for longitudes spaced 30 • apart (-180• , -150 • ,. . . , 150 • ), and colour coded with a gradual change from the anti-stellar (blue) to substellar point (yellow). Full and dashed lines mark locations west and east of the substellar point. The shaded background colour indicates the parameter space where methane is more abundant than CO in chemical equilibrium(Visscher 2012). The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically).

  Figure5. The relative mean temperature difference between the day-and night side is systematically larger for fast-rotating models (orbiting M5-stars) than slow-rotating models (orbiting F5-stars). Furthermore, it increases with planetary effective temperature. The relative mean temperature difference was computed by subtracting the average night-side temperature from the average day-side temperature, and dividing by the average day-side temperature. To compute the average day-or night-side temperature, we have taken the arithmetic mean of the temperature 𝑇 ( 𝜙, 𝜆, 𝑝) in our GCM model for all pressures 𝑝, latitudes (weighted with cos 𝜙) and longitudes 𝜆 of the day or night side respectively. The three panels show the low-, medium-and high-gravity models.

Figure 6 .

 6 Figure 6. Maps of the vertical wind speed on a 70 mbar isobar for the 𝑔 = 10 m/s 2 models in our grid. Upwelling motions (blue) are mostly present on the day side (longitudes between -90 • and +90 • ), whereas downwelling flow (red) is mostly visible on the night side, often as very localized fronts. The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically). Thus, the rotation rate decreases along rows.

Figure 7 .

 7 Figure 7. The eddy diffusion coefficient profiles, corresponding to all planetary climate simulations, show a clear increasing trend with effective temperature, and a decreasing trend with surface gravity. For each stellar type (columns), models of different gravity are displayed (rows). Dark blue and bright yellow correspond to the coolest and hottest planets in the grid, 400 K and 2600 K, respectively. The black dotted line is a reference profile, corresponding to the parametrization given by Parmentier et al. (2013) for HD 209458 b.

Figure 8 .

 8 Figure 8. The abundances for the main molecular species, computed with the pseudo-2D chemistry code for 𝑔 = 10 m/s 2 , show a dichotomy of the zonal mixing regime. On the one hand, hotter planets (bottom rows) sustain large chemical gradients between the hot day-and the cold night sides, as evidenced by the 12 vertical profiles that are plotted for different longitudes (-180 • , -150 • ,. . . , 150 • ). Colder planets (top rows), on the other hand, exhibit zonally quenched chemistry profiles, without any day-night variation. In these cases, all longitudinally sampled vertical profiles are plotted on top of each other. The dashed lines denote the chemical equilibrium composition at the substellar point.

Figure 9 .

 9 Figure9. The vertical advection time-scale 𝜏 vertical is shorter than the zonal advection time-scale 𝜏 zonal for the upper part of this hot Jupiter model atmosphere (with 𝑇 eff = 1400 K, 𝑔 = 10 m/s 2 and G5-star host). In the deeper atmosphere, the zonal mixing time-scale is much shorter. The pressure range in which the chemical abundances become quenched for this model, has been marked in grey. The time-scale associated with the uniform zonal wind used in the pseudo-2D computation, is marked with a solid red line.

Figure 16 .

 16 Figure 16. Physical diagnostics for an HD 209458 b-like planet. This grid model has 𝑇 eff = 1400 K, 𝑔 = 10 m/s 2 and host star type G5. Plotted are the meridional mean temperature used in the pseudo-2D chemistry code (a), the chemical composition (b), and the synthetic transmission spectrum, in a comparison with observations (c). The observational data (grey) are from Sing et al. (2016). The transmission spectrum is cloud-free, and does not contain Na/K opacities. In the chemistry plot, the equilibrium compositions of the morning and evening limb are plotted as thin dashed lines.

  in this work, figure E.1 in Drummond et al. (2020)), we find that our time-scales are about 2 orders of magnitude shorter near 1 bar, which could explain the discrepancy. The reason for the time-scale difference is not easily unravelled, since Drummond et al. (2020) have estimated it directly based on the vertical wind flow 𝑤 as 𝜏 vertical = 𝐻 |𝑤 | , whereas we have based our estimation on the eddy diffusion coefficient as 𝜏 vertical = 𝐻 2 /𝐾 𝑧𝑧 . Upon closer comparison, it seems the 1 bar-temperatures in our model (∼1400 K) are slightly lower than those in Drummond et al. (

Figure 17 .

 17 Figure 17. Physical diagnostics for an HD 189733 b-like planet. This grid model has 𝑇 eff = 1200 K, 𝑔 = 10 m/s 2 and host star type K5. Plotted are the meridional mean temperature used in the pseudo-2D chemistry code (a), the chemical composition (b), and the synthetic transmission spectrum, in a comparison with observations (c). The observational data (grey) are from Sing et al. (2016). The transmission spectrum is cloud-free, and does not contain Na/K opacities. In the chemistry plot, the equilibrium compositions of the morning and evening limb are plotted as thin dashed lines.

Figure 18 .

 18 Figure 18. Physical diagnostics for an WASP-43 b-like planet. This grid model has 𝑇 eff = 1400 K, 𝑔 = 100 m/s 2 and host star type K5. Plotted are the meridional mean temperature used in the pseudo-2D chemistry code (a), the chemical composition (b), and the synthetic transmission spectrum, in a comparison with observations (c). The observational data (grey) are fromKreidberg et al. (2014). The transmission spectrum is cloud-free, and does not contain Na/K opacities. In the chemistry plot, the equilibrium compositions of the morning and evening limb are plotted as thin dashed lines.

Figure 19 .

 19 Figure 19. Physical diagnostics for an WASP-121 b-like planet. This grid model has 𝑇 eff = 2400 K, 𝑔 = 10 m/s 2 and host star type F5. Plotted are the meridional mean temperature used in the pseudo-2D chemistry code (a), the chemical composition (b), and the synthetic transmission spectrum, in a comparison with observations (c). The observational data (grey) are from Evans et al. (2016). The transmission spectrum is cloud-free, and does not contain Na/K opacities. In the chemistry plot, the equilibrium compositions of the morning and evening limb are plotted as thin dashed lines.

Figure 20 .

 20 Figure 20. Physical diagnostics for an WASP-107 b-like planet. This grid model has 𝑇 eff = 800 K, 𝑔 = 1 m/s 2 and host star type K5. Plotted are the meridional mean temperature used in the pseudo-2D chemistry code (a), the chemical composition (b), and the synthetic transmission spectrum, in a comparison with observations (c). The observational data (grey) are from Spake et al. (2018). The transmission spectrum is cloud-free, and does not contain Na/K opacities. In the chemistry plot, the equilibrium compositions of the morning and evening limb are plotted as thin dashed lines.

Figure 21 .

 21 Figure 21. This transmission spectrum, corresponding to a simulation with 𝑇 eff = 1400 K, 𝑔 = 10 m/s 2 and stellar type K5, shows a distinct opacity bump at 4.3 𝜇m because of absorption by CO 2 . The inset shows the CO 2 absorption band in more detail.

Figure 22 .

 22 Figure 22. In this transmission spectrum, corresponding to a simulation with 𝑇 eff = 1200 K, 𝑔 = 10 m/s 2 and stellar type M5, a clear opacity band of NH 3 at 10 𝜇m can be seen. The inset shows a detail of the NH 3 absorption feature.

Figure A1 .

 A1 Figure A1. The root mean square (RMS) wind speed on isobars of this 3D GCM simulation (𝑇 eff = 1400 K, 𝑔 = 10 m/s 2 , G5) shows a quick evolution in the first 400 days of the numerical integration, before attaining a steadystate. The RMS wind speed is weighted with the cosine of the latitude, and horizontally averaged on isobars. The simulation is evaluated after 1500 days and averaged over 100 days.

Figure A2 .

 A2 Figure A2. The number fractions at the night side of this exoplanet atmosphere model (𝑇 eff = 1400 K, 𝑔 = 10 m/s 2 , G5) show convergence to a steady-state situation. Different species show different evolution speeds with respect to simulation time of the pseudo-2D chemistry code. Note that the colour maps for each molecule are different. White contours are spaced at intervals of 0.5 in log-space.

Figure B1 .

 B1 Figure B1. The pressure-temperature profiles corresponding to locations associated with different irradiation angles 𝜇, but constant intrinsic temperatures 𝑇 int = 200 K, end up on different adiabats in the deep convective part of the atmosphere. The profiles are calculated for a planet with 𝑇 eff = 1400 K, 𝑔 = 10 m s -2 and solar metallicity around a G5 star. The night-side profile (𝑚𝑢=0) is not irradiated, so the temperature-pressure profile is purely determined by the intrinsic temperature 𝑇 int . The bold part of the profiles signifies a convective regime. The red line shows the petitCODE model associated with the planetary averaged distribution of the incident flux.

  Figure D1. The radiative-equilibrium temperature for different angles of incidence, plotted per model as separate lines, increases with effective temperature, except for the coldest profile, which corresponds to the planet night side. Gravity acts to shift the temperature structure vertically to higher/lower pressures. This selection of petitCODE(Mollière et al. 2015) models corresponds to models with effective temperatures of 400 K, 1400 K, and 2400 K, and surface gravities of 1 m/s 2 , 10 m/s 2 , and 100 m/s 2 .

Figure E1 .

 E1 Figure E1. Temperature maps for the 𝑔 = 1 m/s 2 models in our grid, plotted as isobaric slices at a pressure of 70 mbar. Grey arrows show the direction of the local horizontal winds, and the arrow length is scaled by the wind speed. The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically). The substellar point is located at (0 • , 0 • ).

Figure E2 .

 E2 Figure E2. Temperature maps for the 𝑔 = 100 m/s 2 models in our grid, plotted as isobaric slices at a pressure of 70 mbar. Grey arrows show the direction of the local horizontal winds, and the arrow length is scaled by the wind speed. The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically). The substellar point is located at (0 • , 0 • ).

Figure E3 .

 E3 Figure E3. Zonally averaged zonal wind plots for the 𝑔 = 1 m/s 2 models in our grid. Positive values of the wind speed (blue) correspond to prograde wind flow. The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically).

Figure E4 .

 E4 Figure E4. Zonally averaged zonal wind plots for the 𝑔 = 100 m/s 2 models in our grid. Positive values of the wind speed (blue) correspond to prograde wind flow. The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically).

Figure E5 .

 E5 Figure E5. Temperature-pressure profiles from the 𝑔 = 1 m/s 2 GCM models in our grid. Temperatures are meridionally averaged between ±20 • , shown for longitudes spaced 30 • apart (-180 • , -150 • ,. . . , 150 • ), and colour coded with a gradual change from the anti-stellar (blue) to substellar point (yellow). Full and dashed lines mark locations west and east of the substellar point. The shaded background colour indicates the parameter space where methane is more abundant than CO in chemical equilibrium(Visscher 2012). The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically).

Figure E6 .

 E6 Figure E6. Temperature-pressure profiles from the 𝑔 = 100 m/s 2 GCM models in our grid. Temperatures are meridionally averaged between ±20 • , shown for longitudes spaced 30 • apart (-180 • , -150 • ,. . . , 150 • ), and colour coded with a gradual change from the anti-stellar (blue) to substellar point (yellow). Full and dashed lines mark locations west and east of the substellar point. The shaded background colour indicates the parameter space where methane is more abundant than CO in chemical equilibrium(Visscher 2012). The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically).

Figure E7 .

 E7 Figure E7. Maps of the vertical wind speed on a 70 mbar isobar for the 𝑔 = 1 m/s 2 models in our grid. Upwelling motions (blue) are mostly present on the day side (longitudes between -90 • and +90 • ), whereas downwelling flow (red) is mostly visible on the night side, often as very localized fronts. The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically). Thus, the rotation rate decreases along rows. For the sake of comparison, the colour scaling is kept consistent with Fig.6and E8.

Figure E8 .

 E8 Figure E8. Maps of the vertical wind speed on a 70 mbar isobar for the 𝑔 = 100 m/s 2 models in our grid. Upwelling motions (blue) are mostly present on the day side (longitudes between -90 • and +90 • ), whereas downwelling flow (red) is mostly visible on the night side, often as very localized fronts. The plots are ordered by host star type (horizontally) and planetary effective temperature (vertically). Thus, the rotation rate decreases along rows. For the sake of comparison, the colour scaling is kept consistent with Fig.6and E7.

Figure E9 .

 E9 Figure E9. The abundances for the main molecular species, computed with the pseudo-2D chemistry code for 𝑔 = 1 m/s 2 are shown as 12 vertical profiles plotted at different longitudes (-180 • , -150 • ,. . . , 150 • ). The dashed lines denote the chemical equilibrium composition at the substellar point. Note that models with effective temperatures of 2200 K, 2400 K and 2600 K are omitted.

Table 1 .

 1 Overview of the model properties in the grid, ordered according to their effective temperatures 𝑇 eff (vertically) and host star type (horizontally).

	𝑇 eff		M5-star K5-star G5-star F5-star
	400 K	𝑎 (AU) 𝑃 orb (days)	0.049 8.61	0.186 35.55	0.421 104.33	0.859 256.13
	600 K	𝑎 (AU) 𝑃 orb (days)	0.021 2.45	0.080 10.13	0.182 29.73	0.372 72.98
	800 K	𝑎 (AU) 𝑃 orb (days)	0.012 1.03	0.045 4.25	0.102 12.46	0.208 30.59
	1000 K	𝑎 (AU) 𝑃 orb (days)	0.008 0.53	0.029 2.17	0.065 6.37	0.133 15.64
	1200 K	𝑎 (AU) 𝑃 orb (days)	0.005 0.30	0.020 1.26	0.045 3.68	0.092 9.04
	1400 K	𝑎 (AU) 𝑃 orb (days)	0.004 0.19	0.015 0.79	0.033 2.32	0.068 5.69
	1600 K	𝑎 (AU) 𝑃 orb (days)	0.003 0.13	0.011 0.53	0.025 1.55	0.052 3.81
	1800 K	𝑎 (AU) 𝑃 orb (days)	0.002 0.09	0.009 0.37	0.020 1.09	0.041 2.68
	2000 K	𝑎 (AU) 𝑃 orb (days)	0.002 0.07	0.007 0.27	0.016 0.80	0.033 1.95
	2200 K	𝑎 (AU) 𝑃 orb (days)	0.002 0.05	0.006 0.20	0.013 0.60	0.027 1.47
	2400 K	𝑎 (AU) 𝑃 orb (days)	0.001 0.04	0.005 0.16	0.011 0.46	0.023 1.13
	2600 K	𝑎 (AU) 𝑃 orb (days)	0.001 0.03	0.004 0.12	0.010 0.36	0.020 0.89

  has been shown that the transit depth of molecular features

		10	HCN NH 3 CO 2 CH 4 H 2 O CO		1.37	Agnostic 1D average
		10									1.36
											1.35
		10	9	8	7	6	5	4	3		1.37	10 0 3D GCM + equilibrium chemistry	10 1
		10									1.36
											1.35
		10	9	8	7	6	5	4	3		1.37	10 0 3D GCM + 1D chemical kinetics	10 1
		10									1.36
											1.35
	Pressure (bar)	10 10	9 9	8 8 Log number fraction 7 6 5 4 7 6 5 4	3 3	Transit radius (R Jup )	1.35 1.36 1.37	10 0 10 0 3D GCM + pseudo-2D chemical kinetics Wavelength (micron)	10 1 10 1

Figure

13

. The same modelling hierarchy as in Fig.

12

is shown, but here for a cooler planet with 𝑇 eff = 800 K.

  The only change brought about by vertical mixing are the increased CO and CO 2

		10	HCN NH 3 CO 2 CH 4 H 2 O CO		1.38	Agnostic 1D average
		10									1.36
											1.34
		10	9	8	7	6	5	4	3		1.38	10 0 3D GCM + equilibrium chemistry	10 1
		10									1.36
											1.34
		10	9	8	7	6	5	4	3		1.38	10 0 3D GCM + 1D chemical kinetics	10 1
		10									1.36
											1.34
	Pressure (bar)	10 10	9 9	8 8 Log number fraction 7 6 5 4 7 6 5 4	3 3	Transit radius (R Jup )	1.34 1.36 1.38	10 0 10 0 3D GCM + pseudo-2D chemical kinetics Wavelength (micron)	10 1 10 1

Figure 14. The same modelling hierarchy as in Fig. 12 is shown, but here for a hotter planet with 𝑇 eff = 2000 K.

  al. 2016;

	Jup )	0.025		G5 F5
	Max limb difference (R	0.005 0.010 0.015 0.020		M5 K5
		0.000
			400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600
			Effective temperature (K)
		1.37	800 K
		1.36
	Jup ) Transit radius (R	1.35 1.34 1.33 1.34 1.35 1.36 1.37 1.38 1.38 1.39	1400 K 2000 K	evening limb averaged morning limb morning limb evening limb averaged
		1.37
		1.36
		1.32 1.33 1.34 1.35	10 0	10 1 morning limb evening limb averaged
			Wavelength (micron)

MNRAS 000,1-41 (2020)

One way to do this would be to use interpolation on a pre-calculated chemical equilibrium grid, at the pressure and temperature under consideration. See e.g. https://petitradtrans.readthedocs.io/en/ latest/content/notebooks/poor_man.html.MNRAS 000,1-41 (2020)

In this and following plots, 70 mbar was selected because it is a suitable pressure region to diagnose the dynamical-radiative coupling, and to facilitate comparisons with earlier work (e.g.Komacek & Showman 2016;Komacek et al. 2017).MNRAS 000,1-41 (2020)

In determining the best matching model to the physical planetary parameters, we first match the temperature and gravity as close as possible. Next, we match the rotation/host star type. In doing so, we aim to obtain a good match for the planetary rotation rate, rather than the actual stellar type of the host star, since the rotation rate has a much bigger impact on the atmospheric modelling.MNRAS 000,1-41 (2020)

The photosphere, the pressure layer at which the atmosphere becomes opaque to radiation, varies strongly with wavelength for planetary atmo-

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000,1-41 (2020)

ACKNOWLEDGEMENTS

We wish to thank Maria Steinrueck, William Pluriel and Tiziano Zingales for helpful discussions. We also thank the anonymous referee for a very helpful review that led to the improvement of this manuscript. RB is a PhD fellow of the Research Foundation -Flanders (FWO). LD acknowledges support from the FWO research grant G086217N. LC acknowledges support by the DFG grant CA

DATA AVAILABILITY STATEMENT

The data underlying this article will be shared on reasonable request to the corresponding author.

APPENDIX B: COMMON ADIABAT

The intrinsic temperature 𝑇 int used in the computation of our pe-titCODE models is associated with the intrinsic flux coming from the planet's interior. When pressure-temperature profiles for different irradiation angles 𝜇 are computed with a constant 𝑇 int , these profiles will end up at different adiabats in the deep, convective part of the atmosphere (see Fig. B1). However, for the construction of our radiative-equilibrium profiles, we assumed that each profile corresponding to a certain irradiation angle converges to the same common adiabatic temperature gradient in the deep part of the atmosphere. To achieve this convergence, a different intrinsic temperature 𝑇 int was imposed for every irradiation angle. This approach was previously used in 1D models of irradiated giant planets (Barman et al. 2005) and in Newtonian-relaxed 3D GCMs for the hot Jupiters HD 209458 b and HD 189733 b (Showman et al. 2008).

The assumption of a common adiabatic temperature lapse rate for all incidence angles is physically motivated. As these atmospheric layers are expected to be convective and lie well below the photosphere 5 , dynamical time-scales will be shorter than radiative

APPENDIX C: THE ASSUMPTION OF SYNCHRONOUS ROTATION

Synchronous rotation is attained by tidal interaction of the planet with its host star. The efficiency of this interaction depends on the the physical properties of the star and the planet, and the distance between them. Since the semi-major axis in our grid is calculated based on the effective temperature of the planet and the star, it might be that for some combination of model parameters the tidal interaction efficiency is very low, and synchronous rotation may not be a valid assumption.

In order to investigate the validity of synchronous rotation, we computed order-of-magnitude estimations of the synchronization times for each of our grid models, using equation (1) in Showman & Guillot (2002):

Here, 𝑀 𝑝 , 𝑅 𝑝 are the planetary mass and radius, 𝑀 * is the stellar mass, 𝐺 is the gravitational constant, 𝑎 is the semi-major axis, 𝜔 is the starting rotational angular velocity and 𝜔 syn is the synchronized angular velocity, namely equal to the orbital angular velocity. Finally, 𝑄 is the tidal dissipation factor. We followed the approach of Showman & Guillot (2002) in adopting 𝑄 = 10 5 and choosing 𝜔 = 1.76 • 10 -4 s -1 , which is the current rotation rate of Jupiter. Likewise, Jupiter values 𝑀 𝑝 = 1 𝑀 Jup and 𝑅 𝑝 = 1 𝑅 Jup were adopted for each model. Other parameters have been calculated per grid model based on Section 3.2. In computing the synchronization times, we have assumed that increasing and decreasing the rotation rate of the planet due to tidal interaction is equally efficient [START_REF] Murray | Solar system dynamics[END_REF]. In Fig. C1, the synchronization time-scales based on our orderof-magnitude calculation are shown for all grid model configurations. We find that for all hot models (𝑇 eff > 1000 K) and almost all models orbiting cool (M5, K5) stars, the synchronization timescales are shorter than 100 Myr, which is at least one order of mag- nitude below the typical age of exoplanets (e.g. Silva Aguirre et al. 2015;Maxted et al. 2015). Hence, for these models synchronous rotation is a reasonable assumption. On the other hand, the coldest models orbiting G5-and F5-stars have synchronization time-scales higher than 10 Gyr and semi-major axes above 0.1 AU (see Table 1). This makes it quite unlikely for these planets to be found in a 1:1 spin-orbit resonance, and for this subset of our grid synchronous rotation is not a valid assumption. However, since this work is mainly a parameter study, we opt to still include the more improbable model configurations as well, since they can facilitate the derivation of general trends in the parameter space. For detailed studies of the atmospheric circulation on asynchronously rotating planets, we refer to Penn & Vallis (2018)

APPENDIX D: SAMPLE OF RADIATIVE-CONVECTIVE EQUILIBRIUM MODELS

We present a selection of petitCODE radiative-convective equilibrium temperatures as a function of pressure and the angle of inci-