Synergies between Asteroseismology and Three-dimensional Simulations of Stellar Turbulence

TL;DR

This paper compares asteroseismically inferred stellar interior structures with predictions from 3D turbulence simulations to improve understanding of stellar convection and reduce reliance on adjustable parameters in traditional models.

Contribution

It demonstrates promising agreement between 3D turbulence simulation predictions and asteroseismic observations of B stars, advancing the integration of simulations into stellar modeling.

Findings

Good match between simulation and asteroseismic composition profiles

Insights into boundary layer structure and wave generation

Potential to improve stellar evolution models

Abstract

Turbulent mixing of chemical elements by convection has fundamental effects on the evolution of stars. The standard algorithm at present, mixing-length theory (MLT), is intrinsically local, and must be supplemented by extensions with adjustable parameters. As a step toward reducing this arbitrariness, we compare asteroseismically inferred internal structures of two Kepler slowly pulsating B stars (SPB's; $M\sim 3.25 M_\odot$) to predictions of 321D turbulence theory, based upon well-resolved, truly turbulent three-dimensional simulations (Arnett , et al. 2015, Christini, et al. 2016) which include boundary physics absent from MLT. We find promising agreement between the steepness and shapes of the theoretically-predicted composition profile outside the convective region in 3D simulations and in asteroseismically constrained composition profiles in the best 1D models of the two SPBs. The structure and motion of the boundary layer, and the generation of waves, are discussed.

Figures

Tables

Table 1. Table 1 Comparison of composition gradients ( θ , ϕ 𝜃 italic-ϕ \theta,\phi and η 𝜂 \eta in Eq.  6 ) to H profile of two Kepler SPBs, and to C and O profiles from 3D simulations (resolutions ∼ 1024 3 similar-to absent superscript 1024 3 \sim 1024^{3} ).

Source	Burning	$\theta$	$\varphi$	$\eta$
Asteroseismology:
KIC 10526294	core H	0.63	0.71	112
KIC 7760680	core H	0.50	0.71	57
3D simulations:
upper boundary	shell C	0.727	0.828	184
upper boundary	shell O	0.42	0.77	103

Equations

$$\mathcal{L}=\Delta \nabla =\nabla -{\nabla }_{\mathrm{ad}}-{\nabla }_Y,$$

$$\mathcal{S}=\nabla -{\nabla }_{\mathrm{ad}},$$

$$u^2\sim g\ell \Delta \nabla >0,$$

$${\mathcal{M}}^2\sim (\ell /H_P)\Delta \nabla .$$

$${\partial }_t\mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}=\mathcal{B}-\mathcal{D}$$

$$u_{r}d{u}_r/dr=g{\beta }_T\Delta \nabla -u_r|u_r|/{\ell }_d.$$

$$\mathcal{S}=\Delta \nabla =\nabla -{\nabla }_{\mathrm{ad}}$$

$$\mathcal{L}=\Delta \nabla =\nabla -{\nabla }_{\mathrm{ad}}-{\nabla }_Y$$

$$u_{r}d{u}_r/dr=d(u_r^2/2)/dr\sim g{\beta }_T\ell \Delta \nabla .$$

$$d^2{u}_h/d{r}^2(u_h-u_0)<0,$$

$$X(z)=\theta +\frac{\varphi -\theta }{1+\exp [-\eta (z-1)]},$$

Full text

**DRAFT FROM **

1. Introduction

Observational and computational advances over the past decade have placed us in a unique era, in which study of the interiors of stars can be pursued at much higher precision than before. The recent observations from space of pulsating stars (through the MOST, CoRoT, Kepler, K2, and BRITE-constellation missions), and the planned future missions (like TESS and PLATO) supply precise photometry of stars, with near-continuous time sampling for durations of weeks to years. Among all targets, those more massive than $\sim 1.4\,M_{\odot}$ have a critical feature in common: during their main sequence lives, they harbor a convective core and a radiative envelope. The gravity (g) modes propagate in the radiative envelope, are reflected from the convective core, providing valuable information regarding the physical conditions near the boundary. Thus, g-mode pulsating stars are excellent tests of the physics of core convection.

