EDP Sciences
Free Access
Issue
A&A
Volume 594, October 2016
Article Number A48
Number of page(s) 24
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201628606
Published online 10 October 2016

© ESO, 2016

1. Introduction

The detection of atoms and molecules such as Na (Charbonneau et al. 2002), K (Sing et al. 2011a), H2O (Swain et al. 2009), CO (Snellen et al. 2010) and CH4 (Swain et al. 2008) in the atmospheres of hot Jupiter exoplanets have provided a first glimpse into the chemical composition of these planets. One of the major results of observational efforts is the detection of clouds/hazes in exoplanet atmospheres. Two key signatures of atmospheric cloud in transmission spectra are an almost featureless optical spectrum, interpreted as Rayleigh scattering by cloud particles, and a relatively flat infrared spectrum with weakened molecular signatures. Hubble and Spitzer observations of exoplanets such as HD 189733b (Pont et al. 2013), WASP-31b (Sing et al. 2015), WASP-6b (Nikolov et al. 2015) and GJ 1214b (Kreidberg et al. 2014) show evidence for clouds in their atmosphere. Deming et al. (2013) present detailed HST observations of HD 209458b, and show that the muted spectral H2O and Na features can be explained by the additional grey opacity of clouds. Sing et al. (2016)’s study of ten hot Jupiters, incl. HD 189733b, in their Hubble and Spitzer observing program, provides compelling evidence that the muted water lines in the infrared can be attributed to high altitude mineral clouds. They show that the signatures of clouds are present across a diverse range of planet equilibrium temperatures (Teq = 960 ... 2510 K), radii (Rpl = 0.96 ... 1.89 RJ) and mass (Mpl = 0.21 ... 1.40 MJ). Lecavelier Des Etangs et al. (2008) and Wakeford & Sing (2015) have shown that a sub-micron cloud composition of silicate minerals such as MgSiO3[s] can fit the observed Rayleigh slope of HD 189733b. Star spots (McCullough et al. 2014) have been discussed as possible explanation for the HD 189733b spectral features. However, this does not rule out the presence of atmospheric cloud. Evans et al. (2013)’s HST analysis of secondary transit observations for HD 189733b show a blueward slope in the geometric albedo at optical wavelengths, suggesting the presence of wavelength dependent backscattering cloud particles.

In recent years, global circulation models (GCM) and radiative-hydrodynamic (RHD) models have been developed in order to understand the global dynamics of extrasolar atmospheres. To date, most groups have focused on modelling the hot Jupiters, HD 189733b (e.g. Showman et al. 2009; Dobbs-Dixon & Agol 2013), HD 209458b (e.g. Showman et al. 2009; Dobbs-Dixon et al. 2010; Rauscher & Menou 2010; Heng et al. 2011; Mayne et al. 2014) and the warm sub-Neptune GJ 1214b (e.g. Kataria et al. 2014; Charnay et al. 2015a). Codes differ in their implementation and hydrodynamic assumptions (see summary in Mayne et al. 2014) and radiative-transfer schemes (Amundsen et al. 2014) as well as the complexity with which gas and cloud chemistry is implemented. However, the overall global thermal structures and jet patterns remain similar for HD 209458b simulations (Heng & Showman 2015). Hot Jupiter simulations of HD 189733b (e.g. Showman et al. 2009; Dobbs-Dixon & Agol 2013) have consistently reproduced offsets of maximum temperature to the East of the sub-stellar point, in line with Knutson et al. (2009)Spitzer thermal maps of atmospheres. Kataria et al. (2016) modelled the hydrodynamics, thermal structure of nine hot Jupiters from the Sing et al. (2016) observing program. They post-process their simulations and examine the differences in the gas phase chemical equilibrium abundances between each planet. Other simulations have investigated the effect of eccentricity (Lewis et al. 2010, 2014), orbital distance (Showman et al. 2015), rotation (Rauscher & Kempton 2014; Showman et al. 2015) and irradiation (Perna et al. 2012; Komacek & Showman 2016) on the hydrodynamics of the atmosphere.

Thus far, investigations of atmospheric cloud properties of hot Jupiters using GCM/RHD simulations have not included a model for describing the formation of clouds. Parmentier et al. (2013) simulated the mixing of constant μm sized (a = 0.1 ... 10μm) cloud tracer particles for HD 209458b. They showed that sub-micron sized grains are likely to not settle out of the upper atmosphere. Charnay et al. (2015b) simulated the warm sub-Neptune GJ 1214b. They applied a constant particle size approach for KCl[s] and ZnS[s] using the phase equilibrium condensate chemistry from Morley et al. (2012). Oreshenko et al. (2016) used a phase equilibrium scheme for MgSiO3[s], Mg2SiO4[s], Fe[s], TiO[s] and Al2O3[s]. They post-process their GCM to model the optical and infrared wavelength phase curves of Kelper-7b (Demory et al. 2013). Parmentier et al. (2016) simulated a suite of hot Jupiter atmospheric models at equilibrium temperatures ranging from Teq = 1300 ... 2200 K. They applied a phase equilibrium approach to model single, homogeneous species Fe[s], Al2O3[s], MgSiO3[s], Cr[s], MnS[s] and Na2S[s] clouds for prescribed, constant grain sizes. They then post-process their suite of models to investigate optical phase curve offsets in the Kepler bandpass. A common feature of the above cloud modelling approaches is the prescription of a constant grain size and homogeneous mineral composition of cloud particles.

We present the first 3D hot Jupiter RHD model with a coupled kinetic cloud formation, opacity and transport of cloud particles. We combine the HD 189733b atmosphere simulation of Dobbs-Dixon & Agol (2013) with a time-dependent, kinetic non-equilibrium cloud formation model based on Woitke & Helling (2003, 2004), Helling & Woitke (2006), Helling et al. (2008). Section 2 provides a summary of the RHD model and the cloud formation model. Section 3 outlines our cloud formation approach, numerical approach and a discussion on the present state of convergence. In Sect. 4 we present the 3D temperature and horizontal velocity field and local cloud properties such as cloud number density, local cloud particle sizes, material composition and element depletion. Section 6 presents wavelength dependent cloud, gas and total opacities and examines the stellar heating rates resulting from the cloud opacity. Section 7 contains the discussion and Sect. 8 outlines our summary and conclusions.

2. Model description

The 3D modelling of dynamic clouds for hot Jupiter atmospheres requires a coupled, time-dependent hydrodynamic, radiative-transfer and cloud formation model scheme. We time dependently evolve the 3D Navier-Stokes equations and two-stream radiative transfer scheme as in Dobbs-Dixon & Agol (2013) coupled with a time dependent, 3D, kinetic phase-non-equilibrium cloud formation and gas/dust phase chemistry module based on Woitke & Helling (2003, 2004), Helling & Woitke (2006), Helling et al. (2008). Gravitational settling of the cloud particles and element conservation under the influence of cloud formation is taken into account.

2.1. 3D radiative-hydrodynamic model

The radiative-hydrodynamic (RHD) model applied in this study combines the fully compressible 3D Navier-Stokes equations to a two-stream, frequency dependent radiative transfer scheme. Equations (1)(3)represent the continuity, momentum and energy conservation of a local fluid element respectively, (1)where ρgas [g cm-3] is the hydrodynamic gas density and ugas [cm s-1] the hydrodynamic gas velocity. (2)where Pgas [dyn cm-2] is the local gas pressure, g = ger [cm s-2] the gravitational acceleration in the radial direction, Ω [rad s-1] the rotational frequency of the planet, r(r, φ, θ) [cm] the spherical coordinate radial distance vector and ν = η/ρgas [cm2 s-1] the constant kinematic viscosity. (3)where egas [erg cm-3] is the internal energy density of the fluid element, Fr [erg cm-2 s-1] the radiative flux in the radial direction, S [erg cm-3 s-1] the incident stellar energy density per second and Dν [erg cm-3 s-1] the local energy density per second from gaseous viscous heating. For the radiation field, the heating due to stellar irradiation S [erg cm-3 s-1] is given by (4)where R [cm] and a [cm] are the stellar radius and semi-major axis, respectively, r [cm] the radial coordinate, μ = cosθ the cosine of the angle between the normal and the incident stellar photons, accounting for the slant path from a vertically defined optical depth and τb, the optical depth of stellar photons for wavelength band b given by (5)Where 1.66 is the diffusivity factor, an approximation that accounts for an exponential integral that arises when taking the first moment of the intensity to calculate the flux (Elsasser 1942) and κb [cm2 g-1] the opacity for wavelength band b (Eq. (27)). In the two-stream approximation, the net radial flux, [erg cm-2 s-1], for wavelength band b in the upward or downward direction is (6)The propagation of radiation intensity in the downward F,b and upward F,b direction at each cell (two-steam approximation) is given by (7)and (8)respectively, where Sb [erg cm-2 s-1] is the source function given by the wavelength integrated Planck function at the local temperature (T = Tgas), (9)Sb(Tbot) in Eq. (8)is the contribution of energy emitted from the interior of the planet, below which the planet is assumed to emit as a blackbody. This interior contribution to the upward flux is fixed to match the observed radius (Dobbs-Dixon & Agol 2013). Stellar heating is fully accounted for in Eq. (4).

More details on the RHD model implementation of the 3D Navier-Stokes equations and the two-stream approximation for radiative transfer can be found in Dobbs-Dixon et al. (2010) and Dobbs-Dixon & Agol (2013), respectively.

2.2. Cloud formation model

On Earth, water clouds are formed when particles made of sand, ash or ocean spray salt are lofted up into atmospheric layers where water vapour is supersaturated. These cloud condensation nuclei or aerosols are the start of the formation of water droplets due to the surface area they provide. Only a small supersaturation ratio (S ~ 1.001; Korolev & Mazin 2003) is required to start water condensation on these surfaces in the Earth’s atmosphere. Since no solid crust exists on hot Jupiters, cloud formation has to start from the gas phase by the nucleation of seed particles. Possible seed formation species have been reviewed by Helling & Fomins (2013) which found TiO2 to be an efficient nucleating gas species for Brown Dwarf and hot Jupiter atmospheres. TiO2[s] is a high-temperature condensate which forms seed particles from a chemical pathway of successive (TiO2)N stable clusters through homogeneous nucleation. Based on the results of Jeong et al. (2003) and Lee et al. (2015b), we are in the position to apply a modified classical nucleation theory where we use the thermodynamical data for the individual cluster sizes (TiO2)N to derive the seed formation rate dependent on the local gas temperature and gas density, and hence calculate the number of locally formed seed particles. After the first seed particles are nucleated from the gas phase, chemical surface reactions grow the grain bulk. This surface can grow or evaporate while the particles are transported across the globe, depending on the local thermodynamic and chemical conditions. The composition of the grain bulk changes when various materials become thermally stable or unstable. This process results in a solid mantle of mixed mineral composition where the seed particle contributes a negligible amount of volume compared to the total volume of the whole cloud particle. The local growth and evaporation chemical surface processes occur at second-to-minute timescales (Helling et al. 2001, 2004) meaning the properties (size, composition) of the cloud particles adapt quickly to the local gaseous thermodynamic properties.

