Three-Dimensional Simulations of Massive Stars: I. Wave Generation and Propagation

TL;DR

This study uses 3D hydrodynamic simulations to explore wave generation and propagation in massive stars, revealing insights into internal gravity waves, their spectra, and surface fluctuations, with implications for stellar oscillation observations.

Contribution

First 3D simulations of intermediate-mass star convection zones capturing wave generation, propagation, and comparison with theoretical and observational spectra.

Findings

Convection dominated by plumes with a double power law spectrum.

Wave spectra consistent with plume penetration excitation.

Surface fluctuations linked to observable brightness variations.

Abstract

We present the first three-dimensional (3D), hydrodynamic simulations of the core convection zone (CZ) and extended radiative zone spanning from 1% to 90% of the stellar radius of an intermediate mass (3 $\mathrm{M}_\odot$) star. This allows us to self-consistently follow the generation of internal gravity waves (IGWs) at the convective boundary and their propagation to the surface. We find that convection in the core is dominated by plumes. The frequency spectrum in the CZ and that of IGW generation is a double power law as seen in previous two-dimensional (2D) simulations. The spectrum is significantly flatter than theoretical predictions using excitation through Reynolds stresses induced by convective eddies alone. It is compatible with excitation through plume penetration. An empirically determined distribution of plume frequencies generally matches the one necessary to explain a large part of the observed spectrum. We observe waves propagating in the radiation zone and excited standing modes, which can be identified as gravity and fundamental modes. They show similar frequencies and node patterns to those predicted by the stellar oscillation code GYRE. The continuous part of the spectrum fulfills the IGW dispersion relation. A spectrum of tangential velocity and temperature fluctuations close to the surface is extracted, which are directly related to observable brightness variations in stars. Unlike 2D simulations we do not see the high frequencies associated with wave breaking, likely because these 3D simulations are more heavily damped.

Tables

Table 1. Table 1: List of 3D simulations

ModelaaModel names are built up from the luminosity boosting factor, rotation rate, and angular resolution.	spectral modes	$\nu /{10}^{13}{\mathrm{cm}}^2{s}^{-1}$	$\kappa$	$\overline{Q}$	$\Omega /{10}^{-6}\mathrm{rad}{s}^{-1}$	RaffThe definitions are given in Eqs. (8) – (10).	ReffThe definitions are given in Eqs. (8) – (10).	PrffThe definitions are given in Eqs. (8) – (10).	$t_{\mathrm{sim}}$ggThis is the total physical runtime of the simulation./d
H6R5	3741bbThis corresponds to angular resolution of 128 ($\theta$) and 256 ($\varphi$).	$8$ (risingddThe value given is the one in the core. In the envelope $\nu$ rises with $\nu \propto {\rho }^{-1/4}$. The profile is continuous.)	${10}^5{\kappa }_{\star }$	${10}^6{\epsilon }_{\star }$	$5$	$8\text{×}{10}^{11}$	126	60 – $0.02$eeThe first value applies to the core, the second to the envelope.	$38.7$
H6R10	3741bbThis corresponds to angular resolution of 128 ($\theta$) and 256 ($\varphi$).	$8$ (risingddThe value given is the one in the core. In the envelope $\nu$ rises with $\nu \propto {\rho }^{-1/4}$. The profile is continuous.)	${10}^5{\kappa }_{\star }$	${10}^6{\epsilon }_{\star }$	$10$	$8\text{×}{10}^{11}$	126	60 – $0.02$eeThe first value applies to the core, the second to the envelope.	$59.3$
H5	3741bbThis corresponds to angular resolution of 128 ($\theta$) and 256 ($\varphi$).	$1$	${10}^5{\kappa }_{\star }$	${10}^5{\epsilon }_{\star }$	0	$6\text{×}{10}^{11}$	468	7 – ${10}^{-6}$eeThe first value applies to the core, the second to the envelope.	$13.3$
H6E	3741bbThis corresponds to angular resolution of 128 ($\theta$) and 256 ($\varphi$).	$10$	${10}^5{\kappa }_{\star }$ – $50{\kappa }_{\star }$eeThe first value applies to the core, the second to the envelope.	$\approx {10}^6{\epsilon }_{\star }$ (exp.)	0	${10}^{12}$	101	100 – 2eeThe first value applies to the core, the second to the envelope.	$61.7$
H6LD	3741bbThis corresponds to angular resolution of 128 ($\theta$) and 256 ($\varphi$).	$10$	${10}^5{\kappa }_{\star }$ – $50{\kappa }_{\star }$eeThe first value applies to the core, the second to the envelope.	${10}^6{\epsilon }_{\star }$	0	${10}^{12}$	100	100 – 2eeThe first value applies to the core, the second to the envelope.	$58.8$
H6LD-HR	14706ccThis corresponds to angular resolution of 256 ($\theta$) and 512 ($\varphi$).	$10$	${10}^5{\kappa }_{\star }$ – $50{\kappa }_{\star }$eeThe first value applies to the core, the second to the envelope.	${10}^6{\epsilon }_{\star }$	0	${10}^{12}$	100	100 – 2eeThe first value applies to the core, the second to the envelope.	$6.5$
H7E	3741bbThis corresponds to angular resolution of 128 ($\theta$) and 256 ($\varphi$).	$8$ (risingddThe value given is the one in the core. In the envelope $\nu$ rises with $\nu \propto {\rho }^{-1/4}$. The profile is continuous.)	$5\text{×}{10}^{13}{\mathrm{cm}}^2{s}^{-1}$	$\approx {10}^7{\epsilon }_{\star }$ (exp.)	0	$2\text{×}{10}^9$	272	1 – 39eeThe first value applies to the core, the second to the envelope.	$16.4$
H7E-HR	14706ccThis corresponds to angular resolution of 256 ($\theta$) and 512 ($\varphi$).	$8$ (risingddThe value given is the one in the core. In the envelope $\nu$ rises with $\nu \propto {\rho }^{-1/4}$. The profile is continuous.)	$5\text{×}{10}^{13}{\mathrm{cm}}^2{s}^{-1}$	$\approx {10}^7{\epsilon }_{\star }$ (exp.)	0	$2\text{×}{10}^9$	272	1 – 39eeThe first value applies to the core, the second to the envelope.	$33.2$

Equations

$$\nabla \cdot \overline{\rho }\text{mathitbf}v=0,$$

$$\frac{\partial \text{mathitbf}v}{\partial t}$$

$$+\overline{\nu }\left({\nabla }^2\text{mathitbf}v+\frac{1}{3}\nabla (\nabla \cdot \text{mathitbf}v)\right),$$

$$\frac{\partial T}{\partial t}$$

$$-v_r\left(\frac{\partial \overline{T}}{\partial r}-(\gamma -1)\overline{T}{h}_{\rho }\right)+\frac{\overline{Q}}{c_v\overline{\rho }}$$

$$+\frac{1}{c_v\overline{\rho }}\nabla \cdot (c_p\overline{\kappa }\overline{\rho }\nabla T)+\frac{1}{c_v\overline{\rho }}\nabla \cdot (c_p\overline{\kappa }\overline{\rho }\nabla \overline{T}).$$

$$C=-\frac{1}{\overline{T}}\left(T+\frac{1}{\overline{g}\overline{\rho }}\frac{\partial \overline{T}}{\partial z}p\right).$$

$$\overline{\rho }\text{mathitbf}v=\nabla \times \nabla \times W\text{mathitbf}\hat{r}+\nabla \times Z\text{mathitbf}\hat{r}.$$

$$T(r,\theta ,\varphi ,t)=\sum _{m=-m_{\max }}^{m_{\max }}\sum _{l=|m|}^{l_{\max }}{T}_{l,m}(r,t)Y_{l,m}(\theta ,\varphi ),$$

$$T_{l,-m}=(-1{)}^m{T}_{l,m}^{\ast },$$

$$\mathrm{Ra}=\frac{g\overline{Q}{D}^5}{c_v{\kappa }^2\nu \overline{T}},$$

$$\mathrm{Re}=\frac{v_{\text{rms}}D}{\nu },$$

$$\mathrm{Pr}=\frac{\nu }{\kappa },$$

$$\overline{Q}=A{c}_v\overline{\rho }\exp (-r/r_{\min })(r-r_{\min })/R_{\star },$$

$$\hat{E}(f)=\frac{1}{t_b-t_a}{\int }_{t_a}^{t_b}E(t)e^{-2\pi ift}dt.$$

$$\hat{E}(f_j)=\frac{1}{n}\sum _{m=0}^{n-1}E(t_m)e^{-2\pi i\frac{mj}{n}},$$

$$L\propto v_{\text{rms}}^3.$$

$$E_l(r)$$

$$\times \left(\frac{l(l+1)}{r^2}|W_l^m{|}^2+{\left|\frac{\partial W_l^m}{\partial r}\right|}^2+|Z_l^m{|}^2\right),$$

$${\hat{E}}_{\mathrm{kin}}=\frac{1}{2}\bar{\rho }\left({\hat{v}}_r^2+{\hat{v}}_{\theta }^2+{\hat{v}}_{\varphi }^2\right),$$