In parallel, two- and three-dimensional (2D/3D) implicit large eddy simulations (ILES) of turbulent convection at different evolutionary phases (e.g., Bazàn & Arnett (1994); Nordlund & Stein (1995); Freytag, Ludwig, & Steffen (1996); Meakin & Arnett (2006, 2007b); Miesch (2005); Nordlund, Stein, & Asplund (2009); Viallet, et al. (2013); Arnett, et al. (2015), and many more) have shed light on the behavior of stellar convection, and in particular on the interface between convective and radiative zones. This allows the development of non-local time-dependent convective theories which are consistent with the 3D simulations. *The numerical data provides closure to the Reynolds-averaged Navier-Stokes (RANS) equations (Arnett, et al., 2015), converting numerical experiments into theory. * This approach is called 321D, and aims to provide alternatives (of increasing sophistication) to the classical Mixing Length Theory of convection (MLT, Böhm-Vitense, 1958), for use in one-dimensional (1D) stellar evolution codes.

Despite this progress in computation, simplifications are necessary. Attainable numerical resolution allows numerical Reynolds numbers $Re>10^{4}$, which are definitely turbulent, but stars have far higher Reynolds numbers (Arnett & Meakin, 2016). To attain highly turbulent simulations, only a fraction of the star is computed (a “box-in-star” approach) and rotation and magnetic fields are ignored. The simulations extend from the integral scale down into the Richardson-Kolmogorov cascade, and the subgrid dissipation merges with the Kolmogorov “four-fifths” law (Frisch, 1995; Arnett, Meakin & Young, 2009).

In addition to the Reynolds number issue, the simulations have negligible radiation diffusion (infinite Péclet number) because of vigorous neutrino cooling, rather than large but finite Péclet numbers found during hydrogen and helium burning (Viallet, et al., 2015).

These simulations represent turbulent solutions to the Navier-Stokes equations, and a step beyond mixing-length theory: they can resolve boundaries.

It is timely to compare the results of observations, asteroseismic modeling and 2D/3D simulations(e.g. Aerts & Rogers, 2015). Although these are independent approaches, it is possible to understand many underlying similarities between the state-of-the-art simulations and modeling. Such synergies shall allow us improve our treatment of turbulent convection in 1D models, and account for the convectively-induced mixing in the radiative interior (through overshooting and internal gravity waves) in a more consistent way.

We review the recent asteroseismic modeling of two Kepler slowly pulsating B stars (SPB’s) in §2. In §3 we compare the theoretical descriptions in MLT and 321D, region by region. In §4 we compare the shapes of the abundance profiles at the convective boundaries, as inferred from astreroseismology and from the 3D simulations of O-burning shell in a 23 $M_{\odot}$ model (Arnett, et al., 2015) and of C-burning shell in a 15  $M_{\odot}$ model (Cristini, et al., 2015, 2016). Our conclusions are summarized in §5.

2. Input Models

Moravveji, et al. (2015, 2016) recently did in-depth forward seismic modeling of two SPB stars having the richest seismic spectra known so far. Both KIC 10526294 (catalog ) (Papics et al., 2014), and KIC 7760680 (catalog ) (Papics et al., 2015) are of spectral class B8V (3.25 $M_{\odot}$), and exhibit long and uninterrupted series of dipole ($\ell=1$) g-modes. In both cases, the relative frequency differences between the observations and models is less than 0.5%. In addition to placing tight constraints on the extent of overshooting beyond the core, and additional diffusive mixing in the radiative interior, the authors concluded that the exponentially-decaying diffusive mixing profile for core overshoot outperforms the step-function prescription. Thus, convectively-induced mixing beyond the formal core boundary seems to have a radial dependence, and decays outwards.

Figure 1 shows the internal mixing profile (colored regions) of the best model of KIC 7760680, with selected positions in the enclosed mass coordinate labeled O, S, M, and W. We will discuss each region, comparing and contrasting the MLT, and a 321D theory based upon well-resolved numerical simulations (see Arnett, et al. (2015); Cristini, et al. (2016) and extensive references and discussion therein).

3. Comparing MLT and 321D

There are four distinct regions in Fig. 1, which we elaborate below: (1) a Schwarzschild core, (2) a braking (overshoot) region, (3) a composition gradient, and (4) a radiative envelope.