To describe the cloud formation process we use the set of dust moment equations derived by Woitke & Helling (2003, 2004), Helling & Woitke (2006), Helling et al. (2008). The dust moments Lj(r) [cmj g-1] (j = 0, 1, 2, 3) are the local integrated particle size distribution, weighted by a power of the grain volume Vj/ 3, defined as (10)where f(V,r) [cm-6] is the distribution of particles in volume space and Vl [cm3] the volume of a seed particle. The conservation equation of dust moments is given by (Woitke & Helling 2003) (11)which has the same conservation transport equation structure as the Navier-Stokes equations. The local time evolution of cloud properties are represented by the source/sink terms (r.h.s.) describing the seed formation and the grain growth/evaporation chemical processes. The second l.h.s. term denotes the advective flux of the conserved quantity (ρgasLj) through space, given by the dust phase hydrodynamic fluid velocity ud(r) [cm s-1]. The relative velocity of gas and dust is given by the drift velocity vdr(r) [cm s-1] defined as (Woitke & Helling 2003) (12)where ugas(r) [cm s-1] is the gas hydrodynamic velocity. Following the analysis in Woitke & Helling (2003) the mean equilibrium drift velocity ⟩ (r) [cm s-1] in the large Knudsen number regime (Kn 1) is given by (13)where ρd [g cm-3] is the bulk (material) dust density, cT [cm s-1] the speed of sound, ger [cm s-2] the gravitational acceleration in the radial direction and a ⟩ (r) [cm] the local mean grain size.

The r.h.s. terms in Eq. (11) denote the local chemical processes that alter the local dust moments; nucleation: J(r) ≥ 0 cm-3 s-1, growth: χnet(r) > 0 cm s-1 and evaporation: χnet(r) < 0 cm s-1. The nucleation rate J(r) [cm-3 s-1] for homomolecular homogeneous nucleation of seed particles, applying modified classical nucleation theory (Jeong et al. 2003; Woitke & Helling 2003; Helling & Fomins 2013; Lee et al. 2015b) is given by (14)where f(1,t) [cm-3] is the number density of the seed forming gas species, τgr [s] the growth timescale of the critical cluster size N, Z(N) the Zeldovich factor (the contribution from Brownian motion to the nucleation rate), ΔG(N) [kJ mol-1] the change in Gibbs free energy from the formation of cluster size N and S(T) the supersaturation ratio defined as (15)where ns is the gas phase number density of species s and (T) the vapour pressure of species s, which is a function of local gas temperature.

The net growth/evaporation velocity χnet(r) [cm s-1] of a grain due to chemical surface reactions (Gail & Sedlmayr 1986; Helling & Woitke 2006) is (16)where r is the index for the chemical surface reaction, the volume increment of the solid s by reaction r, the particle density of the key reactant in the gas phase, the relative thermal velocity () of the gas species taking part in reaction r, αr the sticking coefficient of reaction r and the stoichiometric factor of the key reactant in reaction r. Sr is the reaction supersaturation ratio (Helling & Woitke 2006) and 1//Vtot the volume of solid s, Vs, to the total grain volume containing all species Vtot = ∑ sVs. No prior assumptions about the particular grain size distribution or grain sizes are required to compute the number density of cloud properties. Should detailed analysis using a distribution be required, it can be reconstructed by assuming a distribution shape from the dust moments a posteriori (Helling et al. 2008; Stark et al. 2015).

The local dust number density nd(r) [cm-3] and mean grain radius a ⟩ (r) [cm] are calculated from the dust moments by (17)and (18)respectively.

The composition of material forming on the grain mantle changes as a result of local chemical and thermodynamic conditions and the thermal stability of each material in those conditions. The volume of a specific solid mineral s depends on the growth/evaporation rate of that material. The volume of each material s can be described by a separate moment conservation equation for the third dust moment, L3,s(r) [cm3 g-1], (Helling et al. 2008) (19)where the growth velocity of the solid s, , is given by (20)The local volume fraction Vs/Vtot of each species s is calculated from the L3,s(r) dust moment using the identity (Woitke & Helling 2004; Helling et al. 2008) (21)

2.3. Element abundance and gas phase chemistry

The local gas-phase composition is an input for the cloud formation process as it determines (along with temperature) the cloud formation nucleation and growth/evaporation rates. The elements involved in the cloud formation process are altered by depletion or replenishment depending on the dominating cloud formation processes (nucleation, growth, evaporation). The depletion/enrichment of gas phase elements εi(r) (abundance ratio of element i to hydrogen; εi(r) = ni/nH) by cloud particle growth/evaporation (Helling & Woitke 2006) is given by (22)where i is the index of the element that contributes to the cloud formation process and nH the local total number density of hydrogen. The first r.h.s. terms describes the consumption of elements from the nucleation process. The second r.h.s. term denotes the source (evaporation) and sink (growth) of elements as a result of cloud particle chemical surface reactions. The second l.h.s. term describes the advection of εi(r) through space at the local gas velocity ugas(r).

The composition of the local gas phase is calculated assuming chemical equilibrium based on the element abundances derived from Eq. (22). Using these local element abundances, we use a computationally efficient, hierarchical chemical abundance calculation introduced in Helling et al. (2001).

2.4. Cloud and gas opacity

Cloud particles are a large source of opacity, absorbing and scattering photons at optical and infrared wavelengths. We apply spherical particle Mie theory (Mie 1908) in combination with effective medium theory to calculate the mixed material cloud opacity. Extinction efficiencies, Qext(λ,a), for the local mean grain sizes a are calculated based on the Bohren & Huffman (1983) BHMie routines. We use two approximations for the size parameter (x = 2πa/λ) limits of Mie theory. For large size parameters x ≥ 1000 we use the large particle, hard sphere scattering approximation where the absorption efficiency asymptotically tends towards zero, Qabs = 0, and all extinction is assumed to be from scattering. In the large particle limit, extinction efficiency Qext is then (23)For small size parameters x< 10-6 we use the small metallic sphere particle limit approximation, as, for example, outlined by Gail & Sedlmayr (2014), where (24)and (25)where ϵ is the complex dielectric function and α = (ϵ−1)(ϵ + 2). The second term in Qabs contains the effect due to induction of eddy currents on the grain surface by the electromagnetic field of the photons. The Qsca equation calculates the contribution to the total extinction from Rayleigh scattering. This approximation has been shown to produce similar results to Mie theory for very small size parameters (Gail & Sedlmayr 2014). For all other size parameters the full Mie calculation is carried out. The mass extinction coefficient κλ,cloud [cm2 g-1] at wavelength λ is then given by (26)where a = ⟨ a is the mean grain size from Eq. (18). We calculate the cloud opacity at wavelengths corresponding to the 31 wavelength opacity bin edges used in Showman et al. (2009). The cloud opacity is then averaged across the wavelength range of the bins to calculate the cloud opacity for that band. For gas opacities, we use the temperature, density tabulated frequency dependent results from Sharp & Burrows (2007) Planck averaged over the wavelength bins of Showman et al. (2009). The total band opacity from the gas and cloud κb,total [cm2 g-1] is then given by (27)This local total opacity is treated as a purely absorptive extinction in the radiative transfer scheme.

2.4.1. Effective medium theory

Each mineral material present in the cloud particle will contribute their specific optical (n,k) constants, weighted by their local volume fraction Vs/Vtot (Eq. (21)) of the material in each grain. Effective optical constants for the material mixtures are calculated using effective medium theory. We primarily use the numerical Bruggeman method (Bruggeman 1935) given by (28)where Vs/Vtot is the volume fraction of solid species s, ϵs the dielectric function of solid species s and ϵav the effective, average dielectric function over the total cloud particle volume. A Newton-Raphson minimisation scheme is then applied to solve for ϵav. In cases of rare non-convergence at the furthest infrared bands (λ = 20.00−46.00, 46324.68 μm) where the (n,k) experimental values for materials are most uncertain; we apply the analytic Landau-Lifshitz-Looyenga [LLL] method of Looyenga (1965) given by (29)Values for the optical constants are the same as those used in Lee et al. (2015a). Most solid materials published experimental optical constants do not cover the wavelength range (λ = 0.261...324.68 μm) required for the radiative-transfer scheme. To overcome this, we extrapolate the available data for each material to shorter and longer wavelengths. For wavelengths shorter than published data we assume that (n,k) remains constant. At longer wavelengths, for the non-conducting material considered in this study, we assume that n remains constant while k is reduced linearly from the last data point to the longer wavelengths. This can make the material effective optical constants for the infrared wavelength bins, λ = 20.00−46.00 μm and λ = 46.00−324.68 μm uncertain depending on the volume fraction and available data of each species.

3. Approach

In this section we outline our cloud formation approach, numerical approach, initial conditions and convergence properties for our RHD and cloud formation model. The addition of our cloud formation model the RHD model adds additional costs to the simulation times. For example, the ~60 simulated days presented here took ~20 Earth days using 64 cores1.

3.1. Cloud formation and element abundance

We consider the homogeneous nucleation of TiO2 seed particles based on (TiO2)N cluster data from Jeong et al. (2003), Lee et al. (2015b). We consider five simultaneous surface materials, TiO2[s], SiO[s], SiO2[s], Mg2SiO4[s], MgSiO3[s] with 22 of the corresponding surface chemical reactions found in Helling et al. (2008). The local cloud properties are locally time-dependently computed for each computational domain while the flux of the moments through 3D space can be calculated using a normal advection scheme. This is significant improvement over our previous 1D methods which rely on mixing timescale arguments to calculate clouds properties and do not consider transport in horizontal directions. We solve Eq. (22)for four element abundances: Ti, O, Si and Mg and assume a constant solar element abundance for all other elements. We assume horizontal and meridional frictional coupling of the dust and gas phase (ud,h,m = ugas,h,m). Vertical decoupling between the dust and gas phase is applied given by Eq. (12)(ud,vert = ugas,vert + vdr).

3.2. Numerical approach

We evolve the Navier-Stokes equations (Eqs. (1)(3)) over a spherical grid with a resolution of (Nr, Nφ, Nθ) = (100, 160, 64), where r is the radial direction, φ the longitude and θ the latitude. The upper radial boundary is allowed to vary between ~10-5 and 10-4 bar depending on dynamical properties and the lower boundary is fixed at ~500 bar. Vertical velocity dampening (sponge zone) is implemented near the upper boundaries of the simulation, common to GCM models (e.g. Mayne et al. 2014). We account for flow over polar regions by the method found in Dobbs-Dixon et al. (2012).

One integration of the kinetic cloud formation chemistry can take ~2030 times the computational time of a single hydrodynamic timestep. Therefore, the cloud formation chemistry (r.h.s. Eq. (11)) is integrated every 10 hydrodynamic timesteps to update the local cloud particle properties. The drift velocity vdr(r) and cloud opacity κb,cloud are updated after the cloud formation chemistry. At every hydrodynamic timestep, cloud moments; Lj(r), material volume composition; L3,s(r) and gaseous element abundances; εi(r) are advected around the globe by calculation of the advective term in Eqs. (11), (19)and (22).

Additionally, we take a number of physically based assumptions to reduce the number of cloud chemistry iterations required and to ensure physical solutions to the cloud properties. During an evaporation event, the maximum integration step size is reduced by half to capture the evaporation process more consistently. We limit evaporation of grains to the seed particle size (~0.001 μm); should the integrator attempt to evaporate the moment solutions below seed particle sizes, the end solution is assumed to be at the local seed particle values. At this point, should the seed particles be thermally unstable ( cm s-1) then all seed particles are assumed to be evaporated (i.e. Lj(r) = 0) and the Ti and O elements returned to the gas phase by Eq. (30). This condition is only met at the hottest, deepest parts of the atmosphere (Tgas> 2300 K) in our simulations. Thermal stability may prevail to higher gas temperatures if other high-temperature condensates are included.

