Onset of Nonlinear Internal Gravity Waves in Intermediate-Mass Stars

TL;DR

This paper investigates how internal gravity waves in stars' radiative zones become nonlinear, influenced by stellar parameters, and explores their potential role in angular momentum transport and mixing.

Contribution

It introduces a criterion to determine when IGW become nonlinear in non-rotating stars, considering effects of radiative diffusion and density stratification based on stellar models.

Findings

Nonlinear wave energies increase with stellar mass.

Generation spectrum, convective velocities, and stratification influence wave nonlinearity.

Trends vary with metallicity and age depending on the wave generation spectrum.

Abstract

Internal gravity waves (IGW) propagate in the radiation zones of all stars. During propagation, their amplitudes are affected by two main features: radiative diffusion and density stratification. We have studied the implications of these two features on waves traveling within the radiative zones of non-rotating stars with stellar parameters obtained from the one dimensional stellar evolution code, MESA. As a simple measure of induced wave dynamics, we define a criterion to see if waves can become nonlinear and if so, under what conditions. This was done to understand the role IGW may play in angular momentum transport and mixing within stellar interiors. We find that the IGW generation spectrum, convective velocities and the strength of density stratification all play major roles in whether waves become nonlinear. With increasing stellar mass, there is an increasing trend in nonlinear wave energies. The trends with different metallicities and ages depend on the generation spectrum.

Tables

Table 1. Table 1: The table shows how m 𝑚 m and n 𝑛 n in Eq. ( 5 ), are defined for different works. The flat spectrum represents one with no dependence on ω 𝜔 \omega or ℓ ℓ \ell . The letters shown in square brackets will be used to represent its respective spectrum.

Spectra	m	n
Flat[F]	0	0
Kumar et al. (1999)[K]	-2.17	1
Lecoanet & Quataert (2013)[LD]	-4.25	2.5
Rogers et al. (2013)[R]	-0.6	-0.9

Table 2. Table 2: The table shows percentages of the radiative zone (in terms of the total stellar radius) that were cut out close to the inner convective zone (Inner Cut-out) and close to the surface convection zone (Outer Cut-out) for the different stellar models we used.

Stellar Mass	Stellar Age	Inner Cut-out	Outer Cut-out
3	ZAMS	0.048	0.0
	midMS	0.020	0.001
	TAMS	0.0	0.016
7	ZAMS	0.103	0.007
	midMS	0.059	0.031
	TAMS	0.011	0.311
20	ZAMS	0.394	0.004
	midMS	0.117	0.035

Table 3. Table 3: The table shows the energy ratios of nonlinear waves for different stellar metallicities at ZAMS.

Spectra	$u_c/u_{\mathrm{MLT}}$	$Z$ = 0.06	$Z$ = 0.02	$Z$ = 0.006
F	1	0.968	0.968	0.959
	3	0.971	0.9705	0.9635
R	1	0.0	0.0	$5.54\times {10}^{-2}$
	3	0.234	0.223	0.182
K	1	0.0	0.0	0.0
	3	0.481	0.480	0.606
LD	1	0.0	0.0	0.0
	3	0.0	0.0	0.0

Equations

$$F\sim u_c^3\frac{{\rho }_c{k}_h^3}{r_c^2{N}_c}{\left(\frac{u_c}{{\omega }_c}\right)}^5{\left(\frac{\omega }{{\omega }_c}\right)}^{-10/3}\exp \left[-k_h^2{\left(\frac{u_c}{{\omega }_c}\right)}^2{\left(\frac{\omega }{{\omega }_c}\right)}^{-4/3}\right],$$

$$k_h=\frac{\sqrt{\ell (\ell +1)}}{r}.$$

$$u_{h,K}\sim u_c\frac{k_h^2}{r_c\sqrt{\ell (\ell +1)}}{\left(\frac{u_c}{{\omega }_c}\right)}^3{\left(\frac{\omega }{{\omega }_c}\right)}^{-13/6}$$

$$\exp \left[-\frac{k_h^2}{2}{\left(\frac{u_c}{{\omega }_c}\right)}^2{\left(\frac{\omega }{{\omega }_c}\right)}^{-4/3}\right],$$

$$u_{h,LD}\sim u_c{\left(\frac{k_h{u}_c}{{\omega }_c}\right)}^{5/2}{\left(\frac{\omega }{{\omega }_c}\right)}^{-17/4},$$

$$u_0(\omega ,\ell )\propto {\omega }^m{\left(\sqrt{\ell (\ell +1)}\right)}^n$$

$$\tau (\omega ,\ell ,r)={\int }_{r_{\text{interface}}}^r\frac{\gamma [\omega ,\ell ,r^{\prime }]}{\left|v_g[\omega ,\ell ,r^{\prime }]\right|}d{r}^{\prime }$$

$$\gamma (\omega ,\ell ,r)=K{k}_v^2$$

$$k_v^2=k_h^2\frac{N^2-{\omega }^2}{{\omega }^2},$$

$$N^2=\frac{g^2\rho }{P}\frac{{\chi }_T}{{\chi }_{\rho }}({\nabla }_{\mathrm{ad}}-{\nabla }_T+B),$$

$$v_g(r)=\frac{\partial \omega }{\partial k_v}=-\frac{(N^2-{\omega }^2{)}^{1/2}{\omega }^2}{k_h{N}^2}.$$

$$K(r)=\frac{16\sigma T^3}{3{\rho }^2\kappa c_p},$$

$$\frac{\gamma }{v_g}=\frac{16\sigma T^3}{3{\rho }^2\kappa c_p}\left(\frac{{\left(\ell (\ell +1)\right)}^{3/2}{N}^3}{r^3{\omega }^4}\right){\left(1-\frac{{\omega }^2}{N^2}\right)}^{1/2}.$$

$$u_h(\omega ,\ell ,r)=u_0(\omega ,\ell ){\left(\frac{r_0}{r}\right)}^{3/2}{\left(\frac{{\rho }_c}{\rho (r)}\right)}^{1/2}{\left(\frac{N^2-{\omega }^2}{N_c^2-{\omega }^2}\right)}^{1/4}{e}^{-\tau /2},$$