3.1. Region OS: Schwarzschild core

3.1.1 MLT

In Figure 1, the region extending from the origin at O to S is well mixed; at S the Schwarzschild criterion changes sign, so that within OS buoyancy drives convection.

For a composition gradient of zero, $-\nabla_{\mu}=\nabla_{Y}\equiv\Sigma_{i}X_{i}/A_{i}\rightarrow 0$, where $X_{i}$ and $A_{i}$ are the mass fraction and mass number of each nucleus $i$. The Ledoux discriminant,

[TABLE]

reduces to the the Schwarzschild discriminant,

[TABLE]

which is positive. $\nabla$ and $\nabla_{\rm ad}$ are actual and adiabatic temperature gradients, respectively, and $\nabla_{\mu}$ is the corresponding dimensionless gradient in mean molecular weight.

The convective velocity is approximately

[TABLE]

where $\ell$ is the free parameter, the mixing length, the temperature excess is $\Delta\nabla=\nabla-\nabla_{a}-\nabla_{Y},$ and $g$ is the gravitational acceleration. For the turbulent velocity, the MLT as given by Eq. 1 is an adequate approximation over region OS, and is the commonly used choice. This is a steady-state approximation which may be inferred from the turbulent kinetic energy equation (Meakin & Arnett, 2007b) or alternatively, from the balance of buoyant acceleration and deceleration by turbulent friction ($u|u|/\ell$); see Kolmogorov (1941), Kolmogorov (1962), Smith & Arnett (2014), and Arnett, et al. (2015).

Using the sound speed $s^{2}=\gamma PV$, the Mach number ${\cal M}=u/s$, and the pressure scale height $H_{P}=PV/g$, we may write Eq. 1 as

[TABLE]

For quasi-static stellar evolution, the convective Mach numbers are small (${\cal M}\leq 0.01$). Since $\ell/H_{P}\sim 1$, $\Delta\nabla\ll 1$ giving $\nabla\sim\nabla_{\rm ad}$ for a well-mixed convection zone in a stellar interior.

3.1.2 321D

Fully turbulent 3D simulations of convection may be represented as solutions to a simple differential equation for either the turbulent kinetic energy, or (as shown here) for the turbulent velocity $\bf u$ (Arnett, et al., 2015),

[TABLE]

where the buoyancy term111For strongly stratified media, there is an additional driving term due to “pressure dilatation” (Viallet, et al., 2013). is ${\cal B}\approx g\beta_{T}\Delta\nabla$ ($\beta_{T}=(\partial\ln\rho/\partial\ln T)_{P}$ is the thermodynamic factor222Sometimes denoted $\delta$ or $Q$. to convert temperature excess to density deficit at constant pressure), and the drag term (chosen to be consistent with the Kolmogorov cascade) is ${\cal D}\approx{\bf u}/\tau$, with the dissipation time $\tau=\ell_{d}/|u|$. Here $\ell_{d}$ is the dissipation length, a property of the turbulent cascade, and may differ from the mixing length used in MLT. We average over angles and take the steady state limit, which gives for the radial convective speed $u_{r}$,

[TABLE]

Away from the convective boundaries the gradient of $u$ is small, and this equation resembles Eq. 1 for appropriately chosen $\ell_{d}$. Eq. 3 implies a heating rate due to the dissipation of flow at small scales (Kolmogorov, 1941, 1962), $\epsilon_{K}=u^{2}/\tau=u^{3}/\ell$, a feature ignored in MLT. *This is the frictional cost to the star for moving enthalpy by turbulence. * In practice this term is small but not negligible. The 321D algorithm accounts for this additional heating term in the computation of 1D models (see §4.1.1 in Arnett, et al., 2015).

3.2. Region SM: Braking (overshoot)

3.2.1 MLT

The Schwarzschild discriminant,

[TABLE]