We limit the calculation of the dust chemistry to regions where the number density of grains is nd(r) > 10-10 cm-3 unless the local nucleation rate is J(r) > 10-10 cm-3 s-1. This condition ensures that only regions that are efficiently nucleating or already contain a significant number of cloud particles are integrated. This criteria ensures that deeper atmospheric regions (pgas> 10 bar) which have a negligible nucleation but large growth rate, where a small number of cloud particles grow rapidly large (>1 cm), do not occur. Without this criteria, these regions become computationally challenging as the drift velocity (Eq. (13)) becomes large (~speed of sound) for these particles. The hydrodynamic Courant-Friedrich-Levi [CFL] timestep condition then limits the hydrodynamic timestep to unfeasibly low values. The cloud opacity and drift velocity in regions which contain very little cloud particles (nd(r) < 10-10 cm-3) are assumed to be zero. For regions where nd(r) > 10-10 cm-3, a lower bound cloud opacity of κcloud = 10-7 cm2 g-1 across all wavelength bands is implemented to aid numerical stability. This is required since the results of Mie theory in certain seed particle regions can approach floating point limits. This corresponds to a lower bound of TiO2[s] seed particle opacity (κcloud(λ, aseed) 10-7 cm2 g-1) at optical wavelengths.

Furthermore, cells with very small growth/evaporation rates of | χnet(r) | < 10-20 cm s-1 are assumed to remain constant with respect to the dust moments. Only cells which have local conditions that are significant departures from the equilibrium solution (χnet(r) = 0) are integrated in time. We found integrating cells with | χnet(r) | < 10-20 cm s-1 was computationally prohibitive and did not produce significantly different results.

The need to only update the r.h.s. of the dust moment equation every 10th hydrodynamic timestep may lead to a fast transport of cloud particles into high temperature regions where they will evaporate. This volatile material would evaporate rapidly at some material dependent “evaporation window” as it passed into these unstable regions. The dust moment and element conservation equations become numerically stiff for such an intense evaporation, hence, the dust moments would approach phase equilibrium (i.e. when evaporation stops) in very small timesteps. Small timesteps for the intense evaporation regions are necessary in order to solve the element conservation equation (Eq. (22)) to the best possible precision. To overcome this numerical challenge, we introduce a scheme where, should the integration timestep be too low at the beginning of the cloud chemistry integration, a fraction of the volatile materials are instantly evaporated back into the gas phase. This process is repeated until integration of cloud properties can be computed in a reasonable time. Growing or (near-)stable (| χs(r) | > 10-30 cm s-1) material are not altered and assumed to remain constant. The return of elements to the gas phase from each evaporation species s due to an instantaneous evaporation event is given by (Woitke 2006) (30)where is the element abundance before the instant evaporation step, νi,s the stoichiometric coefficient of element i in species s and 1.427 amu the conversion factor between gaseous mass density ρgas and Hydrogen nuclei density n⟨ H ⟩ and ΔL3,s(r) the difference in grain volume of species s before and after instant evaporation.

This scheme has the added benefit of evenly spreading computational load, since each evaporating cell has a more equal work load. Additionally, since the surface chemical growth of the particles occurs in second-minute timescales (Helling et al. 2001, 2004); if too much of a material is instantly evaporated off, the material can quickly grow back to an equilibrium solution before the end of one hydrodynamical time-step.

3.3. Initial conditions

For initial conditions, we use a well converged RHD model (total simulated time: ~420 Earth days) of Dobbs-Dixon & Agol (2013) which included a parameterised cloud opacity (Eq. (31)). This parameterised opacity is switched off in our simulations. The initial 3D (Tgas, pgas, vgas) structures do not vary significantly from the published results. As in Dobbs-Dixon & Agol (2013) we use a tidally locked HD 189733b set-up with the rotation rate equal to the orbital period (2.22 days). We assume an initial solar abundance of elements throughout the globe given by Asplund et al. (2009).

To set up the 3D RHD simulation, the dust properties are integrated each hydrodynamic timestep for the first 100 iterations and the effects of cloud opacity are neglected until ~5.5 Earth days into the simulation. During this time, larger sized grains with a ⟩ > 1 μm will have gravitationally settled from the upper atmosphere to their pressure supported regions (~1 bar). After these first steps, the opacity of the cloud particles at all positions is accounted for in the radiative transfer scheme.

3.4. Convergence tests

We investigate the present state of the model convergence by examining the the horizontal kinetic energy density Ekin,h = ρgas/2 [erg cm-3] and root-mean-square (rms) horizontal + meridional velocity [m s-1] zonal and meridional averaged at pressure iso-bars. These two quantities show how the global hydrodynamic velocity structure of the atmosphere is changing with atmospheric pressure and with time to check the state of the simulation with respect to a possible statistical steady state. We introduce two properties for the cloud structure, the cloud number density nd [cm-3] (Eq. (17)) and the equilibrium vertical drift velocity vdr [m s-1] (Eq. (13)), which are zonally and meridionally averaged at each iso-bar. Together, these two quantities take into account the time evolution of the global density structure of cloud material as well as the grain size and composition due to their respective dependences.

Overall, the horizontal and meridionial mean gas kinetic energy density, rms velocity and cloud property contours remain reasonably constant throughout our study integration period. This suggests that the horizontal and meridional gas and cloud structures are not significantly varying in time during the ~60 Earth days simulated here, and so further changes to the cloud structure are likely to come from the longer timescale vertical motions. An integration time of ~60 days does not capture the longer vertical settling timescales (>1000 days; Parmentier et al. 2013) of small particles (<0.1 μm) in the upper atmosphere, nor the larger (~1 μm) particles at the clouds base at ~1 bar (Sect. 7).

4. Results

This section presents our results regarding the combined modelling of cloud formation and radiative hydrodynamics for the giant gas planet HD 189733b. We use snapshot results of our simulation at 65 Earth days to illustrate the global cloud formation structures. We split our results into two broad areas, the gas phase properties such as temperature profiles (Sect. 4.1) and horizontal velocity (Sect. 4.2); and the cloud properties such as cloud particle number density structures (Sect. 5.1), mean cloud particle grain sizes (Sect. 5.2) and material composition (Sect. 5.3). Section 5.4 presents the depletion/replenishment of gas phase elements due to cloud formation processes and examines the hydrodynamic transport of elements from dayside to nightside. A global summary of the results and cloud formation physics is presented in Sect. 5.5. Section 6 presents the band by band gas and dust wavelength dependent opacity of the model in order to examine the radiative effects of cloud opacity on the temperature structure of the simulation.

4.1. Global temperature profiles

The local thermodynamic quantities like gas temperature and gas pressure, the local element abundances, determine the local cloud formation processes. The resulting cloud particles affects the local temperatures through their opacity, which in turn is coupled to the local pressure and density through the equation of state and Navier-Stokes equations. We therefore study the local gas temperature which will allow us to develop a global picture of the atmospheric temperature of hot Jupiters like HD 189733b under the influence of cloud formation.

A variation between dayside and nightside is present, with largest gradients in temperatures typically occurring near the terminator regions, most apparent at 0.1, 1 and 10 mbar (Fig. 2). The larger hydrodynamic velocities (super-sonic jet streams) at equatorial regions advects energy density Eastward, resulting in a longitude offset of the temperature maximum by φ ~ 20−40° East compared to the sub-stellar point φ = 0° where the planet receives maximum irradiation from the host star. This is most apparent at 100 mbar and 1 bar in Fig. 2. Differences in temperature (>100 K) between equatorial regions and mid-high latitudes are present in upper atmospheric regions. At depths >5 bar the local gas temperature starts to become more uniform in longitude and latitude. Hottest upper atmosphere regions (Tgas> 1500 K) occur on the dayside at the peak of the temperature inversions (~10 mbar). The coolest regions occur at nightside mid-latitudes with temperatures of ~400 K. These low temperature regions correspond to the large scale, nightside vortex regions where global hydrodynamical motions do not efficiently transport dayside hot gas towards.

Figure 3 displays the 1D (Tgas, pgas) trajectories and shows that the atmosphere contains steep dayside temperature inversions at ~10 mbar. A second, smaller temperature inversion occurs at higher gas pressures ~1 bar for all longitudes at the equator (θ = 0°; Fig. 3, top left). This temperature bump of 100...200 K bump initially develops on the dayside due to a back-warming effect of the larger cloud opacity at ~1 bar (see also Sect. 6). Emission of radiation from hot gas at lower pressure is absorbed by the cloud layer at ~1 bar, resulting in localised heating of the gas. The hot gas is then transported to the nightside by the horizontal winds at this pressure (Fig. 4), resulting in a bump for all equatorial profiles at ~1 bar.

thumbnail Fig. 1

Mean global hydrodynamic and cloud properties as function of local gas pressure during ~60 Earth days of simulation. Top: kinetic energy density log 10Eh,kin [erg cm-3] (left) and global mean rms velocity vrms [m s-1] (right). Bottom: global mean cloud particle number density log 10nd [cm-3] (Eq. (17), left) and drift velocity log 10| vdr | [m s-1] (Eq. (13), right). Parallel contour lines indicate that the global properties of the chosen value in the atmosphere, at that pressure iso-bar, are not changing significantly with time. Note: there is no spin-up period since we use of a well converged continuation simulation of Dobbs-Dixon & Agol (2013) as initial conditions.

Open with DEXTER

Figure 3 also shows the zonal and meridional mean gas temperature Tgas [K] as a function of depth. Zonal mean temperature show how the global temperature structure is changing with latitude and depth. This is useful in order to show global differences between equatorial and mid-high latitude regions for atmospheric properties. Horizontal contours indicate a more uniform variation in temperature in latitude, while vertical contours indicate a greater variation with latitude. Meridional mean temperatures show how the global temperature structure is changing with longitude and depth. This is useful in order to show atmospheric differences between dayside and nightside regions. Horizontal contours indicate a more uniform variation in temperature in longitude, while vertical contours indicate a grater variation with longitude. The highest temperature regions at the upper atmosphere are concentrated at the equator, while a larger (>100 K) difference occurs between equatorial and higher latitude regions. From the meridional mean plot (Fig. 3, bottom right), a stream of hotter gas is present past the φ = 90° terminator at 100 mbar due to hydrodynamic flows advecting energy to the nightside and into deeper regions of higher pressure. The temperature becomes more uniform in longitude and latitude at deeper atmospheric regions >5 bar.

4.2. Atmospheric velocity field

