LOCAL RADIATION HYDRODYNAMIC SIMULATIONS OF MASSIVE STAR ENVELOPES AT THE IRON OPACITY PEAK

Yan-Fei Jiang(姜燕飞); Matteo Cantiello; Lars Bildsten; Eliot Quataert; Omer Blaes

doi:10.1088/0004-637X/813/1/74

1. INTRODUCTION

Massive stars play an important role in many astrophysical environments. The radiative and mechanical energy output from massive stars are important feedback mechanisms regulating star formation and the structure of the interstellar medium in galaxies (Kennicutt 2005; Smith 2014). Radiation from massive stars also plays an important role in the reionization of the universe (Bromm & Larson 2004). When massive stars explode, they produce various types of supernovae and high energy transients, while leaving behind black holes and neutron stars whose properties depend on the previous history of massive stellar evolution (Heger et al. 2003).

Because of the steep mass–luminosity relation of main sequence stars, the luminosity approaches the Eddington-limit for electron scattering in stars with masses larger than $\sim 50-100{M}_{\odot }$ (Crowther 2007; Maeder et al. 2012; Sanyal et al. 2015). In massive star envelopes with non-zero metallicity, the opacity increases significantly compared with electron scattering due to iron opacity at temperatures $\approx 1.8\times {10}^{5}$ K (e.g., Paxton et al. 2013 and references therein). As a result, the actual radiation force can locally exceed the gravitational force near the iron opacity peak. This same physics can operate at somewhat lower temperatures due to the opacity produced by helium and hydrogen (e.g., Maeder 1981).

The locally super-Eddington luminosity near the iron opacity peak (and related opacity peaks) has been a long standing challenge for one-dimensional (1D) stellar evolution models of massive stars. Hydrostatic and thermal equilibrium models with a super-Eddington luminosity require density and gas pressure inversions (Joss et al. 1973; Gräfener et al. 2012; Paxton et al. 2013; Owocki 2015), which are usually convectively unstable. When the iron opacity peak is sufficiently deep in the star, convection is efficient and can rearrange the stellar structure to have constant entropy, removing the density inversion predicted by radiative equilibrium models. When the iron opacity peak is close to the surface, however, convection is inefficient and cannot carry the necessary super-Eddington flux. It is unclear how the stellar structure adjusts in this limit, and in particular if the density inversions are stable (Maeder 2009; Sanyal et al. 2015). This results in very uncertain stellar radii in the very massive star regime (Yusof et al. 2013; Köhler et al. 2015). To complicate the situation even further, stellar evolution calculations of massive stars become numerically difficult when the opacity near the stellar photosphere significantly exceeds the electron scattering opacity. In this situation the convective energy transport is often artificially enhanced in 1D stellar evolution models (Maeder 1987; Paxton et al. 2013), which leads to the disappearance of the density inversions (Stothers & Chin 1979).

The density and temperature dependence of the opacity near composition-specific opacity peaks can also drive radiation hydrodynamic instabilities distinct from convection (Kiriakidis et al. 1993; Blaes & Socrates 2003; Suárez-Madrigal et al. 2013). In many cases, these cannot be captured by 1D stellar evolution models. These instabilities have been suggested as the origin of the large mass loss rates from massive stars that cannot be explained by line-driven winds (Glatzel & Kiriakidis 1993). The amount of mass loss by line driven winds (Puls et al. 2008) also depends on the structure of the stellar photosphere (e.g., the amplitude of density and velocity fluctuations, potentially affecting the degree of wind clumping).

In this paper, we present three-dimensional (3D) radiation hydrodynamic simulations of stellar envelopes near the iron opacity peak and quantify the impact of convective and radiation hydrodynamic instabilities on the stellar envelopes' properties and dynamics. These calculations provide the first direct calibration of classical mixing length theory (MLT; Cox 1968) in the radiation pressure-dominated regime. In addition, we quantify the extent to which convection and/or other hydrodynamic instabilities generate a "porous" atmosphere that substantially increases the radiation flux relative to that predicted by one-dimensional models. This is critical for determining whether locally super-Eddington regions in a stellar envelope can initiate a large-scale stellar wind via continuum radiation pressure or whether the stellar radiation freely escapes through low density channels (Shaviv 1998; van Marle et al. 2008).

Our work utilizes numerical techniques for accurate 3D radiation hydrodynamic simulations based on a variable Eddington tensor (VET) formalism (Davis et al. 2012; Jiang et al. 2012). These techniques have been successfully used to study a wide range of astrophysical problems (Jiang et al. 2013a, 2013b, 2014b) and have significant advantages relative to ad-hoc closures such as the flux-limited diffusion (FLD) approximation (Davis et al. 2014; Proga et al. 2014). The direct solution of the radiation transfer equation (via VET or Monte-Carlo methods) is particularly important for studying mildly optically thick regions near the stellar photosphere. Indeed, recent studies have shown that approximate radiation transfer schemes such as FLD or M1 can reach completely opposite conclusions compared with VET or Monte Carlo methods (Davis et al. 2014; Tsang & Milosavljevic 2015).

The remainder of this paper is organized as follows. In Section 2 we identify the key dimensionless parameter that governs the physics of near surface convection zones in stellar envelopes. We also describe the 1D stellar models that motivate the specific conditions in our 3D radiation hydrodynamic simulations. In Section 3 we describe our numerical methods while in Section 4 we present our results. Finally, in Section 5 we summarize our key results and discuss their implications.

2. THE MODEL PARAMETERS

The relative importance of convective and radiative energy transport can be characterized by a critical optical depth ${\tau }_{{\rm{c}}}$ over the pressure scale height H₀ near the iron opacity peak. As we will show based on our simulations, for conditions typical of massive stars, the maximum turbulent velocity reaches the isothermal sound speed when convection starts to become inefficient. Then the critical optical depth can be estimated by balancing the diffusive radiation flux $\sim {a}_{{\rm{r}}}{T}_{0}^{4}c/\tau$ and the maximal convective flux in the radiation pressure-dominated regime $\sim {c}_{{\rm{g}},0}{a}_{{\rm{r}}}{T}_{0}^{4}:$

$\begin{eqnarray}&&{\tau }_{{\rm{c}}}=c/{c}_{{\rm{g}},0}.\end{eqnarray} \tag{ 1 }$

We define the gas isothermal sound speed ${c}_{{\rm{g}},0},$ radiation sound speed ${c}_{{\rm{s}},0},$ pressure scale height H₀, optical depth ${\tau }_{0}$ and typical dynamic timescale t₀ based on the density ${\rho }_{0}$ and temperature T₀ at the iron opacity peak

$\begin{eqnarray}{c}_{{\rm{g}},0} & = & \sqrt{\displaystyle \frac{{P}_{0}}{{\rho }_{0}}},\\ {c}_{{\rm{s}},0} & = & \sqrt{\displaystyle \frac{{a}_{{\rm{r}}}{T}_{0}^{4}}{3{\rho }_{0}}},\\ {H}_{0} & = & \displaystyle \frac{{c}_{{\rm{s}},0}^{2}}{g}=\displaystyle \frac{{a}_{{\rm{r}}}{T}_{0}^{4}}{3{\rho }_{0}g},\\ {\tau }_{0} & = & {\kappa }_{t}({\rho }_{0},{T}_{0}){\rho }_{0}{H}_{0},\\ {t}_{0} & = & \displaystyle \frac{{H}_{0}}{{c}_{{\rm{s}},0}}.\end{eqnarray} \tag{ 2 }$

The characteristic radiation pressure is ${P}_{{\rm{r}},0}={a}_{{\rm{r}}}{T}_{0}^{4}/3$ while the gas pressure is P₀.

A related way to interpret the critical optical depth ${\tau }_{{\rm{c}}}$ in Equation (1) is that the typical thermal timescale due to photon diffusion t_c near the iron opacity peak is given by

$\begin{eqnarray}&&{t}_{{\rm{c}}}=\displaystyle \frac{{\tau }_{0}{H}_{0}}{c}=\displaystyle \frac{{\tau }_{0}}{{\tau }_{{\rm{c}}}}\displaystyle \frac{{c}_{{\rm{s}},0}}{{c}_{{\rm{g}},0}}{t}_{0}.\end{eqnarray} \tag{ 3 }$

When ${\tau }_{0}\gg {\tau }_{{\rm{c}}},$ the diffusion timescale is longer than the local dynamic timescale and we expect convection to be an efficient energy transport mechanism. When ${\tau }_{0}\ll {\tau }_{{\rm{c}}},$ however, radiative diffusion is rapid enough that convection may behave very differently compared with classical MLT. This can also be related to the Péclet number P_e, which is the ratio of advective to diffusive flux. When ${\tau }_{0}\gg {\tau }_{{\rm{c}}},$ we expect ${P}_{e}\gtrsim 1,$ while when ${\tau }_{0}\ll {\tau }_{{\rm{c}}},$ ${P}_{e}\lesssim 1.$ Note that there is some ambiguity as to whether ${\tau }_{{\rm{c}}}$ should be defined in terms of the gas isothermal sound speed ${c}_{{\rm{g}},0}$ or the radiation sound speed ${c}_{{\rm{s}},0}.$ As we shall show, we find that the former better characterizes the maximum velocities reached in our simulations.

In order to identify stellar conditions spanning the different regimes of radiation-dominated convection, we use Modules for Experiments in Stellar Evolution (MESA, Paxton et al. 2011, 2013, 2015 release r7385) to evolve stars with initial mass $40{M}_{\odot }$ and $80{M}_{\odot }.$ The models have an initial metallicity Z = 0.02 with a mixture taken from Asplund et al. (2005) and are calculated using the OPAL opacity tables (Iglesias & Rogers 1996). Convection is accounted for using MLT in the Henyey et al. (1965) formulation with ${\alpha }_{\mathrm{MLT}}=1.5.$ Convective boundaries are determined using the Ledoux criterion and we use semiconvection with an efficiency ${\alpha }_{\mathrm{sc}}$ = 0.01 (Langer et al. 1983, 1985). An exponentially decaying overshooting (Herwig 2000) has been included with ${\alpha }_{\mathrm{ov}}=0.001$ above and below convective boundaries. Stellar winds are included according to the mass loss scheme described in Glebbeek et al. (2009; Dutch wind scheme) with an efficiency parameter $\eta =0.8.$

We have identified specific times in the MESA massive star models with conditions that range from ${\tau }_{0}/{\tau }_{{\rm{c}}}\ll 1$ to ${\tau }_{0}/{\tau }_{{\rm{c}}}\gg 1.$ We label these models as StarDeep, StarMid and StarTop; they have ${\tau }_{0}/{\tau }_{{\rm{c}}}=15.2,0.97$ , and 0.025, respectively. Figure 1 shows the positions of these models in the HR diagram. Model StarDeep has an initial mass of $40{M}_{\odot }.$ It is evolved to an age of $4.3\times {10}^{6}$ year in MESA at which point it has an effective temperature $7.13\times {10}^{3}$ K, radius $431{R}_{\odot }$ and luminosity $8.96\times {10}^{5}{L}_{\odot }.$ Model StarMid starts from a $80{M}_{\odot }$ star. At an age of $3.02\times {10}^{6}$ year it has an effective temperature $1.17\times {10}^{4}$ K, radius $274{R}_{\odot }$ and luminosity $1.26\times {10}^{6}{L}_{\odot }.$ Model StarTop also starts from an $80{M}_{\odot }$ star but is only evolved to an age of $9.15\times {10}^{4}$ year, at which point it has an effective temperature $4.75\times {10}^{4}$ K, radius $14{R}_{\odot }$ and luminosity $8.96\times {10}^{5}{L}_{\odot }.$