$$ϵ={\xi }_h{k}_h=\frac{u_h}{\omega }{k}_h,$$

$$\frac{k_h}{k_v}=\sqrt{\frac{{\omega }^2}{N^2-{\omega }^2}},$$

$$k_h{\xi }_h\sim \left(\sqrt{\frac{{\omega }^2}{N^2-{\omega }^2}}\right)k_v\left(\sqrt{\frac{N^2-{\omega }^2}{{\omega }^2}}\right){\xi }_v\sim k_v{\xi }_v,$$

$$N(r{)}^2>0$$

$$N(r{)}^2<0$$

$$u_c=u_{\mathrm{MLT}}=\frac{1}{{\mathrm{N}}_g/2}\sum _{{\mathrm{N}}_g/4}^{3{\mathrm{N}}_g/4}{u}_i$$

$${\omega }_c=\frac{u_c}{r_c}.$$

$$E\propto u_h^2$$

$$E$$

$$E\propto {\left(\sqrt{\ell (\ell +1)}\right)}^4{\omega }^{-8.50}.$$

$$E\propto {\left(\sqrt{\ell (\ell +1)}\right)}^2{\omega }^{-4.34}\exp \left[-k_h^2{\left(\frac{u_c}{{\omega }_c}\right)}^2{\left(\frac{\omega }{{\omega }_c}\right)}^{-4/3}\right].$$

$$\frac{E_{ϵ⩾0.1}}{E_{\mathrm{Total}}}=\frac{{\left({\sum }_{l,\omega }{u}_h^2\Delta \omega \right)}_{ϵ⩾0.1}}{{\sum }_{l,\omega }{u}_h^2\Delta \omega },$$

Full text

††pubyear: 2018††pagerange: Onset of Nonlinear Internal Gravity Waves in Intermediate-Mass Stars–Onset of Nonlinear Internal Gravity Waves in Intermediate-Mass Stars

1 Introduction

Internal gravity wave propagation in stably stratified environments is a widely studied subject as these waves occur naturally in Earth’s atmosphere and oceans (Staquet & Sommeria, 2002; Fritts & Alexander, 2003; Sutherland, 2010), and are an important factor in the transport of material, angular momentum and energy. For instance, atmospheric gravity waves are formed when a uniform layer of air encounters an obstacle such as a mountain, in a stratified atmospheric region. These waves can then carry angular momentum and energy upwards towards the ionosphere, leading to the formation of wave clouds. Internal gravity waves also impose themselves on ocean surfaces, leading to formation of periodic calm and rough patches which are visible in satellite images.

In stars, IGW propagate only in radiation zones and as they propagate, some of them are reflected near stellar surfaces, leading to constructive interference and the formation of standing waves. The detection of these g-modes has been one of the central aims of asteroseismology (Aerts et al., 2010). All propagating gravity waves undergo radiative damping (Schatzman, 1993; Zahn et al., 1997). Some of them eventually break (Barker & Ogilvie, 2010) or encounter critical layers. Damping occurs due to radiative diffusion of photons across the IGW wavelength where perturbation temperatures have opposite signs. Wave breaking, on the other hand, is a non-linear process that occurs when the wave amplitude reaches a critical value, which causes the energy in waves to be converted to turbulent kinetic energy. These processes lead to energy, angular momentum transport and chemical mixing within the radiation zone (Press, 1981; Garcia Lopez & Spruit, 1991; Montalban, 1994; Charbonnel & Talon, 2005; Rogers & McElwaine, 2017). In most cases, IGWs are excited stochastically by turbulent convection within the stellar convection zone, leading to a predominantly inward (towards the center) propagation (Brummell et al., 2002; Dintrans et al., 2003) for stars with convective envelopes ($M<1.5\,M_{\odot}$) and predominantly outward propagation (towards the surface) for stars with convective cores ($M>1.5\,M_{\odot}$).

Internal gravity waves have been proposed to solve a host of unexplained observations in stars. Press (1981) first invoked IGW to explain what was then the neutrino problem. Garcia Lopez & Spruit (1991) and a series of papers by Montalban and Schatzman (Montalban, 1994; Montalban & Schatzman, 1996; Montalbán & Schatzman, 2000a) suggested IGW could produce the enhanced mixing (Press & Rybicki, 1981) needed to explain the observed lithium depletion in F stars. Zahn et al. (1997) and later, Kumar et al. (1999) suggested that IGW could explain the uniform rotation of the solar radiative interior. Denissenkov & Tout (2003) showed that turbulent mixing induced by IGWs can mix protons into helium and carbon rich zones causing the formation of carbon-12 pockets in low-mass asymptotic giant branch (AGB) stars. Rogers et al. (2012) suggested IGW could change the surface rotation rate of intermediate-mass stars, thus explaining the discrepancy of obliquity between hot Jupiters around hot and cool stars. Aerts et al. (2009) and Aerts & Rogers (2015) suggested IGW could explain the observed macroturbulence seen across the upper main sequence. Fuller et al. (2014) suggested IGW could cause the enhanced angular momentum transport needed to explain the observed differential rotation in evolved stars and may be responsible for pre-supernovae mass ejection (Fuller et al., 2015). Although IGW may be responsible for a variety of effects in stellar and planetary interiors, the details of their influence depends sensitively on the spectrum of waves generated by turbulent convection, which is still uncertain.

Neglecting the details of wave generation and assuming that wave propagation is linear and no critical layers form, wave amplitudes depend only on density stratification and thermal diffusivity. While very few theories have considered the possibility of wave breaking, it may be possible in some circumstances (Barker & Ogilvie, 2010; Rogers et al., 2013). In stars with convective cores, IGWs propagate outward, into a region of decreasing density, causing the wave amplitudes to increase (in the absence of thermal diffusion), which ultimately means these waves are likely to be more dynamically relevant in intermediate-mass stars than lower-mass stars. Hence, we focus our study here on intermediate-mass stars. While numerical simulations, such as those carried out in Rogers et al. (2013) can simulate the convective-radiative interface and the self-consistent generation of waves, they do so at artificially high viscosities, high thermal diffusivities and/or reduced dimensionality. Therefore, it is currently impossible to self-consistently simulate the generation and propagation of IGWs within a realistic stellar interior.