changes sign at boundary S, and $\nabla_{Y}\approx 0$ there, so Eq. 1 implies an imaginary convective speed. To deal with this singular behavior, different physics is traditionally introduced. A region SM is defined as the “overshoot” region, in which the luminosity is presumably carried entirely by radiative diffusion ($\nabla=\nabla_{\rm rad}$) and a new algorithm is defined, replacing the variable $u$ by a new variable, the effective diffusion rate $D_{ov}$; see Moravveji, et al. (2016) for a short and clear discussion. Inside SM, ${\nabla_{\rm rad}-\nabla_{\rm ad}\leq\cal S}\leq 0$. The coordinate M designates the layer at which this effective “convective diffusion” is no longer able to destroy the composition gradient, so over region SM we still have $\nabla_{Y}\approx 0$.

3.2.2 321D

Because the Ledoux discriminant,

[TABLE]

changes sign at S, ${\cal B}<0$, and the region SM is subjected to buoyant deceleration. Mixing continues over SM so that all of the region OSM is mixed, even though $\cal L$ is negative over SM. Thus SM is the overshoot region, in which the flow turns back to complete its overturn. The vector velocity $\bf u$ has different signs in upflow and downflow; it becomes horizontal at coordinate M, not zero; (see §3.8 in Arnett, et al., 2015). The coordinate M is a shear layer.

Here $\nabla_{Y}\rightarrow 0$, so we have $\nabla-\nabla_{\rm ad}<0$, and $g\ell\Delta\nabla<0$. Eq. 1 is not possible. Near the boundaries the velocity gradient terms dominate over the Kolmogorov term in Eq. 2, so

[TABLE]

For negative buoyancy ($\Delta\nabla<0$), the buoyant deceleration acts to decrease the radial component of the turbulent kinetic energy. Arnett, et al. (2015) show that this is essentially the same as defining the boundary by the gradient Richardson criterion (${\rm Ri}=N^{2}/(\partial u/\partial r)^{2}>1/4$).

The Schwarzschild criterion is only a linear instability condition, derived by assuming infinitesimal perturbations. The Richardson criterion is a nonlinear condition, which indicates whether the turbulent kinetic energy can overcome the potential energy implied by stable stratification. For a well-mixed region ${\cal L}\rightarrow{\cal S}$. Use of $\cal L$ in a stellar code can give fictitious boundaries due to small abundance gradients, which are blown away in a 3D fluid dynamic simulation; the Richardson criterion insures against these.

The boundary of a convective region has a negligible radial velocity of turbulence; Eq. 4 indicates where this occurs, so that by a solution of Eq. 3, 321D automatically determines where the boundary of convection is. This contrasts with MLT for which boundaries are undefined and thus a thorny issue.

3.2.3 Radiation diffusion effects

In order to compare the asteroseismic models from hydrogen burning to simulations from later burning stages, allowance must be made for the difference in strength of radiative diffusion in the two cases. The duration of neutrino-cooled stages is much shorter than the radiative leakage time from the core, while core hydrogen burning takes much longer than its corresponding radiative leakage time. In the terminology of fluid mechanics, the Péclet number is significantly different in these two cases (Viallet, et al., 2015), so their flows may be significantly different too. This is especially important for thin layers, for which the radiative diffusion time is shorter.

For oxygen burning and carbon burning, the flow follows a nearly adiabatic trajectory, having an entropy deficit in the braking region (Fig. 4 in Arnett, et al. (2015)). If radiative diffusion is not negligible, heat will flow into this region of entropy deficit, reducing the strength of the buoyancy braking and thus widening the overshoot layer. This causes the actual gradient to deviate from the adiabatic gradient ($\nabla_{\rm ad}$) and approach the radiative one ($\nabla_{\rm rad}$). The temperature gradient in the overshooting region is expected to lie between the adiabatic and radiative ones $\nabla_{\rm rad}<\nabla<\nabla_{\rm ad}$ (Zahn, 1991; Zhang, 2016).

Calibration of overshooting algorithms using observations of hydrogen-burning stars may be inadequate for later burning stages because of the differences in Péclet number (Viallet, et al., 2015). This could be a problem for helium burning (Mosser et al., 2014; Schindler, et al., 2015; Constantino, et al., 2015; Ghasemi, et al., 2016), and will certainly get worse for the neutrino-cooled stages.

3.3. Region MW: Composition Gradient