**Figure 1.** Stellar evolution tracks of two solar metallicity stars with initial mass of $40{M}_{\odot }$ and $80{M}_{\odot }$ calculated with the 1D stellar evolution code MESA. The stellar symbols identify the locations corresponding to the three initial conditions of the 3D radiation hydro calculations with Athena listed in Table 1. The Athena calculations are effectively local plane-parallel atmosphere calculations motivated by the MESA models.
Download figure:
Standard image High-resolution image

**Figure 1.** Stellar evolution tracks of two solar metallicity stars with initial mass of $40{M}_{\odot }$ and $80{M}_{\odot }$ calculated with the 1D stellar evolution code MESA. The stellar symbols identify the locations corresponding to the three initial conditions of the 3D radiation hydro calculations with Athena listed in Table 1. The Athena calculations are effectively local plane-parallel atmosphere calculations motivated by the MESA models.
Download figure:
Standard image High-resolution image

As described in Section 3, we use plane-parallel 3D radiation hydrodynamics simulations to study the physics of massive stellar envelopes. The initial radial location of the iron opacity peak r₀ and the density ${\rho }_{0},$ temperature T₀ and luminosity L₀ at this location are taken from the corresponding MESA models. The fiducial parameters of our initial conditions are listed in Table 1.

Table 1. Simulation Parameters

Variables/Units	StarDeep	StarMid	StarTop
${r}_{0}/{R}_{\odot }$	125.7	65.6	13.6
T₀/K	$1.87\times {10}^{5}$	$1.67\times {10}^{5}$	$1.57\times {10}^{5}$
${\rho }_{0}/$ (g cm⁻³)	$3.10\times {10}^{-8}$	$3.60\times {10}^{-9}$	$5.52\times {10}^{-9}$
g/(cm s⁻²)	63.59	$3.56\times {10}^{2}$	$1.17\times {10}^{4}$
${F}_{{\rm{r}},i}/$ (erg cm⁻² s⁻¹)	$1.24\times {10}^{12}$	$9.94\times {10}^{12}$	$3.06\times {10}^{14}$
H₀/cm	$1.56\times {10}^{12}$	$1.52\times {10}^{12}$	$2.37\times {10}^{10}$
t₀/s	$1.56\times {10}^{5}$	$6.54\times {10}^{4}$	$1.42\times {10}^{3}$
${\tau }_{{\rm{c}}}\equiv c/{c}_{{\rm{g}},0}$	$5.99\times {10}^{3}$	$6.62\times {10}^{3}$	$6.54\times {10}^{3}$
${\tau }_{0}$	$9.12\times {10}^{4}$	$6.41\times {10}^{3}$	166.5
${P}_{{\rm{r}},0}/{P}_{0}$	3.96	26.32	13.22
${L}_{x,y}/{H}_{0}$	0.89	0.46	1.17
${L}_{z}/{H}_{0}$	6.26	3.71	4.60
${N}_{x,y}$	128	128	128
N_z	768	1024	512

Note. The gravitational acceleration g and radiation flux ${F}_{{\rm{r}},i}$ coming from the bottom of the simulation box are fixed. The optical depth ${\tau }_{0},$ temperature T₀, density ${\rho }_{0},$ radiation, and gas pressure ${P}_{{\rm{r}},0},$ P₀ are the fiducial parameters for each run. They are also the initial conditions around the iron opacity peak. The number of cells are ${N}_{x},{N}_{y}$ , and N_z.

Download table as: ASCII Typeset image

3. NUMERICAL METHOD

3.1. Equations

We model the envelopes of massive stars in plane parallel geometry, neglecting the global stellar geometry. This is suitable for studying the effects of 3D radiation hydrodynamics on energy transport in the stellar envelope, but is not appropriate for studying global effects such as stellar outflows. We solve the frequency-independent (gray) radiation hydrodynamic equations in cartesian coordinates $(x,y,z)$ with unit vectors ( $\hat{x},$ $\hat{y},$ $\hat{z}$ ). The gravitational acceleration g is assumed to be a constant and along the $-\hat{z}$ direction. As the mass of the envelope is less than 1% of the core mass, this is a reasonable assumption (the slight change in $1/{r}^{2}$ in the stellar envelope is also neglected). The equations are (e.g., Jiang et al. 2012, 2013a)

$\begin{eqnarray}\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}) & = & 0,\\ \displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{vv}}+P{\mathsf{I}}) & = & -{S}_{{\rm{r}}}({\boldsymbol{P}})-\rho g\hat{z},\\ \displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\cdot \left[(E+P){\boldsymbol{v}}\right] & = & -{{cS}}_{{\rm{r}}}(E)-g\rho {\boldsymbol{v}}\cdot \hat{z},\\ \displaystyle \frac{\partial {E}_{{\rm{r}}}}{\partial t}+{\rm{\nabla }}\cdot {{\boldsymbol{F}}}_{{\rm{r}}} & = & {{cS}}_{{\rm{r}}}(E),\\ \displaystyle \frac{1}{{c}^{2}}\displaystyle \frac{\partial {{\boldsymbol{F}}}_{{\rm{r}}}}{\partial t}+{\rm{\nabla }}\cdot {{\mathsf{P }}}_{{\rm{r}}} & = & {{\boldsymbol{S}}}_{{\rm{r}}}({\boldsymbol{P}}),\end{eqnarray} \tag{ 4 }$

where the radiation source terms are

$\begin{eqnarray}{{\boldsymbol{S}}}_{{\rm{r}}}({\boldsymbol{P}}) & = & -\rho \left({\kappa }_{\mathrm{aF}}+{\kappa }_{\mathrm{sF}}\right)\left[{{\boldsymbol{F}}}_{{\rm{r}}}-\left({\boldsymbol{v}}{E}_{{\rm{r}}}+{\boldsymbol{v}}\cdot {{\mathsf{P }}}_{{\rm{r}}}\right)\right]/c\\ & & +\rho {\boldsymbol{v}}({\kappa }_{\mathrm{aP}}{a}_{{\rm{r}}}{T}^{4}-{\kappa }_{\mathrm{aE}}{E}_{{\rm{r}}})/c,\\ {S}_{{\rm{r}}}(E) & = & \rho ({\kappa }_{\mathrm{aP}}{a}_{{\rm{r}}}{T}^{4}-{\kappa }_{\mathrm{aE}}{E}_{{\rm{r}}})\\ & & +\rho ({\kappa }_{\mathrm{aF}}-{\kappa }_{\mathrm{sF}})\displaystyle \frac{{\boldsymbol{v}}}{{c}^{2}}\cdot \left[{{\boldsymbol{F}}}_{{\rm{r}}}-\left({\boldsymbol{v}}{E}_{{\rm{r}}}+{\boldsymbol{v}}\cdot {{\mathsf{P }}}_{{\rm{r}}}\right)\right].\end{eqnarray} \tag{ 5 }$

Notice that local thermal equilibrium is assumed for the emission term.

In the above equations, $\rho ,P,{\boldsymbol{v}},c$ are the gas density, pressure, flow velocity and speed of light respectively. The total gas energy density is $E={E}_{{\rm{g}}}+\rho {v}^{2}/2,$ where ${E}_{{\rm{g}}}=P/(\gamma -1)$ is the internal gas energy density with a constant adiabatic index $\gamma =5/3.$ The gas pressure is $P=\rho {k}_{{\rm{B}}}T/\mu ,$ where ${k}_{{\rm{B}}}$ is Boltzmann's constant and $\mu =0.62{m}_{p}$ is the mean molecular weight for nearly fully ionized gas with proton mass m_p. The radiation constant is ${a}_{{\rm{r}}}=7.57\times {10}^{15}$ erg cm⁻³ K⁻⁴, while ${E}_{{\rm{r}}}\mathrm{and}{{\boldsymbol{F}}}_{{\rm{r}}}$ are the radiation energy density and flux. The Rosseland mean absorption and scattering opacities are denoted by ${\kappa }_{\mathrm{aF}}$ and ${\kappa }_{\mathrm{sF}},$ while ${\kappa }_{\mathrm{aP}}$ and ${\kappa }_{\mathrm{aE}}$ are the Planck and energy mean absorption opacities. We describe our choice of these opacities in the following section.

Equations (4) and (5) are solved in the mixed frame, meaning that the radiation flux ${{\boldsymbol{F}}}_{{\rm{r}}}$ and energy density E_r are Eulerian variables, while the material–radiation interaction terms in Equations (5) are written in the co-moving frame (e.g., Lowrie et al. 1999). The Eulerian and co-moving flux ${{\boldsymbol{F}}}_{{\rm{r}},0}$ are related through ${{\boldsymbol{F}}}_{{\rm{r}},0}={{\boldsymbol{F}}}_{{\rm{r}}}-\left({\boldsymbol{v}}{E}_{{\rm{r}}}+{\boldsymbol{v}}\cdot {{\mathsf{P }}}_{{\rm{r}}}\right).$ Notice that only the co-moving radiation flux contributes to the radiation acceleration while the advection flux does not.

We relate the radiation pressure ${{\mathsf{P }}}_{{\rm{r}}}$ and radiation energy density through a VET, ${{\mathsf{P }}}_{{\rm{r}}}={\mathsf{P}}{E}_{{\rm{r}}}.$ We do not use an assumed closure relation (such as FLD or M1) to specify ${\mathsf{P}}.$ Instead, we compute it directly from angular quadratures of the specific intensity I_r for each angle ${\boldsymbol{n}},$ which is calculated from the time-independent radiation transfer equation

$\begin{eqnarray}&&\displaystyle \frac{\partial {I}_{{\rm{r}}}}{\partial s}={\kappa }_{t}\left(S-{I}_{{\rm{r}}}\right),\end{eqnarray} \tag{ 6 }$

where S is the radiation source term and ${\kappa }_{t}$ is the total specific opacity. Because the typical dynamic timescale is much longer than the light crossing time, f barely changes between each time step. Our use of the time independent transfer equation to calculate the angular distribution of the radiation field is thus very accurate. We solve the transfer equation in Equation (6) using the method of short characteristics, as described in Davis et al. (2012). Then f is calculated according to its definition

$\begin{eqnarray}&&{\mathsf{P}}=\displaystyle \frac{\int {I}_{{\rm{r}}}{\boldsymbol{nn}}d{\rm{\Omega }}}{\int {I}_{{\rm{r}}}d{\rm{\Omega }}},\end{eqnarray} \tag{ 7 }$

where Ω is the solid angle. We emphasize that we only use the time-independent transfer Equation (6) to calculate the Eddington tensor ${\mathsf{P}},$ which represents the angular distribution of the radiation field. The moments of I_r are not used to calculate the radiation-material interactions. Instead, the radiation moments E_r, ${{\boldsymbol{F}}}_{{\rm{r}}}$ and the radiation source terms are calculated self-consistently through the time dependent radiation moment Equations (4).

