On the Chemical Mixing Induced by Internal Gravity Waves (IGW)

Detailed modeling of stellar evolution requires a better understanding of the (magneto-)hydrodynamic processes which mix chemical elements and transport angular momentum. Understanding these pro- cesses is crucial if we are to accurately interpret observations of chemical abundance anomalies, surface rotation measurements and asteroseismic data. Here, we use two-dimensional hydrodynamic simula- tions of the generation and propagation of internal gravity waves (IGW) in an intermediate mass star to measure the chemical mixing induced by these waves. We show that such mixing can generally be treated as a diffusive process. We then show that the local diffusion coefficient does not depend on the local fluid velocity, but rather on the wave amplitude. We then use these findings to provide a simple parametrization for this diffusion which can be incorporated into stellar evolution codes and tested against observations.


INTRODUCTION
Accounting for hydrodynamic processes in stellar interiors over stellar evolution times has been the biggest source of uncertainty when comparing theoretical results with observations. While Mixing Length Theory (MLT) has proven extremely useful for characterizing mixing within convection zones (Bohm-Vitense 1958; Kippenhahn et al. 2012), there remain many uncertainties dealing with this mixing at convective-radiative interfaces (Renzini 1987;Zahn 1991) and within radiative regions (Pinsonneault 1997;Heger et al. 2000). Nearly all stars host radiative regions so it is critical we develop methods for accurately parametrizing chemical mixing (and angular momentum transport) in these regions.
Numerous theoretical models have been proposed for incorporating mixing by (magneto-)hydrodynamical processes in stellar radiative zones into stellar evolution models. A myriad of hydrodynamic instabilities (Heger et al. 2000), rotationally induced mixing (Eddington-Sweet circulation, Eddington (1925); Vogt (1925); Sweet (1950)) and magnetically induced mixing (Taylor- Spruit dynamo Spruit (2002)) have been included in modern stellar evolution codes. Still, many questions remain about mixing within stellar radiative interiors. For example, it is typically assumed that rotationally induced mixing is dominant in massive stars, yet the observed lack of correlation between Nitrogen abundance and rotation rate for some stars indicates additional mixing is needed (Brott et al. 2011) and observations in multivariate parameter space suggest pulsational mixing is dominant for slow to moderately rotating OB stars (Aerts et al. 2014). Similarly, differential rotation at late stages of stellar evolution is lower than expected even when all of these mechanisms are considered (Eggenberger et al. 2017).
In general, these multi-dimensional hydrodynamic effects are parametrized as a diffusion coefficient within stellar evolution codes. Each physical process has an instability criteria based on local properties (e.g. shear). Once instability is confirmed a diffusion coefficient is constructed from the lengthscale and growth rate of the instability. This diffusion coefficient is then included locally (in space and time) in the stellar evolution calculation. While this procedure is rather rudimentary it is clear from observations that additional mixing within stellar radiative regions is required. It is only recently, and in limited circumstances, that such prescriptions are being tested against hydrodynamic calculations which self consistently calculate the development of the instability and the subsequent mixing induced (Edelmann et al. 2017).
Internal gravity waves (IGW) are known to propagate and dissipate in radiative regions which could lead to chemical mixing and angular momentum transport. However, the parametrization of IGW in one-dimensional (1D) stellar evolution codes is complex. The transport of angular momentum by these waves can not be treated as a diffusive process, indeed IGW have an anti-diffusive nature. That is, they drive, rather than dissipate, shear flows (Buhler 2009). For this reason, IGW transport has generally not been treated in 1D stellar evolution codes (except in the Geneva code , in which their treatment is complex). Yet, while it is clear that angular momentum transport by IGW can not be parametrized with a diffusion coefficient (Rogers 2015), it is unclear whether the chemical mixing induced by waves could be treated diffusively as previously suggested (Press & Rybicki 1981;Garcia-Lopez & Spruit 1991). The purpose of this letter is to first determine whether wave mixing can be treated diffusively and, if so, determine how efficient that mixing is and whether it could be reasonably parametrized in 1D stellar evolution models.