In this paper, we take a different approach. Assuming a given generation spectrum and that the waves propagate linearly through the radiation zones of stars, we determine what the likelihood is that these waves could become nonlinear. We use realistic thermal diffusivities and other stellar properties generated by the one-dimensional stellar evolution code, Modules for Experiments in Stellar Astrophysics (MESA, version 8845; Paxton et al. (2011)). To investigate how the wave amplitude changes with radius in a star, we consider two main effects: amplitude enhancement due to density stratification and amplitude reduction due to radiative diffusion. We then calculate the nonlinearity parameter as a function of wavenumber, frequency and initial amplitude to determine the likelihood of wave breaking. The work done in Shiode et al. (2013), which predicted the observability of convectively-driven g-modes under the influence of radiative damping in stars with convective cores, is similar, but here, we focus on the possibility of gravity waves becoming nonlinear. While this approach investigates the amplitude evolution of IGWs, it does not look into the influence of these breaking waves on angular momentum transport and chemical mixing.

This paper is organised as follows. In section 2, we discuss IGW physics such as generation by turbulent convection, wave propagation in the linear regime and IGW nonlinearity. In section 3, we discuss the stellar models generated by MESA and the initial parameters used to set up these models. In section 4, we present the analysis done on wave propagation and nonlinearity for varying masses and ages. Section 5 discusses these results and concludes this work.

2 Internal Gravity Wave Phenomenon

2.1 IGW Generation by Turbulent Convection

The classic prescription for wave generation by convection goes back to Goldreich & Kumar (1990), which was originally developed for pressure waves (p-modes). This prescription was then extended to gravity waves by Kumar et al. (1999). The original prescription treats wave generation by bulk convection due to three sources: changes in entropy at fixed pressure (monopole terms), buoyancy variations (dipole terms) and internal stresses (quadrupole terms). In that work, it was argued that while the monopole and dipole terms generate more acoustic radiation, they very nearly cancel each other and therefore, the only term considered was the quadrupole term, or the term due to internal Reynolds stresses. This assumption was carried forward in Kumar et al. (1999), yet it is unclear that this approximation is valid for IGW. Although there is evidence that treating IGW generation as an inhomogeneous wave equation with a source term is appropriate (Lecoanet & Quataert, 2013), it is still unclear what this source term should be or whether the theoretical picture is an accurate one.

The most obvious deficiency of this model is that it does not account for wave generation by overshooting plumes (Rieutord & Zahn, 1995) while it is well known from experiments (Ansong & Sutherland, 2010) and simulations that plume generation is relevant and possibly dominant. In terms of theoretical work, this was first studied by Townsend (1966) in a terrestrial context. Montalbán & Schatzman (2000b) used a plume model developed by Rieutord & Zahn (1995) to derive the IGW spectrum generated at the bottom of the convection zone of a solar-type star. Later, Pinçon et al. (2016) developed a semi analytical model for the same situation to be used in 1D stellar evolution codes. A common result of all these models is an exponential frequency spectrum, proportional to $\exp\left(-\omega^{2}t_{\mathrm{b}}^{2}\right)$ with a characteristic plume incursion timescale $t_{\mathrm{b}}$ and $\omega$ being the generated wave frequency. The wave number dependence is similarly exponential. Then, the total spectrum results from a combination of plumes of different sizes and timescales. Although we believe this generation spectra to be important, we do not include it here but note that the numerical generation spectra from Rogers et al. (2013) is likely to be similar.

To date, no numerical simulations have been able to reproduce any of the theoretically predicted wave spectra (Rogers et al., 2013; Alvan et al., 2014; Alvan et al., 2015) and similarly, neither have any laboratory experiments (Ansong & Sutherland, 2010). However, a recent work by Couston et al. (2018) has confirmed one of theoretical generation spectra investigated in this work through numerical simulations involving the Boussinesq approximation. Also, there is new evidence from 3D simulations that numerical solutions are similar to plume spectra for a variety of plume length and timescales (Edelmann at al (submitted)). Whether this discrepancy between simulations/experiments and theory is due to inadequate assumptions of the theory or unrealistic parameters of the simulations, or both, is still not known. Therefore, the spectrum of waves generated by convection remains a fundamental problem in all theories of wave transport (of angular momentum, energy or species). The first wave generation spectrum explored in this work is that from Kumar et al. (1999). Using flux equations calculated for p-modes from previous work (Goldreich & Kumar, 1990; Goldreich et al., 1994), Kumar et al. (1999) assumed continuity of radial velocities across the convective-radiative interface and derived a flux equation for IGW generated by Reynolds stresses as

[TABLE]

where $\ell$ is the spherical harmonic degree (we will be referring to $\ell$ as wave number, unless stated otherwise) while $u_{c}$, $\omega_{c}$, and $\rho_{c}$ are the bulk convective velocity, convective turnover frequency (in units of Hz) and density at the convective-radiative boundary, respectively. The term, $N_{c}$ is the Brunt–Väisälä frequency at the bottom of the radiation zone while $r_{c}$ is the thickness of the convection zone. The horizontal wavenumber, $k_{h}$ is defined as

[TABLE]

From here on, the subscript $h$ will represent horizontal quantities while $v$ will represent vertical quantities. We convert Eq. (1) into a displacement equation by equating it to $\rho_{c}\left(\omega\xi_{h}\sqrt{\ell(\ell+1)}\right)^{2}v_{g}$ as given in Kumar et al. (1999), where $v_{g}$ is the vertical IGW group velocity (see Eq. (10)), $\xi_{h}$ is the horizontal wave displacement and we have assumed that $\xi_{h}\gg\xi_{v}$. Using $u_{h}=\xi_{h}\omega$ gives the following relation,