The Eddington tensor in Equation (7) closes the radiation moment equations (Equations (4)). We solve these equations using the Godunov radiation MHD code Athena, as described and tested in Jiang et al. (2012; and with additional improvements described in Jiang et al. 2013c). We solve the radiation subsystems as well as the radiation source terms implicitly, while the hydro equations are solved explicitly as in the original Athena code (Stone et al. 2008), with appropriate modifications due to the stiff radiation source terms (Jiang et al. 2012).

3.2. The Opacity

In order to correctly capture the temperature and density dependencies of the opacity, particularly near the iron opacity peak, we directly use the opacity table from MESA (Paxton et al. 2011, 2013, 2015) in our Athena simulations. We assume a metallicity Z = 0.02 and hydrogen fraction X = 0.6. For the density ρ and temperature T in each cell, the total specific opacity ${\kappa }_{t}$ is calculated with bilinear interpolation based on the table. This is used as the total specific flux mean opacity ${\kappa }_{t}={\kappa }_{\mathrm{aF}}+{\kappa }_{\mathrm{sF}}.$ In order to split this into scattering and absorption opacities, we assume a constant electron scattering opacity ${\kappa }_{\mathrm{sF}}=0.32$ g⁻¹ cm². All of the absorption opacities ( ${\kappa }_{\mathrm{aF}},{\kappa }_{\mathrm{aP}},{\kappa }_{\mathrm{aE}}$ ) are assumed to be the same for simplicity.

Figure 2 shows the actual opacity ${\kappa }_{t}$ we use for a range of temperatures and densities appropriate to stellar envelopes. The opacity peak between 10⁵ and 2 × 10⁵ K is due to Fe, which is the region we focus on. The increase in opacity at lower temperatures is due to He and H. The iron opacity peak has a sensitive dependence on temperature but a weak dependence on density. Particularly when $\rho \lesssim 3.6\times {10}^{-11}$ g cm⁻³, ${\kappa }_{t}$ does not depend on density for $T\gtrsim {10}^{5}$ K. For the stellar envelopes we study, the temperature usually does not drop below 4 × 10⁴ K. Moreover, when this does happen, the density is always smaller than 10⁻¹⁰ g cm⁻³. As a result, the He and H opacity peaks usually do not show up in our calculations.

3.3. Initial and Boundary Conditions

We set our initial conditions in the Athena simulations by using the MESA models described in Section 2 to determine the initial radial location r₀, density ${\rho }_{0}$ and temperature T₀ at the iron opacity peak. The gravitational acceleration g is also chosen to be the value at r₀ in the MESA models. We then apply a constant radiation flux ${F}_{{\rm{r}},i}$ through the whole simulation box, which represents the radiation coming from the center of the massive star. The initial radiation flux is determined by the total luminosity L₀ in the MESA models at r₀: ${F}_{{\rm{r}},i}={L}_{0}/(4\pi {r}_{0}^{2}).$

With these fiducial parameters and the opacity given in Section 3.2, we construct our initial conditions by solving for a hydrostatic envelope in thermal equilibrium. Starting from $z={r}_{0},$ all the fluid and radiation quantities are integrated as a function of height using

$\begin{eqnarray}\displaystyle \frac{\partial P}{\partial z} & = & \displaystyle \frac{\rho {\kappa }_{t}}{c}{F}_{{\rm{r}},i}-\rho g,\\ {E}_{{\rm{r}}} & = & {a}_{{\rm{r}}}{T}^{4},\\ \displaystyle \frac{1}{3}\displaystyle \frac{\partial {E}_{{\rm{r}}}}{\partial z} & = & -\displaystyle \frac{\rho {\kappa }_{t}}{c}{F}_{{\rm{r}},i}.\end{eqnarray} \tag{ 8 }$

We integrate above and below r₀ to capture the whole profile of the iron opacity peak region in our simulation box. We then add random perturbations with $0.01\%$ amplitude to the density in the simulation box of dimensions ${L}_{x},{L}_{y},{L}_{z}.$ The initial conditions constructed in this way are slightly different from the 1D profiles returned by MESA models around the iron opacity peak. We only guarantee that the characteristic physical parameters we choose are consistent with the MESA stellar models. In particular, our initial condition does not have any mixing length model of convective or advective energy transport. We then calculate how the energy will be transported using self-consistent 3D radiation hydrodynamic simulations.

At the bottom of the box, we use a reflecting boundary condition for gas quantities. This means that the density, opacity, and the horizontal components of the velocity are copied from the last active zones to the ghost zones used for boundary conditions. The vertical component of velocity in the ghost zones is set to be the opposite of the velocity in the last active zones. The radiation flux is fixed at the constant value ${F}_{{\rm{r}},i}$ while the radiation energy density E_r is integrated to the ghost zones using the diffusion equation. For the top boundary, all gas quantities and the radiation flux ${{\boldsymbol{F}}}_{{\rm{r}}}$ are copied from the last active zones to the ghost zones. When the photosphere is included in the simulation box (StarTop & StarMid), the radiation energy density E_r in the ghost zones is also copied from the last active zones. When the whole simulation box is optically thick (StarDeep), E_r is extended to the ghost zones using the diffusion equation. Periodic boundary conditions are used for all quantities in the horizontal direction.

For the short characteristic module we use to calculate the VET, the incoming specific intensities from the bottom of the domain are chosen such that angular quadratures of the intensities give a constant vertical flux ${F}_{{\rm{r}},i}$ and the same E_r as in the moment equations at the bottom. The incoming specific intensities at the top of the domain are set to zero when the photosphere is inside the simulation box. By contrast, when the whole simulation box is optically thick, the incoming specific intensities from the top of the domain are set to be the same as the outgoing specific intensities.

4. RESULTS

4.1. The Run StarDeep when ${\tau }_{0}\gg {\tau }_{c}$

The run StarDeep is in the regime where the gas and radiation are tightly coupled and we expect efficient convection. Below, we first show how the simulated envelope evolves, then summarize its time averaged structure, and finally compare the energy transport in our 3D radiation hydrodynamic simulation with MLT.

4.1.1. Simulation History

Three snapshots of density in Figure 3 show the evolution of the envelope. The initial density varies only vertically. At the bottom where the temperature is 30% larger than T₀, the radiation acceleration is smaller than the gravitational acceleration and both density and temperature decrease with height. Around $z=130{R}_{\odot },$ where we enter the temperature range of the iron opacity peak, the radiation acceleration becomes larger than the gravitational acceleration for the fixed radiation flux ${F}_{{\rm{r}},i},$ and the density increases with height, peaking around $z=160{R}_{\odot }.$ Because the temperature always decreases with height according to the diffusion equation, the radiation acceleration eventually decreases after the iron opacity peak. When the radiation flux again becomes sub-Eddington, the density drops with height. The end result is that the initial condition has the density inversion associated with hydrostatic envelopes with a super-Eddington radiation flux (Joss et al. 1973). This profile is our initial condition for the simulation, and does not represent the profile in the MESA 1D calculation.

**Figure 3.** Snapshots of density for the run StarDeep at time 0 (left), $21.4{t}_{0}$ (middle), and $97.8{t}_{0}$ (right). The box sizes labeled in the plots are in units of solar radius ${R}_{\odot }.$ The horizontal box size is ${L}_{x}={L}_{y}=20{R}_{\odot }.$ Density is in units of ${\rho }_{0},$ which is given in Table 1 for this run.
Download figure:
Standard image High-resolution image

The middle panel of Figure 3 shows that the initial density inversion is unstable. By $t\approx 20{t}_{0},$ there are large rising and sinking plumes characteristic of convection and/or the Rayleigh–Taylor instability. We now show how this instability can be quantitatively understood as arising from convection in the radiation-dominated regime.

The gas and radiation entropy per unit mass ${S}_{{\rm{g}}},{S}_{{\rm{r}}}$ are given by

$\begin{eqnarray}{S}_{{\rm{g}}} & \equiv & \displaystyle \frac{{k}_{{\rm{B}}}}{\mu (\gamma -1)}\mathrm{ln}\left[\displaystyle \frac{P/{P}_{0}}{{\left(\rho /{\rho }_{0}\right)}^{\gamma }}\right]\ \ \ \mathrm{and}\\ {S}_{{\rm{r}}} & \equiv & \displaystyle \frac{4{E}_{{\rm{r}}}}{3\rho T}.\end{eqnarray} \tag{ 9 }$

In our simulations, the total entropy is dominated by the radiation. This increases with height near the bottom of the simulation box in the initial condition, but necessarily decreases with height near the density inversion where T decreases but ρ increases (see Figure 7 discussed below). We thus expect that our initial condition is convectively stable near the bottom of the domain, but convectively unstable near the density inversion region in 3D. This is exactly what we find in the simulation (see Figure 3). Convection happens around $160{R}_{\odot }$ within $\approx 20{t}_{0}$ (Equation (3)) and mixes the initial density inversion. High density fingers sink and low density gas rises upwards. By $t=97.8{t}_{0}$ (third panel of Figure 3), the bulk of the stellar envelope has reached a statistically steady turbulent state driven by convection.

The initial radiation Brunt–Väisälä frequency $| {N}_{{\rm{r}}}|$ (e.g., Equation (52) of Blaes & Socrates (2003)) at the iron opacity peak is $0.8/{t}_{0},$ but it takes $\approx 20{t}_{0}$ for convection to begin destroying the initial density inversion. The wavelength of the fastest growing initial mode is about half of the horizontal box size, which also differs from the normal Rayleigh–Taylor instability with a density discontinuity (Jiang et al. 2013a) where the shortest wavelengths grow fastest. To understand the linear growth phase of the instability, we first calculate the two-dimensional (2D) discrete Fourier transform in the $x,y$ plane of the density at height $z=160{R}_{\odot }$ as a function of time. The 2D Fourier coefficient is then binned into a 1D power spectrum ${kf}(\rho ),$ where k is the magnitude of the horizontal wavenumber. The time evolution of two modes with horizontal wavelengths $0.41{H}_{0}=9.2{R}_{\odot }$ and $0.18{H}_{0}=4.0{R}_{\odot }$ are shown in Figure 4. The decrease of ${kf}(\rho )$ within t₀ probably originates from an initial vertical oscillation that is present due to the imperfect balance between gravity and pressure gradients in the initial condition. After this point, the long wavelength mode grows exponentially with an e-folding time $3.85{t}_{0}$ until $15{t}_{0}.$ The shorter wavelength mode also grows exponentially, but at a rate that is smaller by a factor of 3.44 initially. However, it then grows quickly after 15t₀ as convection becomes nonlinear and power from long wavelength modes cascades to the short wavelength modes.

**Figure 4.** Time evolution of the density power spectrum ${kf}(\rho )$ for the run StarDeep at height $z=160{R}_{\odot }$ during the linear growth phase of convection. The black and red lines are for the modes with wavelength $0.41{H}_{0}$ and $0.18{H}_{0},$ respectively. The dashed black line is a fit of an exponential function to the solid black line, which indicates an e-folding time $3.85{t}_{0}$ for this mode. The power spectrum is normalized by the mean density at $z=160{R}_{\odot },$ such that the integral of ${k}^{2}{f}^{2}(\rho )$ over all wavenumbers is unity.
Download figure:
Standard image High-resolution image