Hydrodynamic Simulations
In order to measure mixing by IGWs in stellar interiors, we solve the Navier-Stokes equations in the anelastic approximation (Gough 1969;Rogers & Glatzmaier 2005). The equations are solved in two dimensions (2D) representing an equatorial slice of the star. Our reference state model is that of a 3M ⊙ star with a core hydrogen content of 0.35, calculated using Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. (2010Paxton et al. ( , 2015). We solve the equations from 0.03R * to 0.70R * , the initial rotation rate is uniform and equal to 10 −6 rad/s. The simulation is run a total of 4×10 7 s, or approximately 40 convective turnover times. The details of the equations and numerical methods used can be found in Rogers et al. (2013). Time snapshots of vorticity are shown in Fig. 1. Ideally, we would solve the equations in 3D as in Alvan et al. (2015), but few simulations such as these have been done for massive stars (Browning et al. 2004;Augustson et al. 2016) and none which include extended radiative zones. We discuss the role this reduced dimensionality might play in the Discussion.
Like all hydrodynamic simulations our viscous and thermal diffusion coefficients are larger than is typical in real stars. Therefore, our waves are damped more than they would be in actual stellar interiors. Unlike the simulations presented in Rogers et al. (2012Rogers et al. ( , 2013, these models are not "over-forced" so their root-meansquared velocities in the convection zone are similar to that expected from mixing length theory (see Fig. 3).

Tracer Particles and Diffusion
To determine whether mixing by IGW behaves like diffusion we introduce tracer particles within the simulation. We use a tracer particles instead of solving a compositional advection-diffusion equation because such a method would require an explicit diffusion that would dominate any mixing coefficient. For simplicity we regard time as continuous and only measure radial diffusion. At some time of interest we introduce N particles, distributed uniformly in space, and track them for a time T to produce N particle trajectories (particle positions in time, R i (t)). We then consider all the sub-trajectories of duration τ and use a cubic spline function w to interpolate between the start of each sub-trajectory (R i (t)) and the grid position. The length of a sub-trajectory (displacement of a particle) is then R i (t + τ ) − R i (t). We then calculate the following profiles: Here n(r, τ ) is the number of sub-trajectories starting at r of duration τ . P (r, τ ) is the sum of lengths of these subtrajectories and Q (r, τ ) is the sum of the lengths squared of these sub-trajectories. If there is a mean velocity field u(r) then Lower values of τ can result from many time differences, while larger values of τ can only result for long timeline data. 1 Therefore, low values of τ represent more data for a given T , while larger values represent fewer data points. If particle motion is purely diffusive with zero mean then, provided that the distance moved by particles in time τ is smaller than the distance over which the diffusion coefficient varies, we would have where D(r, τ ) is the diffusion coefficient at r. If there is a mean velocity field in addition to diffusion then we would have If the motion is purely diffusive then D will not depend on τ , but if there is a more complicated background velocity field the situation is more complicated as we see in the next section.

Wave Motion with Diffusion
In order to under understand the motion of a particle in a wave acted on by diffusion we consider a particle moving in a wave field with velocity given by ωA cos(φ + ωt) and with random fluctuations given by N (t) so that its equation of motion is This can be integated to give Now we assume that N (t) is Gaussian white noise with standard deviation σ so that t 0 N (s)ds = σW (t), where W (t) is a Wiener process (Doob 1953). We consider an ensemble of trajectories and average uniformly over the phases φ and the Wiener process using W (t) = 0 and W (t)W (s) = min(s, t). Then we have, the expectation (denoted with · · · ) of the position and This function of τ is the sum of a linear function (A 2 + σ 2 τ ) and an an oscillatory function −A 2 cos(ωτ )). At large times (τ ω ≫ 1) there is a linear trend with gradient σ 2 . From this gradient we can then extract the effective diffusion coefficient D = σ 2 /2.

Application to Numerical Simulations
We carry out this procedure in post-processing. That is, given saved velocity data from our hydrodynamic simulations at time intervals t, we introduce N particles into our numerical simulation and measure D using the above procedure. Since the procedure is done in post-processing, we can include an arbitrary number of particles to check numerical convergence, we can also vary the type of particle interpolation and the number of timesteps over which we track the particles. Fig. 1 shows time snapshots of the vorticity and particle positions (only 15000 shown) at four different times.
We checked our procedure against a hydrodynamic simulation with particles run within the simulation in order to confirm our velocity data was taken at fine enough time resolution. We find that in order to measure a diffusion coefficient within the radiation zone, we need long timeline data, but not very finely spaced. However, to smoothly resolve the convective-radiative interface we need finely spaced time data over shorter timeline. In this letter we are concerned with the diffusion profile within the radiation zone, so we integrate over long timelines with longer time steps and note that our convective and overshoot profile is not well resolved nor well described as a diffusive process.

PARTICLE TRAJECTORIES AND DIFFUSION
The first question to address is whether the movement of particles due to the action of waves can be treated as a diffusive process. To test this we plot the mean squared There we see that at low τ the particles indeed behave like diffusion (the mean squared displacement is a linear function of τ ). However at larger τ , there is not enough data to confirm the diffusive nature. We also see that different radial levels have different slopes.
Using low values of τ (2.5 ×10 6 s, for which we have sufficient data), we compute the diffusion coefficient as a function of radius, which is shown in Fig. 2(b), again with different line types showing different numbers of particles. The vertical black solid line shows the convectiveradiative interface, while the vertical black, blue and red dashed lines show the radial positions shown in (a). There we see that the overall radial profile of diffusion coefficient within the radiation zone is robust to variations in particle number, thus demonstrating convergence. We also see that there is a rapid transition from the behavior in the convection zone to that in the radiative region. We note that, while there is a diffusion coefficient (plotted and measured) within the convection zone and overshoot region, the behavior in those regions is generally not diffusive, particularly in the convection zone. The mean squared displacement as a function of τ does not lie on a straight line. Therefore, applying a diffusion coefficient in these regions is not appropriate. However, we can see from Fig. 1 that particles mix much more rapidly within the convection zone than in the radiative region, as ex- Figure 3. Root mean squared velocity (black asterisks) and diffusion coefficients versus time. Diffusion coefficients are at 0.375 R * (red diamonds) and at 0.675 R * (blue asterisks) and are calculated using Equation (6) and a particle tracking time, T=10 6 s.

pected.
In general, within the radiation zone, the amplitude of the diffusion coefficient rises with increasing distance from the convection zone. There is decay just outside the convection zone due to the fact that our thermal diffusivity is constant, rather than a function of radius (we discuss this in Section 4). The radial increase is proportional to a factor between ρ −1 and ρ −1/2 which we show in black dashed lines and which we also discuss in Section 4. Therefore, we conclude that: 1) IGW mixing in the radiation zone can be treated as a diffusive process and 2) the radial profile of diffusion is robust and is proportional to a function between ρ −1 and ρ −1/2 .