3D RHD simulations provide information about the local and global hydrodynamical behaviour. The dominant hydrodynamical feature of GCM/RHD models of hot Jupiter atmospheres is the formation of an equatorial jet. The super-rotating, circumplanetary jet (Tsai et al. 2014) efficiently advects energy density from day to night near the equatorial regions. This jet forms from the planetary rotation (Rossby waves) coupled with eddies which pump positive angular momentum toward the equator (Showman & Polvani 2011). This section studies the local and global velocity profiles of an atmosphere where cloud formation takes place. 1D profiles of zonal/horizontal velocity at the equator and mid-latitudes are presented in Fig. 4. These show that an upper atmosphere super-sonic jet of velocity >4000 m s-1 at equatorial regions. A significant slow down of horizontal velocity at φ ~ 315° longitude occurs West of the sub-stellar (φ = 0°) point. The maximum zonal velocities of >7000 m s-1 occur on the nightside near the night-day φ = 315° terminator. At mid-latitudes (θ ~ 45°), nightside (φ = 135 ... 270°) regions contain super-sonic counter rotating flows with a velocity of <2000 m s-1. Figure 4 also shows the zonal mean gas velocity vh [m s-1] as at different latitudes and atmospheric pressures. This shows that there is supersonic jet flow at the equator in the West-East direction. Counter rotating flows occur at mid-latitudes with lower horizontal velocities than equatorial regions. Below ~1 bar, the horizontal motions are slower and longitude, latitude uniform until reaching the inner boundary of our computational domain. The overall structural features remain similar to Dobbs-Dixon & Agol (2013).

5. Dynamic mineral clouds in HD 189733b

Giant gas planets like HD 189733b form clouds in their atmospheres from a chemically very rich gas phase. Lee et al. (2015a) have shown that the local thermodynamic conditions suggest that clouds form throughout the whole atmosphere of HD 189733b, although this result is limited by the non-global, 1D mixing approach. A similar conclusions was reach for HD 209458b in a comparison study of both planets (Helling et al. 2016). While Lee et al. (2015a) and Helling et al. (2016) present their results for stationary cloud structures, we now discuss the formation of clouds in a dynamic, time-dependent atmosphere in combination with the 3D atmospheric temperature and velocity fields. The following section shows how cloud properties like number density of cloud particles (Sect. 5.1), cloud particle sizes (Sect. 5.2) and the material composition (Sect. 5.3) develop in and form a dynamic cloud structure in an atmospheres with hydrodynamic jets and temperature inversions.

5.1. Seed formation and cloud particle density

The resultant number density structure of the cloud particles is a combination of the initial nucleation of seed particles and the hydrodynamic motions that transport cloud particles across the globe. Early in the simulation, nucleation begins the cloud formation process with the most efficient nucleation occurring at the colder nightside mid-latitude regions (Fig. 5). Nucleation is a quick processes, and a few minutes/hours into the simulation atomic Ti is too depleted, limiting the nucleation of further cloud particles. The nucleation source term, J(r), for the dust moment equations presented in Sect. 2.2 becomes negligible across the globe. Further evolution of the number density structure of cloud particles is then determined by the hydrodynamical and particle settling motions, rather than further nucleation. Hotter dayside temperature regions at ~10 mbar do not nucleate cloud particles at any point in time during the presently simulated epoch, but hydrodynamic motions transport cloud particles into these regions. Seed particles remain thermally stable in these regions throughout the whole time-span of the present simulation.

thumbnail Fig. 2

Each panel shows Tgas[K] (colour bar) and the velocity field (given by the velocity vector | u | = ) (arrows) at different atmospheric pressure iso-bars at pgas = 0.1, 1, 10, 100 mbar and 1, 10 bar for different φ (longitudes) and θ (latitude). Note: the colour bar scale is different for each plot. The sub-stellar point is located at φ = 0°, θ = 0°.

Open with DEXTER

thumbnail Fig. 3

Top: 1D (Tgas, pgas) trajectories at latitudes θ = 0° (left) , 45° (right). Dayside profiles (solid) show steep temperature inversions at ~10 mbar, especially at higher latitudes. Nightside (dashed) and terminator (dotted) profiles have smaller inversions. Bottom: zonal (left) mean Tgas [K] as a function of atmospheric pressure and meridional (right) mean Tgas [K] as a function of atmospheric pressure. The largest differences in latitudinal temperature contrasts occur from 10 mbar to 1 bar. The temperature is generally isothermal at atmospheric regions at pressures >5 bar.

Open with DEXTER

thumbnail Fig. 4

Top: 1D zonal/horizontal velocity vh [m s-1] trajectories at latitudes θ = 0° (left), 45° (right). Equatorial regions show positive super-sonic flow confined to θ = ± 30° latitudes, with maximum velocities greater than 7000 m s-1 at the upper nightside atmosphere. Negative direction velocities occur at higher latitudes (|θ | > 40°). Bottom: zonal (left) and meridional (right) mean horizontal velocity vh [m s-1] as a function of atmospheric gas pressure. Note: the colour scale bar has been normalised to 0 m s-1. The strongest zonal velocities occur at the equator. Negative flows can be found at latitudes of θ ± 40...80°.

Open with DEXTER

thumbnail Fig. 5

Top: 1D nucleation rate log 10J [cm-3 s-1] trajectories ~12 s into the simulation for latitudes θ = 0° (left), 45° (right). The initial, most efficient nucleation regions occur at pgas< 100 mbar for all atmospheric profiles. Dayside equatorial profiles at φ = 0°, 45° have no nucleation occurring from ~10100 mbar, where the gas temperature is too high for the nucleation process. The greatest magnitude of nucleation of seed particles is at nightside, mid-high latitude regions at ~1 mbar.

Open with DEXTER

thumbnail Fig. 6

Cloud particle number density of grains log 10nd [cm-3] (colour bar) and velocity field () at 0.1, 1, 10, 100 mbar and 1, 10 bar for different φ (longitudes) and θ (latitude). Note: the colour bar scale is different for each plot. The sub-stellar point is located at φ = 0°, θ = 0°. Grains are typically more concentrated at equatorial nightside regions. The number density increases until reaching a maximum near 1 bar, which then gradually falls until the lower computational boundary at ~100 bar. The grains follow the flow patterns in the upper atmosphere, showing a preference to transport cloud particles to nightside equatorial regions. Regions deeper than ~1 bar show a more uniform distribution of cloud particles in latitude and longitude.

Open with DEXTER

thumbnail Fig. 7

Top: 1D cloud particle number density log 10nd [cm-3] trajectories at latitudes θ = 0° (left), 45° (right). The density structure is similar across the longitude and latitude range. Density rises to a maximum of ~105 cm-3 at ~1 bar which contains the thickest and most opaque cloud regions. Bottom: zonal mean (left) and meridional mean (right) of the number density of grains log 10nd [cm-3]. Note: the colour bar scale is different for each plot. The most cloud dense region is from ~100 mbar10 bar which is uniform across the globe. The thinest cloud layers are found at the simulation upper (~0.05 mbar) and lower (~500 bar) computational boundaries.

Open with DEXTER

Figure 6 shows the number density nd [cm-3] of cloud particles after the ~60 Earth simulated days at isobars from 0.1 mbar10 bar. The vertical frictional coupling between dust particles and gas is large enough that cloud particles move with the 3D gas flow efficiently. This means, because we assume horizontal coupling, the highest cloud particle number density occurs near and at the equatorial regions for all atmospheric pressures. Drift velocities are small (vdr ~ 10-4... 0.3 m s-1) throughout the atmosphere (Fig. 1), generally <10% the local vertical gas velocity. This is a purely hydrodynamical effect and implies that the local grain sizes are not large enough to cause a significant de-coupling in the vertical direction. Cloud particle motion is therefore dominated by the horizontal gas velocity (i.e. vdrugas,v<ugas,h). The cloud structure predominantly follows the horizontal velocity structures in the atmosphere. For example, cloud particles entering the equatorial jet stream will typically spend a longer time circulating in these regions due to the lower flux of particle out from the central jet, either from meridional motions or vertical settling. Therefore, after a few days of simulation, after seed particle nucleation has become inefficient, equatorial regions are typically denser by ~0.5 mag compared to mid-high latitudes. This is despite the majority of the seed particle nucleation taking place on nightside high-latitude regions, where the flow speed is typically slower, in the first hours of the simulation. There is a build up of material on approach to the φ = 270° terminator corresponding to regions of highest horizontal velocity (uh ~ 7000 m s-1). This is due to the equatorial jet transporting material quickly around to the nightside and slowing down significantly (Fig. 4) when reaching the night-day terminator. This build up of cloud particles means that the overall flux of particles entering the dayside regions is reduced. Additionally, particles entering the equatorial dayside regions are also transported quickly back onto the nightside by the equatorial jet, increasing the flux of particles towards the nightside. This leads to differences in clouds number density, nd, between the dayside and nightside. The cloud number density at 0.1, 1 and 10 mbar is also slightly (~1%) reduced by downward flows dragging the cloud particles to deeper depths near the day-night φ = 90° terminator and replenished by similar amounts at the upwelling night-day φ = 270° terminator. At deeper regions (>100 mbar), due to the lower vertical gas velocities, the effect of vertical velocity on the number density structure is negligible. At depths greater than 100 mbar there is less difference in number density between nightside-dayside and equatorial-high latitude regions. Figure 7 show 1D trajectory plots and the zonal and meridional mean number density as a function of pressure. These clearly show a thicker cloud layer from ~100 mbar10 bar which is relatively uniform throughout the planet which contains micron sized grains (Fig. 9). On average, the atmosphere is fairly uniform in number density with 1020% differences between equatorial and mid-high latitude regions and comparable difference between nightside and dayside regions. This finding has implications for the cloud opacity which will therefore mainly be affected by the size of the cloud particles and their mixed material composition.

The large scale hydrodynamical motions explain the variety of the cloud number density seen on the dayside/nightside and show that cloud particle structures closely follow the horizontal, meridional and vertical gas dynamics at each atmospheric layer. Efficient nucleation of seed particles occurs at mid-high latitudes on the nightside early in the simulation. The gas dynamics then transports them, over time, to the equatorial regions where most of the cloud particles can be found by the end of the simulation. This result may change if frictional coupling of the cloud particles with the atmospheric gas is altered by a force that specifically acts on cloud particles and causes them to move with udugas in horizontal/meridional directions.

5.2. Cloud particle sizes

thumbnail Fig. 8

Mean grain size a [μm] (colour bar) and velocity field () at 0.1, 1, 10, 100 mbar and 1, 10 bar for different φ (longitudes) and θ (latitude). Note: the colour bar scale is different for each plot. The sub-stellar point is located at φ = 0°, θ = 0°. The smallest grains at each layer typically reside at the equatorial regions. The largest grains are typically found on the nightside and at higher latitudes. Deep purple/blue coloured regions at 1, 10 and 100 mbar contain seed particles of sizes ~0.001 μm

Open with DEXTER

thumbnail Fig. 9

Top: 1D mean grain size log 10a [μm] trajectories at latitudes θ = 0° (left), 45° (right). Cloud particles are typically ~0.1 μm from 0.011 mbar. Seed particles are present between 1100 mbar for many profiles. Bottom: zonal mean (left) and meridional mean (right) average grain size log 10a [μm] as a function of pressure. The smallest grains are found on the dayside at pressures 1100 mbar at the equatorial regions. The grain sizes are typically sub-micron above ~5 bar and micron sizes below this pressure level.

Open with DEXTER