**Figure 4.** Time evolution of the density power spectrum ${kf}(\rho )$ for the run StarDeep at height $z=160{R}_{\odot }$ during the linear growth phase of convection. The black and red lines are for the modes with wavelength $0.41{H}_{0}$ and $0.18{H}_{0},$ respectively. The dashed black line is a fit of an exponential function to the solid black line, which indicates an e-folding time $3.85{t}_{0}$ for this mode. The power spectrum is normalized by the mean density at $z=160{R}_{\odot },$ such that the integral of ${k}^{2}{f}^{2}(\rho )$ over all wavenumbers is unity.
Download figure:
Standard image High-resolution image

The measured growth rates agree very well with the short wavelength WKB linear analysis done by Blaes & Socrates (2003; the last factor of their Equation (59)). Plugging in numbers for the long wavelength mode in Figure 4, and assuming that the vertical wavelength of this mode is about equal to its horizontal wavelength, we get a growth rate of $0.27/{t}_{0},$ very close to the value measured in the simulation. This decline should only be valid for wavelengths with growth rates less than the magnitude of the Brunt–Väisälä frequency $| {N}_{{\rm{r}}}| ,$ giving an expected transition wavelength as defined in Blaes & Socrates (2003) of roughly $0.5{H}_{0}$ for this simulation. Hence we expect the longest wavelength convective modes to be only weakly affected by radiative damping, which explains why the measured long wavelength growth rate of $0.26/{t}_{0}$ is of the same order as $| {N}_{{\rm{r}}}| .$ Shorter wavelength modes will be more strongly affected, and grow much more slowly. The measured growth rate of the shorter wavelength mode here is consistent with a reduction factor $\propto {k}^{-1.5}.$ This is consistent with the expected behavior ${k}^{-2},$ given that we are not too far from the transition wavelength, and that the short wavelength mode could have a more horizontal wave vector orientation.

The envelope of the StarDeep simulation shows significant vertical oscillatory motion. We calculate two measures of the time-dependent vertical velocity, weighted by energy and by density,

$\begin{eqnarray}&&{V}_{z,{E}_{{\rm{r}}}}=\displaystyle \frac{\int {v}_{z}{E}_{{\rm{r}}}{dV}}{\int {E}_{{\rm{r}}}{dV}},{V}_{z,\rho }=\displaystyle \frac{\int {v}_{z}\rho {dV}}{\int \rho {dV}},\end{eqnarray} \tag{ 10 }$

where the integrals are done over the whole simulation box. Histories of these two vertical velocity measures are shown in Figure 5. The averaged density, or equivalently the total mass, drops by 9% around $20{t}_{0}.$ At the same time, ${V}_{z,{E}_{{\rm{r}}}}$ and ${V}_{z,\rho }$ reach their peak values. This corresponds to the time when the initial density inversion is destroyed due to convection. After the initial transient, the envelope slowly settles down to a steady state, and both ${V}_{z,{E}_{{\rm{r}}}}$ and ${V}_{z,\rho }$ show oscillations with amplitudes of 1%–2% of the radiation sound speed. The average density also slowly decreases with time because of the mass loss through our open top boundary. We have run the simulation for more than five thermal times to reach the thermal equilibrium.

**Figure 5.** History of the radiation energy density weighted ( ${V}_{z,{E}_{{\rm{r}}}},$ solid line) and density weighted ( ${V}_{z,\rho },$ dashed line) vertical velocity for the run StarDeep.
Download figure:
Standard image High-resolution image

4.1.2. Time Averaged Vertical Structures

In the 3D radiation hydrodynamic simulations, energy can be transported vertically by either photon diffusion or turbulent transport. We quantify the turbulent transport using the advection flux ${F}_{\mathrm{adv}},$ which is

$\begin{eqnarray}&&{F}_{\mathrm{adv}}={v}_{z}{E}_{{\rm{r}}}.\end{eqnarray} \tag{ 11 }$

In principle, photon advection can be due to either a mean flow (e.g., an outflow) or turbulence. In our simulations, the latter dominates.

The solid line in the top panel of Figure 6 shows the time averaged advective radiation flux ${F}_{\mathrm{adv}}$ as a function of height where the time average is done between 31t₀ and $107{t}_{0}.$ Most of the flux at $z\lt 100{R}_{\odot }$ is carried by radiative diffusion. In the turbulent region between $z=110{R}_{\odot }$ and $160{R}_{\odot },$ however, the averaged advection flux is more than 40% of ${F}_{{\rm{r}},i}.$ We compare this to MLT in the next section.

**Figure 6.** Top: horizontally and time averaged advection flux ${F}_{\mathrm{adv}}={v}_{z}{E}_{{\rm{r}}}$ (black line) for the run StarDeep. The red line is the convection flux calculated according to Equation (17) based on the horizontally and time averaged vertical profiles of temperature and pressure in the same run with mixing length parameter $\alpha =0.55.$ Middle: horizontally and time averaged turbulent velocity $\bar{V}$ (Equation (12)), scaled with radiation sound speed c_s (red line) and isothermal sound speed c_g (black line) at each height. Bottom: vertical profile of $\tilde{{\rm{\nabla }}}\equiv ({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}})/{{\rm{\nabla }}}_{\mathrm{ad}}$ for the run StarDeep. The actual gradient ${\rm{\nabla }}\equiv d\mathrm{ln}T/d\mathrm{ln}(P+{P}_{{\rm{r}}})$ is calculated based on the horizontally and time averaged temperature and pressure profiles. The adiabatic gradient is close to 0.25 for the radiation-dominated envelope. The fact that $\tilde{{\rm{\nabla }}}\ll 1$ implies that the convection is almost adiabatic.
Download figure:
Standard image High-resolution image

The turbulent velocity in the envelope can be quantified as

$\begin{eqnarray}&&\bar{V}=\sqrt{\displaystyle \frac{\int \rho {\boldsymbol{v}}\cdot {\boldsymbol{v}}{dV}}{\int \rho {dV}}}.\end{eqnarray} \tag{ 12 }$

The middle panel of Figure 6 shows the time averaged vertical profile of this turbulent velocity relative to both the isothermal sound speed ${c}_{{\rm{g}}}=\sqrt{P/\rho }$ and the radiation sound speed ${c}_{{\rm{s}}}=\sqrt{{E}_{{\rm{r}}}/(3\rho )}.$ The averaged turbulent velocity is only 6%–20% of the radiation sound speed in the turbulent region. Even compared with the isothermal sound speed alone, $\bar{V}/{c}_{{\rm{g}}}$ is smaller than 0.4. This implies that in the efficient convection regime, the turbulent kinetic energy density is much smaller than the thermal energy density ( $\lesssim 2\%$ ), as expected.

Figure 7 shows the quasi-steady envelope structure compared with the hydrostatic initial condition. Note that the density and temperature are relatively unaffected near the bottom of the box indicating that we have placed the bottom boundary deep enough. In the final quasi-steady envelope, the iron opacity peak moves slightly deeper relative to the initial condition, but the density inversion present in the initial condition is gone. This is consistent with the 1D MESA models which also find that for the conditions in StarDeep MLT predicts that convection is efficient and is able to remove the density inversion. Note that the location where the advection flux is largest ( $z\approx 120{R}_{\odot }$ in Figure 6) corresponds to the location of the iron opacity peak in the final state. This is also where the radiation entropy S_r decreases with height, consistent with convection continuing to drive the turbulent motions and energy transport. As expected, the radiation entropy profile is also flatter in the turbulent state after the density inversion is removed.

**Figure 7.** Horizontally and time averaged vertical profiles of various quantities for the run StarDeep. Panel (a): density (solid line at top), gas temperature (solid black line at bottom) and radiation temperature (solid red line at bottom). The gas and radiation temperature profiles overlap as they are strongly thermally coupled. Panel (b): gas entropy S_g (solid black line at top) and radiation entropy S_r (solid black line at bottom). These entropies are calculated according to Equation (9), and both are in units of ${k}_{{\rm{B}}}/\mu .$ Panel (c): total opacity ${\kappa }_{t}$ in unit of $1/({\rho }_{0}{H}_{0})$ (the solid line). The dashed black lines in these three panels are the initial conditions. Panel (d): volume averaged (solid black line) and density weighted (dashed black line) radiation accelerations (a_r, ${\tilde{a}}_{{\rm{r}}}$ ) as well as the acceleration due to the gas pressure gradient (a_g, green line). The red line is the sum of ${\tilde{a}}_{{\rm{r}}}$ and a_g. In this run with ${\tau }_{0}\gg {\tau }_{{\rm{c}}},$ the density inversion in the initial condition is gone.
Download figure:
Standard image High-resolution image

Because part of the energy flux in the turbulent state is carried by advection, the averaged diffusive flux also drops below the initial value ${F}_{{\rm{r}},i}.$ The significant contribution of the advection flux reduces the radiation acceleration so that it is smaller than the gravitational acceleration. The accelerations due to gas a_g and radiation a_r can be calculated as

$\begin{eqnarray}&&{a}_{{\rm{g}}}=-\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial P}{\partial z},\\ &&{a}_{{\rm{r}}}=\displaystyle \frac{{\kappa }_{t}}{c}{F}_{{\rm{r}},0z}.\end{eqnarray} \tag{ 13 }$

Notice that only the diffusive radiation flux ${F}_{{\rm{r}},0z}$ contributes to the radiation acceleration, which is the radiation pressure gradient in the optically thick regime. The last panel of Figure 7 shows that gravity is roughly balanced by the sum of radiation and gas pressure gradients, with a_r being at most 90% of g. The fact that ${a}_{{\rm{r}}}\lesssim g$ is consistent with the absence of a density inversion in the turbulent state.

4.1.3. Time Averaged Horizontal Structures

Figure 8 shows horizontal slices of density ρ and the vertical component of Eulerian radiation flux ${F}_{{\rm{r}},z}$ at $z=140{R}_{\odot }$ and time $97.8{t}_{0}$ for the run StarDeep. There is a clear anti-correlation between fluctuations of ρ and ${F}_{{\rm{r}},z}.$ When ρ is larger (smaller) than the horizontally averaged density, ${F}_{{\rm{r}},z}$ is negative (positive). This is a natural outcome of convection, where the Eulerian radiation flux ${F}_{{\rm{r}},z}$ is dominated by the advection part ${v}_{z}{E}_{{\rm{r}}}$ at each location in the regime ${\tau }_{0}\gg {\tau }_{{\rm{c}}}.$ Figure 8 thus shows that low density regions buoyantly rise while high density fluid elements fall.

To quantify the level of density and temperature fluctuations, the top panel of Figure 9 shows the vertical profiles of standard deviations of ρ and T scaled with the horizontally averaged mean values at each height at time $97.8{t}_{0}.$ The scaled standard deviations peak at $150{R}_{\odot },$ where the iron opacity peak is located at this time. Beyond the convective region, fluctuations drop with height. The scaled density standard deviation is always smaller than 10% in this case. Temperature fluctuations ( $\lesssim 0.4\%$ ) are much smaller. Because radiation pressure ${P}_{{\rm{r}}}\propto {T}^{4},$ the pressure fluctuations are only $\lesssim 1.6\%.$

**Figure 9.** Top: vertical profiles of standard deviations of density ρ and temperature T at time $97.8{t}_{0}$ for the run StarDeep. The standard deviations are scaled with the mean density and temperature at each height. Bottom: vertical profiles of cross correlation coefficients between ρ, T and ρ, ${F}_{{\rm{r}},z}$ at the same time. The anti-correlation between ρ and ${F}_{{\rm{r}},z}$ in this run is due to the buoyancy instead of the "porosity" effect, as ${F}_{{\rm{r}},z}$ is dominated by the advection flux ${v}_{z}{E}_{{\rm{r}}}.$
Download figure:
Standard image High-resolution image