$$f_{\text{TO}}=\frac{v_{\text{rms}}}{\pi r_{\text{CZ}}},$$

$$s=$4.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$\cdot l.$$

$$f_{\text{b}}\sim \frac{v_{\text{pl}}}{{\Delta }_{\text{p}}}.$$

$$E(f)\propto \exp \left(-(f/f_{\text{b}}{)}^2\right),$$

$$A{\int }_c^{\infty }{e}^{-f_{\text{b}}/\alpha }{e}^{-(f/f_{\text{b}}{)}^2}d{f}_{\text{b}}$$

$${\Delta }_{\text{p}}^2=\frac{3}{5}{H}_P{H}_{\kappa }f\frac{\rho v_{\text{pl}}^3}{F_{\text{tot}}},$$

$$v_r\propto {\left(\frac{r_0}{r}\right)}^{3/2}\sqrt{\frac{{\rho }_0}{\rho }}{\left(\frac{N^2-{\omega }^2}{N_0^2-{\omega }^2}\right)}^{1/4}\exp (-\tau /2),$$

$$\tau ={\int }_{r_0}^{r}dr\frac{\kappa {\left[l(l+1)\right]}^{3/2}{N}^3}{r^3{\omega }^4}\sqrt{1-\frac{{\omega }^2}{N^2}},$$

$$\frac{k_{\perp }}{k}=\frac{\omega }{N}.$$

$$k_{\perp }=\frac{\sqrt{l(l+1)}}{r}.$$

$$\max (\overline{\kappa },\overline{\nu })\sim \frac{{\lambda }^2}{{\tau }_{\mathrm{d}}}=\frac{(2\pi {)}^5{r}^2{f}^4}{N^{3}l(l+1)},$$

$$v_g=\frac{\partial \omega }{\partial k}=\frac{r{\omega }^2}{N\sqrt{l(l+1)}}.$$