[TABLE]

where $u_{h}$ is the horizontal fluid velocity and the subscript K refers to “Kumar". The work by Kumar et al. (1999) was revisited by Lecoanet & Quataert (2013), who also produced theoretical predictions for the IGW generation flux, but for three different stratification profiles at the convective-radiative transition. All three profiles have various power laws on $k_{h}$ and $\omega$, which are determined by the Brunt–Väisälä frequency profile across the convective-radiative interface (i.e. discontinuous, smooth continuous or piecewise linear). Assuming $N\gg\omega$, the horizontal velocity perturbation can be written as

[TABLE]

where the subscript LD refers to “Lecoanet Discontinuous". We only work with the discontinuous case here, as this scenario is currently considered the most probable in stars (Lecoanet, private communication).

In addition to theoretical work, there has been significant numerical work done to investigate the transition between convective and radiative regions and convective overshoot (Hurlburt et al., 1986; Brummell et al., 2002) with a few having focused on IGW generation (Rogers & Glatzmaier, 2005; Brun et al., 2011; Rogers et al., 2013; Alvan et al., 2014; Alvan et al., 2015). For instance, Rogers et al. (2013) predicts that for a star with zero rotation, the energy carried111Approximately 79% of the energy in IGWs were fit by a separable function, where E($\omega,l$) $\propto$ g($\omega$)*f($\ell$). by IGW scales as $E\propto\omega^{-1.2}\left(\sqrt{\ell(\ell+1)}\right)^{-1.8}$ in the low wave number ($\ell\lesssim 10$) and low frequency ($\omega\lesssim 10\;\mu$Hz) regime, just outside the convection zone. While it is possible that the negative exponent associated with the wavenumber dependence is due to the reduced dimensionality of those simulations, more recent three-dimensional simulations also show reduced efficiency for shorter wavelengths/larger wavenumbers (Edelmann at al (submitted), Alvan et al. (2014)). Furthemore, these simulations show a similar frequency dependence, so we take the 2D simulations to be representative. Using $E\propto\left(\omega\xi_{h}\right)^{2}\propto u_{h,R}^{2}$ where $\xi$ is the fluid displacement amplitude and the subscript R refers to “Rogers", $u_{h,R}$ is found to be proportional to $\omega^{-0.6}\left(\sqrt{\ell(\ell+1)}\right)^{-0.9}$.

The IGW generation spectrum, which represents the collection of waves generated at the convective-radiative interface can be summarised as power laws in both wave number and frequency and described by:

[TABLE]

where the values of $m$ and $n$ for different prescriptions are shown in Table 1. An additional IGW generation spectrum, the flat spectrum, has been introduced in Table 1, which represents a spectrum with no dependence on $\ell$ and $\omega$. This spectrum will act as a test case and will be hereon referred to as spectrum F. Also, the table does not include the extra dependences of spectrum K on $k_{h}$ and $\omega$ in the exponential term. However, we do consider the effect of this exponential term in our analysis.

2.2 Radiative Damping and Density Stratification

As IGW generated at the convective-radiative interface propagate outward, their amplitudes are affected by two main features; radiative damping and density stratification. We follow the work of Kumar et al. (1999), which built on the work of Press (1981), to investigate the effect of radiative damping on wave amplitudes. The damping opacity, $\tau$ is defined as

[TABLE]

where $\gamma$ is the damping rate, given by

[TABLE]

and $k_{v}$, the vertical wavenumber, is defined as

[TABLE]

where $N$ is the Brunt–Väisälä frequency. In MESA, $N$ is defined using the Ledoux criterion as

[TABLE]

where $g$ is the gravitational acceleration, $\chi_{T}$ is $(\partial\ln P/\partial\ln T)_{\rho}$, $\chi_{\rho}$ is $(\partial\ln P/\partial\ln\rho)_{T}$, $\nabla_{\mathrm{ad}}$ is the adiabatic temperature gradient, $\nabla_{T}$ is the true temperature gradient and $B$ is the compositional gradient. In the case of ideal gas, ${\chi_{T}}/{\chi_{\rho}}=1$. Note that it is common to work in the asymptotic regime and assume $N\gg\omega$. However, we keep the $\sqrt{N^{2}/\omega^{2}-1}$ term in our propagation equation, mainly because although the asymptotic assumption applies in most of the radiation zone, it fails to do so close to the convective-radiative boundaries as $N$ approaches 0. The vertical group velocity, $v_{g}$ is defined as

[TABLE]

The thermal diffusivity, $K$, is given by

[TABLE]

where $\sigma$ is the Stefan-Boltzmann constant while $T$, $\rho$, $\kappa$ and $c_{p}$ are stellar temperature, density, opacity and specific heat capacity, respectively, which are all functions of stellar radius. The full expression of terms inside the integral of Eq. (6) then, becomes

[TABLE]

Using stellar parameters and Brunt–Väisälä frequencies obtained directly from MESA, $\gamma(r)/v_{g}(r)$ was calculated and integrated using the cumulative trapezoidal method. Values of $\tau$ (see Eq. (6)) for different propagation distances, ($r-r_{\text{interface}}$), were then inserted into the following equation to obtain wave amplitude as a function of radius:

[TABLE]

where $u_{0}$ is the velocity perturbation amplitude (which takes the magnitude of convective velocities, $u_{c}$, into account) from different generation spectra described in Section 2.1 and $r_{0}$ is the radius at which the waves are generated. This expression, excluding the damping term, is within the WKB approximation. It can be derived from the anelastic, linearised continuity, momentum and energy equations, as done in Press (1981) and Zahn et al. (1997), with the assumptions of adiabatic perturbations and zero source terms in a medium with density stratification.

2.3 Nonlinearity Parameter

Phillips (1966), Press (1981) and Barker & Ogilvie (2010) have shown that the nonlinearity of waves can be represented by the ratio of wave displacement to wavelength in either horizontal or vertical directions. This leads to the definition of the nonlinearity parameter, $\epsilon$, as follows:

[TABLE]

where we have replaced the horizontal displacement with horizontal perturbation velocity over frequency. An IGW is said to be nonlinear, when its displacement becomes comparable to or larger than its wavelength, which is equivalent to the strong condition, $\epsilon\geqslant 1$. However, one also expects that as $\epsilon$ approaches 1, nonlinear effects can become important.

To consider only the horizontal components, we can first show that

[TABLE]

from rearranging Eq. (8). In the Boussinesq approximation which applies locally, the divergence of the velocity perturbations is zero and thus, $k_{h}u_{h}=k_{v}u_{v}$. So, the following relation is expected to be satisfied in the radiation zone:

[TABLE]

which shows that considering wave nonlinearity in the horizontal and vertical directions leads to the same result.

3 Model Setup

Models of stars with masses ranging from 3 $M_{\odot}$ to 20 $M_{\odot}$ were constructed using the MESA stellar evolution code from ZAMS, when the core hydrogen mass fraction ($X_{c}$) is approximately 0.7 to middle of main sequence (midMS), when $X_{c}\approx$ 0.35, followed by terminal age main sequence (TAMS), when $X_{c}\approx$ 0.01. Mass loss was implemented through a fixed stellar wind scheme. Stellar metallicity, $Z$, was set to be equal to the solar value of $Z$ = 0.02, although we do test the effect of varying metallicities in Section 4.4. The mixing length parameter was set 1.8. The sign of the Brunt–Väisälä frequency was used to determine the extent of the radiation zone which means

[TABLE]

refers to a radiation zone and

[TABLE]

refers to a convection zone. However, to account for the fact that the IGW frequencies must be less than the Brunt–Väisälä frequency, the extent of the radiation zone was further limited to regions where the Brunt–Väisälä frequency exceeds 100 $\mu$Hz, which was the highest IGW wave frequency we tested. It was found that in all our stellar models, the Brunt–Väisälä frequency increases rapidly from zero at the core convective-radiative interface and decreases rapidly back to zero at the surface radiative-convective interface. Therefore, at both the top and bottom convective-radiative boundaries, the regions of the radiative zone cut out by imposing the above condition were small, being less than 0.5% of the star’s total radius (see Table 2). Thermal diffusivities (see Eq. (11)) of 8 stellar models were calculated using opacities from MESA as shown in Fig. 1. The 20 $M_{\odot}$ TAMS model has not been included as it was found to develop an intermediate convection zone. It can be seen that for any particular stellar mass, the convection zone shrinks as the star evolves while the radiation zone expands. This is an important factor in IGW propagation which will be addressed further in the Section 4.5.

Figure 2 is an example of how wave amplitude is affected by density stratification and radiative diffusion for a 3 $M_{\odot}$ star at $X_{c}$ = 0.7. The top panel shows how wave amplitude changes due to density stratification and geometric effects within the radiation zone of the chosen stellar model. Across the radiation zone, density decreases by almost eight orders of magnitude while radius increases by approximately one order of magnitude. This leads to the growth of wave amplitudes as the density stratification term dominates over the geometric term. The trend also shows a sharp increase close to the convective core and a sharp decrease close to the surface convection zone caused by the Brunt–Väisälä frequency changing rapidly (by more than three orders of the magnitude). The middle panel shows the effect of radiative diffusion on IGW with an initial wave amplitude of 1 and wave number of 1. It can be observed that waves of higher frequencies are damped less than waves of lower frequencies due to the damping opacity being proportional to $\omega^{-4}$ (see Eq. (12)). In stellar interiors, IGWs are affected by both radiative damping and density stratification simultaneously and this can be seen in the bottom panel of Fig. 2. Waves that do not succumb to damping rapidly (i.e. higher frequency waves) undergo an increase in amplitude due to density stratification as they propagate through the radiation zone of a star towards the surface. Close to the surface convection zone, both increasing thermal diffusivity and rapidly decreasing Brunt–Väisälä frequency cause IGW amplitudes to decrease rapidly as seen for the case when $\omega$ = 5.0 $\mu$Hz (red line). As the Brunt–Väisälä frequency decreases towards the surface, more waves are likely to be internally reflected at radii where their frequencies are equal to the background Brunt–Väisälä frequency and possibly set up standing modes. Our analysis does not consider this scenario. Looking at the effects of density stratification and radiative damping alone, one may ask if the wave amplitudes increase enough to become nonlinear before reaching the surface.

4 Results

We calculated nonlinearity parameters using Eq. (13) and Eq. (14) for a chosen set of frequencies and wave numbers dictated by the generation spectra listed in Table 1 and present them in Fig. 3. The convective velocity ($u_{c}$), which is required to compute IGW generation wave amplitudes (see Section 2.1), was calculated by taking the average of all MESA-generated convective velocities in the grid points between the first and third quartile (in radius) of the convection zone, as shown below:

[TABLE]

where $i$ is grid space index inside the convection zone and $\mathrm{N}_{g}$ is the last grid index within the convection zone. The term $u_{i}$ refers to the convective velocity at the grid with index $i$. We label $u_{c}$ as $u_{\mathrm{MLT}}$, since the convective velocities in the convection zone are computed using mixing length theory (MLT). The convective turnover frequency, $\omega_{c}$, was set as

[TABLE]

In Fig. 3, solid black horizontal lines in both plots represent $\epsilon=1$ (strong nonlinearity condition) and $\epsilon=0.1$ (weak nonlinearity condition). It is certainly the case that waves with $\epsilon>1$ (Press, 1981) become nonlinear but nonlinear effects may also become important at lower $\epsilon$ (Barker & Ogilvie, 2010). Note that from here on, our analysis of nonlinear waves (by either the strong or weak condition) will be only on those that become nonlinear beyond 0.5 $H_{P}$ from the convective-radiative interface, which is the nominal convective overshoot depth (Rogers et al., 2013). We will refer to this as $r_{\mathrm{min,break}}$.