The cross correlation between two quantities a and b can be quantified by the correlation coefficient ${C}_{a,b},$ normalized by their standard deviations ${\sigma }_{a},{\sigma }_{b}.$ The bottom panel of Figure 9 shows the correlation coefficients between $\rho ,T$ and $\rho ,{F}_{{\rm{r}},z}.$ Consistent with the slice shown in Figure 8, there is an anti-correlation between ρ and ${F}_{{\rm{r}},z}$ in the convective region due to buoyancy. Temperature fluctuations also anti-correlate with density fluctuations, which means an anti-correlation between density ρ and radiation pressure P_r. In this very optically thick regime ${\tau }_{0}\gg {\tau }_{{\rm{c}}},$ density fluctuations do not reduce the effective radiation acceleration, which can be confirmed by comparing the volume averaged a_r (Equation (13)), and the density weighted radiation acceleration

$\begin{eqnarray}&&{\tilde{a}}_{{\rm{r}}}=\displaystyle \frac{\langle \rho {\kappa }_{t}{F}_{{\rm{r}},0z}\rangle }{c\langle \rho \rangle }.\end{eqnarray} \tag{ 14 }$

Here the average $\langle \cdot \rangle$ is done horizontally at each height. Panel (d) of Figure 7 shows vertical profiles of the time averaged radiation accelerations calculated in these two different ways. They are almost identical. This explictly demonstrates that there is no reduction of the radiation acceleration due to density fluctuations in this regime.

The typical size of the turbulent structures can be estimated by calculating the normalized cross correlation ${C}_{\rho ,\tilde{\rho }}$ between $\rho (x,y)$ and $\rho (x+{l}_{\rho },y+{l}_{\rho }),$ where $\rho (x+{l}_{\rho },y+{l}_{\rho })$ is the density shifted horizontally with distance ${l}_{\rho }.$ The correlation length ${l}_{\rho }$ can be defined to to be the shifted distance when ${C}_{\rho ,\tilde{\rho }}$ drops to zero. For the run StarDeep, ${l}_{\rho }$ is typically $\sim 0.3-0.4{H}_{0},$ with corresponding optical depth at the iron opacity peak comparable to $c/{c}_{{\rm{s}},0}$ given in Table 1.

4.1.4. Comparison with MLT

In the radiation-dominated regime, the radiation advection flux ${F}_{\mathrm{adv}}$ from the simulations can be directly compared with the convection flux from MLT. We define the ratio of gas pressure to total pressure (β), the first ( ${{\rm{\Gamma }}}_{1}$ ) and second ( ${{\rm{\Gamma }}}_{2}$ ) adiabatic indices, and specific heat at constant pressure (c_p) for a mixture of gas and radiation (Chandrasekhar 1967), and the logarithmic derivative of density with respect to temperature at constant total pressure to be

$\begin{eqnarray}\beta & \equiv & \displaystyle \frac{P}{P+{P}_{{\rm{r}}}},\\ {{\rm{\Gamma }}}_{1} & \equiv & \beta +\displaystyle \frac{{(4-3\beta )}^{2}(\gamma -1)}{\beta +12(\gamma -1)(1-\beta )},\\ {{\rm{\Gamma }}}_{2} & = & \displaystyle \frac{(4-3\beta ){{\rm{\Gamma }}}_{1}}{3(1-\beta ){{\rm{\Gamma }}}_{1}+\beta }\\ {c}_{p} & = & \displaystyle \frac{1}{\gamma -1}\displaystyle \frac{{{\rm{\Gamma }}}_{1}}{{\beta }^{2}}\left[\beta +12(\gamma -1)(1-\beta )\right],\\ Q & \equiv & -{\left.\displaystyle \frac{d\mathrm{ln}\rho }{d\mathrm{ln}T}\right|}_{P+{P}_{{\rm{r}}}}=\displaystyle \frac{4-3\beta }{\beta },\end{eqnarray} \tag{ 15 }$

where P_r is the zz component of the radiation pressure tensor ${{\mathsf{P }}}_{{\rm{r}}}.$ The adiabatic temperature gradient with respect to total pressure can be calculated as

$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{ad}}\equiv {\displaystyle \frac{d\mathrm{ln}T}{d\mathrm{ln}(P+{P}_{{\rm{r}}})}|}_{\mathrm{ad}}=\displaystyle \frac{{{\rm{\Gamma }}}_{2}-1}{{{\rm{\Gamma }}}_{2}}.\end{eqnarray} \tag{ 16 }$

The actual gradient ${\rm{\nabla }}\equiv d\mathrm{ln}T/d(\mathrm{ln}(P+{P}_{{\rm{r}}}))$ can be calculated based on the mean density and temperature profiles shown in Figure 7. If we assume that the mixing length is related to the pressure scale height via the mixing parameter α, then the convective flux based on non-adiabatic MLT including the effects of radiative cooling is (Henyey et al. 1965)

$\begin{eqnarray}&&{F}_{\mathrm{conv}}=\displaystyle \frac{{Q}^{1/2}{c}_{p}}{4\sqrt{2}}\displaystyle \frac{{k}_{{\rm{B}}}\rho T}{\mu }{\left(\displaystyle \frac{P+{P}_{{\rm{r}}}}{\rho }\right)}^{1/2}{\alpha }^{2}{\left({\rm{\nabla }}-{\rm{\nabla }}^{\prime} \right)}^{3/2},\end{eqnarray} \tag{ 17 }$

assuming that the radiation entropy S_r decreases with height. Here the gradient ${\rm{\nabla }}^{\prime}$ is calculated based on Equations (32) and (42)⁶ in Henyey et al. (1965) given the mean vertical profiles of the envelope. In the adiabatic limit, ${\rm{\nabla }}^{\prime}$ is reduced to ${{\rm{\nabla }}}_{\mathrm{ad}}.$ Then Equation (17) is the same as the formula given by the adiabatic MLT (Cox 1968).

The convective flux calculated in this way with $\alpha =0.55$ is shown as the red line in the top panel of Figure 6, which agrees with the time averaged advection flux in the simulation. Note that the value of α required to explain the flux in the simulations is also quite close to the ratio between the density correlation length and pressure scale height H₀ (see Section 4.1.3), which is the definition of α in MLT. Finally, Figure 6 shows that the difference between the actual gradient and the adiabatic gradient $\tilde{{\rm{\nabla }}}\equiv ({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}})/{{\rm{\nabla }}}_{\mathrm{ad}}$ is quite small, consistent with the expectations of efficient convection.

Taken together, these comparisons demonstrate that the statistical properties of the turbulent energy transport in our radiation hydrodynamic simulation with ${\tau }_{0}\gg {\tau }_{{\rm{c}}}$ are reasonably consistent with the expectations of MLT, although the value of α is somewhat smaller than traditionally assumed in models of massive stars (typical values $\alpha =1.0,...,1.5,$ see, e.g., Yusof et al. 2013; Köhler et al. 2015). This is the first numerical confirmation of the adequacy of MLT in a radiation-dominated regime.

4.2. The Run StarMid when ${\tau }_{0}\gtrsim {\tau }_{c}$

In this regime, rapid radiative diffusion starts to have a non-negligible impact on the properties of convection in the radiation-dominated stellar envelope (the typical thermal timescale near the iron opacity peak is ∼5 dynamical times). Compared with the previous run, we shall see the effects of the reduced optical depth. Note that unlike StarDeep, the photosphere is captured within the domain of this simulation.

4.2.1. Simulation History

The initial condition for StarMid has the iron opacity peak and density inversion at $78{R}_{\odot }.$ The density inversion region is, however, convectively unstable and starts to mix at time $9.4{t}_{0}$ as shown in the left panel of Figure 10. At time $18.8{t}_{0},$ the whole envelope is turbulent (right panel of Figure 10). This temporal evolution is very similar to the run StarDeep, but the properties of the turbulence are different as we now describe.

**Figure 10.** Snapshots of density for the run StarMid at time $9.4{t}_{0}$ (left) and $18.8{t}_{0}$ (right). The box size is in units of ${R}_{\odot }$ while the density unit ${\rho }_{0}$ is given in Table 1. The horizontal box sizes are ${L}_{x}={L}_{y}=10{R}_{\odot }.$
Download figure:
Standard image High-resolution image

**Figure 10.** Snapshots of density for the run StarMid at time $9.4{t}_{0}$ (left) and $18.8{t}_{0}$ (right). The box size is in units of ${R}_{\odot }$ while the density unit ${\rho }_{0}$ is given in Table 1. The horizontal box sizes are ${L}_{x}={L}_{y}=10{R}_{\odot }.$
Download figure:
Standard image High-resolution image

Figure 11 shows the radiation energy density and density weighted vertical velocities ( ${V}_{z,{E}_{{\rm{r}}}},{V}_{z,\rho }$ ) as a function of time. After the initial density inversion is broken around $10{t}_{0},$ $2.4\%$ of the mass is lost through the open top boundary. The averaged vertical velocity also reaches the peak at the same time. The whole envelope oscillates with a period of $4{t}_{0}.$ The vertical velocities can reach 10% of the typical radiation sound speed, which is much larger than the oscillation amplitude in StarDeep (Figure 5). When the envelope expands, it is accelerated by the radiation without the density inversion. But because of the sensitive dependence of the iron opacity on temperature and the temperature decrease with height, the radiation acceleration becomes smaller than g at some height and the acceleration stops. The whole envelope starts to fall back when the vertical velocity reverses.

4.2.2. Time Averaged Vertical Structure and Turbulent Properties

Figure 12 (solid line) shows that there is significant turbulent transport of energy below $65{R}_{\odot },$ but the turbulent transport falls off rapidly at larger heights. At about the same location $65{R}_{\odot },$ the averaged turbulent velocity $\bar{V}$ reaches the isothermal sound speed (middle panel of Figure 12), but is still smaller than the radiation sound speed.

**Figure 12.** Top: time and horizontally averaged vertical profiles of advection flux (black line) and convection flux predicted by mixing length theory (red line). Middle: the averaged turbulent velocity $\bar{V}$ scaled with the radiation sound speed c_s (red line) and isothermal sound speed (black line). Bottom: vertical profile of the superadiabatic parameter $\tilde{{\rm{\nabla }}}$ (black line) and the ratio between the local optical depth per pressure scale height τ and the critical value ${\tau }_{{\rm{c}}}$ (red line). These profiles are for the run StarMid and the time average is done between $17.6{t}_{0}$ and $29.2{t}_{0}.$
Download figure:
Standard image High-resolution image

Figure 13 shows the time and horizontally averaged vertical profiles of density, temperature, opacity, entropy and accelerations for this run. The iron opacity peak moves inward somewhat in radius, but the density inversion in the initial condition is largely gone, replaced with an extended region below $65{R}_{\odot }$ over which the density changes very slowly with height. Gas and radiation are thermally coupled in this region. However, above $65{R}_{\odot },$ the density drops quickly with height and the gas temperature starts to deviate from the radiation temperature. Here the local optical depth per pressure scale height τ becomes smaller than the critical value ${\tau }_{{\rm{c}}}$ and the photon diffusion time becomes smaller than the local dynamic time. This transition to $\tau \lesssim {\tau }_{{\rm{c}}}$ is why the advection flux drops significantly in this region (Figure 12).

