Synergies between Asteroseismology and Three-dimensional Simulations of Stellar Turbulence
W. David Arnett E. Moravveji






1Steward Observatory, University of Arizona, Tucson, AZ 85721, darnett@as.arizona.edu

2Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium


DRAFT FROM December 30, 2023
 

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∼ 3.25 M_⊙) to predictions of 321D turbulence theory, 
based upon well-resolved, truly turbulent three-dimensional simulations <cit.> 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.









§ 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 ∼1.4 M_⊙ 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., <cit.>, 

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 <cit.>, 
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 <cit.>, 
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 <cit.>. 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 <cit.>.

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 <cit.>.


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
<cit.>. 

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 <ref>.



In <ref> we compare the theoretical descriptions in MLT and 321D, region by region.
In <ref> 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_⊙ model <cit.> and of C-burning shell in a 15  M_⊙ model <cit.>. Our conclusions are summarized in <ref>.

 


§ INPUT MODELS
 
 <cit.> recently did 
 in-depth forward seismic modeling of two SPB stars 
having the richest seismic spectra known so far.
Both KIC 10526294 <cit.>, and KIC 7760680 <cit.> are of 
spectral class B8V (3.25 M_⊙), and exhibit long and uninterrupted series of dipole (ℓ=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 <ref> 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 <cit.> and extensive references and discussion therein).







§ COMPARING MLT AND 321D


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



 §.§ Region OS:  Schwarzschild core




  §.§.§ MLT

In Figure <ref>, 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,  -∇_μ = ∇_Y ≡Σ_i X_i/A_i → 0, where X_i and A_i are the mass fraction and mass number of each nucleus i.
The Ledoux discriminant,
 
    L = Δ∇ = ∇ - ∇_ ad - ∇_Y,
 
reduces to the the Schwarzschild discriminant,
 
    S  = ∇ - ∇_ ad,
 
which is positive. 
∇ and ∇_ ad are actual and adiabatic temperature gradients, respectively,  and ∇_μ is the corresponding dimensionless gradient in mean molecular weight.


 
 The convective velocity is approximately

    u^2 ∼ g ℓΔ∇ >0,

where ℓ is the free parameter, the mixing length,  
the temperature excess is Δ∇ = ∇ - ∇_a -∇_Y, and g is the gravitational acceleration. 
For the turbulent velocity, the MLT as given by Eq. <ref> 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 
<cit.> or alternatively, from the balance of buoyant acceleration and deceleration by turbulent 
friction (u|u|/ℓ); see <cit.>, <cit.>, <cit.>, and <cit.>.


Using the sound speed s^2 =γ PV, the Mach number M = u/s, and the pressure scale height H_P = PV/g,  we may write Eq. <ref> as

    M^2 ∼ (ℓ / H_P) Δ∇ .
 
For quasi-static stellar evolution, the convective Mach numbers are small (M≤ 0.01). Since ℓ/H_P ∼ 1,  Δ∇≪ 1 giving ∇∼∇_ ad for a well-mixed convection zone in a stellar interior.





  §.§.§ 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 u <cit.>,

    ∂_t  u + ( u ·∇) u =  B  -  D

