3D mixing in hot Jupiters atmospheres.
I. Application to the day/night cold trap in HD 209458b
Key Words.:
Planets and satellites: atmospheres  Methods: numerical  DiffusionAbstract
Context:Hot Jupiters exhibit atmospheric temperatures ranging from hundreds to thousands of Kelvin. Because of their large daynight temperature differences, condensable species that are stable in the gas phase on the dayside—such as TiO and silicates—may condense and gravitationally settle on the nightside. Atmospheric circulation may counterbalance this tendency to gravitationally settle. This threedimensional (3D) mixing of condensable species has not previously been studied for hot Jupiters, yet it is crucial to assess the existence and distribution of TiO and silicates in the atmospheres of these planets.
Aims:We investigate the strength of the nightside cold trap in hot Jupiters atmospheres by investigating the mechanisms and strength of the vertical mixing in these stably stratified atmospheres. We apply our model to the particular case of TiO to address the question of whether TiO can exist at low pressure in sufficient abundances to produce stratospheric thermal inversions despite the nightside cold trap.
Methods:We modeled the 3D circulation of HD 209458b including passive (i.e.radiatively inactive) tracers that advect with the 3D flow, with a source and sink term on the nightside to represent their condensation into haze particles and their gravitational settling.
Results:We show that global advection patterns produce strong vertical mixing that can keep condensable species aloft as long as they are trapped in particles of sizes of a few microns or less on the night side. We show that vertical mixing results not from smallscale convection but from the largescale circulation driven by the daynight heating contrast. Although this vertical mixing is not diffusive in any rigorous sense, a comparison of our results with idealized diffusion models allows a rough estimate of the effective vertical eddy diffusivities in these atmospheres. The parametrization , valid from 1 bar to a few bar, can be used in 1D models of HD 209458b. Moreover, our models exhibit strong spatial and temporal variability in the tracer concentration that could result in observable variations during either transit or secondary eclipse measurements. Finally, we apply our model to the case of TiO in HD 209458b and show that the daynight cold trap would deplete TiO if it condenses into particles bigger than a few microns on the planet’s nightside, keeping it from creating the observed stratosphere of the planet.
Conclusions:
1 Introduction
The year 1988 was marked by the discovery of the first substellar object outside our solar system (Becklin & Zuckerman 1988). This discovery was followed by numerous other discoveries of cool brown dwarfs, whose spectra differ significantly from stars. This led to the definition of three new spectral classes beyond the stellar classes O to M: the L, T, and Y dwarfs, ranging from temperatures of to (Cushing et al. 2011). To understand of these brown dwarfs and the mechanisms that cause the transition from one spectral class to another, it is fundamental to take condensation processes into account: in cold atmospheres molecules can form, condense, and rain out, affecting the observed spectral features of these objects. In particular, the M to L transition is marked by the disappearance of the titanium oxide (TiO) bands, which is understood by its transformation into titanium dioxide (TiO) and condensation into perovskite (CaTiO) (Lodders 2002). The L to T transition is understood as the switch from a cloudy atmosphere to a cloudfree atmosphere as the cloud layer migrates into the deep, unobservable regions of the star (see Kirkpatrick (2005) for a review). Because brown dwarfs produce their own light, highquality spectra of their atmosphere can be measured. Numerous atmospheric models have been built to fit the data, from simple models with a reduced set of free parameters (Ackerman & Marley 2001) to sophisticated models including the detailed physics of condensation, growth, and settling of particles (Woitke & Helling 2004).
The discovery and characterization of exoplanets followed close behind. The first detection of a planet around a main sequence star by Mayor & Queloz (1995) opened the trail for discovering hundreds of exoplanets. Some years later, atmospheric characterization of these objects became available, both from transit spectroscopy (Charbonneau et al. 2002) and from direct detection of the planet’s thermal emission (Charbonneau et al. 2005). Although the globalmean effective temperatures of these planets are similar to those on brown dwarfs, a major difference is that hot Jupiters are strongly irradiated. Depending on the incident stellar flux, planetary rotation rate, and other factors, atmospheric circulation models show that this daynight heating gradient can lead, at low pressures, to nightside temperatures that are at least 1000 K colder than dayside temperatures (e.g., Showman et al. 2008, 2009; DobbsDixon et al. 2010; Rauscher & Menou 2010, 2012a; Heng et al. 2011a, b; Perna et al. 2012). As a result, a wide variety of condensable species that are stable in the gas phase on the dayside may condense on the nightside, leading to the formation of particles there. In the absence of atmospheric vertical mixing, such particles would gravitationally settle, depleting the atmosphere of these species on both the dayside and the nightside. Sufficiently strong vertical mixing, anywhere on the planet, however, may keep these particles suspended in the atmosphere, allowing them to sublimate into the gas phase in any air transported from nightside to dayside. Thus, the existence of this nightside “cold trap” is crucial for understanding not only the existence of hazes on hot Jupiters (e.g., Pont et al. 2013) but also the gasphase composition of the atmosphere on both the dayside and the nightside.
These arguments are relevant to a wide range of titanium and vanadium oxides, silicate oxides, and other species. In particular, chemicalequilibrium calculations show that, at temperatures of 1000–2000 K, there exist a wealth of condensates including TiO, TiO, and TiO (among other titanium oxides), MgAlO, MgSiO, MgSiO, NaAlSiO, KAlSiO, and several phosphorus oxides (Burrows & Sharp 1999; Lodders 2002). Understanding the possible existence of gasphase TiO, Na, K, and other species on the dayside therefore requires an understanding of the nightside cold trap.
A particularly interesting problem in this regard is the existence of exoplanet stratospheres. For some transiting hot Jupiters, Spitzer IRAC secondaryeclipse observations indicate the presence of thermal inversions (stratospheres) on these planets’ daysides (Knutson et al. 2008). These stratospheres are generally thought to result from absorption of starlight by strong visible/ultraviolet absorbers, but debate exists about the specific chemical species that are responsible. Hubeny et al. (2003) and Fortney et al. (2008) showed that, because of their enormous opacities at visible wavelengths, the presence of gaseous titanium and vanadium oxides can lead to stratospheres analogous to those inferred on hot Jupiters. This hypothesis is supported by the possible detection of TiO by Désert et al. (2008) in the atmospheric limb of HD 209458b. Then, a crucial question is whether the nightside cold trap would deplete the atmosphere of TiO, preventing this species from serving as the necessary absorber (Showman et al. 2009; Spiegel et al. 2009).
As pointed out by several authors, there exists another possible cold trap. 1D radiativetransfer models suggest that, on some hot Jupiters, the globalmean temperaturepressure profile becomes sufficiently cold for condensation of gaseous TiO to occur at pressures of tens to hundreds of bars (e.g., Fortney et al. 2008). Even though a stratosphere on such a planet would be sufficiently hot for TiO—if present—to exist in the gas phase, the condensation of TiO and downward settling of the resulting grains at 10–100 bars might prevent the existence of TiO in the atmosphere (Showman et al. 2009; Spiegel et al. 2009). The strength of these two cold traps is given by a competition between gravitational settling and upward mixing. The vertical mixing results from complex 3D flows and can be inhomogeneous over the planet (Cooper & Showman 2006; Heng et al. 2011a). When using a 1Dmodel, the vertical mixing is usually considered to be diffusive only, and parametrized by a diffusion coefficient . Several studies give an estimate for this vertical diffusion coefficient in hot Jupiters atmospheres. Heng et al. (2011a) uses the magnitude of the Eulerian mean streamfunction as a proxy for the strength of the vertical motions and derived a vertical mixing coefficient of the order of . Other authors used an estimate based on the root mean squared vertical velocity either the local value (Cooper & Showman 2006) or the planetaveraged value (Moses et al. 2011). However these estimates are crude and there is a need for theoretical work to more rigorously characterize the vertical mixing rates in hot Jupiters atmospheres.
In this study we focus in the dynamical mixing of condensable species. We therefore neglect all the potential feedback of the condensable species on the atmospheric flow such as the radiative effects of the condensates (Heng et al. 2012; DobbsDixon et al. 2012), nonequilibrium chemistry due to the depletion of a particular species, latent heat release during condensation among others. We present threedimensional general circulation model (GCM) experiments of HD 209458b to model a chemical species that condenses and settles on the nightside of the planet. We show that mixing in hot Jupiters atmospheres is dominated by largescale circulation flows resolved by the GCM. Finally, from the 3D model, we derive the values of an effective vertical mixing coefficient, representing the averaged vertical mixing in the planet atmosphere. We also compare these 3D models to an idealized 1D model parameterizing the mixing using an eddy diffusivity, with the goal of estimating an effective eddy diffusivity for the mixing rates in the 3D models. These are the first circulation models of hot Jupiters to include the influence of the dynamics on condensable species.
2 3D model
Here, we use a stateoftheart 3D circulation model, coupled to a passive tracer representing a condensable species, to determine how the interplay between dynamical mixing and vertical settling controls the spatial distribution of condensable species on hot Jupiters. Although the day/night cold trap should be present in most hot Jupiters, for concreteness, we must select a particular system to investigate. HD 209458b is among the best studied hot Jupiters. It is believed to harbor a stratosphere (Knutson et al. 2008) and a strong daynight temperature contrast (Showman et al. 2009). We decided to use it as our reference model in this study, keeping in mind that most of the mechanisms discussed here should apply to all hot Jupiters. To model the atmosphere of HD 209458b we use the 3D Substellar and Planetary Atmospheric Radiation and Circulation (SPARC/MITgcm) model of Showman et al. (2009), which couples the planeparallel, multistream radiative transfer model of Marley & McKay (1999) to the MITgcm (Adcroft et al. 2004).
2.1 Dynamics
To model the dynamics of the planet we solve the global, threedimensional primitive equations in spherical geometry using the MITgcm, a general circulation model for atmosphere and oceans developed and maintained at the Massachusetts Institute of Technology. The primitive equations are the standard equations used in stably stratified flows where the horizontal dimensions greatly exceed the vertical ones. In hot Jupiters, the horizontal scales are whereas the vertical scale height fall between and leading to an aspect ration of to . In order to minimize the constraints on the timestep by the CFL criterion, we solve the equations on the cubedsphere grid as described in Adcroft et al. (2004). The simulations do not contain any explicit viscosity nor diffusivity. However, in order to smooth the grid noise and ensure the stability of the code we use a horizontal fourthorder Shapiro filter (Shapiro 1970). In the vertical direction, no filtering process is applied. Kalnay (2003) provide a more detailed description of the equations (see Showman et al. (2009) for the numerical method used to solve them).
We use a gravity of , a planetary radius of , and a rotation rate of (implying a rotation period of 3.5 days). The average pressure ranges from 200 bars at the bottom of the atmosphere to at the top of the secondhighest level, with the uppermost level extending from a pressure of to zero. In most models, is with 53 vertical levels. In some models—particularly those with the largest cloud particle size—we adopt of bar with 47 vertical levels. In either case, this leads to a resolution of almost three levels per scale height. We use a horizontal resolution of C32, equivalent to an approximate resolution of 128 cells in longitude and 64 in latitude and a timestep of . We reran some models at C64 resolution (equivalent to an approximate resolution of in longitude and latitude) to check convergence.
It is worth mentioning that, although the primitive equations are hydrostatic in nature, it does not imply an absence of vertical motion (see section 3.5 of Holton 1992). In the primitive equations, the horizontal divergence is generally nonzero, which requires the presence of vertical motions. Indeed, the dayside radiative heating and nightside radiative cooling are balanced by a combination of horizontal and vertical thermal advection. In many cases, including Earth’s tropics, the vertical advection dominates; in such cases, vertical motions play a crucial, zerothorder role in the thermodynamic energy balance. For conditions relevant to hot Jupiters, GCMs and orderofmagnitude calculations show that vertical motions of or more are expected, despite the fact that the atmosphere is stably stratified. Further discussion can be found in Showman & Guillot (2002) and in Sect. 6.2 of Showman et al. (2008).
2.2 Radiative transfer
The radiative transport of energy is calculated with the planeparallel radiative transfer code of Marley & McKay (1999). The code was first developed for Titan’s atmosphere (McKay et al. 1989) and since then has been extensively used for the study of giant planets (Marley et al. 1996), brown dwarfs (Marley et al. 2002; Burrows et al. 1997), and hot Jupiters (e.g., Fortney et al. 2005, 2008; Showman et al. 2009). We use the opacities developed by Freedman et al. (2008), including more recent updates, and the molecular abundances described by Lodders & Fegley (2002) and Visscher et al. (2006).
As in the models of HD 209458b presented by Showman et al. (2009), our opacity tables include gasphase TiO and VO whenever temperatures are locally high enough for TiO to reside in gas phase. Note that, because of the assumption of local chemical equilibrium, we do not consider the effect of cold traps on the atmospheric composition, and thus on the opacities. TiO opacities are always taken into account where temperatures are high enough for TiO to exist in the gas phase. This causes a warm stratosphere on the dayside of our modeled planet, regardless of the 3D distribution of our tracers (to be described below).
Metallicities are solar, as given by Freedman et al. (2008). A metal enhanced atmosphere should modestly change the daynight temperature difference (at a given pressure) and the pressure of the photosphere itself, as shown in Showman et al. (2009). Cloud opacities are ignored in this study and might alter the flow more significantly (DobbsDixon & Agol 2012).
Opacities are described using the correlated method (e.g., Goody 1961). We consider 11 frequency bins for the opacities ranging from to ; within each bin, opacity information from typically 10,000 to 100,000 frequency intervals is represented statistically over 8 coefficients. This is radically different from other methods in the literature. In particular, DobbsDixon & Agol (2012) uses a multibin approach for the opacity, but inside each of their 30 frequency bins they use a single average opacity. However, inside each bin of frequency, the line by line opacity varies by several orders of magnitudes and no mean, neither the Planck mean nor the Rosseland mean opacity (see Parmentier & Guillot (2013) and Parmentier et al. (2013) for more details) can take into account this large variability. Thus, as stated in DobbsDixon & Agol (2012), our correlated radiative transfer is today’s most sophisticated radiative transfer approach implemented in a hot Jupiters GCM. Showman et al. (2009) provide a detailed description of the radiativetransfer model and its implementation in the GCM.
2.3 Tracer fields
Our target species is represented by a passive tracer field. Tracer fields are often used in GCMs to follow the concentration of a chemical species such as water vapor or cloud amount in the atmosphere or salinity in the oceans. Cooper & Showman (2006) were the first to include a passive tracer in a circulation model of hot Jupiters, in their case to investigate the quenching of CO and CH due to atmospheric mixing.
The tracer field is advected by the flow calculated in the GCM. Thus the tracer abundance is given by the continuity equation:
(1) 
Here, represents the mole fraction of the tracer, i.e., the number of molecules of the species in a given volume, either in gaseous phase or trapped in condensed particles, with respect to the total number of atmospheric molecules in that volume. For simplicity the value of is normalized to its initial value in the deep layers of the planet. In the equation, is a source term and the total derivative is defined by , where is time, is the horizontal velocity, is the horizontal gradient operator on the sphere, is the vertical velocity in pressure coordinates, and is pressure. We model the simplest possible horizontal cold trap in hot Jupiters atmospheres : a situation where the day/night temperature contrast is so strong that our target species is gaseous in the dayside of the planet but condenses and is incorporated in particles of size in the nightside. Thus the source term represents the gravitational settling of these particles and is given by
(2) 
where is the height, increasing upward, is the density of the air and is the settling velocity of the particles defined by Eq. (3). This velocity depends on the size of the particles, which is determined by the complex microphysics of condensation, out of the scope of this study (see Woitke & Helling 2003). Thus we treat as a free parameter in our model. We model spherical particles of radii , , , , , and . Equation (2) describes a simple, bimodal mechanism for the condensation of chemical species on the nightside. Although this scheme is highly simplified (ignoring the detailed temperature and pressure dependence of the condensation curves of possible condensates in hot Jupiters atmospheres), it reflects the fact that a wide range of species will reside in gaseous form on the dayside yet condensed form at low pressures on the nightside.
We consider several independent passive tracers, representing species that condenses in different particle size. Thus, they do not influence either the dynamics or the radiative transfer of the simulation and does not interact with each others. This ignores a role for possible radiative feedback mechanisms – discussed in Sect. 6 – or particle growth but represents a necessary first step toward understanding how dynamics controls the 3D distribution of a condensable species.
2.4 Settling velocity
We assume that the target chemical species condenses, in the nightside only, into spherical particles of radius that reach immediately their terminal fall speed, which is given by (Pruppacher & Klett 1978):
(3) 
where is the viscosity of the gas, is the gravitational acceleration of the planet, is the density of the particle, and the density of the atmosphere. is positive when the particles goes downward. The Cunningham slip factor, , accounts for gas kinetic effects that become relevant when the mean free path of the atmospheric molecules is bigger than the size of the falling particle. This factor has been measured experimentally by numerous experiments. We adopt the expression from Li & Wang (2003) as done by Spiegel et al. (2009).
(4) 
where the Knudsen number is the ratio of the mean free path to the size of the particle :
(5) 
For a perfect gas, the mean free path can be expressed as (Chapman & Cowling 1970):
(6) 
with the diameter of the gas molecules, and the pressure and temperature of the gas, and the Boltzmann constant. In the limit of a highdensity atmosphere, and the terminal speed becomes the Stokes velocity .
For low density gases, the dynamical viscosity is independent of pressure and can be expressed as a power law of the local temperature with an exponent varying between 1/2 in the hardsphere model to near unity, depending on the strength of the interactions between the molecules. Following Ackerman & Marley (2001), we use the analytical formula given by Rosner (2000) for the viscosity of hydrogen :
(7) 
with the molecular diameter, the molecular mass, and the depth of the LennardJones potential well (for H we use \unit2.827 ×10^10\meter and \unit59.7k_B\kelvin respectively). The power law behavior of the viscosity remains valid for temperatures ranging from to and for pressures less than (Stiel & Thodos 1963). At higher temperature, ionization of hydrogen becomes relevant and the viscosity reaches a plateau. However the temperatures of the model are everywhere less than and so Eq. (7) remains valid.
Figure 1 displays the resulting terminal velocity as a function of pressure and particle size. Two different regimes are observed. For Knudsen numbers smaller than unity, the terminal velocity is independent of pressure, whereas for Knudsen numbers exceeding unity, the terminal velocity is inversely proportional to the pressure. At low pressure and for particles bigger than a few tens of micrometers, the Reynolds number becomes higher than unity, and Eq. (3) is no longer valid. However, we show in Appendix A that these differences remain smaller than one order of magnitude and confined to a small parameter space thus we decided to neglect them for this work.
2.5 Integration time—limitation of the study
A challenge for any 3D numerical integrations of hot Jupiters dynamics is the wide range of timescales exhibited by these atmospheres. This is true for the radiative time constant, which varies significantly from low to high pressure (Iro et al. 2005; Showman et al. 2008). Nevertheless, Showman et al. (2009) showed that the computed lightcurve of HD189733b changes little for integration time longer than hundreds of days, indicating that the dynamics in the millibar regime has stabilized. As we integrate 1400 days, we consider the dynamics of the planet to be spun up at pressures less than 100 mbar. For the settling of particles, another timescale must be taken into account. The particle settling timescale can be defined as the time for the particles to fall one atmospheric scale height. To obtain a correct picture of the problem at a given location of the planet, we need to integrate the simulation for at least several times longer than the settling timescale at this location. As we can see in Fig. 2 the settling timescale ranges from tens of seconds for big particles at low pressure to tens of years for small particles at high pressure. Due to computational limitation it is not possible to run the simulation long enough to ensure that every considered tracer field has reached a statistical steady state for the full range of particle sizes we consider. The integration during 1400 days allows us to calculate the steady state for every particle size at pressure lower than and at every pressure for nightside condensates bigger than .
We define the advective timescale as the time for a parcel in the main jet stream (see Sec. 3.1) to cross one hemisphere of the planet. If the advective timescale exceeds the settling timescale at a given level, the particles at these levels will fall several scale heights while on the nightside. We thus expect that these levels will be depleted. Conversely when the advective timescale is shorter than the settling timescale, the coupling between the flow and the particle is essential and particles’ behavior cannot be predicted easily. Typical advective times in our models are 24 hours; this is marked by the upper thick black curve in Fig. 2. We expect depletion to occur above this line.
Our initial conditions for the tracer field correspond to a tracer that is spatially constant everywhere at a normalized mole fraction of 1. The final tracer abundance in our models at low pressures (less than ) is qualitatively insensitive to the initial abundance at those low pressures. This results from the short vertical mixing and settling timescales at low pressures (e.g., Fig. 2).
3 Results
3.1 Dynamical regime
We run the simulations for 1400 days and calculate the time average for all the variables over the last 400 days only, once the simulation reaches a statistical steady state at upper levels. Understanding the flow structure is essential to understanding the Lagrangian advection of particles, thus we first present a brief description of the dynamics. A more complete description of the circulation in our simulations is presented in Showman et al. (2009). The temperature structure is shown in Fig. 3. A hot stratosphere is visible at altitudes above the 10mbar level, due to the strong visiblewavelength absorption by titanium oxide present in solar abundances in the simulation. Titanium oxide abundances should be affected by the horizontal cold trap. However, as explained in Sect. 2.2, we consider its radiative effects as if it was in local chemical equilibrium, regardless of the behavior of our passive tracers. Temperatures reach at low pressures near the substellar point. By contrast, the temperatures deeper than 10 mbar in the dayside are relatively temperate. This could lead to the presence of a vertical cold trap, not considered in this study. The daynight temperature contrast becomes significant at pressures less than 100 mbar, reaching 1600 K near the top of the model. This large temperature difference results from the short radiative timescale at low pressures in comparison to dynamical timescales (Iro et al. 2005; Cooper & Showman 2005; Showman et al. 2008).
The horizontal flow on isobars comprises an eastward (superrotating) jet close to the equator and a daytonight flow pattern at higher latitude. As the pressure increases, the radiative time constant increases and the jet extends to higher latitudes. As seen in Fig. 3, the mean vertical velocities exhibit planetwide variations. The highest velocities coincide with strong horizontal convergence of the flow and occur mostly at the equator. West of the antistellar point, the convergence of the daytonight circulation forces strong downwelling motions. As described in Rauscher & Menou (2010), this convergence point appears in a range of hot Jupiters GCMs of varying complexity. It is usually associated with a shocklike feature (Heng 2012). Our simulations assume local hydrostatic equilibrium and thus can not treat properly the physics of shocks. To date, no global model of hot Jupiters atmospheric dynamics can handle shocks properly. Yet, a similar wind convergence pattern appears when considering the nonhydrostatic nature of the flow as can be seen in DobbsDixon et al. (2010). Between the substellar point and the west terminator, there are additional convergence/divergence points associated with the jet, leading to a region of strong ascending motion of longitude west of the substellar point, and a broad region of descending motion west of that. These vertical flows remain coherent over several orders of magnitude in pressure, giving them the potential to transport vertically large quantities of material. Outside of these points of strong vertical motions, the vertical velocities are more than one order of magnitude smaller. They are mostly upward on the dayside and downward on the nightside.
3.2 Spatial distribution of condensable species
Our simulations show that, as expected, condensation and particle settling on the nightside deplete the tracer from upper levels relative to the abundances at depth—an effect that is stronger for larger particles. This is illustrated in Fig. 4 (solid curves), which shows the globalmean tracer abundance (averaged horizontally on isobars) versus pressure for simulations with nightside condensates sizes ranging from to . In all cases, the horizontally averaged tracer abundance decreases with altitude. The depletion is modest for the smallest particle size (), but for the largest particle sizes, tracer abundances at the top are almost two orders of magnitude smaller than abundances at the bottom. A useful metric is the “ depletion pressure,” that is, the pressure above which the globalmean tracer abundance is less than of the deep abundance. This pressure is only for a particle size of but is for a particle size of . We also note that, for all the models shown in Fig. 4, the particles settling times near the top of the model are much less than our integration times; at low pressures, the tracer abundances have reached a statistical equilibrium where downward transport of tracer due to particle settling is balanced by upward mixing of tracer by the largescale dynamics.
The tracer abundance on isobars exhibits a strong spatial variation as can be seen in Fig. 6. Although a day/night pattern is imposed in the tracer source/sink (with particle settling on the nightside but not the dayside), the threedimensional advection of the tracer field by the atmospheric winds leads to a complex tracer distribution that does not exhibit an obvious daynight geometry. The main pattern appears to be an equatortopole gradient, with large zonalmean abundances at the poles, and smaller zonalmean abundances at the equator. This is particularly true around . Significant longitudinal tracer variations also occur; at 1 mbar (Fig. 6), these variations are particularly prominent at high latitudes. Interestingly, these variations are phase shifted in longitude relative to the daynight pattern with maximum (minimum) peak tracer abundances occurring 60–80 of longitude east of the substellar (antistellar) point.
On top of these main patterns, we clearly see two points depleted in tracers at the equator. These two points correspond to points of horizontal convergence of the flow and high downwelling motions as discussed in the previous section. This correlation between strong downwelling motions and low tracer abundances arises naturally in the presence of a background vertical gradient of tracer abundances. Due to their settling on the nightside, the local tracer abundance generally decreases with height and thus any downwelling motion would carry parcels of gas depleted in tracers whereas any upwelling motion should carry parcels of gas enhanced in tracers.
3.3 Geometry of the mixing
As discussed in the previous section, the tracer abundance is not null everywhere in the planet, which implies that tracers did not rain out during the simulation time. Yet the integration time greatly exceeds the fall times at low pressure for all particle sizes considered, and everywhere throughout the domain for particle sizes exceeding a few . Therefore vertical mixing must happen in order to keep these particles aloft. This vertical mixing is characterized by an upward dynamical flux of tracers that balances the downward flux due to the gravitational settling in the nightside. The upward dynamical flux of tracers across isobars can be calculated as where is the vertical velocity in pressure coordinates, the tracer abundance, the gravity of the planet and the brackets denote the average on isobars. As seen in Fig. 5, the upward flux of tracers due to the dynamics (solid lines) balances nicely the downward flux due to settling (dashed lines), showing that the simulation did reached a quasi steadystate.
Hot Jupiters atmospheres are heated from above and thus are believed to be stably stratified. Then, vertical mixing cannot be driven by small scale convection as it is the case in the deep atmosphere of brown dwarfs (Freytag et al. 2010) and the gas giants of the solar system. As explained in Sect. 2.1, our model does not use any parametrization of subgrid scale mixing. Thus we do not account for mixing induced by smallscale turbulence and gravity wave breaking, two mechanisms that are believed to dominate the mixing in the radiative part of brown dwarfs atmospheres (Freytag et al. 2010). Rather, the upward flux of tracers in our model is due to the largescale, resolved flow of the simulation.
Given that mass is conserved, any upward flux of gas is compensated by a downward flux of gas. Thus, if the tracers concentrations were horizontally homogeneous on isobars, there would be no net upward flux of tracer through that isobar. For a net upward flux of tracers across isobars to occur, there must be a correlation between the horizontal distribution of the tracers and the vertical velocities. Dynamics will produce an upward flux if—on an isobar—ascending regions exhibit greater tracer abundance than descending regions. In other words, an upward tracer flux due to dynamics will occur only if where is the upward velocity in pressure coordinates, is the tracer abundance and the brackets are the mean over one isobar (note that negative implies upward motion). Given a vertical gradient of such that the abundance of decreases upward, an upward flux of gas will naturally bring enhanced material whereas a downward flux will naturally advect parcels of gas depleted in tracers, thereby creating the correlation between and favorable for upward tracer transport. Figure 7 shows the relative contribution to the upward mixing versus longitude and latitude at a given isobar, which can be estimated by the quantity :
(8) 
Mass conservation in the primitive equations implies that : the advection of does not contribute to the net (horizontally averaged) upward flux of material and so we remove the contribution of this term when defining in Eq. (8). The quantity is the mean upward flux of material across isobars, thus the quantity represents the local contribution to the total upward flux on isobars. It is normalized such that .
The strength of the mixing varies significantly with longitude and latitude. Both upward and downward fluxes are one order of magnitude greater than typical values in a handful of specific small areas across the planet—particularly at the two points of horizontal convergence and strong vertical velocities described in Sect. 3.1. This vertical flow remains coherent over several order of magnitude in pressure (see Fig. 3), acting like a vertical “chimney” where efficient transport of material can be achieved.
In summary, the mechanism by which the largescale, resolved atmospheric circulation transports tracers upward is extremely simple and straightforward. The settling of particles leads to a mean vertical gradient of tracer abundance, with, on average, small tracer mixing ratios aloft and large tracer mixing ratios at depth. Given this background gradient, advection by vertical atmospheric motions—whatever their geometry—autmomatically produces a correlation between and () on isobars, with ascending regions exhibiting larger values of than descending regions. In turn, this correlation automatically causes an upward dynamical net flux of tracers when averaged globally over isobars. In statistical steady state, this upward dynamical flux balances the downward transport due to particle settling and allows the atmospheric tracer abundance to equilibrate at finite (nonzero) values despite the effect of particle settling. The mechanism does not require convection, and indeed, the vertical motions that cause the upward transport in our models are resolved, largescale motions in the stably stratified atmosphere. These vertical motions are a key aspect of the globalscale atmospheric circulation driven by the daynight heating contrast.
3.4 Time variability
Besides the spatial variability at a given time, the model exhibits significant temporal variability, both in the 3D flow and especially in the tracer field. The equatorial jet exhibits an oscillation pattern at planetary scale as seen in Fig. 9. At the convergence point west of the substellar point, the jet orientation can be toward the north, the south or well centered on the equator. Whereas DobbsDixon et al. (2010) described a variation in longitude with time of the convergence point of the flow in the nightside, we see a variation in latitude of this convergence point and interpret it as a result of the larger oscillation of the jet itself.
The tracers abundances at specific locations on the planet exhibit strong temporal variability. This is illustrated in Fig. 9, which shows the tracer abundance at and over the globe at several snapshots in time for a model where the radius of particles on the nightside is m. Significant variations in tracer abundance are advected by the equatorial jet and, at high latitudes, by the daytonight flow, leading to large local variations in time. In many cases, the strongest tracer variability seems to involve regionalscale structures with typical sizes of 1–3km but also includes hemisphericscale fluctuations (e.g., in the abundance averaged over the day or night) and between the northern and southern hemispheres. Around the substellar point, the tracer abundance can vary by up to , whereas along the terminator, this temporal variation can reach relative to the mean value. Such variability—if it occurs in radiatively active species like TiO—has important implications for secondaryeclipse and transit observations, which probe the dayside and terminator, respectively.
Figure 8 sheds light on the different timescales at which this variability occurs. The top panel shows the tracer abundance averaged vertically between and and horizontally over a circular patch of centered on the substellar point; this gives a sense of how the tracer abundance would vary in secondaryeclipse measurements probing the dayside. The bottom panel shows the tracer abundance averaged vertically between 0.1–1 mbar and horizontally around the terminator, including all regions within of the terminator itself; this gives a sense of how the tracer might vary in transit measurements. The variability exhibits two characteristic timescales: a short (fast) timescale of order of days and a long (slower) timescale of 50 to 100 days. The bigger the particles on the nightside, the bigger the amplitude of the variations. This comes from the smaller settling timescale of bigger particles. In these models, the amplitude of the longperiod variations exceeds those of the shortperiod variations by a factor of 2–3. The longtimescale variations exhibit similar amplitudes in the dayside and terminator time series. The shortperiod oscillations exhibit stronger amplitude at the terminator and seem related to variations of the flow itself, such as the oscillation of the jet described previously. Figure 8 suggests that radiatively active tracer species that can condense on the nightside, such as TiO or silicates, could lead to detectable time variations in transit or secondary eclipse spectra. The amplitude of this variability will depend on the type of tracer being considered (see Sect. 6) and may vary from planet to planet depending on the availability of the considered species. However we can predict the expected period of these variations: some days for the small amplitude ones and fifty to one hundred days for the biggest ones.
3.5 Limb profile
Transit observations are sensitive to atmospheric composition near the terminator, thus we want to characterize the distribution of our tracer species at the terminator. The possibility of variations in chemical composition between the leading and trailing limbs (as seen during transit) has been discussed in a variety of studies (e.g., Iro et al. 2005; Fortney et al. 2010). However, these studies did not investigate the particular depletion of species due to the interaction between their condensation and the atmospheric dynamics. Our model leads to the first quantitative estimate of how dynamics affects the spatial distribution of condensable species at the daynight terminators, relevant to the interpretation of transit observations. Figure 11 shows the tracer abundance at the terminator at a snapshot in time for our models with particle sizes of 0.1, 0.5, 1, 2.5, 5, and m. Angle represents angle around the terminator and the radial coordinate represents the logarithm of the pressure. In agreement with Fig. 4, the tracers tend to be depleted from upper levels, particularly in models where the particles on the nightside are larger. Moreover, Fig. 11 demonstrates that significant spatial variations occur along the terminator. Depletion occurs first at the equator along the leading limb, corresponding to the terminator west of the substellar point. The superrotating jet carries air depleted in tracer from the nightside directly to this region of the terminator, explaining why abundances are particularly depleted there. In contrast, air along most of the remainder of the terminator has arrived from the dayside, where no particle settling occurs, so depletion is less strong—particularly for particle sizes m. Once the particle size becomes sufficiently large, however, depletion occurs everywhere along the terminator at upper levels regardless of whether the air arrived there from the dayside or nightside.
Considering now the depth dependence of the terminator abundances, our results suggest two different zones (see Figs. 10 and 11) :