**Figure 13.** Time and horizontally averaged vertical profiles of density (top panel of (a)), temperature (bottom panel of (a), solid black line for gas temperature and red line for radiation temperature), entropy (panel (b)), opacity (panel (c)), as well as accelerations (panel (d)) for the run StarMid. As in Figure 7, the dashed lines in panels (a), (b), (c) are initial conditions while other lines are time averaged profiles. The initial condition above $85{R}_{\odot }$ is the region above the photosphere where we set it to be isothermal.
Download figure:
Standard image High-resolution image

The top panel of Figure 12 shows the convective flux calculated based on the time averaged vertical profiles and MLT (Equation (17)) with $\alpha =0.45,$ which is again close to the ratio between the correlation length ${l}_{\rho }$ and scale height H₀. Below $60{R}_{\odot }$ where $\tau \gt {\tau }_{{\rm{c}}},$ the time averaged advection flux is comparable to the MLT convection flux, which is consistent with Figure 6 for the run StarDeep. However, above $65{R}_{\odot }$ the advection flux drops significantly in the simulations, but Equation (17) still predicts a very large convection flux because of the decreasing radiation entropy (see panel (b) of Figure 13). This demonstrates that for $\tau \lesssim {\tau }_{{\rm{c}}},$ convection becomes inefficient and Equation (17) significantly overestimates the amount of advection flux with a constant α.

The vertical profile of $\tau /{\tau }_{{\rm{c}}},$ where $\tau =\rho {\kappa }_{t}{H}_{0}$ is the local optical depth per pressure scale height, shown in the bottom panel of Figure 12 confirms that the transition to inefficient convection happens when $\tau /{\tau }_{{\rm{c}}}\lt 1.$ It also happens roughly when the turbulent velocity $\bar{V}$ becomes larger than the isothermal sound speed (the middle panel of Figure 12), and the superadiabatic parameter $\tilde{{\rm{\nabla }}}$ increases significantly to ∼0.4 (bottom panel of Figure 12). Notice that throughout the whole envelope, $\tilde{{\rm{\nabla }}}$ is comparable to $\bar{V}/{c}_{{\rm{s}}}$ as in Figure 6 for the case StarDeep.

Panel (d) of Figure 13 shows that the volume averaged radiation acceleration a_r is larger than the gravitational acceleration between $50{R}_{\odot }$ and $80{R}_{\odot },$ which is also the region where radiation entropy decreases with height. It is thus perhaps surprising that when convection becomes inefficient above $65{R}_{\odot }$ and there is little turbulent transport of energy, the envelope does not develop a density inversion again. This is because the density weighted radiation acceleration ${\tilde{a}}_{{\rm{r}}}$ (Equation (14), shown as the dashed black line in panel (d) of Figure 13) is significantly smaller than the volume averaged radiation acceleration a_r above $65{R}_{\odot }.$ The sum ${\tilde{a}}_{{\rm{r}}}+{a}_{{\rm{g}}}$ is comparable to g at the iron opacity peak around $60{R}_{\odot }.$ Beyond $65{R}_{\odot },$ it drops below g significantly.

Figure 14 clarifies the origin of the reduced density weighted radiation acceleration. In particular, it shows that density fluctuations reduce the effective radiation acceleration above $65{R}_{\odot }$ because there is an anti-correlation between fluctuations of density and radiation flux. The correlation coefficient ${C}_{\rho ,{F}_{{\rm{r}},z}}$ in Figure 14 shows this explicitly and is always negative. However, the nature of the anti-correlations changes at different heights. Below $65{R}_{\odot },$ when $\tau \gt {\tau }_{{\rm{c}}},$ the situation is similar to the run StarDeep, in that buoyancy causes the anti-correlation between ρ and ${F}_{{\rm{r}},z},$ as ${F}_{{\rm{r}},z}$ is dominated by the advection flux. Above $65{R}_{\odot },$ however, when $\tau \lesssim {\tau }_{{\rm{c}}},$ the situation changes. Convection becomes inefficient and ${F}_{{\rm{r}},z}$ is dominated by the diffusive flux (the advection flux drops as shown in Figure 12). The anti-correlation between ρ and ${F}_{{\rm{r}},z}$ above $65{R}_{\odot }$ is thus not caused by buoyancy. Instead, it is because the "porosity" of the envelope created by the traveling strong shocks allows radiation to propagate through low density regions (Shaviv 1998). This reduces the effective radiation acceleration significantly (see panel (d) of Figure 13). Note also that in StarMid, the standard deviations of density and temperature (scaled with the horizontally averaged quantities) are significantly larger than the corresponding values in the run StarDeep (Figure 9).

**Figure 14.** Top: vertical profiles of density and temperature standard deviations scaled with the horizontally averaged mean values. Bottom: vertical profiles of cross correlation coefficients between $\rho ,T$ and $\rho ,{F}_{{\rm{r}},z}.$ These profiles are for the run StarMid at time $18.8{t}_{0}.$
Download figure:
Standard image High-resolution image

4.3. The Run StarTop when ${\tau }_{0}\ll {\tau }_{c}$

This is the regime when gas and photons are loosely coupled through the whole iron opacity peak region. Convection will be inefficient as in the top region of StarMid. The photosphere is now closer to the region with the iron opacity peak and accurate radiation transfer is critical. For this run, we use a density floor of ${10}^{-10}{\rho }_{0}$ to avoid numerical difficulties.

4.3.1. Simulation History

The iron opacity peak and density inversion are located at $13.7{R}_{\odot }$ initially. This is again convectively unstable and the initial structure is completely destroyed by $50{t}_{0}.$ Just as we found in the simulation StarDeep, this timescale is significantly longer than the buoyancy timescale $1/| {N}_{{\rm{r}}}| ,$ which is ≈t₀ at the iron opacity peak. We therefore again explored the linear growth phase for this run by calculating the density power spectrum at the location of the initial density inversion. The time evolutions of the power ${kf}(\rho )$ in two modes with horizontal wavelengths H₀ and $0.26{H}_{0}$ are shown in Figure 15. As in StarDeep (Figure 4), the long wavelength mode grows exponentially from close to the beginning, with an e-folding time of $6.25{t}_{0}.$ However, in contrast to the StarDeep case, the short wavelength mode is damped initially and only starts to grow quickly after 20t₀ when the e-folding time is $3.33{t}_{0}.$ This later rapid growth is almost certainly due to a nonlinear cascade.

**Figure 15.** Time evolution of the density power spectrum ${kf}(\rho )$ for the run StarTop at height $z=13.7{R}_{\odot }$ during the first $50{t}_{0}.$ The solid black and red lines are for the modes with wavelength H₀ and $0.26{H}_{0}$ , respectively. The long wavelength mode has an e-folding time $6.25{t}_{0}$ as indicated by the dashed black line.
Download figure:
Standard image High-resolution image

We have examined the behavior of other horizontal wavelengths as well, and found that all the modes with horizontal wavelengths larger than $\approx 0.4{H}_{0}$ have similar growth rates after $12{t}_{0},$ while the modes with horizontal wavelengths smaller than $\approx 0.3{H}_{0}$ are damped until the nonlinear cascade sets in. Based on the analysis of Blaes & Socrates (2003), we estimate that modes with wavelengths less than about $2.5{H}_{0}$ are in the regime where radiative diffusion reduces the growth rate of the convective instability. This is larger than the size of the box, so that all convective modes in the simulation are heavily modified by radiative diffusion (consistent with the fact that ${\tau }_{0}\lesssim {\tau }_{{\rm{c}}}$ ). From Equation (59) of Blaes & Socrates (2003), we find that the long wavelength mode in Figure 15 should have a growth rate of $\simeq 0.17/{t}_{0},$ assuming a vertical wavelength that is comparable to the horizontal wavelength (which is also comparable to the vertical width of the density inversion). This agrees very well with the measured growth rate. However, shorter wavelength modes should have smaller growth rates $\propto {k}^{-2},$ in contrast to the behavior measured in the simulations. It is not clear what is causing this discrepancy, although we suspect that the damping of modes with wavelength $\lesssim 0.3{H}_{0}$ is due to the numerical damping rate being larger than the very small growth rate for these modes. Interestingly, traveling acoustic waves should also be unstable according to Equation (62) of Blaes & Socrates (2003), with comparable growth rates to the buoyancy instability (such waves should not be unstable in the StarDeep simulation). Vertically propagating acoustic waves (with infinite horizontal wavelength) should grow fastest, but these would not necessarily show up as growing perturbations at a fixed height. We have not further investigated this possibility, as it is clear that the unstable buoyancy modes dominate the evolution observed in the simulation.

Figure 16 shows snapshots of the density structure of the envelope at time $49.1{t}_{0}$ and $138.9{t}_{0}$ while the histories of the vertical velocities ${V}_{z,{E}_{{\rm{r}}}},{V}_{z,\rho }$ are shown in Figure 17. The whole envelope becomes turbulent after the initial density inversion is destroyed around $50{t}_{0}.$ Unlike the simulations StarDeep and StarMid, the turbulent envelope is not a relatively quiescent convective structure, but instead shows violent, irregular, large amplitude oscillations. The oscillation period varies from ∼5 to $10{t}_{0}.$ The amplitude of ${V}_{z,{E}_{{\rm{r}}}}$ is $0.05{c}_{{\rm{s}},0},$ which is already $\approx 18\%$ of the isothermal sound speed. At time $138.9{t}_{0}$ when ${V}_{z,{E}_{{\rm{r}}}}$ reaches a minimum, the right panel of Figure 16 shows that there is even a low density hole extending over one pressure scale height H₀. The mass loss through the open top boundary is negligible for this run.

**Figure 17.** History of the volume averaged vertical velocity ${V}_{z,{E}_{{\rm{r}}}}$ (solid line at bottom) and ${V}_{z,\rho }$ (dashed line at bottom) for the run StarTop. The volume average is done for the whole simulation box.
Download figure:
Standard image High-resolution image

According to Equation (3), the typical thermal timescale is ${t}_{{\rm{c}}}\sim 0.09{t}_{0}$ for this run. However, for this ${\tau }_{0}\ll {\tau }_{{\rm{c}}}$ regime, t_c is no longer the most relevant timescale for some of the dynamics. Instead, the radiation flux behaves like a roughly constant force reducing the effective gravity (Jiang et al. 2013a). The more relevant timescale is the effective dynamic timescale

$\begin{eqnarray}&&{\tilde{t}}_{0}\equiv \displaystyle \frac{{c}_{{\rm{g}},0}}{| \tilde{g}| },\end{eqnarray} \tag{ 18 }$

where ${c}_{{\rm{g}},0}$ is the isothermal sound speed and $\tilde{g}\equiv g-{\tilde{a}}_{{\rm{r}}}$ is the effective acceleration due to gravity and the radiation force. For an effective acceleration $\tilde{g}=-0.06\;g$ according to the ${\tilde{a}}_{{\rm{r}}}$ shown in panel (d) of Figure 19 (discussed below), the effective dynamic timescale is ${\tilde{t}}_{{\rm{g}}}=4.6{t}_{0},$ which is consistent with the oscillation period of the envelope in Figure 17. Note that the effective scale height is then ${\tilde{H}}_{0}\equiv {c}_{{\rm{g}},0}^{2}/| \tilde{g}| =1.26{H}_{0},$ which happens to be similar to H₀.

