Issue 
A&A
Volume 635, March 2020



Article Number  A159  
Number of page(s)  9  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/201937063  
Published online  30 March 2020 
The effect of internal gravity waves on cloud evolution in substellar atmospheres
^{1}
Division of Computing and Mathematics, Abertay University,
Kydd Building,
Dundee DD1 1HG,
UK
email: 1303985@abertay.ac.uk; c.stark@abertay.ac.uk
^{2}
Mathematics, School of Science & Engineering, University of Dundee,
Nethergate,
Dundee DD1 4HN, UK
^{3}
Atmospheric, Oceanic and Planetary Physics, Department of Physics, University of Oxford,
Oxford OX1 3PU, UK
Received:
5
November
2019
Accepted:
23
February
2020
Context. Substellar objects exhibit photometric variability, which is believed to be caused by a number of processes, such as magneticallydriven spots or inhomogeneous cloud coverage. Recent substellar models have shown that turbulent flows and waves, including internal gravity waves, may play an important role in cloud evolution.
Aims. The aim of this paper is to investigate the effect of internal gravity waves on dust nucleation and dust growth, and whether observations of the resulting cloud structures could be used to recover atmospheric density information.
Methods. For a simplified atmosphere in two dimensions, we numerically solved the governing fluid equations to simulate the effect on dust nucleation and mantle growth as a result of the passage of an internal gravity wave. Furthermore, we derived an expression that relates the properties of the waveinduced cloud structures to observable parameters in order to deduce the atmospheric density.
Results. Numerical simulations show that the density, pressure, and temperature variations caused by gravity waves lead to an increase of the dust nucleation rate by up to a factor 20, and an increase of the dust mantle growth rate by up to a factor 1.6, compared to their equilibrium values. Through an exploration of the wider substellar parameter space, we show that in absolute terms, the increase in dust nucleation due to internal gravity waves is stronger in cooler (T dwarfs) and TiO_{2}rich substellar atmospheres. The relative increase, however, is greater in warm (L dwarf) and TiO_{2}poor atmospheres due to conditions that are less suited for efficient nucleation at equilibrium. These variations lead to banded areas in which dust formation is much more pronounced, similar to the cloud structures observed on Earth.
Conclusions. We show that internal gravity waves propagating in the atmosphere of substellar objects can produce banded clouds structures similar to that observed on Earth. We propose a method with which potential observations of banded clouds could be used to estimate the atmospheric density of substellar objects.
Key words: brown dwarfs / stars: atmospheres / waves / planets and satellites: atmospheres
© ESO 2020
1 Introduction
Brown dwarfs are lowmass, substellar objects below the hydrogen burning limit, with masses between 13 M_{Jup} and 70 M_{Jup}. As a consequence, their atmospheres are sufficiently cool for the formation of dust clouds. The clouds are observed through their effect on spectral features and their association with infrared spectroscopic variability, which is believed to be caused by patchy clouds. Numerous observations (see for example Buenzli et al. 2014) show that a large portion of known brown dwarfs exhibit photometric variability. According to Biller (2017), over 10% of known brown dwarfs show a variation of 1% or more, and over 50% exhibit a variation of 0.1–0.5% or more.
Explaining spectral variability is necessary to understanding the L/T transition (Vos et al. 2019). The L and T components of the Luhman 16AB system, for example, show vastly different patterns (Gillon et al. 2013): Luhman 16B, a T dwarf, exhibits strong, fastchanging periodic spectral variations, while the L dwarf Luhman 16A exhibits no periodic pattern. A model from Saumon & Marley (2008) and Marley et al. (2010) proposes the sinking of parts of the cloud deck, creating thinner, patchy clouds as an explanation for the change in variability patterns around the L/T transition. A study by Stark et al. (2015) proposes the electrostatic disruption of charged cloud particles as a mechanism through which inhomogeneous coverage could be caused. While inhomogeneous dust cloud coverage (Helling & Casewell 2014) is believed to be the main cause for brown dwarf variability, other theories, such as temperature variations (Robinson & Marley 2014) or fingering convection (Tremblin et al. 2016), propose cloudfree models as an explanation.
Internal gravity waves have been simulated in main sequence stars (Alvan et al. 2014). Internal gravity waves triggered by fingering convection have been modelled in objects ranging from mainsequence stars to brown dwarfs (Garaud et al. 2015), and they can reach wavelengths much larger than the source perturbation. Simulations of atmosphere patches by Freytag et al. (2010) show that internal gravity waves are also present in substellar atmospheres, triggered by downdrafts caused by convection patterns, and they are theorised to be one of the main phenomena responsible for transporting dust in the upper atmospheric layers.
Further characterising the impact of gravity waves on dust cloud evolution can advance the understanding of cloud structures in brown dwarfs. Helling et al. (2001) showed that higher frequency acoustic waves, triggered by turbulent flow in brown dwarfs atmospheres, can have a strong impact on cloud formation: by carrying lower temperature perturbations, the passage of waves can temporarily create favourable conditions for dust nucleation in otherwise dusthostile environments, leading to the formation of dust over time.
Internal gravity waves are a type of fluid wave that occurs in atmospheres and oceans. Their defining characteristic is that gravity, in the form of buoyancy, is the restoring force that allows disturbances to propagate. Internal gravity waves can be observed in the Earth’s atmosphere through their effect on clouds. On Earth, a wave cloud is formed from the passage of an internal gravity wave, triggered by stable air flowing over relief. The vertical displacement of the air forces it to oscillate as the buoyancy force tries to restore equilibrium. As the wave propagates, at the wave peaks the displaced air rises and cools resulting in water vapour condensing, forming droplets and clouds; at the wave troughs, the clouds evaporate due to adiabatic heating, leading to clouds that have a distinct banded structure.
In the case of a gas giant planet or brown dwarf, dust clouds are formed instead of water clouds but an analogous process can occur. In this context, wave clouds can be induced as a result of external fluid motion triggering turbulent flow (such as fingering convection in deeper layers of the atmosphere), whereas gas flow over relief would be expected to be the main cause in the case of a rocky terrestrial exoplanet (Roeten et al. 2019). In a substellar atmosphere, oscillating parcels of gas can trigger the nucleation of seed particles and enhanced surface mantle growth, forming banded cloud structures. The nucleation rate is a function of density of the nucleating species and the atmospheric temperature. If the passage of the internal gravity wave perturbs the local thermodynamic structure of the atmosphere it can give enhanced nucleation in localised regions.
The aim of this paper is to investigate and characterise the effect of internal gravity waves on the evolution of dust clouds in the atmospheres of substellar objects and its consequences for cloud variability. This paper presents a novel mechanism for potentially diagnosing the gas density of substellar atmospheres from observations of the resulting cloud structures formed from the passage of internal gravity waves. In Sect. 2 the basic atmos pheric model of internal gravity waves, nucleation and mantle growth is described; in Sect. 3 the numerical methods used to simulate internal gravity waves are presented; in Sect. 4 the results of the simulations are presented and discussed; Sect. 5 summarises and discusses the consequences of the results including a possible way of connecting observations to the wave dispersion relation to diagnose the atmospheric density.
2 Substellar internal gravity waves
For a vertical slice of a substellar atmosphere in hydrostatic equilibrium, the coupled equations of fluid dynamics governing the evolution of the fluid velocity u, the fluid density ρ, and the pressure p, of an atmospheric parcel, under the influence of gravity g are:
where γ_{a} is the ratio of specific heats (for a diatomic gas γ = 7∕5). In static equilibrium u_{0} = 0 and ∂∕∂t = 0, giving the following equilibrium relationship: (4)
where the subscript “0” denotes an equilibrium quantity. For simplicity, in order to capture the fundamental physics, this paper focuses on the effect of gravity waves in the linear regime. We can linearise Eqs. (1)–(3) by decomposing each variable Q into its equilibrium and perturbed value, so that Q = Q_{0} + Q_{1}. In the nonlinear regime Q_{0} ≫ Q_{1}, and powers of Q_{1} higher than 1 can be discarded. For clarity, the subscript “1” is omitted in further equations. Linearisation yields the final system of fluid equations:
To model internal gravity waves, where we deal with incompressible flows ∇ ⋅u = 0, we adopted a vorticitystream function formulation by introducing the vorticity ζ and stream function ψ defined by,
Therefore, the governing fluid equations for incompressible flows in the linear regime become,
where the baroclinic term (∇ρ ×∇p∕ρ^{2}) vanishes since the propagation of internal gravity waves is considered to be an adiabatic process. In an atmospheric vertical plane (x, y):
This yields the system of equations
To obtain the gravity waves’ dispersion relation, Eq. (18) is derived with respect to time, and Eq. (16) substituted for ∂ζ∕∂t: (20)
Differentiating Eq. (20) with respect to time, and then substituting Eq. (17) for ∂ρ∕∂t yields: (21)
At equilibrium, ρ_{0} does not vary along x, therefore: (22)
is the BruntVäisälä buoyancy frequency. In the case of a uniformly stratified atmosphere, assuming a solution of the form ψ ≈ exp[−i(ωt + k_{x}x + k_{y}y)] the dispersion relation for internal gravity waves becomes (Sutherland 2010; Vallis 2017),
Therefore, when the atmospheric density (or equivalently velocity) is perturbed, corresponding oscillations in ρ, ψ, and ζ are triggered, that occur at the BruntVäisälä buoyancy frequency N. As the wavepropagates through the atmosphere, the density variations can affect the resulting nucleation and mantle growth rates.
2.1 Dust nucleation
To quantify the impact of passing waves on dust formation, we used the equation of modified classical nucleation theory presented by Gail et al. (1984); Helling et al. (2001), which defines the nucleation rate J_{*}, the numberof nucleating centres formed per second per unit volume [m^{−3} s^{−1}], as: (26)
where T is the temperature; τ is the seed growth time scale for the gaseous nucleation species x; N_{*} is the size of the critical cluster; n_{x} is the number density of the nucleating species; Z is the Zeldovich factor; S is the supersaturation ratio; defined as follows:
and p_{x} is the partial pressure of the nucleating species x; p_{sat,x} is the saturation vapour pressure of the nucleating species x; σ is the surface tension of the nucleating species; r_{0} is the hypothetical monomer radius; is the hypothetical monomer surface area; m_{x} is the mass of a monomer particle; m_{p} is the proton mass; A is the atomic weight of the monomer; ρ_{m} is the dust material density; and v_{rel} is the thermodynamic equilibrium velocity for the nucleating species studied (TiO_{2} for this paper).
For this paper, we assumed a constant value for the surface tension of TiO_{2}, (see Helling et al. 2001; Lee et al. 2015), and a density ρ_{m} ≈4230 kg m^{−3} for TiO_{2}. In the context of an internal gravity wave propagating through a substellar atmosphere, the wave perturbs the local atmospheric gas density and hence the temperature (via Tρ^{1−γa} = const.) in an adiabatic process. As a result, the passage of the wave perturbs the nucleation rate J_{*}.
2.2 Mantle growth
Once nucleation has established a material surface onto which material can accumulate, dust growth occurs via gasphase surface chemistry (Eq. (24) in Helling & Woitke 2006). We consider a spherical dust grain of radius a, of mass m_{d}, and let n be the number density of the gas phase. The dust grain absorbs gas molecules at a rate απa^{2} n⟨v⟩, where ⟨v⟩ is the mean gas molecular speed and α is the sticking probability that a molecule is absorbed (the sticking factor). Therefore, the mass of the dust grain m_{d} evolves in time as (35)
where ρ = nm is the gas mass density. The mass of a dust grain can be written as , where ρ_{m} is assumed to be constant. Therefore, the time evolution of the radius of a dust grain is (36)
where γ is the growth rate from absorption in units of [ms^{−1}]. Eq. (36) is the archetypal equation describing the absorption of material onto the surface of a dust grain. It is consistent with the dust growth equations presented in Helling et al. (2001), albeit in a much simplified form but still encapsulating the fundamental underlying physics. Furthermore, Eq. (36) is also consistent with mantle growth via ion accretion when dust grains are immersed in a plasma (Eq. (18) in Stark & Diver 2018). Without loss of generality, to investigate the effect of internal gravity waves on the mantle growth rate we simplify the expression by introducing ρ_{s}, the density of the gasphase accreting species, defined as ρ_{s} = f_{s} ρ, where 0 ≤ f_{s} ≤ 1 is the fraction of the surrounding gas composed of the accreting species, (37)
where we have set α = 1 to obtain the optimal growth rate; and ⟨v⟩ = v_{rel,s}. Introducing f_{s} allows us to generalise the effect of different gasphase species, with varying relative abundances in the gasphase, participating in surface chemistry leading to mantle growth.
3 Numerical simulations
3.1 Model equations
The linearised governing equations can be cast in nondimensional form, defining
where L, T, and M are characteristic values of length, time, and mass respectively. Therefore, in nondimensional form, Eqs. (16)–(18) become
Casting the model equations in nondimensional form lets us observe the characteristic behaviour of internal gravity waves without lossof generality.
3.2 Methods
To solve the system of fluid Eqs. (45)–(47) numerically, we used a combination of the leapfrog method for Eqs. (45) and (46), and Successive OverRelaxation (SOR) for Eq. (47) (Vetterling et al. 1992; Mittal 2014). SOR requires that the value of ζ from Eq. (18) be known at the boundaries of the domain. In order to minimise the artefacts caused by boundaries, we used a domain large enough that over short timescales, the waves’ perturbations do not reach the boundaries. Additionally, values of ρ, ψ, ζ, and their spatial firstorder derivatives were interpolated using thirdorder polynomial interpolation. The internal gravity waves were triggered bycreating a Gaussian density perturbation initial condition in the centre of the numerical domain, that was modulated in time by a sine wave with a period equal to a multiple of the local buoyancy frequency, (48)
where ρ_{A} is the maximum amplitude of the perturbation; ς is the spread parameter; and r is the distance to the centre of the numerical domain. The resulting wave solutions were used to calculate the nucleation (Eq. (26)) and dust mantle growth (Eq. (37)) rates as a result of the propagating internal gravity waves. The simulation parameters are presented in Table 1.
3.3 Substellar atmospheric model
The aim of this paper was to investigate the effects of internal gravity waves on the evolution of dust clouds in substellar atmospheres. To this end, we considered a brown dwarf atmosphere characterised by T_{eff}= 1500 K and logg = 5.0 (see Fig. 1) as an exemplar substellar atmosphere, and we used data published in Stark et al. (2013) as input for our numerical simulations. Figure 1 and Table 2 also give the profiles of typical and lowgravity L and T dwarfs for context, to place the simulations in a wider range of brown dwarfs examples. The atmosphericdata was generated by the DRIFTPHOENIX model atmosphere and cloud formation code (Hauschildt & Baron 1999; Helling et al. 2004; Helling & Woitke 2006; Witte et al. 2009, 2011). We note that the atmospheric extension y (Fig. 1; middle panel) is measured from the top of the atmosphere.
The atmosphere model used as input is onedimensional; for this study, we assumed a horizontally uniform atmosphere at equilibrium, and expanded the model to two dimensions. We chose the characteristic length L (see Sect. 3.2) to be the height of the simulation domain (10^{4} m). The numerical domain was defined by the geometric extension y ∈ [25 km, 35 km]; the pressure p_{0} ∈ [3 × 10^{−5} bar, 3 × 10^{−3} bar]; the temperature T_{0} ∈ [700 K, 850 K]; and the density ρ_{0} ∈ [10^{−7} kg m^{−3}, 10^{−5} kg m^{−3}]. For these conditions, (see Fig. 2). The perturbation length scale was of the order of 10^{3} m.
Parameters used for numerical simulations.
4 Results
In a stratified substellar atmosphere, perturbing the background density by vertically displacing a fluid parcel triggers vertical oscillations as the buoyancy force tries to restore equilibrium. The resulting density variations propagate through the atmosphere at the BruntVäisälä buoyancy frequency. Figure 3 shows the buoyancy period T_{N} = 2π∕N as a function of pressure p for the substellar atmosphere characterised by T_{eff} = 1500 K and logg = 5.0. The buoyancy period is of the order of 10 s across the extent of the atmosphere, with the maximum period occurring at high atmosphericpressures. Also shown in Fig. 3, for context, are the buoyancy periods profiles of typical L and T dwarfs (periods of the order of 10–100 s), and lowgravity L and T dwarfs (longer buoyancy periods, of the order of 1000 s). In comparison, Buenzli et al. (2014) observe spectroscopic variations on timescales of 100–1000 s. The key parameter distinguishing between the buoyancy period for different atmospheric models is the surface gravity; since T_{N} ∝ g^{−1∕2}, objects withlower surface gravity will have a longer buoyancy periods. In the case of a neutrally stratified atmosphere (i.e. N = 0), the potential temperature is constant with atmospheric height and internal gravity waves cannot propagate (Sutherland 2010). This scenario is expected deeper in the atmosphere at higher gas pressures (for example, see Tremblin et al. 2015, 2019).
The wavelength of the internal gravity wave is determined by the spatial lengthscale of the instigating perturbation driving the oscillation, and consequently sets the speed of propagation. If the initial density variation is driven by an external source, such as convective motions or largescale turbulent motions, the frequency of the wave is set by thefrequency of the source for frequencies below the buoyancy frequency. In substellar atmospheres the length scale of convective motions deep in the atmosphere is related to the atmospheric pressure scale height by a factor between 1 and 10^{2} (Marley & Robinson 2015; Tremblin et al. 2019)– for example, l_{conv} ≈ 10^{3} km in Freytag et al. (2010) – varying with timescales of the order of 10^{3} s (Tremblin et al. 2019).
Figure 4 shows the characteristic St Andrews cross pattern emanating from density perturbation located at the centre of the numerical domain, for the example driving frequency ω =0.25N and a driving perturbation amplitude ρ_{A} = 0.2ρ_{0}. In an adiabatic process the density variations are accompanied by sympathetic variations in temperature (T ∝ ρ^{γa−1} = ρ^{2∕5}) that go on to affect the local nucleation rate and surface mantle growth rate. Figure 4 shows that areas where ρ is reduced by the passage of internal gravity waves are the areas with the largest increase of dust nucleation. The number density n_{x} of the nucleating species is a component of the total gas density ρ = ∑_{s}m_{s}n_{s} (see Fig. 2, showing the relationship between n and ). Variations of ρ propagated by passing gravity waves are therefore reflected in the partial density of the nucleating species , which in turns leads to an adiabatic change in temperature both of the total gas phase and its components. The primary effect that drives the variation in nucleation rate is the temperature variation: as the density increases (decreases) there is a corresponding adiabatic increase (decrease) in the temperature. Increasing (decreasing)the temperature of the nucleating species in this way, decreases (increases) the supersaturation ratio, S. Efficient nucleation is only possible for temperatures T_{0} such that S(T_{0}) ≫ 1 (Helling et al. 2001); therefore, increasing (decreasing) the temperature can yield a corresponding decrease (increase) in the supersaturation ratio, S, and hence a decrease (increase) in the nucleation rate. The nucleation rate is strongly dependent on the supersaturation ratio and hence the temperature; therefore, local overdensities in the atmosphere as a result of an internal gravity wave can decrease the local nucleation rate; whereas, local underdensities can increase the nucleation rate.
If there is an established particulate onto which material from the atmospheric gas can be absorbed, the growth rate can be also be affected by the passage of an internal gravity wave. In this scenario, an increase (decrease) in the local density, then the number of gas particles passing through a target area per unit time also increases (decreases), and so does the number of interactions per unit time that occurs between the seed particulates and the atmospheric species.
To quantify the impact of the waves on dust nucleation and growth, we ran simulations for a range of driving frequencies and perturbation amplitudes. The density, nucleation and growth responses are presented in Fig. 5. To obtain results comparable across varying frequencies, the values used for plotting were measured for τ = T_{ω}, where T_{ω} is the period of the driving oscillations, and normalised to their values at static equilibrium in the centre of the numerical domain. In order to obtain results comparable across a range of frequencies, the measurements were taken as the maximum values for ρ_{1}, J_{*,1}, and γ_{1} along a vertical slice of the numerical domain, located at a distance X(ω) from the centre of the domain so that X(ω) = X_{ref}ω_{ref}∕ω, where X_{ref} is the location of the slice picked for a reference case at ω_{ref}.
As the perturbation amplitude increases the density amplitude of the resulting wave response linearly increases by up to a factor 1.5, consistent with the linear regime assumed (top panel, Fig. 5). In contrast, the resulting nucleation rate (middle panel, Fig. 5) increases nonlinearly with the perturbation amplitude by up to a factor 20 in the most favourable scenario (ω = 0.25N, ρ_{A} = 0.2ρ_{0}). This is a result of the complex nonlinear dependence of the nucleation rate (Eq. (26)), and not as a consequence of a nonlinear evolution of the internal gravity wave, since our simulations are conducted in the linear regime only. This implies that if the internal gravity wave were to evolve nonlinearly, the corresponding nucleation rate could exhibit an enhanced nonlinear response, giving greater nucleation rate values. A similar assertion can be made regarding the mantle growth rate (bottom panel, Fig. 5): the mantle growth rate increases by up to a factor 1.6 as a nonlinear function of the perturbed density, γ ∝ ρT^{1∕2} = ρ^{(γa−1)∕2}; however, the nonlinear dependence is less pronounced than that for the nucleation rate. The normalised growth rate γ_{1} ∕γ_{0} is independent of the fraction of the surrounding gas composed of the accreting species f_{s}. Varying the spatial scale of the initial density perturbation does not affect the amplitude of the wave response.
When ω > N, internal gravity waves cannot propagate since the system cannot respond quick enough to the imposed driven perturbation and the waves are evanescent. When ω ≤ N, internal gravity waves can freely propagate, where the interplay between the driving frequency and the local buoyancy frequency results in greater response amplitudes for lower driving frequencies than for frequencies approaching the natural buoyancy frequency of the system (see Fig. 5). As the generated wave propagates away from the oscillation source, it encounters regions of differing background density and hence local buoyancy frequency. In response, the wave amplitude, speed, and wavelength change in harmony to conserve wave energy. For example, if the wave propagates into regions of lowerdensity (higherdensity), the amplitude of the wave increases (decreases), the wave speed increases (decreases), and the wavelength decreases (increases) in sympathy. If the wave encounters a region where its frequency is greater than the local buoyancy frequency, the wave ceases to propagate.
The numerical simulations presented in Fig. 5 are normalised to equilibrium reference values, aiding in their generalisation to other brown dwarf and exoplanet models. For example, the density response (top plot, Fig. 5) and the dust growth rate response (bottom plot, Fig. 5) are indicative of the typical response expected in atmospheric models beyond the exemplar presented. In contrast, the nucleation rate response is more complex and is driven by three main parameters: the ambient gas temperature T, the total gas pressure p, and the number density n_{x} of the nucleating species x. These interdependent parameters depend upon a number factors, including the chemical composition of the atmosphere and the chemical processes involved (eg. see Helling et al. 2017; Lee et al. 2015, 2018). To investigate the nucleation rate response, we explored the parameter space of the key atmospheric variables, in order to contextualise the results beyond the exemplar atmospheric model considered. Varying T_{eff} and log g explicitly instead can obfuscate the physical picture leading to the underlying cause of the nucleation rate variations and can be misleading when dealing with microphysical processes. Figure 6 shows the contour of the change in nucleation rate J_{*,1} (normalised to the total nucleation rate J_{*}) for a range of background temperatures T_{0} ∈ [500 K, 1200 K], pressures p_{0} ∈ [10^{−12} bar, 10^{−2} bar], and to cover a wide range of values from contemporary models (see Fig. 4, Lee et al. 2015; Fig. 4, Helling et al. 2008a; Fig. 2, Stark et al. 2013).
For the top plot of Fig. 6, we set p_{0} ≈ 4 × 10^{−5} bar (taken fromthe centre of the domain shown in Fig. 4) and vary and T_{0}. For the bottom plot, we held (see Fig. 2). We computed J_{*,1} using perturbed values of temperature T_{1} and pressure p_{1}. We obtained those values using the adiabatic equations of state, assuming a density perturbation resulting from the passage of an internal gravity wave in the most favourable scenario ρ_{1} ∕ρ_{0} = ±0.48 as shown in Fig. 5. We note that as a result of the wide exploration of parameter space, not all points in Fig. 6 correspond to a selfconsistently calculated pT atmosphericequilibrium state.
Figure 6 shows that the increase in nucleation rate J_{*,1} makes up a larger portion of the total nucleation rate at higher temperatures, lower background pressures (top plot), and lower (bottom plot). These conditions are unfavourable to efficient nucleation, resulting in very small values of J_{*,0}. However, the passage of an interval gravity wave can produce temporary, localised conditions allowing nucleation to occur at an enhanced rate. While that increase in nucleation may not be large in absolute terms, it is much larger than the background values, and leads toa large relative increase. This is consistent with results obtained by Helling et al. (2001) for simulated sound waves.
To elucidate this point further, we consider a slice of constant in the top plot of Fig. 6. As the temperature increases the equilibrium nucleation rate decreases which inhibits the growth of TiO_{2} clusters since nucleation favours cooler temperatures. As a result, the nucleation rate enhancement J_{*,1} due to the passage of an internal gravity wave relative to the equilibrium J_{*,0} is diminished leading to an increase in J_{*,1}∕J_{*}. The opposite occurs if the temperature decreases, leading to an increase in J_{*,1} ∕J_{*}. A similar response is evident when considering a slice of constant p_{0} in the bottom plot of Fig. 6.
Further to this, we consider a slice of constant background temperature in the top plot of Fig. 6. As increases, the equilibrium nucleation rate increases due to the increased supersaturation ratio S, since the gasphase TiO_{2} molecules are more likely to cluster and nucleate. As a result, J_{*,1} is diminished relative to the equilibrium J_{*,0}, leading to a decrease in J_{*,1}∕J_{*}. Similarly, consider a constant slice of temperature T in the bottom plot of Fig. 6: increasing the background pressure leads to a greater absolute TiO_{2} density , resulting in an increased background nucleation rate and lower relative increase due to the passage of an internal gravity wave J_{*,1} ∕J_{*}. If decreases the opposite case occurs, leading to an increase in J_{*,1}∕J_{*}.
To place the relative nucleation rate J_{*,1}∕J_{*} in a wider context, in the bottom plot of Fig. 6 we overplot the pT profiles for the brown dwarf model used for the numerical simulations and of typical and lowgravity L and T dwarfs (see Fig. 1 and Table 2). Furthermore, in the top plot of Fig. 6 we plot the line of constant temperature corresponding to an atmospheric pressure of p_{0} ≈ 10^{−5} bar for each of the model atmospheres considered. We note that each model profile has selfconsistently calculated nucleation rate as a function of atmospheric pressure p_{0}, temperature T_{0}, and atmospheric chemistry that may not correspond to a singular point on the parameter space contours. Therefore, these profiles give a helpful indication of the impact of internal gravity waves on the nucleation rate for the variety of substellar objects considered. For example, we consider the solid line (T_{eff} = 1500 K, logg = 5.0) in the top plot of Fig. 6. For , we can deduce J_{*,1}∕J_{*}≈ 0.95, equivalent to J_{*,1} ≈ 20J_{*,0}, which is consistent with the data shown in Fig. 2. The same result can be obtained using the solid line in the bottom plot, taking p_{0} ≈ 4 × 10^{−5}. The intersection of the constant pressure line and the model line yields J_{*,1}∕J_{*}≈ 0.95 (J_{*,1} ≈ 20J_{*,0}).
In general, Fig. 6 demonstrates that the impact of internal gravity waves on the nucleation rate is significant across the substellar objects considered. The profiles in Fig. 6 show that the strongest relative increases in nucleation are obtained when the conditions at equilibrium are less suited for efficient nucleation, which leads to any increase caused by the passage of a gravity wave to be comparatively large. This is visible on the top plot, where the warmer L dwarf models (LD: T_{eff}= 2000 K, logg = 4.5; LGLD: T_{eff}= 2000 K, logg = 3.0) are linked to strong increase in nucleation for a wider range of than the cooler T dwarf models (TD: T_{eff} = 1200 K, logg = 4.5; LGTD: T_{eff} = 1200 K, logg = 3.0). Similarly, the pT profiles on the bottom plots show that for a constant , the warmer L dwarwillf models exhibit a stronger relative increase in nucleation than T dwarf models.
Fig. 1
Atmospheric diagram (p, T; top panel), (p, y; middle panel), and (p, ρ; bottom panel) for the DRIFTPHOENIX brown dwarf atmosphere model used in numerical simulations (BD: T_{eff} = 1500 K, logg = 5.0, published by Stark et al. 2013 and RodríguezBarrera et al. 2018). Additional profiles for typical (LD, TD: log (g) = 4.5) and lowgravity (LGLD, LGTD: log (g) = 3.0) L dwarf (T_{eff} = 2000 K) and T dwarf (T_{eff} = 1200 K) models are shown for context (see also Table 2.) 
Fig. 2
Plot of the partial number density of TiO_{2} at equilibrium for the model used in simulations (logg = 5.0, T_{eff} = 1500 K). In the range of pressures used in the simulations (10^{−5} bar), is constant around 10^{−17}. 
Fig. 3
Buoyancy period T_{N} = 2π∕N for the DRIFTPHOENIX brown dwarf (logg = 5.0, T_{eff} = 1500 K) atmosphere model used in numerical simulations, as well as typical and lowgravity L and T dwarfs (see Fig. 1, Table 2). The range of periods for higher surface gravity levels (10^{1} –10^{2}s) is consistent with the atmospheric models and simulation results published by Freytag et al. (2010). In lowerdensity brown dwarfs, the buoyancy period rises to the order of 10^{3} s. 
Fig. 4
2D map of ρ_{1}∕ρ_{0} (top panel), J_{*,1}∕J_{*,0} (middle panel) and γ_{1}∕γ_{0} (bottom panel) after one oscillation period. For this example, waves are driven at ω = 0.25 N, and the amplitude of the driving perturbation ρ_{A} = 0.2 ρ_{0}. The vertical line shows the location of the slice used to measure the wave response, at a distance X(ω) from the centre of the domain so that X(ω) = X_{ref}ω_{ref}∕ω. 
Fig. 5
Plots of ρ_{1}∕ρ_{0} (top panel), J_{*,1}∕J_{*,0} (middle panel), and γ_{1}∕γ_{0} (bottom panel) after one period of an internal gravity wave, as a function of the amplitude of the density perturbation used to drive the waves. The nucleation plot shows a strong nonlinear increase in response to increased perturbation amplitude, more pronounced with low driving frequencies. The impact on dust growth is much weaker, and growth increases with the driving frequency. The measurements are taken as the maximum values along a vertical slice of the numerical domain, located at a distance X(ω) from the centre of the domain so that X(ω) = X_{ref}ω_{ref}∕ω. 
Fig. 6
Contour of the change in nucleation rate J_{*,1} (normalised to the total nucleation rate J_{*}), as a functionof equilibrium temperature T_{0}, (top panel) and equilibrium gas pressure p_{0} (bottom panel). The atmospheric profiles for the brown dwarf model used for simulations (BD: T_{eff} = 1500 K, logg = 5.0), as well as typical (LD,TD) and lowgravity (LGLD, LGTD) L and T dwarfs are plotted for context (see model details in Table 2; Fig. 1). The circle on both panels indicate the conditions present in the model used for simulations at the centre of the domain. 
5 Discussion
This paper has investigated and characterised the effect of linear internal gravity waves on the evolution of dust clouds in substellar (brown dwarf and gas giant exoplanetary) atmospheres for the first time. We have shown that in numerical fluid simulations, the passage of an internal gravity wave leads to an increase of dust nucleation by up to a factor 20, and an increase of dust mantle growth rate by up to a factor 1.6. Through an exploration of the wider substellar parameter space, we have shown that, in absolute terms, the increasein dust nucleation due to internal gravity waves is stronger in cooler (T dwarfs) and TiO_{2}rich substellar atmospheres. The relative increase, however, is greater in warm (L dwarf) and TiO_{2}poor atmospheresdue to conditions that are less suited for efficient nucleation at equilibrium. This latter point is important since the stronger thecontrast between the perturbed and equilibrium values, the better the chance of detecting an observable signal. Recent observations (Marocco et al. 2014) and models (Hiranaka et al. 2016) suggest that the extreme reddening of some L dwarfs could be due to a dust haze layer high up in their atmosphere, which could potentially be impacted by internal gravity waves.
The presence of the signature of an internal gravity wave in the spectra of a brown dwarf could indicate the presence of convection deep in the atmosphere or, in the case of a terrestrial exoplanetary atmosphere, that the body has a rocky, solid surface with relief, since such features are known to trigger the buoyancy oscillation required to generate the waves (Roeten et al. 2019). In such a scenario, the wavelength of the resulting wave could give an indication of the scale of the perturbing feature. Further investigation into the nonlinear evolution of internal gravity waves could potentially yield greater variations in atmospheric density and nucleation rate. Moreover introducing additional effects, such as the Coriolis effect and dynamical equilibria, and investigating their impact on the evolution of internal gravity waves and the resulting cloud cover, might yield further insight into inhomogeneous cloud coverage in substellar atmospheres.
Additionally, observations of the photometric variability resulting from the propagation of an internal gravity wave could provide a novel way of diagnosing the atmospheric gas density. We consider two identical, adjacent vertical atmospheric profiles. We assume that one profile contains a gas overdensity of amplitude ρ_{1}, and a corresponding dust overdensity ρ_{d1}, that occurs over a spatial length scale L_{0}, as the result of a propagating internal gravity wave. The ratio of the spectral flux density S from both can be expressed in terms of their respective optical depths, (49)
and we have assumed that the solid angle subtended by the features is the same. In a simple approach, we only consider absorption contributions from the gas and the dust in the atmosphere. To simplify matters further, we assume a total mean opacity to represent the contributions from the gas κ_{g} and from the dust κ_{d} respectively. The nonzero contributions of the optical depth integral over the extent of the atmosphere D, can be approximated to give, (51)
In the linear regime, Fig. 5 maps the normalised amplitude of the driving density perturbation, ρ_{A} ∕ρ_{0}, to the normalised amplitude of the resulting density, ρ_{1}∕ρ_{0}, and nucleation rate, J_{*1}∕J_{*0}, wave response,
χ and δ are functions of ρ_{A}∕ρ_{0}. Therefore, the ratio of flux densities can be expressed as, (54)
where is the mass of a dust particle, and a_{d} is the radiusof a dust grain. To relate the amplitude of the initial density perturbation to the equilibrium atmospheric density, we solve the differential equation for the buoyancy frequency (Eq. (23)) over the length scale of the perturbation L_{0}, giving (56)
where we have assumed that the N is approximately constant across L_{0}, and qN is the wave frequency with q ∈ (0, 1] (we assume q = 1 for simplicity). Therefore, rearranging Eq. (54), we obtain an expression for the equilibrium atmospheric density: (57)
This expression allows us to estimate the density of a substellar atmosphere based on potential observations of S_{0} ∕S_{1} on a timescale consistent with the buoyancy frequency. As an example, for the order of magnitude values list in Table 3, Eq. (58) gives an estimation of ρ_{0} ≈ 10^{−5} kg m^{−3}, which is consistent with contemporary atmospheric numerical models. This demonstrates that from observations of S_{0}∕S_{1}, L_{0}, and the timescale of variation, an estimation of the atmospheric density of a substellar atmosphere can be made. A more indepth analysis could involve calculating synthetic spectra showing the impact of the gravity wave that could be expected from observations, and will be considered in a further paper.
Acknowledgements
The authors are grateful to the anonymous referee for constructive comments and suggestions that have improved this paper. A.P. is grateful for funding and support received from Abertay University as part of the RLINCS studentship programme. C.R.S. is grateful for funding from the Royal Society via grant number RG160840 and from the Carnegie Trust for the Universities of Scotland via research incentive grant number RIG007788. E.K.H.L. acknowledges support from the University of Oxford and CSH Bern through the Bernoulli fellowship and support from the European community through the ERC advanced grant project EXOCONDENSE (PI: R.T. Pierrehumbert).
References
 Alvan, L., Brun, A. S., & Mathis, S. 2014, A&A, 565, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Biller, B. 2017, Astron. Rev., 13, 1 [Google Scholar]
 Buenzli, E., Apai, D., Radigan, J., Reid, I. N., & Flateau, D. 2014, ApJ, 782, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Freytag, B., Allard, F., Ludwig, H. G., Homeier, D., & Steffen, M. 2010, A&A, 513, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gail, H.P., Keller, R., & Sedlmayr, E. 1984, A&A, 133, 320 [NASA ADS] [Google Scholar]
 Garaud, P., Medrano, M., Brown, J. M., Mankovich, C., & Moore, K. 2015, ApJ, 808 [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Jehin, E., et al. 2013, A&A, 555, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hauschildt, P. H., & Baron, E. 1999, J. Comput. Appl. Math., 109, 41 [Google Scholar]
 Helling, C., & Casewell, S. 2014, A&Ar, 22, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Helling, C., & Woitke, P. 2006, A&A, 455, 325 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Helling, C., Oevermann, M., Lüttke, M. J. H., Klein, R., & Sedlmayr, E. 2001, A&A, 376, 194 [Google Scholar]
 Helling, C., Klein, R., Woitke, P., Nowak, U., & Sedlmayr, E. 2004, A&A, 423, 657 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Helling, C., Ackerman, A., Allard, F., et al. 2008a, MNRAS, 391, 1854 [Google Scholar]
 Helling, C., Woitke, P., & Thi, W. F. 2008b, A&A, 485, 547 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Helling, C., Tootill, D., Woitke, P., & Lee, G. 2017, A&A, 603, A123 [Google Scholar]
 Hiranaka, K., Cruz, K. L., Douglas, S. T., Marley, M. S., & Baldassare, V. F. 2016, ApJ, 830, 96 [Google Scholar]
 Lee, G., Helling, C., Giles, H., & Bromley, S. T. 2015, A&A, 575, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lee, G., DobbsDixon, I., Helling, C., Bognar, K., & Woitke, P. 2016, A&A, 594, A48 [Google Scholar]
 Lee, G. K. H., Blecic, J., & Helling, C. 2018, A&A, 614, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marley, M. S., & Robinson, T. D. 2015, ARA&A, 53, 279 [Google Scholar]
 Marley, M. S., Saumon, D., & Goldblatt, C. 2010, ApJ, 723, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Marocco, F., DayJones, A. C., Lucas, P. W., et al. 2014, MNRAS, 439, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Mittal, S. 2014, Int. J. High Perform. Comput. Netw., 7, 292 [Google Scholar]
 Robinson, T. D., & Marley, M. S. 2014, ApJ, 785, 158 [NASA ADS] [CrossRef] [Google Scholar]
 RodríguezBarrera, M. I., Helling, C., & Wood, K. 2018, A&A, 618, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roeten, K. J., Bougher, S. W., Benna, M., et al. 2019, J. Geophys. Res., 124, 3283 [CrossRef] [Google Scholar]
 Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327 [Google Scholar]
 Stark, C. R., & Diver, D. A. 2018, A&A, 611, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stark, C. R., Helling, C., Diver, D. A., & Rimmer, P. B. 2013, ApJ, 776, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Stark, C. R., Helling, C., & Diver, D. A. 2015, A&A, 579, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sutherland, B. R. 2010, Internal Gravity Waves (Cambridge: Cambridge Univerisity Press) [Google Scholar]
 Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJ, 804, L17 [CrossRef] [Google Scholar]
 Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJ, 817, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblin, P., Padioleau, T., Phillips, M. W., et al. 2019, ApJ, 876, 144 [Google Scholar]
 Vallis, G. K. 2017, Atmospheric and Oceanic Fluid Dynamics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Vetterling, W. T., Teukolsky, S. A., Press, W. H., & Flannery, B. P. 1992, Numerical recipes in C, 2nd edn. (Cambridge: Cambridge University Press) [Google Scholar]
 Vos, J. M., Allers, K., Apai, D., et al. 2019, ArXiv eprints [arXiv:1903.06691] [Google Scholar]
 Witte, S., Helling, C., & Hauschildt, P. 2009, A&A, 506, 1367 [Google Scholar]
 Witte, S., Helling, C., Barman, T., Heidrich, N., & Hauschildt, P. 2011, A&A, 529, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1
Atmospheric diagram (p, T; top panel), (p, y; middle panel), and (p, ρ; bottom panel) for the DRIFTPHOENIX brown dwarf atmosphere model used in numerical simulations (BD: T_{eff} = 1500 K, logg = 5.0, published by Stark et al. 2013 and RodríguezBarrera et al. 2018). Additional profiles for typical (LD, TD: log (g) = 4.5) and lowgravity (LGLD, LGTD: log (g) = 3.0) L dwarf (T_{eff} = 2000 K) and T dwarf (T_{eff} = 1200 K) models are shown for context (see also Table 2.) 

In the text 
Fig. 2
Plot of the partial number density of TiO_{2} at equilibrium for the model used in simulations (logg = 5.0, T_{eff} = 1500 K). In the range of pressures used in the simulations (10^{−5} bar), is constant around 10^{−17}. 

In the text 
Fig. 3
Buoyancy period T_{N} = 2π∕N for the DRIFTPHOENIX brown dwarf (logg = 5.0, T_{eff} = 1500 K) atmosphere model used in numerical simulations, as well as typical and lowgravity L and T dwarfs (see Fig. 1, Table 2). The range of periods for higher surface gravity levels (10^{1} –10^{2}s) is consistent with the atmospheric models and simulation results published by Freytag et al. (2010). In lowerdensity brown dwarfs, the buoyancy period rises to the order of 10^{3} s. 

In the text 
Fig. 4
2D map of ρ_{1}∕ρ_{0} (top panel), J_{*,1}∕J_{*,0} (middle panel) and γ_{1}∕γ_{0} (bottom panel) after one oscillation period. For this example, waves are driven at ω = 0.25 N, and the amplitude of the driving perturbation ρ_{A} = 0.2 ρ_{0}. The vertical line shows the location of the slice used to measure the wave response, at a distance X(ω) from the centre of the domain so that X(ω) = X_{ref}ω_{ref}∕ω. 

In the text 
Fig. 5
Plots of ρ_{1}∕ρ_{0} (top panel), J_{*,1}∕J_{*,0} (middle panel), and γ_{1}∕γ_{0} (bottom panel) after one period of an internal gravity wave, as a function of the amplitude of the density perturbation used to drive the waves. The nucleation plot shows a strong nonlinear increase in response to increased perturbation amplitude, more pronounced with low driving frequencies. The impact on dust growth is much weaker, and growth increases with the driving frequency. The measurements are taken as the maximum values along a vertical slice of the numerical domain, located at a distance X(ω) from the centre of the domain so that X(ω) = X_{ref}ω_{ref}∕ω. 

In the text 
Fig. 6
Contour of the change in nucleation rate J_{*,1} (normalised to the total nucleation rate J_{*}), as a functionof equilibrium temperature T_{0}, (top panel) and equilibrium gas pressure p_{0} (bottom panel). The atmospheric profiles for the brown dwarf model used for simulations (BD: T_{eff} = 1500 K, logg = 5.0), as well as typical (LD,TD) and lowgravity (LGLD, LGTD) L and T dwarfs are plotted for context (see model details in Table 2; Fig. 1). The circle on both panels indicate the conditions present in the model used for simulations at the centre of the domain. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.