Dust in brown dwarfs and extra-solar planets. VII. Cloud formation in diffusive atmospheres
AAstronomy & Astrophysics manuscript no. Di ff uDrift c (cid:13) ESO 2019November 12, 2019
Dust in brown dwarfs and extra-solar planets
VII. Cloud formation in diffusive atmospheres
Peter Woitke , , Christiane Helling , , , and Ophelia Gunn Centre for Exoplanet Science, University of St Andrews, St Andrews, UK SUPA, School of Physics & Astronomy, University of St Andrews, St Andrews, KY16 9SS, UK SRON Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, NL SUPA, School of Physics and Astronomy, University of Edinburgh, Edinburgh, EH9 3JZ, UKReceived 10 / / / / ABSTRACT
The precipitation of cloud particles in brown dwarf and exoplanet atmospheres establishes an ongoing downward flux of condensableelements. To understand the e ffi ciency of cloud formation, it is therefore crucial to identify and to quantify the replenishment mech-anism that is able to compensate for these local losses of condensable elements in the upper atmosphere, and to keep the extrasolarweather cycle running. In this paper, we introduce a new cloud formation model by combining the cloud particle moment methodof Helling & Woitke with a di ff usive mixing approach, taking into account turbulent mixing and gas-kinetic di ff usion for both gasand cloud particles. The equations are of di ff usion-reaction type and are solved time-dependently for a prescribed 1D atmosphericstructure, until the model has relaxed toward a time-independent solution. In comparison to our previous models, the new hot Jupitermodel results ( T e ff ≈ g =
3) show fewer but larger cloud particles which are more concentrated towards the cloud base. Theabundances of condensable elements in the gas phase are featured by a steep decline above the cloud base, followed by a shallower,monotonous decrease towards a plateau, the level of which depends on temperature. The chemical composition of the cloud particlesalso di ff ers significantly from our previous models. Due to the condensation of specific condensates like Mg SiO [s] in deeper layers,certain elements, such as Mg, are almost entirely removed from the gas phase early. This leads to unusual (and non-solar) elementratios in higher atmospheric layers, which then favours the formation of SiO[s] and SiO [s], for example, rather than MgSiO [s].Such condensates are not expected in phase-equilibrium models that start from solar abundances. Above the main silicate cloud layer,which is enriched with iron and metal oxides, we find a second cloud layer made of Na S[s] particles in cooler models ( T e ff (cid:47) Key words. planets and satellites: atmospheres – planets and satellites: composition – brown dwarfs – astrochemistry – di ff usion
1. Introduction
The number of confirmed extrasolar planets has reached morethan 4000, but only a hand-full of them can be studied in detail(see e.g. Nikolov et al. 2016; Huitson et al. 2017; Birkby et al.2017; Arcangeli et al. 2018). Indirect observations, like trans-mission spectroscopy, have demonstrated the presence of clouds(Sing et al. 2016; Nikolov et al. 2016; Pino et al. 2018; Gibsonet al. 2017; Kirk et al. 2018; Tregloan-Reed et al. 2018). Fareasier targets for atmosphere studies are brown dwarfs, whichare very similar to planets with respect to their physical param-eters and atmospheric processes. The coolest brown dwarfs (Ydwarfs) reach e ff ective temperatures as low as 250 K (Leggettet al. 2017; Luhman 2014). The observation of brown dwarfsallows us to identify the vertical cloud structures (Apai et al.2013; Buenzli et al. 2015; Yang et al. 2016; Helling & Casewell2014). To date, between 1500 and 2000 brown dwarfs are known(depending on whether late M-dwarfs and / or early L-dwarfs areincluded; Gagné et al. 2015; Best et al. 2018) and are rela-tively well-studied compared to the ∼ ff ects, but cloud parti-cles will also a ff ect the ionisation state of the atmosphere, whichis well known for solar system objects (Helling et al. 2016a,b).E ff orts are therefore ongoing to construct physical models de- scribing the formation of clouds in exoplanet and brown dwarfatmospheres. Such detailed models are necessary tools to pro-vide the context for observations and to uncover processes notdirectly accessible by observations. Part of this e ff ort is the con-sistent coupling of cloud formation with 1D atmosphere modelswith radiative transfer and convection (Tsuji et al. 1996; Acker-man & Marley 2001; Tsuji 2002; Witte et al. 2009; Allard et al.2012; Juncher et al. 2017; also Helling et al. 2008a), but alsoin 3D in order to study the time-dependent climate of extrasolarplanets (Lee et al. 2016; Lines et al. 2018a) and to understandobservational implications beyond 1D (Lee et al. 2017; Lineset al. 2018a).As our understanding of cloud formation progresses (e.g.Lee et al. 2015b; Krasnokutski et al. 2017; Hörst et al. 2019),including its implication for habitability (Narita et al. 2015), westart to refine our approaches. One long-standing discussion ishow to model the element replenishment in 1D cloud formingatmospheres, because without replenishment, a quasi-static at-mosphere must be cloud free (Appendix A in Woitke & Helling2004). Parmentier et al. (2013) utilised passive tracers to studythe atmospheric mixing in 3D (shallow water approximation)simulations for irradiated, dynamic but convectively stable atmo-spheres of (giant gas) planets. They observe that cloud particlesare distributed throughout the whole atmosphere.Parmentier et al. (2013) state: “ In statistical steady state, thisupward dynamical flux balances the downward transport due to
Article number, page 1 of 19 a r X i v : . [ a s t r o - ph . E P ] N ov & A proofs: manuscript no. Di ff uDrift particle settling and allows the atmospheric tracer abundance toequilibrate at finite (non-zero) values despite the e ff ect of par-ticle settling. The mechanism does not require convection, andindeed, the vertical motions that cause the upward transport inour models are resolved, large-scale motions in the stably strat-ified atmosphere. These vertical motions are a key aspect ofthe global-scale atmospheric circulation driven by the day-nightheating contrast. ” This assessment confirms our conclusion thatthe upward transport of condensable elements through the at-mosphere by mixing is indeed the key to understand cloud for-mation. However, challenges arise from the choice of the in-ner boundary condition (Carone et al. 2019), chemical gradi-ents (Tremblin et al. 2019), and the need to include cloud par-ticle feedback in order to test mixing parameterisations. A par-ticular interesting case will be the ultra-hot Jupiters where dayand night-sides can be expected to have very distinct (verti-cal) mixing patterns and scales. In this paper, we consider self-luminous giant gas planets, for which the irradiation from theirhost stars is negligible, such as young giant gas planets andbrown dwarfs. Brown dwarfs atmospheres are by now under-stood to be rather similar to giant gas planets, in particular at-mospheres from low-gravity brown dwarfs and young gas giants(Charnay et al. 2018).Moses et al. (2000) point out that large-scale mixing helpsto homogenise a gravitationally stratified atmospheres consist-ing of di ff erent kinds of molecules. This, however, only prevailsup to a certain altitude above which gas-kinetic di ff usion starts todominate over mixing (Zahnle et al. 2016). Di ff erent approacheshave been chosen to represent this vertical mixing in 1D atmo-sphere models (Ackerman & Marley 2001; Woitke & Helling2004; Helling et al. 2008a; Allard 2014; Juncher et al. 2017) andin 3D models (Lee et al. 2015a; Lines et al. 2018b) by measur-ing vertical velocity fluctuations and deriving mixing parame-terisations from 2D or 3D radiation-hydrodynamics simulations(Ludwig et al. 2002a; Freytag et al. 2010; Parmentier et al. 2013;Zhang & Showman 2018).Cloud formation modelling becomes an increasingly impor-tant part also of exoplanet / brown dwarf retrieval approaches forwhich, however, computational speed is an essential limitation.As part of the ARCiS retrieval platform, Ormel & Min (2019)presented a fast forward model that consistently solves di ff usivemixing and cloud particle growth for exoplanet atmospheres.In this paper, we present a new theoretical approach thatconsistently combines cloud formation modelling with di ff usivetransport for element replenishment. After presenting the mainformula body of our model in Sects. 2.1 to 2.3, we summariseour ansatz for handling the di ff usion coe ffi cient in Sect. 2.4, be-fore we present our main results in Sects. 4 and 5. We concludein Sect. 6. An overview of quantifying di ff usion coe ffi cients inthe literature is provided in Appendix D.
2. Cloud formation with diffusive transport of gasand cloud particles
Cloud formation involves at least seed particle formation (nucle-ation), surface growth and evaporation, element depletion, grav-itational settling and element replenishment. During their decentthrough the atmosphere, cloud particles may change phase or,more general, chemical composition, and may collide with eachothers leading to further growth. These cloud formation pro-cesses have been described previously (Woitke & Helling 2003;Helling & Woitke 2006; Helling & Fomins 2013) and di ff erent ARtful modelling Code for exoplanet Science cloud formation models have been compared by Helling et al.(2008a) with an update by Charnay et al. (2018). We thereforeonly provide a short summary here, a recent review can be foundin Helling (2019).Clouds are made of particles (aerosols, droplets, solid par-ticles). The formation of these particles requires condensationseeds, which are produced, for the case of the Earth atmosphere,by volcano eruptions, ocean sprays and wild fires. In absence ofthese crucial processes, which all require the existence of a solidplanet surface, cloud formation needs to start with the formationof seed particles through chemical reactions in the gas phase,involving the formation of molecular clusters. The formation ofseed particles requires a highly supersaturated gas. Once suchseed particles become available, many materials are already ther-mally stable and can condense on these surfaces simultaneously.Nucleation and growth reduce the local element abundances andhave a strong feedback on the local composition of the atmo-spheric gas. As macroscopic cloud particles form, they displaya spectrum of sizes as well as a mixture of condensed mate-rials. The local particle size distribution and the material mix-ture change as the cloud particles move through the atmosphere(hence, encounter di ff erent thermodynamic conditions), for ex-ample by gravitational settling (rain). Particle-particle collisionwill continue to shape the size distribution function. Cloud par-ticles may break up into smaller units (shattering) or stick to-gether to form even bigger units (coagulation). Cloud particlesmay also be transported upward and downward by macroscopicmixing processes. Particle-particle processes are not part of ourpresent model which focuses on the formation of cloud parti-cles and their feedback on the local chemistry through elementdepletion / enrichment. We note that the surface growth does shuto ff the nucleation process due to e ffi cient element depletion (Leeet al. 2015b) such that a simultaneous treatment of nucleationand growth is required in order to calculate the number of cloudparticles forming in the first place. As introduced in Woitke & Helling (2003), we consider the evo-lution of the size distribution function f ( V ) [cm − ] of cloud par-ticles in the particle volume interval V ... V + dV as a ff ected byadvection, settling, surface reactions and (new) by di ff usion ac-cording to the following master equation ∂ (cid:0) f ( V ) dV (cid:1) ∂ t + ∇ (cid:16) v ( V ) f ( V ) dV (cid:17) = (cid:88) k R k dV − ∇ φ d dV . (1) R k are the various gain and loss rates due to surface chemicalreactions, which lead to growth and evaporation of the particles(see Eqs. 59-62 for large Knudsen numbers, and Eqs. 68-71 forsmall Knudsen numbers in Woitke & Helling 2003). The vol-ume of the particles V is chosen as size variable to formulate thematerial deposit by surface reactions in the most straightforwardway. The last term in Eq. (1) accounts for the additional gainsand losses due to di ff usive mixing. The cloud particle velocity v ( V ) is assumed to be given by the hydrodynamical gas velocity v gas plus a vertical equilibrium drift velocity v ◦ dr ( V ) v ( V ) = v gas + v ◦ dr ( V ) . (2)Applying Fick’s first law (see e.g. Bringuier 2013), the di ff u-sive flux φ d of the cloud particles in volume interval V ... V + dV ( φ d dV has units [cm − s − ]) is given by the concentration gradi- Article number, page 2 of 19eter Woitke et al.: Dust in brown dwarfs and extra-solar planets ent of those particles φ d dV = − ρ D d ∇ (cid:32) f ( V ) dV ρ (cid:33) , (3)where D d [cm s − ] is the di ff usion coe ffi cient for those cloudparticles and ρ [g / cm ] the gas density. We introduce momentsof the cloud particle size distribution as ρ L j = (cid:90) ∞ V (cid:96) f ( V ) V j / dV . (4)Multiplying Eq. (1) with V j / and integrating over volume, weobtain the following system of moment equations for largeKnudsen numbers (see details in Woitke & Helling 2003) ∂ ( ρ L j ) ∂ t + ∇ (cid:0) v gas ρ L j (cid:1) = V j / (cid:96) J (cid:63) (cid:124) (cid:32) (cid:123)(cid:122) (cid:32) (cid:125) nucleation + j χ ρ L j − (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) growth − ∇ (cid:90) ∞ V l v dr ( V ) f ( V ) V j / dV (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) drift − ∇ (cid:90) ∞ V l φ d V j / dV (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) di ff usion , (5)where J (cid:63) [s − cm − ] is the nucleation rate and χ [cm / s] the netgrowth velocity. For large Knudsen numbers and subsonic veloc-ities (Epstein regime), the equilibrium drift velocity, also calledfinal fall speed, is given by Schaaf (1963) v ◦ dr = − √ π g ρ d a ρ c T (cid:98) r , (6)where a is the particle radius, ρ d the cloud particle material den-sity, (cid:98) r the unit vector pointing away from the centre of gravity,and g the gravitational acceleration. c T = (cid:112) kT / ¯ µ is an abbre-viation, T the temperature, k the Boltzmann constant, and ¯ µ themean molecular weight of the gas particles.Using Eq. (6) with a = (cid:16) V π (cid:17) / and assuming that the particledi ff usion coe ffi cient D d is independent of size, we can carry outthe integrations in Eq. (5). The final result is ∂ ( ρ L j ) ∂ t + ∇ ( v gas ρ L j ) = V j / (cid:96) J (cid:63) + j χ ρ L j − + ∇ (cid:32) ξ ρ d c T L j + (cid:98) r (cid:33) + ∇ (cid:16) D d ρ ∇ L j (cid:17) (7)with abbreviation ξ = √ π (cid:32) π (cid:33) / g . (8)A size-dependent di ff usion coe ffi cient, D d ( V ), would lead to anopen set of moment equations as discussed by Helling & Fomins(2013). We assume that all cloud particles are perfect spheres with well-mixed material composition which is independent of size, butdepends on time and location in the atmosphere (Helling &Woitke 2006; Helling et al. 2008c). Using the index s = ... S to distinguish between the di ff erent solid materials, for exampleAl O [s], TiO [s], Mg SiO [s] and Fe[s], we write V = (cid:88) s V s , L = (cid:88) s L s , J (cid:63) = (cid:88) s J s (cid:63) , χ = (cid:88) s χ s , b s = L s L , (9) Table 1.
Variable definitions and units symbol description unit z vertical coordinate cm n (cid:104) H (cid:105) hydrogen nuclei density cm − ρ = µ H n (cid:104) H (cid:105) gas mass density g cm − T gas temperature K a cloud particle radius cm V = π a volume of a cloud particle cm V s volume occupied by material s cm V (cid:96) minimum volume of a cloud particle cm f ( V ) size distribution function cm − ρ L j cloud particle moments cm j − L j j th moment cm j g − J (cid:63) nucleation rate cm − s − χ net growth speed cm / s ρ d mean cloud particle material density g cm − φ di ff usive flux cm − s − D d cloud particle di ff usion coe ffi cient cm s − D gas gas di ff usion coe ffi cient cm s − s index for di ff erent solid materials 1 ... Sr index for the surface reactions 1 ... R where b s is the volume fraction of material s in the cloud par-ticles . The mean cloud particle material density is given by ρ d = (cid:80) s b s ρ s where ρ s is the mass density of a pure material s . Most materials will not nucleate themselves ( J s (cid:63) = ∂ ( ρ L s ) ∂ t + ∇ ( v gas ρ L s ) = V (cid:96) J s (cid:63) + χ s ρ L + ∇ (cid:32) ξ ρ d c T b s L (cid:98) r (cid:33) + ∇ (cid:16) D d ρ ∇ L s (cid:17) . (10)Adding up Eqs. (10) for all solids s again yields Eq. (7) for j = ff erent materials grow at di ff erent speeds which dependon the amount of available atoms and molecules in the gas phaseand on the supersaturation ratio. Islands of some materials maygrow whereas others are thermally unstable and shrink. This be-haviour is obtained by summing up the contributions of all sur-face reactions r = ... R (for examples see Table 1 in Hellinget al. 2008c) χ = (cid:88) s χ s = (cid:88) s (cid:88) r c sr V s , (11)where V s [cm ] is the volume of one unit of material of kind s inthe solid state and c sr [cm − s − ] is the e ff ective surface reactionrate c sr = √ π ν sr n key r v rel r α r ν key r (cid:32) − S r (cid:33) × (cid:40) S r ≥ b s if S r < , (12)where n key r is the particle density [cm − ] of the key species of sur-face reaction r , α r the sticking probability, ν key r its stoichiometricfactor in that reaction, v rel r = (cid:112) kT / (2 π m key ) its thermal relative In Helling et al. (2008c), we have used the notation b s = V s / V tot where V s = ρ L s = (cid:82) ∞ V (cid:96) f ( V ) V s dV is the volume occupied by solid s pervolume of stellar atmosphere and V tot = ρ L = (cid:82) ∞ V (cid:96) f ( V ) V dV is the totalsolid volume per volume of stellar atmosphere.Article number, page 3 of 19 & A proofs: manuscript no. Di ff uDrift velocity and m key its mass. These growth rates are derived froma simple hit-and-stick model where we usually assume α r = α r (cid:44) S r is the reaction supersat-uration ratio as introduced in (Helling & Woitke 2006, see theirApp. B). For example, in the reaction2 CaH + + O ←→ [s] + (13)the key species is either CaH or Ti, depending on which speciesis less abundant. ν key r = s = CaTiO [s] is thesolid species, and ν sr = [s] are produced by onereaction.We note that Eq. (12) di ff ers slightly from our previ-ous definition (Eq. 4 in Helling et al. 2008c). The newgrowth / evaporation rates now always change sign at S r = b s . When supersatu-rated ( S r > s (see Fig. 1 in Helling & Woitke 2006). But forunder-saturation ( S r < b s . An integral part of our cloud formation model is the element con-servation. We must identify a replenishment mechanism whichis able to compensate for the losses of elements due to cloudparticle formation and settling in the upper atmosphere. In thispaper, we include di ff usion of gas particles along their concen-tration gradients. As cloud particles form, they consume certainelements in the upper atmosphere, creating a negative elementabundance gradient. Thus, gas particles containing those ele-ments will ascent di ff usively in that atmosphere. Analogous tothe formulation of the master equation for the dust particles, weformulate the element conservation as di ff usion-reaction system ∂ ( n (cid:104) H (cid:105) (cid:15) k ) ∂ t + ∇ ( v gas n (cid:104) H (cid:105) (cid:15) k ) = − (cid:88) s ν sk N (cid:96) J s (cid:63) − ρ L (cid:88) s (cid:88) r ν sk c sr + ∇ (cid:16) D gas n (cid:104) H (cid:105) ∇ (cid:15) k (cid:17) , (14)where (cid:15) k is the abundance of element k with respect to hydro-gen. The chemical reactions leading to nucleation and growthappear as negative source terms here, because they consume el-ements. We choose n (cid:104) H (cid:105) as density variable in Eq. (14), the to-tal hydrogen nuclei particle density, which is proportional to ρ in hydrogen-dominated atmospheres. n (cid:104) H (cid:105) (cid:15) k is the total numberdensity [cm − ] of nuclei of element k in any chemical form. ν sk is the stoichiometric factor of element k in solid s , for example ν TiO2[s]Ti =
1, and D gas [cm s − ] is the gas di ff usion coe ffi cient.For simplicity, we assume that all molecules are transported bythe same di ff usion coe ffi cient, which is valid within a factor 2or 3 for gas-kinetic di ff usion (sometimes called the binary di ff u-sion coe ffi cient, see Eq. (16) and App. D), and is entirely justifiedwhen eddy di ff usion dominates. The involved di ff usive gas ele-ment flux φ di ff k [cm − s − ] is given by φ di ff k = − D gas n (cid:104) H (cid:105) ∇ (cid:15) k . (15) The di ff usion coe ffi cients provide the kinetic information tocalculate the transport rates from concentration gradients (e.g. Lamb & Verlinde 2011). In general, gas and cloud particles dif-fuse with di ff erent e ffi ciencies because of their di ff erent inertiaand collisional cross sections with the surrounding gas.The determination of the gas di ff usion coe ffi cient D gas is ofcrucial importance for our model. We include gas-kinetic dif-fusion and large-scale turbulent (eddy) di ff usion as mixing pro-cesses. The gas-kinetic di ff usion coe ffi cient is given by D micro =
13 v th (cid:96) , (16)where (cid:96) = / ( σ n ) is the mean free path, n the total gas parti-cle density and σ ≈ . × − cm a typical cross-section forcollisions between the molecules under consideration with H .The thermal velocity is defined as v th = √ kT / ( π m red ) where m red is the reduced mass for collisions between the molecule andH (Woitke & Helling 2003). This gas-kinetic di ff usion ∝ / n isnegligible in the lower high-density layers of brown dwarf andplanetary atmospheres, where instead mixing by large-scale tur-bulent or convective motions is the dominant mixing process(Ackerman & Marley 2001; Ludwig et al. 2002b; Woitke &Helling 2004; Allard et al. 2012; Parmentier et al. 2013; Leeet al. 2015a). The large-scale (turbulent / convective / eddy) gasdi ff usion coe ffi cient is given by D mix ≈ (cid:104) v z (cid:105) L with L = α H p , (17)where (cid:104) v z (cid:105) is the root-mean-square average of the fluctuating partof the vertical velocity in the atmosphere, considering averagesover su ffi ciently large volumes and / or long integration times, L is the mixing length, and H p is the local pressure scale height. α is a dimensionless parameter of the order of one. We use α = α can be fine-tuned to describe theactual mixing scales as revealed by detailed hydrodynamic mod-elling. Inside the convective part of the atmosphere (cid:104) v z (cid:105) ≈ v conv is assumed, where v conv is the convective velocity, which is anintegral part of stellar atmosphere models, derived from mixinglength theory in 1D models. Above the convective atmosphere,where the Schwarzschild criterion for convection is false, (cid:104) v z (cid:105) decreases rapidly with increasing z , but never quite reaches zerodue to convective overshoot (see e.g. Brandenburg 2016). Weapply a power-law in log p to approximate this behaviourlog (cid:104) v z (cid:105) = log v conv − β (cid:48) · max (cid:110) , log p conv − log p ( z ) (cid:111) (18)with a free parameter β (cid:48) ≈ . ... . ff usion coe ffi cient is then D gas = D mix + D micro . (19)At high altitudes, the gas density n is small and hence D micro islarge, whereas D mix is small when β (cid:48) >
0. Therefore, at somepressure level in the atmosphere, the gas-kinetic di ff usion willstart to dominate. Figure 1 shows a typical structure as assumedin our models. The minimum of D gas around 10 − mbar corre-sponds to the crossover point (called the homopause ), upward ofwhich D micro dominates and the atmospheres is not well-mixed.Moses et al. (2000) draw similar conclusions concerning Sat-urn’s atmosphere. The maximum of D gas around p conv = . (cid:104) v z (cid:105) and D gas are approximately constant. Appendix D sum-marises some of the formulas currently applied in the literaturefor the gas di ff usion coe ffi cient. Article number, page 4 of 19eter Woitke et al.: Dust in brown dwarfs and extra-solar planets − − − − − − − − log p [bar] D [ c m / s ] − − − − − − − − / τ m i x [ / s ] D /τ mix Fig. 1.
The gas di ff usion coe ffi cient D gas (Eq. 19) in the new D iffu -D rift model for a brown dwarf atmosphere model with T e ff = g = / s ] and solar abundances. The grey line shows the inversemixing timescale τ mix as assumed in the previous D rift model. τ mix iscalculated according to Eq. (9) in (Woitke & Helling 2004). Both quan-tities are computed for β = β (cid:48) = y -axes show exactly 8 ordersof magnitude. The di ff usion of solid particles due to turbulent gas fluctua-tions was studied, in consideration of protoplanetary discs, byDubrulle et al. (1995), Schräpler & Henning (2004), Youdin &Lithwick (2007) and Riols & Lesur (2018). All works applymean field theory (also called Reynolds decomposition ansatz),where the densities and velocities of both the particles and thegas are decomposed into a mean component (that depends onlyon z ) and a small fluctuating part.The response of the solid particles to the turbulent gas varia-tions is then determined by comparing two timescales. The stop-ping or frictional coupling timescale is given by τ stop = a ρ d v th ρ with v th = (cid:115) kT π ¯ µ , (20)where a the particle radius. Equation (20) follows from a gen-eral relaxation ansatz τ stop = m (cid:0) ∂ F fric /∂ v dr (cid:12)(cid:12)(cid:12) v ◦ dr (cid:1) − , see Eq. (21) inWoitke & Helling (2003), for the special case of large Knud-sen numbers in a subsonic flow (the so-called Epstein regime),which we assume is valid here.The second timescale is the eddy turnover or turbulence cor-relation timescale τ eddy ( l ) in consideration of a spectrum of dif-ferent turbulent modes associated with di ff erent wave-numbers k or di ff erent spatial eddy sizes l . In a Kolmogorov type of powerspectrum P ( k ) ∝ k − / , any given cloud particle of size a tendsto co-move with all su ffi ciently large and slow turbulent eddieswhereas its inertia prevents following the short-term, small tur-bulence modes.In order to arrive at an e ff ective particle di ff usion coe ffi cient,the advective e ff ect of all individual turbulent eddies has to beaveraged, and thereby transformed into a collective di ff usive ef-fect. This procedure is carried out with di ff erent procedures andapproximations. The result of Schräpler & Henning (2004, seetheir Eq. 27), reads D d = D mix + St with St = τ stop τ eddy . (21) where D d is the size-dependent cloud particle di ff usion coe ffi -cient and St is the Stokes number of the particle in considerationof the largest eddy size L . The eddy turnover timescale of thelargest turbulence mode is given by τ eddy = L (cid:104) v z (cid:105) . (22)Both the size of the largest eddy L and the average of the fluctu-ating part of the vertical velocity (cid:104) v z (cid:105) are assumed to be identicalto the mixing length and velocity appearing in Eq. (17). Combin-ing the above equations we find St = a ρ d D mix v th ρ L . (23)The impact of the size dependence of D d on the cloud parti-cle moments was explored by Helling & Fomins (2013), whoshowed that this leads to an open set of moment equations, whichseems impractical for an actual solution. In the frame of thiswork, we will therefore only explore the two limiting cases ofvery large and very small Stokes numbers throughout the atmo-sphere. For small particles with St (cid:28)
1, we have D d → D mix andfor huge particles St → ∞ , we have D d → D d = D mix if all cloud particles are small,case 2: D d = D mix ( z ) is more flat or even increasingwith height, this might be di ff erent.
3. Static plane-parallel atmosphere
Before we proceed with the numerical solution of the full time-dependent model of cloud formation in di ff usive media, we firstdiscuss the 1D static case in order to better understand theexpected results from these equations. Considering the plane-parallel ( ∇ → d / dz ), static ( v gas =
0) and stationary case ( ∂/∂ t = = V j / (cid:96) J (cid:63) + j χ ρ L j − + ξ ddz (cid:32) ρ d c T L j + (cid:33) + ddz (cid:32) D d ρ dL j dz (cid:33) (25)0 = V (cid:96) J s (cid:63) + ρ L (cid:88) r c sr V s + ξ ddz (cid:32) ρ d c T b s L (cid:33) + ddz (cid:32) D d ρ dL s dz (cid:33) (26)0 = − (cid:88) s ν sk N (cid:96) J s (cid:63) − ρ L (cid:88) s (cid:88) r ν sk c sr + ddz (cid:32) D gas n (cid:104) H (cid:105) d (cid:15) k dz (cid:33) . (27) In the hydrostatic stationary case, the total vertical flux of ele-ments (due to vertical settling of cloud particles and di ff usivetransport) must be zero everywhere in the atmosphere and foreach element. This conclusion can be derived formally by addingtogether Eq. (27) and (cid:80) s (Eq. 26) · ν sk / V s , using V (cid:96) = N (cid:96) V s . Thechemical source terms (nucleation and growth terms) cancel out Article number, page 5 of 19 & A proofs: manuscript no. Di ff uDrift exactly, and in case D d = ddz (cid:32) D gas n (cid:104) H (cid:105) d (cid:15) k dz (cid:33) + ξ ddz ρ d c T L (cid:88) s ν sk b s V s = ⇒ D gas n (cid:104) H (cid:105) d (cid:15) k dz (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) − φ di ff k + ξ ρ d c T L (cid:88) s ν sk b s V s (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) φ settle k = const k (28) φ di ff k [cm − s ] is the upward element flux by di ff usion in thegas phase and φ settle k is the downward flux of elements containedin the settling cloud particles at this point. Equation (28) wouldstill allow for solutions with constant (i.e. time-independent andheight-independent) total element fluxes throughout the atmo-sphere, but this would require matching feeding and removingrates at the bottom and the base of the atmosphere, which doesnot seem to be physically plausible. Thus, const k = φ di ff k = φ settle k ⇒ d (cid:15) k dz = − ξ ρ d L c T D gas n (cid:104) H (cid:105) (cid:88) s ν sk b s V s ≤ . (29)According to Eq. (29), the element abundance gradients incloudy, static ( v gas =
0) and stationary ( ∂/∂ t =
0) atmospheresmust be negative because of the downward transport of elementsvia the precipitation of cloud particles, which must be balancedby an upward directed di ff usive flux of elements in the gas phase,which requires a negative concentration gradient. This conclu-sion is correct whenever cloud particles are present ( L >
0) andgravity is active ( ξ > ff erent elements areproportional to the element composition of the settling cloud par-ticles at this point. The abundances of all elements k involved incloud formation must decrease monotonically toward the top ofthe atmosphere.
4. Numerical solution of the time-dependent cloudformation problem
Equations (25) − (27) form a system of (3 + S + K ) coupled 2 nd order di ff erential equations, which can be transformed into asystem of 2 × (3 + S + K ) 1 st order ordinary di ff erential equa-tions (ODEs). Unfortunately, we have not been able to solve thisODE system directly. The boundary conditions are partly givenat the lower and partly at the upper boundary of the model, seeSect. 4.4. The integration direction must be downward in orderto model the nucleation of new cloud particles. Hence, we trieda shooting method where (cid:15) k ( z max ) is varied at the top of the at-mosphere until (cid:15) k ( z min ) is met, i.e. the given values in the deepatmosphere. We found it impossible to proceed this way. A tinychange of (cid:15) k ( z max ) in the 12 th digit was still observed to change (cid:15) k ( z min ) by a factor of two. The reason for this extreme sensitiv-ity seems to be the nucleation rate with its threshold character asfunction of supersaturation, and hence as function of (cid:15) k .Looking for alternatives, we found that a simulation ofthe time-dependent equations on a given vertical grid can beperformed by means of the operator splitting method as ex-plained in Sect. 4.2. We evolve the atmospheric cloud structure (cid:8) L j ( z , t ) , L s ( z , t ) , (cid:15) k ( z , t ) (cid:9) for a su ffi ciently long time, until it ap-proaches the time-independent case (cid:8) L ◦ j ( z ) , L s , ◦ ( z ) , (cid:15) ◦ k ( z ) (cid:9) , whichis the stationary structure we are interested in. Assuming a plane-parallel ( ∇ → d / dz ) and static ( v gas =
0) atmosphere, Eqs. (7), (10) and (14) read, including the time-dependent terms ddt (cid:16) ρ L j (cid:17) = V j / (cid:96) J (cid:63) + j χ ρ L j − (30) + ξ ddz (cid:32) ρ d c T L j + (cid:33) + ddz (cid:32) D d ρ dL j dz (cid:33) ( j = , , ddt ( ρ L s ) = V (cid:96) J s (cid:63) + χ s ρ L (31) + ξ ddz (cid:32) ρ d c T b s L (cid:33) + ddz (cid:32) D d ρ dL s dz (cid:33) ( s = ... S ) ddt ( n (cid:104) H (cid:105) (cid:15) k ) = − (cid:88) s ν sk N (cid:96) J s (cid:63) − ρ L (cid:88) s (cid:88) r ν sk c sr (32) + ddz (cid:32) D gas n (cid:104) H (cid:105) d (cid:15) k dz (cid:33) ( k = ... K ) . The moment Eqs. (30) and (31) are not closed because L ap-pears twice of the right side, a consequence of larger particlessettling faster (Eq. 6). Therefore, a numerical solution requires aclosure condition as L = F ( L , L , L , L ) . (33)We use the closure condition explained in the appendix A.1 of(Helling et al. 2008c). The idea is to approximate the particlesize distribution f by a double δ -function which has four param-eters. These parameters are determined by matching the givenmoments L , L , L and L , and result in the forth moment L according to the definition of the dust moments (Eq. 4). Figure 2 visualises our numerical approach using the operatorsplitting method (Klein 1995).1. We update L j and L s only according to the settling sourceterms (the terms on the right side of Eqs. (30) and (31) con-taining L j + and L ), applying half a timestep ∆ t /
2, see de-tails in App. B.2. We call the di ff usion solver for half a timestep ∆ t / (cid:15) k and, optionally, L j and L s , if the cloud particles are to bedi ff used as well, see App. A.3. We integrate the chemical source terms (nucleation, growthand evaporation) for a full timestep ∆ t . These equationsare sti ff at high densities and require an implicit integrationscheme. We use the implicit ODE-solver L imex T , density n (cid:104) H (cid:105) and element abundances (cid:15) k , wecall the equilibrium chemistry code GG chem (Woitke et al.2018) to calculate all molecular concentrations; (ii) those re-sults are used to calculate the reaction supersaturation ratios S r ; (iii) the nucleation rates J s (cid:63) and net surface reaction rates c sr are determined.4. We finish the timestep by calling again the di ff usion solverfor ∆ t / ∆ t / Article number, page 6 of 19eter Woitke et al.: Dust in brown dwarfs and extra-solar planets Sett Di ff Di ff Sett OP
Fig. 2.
Operator splitting calling sequence. Sett = gravitational settling,Di ff = di ff usion, CF = cloud formation (nucleation, growth and evapora-tion), and OP = output. 1 / is why we do not split CF (Fig. 2) but put it with a full timestepin the centre of the operator splitting calling sequence. The cloudformation part of the code is parallelised and can be executed forall atmospheric layers independent of each other. In order to produce accurate 2 nd order solutions, the timestepmust be limited to make sure that each operator remains inthe linear regime. For example, the sole application of a cloud-chemistry timestep must not change the amount of dust or the el-ement abundances substantially in any computational cell. In or-der to achieve code stability and accuracy, we limit the timestepas follows:1. The cloud particles must not jump over layers by settling ∆ t < . ∆ z v dr , j (34)where ∆ z is the vertical grid resolution and v dr , j is the meandrift velocity a ff ecting moment ρ L j as given by Eq. (B.3).This is the usual Courant-Friedrichs-Lewy (CFL) criterionto stabilise explicit advection scheme, with an additionalsafety-factor 1 / ∆ t , must not change any of the gas elementabundances by more than a given maximum relative change(default accuracy is 15%).3. The timestep must not exceed the maximum explicit di ff u-sion timestep (Eq. A.27).If one of these criteria becomes false during the simulation, thetimestep is discarded and ∆ t reduced. If, on the contrary, the cri-teria are met easily, ∆ t is increased for the subsequent timestep. As our upper boundary condition, we assume that there areno cloud particles settling into the model volume from abovev dr , j ( z max ) =
0. In the di ff usion solver, we use a zero-flux (closedbox) upper boundary condition, i.e. the gradients of (cid:15) k are as-sumed to be zero at z = z max . The same applies to the cloud parti-cle moments ddz L j ( z max ) = ff used as well.The lower boundary is placed well below the main silicatecloud layer to make sure that the temperature is too high toallow for any cloud particles to exist near the lower boundary L j ( z min ) =
0. We also demand that the element abundances atthe lower boundary equal the given values as present deep in theatmosphere (cid:15) k ( z min ) = (cid:15) k , where the (cid:15) k are considered as freeparameters, for example solar abundances (Asplund et al. 2009). We start all simulations from a cloud-free atmosphere L j ( z , t = =
0. Concerning the element abundances, we have experi-mented with two cases: (1) an ‘empty’ atmosphere (cid:15) k ( z , t = = (cid:15) k ( z , t = = (cid:15) k , where the index k is applied to all elements which can potentially be transformedcompletely into solids ( k = Si, Mg, Fe, Al, Ti, ...), but not H, He,C, N, O, etc. For the latter elements we put (cid:15) k ( z , t = = (cid:15) k inboth cases. We found an identical final structure in both cases(see App. C), which is very reassuring. The models calculatedfrom initial condition (2), however, need much more computa-tional time to complete. In this case, the nucleation rate is ini-tially huge and a very large number of tiny cloud particles arecreated shortly after initialisation, which take a long time to set-tle down in the atmosphere.In order to reach the final relaxed, time-independent state, themodel must be evolved until (i) the atmosphere is completely re-plenished several times by fresh elements ascending di ff usivelyfrom the lower boundary to the very top and (ii) new grainsformed high in the atmosphere have su ffi cient time to settle downto the cloud base and evaporate. In comparison, the chemicalprocesses are typically quite fast. We need to evolve one modelfor about 10 timesteps, which, depending on global parameterslike T e ff and log g , translates into real evolutionary times be-tween a few months to a few tens of years. On a parallel cluster,one can complete one such model in a few days real time whenusing 16 processors (about 500 CPU hours), where the chemicalequilibrium solver GG chem is called a few 10 times.
5. Results
In the previous Helling & Woitke cloud formation models(Woitke & Helling 2004; Helling & Woitke 2006; Helling et al.2008c), henceforth called the D rift models, the replenishment ofelements was treated in a di ff erent way, using a prescribed massexchange timescale τ mix ( z ) to replenish the atmosphere withfresh elements from the deep as n (cid:104) H (cid:105) ( (cid:15) k − (cid:15) k ) /τ mix . The mass ex-change timescale was approximated by a powerlaw log τ mix ( z ) = const − β log p ( z ) with power-law index β = . rift modelwith the new di ff usive model, henceforth called the D iffu D rift model. Both approaches model seed formation, kinetic surfacegrowth / evaporation of cloud particles and gravitational settlingin the same way, but di ff er in the treatment of the mixingwhich enters the cloud formation and the element conservationequations. The underlying temperature / pressure structures forall models discussed in this paper are taken from a the D rift -P hoenix atmosphere grid (Dehn et al. 2007; Helling et al. 2008b;Witte et al. 2009, 2011). In Fig. 3 we have selected a model withe ff ective temperature T e ff = g = Z = rift -P hoenix models solve the complete1D model atmosphere problem including convection, radiativetransfer and hydrostatic structure, coupled to our previous D rift model, where the cloud opacities are calculated by Mie and ef-fective medium theory. The resulting atmospheric structure arefrozen for this study, i.e. the feedback of the new cloud formationmodel on the ( p , T )-structure is not included.The chemical setup for this comparison has 16 elements(H, He, Li, C, N, O, Na, Mg, Si, Fe, Al, Ti, S, Cl, K,Ca), one nucleation species (TiO ), 12 solid species (TiO [s],Al O [s], CaTiO [s], Mg SiO [s], MgSiO [s], SiO[s], SiO [s],Fe[s], FeO[s], MgO[s], FeS[s], Fe O [s]) and 60 surface reac- Article number, page 7 of 19 & A proofs: manuscript no. Di ff uDrift log p [bar] T [ K ] DRIFTDIFFUSION (gas+particles)DIFFUSION (gas) log p [bar] -20-15-10-5 l og J [ / s / c m ] log p [bar] -10-50 l og n d [ c m − ] log p [bar] -4-3-2-1012 l og › a fi / [ µ m ] log p [bar] l og › v d r fi [ c m / s ] -8 -7 -6 -5 -4 -3 -2 -1 log p [bar] -12-10-8-6-4-2 l og du s t / ga s Fig. 3.
Comparison of cloud formation models for T e ff = g =
3, metallicity Z =
1, and β = β (cid:48) =
1. The previous D rift model isshown by the thick grey lines. Two D iffu D rift models are overplottedassuming pure gas di ff usion (dashed lines) and gas + particle di ff usion(black solid lines). dust / gas = ρ d L is the dust-to-gas mass ratio, n d = ρ L the number density of cloud particles, (cid:104) a (cid:105) / = (cid:0) L / (4 π L ) (cid:1) / themass-mean particle radius, and (cid:104) v dr (cid:105) = (cid:104) a (cid:105) / √ π g ρ d / (2 ρ c T ) the cor-responding drift velocity according to Eq. (6). tions. The molecular setup in the new models is not quite identi-cal, since the D rift model uses a previous version of the chemi-cal equilibrium solver GG chem , which has been replaced by thelatest version (Woitke et al. 2018) in the D iffu D rift model. Weuse 189 molecules in D rift and 308 molecules in D iffu D rift to find the molecular concentrations in chemical equilibrium.We do not see any substantial di ff erences in molecular concen-trations caused by this data update, unless the local tempera-ture falls below about 400 K. Also the thermochemical data forthe selected solids is not entirely identical, but these di ff erencesare not substantial either, because the local temperatures remainabove 700 K in this test. We assume the mixing powerlaw index to be β = β (cid:48) = τ mix ( z ) in D rift and D gas ( z ) in D iffu -D rift , see Eq. (18) and Fig. 1, albeit the meaning of β and β (cid:48) isslightly di ff erent. We note that using β (cid:48) > β would likely produceresults that are more similar to each other than those presented inthis paper. The lower volume boundary for the size integrationof the cloud particle moments is set to V (cid:96) = × V TiO where V TiO = . × − cm is the assumed volume of one unit ofsolid TiO [s].The resulting cloud structures, as predicted by our previ-ous D rift and the new D iffu D rift models, are compared inFig. 3. The di ff usive transport of condensable elements up intothe high atmosphere with D iffu D rift is much less e ffi cient thancompared to the assumed replenishment in the D rift model. Asthese elements are slowly mixed upwards by di ff usion, they cancollide with existing cloud particles to condense on, and hencemuch less of these elements reach the high atmosphere where thenucleation takes place. This is the main di ff erence between theD rift and the D iffu D rift models. In the previous D rift models,the mixing process was assumed to take place instantly. Cloud structure:
Consequently, the new D iffu D rift model isfeatured by up to 5 orders of magnitude lower nucleation rates(Fig.3) and less cloud particles high in the atmosphere. At inter-mediate pressures (10 − ... − bar) we find that the fewer cloudparticles in the D iffu D rift model grow quickly and reach parti-cle sizes of about 10 µ m at 1 mbar, wheres in the D rift model,since there are so many of them, the cloud particles remainsmaller, about 0.3 µ m. The growth of the cloud particles is lim-ited by the amount of condensable elements available per parti-cle, and therefore, this e ff ect is expected.The dust-to-gas mass ratio, ρ d /ρ gas , increases more steeplyin the D iffu D rift model, but reaches about the same maximumof order 10 − at 1 mbar as in the D rift model. Thus, overall, thecloud formation process is about equally e ff ective, but the cloudsare spatially more confined in the D iffu D rift model, reaching upjust a few scale heights above the cloud base.Table 2 lists vertically integrated cloud column densities forthe three models. We find values of a few milli-grams of conden-sates per cm , where the D iffu D rift model without cloud parti-cle di ff usion is found to be the most dusty one. Using an order ofmagnitude estimate of cloud particle opacities (see Appendix E),values between several 100 cm / g to several 1000 cm / g at λ = µ m are expected, depending on material and particle size dis-tribution, i.e. a column density of 1 milli-gram of condensate percm roughly corresponds to an optical depth of one at λ = µ m. Table 2.
Comparison of cloud column densities [mg / cm ] (cid:63) for the threemodels shown in Fig. 3 and discussed in the text. condensate D rift D iffu D rift D iffu D rift dust + gas di ff usion gas di ff usionTiO . × − . × − . × − Al O .
57 0 .
47 7 . .
29 0 .
040 0 . SiO .
54 0 .
59 0 . .
48 0 .
034 0 . .
14 6 . × − . × − Fe 1 . .
70 0 . . × − . × − . × − MgO 0 . . × − . × − FeS 2 . × − . × − . × − CaTiO .
040 0 .
022 0 . O . × − . × − . × − total 4.0 1.9 10.1 (cid:63) : column densities are calculated as Σ s = (cid:82) ρ s ρ L s dz .Article number, page 8 of 19eter Woitke et al.: Dust in brown dwarfs and extra-solar planets This implies that all three models discussed here have opticallythick cloud layers.The computation of more realistic cloud particle opacitieswill need to take into account the height-dependent materialcomposition, size and possibly shape distribution of the cloudparticles, as done, for example, by Dehn et al. (2007), Witte et al.(2009, 2011) and Helling et al. (2019). In comparison to the pre-vious D rift models, the particles in the D iffu D rift models arelarger, which is likely to cause the optical depths to be somewhatsmaller, although the cloud mass column densities are similar. Inaddition, molecular opacities need to be added to calculate thespectral appearance of the objects and feedback onto the ( p , T )-structure, which goes beyond the scope of the present paper.The resulting particle properties in the main cloud layer be-low 1 mbar depend not only on the treatment of mixing, but alsowhether or not we switch on the dust di ff usion in the D iffu D rift model. In this region, the cloud particles stepwise purify chem-ically as they decent in the atmosphere (see Fig. 5). Thermody-namically less stable materials like Fe O [s], FeO[s] and FeS[s]sublimate o ff the cloud particles sooner. Subsequently, the abun-dant magnesium-silicates MgSiO [s] and Mg SiO4[s] sublimateas well, which causes the cloud particles to shrink significantlyaround 1 mbar. Close to the cloud base, at about 1800 K in thismodel, only the most refractory materials remain, in particularmetal oxides such as Al O [s], CaTiO [s] and TiO [s], beforeeven these materials eventually sublimate and the cloud particlesevaporate completely.As one material sublimates, the liberated elements may re-condense into di ff erent condensates, which are thermodynami-cally more stable, leading to rapid changes in particle size andmaterial composition. There is also a dynamical e ff ect. When theparticles shrink, their fall speeds decrease which leads to spatialaccumulation, hence the number density of cloud particles n d increases. While these e ff ects and the general behaviour of thecloud particles are similar in all three models, the steps of sub-limation are more pronounced in the D iffu D rift model withoutdust di ff usion. Dust di ff usion tends to smooth out the variationsof particle size and density. Element abundances:
Figure 4 compares the resulting gas ele-ment abundances. We see a strong depletion of condensable el-ements in the main cloud layer in all three models, by up to 5orders of magnitude, concerning elements Ti, Al, Mg, Si, Mgand Fe. However, the details are di ff erent. The previous D rift model is featured by minimums of (cid:15) k that are similar in depth ascompared to the overall decrease of (cid:15) k in the D iffu D rift mod-els. High up in the atmosphere, where cloud particles are virtu-ally absent, there is no surface to condense on, and so the in-stantaneous mixing assumption in the D rift models causes a re-increase of (cid:15) k toward the top of the atmosphere, unless the ele-ment can form nuclei. In the extremely low density gas at theseheights, these nuclei simply fall through the atmosphere with-out much interaction, whereas elements, which cannot nucleate,accumulate.In contrast, in the new D iffu D rift models, the abundances ofall elements involved in cloud formation decrease with height ina monotonic way. This behaviour is expected in the final, time-independent, relaxed state of the atmosphere as discussed inSect. 3.1. In the stationary case, the downward transport of con-densable elements via the falling cloud particles must be com-pensated by an upward di ff usive transport of these elements inthe surrounding gas, which implies negative element abundancegradients throughout the cloudy atmosphere, see Eq. (29). The log p [bar] l og † ga s CNOSSiMgFeAlTi log p [bar] l og † ga s CNOSSiMgFeAlTi
Fig. 4.
The impact of our assumptions about the mixing processes in theatmosphere on the resulting gas element abundances, in models with thesame parameters as in Fig. 3. The upper plot shows (cid:15) gas k for instanta-neous mixing as assumed in the previous D rift model. The lower plot shows the results according to the new D iffu D rift models, where thefull lines show the results for gas and dust di ff usion, and the dashedlines show the results for gas di ff usion only. total drop of element abundances is deepest for Ti, but less deepfor Si and Fe as compared to Mg.Freytag et al. (2010) performed two-dimensional radiationhydrodynamical simulations of substellar atmospheres which in-cluded a time-dependent description for the formation of a singlekind of cloud particles for a fixed concentration of seed particles.The paper discusses substellar objects with T e ff =
900 K − g ) = iffu D rift results, where abundances of condensable elements in the gasphase are decreasing fast in the cloud layers, and stay aboutconstant above the clouds. We note that Freytag et al. have pre-scribed the number of seed particles and considered only onegeneric condensate in their simulations. Cloud particle composition:
Figure 5 shows the correspondingsolid material composition (by volume) of the cloud particles.In all three models we assume that the cloud particles have awell-mixed material composition which is independent of size,but that composition changes as the particles fall through theatmosphere, hence material composition depends on height. Allthree models show the same basic vertical structure.
Article number, page 9 of 19 & A proofs: manuscript no. Di ff uDrift log p [bar] v o l u m e f r a c t i o n l og b s DRIFT
TiO2Al2O3CaTiO3Fe FeOFe2O3FeSMgO Mg2SiO4MgSiO3SiOSiO2 log p [bar] v o l u m e f r a c t i o n l og b s DIFFUSION (gas and particles) log p [bar] v o l u m e f r a c t i o n l og b s DIFFUSION (gas only)
Fig. 5.
The material volume composition of the cloud particles b s = L s3 / L = V s / V tot for the same three models as discussed in Figs. 3 and 4.
1. A layer containing only the most stable metal-oxides at thecloud base, in this model Al O [s], TiO [s] and CaTiO [s].The position of the cloud base, which depends on T e ff andlog g , is located at around 1800 K in this model.2. A thin layer of cloud particles around 1500 K which mainlyconsist of metallic Fe[s].3. Main silicate cloud layer composed of Mg SiO [s],MgSiO [s], MgO[s], SiO[s] and SiO [s], mixed with metal-lic iron, upward of about 1400 K in this model.4. Less stable solid materials such as FeS[s], FeO[s] andFe O [s] are incorporated into the silicate cloud particles attemperatures lower than about 1100 K, 1000 K and 850 K,respectively, in this model.5. Pure nuclei at the top, here TiO [s], which fall through theatmosphere so quickly that they practically do not grow.Further inspection shows, however, that the material composi-tion of the main silicate cloud layer (3) di ff ers substantially be-tween the D rift and the D iffu D rift models. In the new di ff usiveD iffu D rift models, the first magnesium-silicate to form is fos-terite Mg SiO [s], which has a stoichiometric ratio Mg : Si = SiO [s] stops once the reservoirof Mg is exhausted, still leaving about half of the available Siin the gas phase. Since the mixing is di ff usive, very little Mg T eff [K] l og Σ [ g / c m ] TiO2Al2O3MgSiO3 Mg2SiO4SiOSiO2 FeFeSFeO MgONa2S
Fig. 6.
Column densities [g / cm ] of di ff erent condensates in the atmo-sphere along a sequence of models with decreasing T e ff but constantlog g = β (cid:48) =
1. A value of 10 − g / cm roughly corresponds to anoptical depth of one at a wavelength of λ = µ m (see Appendix E). can be mixed upwards through these Mg SiO [s] clouds. Thus,the remaining amount of Si above the Mg SiO [s] clouds prefer-entially forms other silicate materials, in particular SiO [s] andSiO[s], but only very little MgSiO [s]. This is di ff erent in ourprevious D rift model where the depleted elements are assumedto be instantly replenished at similar rates, in which case bothMg SiO [s] and MgSiO [s] are found to be about equally abun-dant condensates in the main silicate cloud layer.Another di ff erence concerns FeS[s] (troilite). FeS[s] is foundto form in large quantities in the previous D rift models, causing (cid:15) S to drop significantly, see upper part of Fig. 4. However, thisdepends on our assumptions about how the elements are replen-ished. In the new di ff usive models, upward mixing of gaseous Feis rather ine ffi cient because the Fe atoms have plenty of oppor-tunity to condense in form of Fe[s] on existing cloud particlesalong their way upwards in the atmosphere. Once the tempera-ture is low enough to allow FeS[s] to form, there is so little Fe leftin the gas phase that the S abundance is more or less una ff ectedby FeS[s] formation, and therefore sulphur remains available forother condensates to form. T e ff In this section, we study the results of a sequence of the new D if - fu D rift cloud formation models with decreasing e ff ective tem-perature T e ff . We are using a slightly di ff erent chemical setuphere that will allow us to discuss secondary cloud layers. Weconsider four nucleation species (TiO ) N , (SiO) N , (KCl) N and(C) N . The nucleation rates of TiO , KCl and C are calculated bymodified classical nucleation theory (Helling et al. 2017; Gailet al. 1984), with a surface tension value for KCl from (Lee et al.2018). The nucleation rate of SiO is calculated according to (Gailet al. 2013). We have 16 elements in this setup (H, He, Li, C, N,O, Na, Mg, Si, Fe, Al, Ti, S, Cl, K, Ca), 14 condensed species(TiO [s], Al O [s], MgSiO [s], Mg SiO [s] ,SiO[s], SiO [s],Fe[s], FeS[s], FeO[s], MgO[s], KCl[s], NaCl[s], Na S[s], C[s]),
Article number, page 10 of 19eter Woitke et al.: Dust in brown dwarfs and extra-solar planets log p [bar] T [ K ] Teff=2800KTeff=2600KTeff=2400KTeff=2100KTeff=1900KTeff=1700KTeff=1500KTeff=1300K log p [bar] l og D ga s [ c m / s ] Teff=2800KTeff=2600KTeff=2400KTeff=2100KTeff=1900KTeff=1700KTeff=1500KTeff=1300K log p [bar] l og du s t / ga s Teff=2800KTeff=2600KTeff=2400KTeff=2100KTeff=1900KTeff=1700KTeff=1500KTeff=1300K log p [bar] l og † S i log p [bar] l og J [ c m − s − ] J total J TiO J SiO log p [bar] l og › a fi / [ µ m ] Teff=2800KTeff=2600KTeff=2400KTeff=2100KTeff=1900KTeff=1700KTeff=1500KTeff=1300K
Fig. 7.
Sequence of cloud forming models with decreasing T e ff at constant log g = β (cid:48) = Top row: gas temperature T and di ff usion coe ffi cient D gas as function of pressure (both assumed). Middle low: resulting dust to gas mass ratio and element abundance ofsilicon in the gas phase (cid:15) Si . Lower row: resulting nucleation rates J (cid:63) and mean particle sizes (cid:104) a (cid:105) / .
308 molecules and 50 surface growth reactions. Molecular equi-librium constants and Gibbs free energies of the condensates areall taken from Woitke et al. (2018). Dust di ff usion is included inall models.Figure 6 shows the total column densities of selected cloudmaterials Σ s [g / cm ] in a series of D iffu D rift models with con-stant log g = β (cid:48) =
1, but decreasing T e ff . Thecolumn densities of the condensed species are computed as Σ s = (cid:90) ρ s ρ L s dz , (35)where ρ s [g / cm ] is the material density of the pure condensateof kind s and ρ L s [cm / cm ] is the volume of condensed kind s per volume of atmosphere. For example, for T e ff = SiO [s], Fe[s] and Al O [s], followed by SiO[s]and MgSiO [s].On the left side of this plot, the first model that shows con-densation appears at T e ff = O [s] and TiO [s]. In the nextfew models down to T e ff = ff ective temperatures,Al O [s] still has the largest column density because the metaloxide layer is situated deeper in the atmosphere where the den-sities are higher. Only for T e ff < T e ff < S[s].Figure 7 shows a few more details from this T e ff -series ofnew cloud formation models. The upper left plot shows the at- Article number, page 11 of 19 & A proofs: manuscript no. Di ff uDrift mospheric density / temperature structures assumed (taken fromthe D rift -P hoenix atmosphere grid (Dehn et al. 2007; Hellinget al. 2008b; Witte et al. 2009, 2011). The kinks in deep layers( T ∼ − ff usion coe ffi cient in theatmosphere, which decreases with T e ff , because the convectivelayer sinks into deeper layers, hence the spatial distance to thesource causing the mixing motions in the atmosphere increases.The left middle plot in Fig. 7 shows the dust-to-gas mass ra-tio, which has its maximum in the main silicate-iron layer, and ashoulder on the right due to the deeper metal-oxide clouds whichare made of the rarer elements with the highest condensationtemperatures, namely aluminium, calcium and titanium. As T e ff decreases, both layers move inward to deeper layers and becomesuccessively more narrow, until finally, for T e ff = S[s]. The right middle plot shows how the silicon abundancein the gas phase is a ff ected. All curves are monotonic decreas-ing towards the surface, with higher Si depletions for lower T e ff where the silicate cloud particle formation is more complete.The nucleation rates of (TiO ) N and (SiO) N particles are de-picted in the lower left plot. A complicated, double-peaked pat-tern shows, which has a minimum around the main peak of thedust-to-gas ratio (at the peak position of the main silicate-ironlayer). (TiO ) N is usually the most significant nucleation species,but cooler models show additional contributions by (SiO) N . Theresulting mean particle sizes are plotted on the lower right, witha tendency to produce larger particles deep in the atmosphere forlower T e ff . An in-between minimum in particle size occurs wherethe main silicate material evaporates. Only the coolest model hasa second minimum where Na S[s] evaporates. Interestingly, thehottest and the coolest model in Fig. 7 show about equally largecloud particles at high altitudes, whereas all other models showsmaller particles.
6. Summary and Discussion
This paper has introduced a new cloud formation model appli-cable to the atmospheres of brown dwarfs and gas giant (exo-)planets. We have combined our previous cloud particle mo-ment method (Woitke & Helling 2004; Helling & Woitke 2006;Helling et al. 2008c) with a di ff usive mixing approach, accordingto which, in the final relaxed time-independent state of the atmo-sphere, fresh condensable elements are di ff usively transportedupwards to replenish the upper atmosphere via a combination ofturbulent (eddy) mixing and gas-kinetic di ff usion. Our formula-tion of the problem arrives at a system of about 30 second orderpartial di ff erential equations of reaction-di ff usion type, wherethe formation and growth of the cloud particles follows from akinetic treatment in phase-non-equilibrium. Model setup:
The new cloud formation model is applied to agiven one-dimensional ( p , T ) atmospheric structure in this pa-per. The model is calculated time-dependently, using an oper-ator splitting technique. All models are found to relax towarda time-independent, stationary solution, where the condensableelements are constantly mixed up di ff usively, cloud particles nu-cleate from the gas phase high in the atmosphere, grow by thesimultaneous condensation of di ff erent solid materials on theirsurface, and then decent through the atmosphere due to gravita-tional settling, before the particles stepwise purify and eventu-ally sublimate completely at the cloud base. Fig. 8.
Concentration of condensed species in a model with T e ff = g = β (cid:48) ==
1, showing a secondary cloud layer almostentirely made of di-sodium sulfide Na S[s]. n s cond = ρ L s / V s mat [cm − ] isthe number density of solid units of condensed species s , V s mat = m s /ρ s isthe volume occupied by one unit of solid s in the pure material, and m s is the mass of one such units, for example 100.4 amu for s = MgSiO . n s cond / n (cid:104) H (cid:105) is directly comparable to element abundances. Timescales:
The real-time simulation time required to reachthat stationary solution varies between a few months to severaltens of years, depending on log( g ) and T e ff . The relaxation isquicker when models are started from an atmosphere that is de-void of any condensable elements at t =
0. These relatively longsimulation times make these models computationally expensive(of order 500 CPU-hours per model), because the intrinsic nu-cleation and growth reactions are very fast, which means thatthe models need to be advanced on short computational timesteps of the order of seconds to guarantee numerical stability.The long physical timescales involved in the simulations are (i)the overall settling time for small particles inserted high in theatmosphere, and (ii) the overall mixing time for gas parcels todi ff usively reach the highest point in the model from the cloudbase. This implies that 3D simulations of cloud formation (GCMmodels, for example Freytag et al. 2010; Lee et al. 2016; Lineset al. 2018a; Powell et al. 2018; Charnay et al. 2018) must be ad-vanced for similar real-time simulation times before a relaxedphysical structure can be expected. However, how long thesephysical timescales actually are will depend on the exact for-mulation of mixing and setting in the GCM models. Cloud density and particle sizes:
In comparison to our previ-ous D rift models, the D iffu D rift models show fewer but largercloud particles, which are more concentrated towards the cloudbase. However, the physical properties of the cloud particles inthe main silicate-iron layer towards the bottom of the clouds(dust to gas mass ratio, particle sizes, optical depth, chemicalcomposition, etc.) are found to be similar to the results of theprevious models. The dust-to-gas ratio in the main silicate-ironlayer reaches a peak value of about 0.002 to 0.003, quite inde-pendent of T e ff , for not too hot models ( T e ff > iffu D rift modelsis that the di ff usive element replenishment is less e ff ective for theupper atmosphere, because the molecules carrying the elementsdi ff usively upwards have a high probability to collide with exist- Article number, page 12 of 19eter Woitke et al.: Dust in brown dwarfs and extra-solar planets ing cloud particles on their way up the atmosphere. This e ff ectwas not accounted for in the previous models. Element abundances:
The concentration of condensable ele-ments in the gas phase shows a steep decline in the D iffu D rift models above the cloud base, followed by a monotonous de-crease towards a plateau which then continues on that level to-ward the upper boundary of the model. This behaviour is ex-pected in the time-independent relaxed state, because the down-ward flux of condensable elements due the falling cloud parti-cles must be compensated for by an upward di ff usive flux of ele-ments in the gas phase, which requires a negative concentrationgradient. We find the abundances of the condensable elementshigh above the cloud layers to be strongly dependent on e ff ec-tive temperature, in agreement with the results of 2D radiationhydro-models by Freytag et al. (2010). For example, the sili-con abundance is reduced by about 2.5 orders of magnitude for T e ff = T e ff = Cloud material composition:
The chemical composition of thecloud particles is characterised by (i) a deep layer with themost stable metal-oxides at the cloud base (Al O [s], TiO [s],CaTiO [s] in the current setup), (ii) the main silicate-iron layer(mainly Mg SiO [s] and Fe[s]) which, with increasing height,is then mixed with other silicates like SiO[s] and SiO [s] andMgSiO [s]. Some less stable condensates are also found insmaller quantities, in particular FeS[s] and FeO[s]. The conden-sation in both these cloud layers leads to a removal of certainelements from the gas phase, and the stoichiometry of the con-densates decides upon which elements remain for further con-densation higher in the atmosphere. In particular, the formationof Mg SiO [s], with stoichiometry Mg : Si = [s] above the Mg SiO [s] layer,rather than the formation of MgSiO [s], which is a rather unim-portant condensate in our new D iffu D rift models. Having solittle MgSiO [s] in the main silicate-iron layer is a result thatdi ff ers from the results obtained with our previous D rift mod-els, and from phase-equilibrium models starting from completesolar abundances (Woitke et al. 2018). Na S[s] clouds:
Our coolest D iffu D rift models show the oc-currence of a secondary cloud layer almost entirely made of di-sodium sulfide Na S[s], see Fig. 8. The presence of Na S-cloudsin brown dwarf atmospheres has been proposed by Morley et al.(2012) to fit the optical appearance of two red T-dwarfs. Theformation of Na S-clouds requires the presence of sulphur andsodium in the gas phase at low temperatures. In phase equilib-rium models starting from solar abundances, such a combinationis prevented by the formation of FeS[s] (troilite), which con-sumes the sulphur. However, in our new di ff usive kinetic cloudformation models, iron is depleted by the formation of metalliciron Fe[s] at high temperatures, so FeS[s] cannot form in largequantities. Consequently, sulphur remains available to eventu-ally form Na S[s] at lower temperatures. The condensation ofNa S[s] then reduces the possibility to form NaCl[s] at evenlower temperatures, and so on. Therefore, the new di ff usive D if - fu D rift models reveal new details about the condensation se-quence in cloudy atmospheres, and we need more experimentswith our selection of condensates during model initialisation toarrive at more distinct conclusions.
7. Conclusions
The physical description of the replenishment mechanism forcondensable elements in planetary atmospheres seems crucialfor realistic cloud formation models. This paper has used aquasi-di ff usive approach in 1D to simulate the turbulent eddy-mixing processes in cloudy atmospheres, using the new D iffu -D rift models. This approach can be considered as the limitingcase of small-scale mixing. On the other extreme, large-scalehydrodynamic motions (convection, Hadley-cells, etc.) may beable to dredge up those elements maybe in a more immediatestraightforward way, which was the idea in our previous D rift models. In reality, there is not only vertical, but also horizon-tal mixing, which is likely to be very e ffi cient for example insuper-rotating horizontal jets as known from Jupiter (Schneider& Liu 2009), assuming that there are considerable horizontalabundance gradients present in the atmosphere. More 3D numer-ical experiments are required to quantify the e ffi ciency of mixingto inform our cloud formation models. Acknowledgements.
ChH thanks Will Best and Jonathan Gagné for the discus-sion on the number of known brown dwarfs. We thank Robin Baeyens and Lud-mila Carone for insightful discussions on mixing regimes in giant gas planets.The computer simulations were carried out on the UK MHD Consortium paral-lel computer at the University of St Andrews, funded jointly by STFC and SRIF.
References
Ackerman, A. S. & Marley, M. S. 2001, ApJ, 556, 872Allard, F. 2014, in IAU Symposium, Vol. 299, Exploring the Formation and Evo-lution of Planetary Systems, ed. M. Booth, B. C. Matthews, & J. R. Graham,271–272Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of theRoyal Society of London Series A, 370, 2765Apai, D., Radigan, J., Buenzli, E., et al. 2013, ApJ, 768, 121Arcangeli, J., Désert, J.-M., Line, M. R., et al. 2018, ApJL, 855, L30Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481Best, W. M. J., Magnier, E. A., Liu, M. C., et al. 2018, ApJS, 234, 1Birkby, J. L., de Kok, R. J., Brogi, M., Schwarz, H., & Snellen, I. A. G. 2017,AJ, 153, 138Brandenburg, A. 2016, ApJ, 832, 6Bringuier, E. 2013, European Journal of Physics, 34, 1103Buenzli, E., Marley, M. S., Apai, D., et al. 2015, ApJ, 812, 163Carone, L., Baeyens, R., Mollière, P., et al. 2019, arXiv e-prints,arXiv:1904.13334Charnay, B., Bézard, B., Baudino, J. L., et al. 2018, ApJ, 854, 172Dehn, M., Helling, C., Woitke, P., & Hauschildt, P. 2007, in IAU Symposium,Vol. 239, Convection in Astrophysics, ed. F. Kupka, I. Roxburgh, & K. L.Chan, 227–229Deuflhard, P. & Nowak, U. 1987, in Large scale scientific computing, ed. P. Deu-flhard & B. Engquist, Prog. Sci. Comp. 7 (Birkhäuser), 37–50Draine, B. T. & Lee, H. M. 1984, ApJ, 285, 89Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237Freytag, B., Allard, F., Ludwig, H.-G., Homeier, D., & Ste ff en, M. 2010, A&A,513, A19Gagné, J., Faherty, J. K., Cruz, K. L., et al. 2015, ApJS, 219, 33Gail, H.-P., Keller, R., & Sedlmayr, E. 1984, A&A, 133, 320Gail, H.-P., Wetzel, S., Pucci, A., & Tamanai, A. 2013, A&A, 555, A119Gibson, N. P., Nikolov, N., Sing, D. K., et al. 2017, MNRAS, 467, 4591Gierasch, P. J. & Conrath, B. J. 1985, in Recent Advances in Planetary Meteo-rology, ed. G. E. Hunt (Cambridge University Press), 121–146Helling, C. 2019, Annual Review of Earth and Planetary Sciences, 47, 583Helling, C., Ackerman, A., Allard, F., et al. 2008a, MNRAS, 391, 1854Helling, C. & Casewell, S. 2014, A&A Rev., 22, 80Helling, C., Dehn, M., Woitke, P., & Hauschildt, P. H. 2008b, ApJL, 675, L105Helling, C. & Fomins, A. 2013, Philosophical Transactions of the Royal Societyof London Series A, 371, 20110581Helling, C., Harrison, R. G., Honary, F., et al. 2016a, Surveys in Geophysics, 37,705Helling, C., Iro, N., Corrales, L., et al. 2019, arXiv e-prints, arXiv:1906.08127Helling, C., Rimmer, P. B., Rodriguez-Barrera, I. M., et al. 2016b, PlasmaPhysics and Controlled Fusion, 58, 074003Helling, C., Tootill, D., Woitke, P., & Lee, G. 2017, A&A, 603, A123Helling, C. & Woitke, P. 2006, A&A, 455, 325 Article number, page 13 of 19 & A proofs: manuscript no. Di ff uDrift Helling, C., Woitke, P., & Thi, W.-F. 2008c, A&A, 485, 547Hörst, S., He, C., Lewis, N., et al. 2019, in European Planetary Science Congress,Vol. 2019, EPSC–DPS2019–1095Huitson, C. M., Désert, J.-M., Bean, J. L., et al. 2017, AJ, 154, 95Juncher, D., Jørgensen, U. G., & Helling, C. 2017, A&A, 608, A70Kirk, J., Wheatley, P. J., Louden, T., et al. 2018, MNRAS, 474, 876Klein, R. 1995, Journal of Computational Physics, 121, 213Krasnokutski, S. A., Goulart, M., Gordon, E. B., et al. 2017, ApJ, 847, 89Lamb, D. & Verlinde, J. 2011, Physics and Chemistry of Clouds (CambridgeUniversity Press)Laor, A. & Draine, B. T. 1993, ApJ, 402, 441Lee, G., Dobbs-Dixon, I., Helling, C., Bognar, K., & Woitke, P. 2016, A&A,594, A48Lee, G., Helling, C., Dobbs-Dixon, I., & Juncher, D. 2015a, A&A, 580, A12Lee, G., Helling, C., Giles, H., & Bromley, S. T. 2015b, A&A, 575, A11Lee, G. K. H., Blecic, J., & Helling, C. 2018, A&A, 614, A126Lee, G. K. H., Wood, K., Dobbs-Dixon, I., Rice, A., & Helling, C. 2017, A&A,601, A22Leggett, S. K., Tremblin, P., Esplin, T. L., Luhman, K. L., & Morley, C. V. 2017,ApJ, 842, 118Lines, S., Mayne, N., Boutle, I., et al. 2018a, A&ALines, S., Mayne, N. J., Boutle, I. A., et al. 2018b, A&A, 615, A97Ludwig, H.-G., Allard, F., & Hauschildt, P. H. 2002a, A&A, 395, 99Ludwig, H.-G., Allard, F., & Hauschildt, P. H. 2002b, A&A, 395, 99Luhman, K. L. 2014, ApJL, 786, L18Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, ApJ, 756, 172Moses, J. I., Bézard, B., Lellouch, E., et al. 2000, Icarus, 143, 244Narita, N., Enomoto, T., Masaoka, S., & Kusakabe, N. 2015, Scientific Reports,5, 13977Nikolov, N., Sing, D. K., Gibson, N. P., et al. 2016, ApJ, 832, 191Ormel, C. W. & Min, M. 2019, A&A, 622, A121Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91Pino, L., Ehrenreich, D., Wyttenbach, A., et al. 2018, A&A, 612, A53Powell, D., Zhang, X., Gao, P., & Parmentier, V. 2018, ApJ, 860, 18Rimmer, P. B. & Helling, C. 2016, ApJS, 224, 9Riols, A. & Lesur, G. 2018, A&A, 617, A117Schaaf, S. A. 1963, Handbuch der Physik, 3, 591Schneider, T. & Liu, J. 2009, Journal of Atmospheric Sciences, 66, 579Schräpler, R. & Henning, T. 2004, ApJ, 614, 960Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature, 529, 59Tregloan-Reed, J., Southworth, J., Mancini, L., et al. 2018, MNRAS, 474, 5485Tremblin, P., Padioleau, T., Phillips, M. W., et al. 2019, ApJ, 876, 144Tsuji, T. 2002, ApJ, 575, 264Tsuji, T., Ohnaka, K., & Aoki, W. 1996, A&A, 305, L1Witte, S., Helling, C., Barman, T., Heidrich, N., & Hauschildt, P. H. 2011, A&A,529, A44Witte, S., Helling, C., & Hauschildt, P. H. 2009, A&A, 506, 1367Woitke, P. & Helling, C. 2003, A&A, 399, 297Woitke, P. & Helling, C. 2004, A&A, 414, 335Woitke, P., Helling, C., Hunter, G. H., et al. 2018, A&A, 614, A1Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103Yang, H., Apai, D., Marley, M. S., et al. 2016, ApJ, 826, 8Youdin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588Zahnle, K., Marley, M. S., Morley, C. V., & Moses, J. I. 2016, ApJ, 824, 137Zhang, X. & Showman, A. P. 2018, ApJ, 866, 1Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73
Article number, page 14 of 19eter Woitke et al.: Dust in brown dwarfs and extra-solar planets
Appendix A: The diffusion solver
We use a self-developed 1D explicit / implicit di ff usion solver inthis paper which has second order accuracy in both the formula-tion of the di ff erential equations and the boundary conditions .In case of a plane-parallel static atmosphere, the di ff usion equa-tion for element k is given by d ( n (cid:104) H (cid:105) (cid:15) k ) dt = ddz (cid:32) n (cid:104) H (cid:105) D gas d (cid:15) k dz (cid:33) (A.1)where the di ff usive element flux is j di ff k = − n (cid:104) H (cid:105) D gas d (cid:15) k dz (A.2) Appendix A.1: Vertical grid and discretisation of derivatives
We introduce an ascending vertical grid z i ( i = , ... , I ). The firstand second derivatives of any quantity f ( z i ) = f i at grid point z i are approximated as d f i dz = d l , i f i − + d m , i f i + d r , i f i + (A.3) d f i dz = d l , i f i − + d m , i f i + d r , i f i + , (A.4)i.e. as linear combinations of the function values on the neigh-bouring grid points, where e.g. d l , i is the coe ffi cient for the firstderivative on the point left of the grid point i , d m , i the same onthe mid point and d r , i the same on the point right of grid point i .Similar, for the second derivative, the coe ffi cients are d l , i , d m , i and d r , i . Using a second order polynomial approximation forfunction f ( z ) the coe ffi cients are given by d l , i = − h ri ( h ri + h li ) h li (A.5) d m , i = + h ri − h li h li h ri (A.6) d r , i = + h li ( h r ( i ) + h li ) h ri (A.7) d l , i = + h ri + h li ) h li (A.8) d m , i = − h ri h li (A.9) d r , i = + h ri + h li ) h ri (A.10)where h li = z i − z i − and h ri = z i + − z i are the l.h.s. and the r.h.s.grid point distances. For the special case of an equidistant grid,we have h = h li = h ri and hence d f i dz = f i + − f i − h (A.11) d f i dz = f i + − f i + f i − h (A.12) The code is available at https://github.com/pw31/Diffusion . The above equations are valid for grid points i = , ... , I −
1. Forthe first derivative at the boundaries we write d f dz = d l , f + d m , f + d r , f (A.13) d f dz = d l , I f I − + d m , I f I − + d r , I f I (A.14)which is also second order accuracy by using the informationon the 3 leftmost or 3 rightmost grid points, respectively. Thecoe ffi cients are given by d l , = − h + h h h (A.15) d m , = h h ( h − h ) (A.16) d r , = − h h ( h − h ) (A.17) d r , I = h I − + h I − h I − h I − (A.18) d m , I = − h I − h I − ( h I − − h I − ) (A.19) d l , I = h I − h I − ( h I − − h I − ) (A.20)where h = z − z , h = z − z , h I − = z I − z I − and h I − = z I − z I − . Appendix A.2: Spatial derivatives
The di ff usion term at grid point z i ( i = ... I −
1) is numericallyresolved, with abbreviation D gas ( z i ) = D i , as ddz (cid:32) n (cid:104) H (cid:105) D gas d (cid:15) k dz (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z i = d (cid:0) n (cid:104) H (cid:105) D gas (cid:1) dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z i · d (cid:15) k dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z i + n (cid:104) H (cid:105) D gas d (cid:15) k dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z i = (cid:16) d l , i n (cid:104) H (cid:105) , i − D i − + d m , i n (cid:104) H (cid:105) , i D i + d r , i n (cid:104) H (cid:105) , i + D i + (cid:17) · (cid:16) d l , i (cid:15) k , i − + d m , i (cid:15) k , i + d r , i (cid:15) k , i + (cid:17) + n (cid:104) H (cid:105) , i D i (cid:16) d l , i (cid:15) k , i − + d m , i (cid:15) k , i + d r , i (cid:15) k , i + (cid:17) (A.21)and the di ff usive fluxes across the lower and upper boundariesare φ k , = − D n (cid:104) H (cid:105) , (cid:16) d l , (cid:15) k , + d m , (cid:15) k , + d r , (cid:15) k , (cid:17) (A.22) φ k , I = − D I n (cid:104) H (cid:105) , I (cid:16) d l , I (cid:15) k , I − + d m , I (cid:15) k , I − + d r , I (cid:15) k , I (cid:17) . (A.23) Appendix A.3: Boundary conditions
As boundary conditions, we have implemented three options, forexample considering the lower boundary:1. fixed concentration: (cid:15) k , is a given constant2. fixed flux: φ k , is a given constant3. fixed outflow rate: The flux through a boundary is assumedto be proportional to the concentration of species k at theboundary, e.g. φ k , = β k n (cid:104) H (cid:105) , (cid:15) k , v k [cm − s − ] (A.24)where the β k is a given probability (fixed value) and v k is thespeed at which the particles of kind k are moving through theboundary (also fixed value). Article number, page 15 of 19 & A proofs: manuscript no. Di ff uDrift Appendix A.4: Explicit integration
A straightforward way to integrate Eq. (A.1), for a timestep ∆ t ,is the following explicit scheme f ( n ) i = f ( n − i + ∆ t d f ( n − i dt (A.25)where f ( n ) i is some quantity on grid point i at time t n and f ( n − i is the quantity on grid point i at time t n − with t n = t n − +∆ t . Inconsideration of Eq. (A.1), this leads to (cid:15) ( n ) k , i = (cid:15) ( n − k , i + ∆ t I (cid:88) j = A i j (cid:15) ( n − k , j , (A.26)where A is a tri-diagonal matrix, the elements A i j of which aregiven by Eq. (A.21). Equation (A.26) applies to the grid points i = , ... , I −
1, but not to the boundaries. On the boundary points,the following equations are applied depending on the choice ofboundary conditions, here for example the lower boundary1. fixed concentration: (cid:15) ( n ) k , = (cid:15) k
2. fixed flux: (cid:15) ( n ) k , = d l , (cid:32) − φ k , D n (cid:104) H (cid:105) , − d m , (cid:15) ( n ) k , − d r , (cid:15) ( n ) k , (cid:33)
3. fixed outflow rate: β k n (cid:104) H (cid:105) , (cid:15) k , v k = − D n (cid:104) H (cid:105) , (cid:16) d l , (cid:15) k , + d m , (cid:15) k , + d r , (cid:15) k , (cid:17) ⇒ (cid:15) ( n ) k , = − d m , (cid:15) ( n ) k , − d r , (cid:15) ( n ) k , d l , + β k v k D . These assignments are applied at time t n , i.e. after an explicitdi ff usion timestep has been completed on grid points i = ... I −
1. To guarantee numerical stability, the explicit timestep must belimited by α ≤ . ∆ t = α · min i = , ... , I ( z i − z i − ) ( D i + D i − ) . (A.27) Appendix A.5: Implicit integration
To avoid the timestep limitations in the explicit solver, and toguarantee numerical stability for much larger timesteps, an im-plicit integration scheme can optionally be applied f ( n ) i = f ( n − i + ∆ t d f ( n ) i dt (A.28)which is a system of linear equations for the unknowns f ( n ) i . Inconsideration of Eq. (A.1), we have (cid:15) ( n ) k , i = (cid:15) ( n − k , i + ∆ t I (cid:88) j = A i j (cid:15) ( n ) k , j (A.29)We re-write this equation more generally, including the bound-ary conditions, by means of another unit-free matrix as B (cid:15) ( n ) k = R k , (A.30)where we have B i j = ( − ∆ t A ) i j and R k , i = (cid:15) ( n − k , i for i = ... I − z [cm] c o n ce n t r a t i o n t = 0.000 st = 0.010 st = 0.025 st = 0.050 st = 0.100 st = 0.200 st = 0.400 s Fig. A.1.
Test problem with fixed concentrations on the left and rightboundaries ( (cid:15) = (cid:15) ( z , t ) = exp( − ω t ) sin( kz ) with k = π and ω = Dk .
1. fixed concentration: B = R k , = (cid:15) k
2. fixed flux: B = B = d m , / d l , , B = d m , / d l , , and R k , = − φ k , / (cid:0) n (cid:104) H (cid:105) , D d l , (cid:1)
3. fixed outflow rate: B = + β k v k / (cid:0) D d l , (cid:1) , B = d m , / d l , , B = d m , / d l , , and R k , = (cid:15) ( n ) k = B − R k (A.32)where B − is the inverse of the matrix B . As long as the spa-tial grid points z i , the densities n (cid:104) H (cid:105) , i and di ff usion constants D i ,the constants involved in the boundary conditions (e.g. φ k , or β k ), and the timestep ∆ t do not change, we need to perform thematrix inversion only once. Successive time steps are then per-formed by simply incrementing n , re-computing the vector R k ,and applying again Eq. (A.32). B − is also usually the same forall elements k to be di ff used.This favourable property of B makes the computation of im-plicit timesteps actually very fast. We note, however, that B − ,in general, is a full I × I matrix where all entries are positive( B − ) i j >
0. This leads to a very stable numerical behaviour forarbitrary time steps. In contrast, the matrix A has positive entriesalong the main diagonal, but negative entries along both semi-diagonals, which leads to numerical instabilities when the timestep is too large. Appendix B: The settling solver
For the 1D vertical settling in the Epstein regime we solve d ( ρ L j ) dt = ddz (cid:90) V (cid:96) v ◦ dr ( V ) f ( V ) V j / dV = ddz (cid:32) ξ ρ d c T L j + (cid:33) , (B.1)according to Eqs. (6) and (8). The settling flux for cloud particlemoment ρ L j is hence φ j = − ξ ρ d c T L j + = − ρ L j v dr , j (B.2) Article number, page 16 of 19eter Woitke et al.: Dust in brown dwarfs and extra-solar planets z [cm] c o n ce n t r a t i o n t = 0.000 st = 0.001 st = 0.001 st = 0.002 st = 0.005 st = 0.010 st = 0.020 st = 0.050 st = 0.100 st = 0.300 s Fig. A.2. Di ff usive evolution of an initial δ -peak with analytic solutionoverplotted. The analytic solution is (cid:15) ( z , t ) = A ( t ) exp (cid:0) − ( z − . / w ( t ) (cid:1) with A ( t ) = √ t / t and w ( t ) = √ Dt . where we have introduced mean drift velocities for the cloudparticle moments ρ L j asv dr , j = ξ ρ d ρ c T L j + L j . (B.3)The cloud particle moments are updated according to the fol-lowing explicit upwind advection scheme. We first calculate allvertical moment fluxes φ j , i = φ j ( z i ) via Eq. (B.2) and then apply ρ L ( n ) j , i = ρ L ( n − j , i + ∆ t ∆ z (cid:16) φ ( n − j , i + − φ ( n − j , i (cid:17) (B.4) ρ L s , i ( n ) = ρ L s , i ( n − + ∆ t ∆ z (cid:16) b si + n − φ ( n − , i + − b si ( n − φ ( n − , i (cid:17) , (B.5)where the notation f ( n ) i is some quantity on grid point i at time t n , ∆ t = t n − t n − the timestep and ∆ z the vertical extension of theconsidered atmospheric cell. Appendix C: Verification tests
We have carefully checked our di ff usion solver against analyticaltest problems and by cross-checking the results from the explicitand implicit solvers. Figures A.1 and A.2 show two test problemson domain z = [0 , n (cid:104) H (cid:105) = D =
1. Theblack thin lines are the overplotted analytic solutions, showingexcellent agreement. The tests use an equidistant z -grid with 101points, and can be computed within less than 1 CPU-sec.The convergence of the full cloud formation model was stud-ied by comparing the results obtained with di ff erent initial condi-tions (see Sect. 4.5). Figure A.3 shows the results of two modelsfor T e ff = K , log g = β (cid:48) = ff ective at high densi-ties and low temperatures, causing transient minimums of (cid:15) k ( z ).The initially ‘empty’ model needs more time to start formingclouds because the condensable elements first need to be trans-ported upwards by di ff usion, resulting in a more gradual onsetof cloud formation. The final states after t = (cid:15) k ( z ) decrease monotonically withheight and have zero gradients at the upper boundary, as it shouldbe in the time-independent case, see Sect. 3.1. Appendix D: Diffusion coefficients in the literature Di ff usion, in principle, is a microscopic process driven by par-ticle concentration gradients ∇ c j , where for example c CO = n CO / n tot for j = CO and n tot = (cid:80) n j is the total particle density.Such gradients can result from gravity (Zahnle et al. 2016), fromchemical processes (Moses et al. 2000) and from cloud conden-sation as shown in this paper. Di ff usion will always counteractthese concentration gradients. Experiments have been conductedto measure di ff usion constants for gases relevant for solar systemplanets. Lamb & Verlinde (2011) provide, for example, the gas-kinetic di ff usion coe ffi cient for water molecules near sea level inthe Earth atmosphere ( D H2O ≈ × − cm / s) and note that itvaries inversely proportional with the atmospheric pressure.The e ff ect of mixing on larger scales has been modelled dif-ferently in di ff erent communities, and terminology is usuallynot unique. Often, a quasi-di ff usive approach is used where thedi ff usion constant is replaced by a function of height or den-sity / pressure. Transport of matter due to turbulent mixing hasbeen termed ’turbulent di ff usion’ in protoplanetary disk mod-elling, describing the averaged e ff ect of advection of the indi-vidual turbulent eddies (e.g. Schräpler & Henning 2004) and as’eddy di ff usion’ in planetary atmosphere modelling.Studying solar system giants, Moses et al. (2000) havedemonstrated that ISO observations of hydrocarbon moleculesin Saturn’s atmosphere can be well fitted by assuming an eddydi ff usion coe ffi cient of D mix = . × cm / s (cid:32) . × cm − n tot (cid:33) β (D.1)with slope β between 0.3 and 0.7, i.e. D mix increases with height.Ackerman & Marley (2001) consider an equilibrium be-tween upward mixing of vapour in the gas phase and gravi-tational settling of particles condensed from the vapour. TheirEq. (4) reads − D mix ddz ( q c + q v ) − f sed w (cid:63) q c = q v and q c are the mixing ratios of vapour and conden-sate, respectively (moles of vapour / condensate per mole of at-mosphere). f sed w (cid:63) represents an average sedimentation velocityof the condensed particles with f sed being adjusted as needed.We note that Eq. (D.2) is very similar to our Eq. (29) for thecase where cloud and gas particles are equally a ff ected by eddy-di ff usion. Their eddy-di ff usion coe ffi cient D mix is defined ac-cording to (Gierasch & Conrath 1985) as D mix = H (cid:32) (cid:96) mix H (cid:33) / (cid:32) RF conv µ ρ c p (cid:33) / (D.3)where H = RT / ( µ g ) is the pressure scale height, (cid:96) mix the mix-ing length, c p = ( f + R / (2 µ ) is the isobaric specific heat, f isthe mean degree of freedom of the gas particles, R the universalgas constant, and F conv is the convective heat flux. Ackerman &Marley (2001) assume F conv = σ T ff , i.e. they assume that the at-mosphere is fully convective, which leaves open the problem ofwhat to do in radiative layers, for example in brown dwarf atmo-spheres. Charnay et al. (2018) note that the factor 1 / l mix is calculated as l mix = H · max { Λ , Γ / Γ adb } (D.4) Article number, page 17 of 19 & A proofs: manuscript no. Di ff uDrift l og † ga s ( a ) t =0 days $ O &