© ESO 2020

1 Introduction

Brown dwarfs are low-mass, sub-stellar 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, fast-changing 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 cloud-free 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 main-sequence 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 sub-stellar 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 dust-hostile 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 sub-stellar 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 sub-stellar objects and its consequences for cloud variability. This paper presents a novel mechanism for potentially diagnosing the gas density of sub-stellar 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 Sub-stellar internal gravity waves

For a vertical slice of a sub-stellar 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 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₀ + Q₁. In the non-linear regime Q₀ ≫ Q₁, and powers of Q₁ 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 vorticity-stream 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∕ρ²) 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

where (19)

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, ρ₀ does not vary along x, therefore: (22)

where (23)

is the Brunt-Väisälä buoyancy frequency. In the case of a uniformly stratified atmosphere, assuming a solution of the form ψ ≈ exp[−i(ωt + k_xx + k_yy)] 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 Brunt-Vä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⁻³ s⁻¹], 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:

where,

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₀ 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₂ for this paper).

For this paper, we assumed a constant value for the surface tension of TiO₂, (see Helling et al. 2001; Lee et al. 2015), and a density ρ_m ≈4230 kg m⁻³ for TiO₂. In the context of an internal gravity wave propagating through a sub-stellar 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 gas-phase 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² 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⁻¹]. 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 gas-phase 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 gas-phase species, with varying relative abundances in the gas-phase, participating in surface chemistry leading to mantle growth.

3 Numerical simulations

3.1 Model equations

The linearised governing equations can be cast in non-dimensional form, defining

where L, T, and M are characteristic values of length, time, and mass respectively. Therefore, in non-dimensional form, Eqs. (16)–(18) become

Casting the model equations in non-dimensional 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 Over-Relaxation (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 first-order derivatives were interpolated using third-order 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 Sub-stellar atmospheric model

The aim of this paper was to investigate the effects of internal gravity waves on the evolution of dust clouds in sub-stellar 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 sub-stellar 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 low-gravity L and T dwarfs for context, to place the simulations in a wider range of brown dwarfs examples. The atmosphericdata was generated by the DRIFT-PHOENIX 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 one-dimensional; 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⁴ m). The numerical domain was defined by the geometric extension y ∈ [25 km, 35 km]; the pressure p₀ ∈ [3 × 10⁻⁵ bar, 3 × 10⁻³ bar]; the temperature T₀ ∈ [700 K, 850 K]; and the density ρ₀ ∈ [10⁻⁷ kg m⁻³, 10⁻⁵ kg m⁻³]. For these conditions, (see Fig. 2). The perturbation length scale was of the order of 10³ m.

4 Results