4.3.2. Time Averaged Vertical Structures

Figures 18 and 19 show vertical profiles of various time and horizontally averaged quantities in the nonlinear state after $55.6{t}_{0}.$ Unlike the previous two runs StarDeep and StarMid, the averaged advection flux is less than $0.5\%$ of ${F}_{{\rm{r}},i}$ in the whole envelope as shown in the top panel of Figure 18. This is because ${\tau }_{0}\ll {\tau }_{{\rm{c}}}$ everywhere and gas is not able to trap photons. The averaged radiation flux ${F}_{{\rm{r}},0z}$ is almost a constant vertically except very near the top of the domain; the radiation flux is nearly constant in time as well, with $\lesssim 8\%$ variation during the whole simulation.

**Figure 18.** Top: time and horizontally averaged vertical profiles of advection flux (black line), and convection flux according to the non-adiabatic MLT with $\alpha =0.4$ (Equation (17)) as given by Henyey et al. (1965; red line). Middle: the averaged turbulent velocity $\bar{V}$ scaled with the radiation sound speed c_s (red line) and isothermal sound speed (black line). Bottom: vertical profile of superadiabatic parameter $\tilde{{\rm{\nabla }}}.$ The time average is done between $55.6{t}_{0}$ and $171.0{t}_{0}.$ These profiles are for the run StarTop.
Download figure:
Standard image High-resolution image

**Figure 18.** Top: time and horizontally averaged vertical profiles of advection flux (black line), and convection flux according to the non-adiabatic MLT with $\alpha =0.4$ (Equation (17)) as given by Henyey et al. (1965; red line). Middle: the averaged turbulent velocity $\bar{V}$ scaled with the radiation sound speed c_s (red line) and isothermal sound speed (black line). Bottom: vertical profile of superadiabatic parameter $\tilde{{\rm{\nabla }}}.$ The time average is done between $55.6{t}_{0}$ and $171.0{t}_{0}.$ These profiles are for the run StarTop.
Download figure:
Standard image High-resolution image

**Figure 19.** Time and horizontally averaged vertical profiles of density (top panel of (a)), temperature (bottom panel of a, black and red lines for gas and radiation temperature), entropy (panel (b)), opacity (panel (c)), as well as accelerations (panel (d)) for the run StarTop. The dashed lines in (a), (b), and (c) are initial conditions while other lines are time averaged profiles between $55.6{t}_{0}$ and $171.0{t}_{0}.$
Download figure:
Standard image High-resolution image

**Figure 19.** Time and horizontally averaged vertical profiles of density (top panel of (a)), temperature (bottom panel of a, black and red lines for gas and radiation temperature), entropy (panel (b)), opacity (panel (c)), as well as accelerations (panel (d)) for the run StarTop. The dashed lines in (a), (b), and (c) are initial conditions while other lines are time averaged profiles between $55.6{t}_{0}$ and $171.0{t}_{0}.$
Download figure:
Standard image High-resolution image

Because the radiation flux dominates over the advection flux, the full super-Eddington flux contributes to the radiation acceleration of the gas. However, because of the large density fluctuations, most photons go through regions that have below average density. This reduces the density weighted radiation acceleration ${\tilde{a}}_{{\rm{r}}}$ (Equation 14), as shown in panel (d) of Figure 19. However, unlike the assumptions made in previous models (Shaviv 1998), the "porosity" of the envelope due to convection cannot reduce the effective radiation acceleration to a value below g. At $z\approx 13.3{R}_{\odot }$ where the iron opacity peak is located now, ${\tilde{a}}_{{\rm{r}}}$ is still about 6% larger than g, while a_r is about 10% larger than g. This is why the average density profile still shows a mild inversion near the iron opacity peak, as shown in panel (a) of Figure 19. We stress that although the density inversion remains, the resulting structure is nonetheless very different from the hydrostatic density inversion in the initial condition. The high density region near the iron opacity peak is continuously reformed and then mixed, triggering violent oscillations, shocks, and large density fluctuations in the stellar envelope. Physically, in the regime ${\tau }_{0}\ll {\tau }_{{\rm{c}}},$ photon diffusion is so rapid that radiation pressure does not respond to the density and velocity fluctuations. As the turbulent velocity becomes larger than the isothermal sound speed, shocks are formed in the envelope (Figure 16), which causes the large density fluctuations. This phenomenon is also observed in radiation magneto-hydrodynamic simulations of black hole accretion disks (Turner et al. 2003; Jiang et al. 2013c).

Figure 19 shows that the radiation and gas entropies S_r and S_g decrease rapidly with height around the iron opacity peak. If we still adopt ${l}_{\rho }/{H}_{0}\approx 0.4$ for the mixing length parameter α, Equation (17) over-predicts the flux relative to the simulations by a factor of more than ∼3, and the location of the peak in the predicted flux also offsets from the peak of time averaged advection flux given by the simulations, highlighting that existing models of inefficient convection in the radiation pressure-dominated regime are not accurate compared to our simulation results. The failure of convection to carry a significant energy flux is intimately connected with the turbulent velocities becoming supersonic, at least relative to the gas (isothermal) sound speed. This is expected on theoretical grounds (Section 2) and is also consistent with the simulation results, as shown in the middle panel of Figure 18. The superadiabatic parameter $\tilde{{\rm{\nabla }}}$ is also comparable to $\tilde{V}/{c}_{{\rm{s}}}$ and similar to the top region in StarMid (Figure 12).

4.3.3. Properties of the Turbulent Structures

Figure 20 quantifies the large density fluctuations and strong anti-correlation between density and radiation flux in simulation StarTop. The correlation coefficient ${C}_{\rho ,{F}_{{\rm{r}},z}}$ is almost $-1$ in the turbulent region between $12.8{R}_{\odot }$ and $13.6{R}_{\odot },$ while ${C}_{\rho ,T}$ only fluctuates around 0. This demonstrates that unlike for StarDeep, the anti-correlation between density and radiation flux is not due to buoyancy in the radiation pressure-dominated regime but is instead due to the photons preferentially diffusing through low density channels. The scaled standard deviations $\delta \rho /\rho$ and $\delta T/T$ are also much larger than the corresponding values in StarDeep and StarMid, which is consistent with the images in Figure 16.

The typical size of the turbulent structures can also be quantified by the correlation length ${l}_{\rho },$ which varies between $0.4{H}_{0}$ and $0.5{H}_{0}.$ This is likely set primarily by geometry as in standard mixing length arguments. In particular, the optical depth across the correlation length ${\tau }_{l}$ is much smaller than ${\tau }_{{\rm{c}}},$ which means that photon diffusion is very rapid for all of the turbulent fluctuations and cannot be important for setting a characteristic scale. The correlation length ${l}_{\rho }$ is also much smaller than the characteristic turnover wavelength of acoustic waves driven unstable by the dependence of opacity on density (Blaes & Socrates 2003). Once again, it is simply the fact that the convective turbulence is supersonic with respect to the isothermal sound speed in the gas, and that radiation pressure cannot respond due to the rapid diffusion, that is producing these nonlinear structures.

5. DISCUSSIONS AND CONCLUSIONS

We have presented three-dimensional (3D) radiation hydrodynamic simulations of the radiation pressure-dominated envelopes of massive stars near the iron opacity peak. The resulting dynamics depends sensitively on the ratio of the photon diffusion time to the dynamical time. This is set by the ratio of the local optical depth ${\tau }_{0}$ to the critical optical depth ${\tau }_{{\rm{c}}}=c/{c}_{{\rm{g}}},$ where c is the speed of light and c_g is the isothermal sound speed. Our simulations cover the regimes ${\tau }_{0}\gg {\tau }_{{\rm{c}}}$ (StarDeep), ${\tau }_{0}\sim {\tau }_{{\rm{c}}}$ (StarMid), and ${\tau }_{0}\ll {\tau }_{{\rm{c}}}$ (StarTop).

The iron opacity peak at $T\approx 1.6\times {10}^{5}$ K leads to locally super-Eddington fluxes in the envelopes of massive stars. In radiative equilibrium this produces a pronounced density inversion in one-dimensional stellar models (Joss et al. 1973; Paxton et al. 2013). The resulting density inversions are convectively unstable, but the properties of the resulting convection zones have until now been very uncertain when ${\tau }_{0}\lesssim {\tau }_{{\rm{c}}}$ and convection is expected to be relatively inefficient.

Figures 21 and 22 compare the time and horizontally averaged turbulent properties of our three runs, showing the gas, radiation, and kinetic energy densities (Figure 21) and the ratio of the convective and diffusive energy fluxes, respectively.

**Figure 21.** Time and horizontally averaged vertical profiles of radiation energy density (E_r, red lines), gas internal energy density (E_g, solid black lines), and kinetic energy density (E_k, green lines). From top to bottom, they are for StarDeep, StarMid and StarTop respectively. The black dashed vertical lines in each panel indicate the time averaged location of the iron opacity peak. The horizontal axis is the distance from the bottom boundary z₀ in units of the scale height H₀ for each run.
Download figure:
Standard image High-resolution image

When ${\tau }_{0}\gg {\tau }_{{\rm{c}}}$ as in StarDeep, convection is efficient and the averaged radiation advection flux is consistent with the convection flux as calculated using the MLT formalism (Equation (17)). We find that the mixing length parameter required to match the advection flux in the 3D simulation is $\alpha =0.55,$ roughly consistent with the typical size of the turbulent structures ${l}_{\rho }/{H}_{0}\approx 0.3-0.4$ (as assumed in MLT). The density inversion, present in the 1D radiative equilibrium initial condition, is absent in the turbulent state. This is because a significant fraction of the luminosity is transported by convection (Figure 22), which reduces the radiation acceleration to be sub-Eddington. For StarDeep, the turbulent velocity is very subsonic with the turbulent kinetic energy less than 1% of the thermal energy (Figure 21). The difference between the actual thermal gradient ∇ and the adiabatic gradient ${{\rm{\nabla }}}_{\mathrm{ad}}$ is smaller than 8% of ${{\rm{\nabla }}}_{\mathrm{ad}}.$ All of these properties are consistent with expectations for efficient convection.

**Figure 22.** Vertical profiles of the ratio between time averaged advection (convection) flux ${F}_{\mathrm{adv}}={v}_{z}{E}_{{\rm{r}}}$ and diffusive flux ${F}_{{\rm{r}},0z}$ for the three runs. Convection transports significant energy in StarDeep and in the deeper layers of StarMid. But in the outer layers of StarMid and throughout StarTop, photon diffusion is so rapid that it strongly suppresses the turbulent transport of energy.
Download figure:
Standard image High-resolution image

When ${\tau }_{0}\lt {\tau }_{{\rm{c}}},$ the advection flux drops significantly as photons cannot be trapped by the fluid due to rapid diffusion. We see this transition in the surface layers of run StarMid. As shown in Figure 12, below $60{R}_{\odot },$ the averaged radiation advection flux is consistent with MLT (Equation (17)). The turbulent velocity is subsonic and $\tilde{{\rm{\nabla }}}\ll 1$ as in StarDeep. However above $60{R}_{\odot },$ Equation (17) significantly overpredicts the average convection flux if we use a constant α, a choice motivated by the fact that the correlation length of the turbulence ( ${l}_{\rho }/{H}_{0}$ ) is almost the same through the whole envelope. Above $60{R}_{\odot }$ the turbulent velocity reaches the isothermal sound speed and $\tilde{{\rm{\nabla }}}$ increases significantly.