Moravveji, et al. (2015, 2016) found it desirable to add an extra diffusive mixing (their $D_{\rm ext}$) beyond the well-mixed region OSM, because the agreement between observed and modeled frequencies improve – in $\chi^{2}$ sense – by a factor 11. This is an important clue regarding the structure of this region.

The physical basis for extra mixing seems to be at least two-fold. Because coordinate M is a shear layer, Kelvin-Helmholtz instabilities will mix matter above and below coordinate M. In Eq. 3 we saw that the radial velocity could be decelerated at a boundary. As the flow turns, a horizontal velocity $u_{h}$ must develop (Arnett, et al., 2015). This variable does not appear in MLT (the “blob disappears back into the environment”). A finite $u_{h}$ implies a shear instability may occur. A necessary condition for instability (due to Rayleigh and to Fjørtoft, see §8.2 in Drazin (2002)) is that

[TABLE]

somewhere as we move in radius through the flow at the boundary. Here $u_{h}$ is the horizontal velocity and $u_{0}$ is that velocity at the radius at which $d^{2}u_{h}/dr^{2}=0$ (the point of inflection). While a stably-stratified composition gradient may tend to inhibit the instability, Eq. 5 illustrates a basic feature: the velocity field drives the instability.

If this horizontal flow is stable, a composition gradient can be maintained. Turbulent fluctuations may cause the stability criterion to be violated, and thus erode the layer until it is again stable. If there is entrainment at the boundary, as the 3D simulations show (Meakin & Arnett, 2007b), then the boundary moves away from the center of the convection zone. The result is that the composition gradient will be left near the margin of instability. This gives the sigmoid shape, found in 3D simulations of oxygen burning (Arnett, et al., 2015) and carbon burning (Cristini, et al., 2016).

At mass coordinate M the average turbulent velocity in the radial direction is zero. Because the flow is turbulent, it is zero only on average, so $\langle u\rangle=0$, but has a finite rms value, $\langle u^{2}\rangle\neq 0$. Over the regions SM and MW, $\Delta\nabla<0$, so the Brunt-Väisälä frequency is real. Thus, this region can support waves, which the velocity fluctuations necessarily excite. The turbulent fluctuations couple best to g-modes at low Mach numbers (Meakin & Arnett, 2007a). The mixing induced here is much slower than in the convection region333Most of the turbulent energy is at large wavelengths (g-modes), so the induced waves are less able to cause mixing than shorter wavelengths might., and allows a composition gradient to be sustained over MW. Because the mixing is slow, the composition structure may be a combination of previous history and slow modification.

3.4. Region W: Radiative envelope

To the extent that nonlinear interactions of waves generate entropy, the circulation theorem implies that slow currents will be induced, leading to more “extra mixing” in radiative regions such as those above coordinate W. In 1D stellar evolution, this type of mixing is treated diffusively ($D_{\rm ext}$), and asteroseismology can tightly constraint it.

4. Comparing Abundance Profiles

Figure 2 shows the inferred profile of hydrogen in the two Kepler SPB’s. This may be well-fitted by a logistic function strikingly similar to the composition profile shapes resulting from 3D simulations. We choose to fit the abundance profiles $X$ (in mass fraction) to

[TABLE]

where $z=r/r_{\rm mid}$ is the normalized radius and $r_{mid}$ is the radius at the mid-point of the abundance gradient. The steepness of the composition profile is encoded in $\eta$. At the “‘core”, $-\eta(z-1)\gg 1$, and $X(z)\rightarrow\theta$. At $r=r_{\rm mid},$ $X_{\rm mid}=(X_{\rm core}+X_{\rm surf})/2$, so $z=1$ and $X_{\rm mid}=0.5(\theta+\phi)$. At the stellar “surface” ($z\gg 1$) and for large $\eta$, $X_{\rm surf}=\phi$. Convective boundaries in stars are often narrow with respect to radius. The mathematical form in Eq. 6 stretches the abundance gradient over the whole $z$ axis. Thus, it is straightforward to fit the sigmoid shape to any abundance profile from stellar models at any evolutionary phase from 1D or 2D/3D (after temporal and angular averaging). The normalization used in Eq. 6 tends to separate shape from absolute scale, so that it may be used for different compositions (burning stages).