PARAMETRIZATION OF THE DIFFUSION PROFILE
In order to find a useful parametrization of IGW to be used in one-dimensional (1D) stellar evolution models, we would like to know both the amplitude and the radial profile of diffusion from our simulations. We find that the amplitude of the diffusion coefficient we obtain within the radiation zone is correlated with the root mean squared (rms) velocity within the convective zone. In the model presented, the root mean squared velocity in the convection zone is decreasing in time. This is an artifact of the fact that we force the convection through a superadiabaticity, which is reduced due to efficient convection. In previous simulations we have a used a forcing term to drive convection in order to avoid this, but here we allowed this to investigate the role of convective velocities. We see how this decay affects the diffusion coefficient in the radiation zone in Fig. 3 -as convective velocities decrease, the diffusion coefficient decreases. As we show later, this is because the rms velocity within the convection zone sets the wave amplitudes within the radiation zone. Within a 1D stellar evolution code the rms velocity could be approximated as the convective velocity given by mixing length theory (MLT).
Independent of the amplitude of the diffusion coefficient, we find that the radial profile of diffusion seen in Fig. 2 is robust across choices of numerical diffusion (the thermal and viscous diffusivity are each varied by an order of magnitude) and convective velocities. Now we Figure 4. Root mean squared velocity (black asterisks, cm 2 s −2 ), wave amplitude squared (red crosses,cm 2 s −2 ) and diffusion coefficient (blue crosses,cm 2 s −1 ) versus radius. Wave amplitude squared is calculated using the larges scale wave (k h =1) and all frequencies. Black lines represent the wave amplitude squared calculated using Eqn. (12). The solid line assumes κ = 2 × 10 12 cm 2 /s as in the numerical simulation and a frequency spectrum at generation of velocities ∝ ω −1 , while the dotted and dashed lines use κ from MESA and a frequency spectrum at generation of velocities ∝ ω −1 , ω −3 , respectively. would like to determine what sets the radial profile of diffusion. In Fig 4 we plot the rms velocity (black asterisks), the wave amplitude squared (red crosses) and the diffusion coefficient (blue crosses) as a function of radius. The wave amplitude is calculated with Eqn. (12) using the largest scale wave and all frequencies we calculate (up to 500 µHz). In general, frequencies above ∼40 µHz do not contribute as their generation amplitudes are too low. We see that the diffusion profile closely correlates with the wave amplitude squared. This numerical finding is consistent with the theoretical prediction that diffusion is the autocorrelation of the lagrangian velocity field. In this case the velocity field is the wave amplitude, which is correlated over long timescales. The diffusion profile is not correlated with the local rms velocity, which is again consistent with with the theoretical expectation that diffusion is the autocorrelation of the velocity field, as the local rms velocities are not correlated over long timescales and hence have negligible autocorrelation.
Due to numerical constraints (dimensionality, diffusion coefficients) it is likely the IGW amplitudes in our simulations are not realistic. Using our numerical result that diffusion is proportional to the wave amplitude squared, we now turn to theoretical models for a parametrization of wave diffusion in stellar interiors. The amplitude of an internal gravity wave (IGW) depends on the wave driving at the convective-radiative interface, the radiative damping of the wave and the density stratification. The wave driving is directly correlated with the rms velocity within the convection zone, v rms−cz (Garcia-Lopez & Spruit 1991). Within the radiation zone the wave is damped by thermal diffusion (Kumar et al. 1999) and amplified by the density stratification (due to conservation of pseudomomentum (Buhler 2009)). Therefore, the wave amplitude is determined from the simple function where ρ tcz is the density at the top of the convection zone. From Fig. 4 we estimate the diffusion due to IGW as where A has units of s and τ (ω, k h , r) is the damping "optical depth" of a wave defined in Kumar et al. (1999) as: where κ is the radiative diffusivity, k h is the horizontal wavenumber of the wave, r tcz is the radius at the top of the convection zone, N is the Brunt-Väisälä frequency and ω is the frequency of the wave. For simplicity, we have assumed no doppler shift in the frequencies of the waves. A is an unknown constant, which is ∼1s in our models. Although the precise value is unknown, we don't expect it to vary significantly from this (see Discussion). The initial decay of diffusion outside the convection zone is due to the fact that we use a constant thermal diffusion coefficient, κ, rather than the stellar radiative value. The numerical value used is 2 × 10 12 cm 2 /s throughout, while the value in the star varies from 10 7 cm 2 /s at the convective-radiative interface to 10 15 cm 2 /s at the surface of the star. Therefore, throughout our computed radiative zone we are damping the waves more than they would be damped in the stellar interior, this is particularly true just outside the convection zone and is the reason for the initial decay seen in the wave amplitude squared (red line, Fig. 4). This is demonstrated by the black lines which show the wave amplitude computed using Eqn. (12) with the value of κ used in our simulation (solid line) and the values of κ from the stellar model (dashed and dotted lines). For this simple calculation we have assumed that the waves are linear, non-interacting and that the wave amplitude at the convective-radiative interface is half the rms velocity. 2 Predictions for the frequency spectra of waves generated by convection at a given wavenumber range from ω −3 (Kumar et al. 1999;Lecoanet & Quataert 2013) to ω −1 (Rogers & MacGregor 2010;Rogers et al. 2013). Therefore, the dashed and dotted lines represent a frequency generation spectra of velocity proportional to ω −3 and ω −1 respectively, to cover that range. We integrate over the same frequency range and scales as for the numerical results (red asterisks in Fig. 4). We see that the initial decay outside the convection zone is purely an artifact of enhanced diffusion for numerical purposes and is not physical.