Since the damping opacity is proportional to $\omega^{-4}$ and $\ell^{3}$ (see Eq. (12)), higher-frequency and lower-$\ell$ waves are expected to be damped less, without consideration of generation spectra. Thus, they are more likely to become nonlinear close to the stellar surface from amplitude growth alone. When all the generation spectra are considered as shown in Fig. 3, we see that not only do the waves need to have higher frequency and lower wavenumber to become nonlinear close to the surface, they must also be generated with sufficient amplitudes. Nonlinear waves from spectrum F satisfy both these criteria for our example. On the other hand, low frequency waves that are generated with high amplitudes have a higher chance of becoming nonlinear close to the convective core. Figure 3 shows that as frequency decreases, the initial $\epsilon$ for all generation spectra increases. This means that before the effect of damping starts to dominate, these low-frequency waves can become nonlinear near the core.

4.1 IGW Generation Spectra Dependence

The different line styles in Fig. 3 represent different IGW generation spectra introduced in Section 2.1. It can be seen that when $l=1$ (top panel), only 1 out of the 14 waves shown become nonlinear under the strong condition, which is the 10 $\mu$Hz wave generated from spectrum F. The flat spectrum allows high frequency waves to be generated with high amplitude. This, coupled with the fact that high-frequency waves experience lower damping, allows the 10 $\mu$Hz wave to become strongly nonlinear close to the stellar surface. Note that zero dependence on IGW frequency or wave number means that all waves generated from spectrum F start with the same initial amplitude. However, since the nonlinearity parameter is $k_{h}u_{h}/\omega$, the initial $\epsilon$ will be different for different frequencies and horizontal wave numbers. Waves from the other three generation spectra do not become strongly nonlinear, either because they are not generated with a large enough amplitude or because they experience large damping.

Looking at waves with $\epsilon>0.1$ at $r>r_{\mathrm{min,break}}$, we can see one additional wave generated with spectrum F satisfying the weak condition. At lower frequencies, gravity waves are more likely to become nonlinear closer to the convective-radiative interface because of higher generation amplitudes and the boost due to the rapid increase of the Brunt–Väisälä frequency close to the interface. However, even though the initial amplitudes of low-frequency waves generated from spectrum K, spectrum LD and spectrum R are higher compared to their higher frequency counterparts, they are not boosted sufficiently to become nonlinear.

When $\ell$ is increased to 3 (bottom panel of Fig. 3), a few differences can be noted in the $\epsilon$ behaviour compared to that for $\ell$ = 1. First, when waves of the same frequency are considered, small-scale waves are generated with higher $\epsilon$. Spectrum LD does this most efficiently, followed by spectrum K, spectrum F and spectrum R. Second, increasing $\ell$ results in increased damping which causes the peaks of $\epsilon$ to be lower in the $l=3$ case compared to the $l=1$ case. This can lead to higher $\ell$ waves not achieving sufficient $\epsilon$ to become nonlinear. Thus, for a fixed frequency, the wave number determines which of the two effects mentioned above dominate.

Fig. 4 shows the nonlinearity parameter with a broader range of IGW horizontal wave numbers ($\ell$ = 1 – 10) and frequencies ($\omega$ = 0.5 $\mu$Hz to 100 $\mu$Hz). The colour scheme represents the maximum $\epsilon$ reached by an IGW of a given wave number and frequency at $r>r_{\mathrm{min,break}}$. The range in frequency of 0.5 $\mu$Hz to 100 $\mu$Hz was chosen to look at the range of IGW thought to be generated by convection. The spacing between wave frequencies was chosen to be equal in the logarithmic scale. Note that in Fig. 4, any nonlinearity parameter above 1.0 or below 0.01 will be red or blue respectively. The lower limit of $\epsilon$ has been set to 0.01 to provide a clear view of how $\epsilon$ varies across the weak condition. It was found that spectrum F has the highest fraction of nonlinear waves (for both strong and weak nonlinearity conditions), while spectrum K and spectrum LD show no nonlinear waves. Spectrum R shows that potentially nonlinear waves are concentrated in the lower $\ell$ regime. However, there is an irregular trend with increasing frequencies. At the lowest frequencies (0.5 $\mu$Hz - 1.0 $\mu$Hz), damping is at its highest causing most waves to be damped rapidly. Across all other tested frequencies at low-$\ell$, we see two peaks in $\epsilon$. The peak at low frequency is due to waves becoming nonlinear close to the convective core. The second peak at higher frequency is due to waves breaking at the surface. In between, waves are of high enough frequencies to not be generated with sufficient amplitude to break near the core but of low enough frequency to still be damped before reaching the surface. This damping is of course wavelength dependent, hence the structure seen in $\ell$.

In conclusion, one can see that with MLT velocities, only waves from spectrum R and spectrum F can become nonlinear.

4.2 Convective Velocity Dependence

The bulk convective velocities, $u_{c}$, were obtained from the formulation of mixing length theory (Biermann, 1951; Böhm-Vitense, 1958) in MESA. While being widely adopted, it is also well known that this prescription is likely not an accurate representation of the three-dimensional nature of turbulent convection. In particular, simulations of turbulent convection generally show significantly more intermittency than this prescription allows for. In addition to the uncertainty related to MLT, the MLT velocity itself, $u_{\mathrm{MLT}}$, can vary significantly within the convection zone. Generally, within the convection zone of a MESA model, the MLT velocity varied by approximately one order of magnitude for all the masses investigated. Furthermore, a recent work (Couston et al., 2018), which looked at 3D simulations of IGW generated by turbulent convection, has shown that at the point of generation, the steep fall-off in frequency dependence as in spectrum K and spectrum LD occurs well away from the traditional convective turnover frequency relation (a factor of 2-5). This introduces further uncertainties to a factor of 2-5 in the value of $u_{c}$. Thus, we investigated the effect of varying convective velocities by varying $u_{c}$ from one-third to three times $u_{\mathrm{MLT}}$. Fig. 5a and Fig. 5b show how the nonlinearity parameter changes with higher and lower convective velocities. In Fig. 5a, where $u_{c}$ was set to one-third of the initial convective velocity given by MESA, only waves from spectrum F show strong nonlinearity. On the other hand, when $u_{c}$ was set to 3 $u_{\mathrm{MLT}}$, waves from all spectra, except from spectrum LD, exhibit strong nonlinearity