$${\hat{v}}_{r,\theta ,\varphi }^{\mathrm{filt}}=\{\begin{array}{ll}{\hat{v}}_{r,\theta ,\varphi }& \text{if NωCcut1<kk⊥<NωCcut,}\\ 0& \text{otherwise.}\end{array}$$

Full text

hydrodynamics - stars: interiors - convection - waves

††journal: ApJ††software: matplotlib (Hunter, 2007), scipy (Jones et al., 2001–), MESA (Paxton et al., 2011, 2013, 2015, 2018), GYRE (Townsend & Teitler, 2013)

\savesymbol

tablenum \restoresymbolSIXtablenum

1 Introduction

In addition to sound waves, fluid dynamical systems can have other wave-like solutions for which the restoring force is not pressure but buoyancy. These waves are commonly referred to as internal gravity waves (IGWs) to distinguish them from surface gravity waves. They occur in many stratified systems, such as atmospheres and oceans, in many of which they have an important impact on the large scale dynamics. IGWs excited by equatorial convection were found to be crucial in driving the quasi-biennial oscillation (QBO) in the Earth’s equatorial stratosphere (Baldwin et al., 2001). In the oceans, IGWs excited through the surface wind or tides cause turbulent mixing when they break (Munk & Wunsch, 1998).

In stars IGWs have been suggested to play an important role in angular momentum transport and chemical mixing in radiative regions, where other mechanisms are not efficient. Press (1981) suggested that IGWs in the sun can cause mixing in the convectively stable interior and affect the effective radiative opacity by a factor of two or more. IGW mixing was also suggested as the cause of lithium depletion in F stars (Garcia Lopez & Spruit, 1991) and in the sun (Schatzman, 1993; Montalbán, 1994; Talon & Charbonnel, 2005).

IGWs are candidates for being the cause of some of the observed properties of stars that are poorly explained by current stellar models, such as the internal rotation structure of stars (Beck et al., 2012; Aerts et al., 2017a), stellar cores counter-rotating to their envelopes (Triana et al., 2015; Rogers, 2015), or the enhanced mass loss needed to explain certain classes of core-collapse supernovae (Quataert & Shiode, 2012). Photometric observations suggest the presence of convectively generated IGWs in, at least, some massive stars since the observed velocity spectrum at the surface compares well to that obtained using numerical simulations of IGWs (Aerts & Rogers, 2015; Bowman et al., 2019).

To understand the role IGWs play in all these physical situations it is important to know what spectrum of waves in frequency and wave number space is excited by convection. Theoretical work characterizing these spectra mostly focuses on two mechanisms, excitation through the Reynolds stresses of convective eddies or through penetration of plumes. The former approach was taken by Lighthill (1952), Goldreich & Kumar (1990), Kumar et al. (1999), and later by Lecoanet & Quataert (2013). All these studies found a power law dependence in frequency, i.e. proportional to $f^{-\alpha}$, with wave frequency $f$ and exponent $\alpha$. The exact value of the exponent depended on the profile of the Brunt–Väisälä frequency at the convective boundary (CB). The spectrum generated by plume penetration was first studied by Townsend (1966) in a terrestrial context and later extended to stars by Montalbán & Schatzman (2000). A recent semianalytical model for the IGW flux caused by plumes at the base of a convection zone, as is the case in the sun, has been developed by Pinçon et al. (2016). The predicted spectrum takes a very similar functional form in all these plume-driven cases, which is proportional to $\exp[-(f/f_{\text{b}})^{2}]$, with wave frequency $f$ and the plume frequency $f_{\text{b}}$.

Multidimensional hydrodynamic simulations generally do not impose a specific IGW generation mechanism, and are able to follow convection, IGW generation and propagation directly from the basic equations. Yet numerical limitations and the extreme scales within stellar interiors often restrict them to a more dissipative regime than is realistic in stars. Nevertheless, careful choice of parameters and interpretation of the results allow us to assess theoretical predictions. Simulations showed that the Li depletion in the sun cannot be explained by IGWs (Rogers et al., 2006). Similarly the uniform rotation of the sun’s radiative interior is not completely caused by IGWs (Rogers & Glatzmaier, 2005; Denissenkov et al., 2008). Rogers et al. (2013) performed two-dimensional (2D) hydrodynamic simulations of IGW generation at the boundary of convective cores of massive stars. They found that the IGW generation spectrum is generally much shallower than theoretical predictions. It shows two frequency regimes with different slopes, suggesting different excitation mechanisms at work. Recent research on breaking of IGWs in the radiative envelopes of massive stars affirmed the importance of the shape of the wave generation spectrum (Ratnasingam et al., 2019).

Browning et al. (2004) performed simulations of the inner 30% in radius of a $2\,\mathrm{M}_{\odot}$ star with methods very similar to the ones used in this work. Their work focused on convective motions in the core, overshooting, and the influence of rotation. They do not study IGWs in detail, but mention their excitation at the convective boundary. Later work by the same group included magnetic fields and specifically studied the dynamo in the convective core (Brun et al., 2005). In contrast, our work specifically studies IGW excitation and propagation and therefore includes a much larger part of the radiation zone (up to 90% in radius).

In their work on IGWs in solar-like stars, Alvan et al. (2014) performed a detailed analysis of wave excitation and propagation, similar to the one carried out in our work. The main difference is their work is based on solar-like stars with a convective envelope and radiative core, while the opposite is the case in our $3\,\mathrm{M}_{\odot}$ star. Propagation through a radiative envelope along a falling density gradient causes wave amplification, which makes nonlinear behavior more likely in intermediate-mass and massive stars.

The remaining parts of this paper are structured as follows: Section 2 describes the hydrodynamic equations solved and their pseudo-spectral discretization. Section 3 discusses the stellar models used as the background state of the simulations and assumptions on heating and dissipation needed for numerical reasons. The general properties of three-dimensional (3D) convection in the core are presented in Sect. 4.1. Frequency spectra of core convection and their implications for the generation of IGWs are discussed in Sect. 4.2. The properties of the overshooting region is the subject of Sect. 4.3. Section 4.4 treats IGW propagation and the modes excited in the radiation zone, while Sect. 4.6 discusses the signature they are expected to leave on the surface, which could be observed by photometry and spectroscopy. We conclude in Sect. 5.

2 Computational Method

The simulations presented here are a logical continuation of those of Rogers et al. (2013). One caveat of their work is the restriction to 2D geometry, which is expected to yield significantly different behavior of turbulence and also altered wave propagation to some degree due to the difference between 2D annulus geometry and a 3D sphere. We extend their method to 3D by using the same type of anelastic approximation, but discretizing the horizontal part of the equation in terms of spherical harmonics instead of $\sin$ and $\cos$ functions.

We solve the following equations for the deviation from the reference state (indicated by a bar, e.g., $\overline{\rho}$) given by a hydrostatic stellar evolution model,

[TABLE]

[TABLE]

Here, $\mathitbf{v}$ is the 3D fluid velocity, $v_{r}$ its radial component, $\overline{\rho}$ is the background density, $\gamma$ is the adiabatic index of the ideal gas equation of state, $\overline{T}$ and $T$ are the temperature background and fluctuation, $\overline{\kappa}$ and $\overline{\nu}$ are the thermal and viscous diffusivities, $\overline{Q}$ is the energy release rate, $c_{v}$ is the specific heat at constant density, $\overline{g}$ is gravitational acceleration, and $h_{\rho}=\partial\ln\overline{\rho}/\partial r$ is the negative inverse of the density scale height. We use a standard spherical coordinate system with radius $r$, colatitude $\theta$, and azimuthal angle $\phi$. The unit vector $\mathitbf{\hat{r}}$ points in radial direction. Rotating stars are set up using a rotating frame of reference with an angular velocity $\Omega$ and rotation axis along $\mathitbf{\hat{z}}$ in direction of the pole at $\theta=0$.

This formulation of the anelastic equations includes self-gravity perturbations $\Phi$ to the reference state gravitational potential $\overline{\Phi}$ by introducing the reduced pressure $P=p/\overline{\rho}+\Phi$ and co-density $C$ (Braginsky & Roberts, 1995; Rogers & Glatzmaier, 2005). This introduces no additional computational effort as long as the thermodynamic pressure is not calculated. The co-density takes the form

[TABLE]

In their comparison of different variants of the anelastic approximation Brown et al. (2012) also investigated this variant of the anelastic equations111They call this set of equations the RG equations.. They find that the equations do not conserve energy for non-isothermal stratifications and suggest the removal of the term proportional to $p$ in Eq. (4) to ensure energy conservation. While all simulations used for the analysis in Sect. 4 did not contain the Brown et al. (2012) modification, we performed several test calculations including it. The modes we find in the radiation zone are not affected by the inclusion or exclusion of this factor.

Brown et al. (2012) also predict that IGW frequencies are larger by a factor of $\sqrt{\gamma}$. When comparing the frequencies generated in our simulations to those generated with the 1D pulsation code GYRE (Townsend & Teitler, 2013), we find small deviations, initially of the order of a few $\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ ($\sim$ 2% relative deviation) getting larger at higher wavenumbers, but this $\sqrt{\gamma}$ factor does not explain the differences. Hence, we are unsure how this factor manifests itself in our simulations or how that work extends to these fully nonlinear simulations.

The numerical solution method we choose is similar to the approach taken by Glatzmaier (1984) and in the ASH code (Clune et al., 1999) with some different choices adapted to the application at hand. To implicitly fulfill Eq. (1) we replace the mass flux $\overline{\rho}\mathitbf{v}$ by its decomposition into a poloidal ($W$) and toroidal stream function ($Z$). These are related to the mass flux by

[TABLE]

Resulting purely from a curl of a vector this is naturally divergence free. Together with temperature $T$ and reduced pressure $P$ these form the four unknown quantities we are solving for. They are expressed as a linear combination of spherical harmonics $Y_{l,m}$ and radius-dependent, complex coefficients. For temperature this is

[TABLE]

and its equivalent for the other quantities. This allows us to compute horizontal derivatives via computationally inexpensive recursion relations and it avoids the singularities at the poles. As the coefficients of real-valued quantities fulfill

[TABLE]

only the components with $m\geq 0$ need to be stored. The choice of triangular truncation ($l_{\mathrm{max}}=m_{\mathrm{max}}$) results in uniform angular resolution. The method we use is pseudo-spectral, i.e. the linear terms are computed in spectral space and the nonlinear terms are computed in grid space. This approach makes it necessary to chose the number of grid points in latitudinal and longitudinal direction, $N_{\theta}$ and $N_{\phi}$, corresponding to the number of spectral modes. To avoid aliasing errors we set $N_{\theta}=(3l_{\mathrm{max}}+1)/2$ and $N_{\phi}=2N_{\theta}$. Details on this kind of spectral discretization can be found in Glatzmaier (2013). We do not use a spectral basis in the radial direction to be able to easily adjust the grid to the underlying stellar model. In the present case we have an increased radial resolution in the convection zone. Radial derivatives are computed using second-order finite differences accounting for the nonuniform grid spacing.

We use the implicit Crank–Nicolson method for the linear diffusion terms to avoid the strict CFL condition that depends quadratically on the step size associated with explicit time-stepping. The nonlinear terms are calculated using the explicit Adams–Bashforth linear multistep method, which makes the method second-order accurate in time. We choose a constant time step of $1\text{\,}\mathrm{s}$, which is well below the CFL condition of the explicit terms and makes later Fourier analysis of the time series easier.

The code is parallelized using a domain decomposition in the radial coordinate only. The communication involves halo updates for computing finite differences in the radial direction and all-to-all communication for solving the linear equations involved in implicit time-stepping. It is implemented using the message passing interface (MPI). To alleviate the problem that domains become small when using many cores we additionally implement thread-based parallelization using OpenMP, which starts to be more efficient than pure MPI when there are less than 3 radial points per task (see Fig. 1). The achieved scaling efficiency from 100 to 1500 cores is 77% on the NASA Pleiades system.

3 Simulations

The equations discussed in Sect. 2 rely on a spherically symmetric reference state for the thermodynamic variables on top of which the evolution of small perturbations is calculated. We use the MESA (Modules for Experiments in Stellar Astrophysics) stellar evolution code222The MESA version used was SVN revision number 10000. (Paxton et al., 2011, 2013, 2015, 2018) to produce the reference state. We use the default settings to generate a nonrotating, $3\,\mathrm{M}_{\odot}$ zero-age main-sequence (ZAMS) star of metallicity $Z=10^{-2}$. The exact code configuration (inlists) and MESA profiles can be obtained at this URL333https://www.mas.ncl.ac.uk/~npe27/igw3d/. No convective overshooting was used. The values of density, temperature, and gravity are adopted unchanged from the model and interpolated onto a grid with 400 cells in the convection zone and 1100 cells in the radiation zone. The total radius of the star is $R_{\star}=$1.42\text{\times}{10}^{11}\text{,}\mathrm{c}\mathrm{m}$=2.05\,\mathrm{R}_{\odot}$.

Figure 2 shows the density and temperature profile of the stellar model. The radial extent of the 3D simulation domain indicated by vertical, dashed lines is limited at 1% of the stellar radius to avoid the coordinate singularity at the core and at 90%, where density drops below $10^{-4}\,\mathrm{g\,cm^{-3}}$, covering nearly six orders of magnitude in density.444Models H7E and H7E-HR were run with an earlier version of the stellar model, which reached 95% of stellar radius. There is no qualitative change in the wave spectra of these models. The Brunt–Väisälä frequency profile, which governs the propagation of IGWs, is plotted in Fig. 3. It shows the convective radiative boundary at 13% of the total radius.

For numerical stability we need to increase the thermal diffusivity $\kappa$ and kinematic viscosity $\nu$ beyond their physical values in the star, $\kappa_{\star}$ and $\nu_{\star}$, respectively. In the $3\,\mathrm{M}_{\odot}$ MESA model $\kappa_{\star}$ ranges from ${10}^{7}\text{\,}\mathrm{c}\mathrm{m}^{2}\,\mathrm{s}^{-1}$ in the CZ to ${10}^{12}\text{\,}\mathrm{c}\mathrm{m}^{2}\,\mathrm{s}^{-1}$ at the top of the simulated region ($r=0.9\,R_{\star}$), and $\nu_{\star}$ ranges from $60\text{\,}\mathrm{c}\mathrm{m}^{2}\,\mathrm{s}^{-1}$ to $5\text{\times}{10}^{4}\text{\,}\mathrm{c}\mathrm{m}^{2}\,\mathrm{s}^{-1}$. As increased diffusivity and viscosity would damp convection too strongly in order to reach a somewhat turbulent state, we increase the luminosity of the star by a similar factor to balance the increased damping. In a series of models we explore the effect of increased forcing and that of using different profiles for viscosity. These are summarized in Table 1. To put this into context we compare several characteristic nondimensional numbers. The Rayleigh number,

[TABLE]

with a typical length scale $D$ (chosen to be the size of the convective core in this case), controls the details of convection and determines if energy transport is mostly through radiation or convection. This particular form of $\mathrm{Ra}$ is also called a flux Rayleigh number. The stellar value is ${10}^{28}$, which is more than six orders of magnitude higher than the values reached in the simulations. This is the rationale for increasing the energy release. If we had used the original value of $\overline{Q}$, $\mathrm{Ra}$ would be approximately ${10}^{6}$, which might even be subcritical. The actual convection in the star is likely even more vigorous and plume dominated than that observed in the simulations.

Flows with a high Reynolds number,

[TABLE]

develop turbulence, while low values of $\mathrm{Re}$ normally result in laminar flow. Due to the extreme length scales $D\approx 14\%\,R_{\star}$ and velocities, $\mathrm{Re}$ is typically extremely large in stellar environments, in the current case $\mathrm{Re}\approx${10}^{12}$$. These parameters are not currently possible in numerical simulations. As can be seen in Table 1 we can only reach values of approximately ${10}^{2}$ in the CZ. However, it is expected that as long as a part of the inertial range of the turbulent cascade is numerically resolved, the energy dissipation rate will not change significantly at higher $\mathrm{Re}$ (e.g., Frisch, 1995, Chapter 5). Yet the small scale velocity field will definitely show differences, which is a caveat of the presented simulations. The Reynolds number is another reason for using an increased convective forcing, as using the original value would result in velocities corresponding to $\mathrm{Re}\approx 1$.

The Prandtl number,

[TABLE]

is the ratio of viscous to thermal diffusion. In stars it is typically extremely low, ranging from ${10}^{-6}$ in the core to ${10}^{-9}$ at the surface. The only way to reach these values in our numerical simulations would be to increase $\kappa$ to very high values, which would damp the waves too much. As a compromise we settle on $\mathrm{Pr}$ around 1 in the envelope and around 100 in the core in most models (see Tab. 1). In a few models $\mathrm{Pr}$ reaches much lower values of $0.02$ or even ${10}^{-6}$ in the envelope, but these are subject to excessive damping due to too much thermal diffusion.

In most models we increase luminosity by setting the heating function $\overline{Q}$ to the nuclear energy generation rate from MESA multiplied by a constant factor. These models are referred to with a name starting with “H$X$” for $10^{X}$ times the stellar luminosity $L_{\star}$. For example, “H6” corresponds to a luminosity of $10^{6}\,L_{\star}$. In a few models we used an exponential heating function,

[TABLE]

with a scaling factor $A$, which is used to adjust it to a boosted stellar luminosity. These are labeled with “H$X$E” for an exponential heating profile corresponding to a luminosity of $10^{X}$ the stellar value, and $r_{\mathrm{min}}$ the innermost radius of the simulation domain.

The thermal diffusivity is treated in a similar way by multiplying the stellar value with a constant factor, which was the lowest value that did not show stability problems. As the increased diffusivity is mainly needed in the convection zone, we also tried a different approach where just the CZ is subject to a value of $10^{5}\kappa_{\star}$, while diffusivity in the radiation zone can be reduced to $50\kappa_{\star}$. Both regions are blended using a hyperbolic tangent function with a width of $5\text{\times}{10}^{9}\text{\,}\mathrm{c}\mathrm{m}$ (3.5% of the stellar radius and 26% of the size of the convection zone). This was used in models H6E and H6LD. Model H6LD is a combination of the low diffusivity of Model H6E with the boosted MESA energy release of model H6.

Important conclusions in this paper are drawn from spectra. To clarify their interpretation we give an exact definition here. For a real function $E(t)$ sampled in an interval $[t_{a},t_{b}]$ the Fourier transform is

[TABLE]

By normalizing with the length of the interval the units of $\hat{E}$ are the same as those of $E$, which makes it easier to interpret the magnitude of components of the spectra.

As data from the simulations is sampled at discrete times $t_{0},\ldots,t_{n-1}$ with equidistant spacing $\Delta t$, we approximate Eq. (12) with a discrete Fourier transform (DFT),

[TABLE]

where $j$ takes values from 0 to $\lfloor\frac{n}{2}\rfloor$. Higher values of $j$ are redundant due the real input data. The corresponding frequencies are $f_{j}=\frac{j}{n\Delta t}$.

4 Results

4.1 Convection Zone

As we have to increase the heating term $\overline{Q}$ (equivalent to an increase in luminosity $L$), thermal diffusivity $\kappa$, and kinematic viscosity $\nu$ for numerical reasons, the convective velocities are higher than those predicted by mixing-length theory (MLT) using quantities from the stellar evolution model. Figure 4 compares the angular average of velocities of the different models to the MLT value. In the convection zone all our simulations have velocities one to two orders of magnitude higher than the MLT value. This causes waves at the convective–radiative boundary (CB) to be excited at higher amplitudes, which is intended to offset the increased dissipation within the RZ with the hope of surface amplitudes being more realistic. The rise of velocity close to the largest radii is related to the outer boundary condition.

The scaling of convective velocities with changing luminosity has been subject of previous studies. Other hydrodynamic simulations of convection zones in stars (Porter & Woodward, 2000; Viallet et al., 2013; Jones et al., 2017) find that,

[TABLE]

This is also the result found using MLT (e.g., Kippenhahn et al., 2012). The scaling relation agrees perfectly with the observations in our simulations, which fit $v_{\text{rms}}\propto L^{0.34}$ (see Fig. 5).

As expected from the stratification of the 1D reference state, convection immediately starts to develop in the core. From early times convection is dominated by large plumes. These plumes often rise until they reach the CB, but are sometimes dissolved by interacting with large eddies. Their disintegration at the convective boundary perturbs the stably stratified radiation zone directly above. This process can be seen in the time series in Fig. 6.

Figure 7 illustrates the correlation of positive radial velocity and temperatures higher than the horizontal average. It shows that in model H6R10 plumes reach typical velocities of $10\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$. Scaling this down to the rms convective velocities of the actual star using Eq. (14) yields a rising speed of $0.1\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$. Figure 8 shows a 3D view of the whole star using the same model.

The series of meridional slices in Fig. 6 show an example of several plumes hitting the convective boundary and triggering wave motion in the region above. Between $t=0$ and $t=$2.8\text{,}\mathrm{h} the large plume in the bottom part of the slice splits into two parts, which subsequently cause small-scale disturbances in the previously much more uniform temperature field of that region. At $t=$5.6\text{\,}\mathrm{h} a larger plume hits the boundary in the upper left corner of the convection zone. It spreads out at the boundary over more than half a hemisphere and causes Kelvin–Helmholtz-like vortices on its inner side. These seem to be the cause of many of the small-scale eddies at the interface, which can themselves drive waves in the RZ.

Turbulent kinetic energy in the CZ shows a typical cascade behavior, where most energy is present at low wavenumbers, i.e. large length scales. Figure 9 shows the kinetic energy spectrum of several models as a function of $l$ mode. The energy contained in a single $l$ mode is computed from the poloidal ($W$) and toroidal ($Z$) decomposition (see Eq. (5)) with the expression (e.g., Glatzmaier, 2013, Sect. 10.6.6),

[TABLE]

where the primed sum means that the $m=0$ terms are multiplied by $1/2$. In all cases most energy is contained in the low-order modes ($l\lesssim 5$), although the actual peak of the spectrum varies between $l=1$ and $l=3$ for the different parameters.

Although numerical diffusivity limits the inertial range in these spectra, we can still obtain a power law slope for each of the models. The slope becomes negatively steeper with increased convective forcing. In Fig. 9 we fit the inertial range of each model with power laws. In the strongly forced models H6LD and H6R10, in which we see a strong influence of rising plumes (see Fig. 6), the kinetic energy spectrum drops with $l^{-2.1}$ or $l^{-2.3}$, respectively. This approaches the value predicted by Bolgiano–Obukhov scaling of $l^{-2.2}$ for buoyancy-driven turbulence (Obukhov, 1959; Bolgiano, 1959). The more strongly forced models H7E and H7E-HR, show an even steeper slope in the inertial range, following $l^{-3.0}$. This is significantly steeper than the $l^{-5/3}$ relation predicted by Kolmogorov (1941), which forms the basis for theoretical spectra using the eddy excitation mechanism. This might explain why our simulations show a different slope in the frequency spectra.

The deviation from theoretically predicated slopes might be due to the relatively low Reynolds numbers reached in the simulations (see Tab. 1). On the other hand, the case of heating concentrated in a small spherical region is quite different from the plane-parallel, Boussinesq convection underlying some theoretical models and the velocity field is not necessarily isotropic in this case.

Comparing with previous hydrodynamic simulations we find that Model H6R10 agrees well with a comparable model from Rogers et al. (2013), who find a broken power law fit with exponents of $-4.8$ and $-1.9$ and a break at $l\approx 10$ in a singular-value decomposition of the frequency and wavenumber spectrum. Augustson et al. (2016) obtain a qualitatively similar spectrum in their simulations of magnetic turbulence, with a low-wavenumber exponent of approximately $-3$ and a steeper power law for higher wavenumbers. Their simulations also have a peak in the kinetic energy spectrum at low spherical harmonic degree, in the range from $l=1$ to $l=10$. We note that the models of Augustson et al. (2016) do not have enhanced forcing and yet show a similar spectrum to those in this work. This indicates that the spectrum is more dependent on the regime that nondimensional numbers like $\mathrm{Ra}$ and $\mathrm{Re}$ are in than the actual value of convective forcing, as expected.

The inset in Fig. 9 shows the comparison of the high-wavenumber tail for the simulations H7E and H7E-HR, which are run with identical parameters except for the number of spectral modes being used. Their spectra are almost identical apart from a small bend at the highest $l$ values. This suggests that enough of the inertial subrange of the turbulent cascade is resolved to get the correct energy dissipation and that the simulations do not suffer from severe anomalous behavior at the smallest resolved length scales. The same is true for H6LD and its high-resolution counterpart H6LD-HR.

As a main concern of this paper are the IGW spectra, we also assess the impact of angular resolution on these. Figure 10 shows the frequency spectrum of radial velocity in the radiation zone for the two simulations H6LD and H6LD-HR, where both are identical except for latter having twice the angular resolution. We see that both simulations are very similar, including the magnitude and position of the modes between $200\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ and $300\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$, the continuous signal between $20\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ and $200\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$, and the low frequency drop due to radiative damping at $20\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$.

4.2 IGW generation

It is controversial which physical mechanism is most important for the excitation of IGWs at the CB. The two common candidates are bulk Reynolds stresses produced by convective eddies (Lighthill, 1952; Goldreich & Kumar, 1990) and plume overshoot (Townsend, 1966; Zahn, 1991). Most theories about the effect of IGWs in stellar interiors (e.g., Talon & Charbonnel, 2005; Fuller et al., 2014) employ the spectrum of IGWs derived from convective eddies (Kumar et al., 1999; Lecoanet & Quataert, 2013). Therefore our analysis focuses on this spectrum, but the plume spectrum is considered later.

To study the spectrum of waves generated, we first investigate the spectrum of motions generated at the CB. We analyze our 3D data by computing the spectrum of kinetic energy density at a radius of 0.07 $H_{P}$555The pressure scale height is defined as $H_{P}=-\partial r/\partial\ln P$. above the top of the convection zone (as defined by the Schwarzschild criterion). This spectrum is given by,

[TABLE]

with the Fourier transforms of the individual velocity components, $\hat{v}_{r}$, $\hat{v}_{\theta}$, $\hat{v}_{\phi}$, according to Eq. (13). Figure 11 shows the spectra for different models. For guidance we show an estimate of the convective turnover frequency given by,

[TABLE]

which assumes the largest eddy extends from the center of the star to the radius of the convection zone $r_{\text{CZ}}$ and it turns at the rms velocity. Panel (b) shows the spectra multiplied by $f$ to account for integration over $d\log f$, which makes it easier to see the regions containing most energy in the logarithmic plot. We see that, while the peak is not too far from $f_{\text{TO}}$, the distribution is almost flat in the low frequency regime.

Clearly, this is a spectrum of motions at this radius and is not necessarily waves (although see Sect. 4.5). However, this motion is what drives the waves and if it has a high-frequency component then high-frequency waves can be efficiently driven.

At this radius the integrated (i.e. including all harmonic degrees $l$) frequency spectrum is nearly flat with a transition to a more steeply declining power law at higher frequencies ($f\gtrsim$20\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$$). The spectrum is not dominated by values at $f_{\text{TO}}$ and indeed it is hard to make out this frequency in the spectrum. However, if we look at the frequency spectrum at particular length scales, by selecting individual values of $l$, we start to see a sharp transition between the power laws at low and high frequency as evidenced in Fig. 12.

In this scale-dependent spectrum the break point between the two power laws depends mostly linearly on angular degree $l$ and can be approximated with,

[TABLE]

The slope is not too far from the estimate for the convective turnover frequency $f_{\text{TO}}=$7.2\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$$ for this model, considering the uncertainty in the estimate of $f_{\text{TO}}$ in Eq. (17). This fits the conjecture by Rogers et al. (2013) that the eddy mechanism efficiently generates waves below this frequency.

It is worth noting that, even in this scale-dependent spectrum in which the break between power laws corresponds to the scale-dependent turnover frequency, the energy is not concentrated at that frequency. This is in stark contrast to the theoretical predictions which posit that the frequency spectrum (within the CZ) is strongly peaked at the convective turnover frequency.

In a more systematic analysis of the broken power law fits to the frequency spectrum at the CB in Fig. 13, we notice that the exponent of the low frequency regime stays relatively constant for $l>3$ at a value of -0.46$\pm$0.07. The high frequency component covers a wider range of exponents from $-3.6$ at $l=2$ to $-6.6$ at $l=33$. Both exponents show very little change at higher values of $l$. The position of the frequency break point $s$ in the power law in Fig. 14 on the other hand is rising with $l$ roughly following the estimate for the convective turnover frequency from Eq. (17) multiplied by $l$ to account for the smaller length scales at higher spherical harmonic degree. For $l\gtrsim 31$ we observe a rise in the exponent of the high frequency range, which is due to the difficulty in fitting an increasingly smaller part of the curve. This makes the determination of $s$ less certain as well. The value of $s$ lies within the range between 12 and $200\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ in our simulations, which is a bit higher than the range of 10 to $80\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ in the 2D simulations of Rogers et al. (2013). This is understandable if the position of the break really depends on $f_{\text{TO}}$, and in turn $v_{\text{rms}}$, because the 2D simulations show a lower convective velocity.

The logarithmic scaling of Fig. 11 makes it hard to see, which frequencies contribute most to kinetic energy at the top of the convection zone. To analyze this we plot the cumulative energy distribution, i.e. the function of energy contained below a certain frequency. It is shown normalized to the full energy of the particular model in Fig. 15. We see that for models with a heating rate increased by a factor of $10^{6}$, roughly 50% of the energy is below $f_{\text{TO}}$ and 50% is above. While the steepest increase is around $f_{\text{TO}}$, the distribution of energy is widely spread in frequency. The shaded area shows the frequency range around $f_{\text{TO}}$ that contains 40% of energy. In model H7E with a heating rate of $10^{7}$ times the stellar value, 80% of the energy is located below $f_{\text{TO}}$, while model $H5$ has almost all energy far above its value of $f_{\text{TO}}$ at $0.5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$. As discussed earlier we believe model H6LD is the best trade-off between an increased heating rate and increased diffusivity.

Comparing these numerical spectra to theoretical spectra is not straightforward as clearly the former include wave motion as well as overshooting motion (although see Sect. 4.5). However, one can trace differences between theoretical IGW spectra to differences in the assumed convective spectra. While the theoretical wavenumber spectrum (following Kolmogorov, 1941) has some observational basis, the frequency spectra supposed in the theoretical analysis of Kumar et al. (1999) and Lecoanet & Quataert (2013) – based on the assumption of Kolmogorov scaling of eddy sizes and their corresponding turnover times – do not. Yet it is this frequency spectrum in the CZ that determines the frequency spectra of excited IGWs.

For example, theoretical spectra do not efficiently generate high-frequency waves because of the assumption that most of the convective energy is concentrated at the convective turnover frequency. If the energy of convection itself is not limited to a narrow band around the convective turnover frequency, there is no reason to suppose that the IGW frequency spectrum would be. Moreover, if the CZ has high frequencies then it can efficiently generate waves of high frequency. Therefore, based on comparisons of CZ spectra one can see two important issues arise between theoretical and numerical results that would affect the IGW spectra: (1) for an integrated spectrum, energy is not concentrated at the convective turnover frequency, but is spread among a wide range of frequencies; and (2) frequencies higher than $f_{\mathrm{TO}}$ are clearly present with significant energy within the convection zone. We also note that while Kolmogorov scaling may be the appropriate description for isotropic turbulence in a Boussinesq box, it is wholly unclear that it is appropriate for spherical configurations with a centrally peaked heating term, such as stars.

Kumar et al. (1999) mention that they deliberately ignore wave excitation by plumes due to limited information on their properties. As the flow pattern we observe in the simulations is obviously dominated by large plumes, it is natural to compare the spectra to theory of IGW excitation by plume penetration. Montalbán & Schatzman (2000) developed expressions for the IGW spectrum generated by plume penetration at the bottom of the solar convection zone. It is based on the plume model by Rieutord & Zahn (1995). They explicitly caution against its use at the top of a convective envelope because of the typical importance of radiative cooling there, characterized by a low Péclet number ($\mathrm{Pe}\approx 1$). This argument does not apply here, where $\mathrm{Pe}\gtrsim${10}^{4}. For comparison, the stellar value in the core is $\mathrm{Pe}\gtrsim${10}^{6}. Therefore, for lack of a dedicated theory, we apply their model to our spectra.

The frequency dependence of the energy spectrum is determined by the plume timescale $t_{\text{b}}$, which is often approximated by the ratio of plume velocity $v_{\text{pl}}$ and plume incursion depth $\Delta_{\text{p}}$. The corresponding frequency is then,

[TABLE]

The predicted energy spectrum takes the form (Montalbán & Schatzman, 2000),

[TABLE]

where we absorbed all factors depending on radius and wavenumber into the proportionality constant. We choose to work directly with the expression for kinetic energy instead of wave flux because the spectrum is taken right at the top of the CZ, where a conversion to flux is not straightforward.

A combination of plumes with different timescales can be fit to the simulation spectra. Figure 16 shows fits with one and three values of $f_{\text{b}}$. This shows that the plume spectrum as described by Eq. (20) generally fits the shallow power law in the low-frequency regime very well, even with a single plume frequency (red line). With just three plume timescales (orange line) it is possible to fit most of the spectrum. Assuming plumes are distributed so that their frequencies follow an exponential function with a low frequency cut-off, allows us to fit the whole spectrum apart from the high frequency turnoff (cyan line). This heuristic expression has the form,

[TABLE]

where the shape of the exponential is given by $\alpha=$11.0\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z} and the low-frequency cut-off is $c=$1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}. It should be noted that similar fits can be obtained with other steeply declining functions (e.g., a power law with a negative exponent) for the plume frequency distribution.