At altitudes above the 1mbar level, the tracer abundance is homogeneous over most of the limb except for the east (trailing) equatorial limb that is strongly depleted.

At altitudes below the 1mbar level, the east and west equatorial limb are rather homogeneous; however, the east/west dichotomy shifts to higher latitudes and the west limb above is more depleted than the equivalent region of the east limb.
Moreover, it should be noted that due to the shift of the hot spot, the temperatures at the east limb are higher than at the west limb, thus a given species is more likely to be gaseous and detectable on the trailing limb than on the leading limb.
4 1D model of the daynight cold trap
Although hot Jupiters atmospheres are inherently threedimensional, 1D models continue to play a useful role for understanding the vertical thermal and chemical structure of these atmospheres. In particular, many groups have explored the chemistry of hot Jupiters using 1D models in which the vertical mixing caused by the largescale dynamics is parameterized by a specified eddy diffusivity (e.g., Spiegel et al. 2009; Zahnle et al. 2009a, b; Youdin & Mitchell 2010; Line et al. 2010, 2011; Madhusudhan & Seager 2011). In these studies, the chosen eddy diffusivity is ad hoc, with no convincing theoretical support. Although the vertical mixing in our 3D models is not diffusive in any rigorous sense, there is merit in comparing the results of our 3D models with 1D models parameterized by eddy diffusivity. This allow us to make approximate estimates of the magnitudes of eddy diffusivity—in the context of a 1D model—that produce similar horizontalmean behavior as our 3D models. Such estimates of eddy diffusivity should guide the parameter choices in 1D chemical models like those cited above. A comparison between our 3D models and 1D diffusive models also allow us to investigate how the horizontalmean tracer depletion relates to the amplitudes of spatial tracer variation.
Therefore, in this section, we present a simple 1D model, including particle settling, with atmospheric mixing represented as an eddy diffusivity.
4.1 System studied
The presence of a superrotating, eastward equatorial jet is a dominant dynamical feature of many 3D circulation models of hot Jupiters. This superrotating jet was first predicted by Showman & Guillot (2002), and later emerges from almost all 3D simulations of hot Jupiters atmospheres (Cooper & Showman 2005; Showman et al. 2008, 2009, 2013; DobbsDixon & Lin 2008; Rauscher & Menou 2010, 2012b, 2012a; Perna et al. 2010, 2012; Heng et al. 2011a, b; Lewis et al. 2010; Kataria et al. 2013) including ours (see Sect. 3.1) and has been theoretically understood (Showman & Polvani 2011). A shift of the hottest point of the planet eastward from the substellar point has been directly observed in several exoplanets (Knutson et al. 2007, 2009, 2012; Crossfield et al. 2010) and interpreted as a direct consequence of this jet. Thus we believe that any study of the day/night cold trap in hot Jupiters atmospheres must account for this feature.
To include the presence of this jet in our model, we choose as a study system a vertical column of gas homogeneously advected around the equator by the superrotating jet. Such a column is transported from day to night and from night to day with a period where is the advective timescale, is the planetary radius and the equatorial jet velocity. is around 48h for HD 209458b.
As in the 3D case, we focus on a hypothetical chemical species which is gaseous on the dayside and trapped in condensates of size on the nightside. This species freely diffuses with a vertical diffusion coefficient on both the dayside and the nightside. Because we envision the species as condensed on the nightside, we additionally include downward settling via StokesCunningham drift on the nightside (Eq. (3)) but not on the dayside. The model includes no horizontal dimensions. Rather, we model a single column of gas that is advected from day to night. The horizontal variations therefore translate into a timedependant settling term. We assume this chemical species to be a minor constituent of a Hatmosphere and neglect the latent heat released during the condensation.
4.2 1D diffusion equation
As before, is the local mole fraction of the target chemical species i.e.the number of moles of tracer species (whether in gaseous or condensed form) to the total moles of air in a given volume. On the dayside, the molecules of the target chemical species can freely diffuse with a diffusion coefficient , according to the equation:
(9) 
with the density of the atmosphere and the vertical coordinate.
On the nightside, the molecules of the target species are trapped into particles that both diffuse and settle with the velocity described in 2.4. Thus follows the same equation as Eq. (9), plus a source term describing the settling:
(10) 
We can define the diffusive time scale as and a reference free fall time scale with the atmospheric scale height and the Stokes velocity (see 2.4). We note that is a reference time scale and is not equal to the effective free fall time scale for high Knudsen numbers. Assuming hydrostatic balance, we can use the pressure as the vertical coordinate and using the perfect gas law the system to solve become:
(11) 
4.3 Timedependent solution
In order to solve the system of equations (11) for , we need two boundary conditions. We assume the species to be well mixed below with a molecular abundance . At the top of the atmosphere, we assume that no molecule crosses the upper boundary (i.e.). Then we can solve the system with an implicit time stepping code using GNU OCTAVE^{1}^{1}1http://www.octave.org, an opensource, free software equivalent to MATLAB. Assuming that the column of gas spends in each hemisphere, we can reach a periodic behavior where the initial condition is forgotten. While on the dayside, there is no settling and, at upper levels, the tracer diffuses upward. Thus at a given pressure the molecular abundance increases with time. This is shown by the red curves in Fig. 12. While on the nightside, the particles both diffuse vertically and settle downward at their terminal velocity and thus, at at given pressure level, the abundance decreases with time. This is shown by the blue curves in Fig. 12. In the upper atmosphere, the settling timescale and the diffusion timescale become smaller than the advective timescale. The particles have time to settle several scale heights during the time required for the air column to cross the dayside. Thus, at upper levels, the molecular abundance vary strongly throughout the diurnal cycle, and the vertical profiles vary widely around the mean value.
4.4 Steadystate solutions
Our system being forced periodically, there is no steadystate solution stricto sensu. However, the mean of the tracer abundance over one period should remain constant. We can thus integrate eqs. (11) over one period and get an equation for , the mean molecular abundance over one period. Then, once the periodic state is reached, the mean over one period of the tracer abundance in our column of gas, , is the same as the mean over longitude at a given time, and we obtain :
(12) 
where the factor appears because the source term is only integrated during the night. To derive Eq. (12) we assumed that the mean tracer abundance on the nightside is close to the globalmean tracer abundance. This is true given than and are much bigger than . This approximation breaks down at low pressure and our analytical solution diverges from the real one as can be seen by comparing the black and the green curves in Fig. 12. However, the discrepancy is small such that the analytic solution is still a good representation of the timemean of the full numerical solution of the 1D model everywhere deeper than the level. is a constant coming from the integration over pressure. When goes to , goes to . Moreover we assume that does not go to infinity. Then the constant must be zero and we obtain :
(13) 
To simplify the problem, we neglect the transitional regime for from Eq. (4) and use the expression :
(14) 
Then choosing a functional form , where is a constant, we can solve equation (13). For and , we get :
(15) 
For the particular case of a constant () the formula becomes :
(16) 
In the case where is inversely proportional to ():
(17) 
Where is the abundance of at and .
4.5 Diffusivities needed to keep the tracers suspended
Because particle settling acts to transport condensates downward, the tracers only exhibit significant abundances in the upper regions of the atmosphere if the eddy diffusion coefficient exceeds some critical value, which depends on the particle size of the condensates. Here we solve for an approximate analytical expression for this critical magnitude as a function of particle size and other parameters.
Assuming a constant vertical diffusion coefficient, we can use Eq. (16) to derive an expression for the needed to achieve a given molecular abundance at a given pressure . Like in Spiegel et al. (2009), we use and . We first note that Eq. (16) is composed of two terms. The first exponential is given by the Stokes regime () whereas the second term is given by the Cuningham regime (). As can be seen in Fig. 1, particles smaller than should be in the Cunhingham regime at . Thus the second term in Eq. (16) is dominant and we can use as a condition :
(18) 
As in the 3D model, we assume that below the tracers are no more trapped into condensates and no more subject to settling. Thus we choose at . Then, for , replacing and by their expressions, using Eq. 7 for the viscosity and Eq. 6 for the mean free path, we obtain a condition on the diffusion coefficient:
(19) 
For big particles, the Stokes regime become dominant thus we neglect the second term of Eq. (16) and the condition turns to be:
(20) 
As the relevant range of particle size span several order of magnitudes, a good approximation of can be obtained by taking the sum of these two coefficient:
(21) 
Again we note that these limits for the diffusion coefficient are independent of the planet considered. Spiegel et al. (2009) stated that an abundance of half the solar composition at 1mbar would be necessary to produce an observable stratosphere. Applying formula (21) with and and assuming a well mixed layer below (i.e. and ) we obtain for , for and for . These results are of the same order of magnitude as the ones found by Spiegel et al. (2009) for the vertical cold trap. This is expected, since we compare two similar mechanisms: settling and diffusion of particles. However, this similarity of the results shows that the daynight cold trap is at least as important as the vertical one in hot Jupiters atmospheres. Moreover, the condition on we derived is independent of the planet studied and holds for very hot Jupiters such as WASP12b or WASP33b where the vertical cold trap could be inefficient or nonexistent.
5 Effective vertical diffusion coefficient
As described previously, 1D models have been extensively used to investigate chemistry and vertical structure of hot Jupiter atmospheres, with the 3D dynamics parameterized as vertical eddy diffusion with a specified diffusivity (e.g Spiegel et al. 2009; Zahnle et al. 2009b, a; Youdin & Mitchell 2010; Line et al. 2010, 2011; Madhusudhan et al. 2011; Moses et al. 2011; Venot et al. 2012).
There is no theoretical reason for the mixing by the large scale flow patterns in hot Jupiters atmospheres to behave like a onedimensional diffusion process. However, deriving an a posteriori effective diffusion coefficient that describes as closely as possible the averaged vertical mixing within the atmosphere can be a useful way to roughly characterize the strength of the vertical fluxes of material and guide 1D modelers in their choice of vertical mixing parameters.
A first way to define a vertical mixing coefficient from our simulation is to choose the that best reproduces the planet averaged tracer profiles. The 1D model developed in Sect. 4 describes the equilibrium between vertical diffusion of tracers with a specified heightdependent and their settling on the nightside. We tune the diffusivities to obtain a good match between the solutions of our 1D model and the horizontalmean tracer abundance versus pressure from the 3D models. To use Eq. (15) we must specify the temperature of the atmosphere, constant in the analytical model. The temperature appears in the expression for the Knudsen number and in the expression for the viscosity of hydrogen. Both quantities are related to the settling of the particles. Thus the temperature to consider is the nightside temperature. In our GCM nightside temperatures range from to and we decided to use a mean temperature of 1000 K. However, we note that the derived value of does not depend strongly on this choice. Then using a value proportional to the inverse square root of the pressure ( in Eq. (15)), we obtain a remarkably good agreement between our 1D model (dotted lines in Fig. 4) and the horizontal average of the 3D model (solid lines of Fig. 4). The resulting value for the vertical mixing coefficient that best fits the different tracer field used in the simulation is :
(22) 
which is valid over a pressure range from 1 bar to a few bar (see Fig. 13). Fundamentally, this represents an adhoc fitting of our simulation results, although we find it a good match to the overall mixing properties of our 3D models over a wide range of particle sizes. Note that the upwardincreasing mixing rates captured in Eq. (22) arise naturally from the fact that the radiative heating rates and vertical velocities tend to increase with decreasing pressure in our models, leading to greater mixing rates at lower pressure. The 3D models adopt the stellar insolation and other properties for HD 209458b, so the results are most germane to that planet; the mixing rates are likely to be higher for hotter planets and lower for cooler planets than implied by Eq. (22).
Another way to define a one dimensional vertical mixing coefficient from the threedimensional simulation is to find the that leads to an upward diffusive flux of material that matches the averaged vertical flux produced by the dynamics. This can be written (e.g., Chamberlain & Hunten 1987, p. 90):
(23) 
where the brackets represent the horizontal average along isobars over the whole planet. This expression does not necessitate any assumption on the functional form of nor on the nightside mean temperature. It also has the advantage that no comparisons or fits to a 1D diffusion model are necessary; the effective values of can be derived directly from the 3D GCM data via Eq. (23). As a tradeoff, this expression depends on the vertical tracer gradients that are affected by long term temporal variability (see Sect. 3.4), which are not smoothed out completely given the limited integration time of the simulation. This leads to some strong vertical variation of . The profile of calculated from Eq. (23) are shown in the black curves of Fig. LABEL:fig:Kzz. Despite the vertical fluctuations, the overall shape of obtained with this method is close to the estimate using the planet averaged tracer profile (see the red curve in Fig. 13). Like before, we note that the derived value does not depend strongly on the particle size, consistent with the fact that Eq. (22) leads to a good fit between 1D models and our 3D model for a wide range of particle sizes.
Several previous studies have attempted to estimate the vertical diffusion coefficient in hot Jupiters atmospheres (Cooper & Showman 2006; Moses et al. 2011; Heng et al. 2011a; Lewis et al. 2010). Cooper & Showman (2006) adopted an estimate for based on the product of a rootmeansquare vertical velocity from their 3D GCMs and an appropriately chosen vertical length scale following the formulation of Smith (1998). Moses et al. (2011) and Lewis et al. (2010) followed a similar procedure but adopted an atmospheric scale height for the vertical length scale. These estimates are crude, although Cooper & Showman (2006) showed that this formulation for allows 1D models to match the full tracer profiles from 3D GCMs reasonably well. More recently, Heng et al. (2011a) used the magnitude of the Eulerian mean streamfunction as a proxy for the strength of the vertical motions and derived a vertical mixing coefficient of . Note, however, that the Eulerianmean velocities are known to be a poor descriptor of tracer advection rates in planetary atmospheres, since mixing by large scale eddies (potentially resolved by GCMs) often dominates over transport due to the Eulerianmean circulation (see, e.g., Andrews et al. 1987, Chapter 9).
Although we confirm that mixing in hot Jupiters atmospheres is strong, we find a value that is significantly smaller than the previous ones. In particular, our value is two orders of magnitude smaller than what is obtained when multiplying the vertical scale height by the root mean square of the vertical velocity (blue curve in Fig. 13), a common estimate for in the literature (e.g., Lewis et al. 2010; Moses et al. 2011).
As stated in Sect. 2.1, the model does not include any subgrid vertical diffusion coefficient. Yet, given the large values for that we derive from the resolved flow, it seems unlikely that subgrid turbulent mixing would contribute significantly to the total mixing. However, the interaction between smallscale turbulence and the global flow might not be trivial and a more detailed study would be required to draw a firm conclusion.
We emphasize that vertical mixing by the global circulation appears to be planetwide and differs from region to region. Although the globally averaged dynamics seem to be reasonably described by a vertical mixing coefficient, that is not the case for the local flow in the simulation. It is therefore difficult to define mixing coefficient values for particular locations in the planet.
Along the equator, where the strong flow efficiently mixes the tracers longitudinally, we expect a good agreement of our 1D model to the 3D flow. Indeed, using the value of derived in Sect. 5 we realize that the spread of the tracer profiles along the equator (Fig. 14) is of the same order of magnitude as the spread predicted by the 1D model (middle panel of Fig. 12). However, in the 1D model, profiles from equally sampled longitude are equally spaced in abundances whereas the profiles obtained from the 3D model are sometimes packed together and sometimes widely spread, denoting an unequal strength of the vertical mixing longitudinally.
6 Applications
6.1 Presence of a stratosphere on hot Jupiters
TiO is a leading hypothesis for the absorber needed to create temperature inversions in hot Jupiters atmospheres (Hubeny et al. 2003; Fortney et al. 2008). However, TiO condenses at temperatures lower than . On most hot Jupiters, while gaseous TiO can be stable on the dayside, it should condense on the nightside. Thus, our results allow us to address the question of whether TiO can remain suspended in the atmospheres of hot Jupiters, and hence whether TiOinduced stratospheres are indeed viable. Spiegel et al. (2009) predicted that for a solar abundance of TiO in the planet, an abundance of 0.5 times the deep abundance at 1 mbar is necessary to maintain a temperature inversion on the dayside. Our simulations suggest that if TiO condenses into particles bigger than several micrometers, the daynight cold trap will be sufficiently efficient to deplete it from the dayside. If, on the contrary, TiO cannot condense in particles bigger than several micrometers, it should remain present on the dayside and produce a stratosphere (see Figs. 4 and 6).
The size of the condensate, a free parameter in our study, results from complex microphysical processes. Once on the nightside, TiO is oversaturated, and thus Tibearing condensates are expected to appear. We used the formalism of Woitke & Helling (2003) to calculate the characteristic growth time scale of particles assuming that all the titanium is contained in and that this last fully saturates the atmosphere. If the condensate growth time scale exceeds the advective time scale – the time for the jet to travel across one hemisphere – the condensate will be back in the dayside before reaching it’s full size and will vaporize. Thus, all the particles above the black line in Fig. 15 are unlikely to form. The low elemental abundance of titanium—solar abundance is compared to H (Lodders 2002)—kinetically inhibits the formation of micrometer size particles above .
However, titanium is not the only element that can form condensates on the nightside of hot Jupiters. Silicates are believed to condense and could incorporate titanium atoms into their grains. Submicron size particles could even be used as seeds for the formation of silicates grains. In that case the relevant time scale is not tied to the growth of grains but rather to the growth of based grains ( for example). We calculated the growth time scale of based grains using the same formalism as for grains, assuming that is fully saturated and that all the silica atoms are in molecules. As can be seen in Fig. 16, particles as big as can form at pressures as low as in that case.
The estimates in Fig. 16 show that the growth time of silicate particles is comparable to the time for the jet to cross an hemisphere in the pressure range where the stratosphere forms (0.01–1mbar). This comparison suggests that more detailed calculations, coupling the 3D dynamics to microphysics that allow selfconsistent prediction of particle growth, may be necessary to obtain a firm conclusion about whether particle growth timescales are sufficiently long to inhibit loss of TiO from the atmosphere.
Spiegel et al. (2009) studied the deep cold trap in the deep layer of the planet, at pressures exceeding tens of bars, where dynamical mixing rates are probably low. In planets that exhibit such a cold trap, mixing TiO upward to altitudes where it could be affected by the daynight cold trap may be difficult. On the other hand, the opacities and therefore temperature structure in these deep regions are rather uncertain; moreover, hot Jupiters that are particularly highly irradiated would exhibit warmer temperatures and would therefore be less likely to exhibit such a vertical cold trap. On these planets, the daynight cold trap at low pressure (bar) would then dominate.
The spatial variations in the gaseous tracer abundance in our 3D models (Fig. 6) suggest that the TiO abundance on the dayside, and hence the stratosphere itself, could be patchy, with some regions of the dayside exhibiting a stronger temperature inversion than others. This would have interesting consequences for the interpretation of dayside infrared spectra.
The time variability described in Sect. 3.4 could affect the presence of the stratosphere, leading to strong temporal variability in the upper atmospheric temperatures. In particular, the tracer abundance averaged over much of the dayside exhibits largeamplitude fluctuations (Fig. 8), particularly for larger particle sizes. This suggests that, at least for some planets, the TiO abundance could fluctuate between values large enough to generate a stratosphere and values too small for a stratosphere to form. The stratosphere itself might then fluctuate episodically in time, leading to variations of a factor two in the thermal flux emitted by the planet at some wavelength (see Fig. 12 of Fortney et al. 2008). Although not included in our current models, there is the possibility of feedbacks with the flow itself, since the presence (or absence) of a stratosphere exerts a significant impact on the flow structure and vertical mixing rates. If the feedback is positive, i.e. if the presence of high temperatures in the upper atmosphere enhances the mixing, then hot Jupiters could oscillate between a state with strong vertical mixing and stratospheric heating by TiO and a state with no stratospheric heating and less vertical mixing. This is a twostate atmosphere analogous to that described by Hubeny et al. (2003). However, this possibility remains speculative and further models that include the feedback of the tracer field on the flow are necessary to draw a firm conclusion.
TiO is thought to be the major Tibearing gas in hot Jupiters atmospheres (see Lodders 2002). However, other Tibearing gases, being the most abundant, are believed to be the condensable Tibearing species. Thus, a parcel of gas experiencing a sudden drop in temperature due to its advection to the nightside might not see all its titanium incorporated into condensates, but rather only the titanium atoms that already reside in other Tibearing molecules such as . The reaction between gaseous species is fast under the conditions relevant to the dayside of hot Jupiters (see Fortney et al. 2008). However, the reaction might be kinetically inhibited on the nightside of the planet where the temperature drops significantly, leading to a smaller depletion of Ti than in the hypothetical case where could condense by itself.
6.2 Clouds in hot Jupiters atmospheres
The huge daynight temperature difference and cold nightside temperatures predicted on many hot Jupiters at low pressure (e.g., Showman et al. 2009) suggest that, in addition to TiO, a wide range of other chemical species, including silicates and iron, will condense on the nightside. Some of them could also stay in a condensed state in part or all of the dayside hemisphere. The Rayleigh scattering slope in the transmission spectrum of HD189733b, first observed by Lecavelier Des Etangs et al. (2008) and later confirmed by numerous observations (see Pont et al. (2013) for a review of the different observations of this planet) is best fitted by models including submicron sized particles. Then the strong spatial and temporal variations observed in our model can also be interpreted as spatial variation of the cloud coverage in the atmospheres of hot Jupiters. This could lead to albedo variations along the dayside of the planet. As the hottest point of the planet is shifted to the east, the western part of the dayside is colder than the eastern one (see Fig. 3). Therefore, the western part of the dayside could be colder than the condensation temperature of some species whereas the eastern part could be at higher temperatures, leading to a more cloudy atmosphere west of the substellar point than east of it. Moreover, due to the eastward superrotating jet, material west of the substellar point arrives from the cold nightside – where condensation is thought to happen — whereas material east of the substellar point arrives from the hot substellar point — where the material is thought to have sublimated. Thus, clouds, if present on the dayside of the planet, should form more easily west of the substellar point than east of it, leading to a longitudinal variation of the albedo that contributes to the spatial variation described before. The large amount of data in the visible from the Kepler space telescope is ideal to search for such a spatial and temporal variability in albedo pattern of tidally locked planets.
6.3 Parameter retrieval
Atmospheric characteristics of hot Jupiters are usually derived from diskintegrated fluxes (for secondary eclipses measurements) or limbintegrated transmission (for transits spectroscopy). Therefore, interpretation of the data is usually done using one dimensional atmospheric models, assuming an homogenous atmosphere, both in term of temperature structure and composition (e.g., Madhusudhan & Seager 2009; Lee et al. 2012; Benneke & Seager 2012, among others). However, given the strong spatial variability in the tracer distribution both on the dayside hemisphere and in the limb profiles observed, some future exoplanet spectra — obtained with more and better data than nowadays — might be better understood by considering spatial variation in the atmospheric profiles and chemical composition along the planet (see also Agúndez et al. 2012, on the longitudinal variability in the chemical composition of HD209458b). For example, an inhomogeneous distribution of TiO would lead to strong brightness differences in the emitted flux from the different locations of the planet, which could affect features such as the apparent eastward offset of the brightest spot of the planet. Hazes could also show a similar behavior, leading to planetary spectra that might be better explained with the combination of two different (cloudy and cloudless) 1D models—as is the case in recent models of brown dwarfs (e.g., Marley et al. 2010; Burrows et al. 2011). The rising technique of secondary eclipse mapping (see Majeau et al. 2012; De Wit et al. 2012), combined with the accuracy of telescopes like EChO (Tinetti et al. 2012) or JWST might soon allow us to constrain better the spatial inhomogeneities in the disk of the planet.
Conclusion
We presented global, threedimensional numerical simulations of the atmospheric circulation of HD 209458b including a passive tracer. This is the first circulation model of a hot Jupiter to include the dynamical mixing of condensable species. We applied our model to chemical species that are gaseous on the dayside but condense on the nightside of the planet and are trapped in particles of a given size. Given the strong day/night contrast present in hot Jupiters, our model applies to a wealth of different chemical species such as titanium, vanadium, and silicate oxides among others.
Prior studies of hot Jupiters circulation demonstrate the presence of 3D circulation patterns including both strong horizontal and vertical velocities that are necessary for mixing to happen (e.g., Showman & Guillot 2002; Cooper & Showman 2005, 2006; Showman et al. 2008, 2013; Rauscher & Menou 2012b, a, 2013; Heng et al. 2011a, b, 2012; Lewis et al. 2010; DobbsDixon & Lin 2008; DobbsDixon et al. 2010, 2012; Perna et al. 2010, 2012). While preliminary attempts have been made to quantify the mixing rates (Cooper & Showman 2006; Heng et al. 2011b), these methods are approximate and do not lead to a rigorous quantification of the mixing rate. Here we demonstrate that, although hot Jupiters atmospheres are believed to be stably stratified (i.e. locally nonconvective), they are strongly mixed. In the presence of a background gradient of chemical species, largescale circulation patterns naturally create upward mixing. This mixing, resolved by the GCM, is strong and likely dominates over molecular, convective or turbulent mixing. In HD 209458b, the mixing is strong enough to keep a condensable species aloft if it condenses into particles smaller than a few microns on the nightside of the planet.
The coupling between 3D flow and particle settling leads to strong spatial and temporal variations in the abundance of a given condensable species. Around , the tracer abundance is homogeneous in longitude but exhibits a large latitudinal variation, the equator being more depleted than the poles. Around , at high latitudes, the day/night contrast becomes important. According to our models, variability of up to 50% in the dayside tracer abundance, and of up to 75% in the tracer abundance along the limb, can occur for sufficiently large particle sizes (m). This variability characteristic periods ranging from days to 50–100 days. The observability of such a variation depends on the radiative properties of the considered species and will be quantified in a future work.
These results can be applied to a wide range of molecules in hot Jupiters atmospheres. Titanium oxide, the best candidate for creating a temperature inversion in the dayside of hot Jupiters, should condense on the nightside of most planets. Our results imply that the day/night cold trap could impede the formation of a stratosphere in the dayside if TiO condenses into particles bigger than a few microns on the nightside. Growing particles to such a size seems difficult when TiO alone is considered due to its small abundance. However, TiO can be incorporated into condensates from more abundant gases such as silicate oxides. In that case the day/night cold trap could be strong enough to impede the formation of a hot stratosphere on the dayside. Spatial variability of TiO could significantly affect the dayside temperature structure and exert interesting effects on infrared spectra and lightcurves. For example, midtohigh latitudes might keep enough TiO to create an inversion whereas the equator, more depleted, might not be able to sustain the inversion. Such a latitudinal contrast could be observed using the secondary eclipse mapping technique (e.g., Majeau et al. 2012; De Wit et al. 2012). The temporal variability observed in the model could lead to the appearance and disappearance of the stratosphere on timescales of 10–100 earth days. Highly irradiated planets have significant thermal emission in the Kepler bandpass (e.g., Spiegel & Burrows 2010, for HatP7b). Using the long photometric series from the Kepler Space Telescope, such a variability might be observable.
Our results also apply to silicate hazes. Temporal and spatial variability in the cloud coverage could strongly affect the albedo and the thermal emission of the planet. For moderately irradiated planets, the Kepler spacecraft observes the reflected light of the star by the planet (e.g., Demory et al. 2011, for Kepler7b). Thus, the time series from Kepler could be used to build albedo, and therefore cloud maps of the planet.
Although there is no theoretical reason for the upward mixing driven by the global circulation to be diffusive, it is interesting to quantify the averaged vertical mixing with a diffusive model. The parameterization value of or , between and , can be used in 1D models of HD 209458b.
This study, the first one to include the influence of the dynamics on condensable species in a GCM of a hot Jupiter, confirms that hot Jupiters atmospheres are strongly mixed and that large scale spatial and temporal variability are expected in any condensable chemical constituents. Today, observers can already detect longitudinal variations in the emitted thermal flux of the planet. In the next decade, both longitudinal and latitudinal variations in thermal emission and albedo of the hot Jupiters will be observable, expanding the study of weather to extrasolar planets.
Acknowledgement
The authors wish to thank the ISIMA program (http://isima.ucsc.edu) where this project was initiated. This project was supported by NASA Planetary Atmospheres and Origins grants to APS. We thank Christiane Helling for discussions on dust formation, Franck Hersant for discussions on the calculations, Nikole Lewis for sharing her Matlab routines and Franck Selsis for useful input. We also acknowledge the anonymous referee whose numerous comments improved the clarity of this paper.
References
 Ackerman & Marley (2001) Ackerman, A. S. & Marley, M. S. 2001, ApJ, 556, 872
 Adcroft et al. (2004) Adcroft, A., Campin, J.M., Hill, C., & Marshall, J. 2004, Monthly Weather Review, 132
 Agúndez et al. (2012) Agúndez, M., Venot, O., Iro, N., et al. 2012, Astron. & Astrophys, 548, A73
 Andrews et al. (1987) Andrews, D. G., Holton, J. R., & Leovy, C. B. 1987, Middle atmosphere dynamics.
 Becklin & Zuckerman (1988) Becklin, E. E. & Zuckerman, B. 1988, Nature, 336, 656
 Benneke & Seager (2012) Benneke, B. & Seager, S. 2012, ApJ, 753, 100
 Burrows et al. (2011) Burrows, A., Heng, K., & Nampaisarn, T. 2011, ApJ, 736, 47
 Burrows et al. (1997) Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856
 Burrows & Sharp (1999) Burrows, A. & Sharp, C. M. 1999, ApJ, 512, 843
 Chamberlain & Hunten (1987) Chamberlain, J. W. & Hunten, D. M. 1987, Theory of planetary atmospheres. An introduction to their physics andchemistry.
 Chapman & Cowling (1970) Chapman, S. & Cowling, T. 1970, The mathematical theory of nonuniform gases, 3rd edn. (Cambridge mathematical library)
 Charbonneau et al. (2005) Charbonneau, D., Allen, L. E., Megeath, S. T., et al. 2005, ApJ, 626, 523
 Charbonneau et al. (2002) Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377
 Cooper & Showman (2005) Cooper, C. S. & Showman, A. P. 2005, ApJ, 629, L45
 Cooper & Showman (2006) Cooper, C. S. & Showman, A. P. 2006, ApJ, 649, 1048
 Crossfield et al. (2010) Crossfield, I. J. M., Hansen, B. M. S., Harrington, J., et al. 2010, ApJ, 723, 1436
 Cushing et al. (2011) Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al. 2011, ApJ, 743, 50
 De Wit et al. (2012) De Wit, J., Gillon, M., Demory, B.O., & Seager, S. 2012, Astron. & Astrophys, 548, A128
 Demory et al. (2011) Demory, B.O., Seager, S., Madhusudhan, N., et al. 2011, ApJ, 735, L12
 Désert et al. (2008) Désert, J.M., VidalMadjar, A., Lecavelier Des Etangs, A., et al. 2008, Astron. & Astrophys, 492, 585
 DobbsDixon & Agol (2012) DobbsDixon, I. & Agol, E. 2012, ArXiv eprints
 DobbsDixon et al. (2012) DobbsDixon, I., Agol, E., & Burrows, A. 2012, ApJ, 751, 87
 DobbsDixon et al. (2010) DobbsDixon, I., Cumming, A., & Lin, D. N. C. 2010, ApJ, 710, 1395
 DobbsDixon & Lin (2008) DobbsDixon, I. & Lin, D. N. C. 2008, ApJ, 673, 513
 Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419
 Fortney et al. (2005) Fortney, J. J., Marley, M. S., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Astronomische Nachrichten, 326, 925
 Fortney et al. (2010) Fortney, J. J., Shabram, M., Showman, A. P., et al. 2010, ApJ, 709, 1396
 Freedman et al. (2008) Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504
 Freytag et al. (2010) Freytag, B., Allard, F., Ludwig, H.G., Homeier, D., & Steffen, M. 2010, Astron. & Astrophys, 513, A19
 Goody (1961) Goody, Y. 1961, Atmospheric Radiation Theoretical Basis (Oxford university press)
 Heng (2012) Heng, K. 2012, ApJ, 761, L1
 Heng et al. (2011a) Heng, K., Frierson, D. M. W., & Phillipps, P. J. 2011a, MNRAS, 418, 2669
 Heng et al. (2012) Heng, K., Hayek, W., Pont, F., & Sing, D. K. 2012, MNRAS, 420, 20
 Heng et al. (2011b) Heng, K., Menou, K., & Phillipps, P. J. 2011b, MNRAS, 413, 2380
 Holton (1992) Holton, J. R. 1992, An introduction to dynamic meteorology
 Hubeny et al. (2003) Hubeny, I., Burrows, A., & Sudarsky, D. 2003, ApJ, 594, 1011
 Iro et al. (2005) Iro, N., Bézard, B., & Guillot, T. 2005, Astron. & Astrophys, 436, 719
 Kalnay (2003) Kalnay, E. 2003, Atmospheric modeling, data assimilation, and predictability (Cambridge, U.K., New York: Cambridge University Press)
 Kataria et al. (2013) Kataria, T., Showman, A. P., Lewis, N. K., et al. 2013, ApJ, 767, 76
 Kirkpatrick (2005) Kirkpatrick, J. D. 2005, ARA&A, 43, 195
 Knutson et al. (2008) Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., & Megeath, S. T. 2008, ApJ, 673, 526
 Knutson et al. (2007) Knutson, H. A., Charbonneau, D., Allen, L. E., et al. 2007, Nature, 447, 183
 Knutson et al. (2009) Knutson, H. A., Charbonneau, D., Cowan, N. B., et al. 2009, ApJ, 690, 822
 Knutson et al. (2012) Knutson, H. A., Lewis, N., Fortney, J. J., et al. 2012, ApJ, 754, 22
 Lecavelier Des Etangs et al. (2008) Lecavelier Des Etangs, A., Pont, F., VidalMadjar, A., & Sing, D. 2008, Astron. & Astrophys, 481, L83
 Lee et al. (2012) Lee, J.M., Fletcher, L. N., & Irwin, P. G. J. 2012, MNRAS, 420, 170
 Lewis et al. (2010) Lewis, N. K., Showman, A. P., Fortney, J. J., et al. 2010, ApJ, 720, 344
 Li & Wang (2003) Li, Z. & Wang, H. 2003, Physical Review E  Statistical, Nonlinear and Soft Matter Physics, 68, 061206
 Line et al. (2010) Line, M. R., Liang, M. C., & Yung, Y. L. 2010, ApJ, 717, 496
 Line et al. (2011) Line, M. R., Vasisht, G., Chen, P., Angerhausen, D., & Yung, Y. L. 2011, ApJ, 738, 32
 Lodders (2002) Lodders, K. 2002, ApJ, 577, 974
 Lodders & Fegley (2002) Lodders, K. & Fegley, B. 2002, Icarus, 155, 393
 Madhusudhan et al. (2011) Madhusudhan, N., Mousis, O., Johnson, T. V., & Lunine, J. I. 2011, ApJ, 743, 191
 Madhusudhan & Seager (2009) Madhusudhan, N. & Seager, S. 2009, ApJ, 707, 24
 Madhusudhan & Seager (2011) Madhusudhan, N. & Seager, S. 2011, ApJ, 729, 41
 Majeau et al. (2012) Majeau, C., Agol, E., & Cowan, N. B. 2012, ApJ, 747, L20
 Marley & McKay (1999) Marley, M. S. & McKay, C. P. 1999, Icarus, 138, 268
 Marley et al. (2010) Marley, M. S., Saumon, D., & Goldblatt, C. 2010, ApJ, 723, L117
 Marley et al. (1996) Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Science, 272, 1919
 Marley et al. (2002) Marley, M. S., Seager, S., Saumon, D., et al. 2002, ApJ, 568, 335
 Mayor & Queloz (1995) Mayor, M. & Queloz, D. 1995, Nature, 378, 355
 McKay et al. (1989) McKay, C. P., Pollack, J. B., & Courtin, R. 1989, Icarus, 80, 23
 Moses et al. (2011) Moses, J. I., Visscher, C., Fortney, J. J., et al. 2011, ApJ, 737, 15
 Parmentier & Guillot (2013) Parmentier, V. & Guillot, T. 2013, submitted to A&A
 Parmentier et al. (2013) Parmentier, V., Guillot, T., J., F., & M., M. 2013, in prep.
 Perna et al. (2012) Perna, R., Heng, K., & Pont, F. 2012, ApJ, 751, 59
 Perna et al. (2010) Perna, R., Menou, K., & Rauscher, E. 2010, ApJ, 719, 1421
 Pont et al. (2013) Pont, F., Sing, D. K., Gibson, N. P., et al. 2013, MNRAS
 Probstein (1968) Probstein, R. F. 1968, Problems of hydrodynamics and continuum mechanics (Society for industrial and applied mathematics)
 Pruppacher & Klett (1978) Pruppacher, H. & Klett, J. 1978, Microphysics of clouds and precipitation (D. Reidel publishing compagny)
 Rauscher & Menou (2010) Rauscher, E. & Menou, K. 2010, ApJ, 714, 1334
 Rauscher & Menou (2012a) Rauscher, E. & Menou, K. 2012a, ApJ, 750, 96
 Rauscher & Menou (2012b) Rauscher, E. & Menou, K. 2012b, ApJ, 745, 78
 Rauscher & Menou (2013) Rauscher, E. & Menou, K. 2013, ApJ, 764, 103
 Rosner (2000) Rosner, D. E. 2000, Transport processes in chemically reacting flow systems (Dover : Mineola)
 Shapiro (1970) Shapiro, R. 1970, Rev. Geophys., 8, 359
 Showman et al. (2008) Showman, A. P., Cooper, C. S., Fortney, J. J., & Marley, M. S. 2008, ApJ, 682, 559
 Showman et al. (2013) Showman, A. P., Fortney, J. J., Lewis, N. K., & Shabram, M. 2013, ApJ, 762, 24
 Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564
 Showman & Guillot (2002) Showman, A. P. & Guillot, T. 2002, Astron. & Astrophys, 385, 166
 Showman & Polvani (2011) Showman, A. P. & Polvani, L. M. 2011, ApJ, 738, 71
 Smith (1998) Smith, M. D. 1998, Icarus, 132, 176
 Spiegel & Burrows (2010) Spiegel, D. S. & Burrows, A. 2010, ApJ, 722, 871
 Spiegel et al. (2009) Spiegel, D. S., Silverio, K., & Burrows, A. 2009, ApJ, 699, 1487
 Stiel & Thodos (1963) Stiel, L. & Thodos, G. 1963, Industrial and Engineering Chemistry Fundamentals, 2, 233
 Tinetti et al. (2012) Tinetti, G., Beaulieu, J. P., Henning, T., et al. 2012, Experimental Astronomy, 34, 311
 Venot et al. (2012) Venot, O., Hébrard, E., Agúndez, M., et al. 2012, Astron. & Astrophys, 546, A43
 Visscher et al. (2006) Visscher, C., Lodders, K., & Fegley, Jr., B. 2006, ApJ, 648, 1181
 Woitke & Helling (2003) Woitke, P. & Helling, C. 2003, Astron. & Astrophys, 399, 297
 Woitke & Helling (2004) Woitke, P. & Helling, C. 2004, Astron. & Astrophys, 414, 335
 Youdin & Mitchell (2010) Youdin, A. N. & Mitchell, J. L. 2010, ApJ, 721, 1113
 Zahnle et al. (2009a) Zahnle, K., Marley, M. S., & Fortney, J. J. 2009a, ArXiv eprints
 Zahnle et al. (2009b) Zahnle, K., Marley, M. S., Freedman, R. S., Lodders, K., & Fortney, J. J. 2009b, ApJ, 701, L20
Appendix A Departure from the Cunningham velocity
The StokesCunningham velocity defined in Eq. (3) is derived under the assumption of low Reynolds number. Therefore it is not valid for turbulent flow and other expressions may be used when the Reynolds number increases. Here we derive better laws for intermediate and large Reynolds number.
a.1 Low Knudsen number
For small Reynolds number and small Knudsen number, the drag force exerted by a fluid on a sphere at rest is considered proportional to the kinetic energy of the fluid and the projected area of the sphere. The coefficient of proportionality, or drag coefficient, is given by :
(24) 
Then, equating gravity and drag forces leads to the settling velocity of a particle in an atmosphere :
(25) 
Where is the density of the particle. For small Reynolds numbers and high Knudsen number, is constant and the settling velocity is the Stokes velocity. When increasing the Reynolds number, the non linear terms of the NavierStokes equation become important and is no longer constant. We used tabulated values of the drag coefficient as a function of the Reynolds number given by Pruppacher & Klett (1978). We assume that when the Reynolds number reaches 1 to stay consistent with Stokes flow and that reaches its asymptotic value, , when and fit the relationship :
(26) 
Then we follow the same method as Ackerman & Marley (2001). Noting that :
(27) 
is independent of the velocity, we use the fit of Eq. 26 and extract the velocity :
(28) 
a.2 High Knudsen number
In the freemolecular regime, calculations have been made by Probstein (1968) leading to an expression for the drag coefficient:
(29) 
where is the ratio of the object velocity over the thermal speed of the gas ( where is the sound speed) and is the error function.
For velocities much smaller than the sound speed, and we can use an equivalent of the error function in 0:
(30) 
Using this equation inside Eq. (29) and taking the limit , the term goes to 1 and the terms proportional to cancels out leading to:
(31) 
For velocities much greater than the sound speed, the limit of Eq. (29) when is:
(32) 
In order to simplify Eq. (29) we use the following expression for the drag coefficient at high Knudsen number :
(33) 
Our approximation fits correctly the exact expression in the limit of low and high velocities. In between the difference to the exact expression is at most . Replacing by its value in Eq. (25) we obtain a second order equation for the velocity :
(34) 
which leads to :
(35) 
with . When the speed becomes small compared to the sound speed () we obtain :
(36) 
which is in good agreement with Eq. (3), derived for high Knudsen number and small Reynolds numbers.
a.3 Comparison with the Cunningham velocity
Figure 17 shows the ratio of the Cunningham velocity (see Eq. (3)) to the ones we just derived. The difference is noticeable only for particles of the order of at pressures less than the level and exceeding the level. This difference is always less than one order of magnitude and concern only a tiny portion of the parameter space which has little relevance to our study (the largest particle sizes considered in our 3D models is m). Thus we decided to neglect this discrepancy in the main study.