Cloud particle sizes are a direct result of our cloud formation model (Sect. 2.2). Each growth species local growth/evaporation rate (Eq. (16)) determines the local grain size, which is a function of the local number of elements and temperature. It also depends on the sinking/settling of cloud particles of different sizes over time as larger grains sink faster to higher pressure regions. Figure 8 shows the mean grain size a [μm] of cloud particles at pressure isobars from 0.1 mbar10 bar. Dayside regions from 0.1 mbar1 bar contain smaller grains while the nightside contains larger grains. This is most evident from the 0.1 and 10 mbar plots in Fig. 8 where larger grains (>0.1 μm) reside on the nightside while smaller grains (<0.1 μm) are found on the dayside. Asymmetry in grain size between equatorial and mid-high latitudes is also seen, with equatorial regions containing the smallest grains at any given atmospheric pressure and larger grains supported at higher latitudes. At ~10 mbar seed particles (~0.001 μm) reside on the dayside, corresponding to the highest upper atmosphere temperature regions where all other growth species are evaporated. At depths ~1 bar particles grow to 1 μm sizes or larger. At such high densities, the frictional coupling to the gas phase is almost complete, resulting in very small drift velocities (vdr< 0.001 m s-1).

In Fig. 9 the 1D trajectories at equatorial and mid-latitudes show a varied depth dependent grain size. Nightside and mid-latitude terminator regions regions typically contain grain sizes of ~0.1 μm or less down to 1 bar where they grow to micron sizes and above. Dayside and the φ = 90°, 135° equatorial profiles show the presence of seed particles from ~0.1100 mbar. Figure 9 shows the zonal and meridional averaged mean particle size (note: log scale) as a function of pressure. The atmosphere typically contains sub-micron particles down to a pressure level of ~5 bar. The equatorial dayside regions contain the smallest particulates from 0.1100 mbar. The nightside mean particle size is also ~0.51 mag larger than the dayside grains but remain sub-micron at these pressure levels. The largest particles (~0.11 mm) reside at the most dense parts of the atmosphere at gas pressures >10 bar. Gradients (up to 1 order of magnitude) in cloud particle size occur near the φ = 90° terminator, while grain sizes at the φ = 270° terminator remain relatively homogenous with longitude. This suggests that transit spectroscopy (Pont et al. 2013) would sample a variety of cloud particle sizes. The lower temperature regions at the center of vortex regions on the upper atmosphere nightside contain the largest cloud particles at each pressure level, suggesting that these vortex regions can efficiency trap and grow larger particles.

In summary, our model produces a variety of cloud particle sizes dependent on the local thermo-chemical conditions of the atmosphere. A large portion of the hot equatorial dayside contains seed particles of nm size in contrast to the cooler, nightside and mid-latitude regions where grain sizes can be >0.1 μm. We, however, note that the present results are limited to the growth of TiO2[s], SiO[s], SiO2[s], MgSiO3[s], Mg2SiO4[s] materials, which suggest that the present grain sizes to be a lower limit (Sect. 7).

5.3. Cloud material composition

thumbnail Fig. 10

Meridional polar slices of atmospheric cloud properties at the equator (θ = 0°). Top: local gas temperature (Tgas [K]), cloud particle number density (log 10nd [cm-3]) and mean cloud particle grain size (log 10a [μm]). Middle: volume fraction (Vs/Vtot [%]) of the cloud particle composition containing TiO2[s], SiO[s] and SiO2[s]. Bottom: volume fraction (Vs/Vtot [%]) of the cloud particle composition containing MgSiO3[s] and Mg2SiO4[s]. Outer circular values denote longitude at intervals of φ = 45° from the sub-stellar point (φ = 0°). Radial values indicate log 10pgas isobars from 0.1 mbar100 bar. The globe is irradiated from the direction of the colour bars. Note: the size of the annulus is not scaled to planetary radius.

Open with DEXTER

thumbnail Fig. 11

Meridional polar slices of atmospheric cloud properties at mid-latitude (θ = 45°). Top: local gas temperature (Tgas [K]), cloud particle number density (log 10nd [cm-3]) and mean cloud particle grain size (log 10a [μm]). Middle: volume fraction (Vs/Vtot [%]) of the cloud particle composition containing TiO2[s], SiO[s] and SiO2[s]. Bottom: volume fraction (Vs/Vtot [%]) of the cloud particle composition containing MgSiO3[s] and Mg2SiO4[s]. Outer circular values denote longitude at intervals of φ = 45° from the sub-stellar point (φ = 0°). Radial values indicate log 10pgas isobars from 0.1 mbar100 bar. The globe is irradiated from the direction of the colour bars. Note: the size of the annulus is not scaled to planetary radius.

Open with DEXTER

Cloud particles form that are made of a mix of materials that are locally thermally stable. The composition of these material mixes changes depending on the local thermo-chemical conditions that a cloud particle may encounter when being advected due to the presence of a velocity field. Depending on the local temperature and element abundance properties, different solid growth species may dominate the bulk composition compared to others. The local composition of our mixed material cloud particles is therefore dependent on the local growth and evaporation rates as well as the transport of cloud particles and gas phase elements. Figures 10 and 11 show meridional slices of temperature, cloud particle number density, mean grain size and composition at latitudes θ = 0°, 45° respectively. These show a complicated, non-uniform composition structure depending on what material is thermo-chemically favourable at each local atmospheric regions. These plots visualise the interplay between the gas temperature, number density, grain size and composition of the cloud particles.

TiO2[s] is a high-temperature condensate which forms stable clusters that become subsequently more stable with size through homogenous nucleation. TiO2[s] is therefore an efficient seed formation species and will also contribute to the material richness of the grain mantle by surface growth processes. TiO2[s] rich grains (Vs/Vtot ≳ 80%) are generally found between pressures of 1100 mbar on the dayside of planet, corresponding to the hottest regions of the upper atmosphere. These regions primarily consist of near seed particle size (~0.001 μm) cloud particles due to the more volatile materials evaporating off the grain surface. These seed particles also appear on the nightside of the planet from φ = 90° to φ = 135° at ~100 mbar due to the equatorial jet efficiently circulating hot gas to the nightside and to greater depths. At mid-latitudes, pure TiO2[s] grains are only found in regions with the highest local gas temperatures at ~10 mbar, also seed particle sized. These seed particles are thermo-chemically stable. Elsewhere in the atmosphere, TiO2[s] constitutes less than 5% of the grain volume. Other materials grow more efficiently due to the greater element abundance of their constituent elements. Deep atmospheric regions near the lower computational boundary (~500 bar) contain pure TiO2[s] seed particles where other material is thermally unstable.

SiO[s] is typically <5% of the volume fraction in most of the atmosphere. However, it is found in significant volume fractions of >33% at the equatorial regions from φ = 90° to φ = ~ 315° at gas pressures of pgas ~ 0.1−10 mbar. This corresponds to regions of lower gas temperature and density where Mg/Si-growth is unfavourable. SiO[s] can also be found at the hotter and denser inner atmosphere from 10100 bar where SiO[s] contributes 10% to the total volume of the cloud particles.

SiO2[s] dominates (Vs/Vtot ≳ 33%) the dayside equatorial upper atmospheric regions from 0.11 mbar. It is especially dominant at the upper hotter regions from 0.11 mbar where Mg/Si-materials are thermodynamically unfavourable with near 100% composition in some regions. Grain sizes at these regions are ~0.1 μm. At mid-latitudes, SiO2[s] contributes with >10% to the total cloud particle volume at all longitudes and pressure levels, with large volume fractions >50% at dayside pressures of 0.11 mbar.

MgSiO3[s] is perhaps the most interesting species included our models since its optical properties have been used to fit transit spectra, Rayleigh slope observations. We find that it comprises a large amount (Vs/Vtot> 20%) of the grain composition at mid-high latitudes and at all depths, excluding seed particle regions. However, at equatorial regions it can only be found at the deeper, denser parts of the atmosphere from 100 mbar. At equatorial longitudes from φ = 45°−90° at ~1 mbar it can be found to be 1020 % of grain volume.

Mg2SiO4[s] is found to be the most abundant Mg/Si material. Mg2SiO4[s] and MgSiO3[s] can be found in the same regions in the atmosphere, and follow similar trends for their thermal stability. However, Mg2SiO4[s] contributes a larger volume fraction when both materials are thermally stable due to its larger monomer volume. It is the most dominant material at pressures greater than ~500 mbar with grain volumes over 50%.

Overall, a complicated longitude, latitude and depth dependence of the cloud composition across the globe with the gas temperature playing a key role. The thermal instability of the silicate materials at the dayside upper atmosphere regions leads to large volumes of the dayside containing thermally stable, nm-sized TiO2[s] seed particles. An equatorial belt of SiO2[s] and SiO[s] forms due to the different thermo-chemical conditions between mid-latitudes and equatorial regions. Silicate materials such as SiO2[s], MgSiO3[s] and Mg2SiO4[s] are abundant at terminator regions (φ ~ 90°, 270°) probed by transit spectroscopy.

5.4. Non-uniform element abundances

thumbnail Fig. 12

Meridional polar slices of gas phase element abundance log 10εi = ni/n⟨ H ⟩ (ratio i to hydrogen abundance) of top: Ti (blue), second row: O (red), third row: Si (purple) and fourth row: Mg (orange/brown) at θ = 0° (left) and θ = 45° (right), respectively. For reference, solar metallicity of the elements from Asplund et al. (2009) are, Ti: −7.05, O: −3.31, Si: −4.49, Mg: −4.40. Lighter coloured regions indicate a depletion of elements due to the cloud formation processes. Darker coloured regions indicate a local element abundance closer to the initial solar values. Outer circular values denote longitude at intervals of φ = 45° from the sub-stellar point (φ = 0°). Radial values indicate log 10pgas isobars from 0.1 mbar100 bar. The globe is irradiated from the direction of the colour bars. Note: the size of the annulus is not scaled to planetary radius.

Open with DEXTER

To complete our understanding of why certain mineral materials are thermally unstable at certain atmospheric regions, we look at the elemental abundances of elements involved in the cloud formation. The local element abundances determine the gas phase chemical composition which are used to determine the composition of the cloud particles. Several materials can be thermally stable at a particular Tgas-pgas, meaning each of the condensing materials (S> 1) must be considered when calculating the element depletion (Appendix B, Helling & Woitke 2006). If one solid material has reached an equilibrium state (S = 1), this does not imply that all the constituent elements of the material have been condensed. If the initial abundance of elements in the gas phase is low for a particular element (e.g. Ti), then any changes in the elemental abundance due to the cloud formation processes will be greater compared to more abundant elements. Elements are depleted from the gas phase where cloud particle formation has occurred and are replenished where the material that the cloud particles are composed of has evaporated from the grain bulk. The “evaporation window” of a material marks the atmospheric longitude (φ) location where a particular cloud particle material becomes thermally unstable. As cloud particles are transported towards this longitude by circulating gas, the elements from the unstable material replenishes the gas phase. More volatile materials evaporate further away from the sub-stellar point (φ = 0°), while more thermally stable material evaporates closer to the sub-stellar point. Figure 12 shows meridional slices of the gaseous elemental abundance of Ti, O, Si and Mg at θ = 0°, 45°. This shows where elements have been depleted in the atmosphere by the formation of cloud particles or replenished by the evaporation of thermally unstable materials.