To acquire an estimate for typical plume length and time scales from our simulations, we study the process of plume incursion in more detail. Figure 17 shows a Cartesian projection of temperature perturbation at the CB in two different cases. In the left panel a large plume hits the boundary exciting waves at a large range of phase angles, including very steep angles. A representation of excitation through eddies is shown in the right panel. It results in much smaller phase angles. Both cases are remarkably similar to 2D simulations (Rogers et al., 2013, Fig. 4).

The large plume of Fig. 17 is studied in greater detail in Fig. 18 to extract its size, velocity, and penetration depth. The Brunt–Väisälä frequency $N^{2}$ (Panel (d)) is significantly reduced in the overshooting region above the original convective boundary. This coincides almost perfectly with the penetration depth of the plume, which can be identified by $v_{r}$ approaching 0 (Panel (b)) and a discontinuity in $T$ (Panel (c)). The penetration depth is $\Delta_{\text{p}}=0.4\,H_{P}$ from the original convective boundary. The maximum plume velocity is $v_{\text{pl}}=$7.8\text{,}\mathrm{k}\mathrm{m},\mathrm{s}^{-1}. This allows us to estimate a plume incursion time $t_{\text{b}}\sim\frac{\Delta_{\text{p}}}{v_{\text{pl}}}=$2.2\text{\,}\mathrm{h}, which corresponds to a frequency $f_{\text{b}}$ of $128\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$. The lateral extent of the plume $b$ can be defined as the region of positive $v_{r}$ in Panel (e). It has a value of $b=$1.4\text{\times}{10}^{10}\text{,}\mathrm{c}\mathrm{m}$$ in this case. This is also sometimes used to compute the plume timescale, which results in a value of $56\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ here. As it is easier to extract from simulations in a systematic way, we stay with the first definition using $\Delta_{\text{p}}$ in the following analysis.