where the buoyancy  term[For strongly stratified media, there is an additional driving term due to “pressure dilatation” <cit.>.] is B≈  gβ_T Δ∇  (β_T = (∂lnρ/ ∂ln  T )_P is the thermodynamic factor[Sometimes denoted δ 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 D≈ u/τ, with the dissipation time τ = ℓ_d / |u|. 
Here ℓ_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,

    u_r d u_r/dr =  gβ_T Δ∇  - u_r|u_r|/ℓ_d.

Away from the convective boundaries the gradient of u is small, and this equation resembles Eq. <ref> 
for appropriately chosen ℓ_d.

Eq. <ref> implies a heating rate due to the dissipation of flow at small scales <cit.>, ϵ_K = u^2/τ = u^3/ℓ, 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 <cit.>.
 


 §.§ Region SM: Braking (overshoot)
 


  §.§.§ MLT

The Schwarzschild discriminant, 
    S = Δ∇ = ∇ - ∇_ ad
 changes sign at boundary  S, 
and ∇_Y ≈ 0 there, so Eq. <ref> 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 (∇  = ∇_ rad) and a new algorithm is defined, replacing the  variable u by a new variable, the effective diffusion rate D_ov; see <cit.> for a short and clear discussion. 
Inside SM, ∇_ rad-∇_ ad≤ S≤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 ∇_Y ≈ 0.




  §.§.§ 321D

Because the Ledoux discriminant, 
    L = Δ∇ = ∇ - ∇_ ad - ∇_Y
 changes sign at S, 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 L is negative over SM. Thus SM is the overshoot region, in which the flow turns back to complete its overturn. The vector velocity u has different signs in upflow and downflow; it becomes horizontal at coordinate M, not zero; <cit.>. The coordinate M is a shear layer.


Here ∇_Y → 0, so we have  ∇ - ∇_ ad < 0, and g ℓΔ∇ <0.  
Eq. <ref> is not possible.
Near the boundaries the velocity gradient terms dominate over the Kolmogorov term in Eq. <ref>, so

    u_r du_r/dr = d (u_r^2/2)/dr ∼ g β_T ℓΔ∇.

For negative buoyancy (Δ∇ <0), the buoyant deceleration acts to decrease the radial 
component of the turbulent kinetic energy.

<cit.> show that this is essentially the same as defining the boundary by the gradient Richardson 
criterion (Ri = N^2/(∂ u /∂ 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 L→ S.
Use of 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. <ref> 
indicates where this occurs, so that by a solution of Eq. <ref>, 321D automatically determines where the 
boundary of  convection is.  
This contrasts with MLT for which boundaries are undefined and thus a thorny issue. 










  §.§.§ 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 <cit.>, 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 <cit.>). 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 (∇_ ad) and approach the radiative one (∇_ rad).

The temperature gradient in the overshooting region is expected to lie between the adiabatic and radiative ones ∇_ rad<∇<∇_ ad <cit.>. 



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 <cit.>. This 
could be a problem for helium burning <cit.>, and will certainly get worse for the neutrino-cooled stages.




 §.§ Region MW: Composition Gradient

<cit.> found it desirable to add an extra diffusive mixing (their D_ ext) beyond the well-mixed region OSM, 
because the agreement between observed and modeled frequencies improve – in χ^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. <ref> we saw that the radial velocity could be decelerated at a boundary. As the flow turns, a horizontal velocity u_h  must develop <cit.>. 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 <cit.>) is that 

    d^2 u_h/dr^2 (u_h-u_0) < 0,
 
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^2u_h/dr^2 = 0 (the point of inflection).
While a stably-stratified composition gradient may tend to inhibit the instability, Eq. <ref> 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 <cit.>, 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 <cit.> and carbon burning <cit.>.


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 ⟨ u ⟩ = 0, but has a finite rms value, ⟨ u^2 ⟩≠ 0.
Over the regions SM and MW, Δ∇ < 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 <cit.>.
The mixing induced here is much slower than in the convection region[Most 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. 




 §.§ 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_ ext), and 
asteroseismology can tightly constraint it.



§ COMPARING ABUNDANCE PROFILES


Figure <ref> 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

    X(z) = θ + ϕ - θ/1 + exp[-η (z-1)],

where z=r/r_ 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 η.
At the “`core”, -η(z-1) ≫ 1, and X(z) →θ.

At r=r_ mid ,  X_ mid= (X_ core + X_ surf)/2, 

so z=1 and X_ mid=0.5(θ+ϕ).
At the stellar “surface” (z ≫ 1) and for large η, X_ surf=ϕ.
Convective boundaries in stars are often narrow with respect to radius.
The mathematical form in Eq. <ref> 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. <ref> 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 <cit.>, 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 η); see discussion in <cit.>. 
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 <cit.>. 
Because of negligible  
radiative diffusion in carbon and oxygen burning,  those 
corresponding values for η in Table <ref> 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 <cit.> and references therein to earlier works. Despite this, the  large value of η  reasonably represents the time- and angle-averaged O-profiles.





In Table <ref> we compare upper boundaries because
core H burning has only an upper boundary. 


The large values of η≫ 1 imply that the composition boundaries are narrow 
relative to the radius, in all cases. 

The η for KIC 10526294 is roughly twice that of KIC 7760680.

The former star is very young (X_ core=0.63) with a
steeper composition jump; the latter is more evolved (X_ core=0.50) with broader
μ-gradient zone.
Further investigation of precisely which physical processes determine the parameter η is underway (see <ref>).

In Fig. <ref>, the maximum and minimum mass fractions of fuel are simply a function of the chosen consumption by burning, fixing ϕ and θ. 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 η) 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. 




§ 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 (<cit.> 40, <cit.>).
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 <cit.>, 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 <ref> 
and Fig. <ref> 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).

















[Aerts & Rogers (2015)]conny Aerts, C., & Rogers, T. M. 2015, , 806, 33
























[Arnett, Meakin & Young(2009)]amy09vel Arnett, W. D., Meakin, C., & Young, P. A., 2009, , 690, 1715 
















[Arnett, et al.(2015)]321D Arnett, W. D., Meakin, C. A., Viallet, M., Campbell, S. W., Lattanzio, J. C. & Mocák, M., 2015, , 809, 30 
[Arnett & Meakin(2016)]ropp Arnett, W. D. & Meakin, C. A., 2016, Reports on Progress in Physics, in press







[Bazàn & Arnett(1994)]ba94 Bazàn, G., & Arnett, D., 
1994, , 433, L41














[Böhm-Vitense(1958)]bv58 Böhm-Vitense, E., 1958, , 46, 108


































































[Constantino, et al.(2015)]tom15 Constantino, T., Campbell, S., Christensen-Dalsgaard, J., Lattanzio, J., Stello, D., 2015, , 452, 123












[Cristini, et al.(2015)]andrea2015 Cristini, A., Hirschi, R., Georgy, C., Meakin, C., Arnett, D., & Viallet, M., 2015, IAUS 307, 459 

[Cristini, et al.(2016)]andrea2016 Cristini, A., Hirschi, R., Georgy, C., Meakin, C., Arnett, D., & Viallet, M., 2016, submitted to 











[Drazin(2002)]drazinDrazin, P. G., 2002, Introduction to Hydrodynamic Stability, Cambridge University Press, Cambridge, U. K.















[Freytag, Ludwig, & Steffen(1996)]fls96 Freytag, B., Ludwig, H.-G., & Steffen, M., 1996, , 313, 497





[Frisch(1995)]frisch Frisch, U., 1995, Turbulence, Cambridge University Press, Cambridge





[Ghasemi, et al.(2016)]ghasemi16 Ghasemi, H., Moravveji, E., Aerts, C., et al. 2017, , 465, 1518



[Gough(1977)]gough77 Gough, D. O., 1977, 38th Coll., Problems of Stellar Convection (Berlin: Springer), 71, 799
 

























[Kitiashvili, et al.(2016)]kitiashvili Kitiashvili, I N., Kosovichev, A. G., Mansour, N. N., & Wray, A. A. 2016, , 821, 17

[Kolmogorov(1941)]kolmg41 Kolmogorov, A. N., 1941, Dokl.  Akad. Nauk SSSR, 30, 299

[Kolmogorov(1962)]kolmg Kolmogorov, A. N., 1962, J. Fluid Mech., 13, 82


[Landau & Lifshitz(1959)]llfm Landau, L. D. & Lifshitz, E. M., 1959, Fluid Mechanics,
Pergamon Press, London. 




























[Meakin & Arnett(2006)]ma06 Meakin, C. A., & Arnett, W. D., 2006, , 637, 53





[Meakin & Arnett(2007a)]ma07a Meakin, C. A., & Arnett, W. D., , 665, 690


[Meakin & Arnett(2007b)]ma07b Meakin, C. A., & Arnett, W. D., 2007b, , 667, 448





[Miesch(2005)]miesch Miesch, M. S., 2005, Living Reviews in Solar Physics, 2, 1









[Moravveji, et al.(2015)]ehsan1 Moravveji, E.,  Aerts, C., Pápics, et al. 2015, , 580, 27

[Moravveji, et al.(2016)]ehsan2 Moravveji, E., Townsend, R. H. D., Aerts, C., Mathis, S., 2016, , 823, 130

[Mosser et al.(2014)]mosser Mosser, B., Benomar, O., Belkacem, K., et al., 2014, , 572, 5







[Nordlund & Stein(1995)]ns95 Nordlund, A., & Stein, R., 1995, Stellar Evolution: What Should Be Done,  32nd Liège Int. Astroph. Coll., 32, 75




[Nordlund, Stein, & Asplund(2009)]nsa Nordlund, A., Stein, R., & Asplund, M., 2009,
  <http://www.livingreviews.org/lrsp-2009-2>






[Papics et al. (2014)]papics1 Pápics, P. I., Moravveji, E., Aerts, C., et al. 2014, , 570, 8

[Papics et al. (2015)]papics2 Pápics, P. I., Tkachenko, A., Aerts, C., et al. 2015, , 803, 25

































[Schindler, et al.(2015)]jt2015 Schindler, J.-T., Green, E. M., & Arnett, W. D., 2015, , 806, 178






[Smith & Arnett(2014)]nathan2014 Smith, Nathan, & Arnett, W. D., 2014, , 785, 82


































[Viallet, et al.(2013)]viallet2013 Viallet, M., Meakin, C., Arnett, D., Mocák, M., 2013, , 769, 1


[Viallet, et al.(2015)]viallet2015 Viallet, Maxime, Meakin, C., Prat, V., & Arnett, D., 2015, , 580, 61 




































[Zahn(1991)]zahn91 Zahn, J.-P. 1991, , 252, 179
[Zhang(2016)]zhang-2016-01 Zhang, Q. S. 2016, , 818, 146