The run StarTop has ${\tau }_{0}\ll {\tau }_{{\rm{c}}}$ and convection is very inefficient. In this case the average radiation advection flux is $\ll 1\%$ of the diffusive flux throughout the envelope (Figure 22). If we adopt the non-adiabatic MLT model (e.g., Henyey et al. 1965; Magic et al. 2015) with $\alpha =0.4,$ which is the ratio of the turbulent correlation length to the pressure scale height, we find that it still significantly over-predicts the convection flux (Figure 18). The difference is probably caused by a combination of the porosity effect and vertical oscillation of the envelope as explained blow.

The turbulent velocity in StarTop is larger than the isothermal sound speed and strong shocks are formed, driving large density fluctuations. The bottom panel of Figure 21 shows, however, that the kinetic energy density only reaches the gas internal energy density, which is still much smaller than the radiation energy density. The large density fluctuations allow radiation to escape somewhat more efficiently (Shaviv 1998) than through a homogenous medium ("porosity"), reducing the effective radiation acceleration. This effect is not, however, large enough to make the radiation acceleration smaller than the gravitational acceleration. The density inversion present in the initial conditions is initially destroyed by convection and the envelope expands, but the gas quickly falls back because of the reduced opacity and radiation acceleration. The density inversion then reforms. This cycle repeats, with the envelope undergoing strong vertical oscillations with a period of a few hours. The density inversion predicted by 1D radiative equilibrium models still exists in the time averaged vertical profile of StarTop, in contrast to StarDeep and StarMid where convection is able to largely eliminate the density inversion.

5.1. Observational Consequences

A signature of inefficient convection is the appearance of supersonic turbulent velocities with respect to the isothermal sound speed (but not the radiation sound speed). This is because rapid diffusion leads the photons and gas to decouple and the radiation pressure cannot respond to density fluctuations. This in turn makes the gas highly compressible, leading to shocks and large density fluctuations.

In the calculations that include the photosphere (StarTop and StarMid) we observe turbulent velocities reaching the isothermal sound speed close to the stellar surface ( $\approx 50\;\mathrm{km}\;{{\rm{s}}}^{-1}$ ). These large velocities can impact the spectroscopic measurements of line profiles, as quantified by the macroturbulence (microturbulence) parameters. These parameters quantify the velocity fields at the stellar photosphere with correlation length larger (smaller) than the line-forming region. The extent of the line forming region is comparable to the local pressure scale-height, which is about the correlation length of the turbulent structures in our calculations. Therefore we expect the velocity fields observed in the simulations to potentially impact both micro- and macro- turbulence measurements. This is important, as the broadening of spectral lines is used to measure the (projected) rotational velocity of stars, which in the presence of extra turbulent broadening could be be systematically higher than the actual rotation rate (Ramírez-Agudelo et al. 2015). The idea that inefficient convection at the iron opacity peak could be responsible for the extra broadening of spectral lines in hot massive stars has been discussed previously (Cantiello et al. 2009), where they focused on the role of gravity waves excited by turbulent convection. Recently the potentially observable signature of turbulent pressure in the inefficient outer convective regions of massive stars has been highlighted by Grassitelli et al. (2015), who showed a correlation between observed macroturbulence and turbulent pressure (as calculated in the context of 1D MLT).

The vertical oscillations of the envelope in StarMid and StarTop also cause the radiation flux coming from the photosphere to vary with time. The standard deviations of the total radiation flux leaving from the photosphere are $\approx 13\%$ and $\approx 2.7\%$ of the time averaged mean values in StarMid and StarTop respectively. This means that this oscillation can be observed in principle. However, direct comparisons with observations require global calculations of the envelope.

Strange mode instabilities are commonly associated with opacity bumps and density inversions in 1D models of massive stars (e.g., Kiriakidis et al. 1993; Saio 2009; Saio et al. 2013). These standing wave instabilities are fundamentally acoustic in nature, and are generally associated with a trapping region around the opacity bump/density inversion. The excitation mechanism of these modes is due to opacity variations in the acoustic wave, and is related to the growth mechanism of traveling acoustic waves discussed by Blaes & Socrates (2003). As we briefly noted above in Section 4.3.1, conditions in StarTop at least are such that these traveling wave instabilities should be driven, and the fact that the density inversion persists in the time-averaged structure may permit the existence of trapped strange modes. It would be interesting to do a global stability analysis on the time-averaged structure of StarMid and StarDeep, to see if standing acoustic wave instabilities could still be trapped even in the absence of density inversions.

5.2. Implications for 1D Stellar Evolution

For ${\tau }_{0}\gtrsim {\tau }_{{\rm{c}}},$ we find that MLT provides a reasonable description of the convective energy transport in the 3D simulations. However, the mixing length parameters we find ( $\alpha \approx 0.5$ ) are noticeably smaller than the canonically adopted values ( $\alpha \gt 1,$ see, e.g., Yusof et al. 2013; Köhler et al. 2015). We also find that as ${\tau }_{0}/{\tau }_{{\rm{c}}}$ decreases, the mixing length parameter decreases. This is exactly the opposite of what has been proposed by different groups dealing with the numerical difficulties associated with radiation-dominated stellar envelopes in 1D. They suggested that convection in this situation should be artificially enhanced (equivalent to increasing the α parameter with respect to the canonical values above), justifying this choice with either the need to remove an "unphysical density inversion" (Maeder 1987; Ekström et al. 2012) or the need to include some (unknown) extra energy transport mechanism (Paxton et al. 2013).

Moreover we find that density inversions are present in the time averaged structure of the stellar envelope when ${\tau }_{0}\ll {\tau }_{{\rm{c}}},$ i.e., when an opacity peak is close to the stellar surface. In the case of the iron opacity peak, this occurs for very massive main sequence stars, though we expect the same physics to occur in cooler objects due to the effect of the H and He opacity bumps (not studied in this work). The behavior of these density inversions is quite different as compared to the steady 1D solution, as they tend to be cyclically destroyed and reformed and are associated with very large density inhomogeneities. One important effect that is currently missing in the 1D models is the inclusion of a "porosity factor," reducing the effective radiation acceleration when ${\tau }_{0}\lt {\tau }_{{\rm{c}}}$ (e.g., Paczyński 1969). The dependency of the porosity on optical depth, resolution and magnetic fields will be the focus of future work.

The turbulent, non-stationary behavior of radiation-dominated envelopes suggests a possible modification of the stellar mass loss relative to spherically symmetric hydrostatic models (e.g., Owocki et al. 2004). Unfortunately, as we only simulate a local patch of the envelope, we cannot determine the fate of the mass that leaves our simulation box. The mass loss rate also cannot be calculated self-consistently based on these simulations. This requires global calculations of the envelope to capture the change of gravitational acceleration and geometric dilution correctly.

The majority of massive stars are found in binary systems (Sana et al. 2012). Depending on their orbital separation and the evolution of their radii, these stars can fill their Roche lobe and strongly interact already during the main sequence (e.g., de Mink et al. 2014). As shown in this paper, the energy transport in the radiation-dominated envelopes of these stars is not properly modeled in current 1D calculations, resulting in very uncertain values for the stellar radii. This could have important implications for the evolution of massive binaries and, for example, the predicted rates of BH–BH binary mergers, which are promising sources of gravitational waves (Belczynski et al. 2014).

One limitation of the simulations presented here is that we have neglected the effects of the stellar magnetic field. Recently, magnetic buoyancy has been found to significantly enhance the radiation advection flux, even with a relatively weak magnetic field (Blaes et al. 2011; Jiang et al. 2014a). A fraction of massive stars do show significant magnetization (Wade et al. 2014), and dynamo generated magnetic fields could be present in the iron convection zone (Cantiello & Braithwaite 2011). Therefore the effects of magnetic field should be studied in future work.

Y.F.J thanks Jeremy Goodman, Jim Stone, Shane Davis and all members of the SPIDER network for helpful discussions. E.Q. thanks Mark Krumholz and Todd Thompson for useful conversations. This work was supported by the computational resources provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center; the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575; and the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Y.F.J. is supported by NASA through Einstein Postdoctoral Fellowship grant number PF-140109 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. E.Q. is supported in part by a Simons Investigator Award from the Simons Foundation. This project was supported by NASA under TCAN grant number NNX14AB53G and the NSF under grants PHY 11-25915 and AST 11-09174.

LOCAL RADIATION HYDRODYNAMIC SIMULATIONS OF MASSIVE STAR ENVELOPES AT THE IRON OPACITY PEAK

Article metrics

Permissions

Author affiliations

Author notes

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION

2. THE MODEL PARAMETERS

3. NUMERICAL METHOD

3.1. Equations

3.2. The Opacity

3.3. Initial and Boundary Conditions

4. RESULTS

4.1. The Run StarDeep when ${\tau }_{0}\gg {\tau }_{c}$

4.1.1. Simulation History

4.1.2. Time Averaged Vertical Structures

4.1.3. Time Averaged Horizontal Structures

4.1.4. Comparison with MLT

4.2. The Run StarMid when ${\tau }_{0}\gtrsim {\tau }_{c}$

4.2.1. Simulation History

4.2.2. Time Averaged Vertical Structure and Turbulent Properties

4.3. The Run StarTop when ${\tau }_{0}\ll {\tau }_{c}$

4.3.1. Simulation History

4.3.2. Time Averaged Vertical Structures

4.3.3. Properties of the Turbulent Structures

5. DISCUSSIONS AND CONCLUSIONS

5.1. Observational Consequences

5.2. Implications for 1D Stellar Evolution

Footnotes

LOCAL RADIATION HYDRODYNAMIC SIMULATIONS OF MASSIVE STAR ENVELOPES AT THE IRON OPACITY PEAK

Article metrics

Permissions

Share this article

Author affiliations

Author notes

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION

2. THE MODEL PARAMETERS

3. NUMERICAL METHOD

3.1. Equations

3.2. The Opacity

3.3. Initial and Boundary Conditions

4. RESULTS

4.1. The Run StarDeep when {\tau }_{0}\gg {\tau }_{c}

4.1.1. Simulation History

4.1.2. Time Averaged Vertical Structures

4.1.3. Time Averaged Horizontal Structures

4.1.4. Comparison with MLT

4.2. The Run StarMid when {\tau }_{0}\gtrsim {\tau }_{c}

4.2.1. Simulation History

4.2.2. Time Averaged Vertical Structure and Turbulent Properties

4.3. The Run StarTop when {\tau }_{0}\ll {\tau }_{c}

4.3.1. Simulation History

4.3.2. Time Averaged Vertical Structures

4.3.3. Properties of the Turbulent Structures

5. DISCUSSIONS AND CONCLUSIONS

5.1. Observational Consequences

5.2. Implications for 1D Stellar Evolution

Footnotes

4.1. The Run StarDeep when ${\tau }_{0}\gg {\tau }_{c}$

4.2. The Run StarMid when ${\tau }_{0}\gtrsim {\tau }_{c}$

4.3. The Run StarTop when ${\tau }_{0}\ll {\tau }_{c}$