We apply this estimate of the plume frequency statistically to all plumes in simulation H6LD. For each longitude and latitude, and each output snapshot (every $1000\text{\,}\mathrm{s}$) we determine if there is a rising plume and in that case compute a plume incursion depth and plume velocity. The criterion for a plume is that $v_{r}$ is positive at the position of the convective boundary. The incursion radius $r_{\text{p}}$ is then defined as the radius at which $v_{r}$ first becomes negative along a line at this particular angle. The penetration depth is calculated as the distance to the convective boundary $\Delta_{\text{p}}=r_{\text{p}}-r_{\text{conv}}$. The plume velocity $v_{\text{pl}}$ is the highest value of $v_{r}$ between the top of the convective boundary and $r_{p}$. Using the estimate for the plume frequency from Eq. (19), we compute the probability density function (PDF) of $f_{\text{b}}$ throughout the simulation (Fig. 19). It rises sharply to its maximum at $15\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ and then drops roughly following an exponential distribution. An exponential fit to the data does not perfectly match the values found when fitting Eq. (21). The parameter $\alpha$ is too high by a factor of 3. Yet considering the simplistic definition of $f_{\text{b}}$, this still makes a strong argument for an exponential distribution of plume frequencies as the explanation of a large part of the kinetic energy spectrum at the top of the convection zone and hence, the IGW frequency spectrum.