In a stratified sub-stellar 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 Brunt-Väisälä buoyancy frequency. Figure 3 shows the buoyancy period T_N = 2π∕N as a function of pressure p for the sub-stellar 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 low-gravity 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 length-scale 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 large-scale turbulent motions, the frequency of the wave is set by thefrequency of the source for frequencies below the buoyancy frequency. In sub-stellar 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² (Marley & Robinson 2015; Tremblin et al. 2019)– for example, l_conv ≈ 10³ km in Freytag et al. (2010) – varying with timescales of the order of 10³ 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ρ₀. 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 ρ = ∑_sm_sn_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₀ such that S(T₀) ≫ 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 over-densities in the atmosphere as a result of an internal gravity wave can decrease the local nucleation rate; whereas, local under-densities 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 ρ₁, J_*,1, and γ₁ 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 non-linearly with the perturbation amplitude by up to a factor 20 in the most favourable scenario (ω = 0.25N, ρ_A = 0.2ρ₀). This is a result of the complex non-linear dependence of the nucleation rate (Eq. (26)), and not as a consequence of a non-linear 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 non-linearly, the corresponding nucleation rate could exhibit an enhanced non-linear 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 non-linear function of the perturbed density, γ ∝ ρT^1∕2 = ρ^{(γ_a−1)∕2}; however, the non-linear dependence is less pronounced than that for the nucleation rate. The normalised growth rate γ₁ ∕γ₀ 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 lower-density (higher-density), 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₀ ∈ [500 K, 1200 K], pressures p₀ ∈ [10⁻¹² bar, 10⁻² 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₀ ≈ 4 × 10⁻⁵ bar (taken fromthe centre of the domain shown in Fig. 4) and vary and T₀. For the bottom plot, we held (see Fig. 2). We computed J_*,1 using perturbed values of temperature T₁ and pressure p₁. 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 ρ₁ ∕ρ₀ = ±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 self-consistently 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₂ 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₀ 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 gas-phase TiO₂ 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₂ 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 low-gravity 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₀ ≈ 10⁻⁵ bar for each of the model atmospheres considered. We note that each model profile has self-consistently calculated nucleation rate as a function of atmospheric pressure p₀, temperature T₀, 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 sub-stellar 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₀ ≈ 4 × 10⁻⁵. 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 sub-stellar 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 DRIFT-PHOENIX brown dwarf atmosphere model used in numerical simulations (BD: Teff = 1500 K, logg = 5.0, published by Stark et al. 2013 and Rodríguez-Barrera et al. 2018). Additional profiles for typical (LD, TD: log (g) = 4.5) and low-gravity (LGLD, LGTD: log (g) = 3.0) L dwarf (Teff = 2000 K) and T dwarf (Teff = 1200 K) models are shown for context (see also Table 2.)

Fig. 2 Plot of the partial number density of TiO2 at equilibrium for the model used in simulations (logg = 5.0, Teff = 1500 K). In the range of pressures used in the simulations (10−5 bar), is constant around 10−17.

Fig. 3 Buoyancy period TN = 2π∕N for the DRIFT-PHOENIX brown dwarf (logg = 5.0, Teff = 1500 K) atmosphere model used in numerical simulations, as well as typical and low-gravity L and T dwarfs (see Fig. 1, Table 2). The range of periods for higher surface gravity levels (101 –102s) is consistent with the atmospheric models and simulation results published by Freytag et al. (2010). In lower-density brown dwarfs, the buoyancy period rises to the order of 103 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(ω) = Xrefω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 non-linear 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(ω) = Xrefωref∕ω.

Fig. 6 Contour of the change in nucleation rate J*,1 (normalised to the total nucleation rate J*), as a functionof equilibrium temperature T0, (top panel) and equilibrium gas pressure p0 (bottom panel). The atmospheric profiles for the brown dwarf model used for simulations (BD: Teff = 1500 K, logg = 5.0), as well as typical (LD,TD) and low-gravity (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 sub-stellar (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 sub-stellar 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₂-rich sub-stellar atmospheres. The relative increase, however, is greater in warm (L dwarf) and TiO₂-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 non-linear 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 sub-stellar 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 over-density of amplitude ρ₁, and a corresponding dust over-density ρ_d1, that occurs over a spatial length scale L₀, 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)

where (50)

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 non-zero 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 ∕ρ₀, to the normalised amplitude of the resulting density, ρ₁∕ρ₀, and nucleation rate, J_*1∕J_*0, wave response,

χ and δ are functions of ρ_A∕ρ₀. Therefore, the ratio of flux densities can be expressed as, (54)

where (55)

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₀, giving (56)

where we have assumed that the N is approximately constant across L₀, 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)

where (58)

This expression allows us to estimate the density of a sub-stellar atmosphere based on potential observations of S₀ ∕S₁ 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 ρ₀ ≈ 10⁻⁵ kg m⁻³, which is consistent with contemporary atmospheric numerical models. This demonstrates that from observations of S₀∕S₁, L₀, and the timescale of variation, an estimation of the atmospheric density of a sub-stellar atmosphere can be made. A more in-depth 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

All Tables

All Figures

Fig. 1 Atmospheric diagram (p, T; top panel), (p, y; middle panel), and (p, ρ; bottom panel) for the DRIFT-PHOENIX brown dwarf atmosphere model used in numerical simulations (BD: Teff = 1500 K, logg = 5.0, published by Stark et al. 2013 and Rodríguez-Barrera et al. 2018). Additional profiles for typical (LD, TD: log (g) = 4.5) and low-gravity (LGLD, LGTD: log (g) = 3.0) L dwarf (Teff = 2000 K) and T dwarf (Teff = 1200 K) models are shown for context (see also Table 2.)

In the text

	Fig. 2 Plot of the partial number density of TiO2 at equilibrium for the model used in simulations (logg = 5.0, Teff = 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 TN = 2π∕N for the DRIFT-PHOENIX brown dwarf (logg = 5.0, Teff = 1500 K) atmosphere model used in numerical simulations, as well as typical and low-gravity L and T dwarfs (see Fig. 1, Table 2). The range of periods for higher surface gravity levels (101 –102s) is consistent with the atmospheric models and simulation results published by Freytag et al. (2010). In lower-density brown dwarfs, the buoyancy period rises to the order of 103 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(ω) = Xrefω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 non-linear 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(ω) = Xrefω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 T0, (top panel) and equilibrium gas pressure p0 (bottom panel). The atmospheric profiles for the brown dwarf model used for simulations (BD: Teff = 1500 K, logg = 5.0), as well as typical (LD,TD) and low-gravity (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