Ti is typically depleted by 220 orders of magnitude depending on the location in the atmosphere. The magnitude of this severe depletion is due to the initial solar metallicity low abundance of Ti in the atmosphere. The highest abundance of gas phase Ti elements occurs at the dayside at ~10 mbar where Ti is replenished by the evaporation of TiO2[s] in these regions. Seed particles are not evaporated, however, indicating that the TiO2[s] seed particles remain thermally stable in these regions. The gas temperature is also too high to nucleate seed particles at these regions. The atmosphere never returns to solar metallicity values due to the thermal stability of the TiO2[s] seed particles throughout the 3D radial extent of the simulation boundaries. The evaporation window for TiO2[s] occurs at approximately φ = 315° longitude at 10 mbar. At θ = 0° the equatorial jet carries uncondensed gas towards the nightside where the local temperatures are low enough to allow an efficient surface growth of TiO2[s]. The result is a decreasing Ti-element abundance by >10 mag on the nightside. The growth efficiency of TiO2[s] is less compared to other materials because less Ti is available to condense into TiO2[s] (compare Helling & Woitke 2006).

O, the most abundant element considered here, is depleted by a maximum of 30% throughout the atmosphere. Replenishment of oxygen at θ = 45°, 10 mbar occurs due to the evaporation of O bearing silicate materials.

Si is the least depleted element on the nightside in comparison to Ti and Mg. Upper cooler atmospheric regions where SiO[s] is present tend to be less depleted. The depletion of Si at the equator, θ = 0°, is from SiO2[s] growth. The evaporation window for Si/O material is φ ≈ 300° longitude at 10 mbar. Equatorial jets resupply the nightside with Si elements from evaporating dayside silicate materials.

Mg is typically more depleted on nightside regions indicating the efficient formation of Mg/Si-material. Gas containing Mg is transported to the nightside by the equatorial jet, where the lower temperatures allow the Mg/Si-material to be thermally stable. Severe depletion occurs at ~110 mbar where Mg/Si-material surface growth occurs. The evaporation window for Mg/Si-material is φ 280° longitude at 10 mbar.

Due to the global dynamics of the upper atmosphere, any replenished elements from evaporating material at mid-latitudes gets funnelled towards the equator. This hydrodynamic preference is how the mid-latitude, nightside regions contain some of the most element depleted areas of the globe with >10 mag depletion. Si/O materials are also dominant at the equator due to these hydrodynamic motions, compared to Mg/Si despite both elements have similar initial solar abundances. At mid-latitude regions the efficient growth of Mg2SiO4[s] removes 1 Mg atom more from the gas phase compared to Si atoms for each Mg2SiO4[s] monomer condensed onto the cloud particle. The meridional motions then funnel gas from these mid-latitudes regions towards the equatorial jet, leading to a greater gaseous abundance of Si at the equator compared to Mg. This is seen in the different ranges of colour bar for Si and Mg in Fig. 12 where Si is typically more abundant by 510 mag. The supersaturation ratio (Eq. (15)) for SiO[s] and SiO2[s] materials are therefore larger than Mg2SiO4[s] and MgSiO3[s], since there are less gas phase Mg atoms to condense when grains are transported to the equator. Mg2SiO4[s] and MgSiO3[s] constitute less of the the grain volume as the supersaturation of SiO[s] and SiO2[s] increases in these regions while Mg2SiO4[s] and MgSiO3[s] decreases. That is, Mg/Si material can become thermo-chemically unstable at equatorial regions while Si/O material remains stable. SiO[s] and SiO2[s] are more efficient growth materials than Mg2SiO4[s] and MgSiO3[s] at the equatorial thermo-chemical conditions and become the most dominant material compositions at the equator. Over time, this leads to an upper atmosphere equatorial band of SiO[s] and SiO2[s]. A longitude, latitude and height inhomogeneous element abundance structure develops as the evaporation and growth windows are different for each species as they travel around the atmosphere.

In summary, the atmosphere is depleted of elements due the cloud formation processes. Local regions of cloud material evaporation on the dayside replenish the gas phase of elements. Replenishment of the nightside of elements is governed through the equatorial jet, which transports uncondensed material from the dayside to the nightside. The gas phase elements are then depleted by growth of material at the cooler regions. A consequence of the dayside thermal instability and hydrodynamic transport of material is that the φ = 90° terminator region is replenished in metallic elements while the φ = 270° terminator is severely depleted. The atmosphere never returns to solar element abundances for those elements involved in cloud formation.

5.5. Summary of dynamic cloud formation results

The previous sections presented each of our cloud formation results individually. In this section, we examine the physics of our results as a whole and explain the cloud formation process on a global scale.

The nucleation of seed particles early in the simulation begins the cloud formation process. After a few minutes/hours of simulation, the rate of seed particle formation becomes negligible throughout the entire atmosphere due to the depletion of Ti elements from the gas phase. Meanwhile, cloud particles that grow 1 μm begin to settle rapidly from the upper atmosphere to deeper layers (~1 bar), to their pressure supported levels. Settling of cloud particles results in a globally uniform cloud layer of maximum number density at ~1 bar. These settling grains deplete elements from the upper atmosphere. This leaves sub-micron cloud particles in the upper atmosphere from pgas< 1 bar. Cloud particles are transported by meridional gas motions towards the equatorial jet where the largest number density nd of cloud particles can be found. Cloud particles then follow a circulation cycle as they are transported across the globe, dependent on the local temperature and element abundance conditions. This large scale cloud formation circulation cycle can be summarised as follows:

  • As cloud particles are transported towards the night-dayterminator(φ = 270°), the increase in gas temperature leads material to become thermally unstable. The most volatile material evaporates first, while the more stable material evaporates closer to the sub-stellar point (φ = 0°). TiO2 seed particles (~0.001 μm) remain thermally stable. The local gas phase is replenished in elements from the evaporating material.

  • The equatorial jet transports the element replenished gas phase and the thermally stable seed particles (~0.001 μm) past the sub-stellar point (φ = 0°) and to the φ = 90° day-night terminator. This replenishes the nightside regions near the day-night terminator of uncondensed elements.

  • From longitudes φ = 90−180°; due to the colder nightside temperatures and the replenishment of gas phase elements and transport of seed particles from the equatorial jet, silicate material is thermally stable and condenses onto the seed particles. This depletes the gas phase of elements on the nightside. Due to this cloud formation, particles are ~0.010.1 μm in these regions.

  • Cloud particles remain thermally stable as they are transported through longitudes φ = 180−270°. The gas phase is severely depleted of elements but the cooler temperatures keep the cloud material thermally stable. The cloud particles are then transported to the night-day terminator (φ = 270°) and the cloud formation cycle repeats.

A characteristic of this element cycle is that the amount of metallic elements returned to the gas phase by dayside grain evaporation is not enough to grow the grains to large sizes (>1 μm) when they are thermally stable on the nightside. Over time, an upper atmosphere equatorial band of SiO2[s] and SiO[s] rich grains develops due to greater number of uncondensed Si atoms available at the equator, while Mg atoms condense at higher latitudes.

6. Cloud/gas opacity and radiative effects of clouds

thumbnail Fig. 13

Cloud (top), gas (middle) and total (bottom) opacities at the sub-stellar point (left) and φ = 0°, θ = 45° (right) at the centre of each of the wavelength bands used in the RHD radiative transfer scheme. The thick dashed line shows the parameterised cloud opacity used in Dobbs-Dixon & Agol (2013), Eq. (31). The addition of the the size and composition dependent cloud opacity results in a more inhomogeneous opacity structure.

Open with DEXTER

thumbnail Fig. 14

Stellar heating rate S [erg cm-3 s-1] (Eq. (4)) at the equator and θ = 45° latitude, φ = 0°, 45°, 315° longitude as a function of pressure. The heating rate steadily rises up to a maximum at ~10 mbar corresponding to the highest temperature regions at the upper atmosphere. The heating rate drops off rapidly at 100 mbar.

Open with DEXTER

From Sect. 4.1 the temperature structure of the planet is affected by the presence of cloud particles. The cloud particles add an inhomogeneous opacity distribution to the atmosphere, altering the local radiation fields. A specific feature due to the presence of cloud particles is the equatorial regions temperature bump at ~1 bar of >100 K. Figure 13 shows the cloud, gas and total opacity at the wavelength midpoints of the wavelength opacity bands in Showman et al. (2009), applied in this paper, at the sub-stellar point and φ = 0°, θ = 45°. Figure 13 also shows the parameterised cloud opacity used in Dobbs-Dixon & Agol (2013; black dashed line), (31)where κgrey = 0.035 cm2 g-1 and κRS = 0.6 cm2 g-1. The size and composition dependent cloud opacity that results from the present simulations is typically lower than the Dobbs-Dixon & Agol (2013) opacity in the upper dayside atmosphere (pgas< 1 bar) but larger and greyer, i.e. less strongly wavelength dependent, at all wavelengths from pgas> 1.

Figure 14 shows the stellar heating rate S (Eq. (4)) at the sub-stellar point, θ = 45° at φ = 45°, 315° longitudes. The peak of energy deposition by stellar photons occurs at ~10 mbar, where the highest dayside temperatures occur. The stellar energy deposition drops off rapidly at 100 mbar where some of the coolest dayside temperatures can be found. Compared to a self-consistent gas phase opacity simulation (Fig. 1, Tsai et al. 2014), the peak of stellar energy deposition is at the same pressure level ~10 mbar. However, the peak of the heating at the sub-stellar point for the gas phase opacity simulation is ~3 erg cm-3 s-1 less and there is a more gradual drop off in heating to ~1 bar. The Tsai et al. (2014) simulation was found to be too cold when compared to the observations (Dobbs-Dixon & Agol 2013). Our microphysical cloud structure maintains similar stellar heating regions seen in Dobbs-Dixon & Agol (2013), which suggests that a cloud opacity (parameterised or microphysical) pushes the stellar energy deposition further upward on the dayside atmosphere. The upper atmosphere gas temperature is typically cooler on the dayside and nightside with the microphysical cloud model compared to the parameterised cloud in Dobbs-Dixon & Agol (2013). This indicates that the lower cloud opacity in the microphysical model allows the gas to cool more efficiently in these regions than in the Dobbs-Dixon & Agol (2013) simulation.

Amundsen et al. (2014) suggest that the Planck averaged gas opacities used in the current study can lead to greater uncertainties in the stellar heating rate compared to other methods. The addition of cloud opacity may reduce this error by muting or washing out the rich molecular lines when the cloud opacity approaches the gas opacity, which can be seen in Fig. 13. However, in regions of low cloud opacity (e.g. seed particle regions) the results of Amundsen et al. (2014) suggest that the maximum of the stellar heating rate may occur in deeper atmospheric layers. We hypothesise that this would lead to a smaller or larger vertical extension of the seed particle region on the dayside, depending if this change decreased or increased the temperature at the seed particle region boundaries. Energy deposited at greater depth would be advected more efficiently in the vertical and horizontal directions, which may impact the overall trends of the cloud particles.

Since the stellar energy deposition is negligible at pgas> 100 mbar, the Tgas> 100 K temperature bump seen at the equator (Fig. 3) cannot be due to stellar heating. A backwarming effect due to the presence of the opaque cloud base at ~1 bar occurs. This backwarming was not seen in Dobbs-Dixon & Agol (2013), suggesting that the increased cloud opacity in these regions (Fig. 13) is responsible for this feature. The gas irradiated by the host star at ~10 mbar radiates with a Planck function B(λ, T) peak at λ ~ 1−2μm. This emitted radiation is then absorbed deeper in the atmosphere where the cloud opacity at ~1 bar is largest at these infrared wavelengths, heating the local gas phase. The backwarming effect is not seen at mid-latitudes due to the larger grain sizes of the cloud at pgas< 1 bar, producing larger infrared opacity (by up to 3 mag) compared to the equator. The remitted infrared at mid-latitude regions is then absorbed closer to the temperature peak at 10 mbar, which results in a slightly flatter temperature inversion. Therefore, the 100200 K temperature bump is not seen for mid-high latitude regions. Similar backwarming effects due to cloud particle and/or gaseous opacity have been investigated in previous studies (e.g. Helling et al. 2000; Tsuji 2002; Burrows et al. 2006; Witte et al. 2009; Heng et al. 2012; Heng & Demory 2013).