To understand the effect of increased forcing and diffusivity we follow the discussion of plume lifetimes of Pinçon et al. (2016). They argue that plume velocity scales with luminosity as $v_{\text{pl}}\propto L^{1/3}$ which is consistent with the scaling of the convective velocities from Eq. 14. A luminosity increased by a factor of $10^{6}$ would thus result in $v_{\text{pl}}$ increased by a factor of 100. The penetration depth $\Delta_{\text{p}}$ is not expected to be strongly affected by the change in forcing (see end of Sect. 4.3 for an estimate). The effect of radiative thermalization, while strongly increased due to the higher value of $\kappa$ in the simulations, is still negligible as the timescale $t_{\text{rad}}\sim\Delta_{\text{p}}/\kappa$ is of the order of two years, much longer than any observed plume lifetime. The turbulent timescale inside the plume $t_{\text{turb}}\sim b/v_{\text{pl}}$, with the lateral plume size $b$. Assuming that $b$ is not strongly affected by increased forcing, similar to $\Delta_{\text{p}}$, this means that plume frequency $f_{\text{b}}=1/t_{\text{b}}$ scales like $v_{\text{pl}}\propto L^{1/3}$. In model H6LD ($L=10^{6}L_{\star}$) this results in $f_{\text{b}}$ being too high by a factor of 100.

4.3 Convective overshoot

The treatment of convective–radiative boundaries (CB) in 1D stellar evolution codes is a long-standing problem. It can have a significant impact on the evolution and nucleosynthetic signature of stars by mixing of species beyond convective regions. Hydrodynamic simulations in two or three dimensions promise insight based on first principles and have been subject of previous work (e.g., Freytag et al., 1996; Rogers et al., 2006; Meakin & Arnett, 2007; Jones et al., 2017; Cristini et al., 2017). There is no single accepted definition of the overshooting depth in terms of angular averages of 3D quantities. For better comparability between different stellar parameters the overshooting depth is usually stated in multiples of the pressure scale height $H_{P}$ above the convective boundary as defined by the Schwarzschild or Ledoux criterion. Both criteria are equivalent in the case studied here because the star is chemically homogeneous.

The statistics of velocities and penetration depth $\Delta_{\text{p}}$ from Sect. 4.2 can also be used to make statements on the size of the overshooting region. Figure 20 shows the PDF of plume penetration depth in model H6LD. The distribution at low $\Delta_{\text{p}}$ is relatively flat until it peaks at $0.54\,H_{P}$. Beyond that it drops following a power law with exponent $-1.4$. This is consistent with the picture in Panel (d) of Fig. 18, where $N^{2}$ is affected by penetration up to a value of approximately $0.5\,H_{P}$. 95% of plumes penetrate no further than $0.695\,H_{P}$, which is the value we will use as the boundary of the overshooting region in Sect. 4.5.

Figure 21 shows the distribution of updrafts and downdrafts in the overshooting region ($0.42\,H_{P}$ above $r_{\text{conv}}$). The PDF is peaked in Lorentzian shape at $-0.2\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$. The inward velocities are distributed in a smaller range, $-10\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$ at most, than the outward velocities, which extend up to $25\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$.

Our use of an increased convective forcing and thermal diffusivity raises the question of the validity of these results for the actual stellar values. In his study of convective penetration in stellar interiors Zahn (1991) found a simple scaling law for the size of the penetrative region (see also, Rogers et al., 2006),

[TABLE]

with the scale height of thermal diffusivity $H_{\kappa}=-d\ln r/d\ln\kappa$, plume filling factor $f$, and total energy flux $F_{\text{tot}}$. As $\kappa$ is only multiplied by a constant in the radiation zone, $H_{\kappa}$ is identical to the stellar value. The same is true for $H_{P}$ and $\rho$. While the simulations have an increased $F_{\text{tot}}$, we found the scaling $F_{\text{tot}}\propto v^{3}$ (Eq. (14) and Fig. 5). This means the penetration depth in the simulations and in the star only vary by a factor of $\sqrt{f_{\text{sim}}/f_{\star}}$, which we expect to be a number not too far from unity.

4.4 IGW propagation

Convective motions in the core generate IGWs at the CB, which propagate through the cavity of positive $N^{2}$ in the radiation zone. To visualize the excited frequencies and the change of the wave spectrum with radius, we compute the frequency spectrum of $v_{r}$ sampled at several longitudes around the stellar equator at all times for all radii and show it as a heat map in Fig. 22 for model H6LD. In the convection zone ($r\lesssim$2\text{\times}{10}^{10}\text{,}\mathrm{c}\mathrm{m}$=0.14\,R_{\star}$) we note the presence of all frequencies with a clear dominance of the range below $50\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$. This is reflective of the large range of timescales of convective motion (see Sect. 4.2). In the radiation zone frequencies up to $300\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ are excited. Low-frequency waves are strongly damped and only frequencies above $40\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ reach the top of the simulation domain. This is qualitatively in agreement with linear theory, which predicts that lower frequency waves experience stronger damping (e.g., Kumar et al., 1999). The exact position of this lower cut-off depends on the value of $\kappa$ as well, which is why the present simulations cannot predict it quantitatively. Figure 23 shows line plots of the same spectrum at different radii. It shows that the spectrum in the RZ at low frequencies ($\lesssim$20\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$$) does not reach the numerical noise level, as would be expected by the excessive numerical diffusion at this frequency, but turns flat at a higher value. The figure also indicates the expected frequency dependence of $v_{r}$ from theoretical work by Lecoanet & Quataert (2013), which is $f^{-3.25}$ for the radial velocity of waves excited at a discontinuous $N$ profile. We see that the simulated spectrum is much flatter than this prediction, following $f^{0.8}$. The steep drop around $200\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ for $r>0.8\,R_{\star}$ is due to the limit imposed by the Brunt–Väisälä frequency at these radii.

Strong vertical features are visible in Fig. 22. These are peaks in the spectrum which are present at the same frequency at all radii in the radiation zone. This identifies them as standing waves. Their frequencies are determined by the cavity they resonate in and can be computed numerically using the stellar oscillation code GYRE (Townsend & Teitler, 2013).