For a larger $u_{c}$, spectrum K can be seen to exhibit less nonlinear waves at high frequencies as these waves are generated with lower amplitudes. At low frequencies and larger $\ell$, generation amplitudes are larger but damping is also larger. Moreover, the exponential term reduces amplitudes further. This causes only a small range of waves with specific frequencies to become nonlinear. Whether waves become nonlinear depends sensitively on the amplitudes of convective velocity which may very well be time-dependent. Hence, it is entirely possible that the occurence of breaking waves is also time-dependent.

4.3 IGW Energies

To better understand the impact breaking waves would have in stellar processes like mixing and angular momentum transfer, we look at the energy in waves that become nonlinear, which can be represented by the following relation:

[TABLE]

For the flat spectrum, the nonlinear wave energies are distributed equally for all wave numbers and frequencies. For spectrum R,

[TABLE]

which means that waves with the lowest wave numbers and lowest frequencies possess the highest energies. On the other hand, for spectrum LD,

[TABLE]

Therefore, the waves that carry the most energy are the high-wave number, low-frequency waves. Finally, for spectrum K,

[TABLE]

The distribution of wave energies for spectrum K is not as direct as for the other spectra due to the exponential term. In Figure 6, which shows the nonlinear wave energy trend for all generation spectra, we can see that for spectrum K and $u_{c}=u_{\mathrm{MLT}}$, the nonlinear wave energies are highest for low-$\ell$ and low-$\omega$ waves.

The ratio of energy in waves that satisfy $\epsilon\geqslant 0.1$ to the total energy of all waves within the frequency and wave number ranges considered, is

[TABLE]

where the numerator on the right-hand side represents the sum over IGW amplitudes that satisfy $\epsilon\geqslant 0.1$, at radii larger than $r_{\mathrm{min,break}}$, while the denominator is the total energy input in waves. The spacing between frequencies is represented by $\Delta\omega$. The spacing between $\ell$ is 1. Both energy terms in the numerator and the denominator are energies supplied by convective processes, so this ratio was calculated at the point of wave generation.

Applying Eq. (21) to the same stellar model for cases of $u_{c}=u_{\mathrm{MLT}}$ and $u_{c}=3\;u_{\mathrm{MLT}}$ gives the values shown in Table 3 at $Z=0.02$. Energy ratios for spectrum F are always close to one as most waves generated from spectrum F become nonlinear as shown in Fig. 4 and Fig. 5. For $u_{c}=u_{\mathrm{MLT}}$, only spectrum F shows any appreciable energy in nonlinear waves (see Fig. 4). When $u_{c}$ is increased to $3\;u_{\mathrm{MLT}}$, all the generation spectra, except spectrum LD, produce nonlinear waves. Nonlinear waves from spectrum R and spectrum K have moderately higher energy ratios (see Table 3). This due to the significant overlaps between the nonlinearity parameter trend (see Fig. 5b) and the wave energy trend (Fig. 6), for both spectrum R and spectrum K.

4.4 Metallicity Dependence

In the results discussed above, we have used a stellar model generated with an initial metallicity of 0.02, equal to the solar metallicity. To investigate the effect of different initial metallicities on the production of nonlinear waves, we generated stellar models at ZAMS with initial metallicities of 0.006 and 0.06, and tested these models with different wave generation spectra at $u_{c}$ = $u_{\mathrm{MLT}}$ and $u_{c}$ = 3 $u_{\mathrm{MLT}}$. As in the previous section, we have chosen a fixed nonlinearity condition of $\epsilon\geqslant 0.1$. We calculate energy ratios based on Eq. (21) and show them in Table 3.

The trend with spectrum F shows a slight decrease with decreasing metallicities due to increasing damping as lower metallicties lead to lower opacities (see Fig 7a). Other than spectrum F, virtually no other spectrum shows nonlinear waves for any metallicity at $u_{c}$ = $u_{\mathrm{MLT}}$. The exception is spectrum R which shows a slight increase in nonlinear wave energy at $Z$ = 0.006 because the product of the density stratification term and geometric term (see Fig. 2) becomes larger due to decreasing total stellar radius (as shown in Fig. 7) and increasing Brunt–Väisälä frequency spike (see Fig. 7b) causing wave amplification close to the stellar surface.

Moving to $u_{c}$ = 3 $u_{\mathrm{MLT}}$, we see that spectrum R and spectrum K have opposite trends in nonlinear energy ratios for decreasing metallicity. Spectrum R shows a decreasing trend because of the same reason spectrum F shows a decreasing trend. Spectrum K shows an increasing trend as convective velocities are found to be higher when stellar metallicity is lower.

4.5 Age and Mass Dependences

To investigate whether there are any trends of waves becoming nonlinear with varying stellar mass and age, we extended our analysis to include stellar masses up to 20 $M_{\odot}$ and $X_{c}$ down to 0.01, where we have used the central hydrogen mass fraction, $X_{c}$, as a proxy for stellar age. A fixed nonlinearity condition of $\epsilon\geqslant 0.1$ was chosen.

Figure 8 shows how these energy ratios change as a function of $X_{c}$, with $u_{c}=$ 3 $u_{\mathrm{MLT}}$, for 3 $M_{\odot}$, 7 $M_{\odot}$ and 20 $M_{\odot}$ main-sequence stars. From midMS onwards, no data is plotted for 20 $M_{\odot}$ (bottom panel of Fig. 8), as the star develops a convection zone in the radiation zone between the convective core and the surface convection zone. A quick analysis on waves generated from the intermediate convection zone showed that no waves become nonlinear for all spectra except spectrum F. This was due to the very small convective velocities generated in this narrow convection zone. Waves generated with such small amplitude never reach sufficient amplitudes to become nonlinear. This intermediate zone also significantly attenuates waves from the core.