7. Discussion

Our results suggest that the atmosphere of HD 189733b contains a silicate mineral cloud component in line with interpretation of transit spectra and albedo observations (e.g. Lecavelier Des Etangs et al. 2008; Sing et al. 2011b, 2016; Pont et al. 2013; Evans et al. 2013; Wakeford & Sing 2015). Lecavelier Des Etangs et al. (2008) suggest that a sub-micron (~0.01...0.1 μm) grain size of silicate composition such as MgSiO3[s] located at local gas pressures of ~10-6...10-3 bar can fit the optical wavelength Rayleigh slope. The results of our model suggest that transit spectra observations would sample a variety of grain sizes at different pressure levels within the size and pressure ranges suggested by Lecavelier Des Etangs et al. (2008) and Wakeford & Sing (2015). Additionally, the composition of the grains changes with longitude and latitude, meaning these observations would also sample differences in cloud particle composition.

We neglect possible condensation of more volatile material such as ZnS[s], KCl[s] and Na2S[s]. These materials have been used for modelling cooler objects such as GJ 1214b (Charnay et al. 2015b) and T Brown Dwarfs (Morley et al. 2012), where the more stable species considered in this study are found deeper (below τ ~ 1) in the atmosphere. However, if the atmosphere can efficiently mix solid/gas material upwards, we can expect the more stable condensates to also be present at high altitudes for these cooler objects. Sing et al. (2015) suggest that Na2S[s] and KCl[s] condensation could be responsible for the sub-solar Na to K ratio observed on hot Jupiter WASP-31b (Teq = 1575 K). SiO[s], the most volatile species in the current set-up, is only abundant at specific regions on the nightside which are thermo-chemically stable for it. The atmosphere is generally highly depleted (>10 orders of magnitude less than solar) of the Ti, Si, Mg elements which take part in the cloud formation, however, longitude, latitude and depth differences in atomic abundance are present. From our model, the dayside-nightside terminator region is replenished of elements by the equatorial jet after material has evaporated at the hottest dayside regions. We suggest from the presented results, that Na2S[s] and/or KCl[s] condensation could occur on the nightside and deplete Na/K on the nightside-dayside terminator boundary for both HD 189733b and WASP-31b. Even if Na2S[s] and/or KCl[s] are a minor component of the total dust volume, the condensation of the materials can cause a large decrease in elemental abundance, as discussed in Helling et al. (2016). The grains would evaporate their Na/K content once they travel to the dayside which would replenish the Na/K atomic abundance for the dayside-nightside terminator. Additionally, the different thermo-chemical kinetics of Na2S[s] and KCl[s] could lead to latitudinal variance, similar to our equatorial band of Si/O and mid-latitude Mg/Si material dominated regions.

thumbnail Fig. 15

Snapshot horizontally and meridionally averaged iso-bars of the time dependent changes in nd (left) and εMg (right) due to the vertical advection. εMg is given in the range 1−100 bar to avoid the skewing of the global averages by dayside particle evaporation in the upper atmosphere. Settling of grains by the drift velocity is taken into account for the nd plot.

Open with DEXTER

As noted in Mayne et al. (2014), for most GCM/RHD modelling of hot Jupiters, the deep atmosphere (p ≳ 1 bar) takes longer to reach a steady state due to the slow (t> 1200 days) momentum exchange between the lower and upper layers. There is evidence from Fig. 1 that the velocity structure is still evolving slowly at these deep regions from the effect of the cloud opacity. Due to the added cost of the microphysical cloud model, the effect of the upper and lower atmospheric cloud opacity on the dynamics of deeper regions may take many more months or years of simulation, beyond the scope of this early investigation.

The mixing of replenished gas material upward from the deeper depths (~100 bar) where cloud particles evaporate their volatile contents is also expected to occur on the momentum exchange timescales suggested by Mayne et al. (2014); as it is this timescale where the information of the gaseous elemental content is exchanged between the deep and upper atmosphere. The replenishment rate over 1 scale height can be approximated by the mixing timescale τmix at these depths. In Lee et al. (2015a) we estimated that the mixing timescale would be on the order of τmix ~ 108 s at ~100 bar. The replenishment rate for Mg abundances at these depths can be estimated from by τmix, /τmix ≈ 4 × 10-13εMg s-1; and τmix = 108 s. This value would be many times smaller as the mixing material travels several atmospheric scale heights before reaching the upper atmosphere. Agúndez et al. (2014) suggest that the mixing timescale may be on the order τmix ~ 109 s at these depths for HD 189733b, calculated from the GCM mixing tracing method of Parmentier et al. (2013). To illustrate this point, Fig. 15 shows snapshot horizontal and meridional mean dnd/dt and dεMg/dt due to vertical advection at gas pressure iso-bars. The small magnitudes of these changes compared to the absolute values suggest that vertical advection may not significantly alter upper atmosphere cloud particle results during the epoch of the simulation discussed here. Longer integration times (>1000 days) will be required to better understand the effect of deeper mixing of gas phase elements and the settling of smaller particles. However, the results of Parmentier et al. (2013) suggest, for HD 209458b, that sub-micron sized cloud particles may remain present in the upper atmospheric (>1bar) layers over longer timescales. Overall, this is in contrast to more convective atmospheres, e.g. Brown Dwarfs, where τmix is estimated to be ~300 s at ~10 bar (Woitke & Helling 2004), which increases the resupply rate of elements to the upper atmosphere.

Through our modelling we have shown that mineral cloud particles can survive for many advective timescales in the atmosphere of HD 189733b. Cloud particles are continually transported between warmer and cooler regions by either horizontal motions from warm dayside to cool nightside, or from vertical flows from the cooler upper atmosphere to the warmer inner atmosphere and back again. Helling & Rietmeijer (2009) show that SiO2[s], MgSiO3[s] and Mg2SiO4[s] can crystallise efficiently in supersaturated atmospheric conditions similar to that of hot Jupiters, especially if grains are transported from hotter to cooler regions and vice-versa. In addition, since the cloud particles can survive for long periods in the atmosphere, the internal structure of the grains have time to rearrange into a crystalline structure. This would change the optical properties of the grains, where they would appear more crystalline compared to amorphous grains. This may increase the wavelength dependent back-scattering of photons.

Our presented cloud modelling also lends weight to a patchy cloud scenario to explain the Westward offsets in optical phase-curves of some Kepler hot Jupiter exoplanets (Heng & Demory 2013; Hu et al. 2015; Esteves et al. 2015; Shporer & Hu 2015). The evaporation window of the reflecting silicate materials in this study occurs φ = 45−90° longitude Westward from the sub-stellar point, leading to smaller grain sizes, a changed material composition and a large drop in cloud opacity across large areas of the mid and Eastern dayside. Due to the different atmospheric properties of the Kepler planets, should this evaporation window occur closer to the sub-stellar point, the difference in reflectivity of the East and West cloud properties would be become more noticeable. Our current study therefore provides a microphysical basis to the analysis in Oreshenko et al. (2016) and Parmentier et al. (2016). By examining the kinetics of the cloud formation with radiative feedback effects for the Kepler planets, a detailed picture of the composition and size of the optical backscattering material can be calculated. The variability in the reflectivity of Kepler-7b has been shown to be small over the Kepler observing period (Demory et al. 2013), suggesting that a long-term, stable cloud particle hydrodynamic circulation and thermo-chemical cycle is a likely possibility for these atmospheres.

7.1. Comparison to 1D results

In this section we compare our 3D coupled model to our previous 1D post-processing, non-global approach in Lee et al. (2015a). In the current study we have investigated 5 cloud species (TiO2[s], SiO[s], SiO2[s], MgSiO3[s], MgSiO4[s]), neglecting 7 materials (Fe[s], FeO[s], FeS[s], Fe2O3[s], CaTiO3[s], MgO[s], Al2O3[s]) which were included in previous modelling efforts (Helling & Woitke 2006; Helling et al. 2008, 2016; Lee et al. 2015a). The previous 1D study Lee et al. (2015a) showed that some of the neglected species can have significant volume fractions in certain parts of the atmosphere. Al2O3[s] and CaTiO3[s] were found to be thermally stable at regions of the hot dayside, this is likely to be the same in the 3D RHD case. This suggests that the grain sizes could be underestimated for these regions in the RHD model. However, Al2O3[s] and CaTiO3[s] are typically not efficient growth species and are unlikely to grow the grain to significantly larger (>0.01 μm) sizes. The addition of Fe[s] will effect the thermal stability of the cloud particles, grain sizes and drift velocity and therefore the local degree of element depletion, especially in the deeper regions where other species are less thermally stable compared to Fe[s]. A richer chemical composition can be expected when additional high-temperature condensates are included.

Lee et al. (2015a) generally reproduces the regions of efficient nucleation, growth and evaporation compared to the 3D RHD model, indicating that the chemical processes are accurately captured by the 1D models. The influence of dynamics on the specific cloud properties is large however, which leads to differences in predicted grain sizes between the two approaches. The time dependent settling of ~1 μm grains to their pressure supported regions near 1 bar results in a higher cloud opacity deeper in the atmosphere for the RHD model compared to the 1D models.

In more general terms, our 1D cloud simulations, Drift, are a valuable analysis tool because they are fast at providing a stationary solution to cloud properties with a substantial degree of chemical details. The 3D RHD simulations with our cloud formation module, are time-consuming but allows us to resolve the time and spacial evolution of a cloud-forming hot Jupiter atmosphere. At present, the time evolution of the current study has focused on the first 60 Earth days and the spacial resolution is limited to that of the RHD/GCM cells. Long-term studies which address the limitations of the current implementation are under development.

8. Summary and conclusions

We have developed a 3D kinetic cloud formation module for exoplanet RHD/GCM simulations. Through our coupled RHD and cloud model we have shown that HD 189733b has extremely favourable thermodynamic and hydrodynamic conditions for efficient cloud formation and growth. We lend weight to previous interpretations of observations of a thick mineral cloud component containing sub-micron sized particles in the upper atmosphere. The interplay between the hydrodynamical motions and the cloud formation produces an inhomogeneous opacity structure which has effects on the global atmospheric conditions. A summary of key results include:

  • Grain sizes are sub-micron at atmospheric pressures andterminator regions probed by transmission spectroscopy.

  • Silicate materials are thermally stable in regions probed by transmission spectroscopy.

  • Grain sizes, number density and opacity of cloud particles are non uniform across the globe with significant differences in longitude, latitude and depth.

  • Due to the global elemental depletion from cloud formation, the equatorial regions are dominated by SiO[s] and SiO2[s] while mid-upper latitudes mostly contain MgSiO3[s] and Mg2SiO4[s].

  • Hydrodynamic motions primarily govern the global distribution of cloud particles, transporting cloud particles towards nightside equatorial regions.

  • The existence of the clouds as well as the particle sizes and material mixes are the result of the kinetic cloud formation processes.

  • The atmosphere is severely depleted (10 orders of magnitude) of elements Ti, Si, Mg used in the cloud formation theory while O is only depleted by ~30%.

  • Mid-high latitude nightside regions are not efficiently replenished in elements and contain the most reduced gas phase elemental abundances.

  • Thermally unstable materials on the dayside replenish elements to the gas phase; these uncondensed elements are then transported to the nightside via the equatorial jet.

  • Maximum thermal stability of the cloud particles are found at the coolest parts of the nightside inside the large scale vortex regions.