It is hard to disentangle individual modes because the contributions of several wave numbers overlap, but due to the horizontal discretization of the simulations using spherical harmonics it is simple to extract a spectrum for particular $l$ and $m$ modes. The panels in Fig. 24 show the frequency spectrum for the modes $l=2$ and $l=4$. Here, the radial order of the individual standing modes can be clearly identified by the number of radial nodes. The strong mode at $210\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ without any nodes is a fundamental mode or f mode. The other visible modes show an increasing number of nodes with decreasing frequency. This identifies them as g modes (e.g., Aerts et al., 2010, Sect. 3.5). We computed expected mode frequencies with GYRE666We used version 5.1 from the GYRE web page. for comparison. They are labeled in the figure using the Eckart-Osaki-Scuflaire-Takata scheme (e.g., Aerts et al., 2010), where negative numbers indicate g modes, 0 is the f mode, and positive numbers are p modes. In the case of $l=2$ we find quite good agreement for the g modes (at least up to $\text{g}_{4}$) and the f mode, especially considering that our 3D simulation has a slightly different resonant cavity due to the different equation of state and outer boundary compared to the 1D MESA model. As expected there are no p modes as the chosen set of equations (see Sect. 2) does not include the physics of sound waves. For $l=4$ (lower panel of Fig. 24) the $\text{g}_{1}$, $\text{g}_{2}$, and $\text{g}_{3}$ modes match very well, while the identified f mode is within a few $\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ of the expected frequency of the $\text{p}_{1}$ mode according to GYRE. This is probably coincidental as the discrepancy between 3D hydrodynamics and GYRE gets even larger at higher wave numbers. The identification of modes by counting the number of nodes in the RZ is illustrated in Figure 25, which shows the radial change of $v_{r}$ amplitude of particular frequency components corresponding to the standing waves. The amplitudes were computed by projection on a complex phase angle of the Fourier transform at the radius with the maximum absolute value in the RZ. The $\mathrm{g_{2}}$ and $\mathrm{g_{3}}$ show nodes at the top of the overshooting region, which is ignored for the mode identification.

Linear theory predicts amplification of waves moving along a decreasing density profile through pseudomomentum conservation (e.g., Bühler, 2009). At the same time thermal diffusivity damps the wave. Ratnasingam et al. (2019) give an expression for the linear wave amplitude based on Press (1981) and Kumar et al. (1999). The amplitude of the radial velocity follows

[TABLE]

with

[TABLE]

using $\omega=2\pi f$ and the starting radius of wave propagation $r_{0}$ with its corresponding density $\rho_{0}$ and Brunt–Väisälä frequency $N_{0}^{2}$. We extract the amplitude of $v_{r}$ at several frequencies for a particular $l$ mode and show it together with the theoretical prediction from Eq. (23) in Fig. 26. We see that the waves generally follow amplification through the $\sqrt{\rho_{0}/\rho}$ term and are hardly affected by radiative damping, except for the low frequency case, as expected.

Generally the match between the GYRE predictions and data extracted from 3D hydrodynamics is quite promising, considering that both approaches make slightly different assumptions about the physics. Even with the high thermal diffusivity needed for the simulation we can see wave amplification. We might be able to observe wave breaking in future simulations which extend to regions closer to the surface at much lower densities.

4.5 Nature of the signal in the radiation zone

Although a visualization of the temperature field in the radiation zone such as in Fig. 8 suggests a wave nature of the flow field, a more rigorous analysis is needed to prove the motions are indeed IGWs excited close to the convective boundary or by nonlinear interaction in the RZ. IGWs are naturally limited to frequencies below the Brunt–Väisälä frequency, i.e. $\omega<N$ with $\omega=2\pi f$. In Fig. 27 we show the high-frequency part of the spectrum of $v_{r}$ at all radii. The white line in the figure indicates the local Brunt–Väisälä frequency. We see that the bulk of the signal in the RZ is constrained to the region $\omega<N$. Beyond this frequency there is a sharp drop in the amplitude which is consistent with IGW nature. A notable exception are the strong f modes (e.g. at $320330340\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$) going beyond that limit, which does not contradict this interpretation because these modes are not subject to the frequency limit.

In the locally Boussinesq but globally anelastic approximation IGWs follow the dispersion relation (e.g., Press, 1981),

[TABLE]

Here, $k_{\perp}$ is the horizontal wave number, $k_{r}$ the radial wave number, and $k=\sqrt{k_{\perp}^{2}+k_{r}^{2}}$ the magnitude of the total wave vector. We verify this relation for individual values of angular degree $l$ and frequencies because the resulting velocity field is a combination of many individual waves. The horizontal wave vector can easily be computed for a given $l$ by

[TABLE]

The radial wave number is not as straightforward to derive because the wave length changes with radius as the Brunt–Väisälä frequency varies. An additional complication is that the wavelength becomes comparable to the stellar radius above $r\sim 0.5\,R_{\star}$, which makes an accurate determination very hard.

We determine the radial wavelength $\lambda_{r}$ for each individual frequency by finding the peaks of $v_{r}$ along a ray in the radial direction and calculating the difference between them. For this we employ the routine signal.find_peaks from the scipy Python package (Jones et al., 2001–), which finds isolated local extrema and is resilient to small numerical noise. As there are only very few wave cycles along the total radius of the star, we use a cubic spline to interpolate the wavelengths at every radial coordinate. Several other methods proved unsuccessful in this particular case: using a radial Fourier transform with a sliding window is inaccurate, as there are only few wave cycles per window; calculating the radial derivative of the phase of the spectrum works well except for regions where the phase is poorly defined when the amplitude is close to 0. This makes this method inapplicable to determining the wavelength in standing modes. The simple method of measuring the distance between peaks and interpolating the found wavelengths is the most robust.

Using this method we calculate the radial wave number $k_{r}=2\pi/\lambda_{r}$ for every frequency at every radius to check how closely Eq. (25) is fulfilled. The two panels in Fig. 28 show this for $l=2,4,10,20$. White regions indicate a match of the dispersion relation, red regions have a too large $k_{\perp}$, blue regions have a too small $k_{\perp}$. All components show no match in the very low frequency range ($\lesssim$10\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$$), which is expected to be totally dominated by diffusion at all radii. Just above the overshooting region at $0.2\,R_{\star}$ we find excellent agreement at higher frequencies, which we interpret as waves being emitted from the convection zone over a large range of frequencies. Due to the increased thermal diffusivity needed for numerical reasons, low-frequency waves cannot propagate far into the RZ. This is evidenced by the increasing size of the non-IGW (red) region at the low-frequency end.

The higher $l$ values show a remarkable phenomenon. At $r\gtrsim 0.4\,R_{\star}$ in the low-frequency region which should be completely dominated by damping ($f\lesssim$70\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$$ for $l=10$) a signal appears which matches the dispersion relation. This cannot be explained by waves originating from the convective boundary because there are no waves of these frequencies present at lower radii. A plausible explanation is that these are generated by nonlinear interaction of low $l$ waves in the middle of the RZ. These secondary waves reach frequencies from $10\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$.

The ability of the discretization to resolve IGWs is checked in this context as well. For given values of $k_{\perp}$ and $N$ we can calculate a frequency below which the radial IGW wavelength would not be resolved by at least 10 grid points in the radial direction. This frequency forms the lower limit for resolving IGWs at a given radius in this simulation. The limit is indicated as a magenta line in Fig. 28. The lower limit is highest close to the convective boundary, where $N$ is lowest.

Another limit on wave resolution is imposed by diffusion. As a rough estimate for the minimum wave length of waves not dissipated by diffusion and viscosity we use,

[TABLE]

with the IGW wavelength $\lambda=2\pi/k$ and the diffusion time $\tau_{\mathrm{d}}=\lambda_{\perp}/v_{\mathrm{g}}$. This uses the magnitude of the group velocity,

[TABLE]

Solving this equation for $f$ leads to the yellow colored line in Fig. 28. As expected, we can see that motions below this frequency largely do not fulfill the dispersion relation. This proves that the radial resolution in our simulations is sufficient to resolve waves with frequencies above $30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ in the most energy bearing wave numbers ($l\lesssim 5$).

While this analysis gives a good overview of the frequencies, wave numbers, and radii where the dispersion relation is fulfilled well, it is also important to see what fraction of the kinetic energy is actually contained in the wave motions. To compute this we filter the Fourier transform of the velocity components $\hat{v}_{r}$, $\hat{v}_{\theta}$, and $\hat{v}_{\phi}$ to only include values at frequencies and radii, where $k_{\perp}/k$ is close to the IGW dispersion relation $\omega/N$. The filtered velocities are,

[TABLE]