In general, we see a decreasing trend of nonlinear wave energies with increasing age for spectrum F, spectrum R and initially for spectrum K. In fact, for spectrum K, this decrease leads to no waves becoming nonlinear with the effect being more prominent at lower stellar masses. To understand this, we look at various factors that affect the wave amplitude. First, Fig. 9 shows that when a star starts to undergo hydrogen burning and evolve from ZAMS, its core contracts slowly leaving behind a steep hydrogen abundance gradient (Dziembowski et al., 1991). This leads to the formation of a sharp “spike" in the Brunt–Väisälä frequency profile near the edge of convective core, which broadens and moves inwards with the convective core surface (see Fig. 9). This spike effectively traps low-frequency, high-wave number waves causing less waves to become nonlinear close to the stellar surface. Moreover, the Brunt–Väisälä frequencies are lower for older stars which means the density stratification term (from Eq. (13)) will be smaller too. Furthermore, as the star ages, stellar density varies over the same magnitude across the radiation zone but the star grows in size. This causes the density gradient to decrease and the geometric term to grow, leading to their product (as shown in Eq. (13)) being smaller in older stars. It was thus observed that despite the increase in convective velocities with age (see Fig. 10), which should lead to an increasing trend of nonlinear wave energies, the trend seen in Fig. 8 shows that the factors mentioned above dominate.

We see an increasing trend in nonlinear wave energy ratios for spectrum K at higher ages. To explain this, we look at Fig. 9, which shows that the peak Brunt–Väisälä frequency close to the stellar surface is lower for older stars. This causes higher frequency waves to experience lower damping close to the surface which allows more waves to cross the $\epsilon>0.1$ threshold.

Figure 11 shows wave energy ratios as a function of mass with the top panel representing ZAMS models, the middle panel representing midMs models and the bottom panel representing TAMS models. As in the case for age dependence, energy ratios of waves from spectrum R are higher than energy ratios of waves from spectrum K. In general, the nonlinear energy ratios for all spectra increase slightly with increasing mass. This can be explained by the increasing convective velocities with mass (see Fig. 10). Furthermore, Fig. 9 shows that for larger stellar masses, the Brunt–Väisälä frequencies are lower, which leads to lower damping.

5 Discussion and Conclusion

In this paper, we have shown that depending on the spectrum and amplitudes of waves generated by convection, some waves may become nonlinear in the radiative regions of intermediate-mass stars. This can occur because the amplitudes of high(er) frequency/longer wavelength waves are amplified by density stratification. Therefore, if such waves are generated with sufficient amplitude, they can become nonlinear along their journey towards the stellar surface. On the other hand, lower frequency waves, which are generated with larger amplitudes, can become nonlinear close to the convective core. These waves that break near the convective core may lead to enhanced mixing beyond a nominal “overshoot depth”.

One crucial element as to whether waves can become nonlinear and break in stellar radiative regions is the wave generation spectrum. Compared to the theoretical generation spectra from Lecoanet & Quataert (2013) and Kumar et al. (1999), with their steep fall off at frequencies higher than the convective turnover frequency, the numerical generation spectra from Rogers et al. (2013) show that more waves may become nonlinear and consequently, more energy is in these nonlinear waves for MLT convective velocities.

The assumed convective velocity also plays an important role in the ability of waves to become nonlinear, with even a factor of three increase causing significant differences in total nonlinear wave energies. One such significant difference occurs when convective velocities are set as three times the MLT values. In this case, nonlinear wave energy ratios from the theoretical generation spectrum of Kumar et al. (1999) became comparable to those from the numerical generation spectrum R. Given the likelihood of variation in convective velocities by at least this order, it is worth trying to understand the role that intermittency plays in determining the wave spectrum. This is important especially after the results from a recent numerical work from Couston et al. (2018) which shows that the steep fall-off in energy with frequency dependence occurs away from the classical theory for convective turnover frequency. Indeed, it may be possible that the wave generating process at the convective-radiative interface and wave breaking process at the surface may be time-dependent.

We found that age, mass and metallicity also affect the fraction of waves becoming nonlinear, but to a lesser degree. Generally, as stars become older the fraction of energy that goes into nonlinear waves decreases because the receding convective core leaves behind a compositional gradient and hence, spike in the Brunt–Väisälä frequency which effectively filters waves from propagating outward. As the mass increases, the fraction of nonlinear waves increases slightly due to increasing convective velocities. As the metallicity increases, the fraction of nonlinear waves behaves differently, depending on generation spectra.

Whether or not waves can become nonlinear in the radiative regions of stars plays an important role on the efficiency with which waves can transport momentum and mix chemical elements. Clearly, wave breaking and the ensuing turbulence that is generated would lead to more efficient mixing and angular momentum transport than simple radiative dissipation. The strong observational evidence for enhanced coupling between convective and radiative regions (Beck et al., 2012; Deheuvels et al., 2012; Tayar & Pinsonneault, 2013) and the possible observational detections of IGW(Aerts & Rogers, 2015; Aerts et al., 2017, 2018; Ramiaramanantsoa et al., 2018) both point to breaking IGW, which in turn, points to a flatter frequency generation spectrum from convection than what theoretical models predict or higher convective velocities than those predicted by MLT. It is therefore vital that more theoretical and numerical investigations into the convective generation spectrum of IGW, and the physics that drives it, are undertaken. In particular, it is worth further studying plume models, which depending on parameter choices, can produce generation spectra with flatter dependences on frequencies (Pinçon et al., 2016).

Acknowledgments

We acknowledge support from STFC grant ST/L005549/1 and NASA grant NNX17AB92G. Computing was carried out on Pleiades at NASA Ames. We thank Conny Aerts, Daniel Lecoanet, May Gade Pedersen and Dominic Bowman for useful conversations leading to the development of this manuscript. Finally, we would like to thank Manoj Previn Ratnasingam and Romy Klaasje Verhoeven for the proofreading support they provided.