We improve on our previous 1D approach by developing an atmospheric RHD model coupled to a cloud formation module, including self consistent opacity feedback, element mixing, cloud particle transport and settling. We emphasise that we do not rely on any mixing parameterisation of gaseous/solid material (e.g. by use of Kzz), assumptions about the grain sizes or particle size distributions. Through this early study, we have demonstrated that our kinetic cloud formation model is well suited to be applied to 3D hydrodynamic studies of exoplanet atmospheres. However, long-term studies are required to address the limitations of the current implementation of the cloud module. Potential future applications of our model are 3D Brown Dwarf atmosphere simulations such as those in Showman & Kaspi (2013), Zhang & Showman (2014). The model can be extended for warm-Neptune studies with the addition of more volatile cloud species such as ZnS[s], KCl[s] and Na2S[s] to the chemical scheme.


1

Using the NYU Abu Dhabi HPC cluster, BuTinah.

Acknowledgments

We thank the anonymous referee for insightful comments and suggestions which improved the manuscript substantially. G.L. and Ch.H. highlight the financial support of the European community under the FP7 ERC starting grant 257431. Our local HPC computational support at Abu Dhabi and St Andrews is highly acknowledged. Most plots were produced using the community open-source Python packages Matplotlib, SciPy and AstroPy. We thank J. Blecic and current/former members of the LEAP team for insightful discussions and suggestions.

References

All Figures

thumbnail Fig. 1

Mean global hydrodynamic and cloud properties as function of local gas pressure during ~60 Earth days of simulation. Top: kinetic energy density log 10Eh,kin [erg cm-3] (left) and global mean rms velocity vrms [m s-1] (right). Bottom: global mean cloud particle number density log 10nd [cm-3] (Eq. (17), left) and drift velocity log 10| vdr | [m s-1] (Eq. (13), right). Parallel contour lines indicate that the global properties of the chosen value in the atmosphere, at that pressure iso-bar, are not changing significantly with time. Note: there is no spin-up period since we use of a well converged continuation simulation of Dobbs-Dixon & Agol (2013) as initial conditions.

Open with DEXTER
In the text
thumbnail Fig. 2

Each panel shows Tgas[K] (colour bar) and the velocity field (given by the velocity vector | u | = ) (arrows) at different atmospheric pressure iso-bars at pgas = 0.1, 1, 10, 100 mbar and 1, 10 bar for different φ (longitudes) and θ (latitude). Note: the colour bar scale is different for each plot. The sub-stellar point is located at φ = 0°, θ = 0°.

Open with DEXTER
In the text
thumbnail Fig. 3

Top: 1D (Tgas, pgas) trajectories at latitudes θ = 0° (left) , 45° (right). Dayside profiles (solid) show steep temperature inversions at ~10 mbar, especially at higher latitudes. Nightside (dashed) and terminator (dotted) profiles have smaller inversions. Bottom: zonal (left) mean Tgas [K] as a function of atmospheric pressure and meridional (right) mean Tgas [K] as a function of atmospheric pressure. The largest differences in latitudinal temperature contrasts occur from 10 mbar to 1 bar. The temperature is generally isothermal at atmospheric regions at pressures >5 bar.

Open with DEXTER
In the text
thumbnail Fig. 4

Top: 1D zonal/horizontal velocity vh [m s-1] trajectories at latitudes θ = 0° (left), 45° (right). Equatorial regions show positive super-sonic flow confined to θ = ± 30° latitudes, with maximum velocities greater than 7000 m s-1 at the upper nightside atmosphere. Negative direction velocities occur at higher latitudes (|θ | > 40°). Bottom: zonal (left) and meridional (right) mean horizontal velocity vh [m s-1] as a function of atmospheric gas pressure. Note: the colour scale bar has been normalised to 0 m s-1. The strongest zonal velocities occur at the equator. Negative flows can be found at latitudes of θ ± 40...80°.

Open with DEXTER
In the text
thumbnail Fig. 5

Top: 1D nucleation rate log 10J [cm-3 s-1] trajectories ~12 s into the simulation for latitudes θ = 0° (left), 45° (right). The initial, most efficient nucleation regions occur at pgas< 100 mbar for all atmospheric profiles. Dayside equatorial profiles at φ = 0°, 45° have no nucleation occurring from ~10100 mbar, where the gas temperature is too high for the nucleation process. The greatest magnitude of nucleation of seed particles is at nightside, mid-high latitude regions at ~1 mbar.

Open with DEXTER
In the text
thumbnail Fig. 6

Cloud particle number density of grains log 10nd [cm-3] (colour bar) and velocity field () at 0.1, 1, 10, 100 mbar and 1, 10 bar for different φ (longitudes) and θ (latitude). Note: the colour bar scale is different for each plot. The sub-stellar point is located at φ = 0°, θ = 0°. Grains are typically more concentrated at equatorial nightside regions. The number density increases until reaching a maximum near 1 bar, which then gradually falls until the lower computational boundary at ~100 bar. The grains follow the flow patterns in the upper atmosphere, showing a preference to transport cloud particles to nightside equatorial regions. Regions deeper than ~1 bar show a more uniform distribution of cloud particles in latitude and longitude.

Open with DEXTER
In the text
thumbnail Fig. 7

Top: 1D cloud particle number density log 10nd [cm-3] trajectories at latitudes θ = 0° (left), 45° (right). The density structure is similar across the longitude and latitude range. Density rises to a maximum of ~105 cm-3 at ~1 bar which contains the thickest and most opaque cloud regions. Bottom: zonal mean (left) and meridional mean (right) of the number density of grains log 10nd [cm-3]. Note: the colour bar scale is different for each plot. The most cloud dense region is from ~100 mbar10 bar which is uniform across the globe. The thinest cloud layers are found at the simulation upper (~0.05 mbar) and lower (~500 bar) computational boundaries.

Open with DEXTER
In the text
thumbnail Fig. 8

Mean grain size a [μm] (colour bar) and velocity field () at 0.1, 1, 10, 100 mbar and 1, 10 bar for different φ (longitudes) and θ (latitude). Note: the colour bar scale is different for each plot. The sub-stellar point is located at φ = 0°, θ = 0°. The smallest grains at each layer typically reside at the equatorial regions. The largest grains are typically found on the nightside and at higher latitudes. Deep purple/blue coloured regions at 1, 10 and 100 mbar contain seed particles of sizes ~0.001 μm

Open with DEXTER
In the text
thumbnail Fig. 9

Top: 1D mean grain size log 10a [μm] trajectories at latitudes θ = 0° (left), 45° (right). Cloud particles are typically ~0.1 μm from 0.011 mbar. Seed particles are present between 1100 mbar for many profiles. Bottom: zonal mean (left) and meridional mean (right) average grain size log 10a [μm] as a function of pressure. The smallest grains are found on the dayside at pressures 1100 mbar at the equatorial regions. The grain sizes are typically sub-micron above ~5 bar and micron sizes below this pressure level.

Open with DEXTER
In the text
thumbnail Fig. 10

Meridional polar slices of atmospheric cloud properties at the equator (θ = 0°). Top: local gas temperature (Tgas [K]), cloud particle number density (log 10nd [cm-3]) and mean cloud particle grain size (log 10a [μm]). Middle: volume fraction (Vs/Vtot [%]) of the cloud particle composition containing TiO2[s], SiO[s] and SiO2[s]. Bottom: volume fraction (Vs/Vtot [%]) of the cloud particle composition containing MgSiO3[s] and Mg2SiO4[s]. Outer circular values denote longitude at intervals of φ = 45° from the sub-stellar point (φ = 0°). Radial values indicate log 10pgas isobars from 0.1 mbar100 bar. The globe is irradiated from the direction of the colour bars. Note: the size of the annulus is not scaled to planetary radius.

Open with DEXTER
In the text
thumbnail Fig. 11

Meridional polar slices of atmospheric cloud properties at mid-latitude (θ = 45°). Top: local gas temperature (Tgas [K]), cloud particle number density (log 10nd [cm-3]) and mean cloud particle grain size (log 10a [μm]). Middle: volume fraction (Vs/Vtot [%]) of the cloud particle composition containing TiO2[s], SiO[s] and SiO2[s]. Bottom: volume fraction (Vs/Vtot [%]) of the cloud particle composition containing MgSiO3[s] and Mg2SiO4[s]. Outer circular values denote longitude at intervals of φ = 45° from the sub-stellar point (φ = 0°). Radial values indicate log 10pgas isobars from 0.1 mbar100 bar. The globe is irradiated from the direction of the colour bars. Note: the size of the annulus is not scaled to planetary radius.

Open with DEXTER
In the text
thumbnail Fig. 12

Meridional polar slices of gas phase element abundance log 10εi = ni/n⟨ H ⟩ (ratio i to hydrogen abundance) of top: Ti (blue), second row: O (red), third row: Si (purple) and fourth row: Mg (orange/brown) at θ = 0° (left) and θ = 45° (right), respectively. For reference, solar metallicity of the elements from Asplund et al. (2009) are, Ti: −7.05, O: −3.31, Si: −4.49, Mg: −4.40. Lighter coloured regions indicate a depletion of elements due to the cloud formation processes. Darker coloured regions indicate a local element abundance closer to the initial solar values. Outer circular values denote longitude at intervals of φ = 45° from the sub-stellar point (φ = 0°). Radial values indicate log 10pgas isobars from 0.1 mbar100 bar. The globe is irradiated from the direction of the colour bars. Note: the size of the annulus is not scaled to planetary radius.

Open with DEXTER
In the text
thumbnail Fig. 13

Cloud (top), gas (middle) and total (bottom) opacities at the sub-stellar point (left) and φ = 0°, θ = 45° (right) at the centre of each of the wavelength bands used in the RHD radiative transfer scheme. The thick dashed line shows the parameterised cloud opacity used in Dobbs-Dixon & Agol (2013), Eq. (31). The addition of the the size and composition dependent cloud opacity results in a more inhomogeneous opacity structure.

Open with DEXTER
In the text
thumbnail Fig. 14

Stellar heating rate S [erg cm-3 s-1] (Eq. (4)) at the equator and θ = 45° latitude, φ = 0°, 45°, 315° longitude as a function of pressure. The heating rate steadily rises up to a maximum at ~10 mbar corresponding to the highest temperature regions at the upper atmosphere. The heating rate drops off rapidly at 100 mbar.

Open with DEXTER
In the text
thumbnail Fig. 15

Snapshot horizontally and meridionally averaged iso-bars of the time dependent changes in nd (left) and εMg (right) due to the vertical advection. εMg is given in the range 1−100 bar to avoid the skewing of the global averages by dayside particle evaporation in the upper atmosphere. Settling of grains by the drift velocity is taken into account for the nd plot.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.