The kinetic energy computed from these velocities is then identified as the energy in IGW motions $E_{\mathrm{IGW}}$. Figure 29 shows the ratio of this energy to the unfiltered kinetic energy integrated over angular degrees $l\leq 10$. We do not include frequencies below the limit definitely dominated by diffusion (Eq. (27)) in this analysis. As expected, we see almost no energy is in IGWs from the center up to the top of the overshooting region at $r=0.2\,R_{\star}$, which matches the previously determined position from Sect. 4.2 as indicated by the vertical dotted line. In the RZ it rises to 90% when applying the error margin $C_{\mathrm{cut}}=1.3$. The fraction of kinetic energy in IGWs drops beyond $r=0.4\,R_{\star}$. The main cause of this is the uncertainty in determining $\lambda_{r}$ at large radii, where $\lambda_{r}$ approaches $R_{\star}$ and our method of measuring the distance between peaks breaks down. Another reason is the growth of the low-frequency, red regions in Fig. 28, which arises because of the limited range of IGWs due to high numerical diffusion. At even higher radii, $r\gtrsim 0.7\,R_{\star}$, we notice an increase in the IGW energy fraction, which is likely due to the increased fraction of secondary waves, as discussed earlier.

This analysis makes us confident that the motions in the RZ are indeed of IGW nature to a significant fraction, with the exception of very low frequency motions dominated by numerical diffusion. We see evidence for secondary generation of waves within the RZ.

4.6 IGW surface signature

The observed brightness variations in O stars have been suggested as signatures of convectively excited IGWs (Aerts & Rogers, 2015). Their spectrum is likely linked to that of tangential velocity close to the surface of the star (De Cat & Aerts, 2002; Tkachenko et al., 2014). Figure 30 shows a spectrum of latitudinal velocities from model H6LD (orange line). These are less affected by numerical influence from the boundary condition than the azimuthal velocities. This is compared to a spectrum obtained from 2D simulations of a $3\,\mathrm{M}_{\odot}$ star from Rogers et al. (2013) (blue line). We see the same low-frequency power excess in range from $2\text{\,}\mathrm{d}^{-1}6\text{\,}\mathrm{d}^{-1}$ and a similar drop in amplitude below $2\text{\,}\mathrm{d}^{-1}$ as in the 3D simulations. We also plot photometric observations (Blomme et al., 2011) with an amplitude ratio of $1\text{\,}\mathrm{m}\mathrm{m}\mathrm{a}\mathrm{g}\,\mathrm{(}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}\mathrm{)}^{-1}$ (De Cat & Aerts, 2002; Aerts & Rogers, 2015) for comparison.

The same low-frequency power excess shows in the spectra of temperature fluctuations in Fig. 31. It is expected to be the dominant cause of photometric variability in observations. In the spectrum integrated over all $l$ values (top panel) this excess makes it hard to distinguish individual excited mode frequencies, expect for one mode at $210\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ ($18.1\text{\,}\mathrm{d}^{-1}$), which is part of the $l=2$ component. In the spectra for individual $l$ values, several modes can be identified, corresponding to those in Fig. 24. This decomposition also makes it clear that the low-frequency power excess is a combination of the power excesses in different $l$ components, each contributing to a small frequency range. The lack of signal at low-frequencies in numerical simulations is due to the high numerical diffusivity and not expected to be physical. The amplitude of the waves is expected to increase as they propagate from $r=0.89\;R_{\star}$, where the spectrum was computed, to the surface. According to pseudomomentum conservation it should increase by a factor of 380.

The low-frequency power excess is also found in observations of stars with a convective core (Bowman et al., 2019). In the simulations it is caused by the high density of high radial order, low-$l$ g modes and because most energy in the CZ is at low $l$ values (see Fig. 9). The drop in amplitude below $2\text{\,}\mathrm{d}^{-1}$ on the other hand is in disagreement with observed photometry (Blomme et al., 2011; Aerts et al., 2017b, 2018; Bowman et al., 2019). This disagreement is likely caused by the increased thermal diffusivity in the simulations (both in 2D and 3D) which damps low-frequency waves more strongly than in stellar interiors. Another possibility is the lack of differential rotation in our 3D model. Rogers et al. (2013) found that differential rotation between core and envelope introduces a significant low frequency component in the spectrum. At higher frequencies above $10\text{\,}\mathrm{d}^{-1}$ the 2D simulations drop more slowly than the 3D simulations and show many excited modes. This is possibly due to the lack of wave breaking brought about by the high thermal diffusivity and viscosity needed in our present set of 3D simulations.

5 Conclusions and Outlook

We showed the first 3D simulations of convection in the core of an intermediate-mass star, with a convective core and radiative envelope, that also include a large part of the radiation zone (RZ). The simulations using the anelastic equations (i.e., removing the physics of sound waves) and a spectral discretization using spherical harmonics were run using a realistic reference state from the stellar evolution code MESA. For numerical reasons the simulations were run with increased thermal and viscous diffusivity. To compensate for the increased wave damping this produces we increased the luminosity of the star causing higher velocities in the convective core. We do this in the hope that wave velocities at the surface of the star are more realistic.

We see wave patterns in the RZ, which are identified to be standing g and f modes with frequencies similar to those predicted by the oscillation code GYRE. Although there are differences, they are not of the sort predicted in Brown et al. (2012) and are dependent on the $l$ and $m$ values of the spherical harmonics. These differences are likely due to slightly different physics (e.g., equation of state, outer boundary condition) and changes in the temperature profile at the top of the convection zone due to overshooting.

Apart from the standing modes the simulations also show a continuous signal in the RZ between frequencies of approximately $20\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$ and $200\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$. An analysis of the dispersion relation (see Fig. 28) identifies the physical mechanism as IGWs. The decline of this continuous spectrum with frequency is markedly smaller than theoretically predicted values for excitation purely due to convective eddies (see Fig. 23).

An analysis of the kinetic energy distribution over spherical harmonic degree $l$ shows a spectrum which peaks at a low value of $l$ and then declines with a power law with an exponent in the range from $-2-3$ in the inertial range. This is closer to the theoretical value for Bolgiano–Obukhov scaling ($-2.2$) of buoyancy dominated convection than to the Kolmogorov value ($-1.6$) of isotropic turbulence. We do not have enough information to get a conclusive answer on the realized scaling in convective stellar cores. This should be studied further in detailed simulations of just the core. These slopes are measured at the top of the CZ, which is subject to convective overshooting, and therefore do not directly match the IGW spectrum. Yet they show what energy is available for wave excitation at a given frequency. The slope in the inertial range is similar to that observed in other 2D or 3D simulations of core convection (Rogers et al., 2013; Augustson et al., 2016).

The broken power law structure of the frequency spectra of kinetic energy above the convection zone is similar to those obtained in the 2D simulations of Rogers et al. (2013), suggesting the mechanism driving the bulk of this spectrum does not fundamentally change with dimensionality. It is likely that bulk Reynolds stresses induced by convective eddies contribute more in higher Reynolds number flows, but this would still only affect the low frequencies ($f<f_{\text{TO}}$) and hence, have little impact on angular momentum transport or mixing within the bulk of the RZ (Shiode et al., 2013; Kumar et al., 1999; Lecoanet & Quataert, 2013).

Excitation by plume penetration is obviously involved as can be seen in the temperature and velocity fields (see Fig. 7). It can explain the excitation of higher frequency waves and the extracted distribution of plume frequencies fits a large part of the simulation spectrum. One may argue that the plume penetration depths, and hence, frequencies generated are too large. However, at least in the theory by Zahn (1991), this penetration depth scales like the velocities cubed divided by the total flux, a number which is the same in the simulations as it is in the star. The production of high-frequency waves is extremely important for explaining the photometrically observed brightness variations at high frequencies (see Fig. 30). They are likely underestimated in our simulations due to high dissipation preventing wave breaking.

Stochastic brightness variations caused by velocity and temperature fluctuations at the stellar surface have been inferred to be caused by IGWs in massive stars (Aerts & Rogers, 2015; Aerts et al., 2017a, 2018; Bowman et al., 2019). We extracted frequency spectra of these quantities from the simulation close to the stellar surface. General features are a low-frequency power excess and the presence of standing modes at low $l$ harmonics. This is in agreement with the findings of 2D simulations (Rogers et al., 2013), which match observations in the power bearing range, but lack both amplitudes at very low frequencies (due to excessive radiative damping) and high frequencies (possibly due to lack of wave breaking from overdamped waves).

The simulations presented in this article show the feasibility of hydrodynamic modeling of convectively excited IGWs and their propagation through a large part of the radiative zone using a consistent numerical treatment. In future work employing more computational resources the limitations forcing us to use unphysically high diffusivities and luminosities can hopefully be overcome to achieve more realistic wave amplitudes throughout the interior and at the surface, and hence more realistic angular momentum transport. More realistic physical parameters in the simulations combined with coverage of a wider range of stellar models will also allow us to make quantitative predictions of the expected signature of IGWs in asteroseismological observations.

The general similarity of our results with those of previous 2D simulations encourage us to consider those results with less reservations due to their dimensionality and to use 2D simulations as a tool in the future to quickly cover a wider parameter range in models than is possible with 3D simulations.

Support for this research was provided by STFC grant ST/L005549/1 and NASA grant NNX17AB92G. MGP and DMB received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No670519: MAMSIE). VP acknowledges support from the European Research Council through ERC grant SPIRE 647383. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. This research made use of the Rocket High Performance Computing service at Newcastle University. The authors thank C. Pinçon and M. Rieutord for helpful comments.