DISCUSSION
2 This amplitude comes from a simple calculation assuming Fw ∼ M Fc (Lecoanet & Quataert 2013), where M is the Mach number of the convection and Fc and Fw are the convective and wave fluxes, respectively. We assume that Fc, Fw ∼ v 3 c , v 3 w , therefore, vw ∼ M 1/3 vc. Assuming M ∼ 0.1 then vw ∼ 0.47vc.
In this letter we have demonstrated that the mixing by IGW within radiative regions can be treated as a diffusive process. We have further shown that the local amplitude of the diffusion coefficient depends on the local wave amplitude. The wave amplitude, in turn, depends on the convective forcing (v rms−cz ), the thermal damping of the wave and the density stratification in a simple way. Therefore, a prescription for the diffusion coefficient due to mixing by IGW can be easily implemented using Eqns. (12-14) assuming MLT velocities for the rms velocity and all other parameters determined by the stellar model. The one parameter that is left is A. While this value is ∼1 in our simulations, its precise value may depend on the stellar viscosity/thermal diffusivity, rotation and the dimensionality. Thermal diffusivity is already accounted for in (13). Since viscosity is enhanced in our simulations one would expect A = 1 to be a lower limit. However, one does not expect viscosity to play a role in the propagation of linear waves, so the prescription in (13) with A = 1 likely still holds. In the case of rotation, fast rotation would likely reduce the wave amplitude and hence, A = 1 would be considered an upper limit, but that is likely a small effect.
It is unclear what effect our reduced dimensionality has on the waves. Assuming that simple one dimensional wave propagation is sufficient, the dimensionality likely only affects the wave generation spectrum. Taking all this into account, the best approach would be to assume A = 1 and vary the wave generation spectrum incorporated through Eqn. (14). Then the one parameter of the model would be the exponent of the frequency spectrum of wave generation.
Given numerical limitations, for the forseeable future the most reliable constraints on A would come from comparisons between theoretical models using this prescription and observations of slowly rotating intermediate mass stars. Asteroseismic inversions could place constraints on near-core mixing (Moravveji et al. 2015(Moravveji et al. , 2016, while spectroscopic observations may place constraints on subsurface mixing. Simultaneous comparisons between theoretical evolution models, spectroscopic and asteroseismic data could provide constraints on the entire diffusion profile and indeed may help place constraints on the wave generation spectrum. Finally, our simulations only extend to 0.7R * . Extending our linear calculations using Eqn. (12) shows that the wave amplitude continues to increase until just beneath the stellar surface. It is likely that these waves break (Rogers et al. 2013) and hence the surface diffusion coefficient may be enhanced (due to turbulent mix-ing) beyond what is expected from linear wave behavior. Numerical simulations attempting to resolve the wave dynamics near the stellar surface will be forthcoming.