Radiative transfer is negligible for carbon burning and oxygen burning, so the Peclét number is essentially infinite (Viallet, et al., 2015), while radiative transfer does affect hydrogen burning. In the braking region, radial deceleration (turning of the flow) is due to negative buoyancy. At the top of a convection zone, rising matter becomes cooler (has lower entropy) than its surroundings. Radiative transfer tends to heat this cooler matter, reducing the negative buoyancy (radial braking) so that the matter must continue further before it is turned. This results in a wider braking layer (smaller $\eta$); see discussion in (Arnett, et al., 2015). Boundary layers are narrower for neutrino-cooled stages, so that calibration of boundary widths from photon-cooled stages are systematically in error. Simulations in 3D of H burning are difficult without artificial scaling of the heating rate (Meakin & Arnett, 2007b). Because of negligible radiative diffusion in carbon and oxygen burning, those corresponding values for $\eta$ in Table 1 should tend to be higher than would be expected in hydrogen burning, as they do.

Simulations of oxygen burning give behavior similar to carbon burning, but are more complicated to interpret because of a significant initial readjustment of convective shell size, and an episode of ingestion of 20Ne; see Arnett, et al. (2015) and references therein to earlier works. Despite this, the large value of $\eta$ reasonably represents the time- and angle-averaged O-profiles.

In Table 1 we compare upper boundaries because core H burning has only an upper boundary. The large values of $\eta\gg 1$ imply that the composition boundaries are narrow relative to the radius, in all cases. The $\eta$ for KIC 10526294 is roughly twice that of KIC 7760680. The former star is very young ($X_{\rm core}=0.63$) with a steeper composition jump; the latter is more evolved ($X_{\rm core}=0.50$) with broader $\mu$-gradient zone. Further investigation of precisely which physical processes determine the parameter $\eta$ is underway (see §3.3).

In Fig. 2, the maximum and minimum mass fractions of fuel are simply a function of the chosen consumption by burning, fixing $\phi$ and $\theta$. The position of the boundary is normalized relative to the center of the gradient $r_{mid}$ (which avoids the important issue of entrainment). What is important here is the narrowness of the composition gradient (the large value of $\eta$) and the shape of the curve joining the high and low fuel abundances, both of which are predicted by the 3D simulations, and independently inferred from asteroseismology, giving an encouraging agreement.

5. Summary

Asteroseismically-inferred composition gradients are strikingly similar to those found in 3D simulations; MLT cannot predict boundary structure. We compare, zone by zone, the predictions of MLT and 321D, from the convective core to the radiative mantle. Advantages of 321D are its time-dependence, non-locality, incorporation of the Kolmogorov cascade and turbulent heating, and resolution of convective boundaries. The properties of the fully convective regions are similar in 321D and MLT. Use of the 321D can avoid imaginary convective velocities in the braking regions, which are related to the development of singularities in boundary layers (Landau & Lifshitz (1959) §40, Gough (1977)). The 321D can provide a continuous description of the convective boundary – from fully convective to radiative – avoiding the awkward patching characteristic of MLT.

The 321D procedure promises a dynamical treatment of the overshooting and wave generation in stellar models. Possible effects on convective flow from radiative diffusion should be further explored for regions with moderate Péclet numbers (Viallet, et al., 2013), between existing deep interior and atmospheric simulations. Better boundaries can provide a possible solution to the behavior of convective helium burning cores, which the MLT fails to represent well.

The sigmoid fits to the hydrogen profiles of the two Kepler targets come from 1D (MESA) including ad-hoc core overshooting and extra diffusive mixing to match the observed g-mode frequencies. In contrast the simulations of C and O burning shells have no free parameters to adjust. Despite dealing with such different burning stages, the similarity shown in Table 1 and Fig. 2 is striking.

Special thanks are due to Simon Campbell and to Andrea Cristini for access to their simulation data prior to publication, and to Raphael Hirschi, Casey Meakin, Cyril Georgy, Maxime Viallet, John Lattanzio and Miro Mocák for helpful discussions. We thank an anonymous referee who helped to improve the manuscript. This work was supported in part by the Theoretical Program in Steward Observatory, and by the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme FP7/2007-2013/ under REA grant agreement N∘ 623303 (ASAMBA).