This paper was written collaboratively, on the web, using Authorea. JF acknowledges partial support from NSF under grant no. AST-1205732 and through a Lee DuBridge Fellowship at Caltech. DL is supported by a Hertz Foundation Fellowship and the National Science Foundation Graduate Research Fellowship under Grant No. DGE 1106400. EQ was supported in part by a Simons Investigator award from the Simons Foundation and the David and Lucile Packard Foundation. This research was supported by the National Science Foundation under grant No. NSF PHY11- 25915 and by NASA under TCAN grant No. NNX14AB53G.
\label{model}
Our stellar models are constructed using the MESA stellar evolution code \citep{paxton:11,paxton:13}, version 6794. The models are run using the following inlist controls file.
\& star_job
kappa_file_prefix = 'gs98'
/ ! end of star_job namelist
\& controls
initial_mass = 12.
initial_z = 0.020
sig_min_factor_for_high_Tcenter = 0.01
Tcenter_min_for_sig_min_factor_full_on = 3.2d9
Tcenter_max_for_sig_min_factor_full_off = 2.8d9
logT_max_for_standard_mesh_delta_coeff = 9.0
logT_min_for_highT_mesh_delta_coeff = 10
delta_Ye_highT_limit = 1d-3
okay_to_reduce_gradT_excess = .true.
allow_thermohaline_mixing = .true.
thermo_haline_coeff = 2.0
overshoot_f_above_nonburn = 0.035
overshoot_f_below_nonburn = 0.01
overshoot_f_above_burn_h = 0.035
overshoot_f_below_burn_h = 0.0035
overshoot_f_above_burn_he = 0.035
overshoot_f_below_burn_he = 0.0035
overshoot_f_above_burn_z = 0.035
overshoot_f_below_burn_z = 0.0035
RGB_wind_scheme = 'Dutch'
AGB_wind_scheme = 'Dutch'
RGB_to_AGB_wind_switch = 1d-4
Dutch_wind_eta = 0.8
include_dmu_dt_in_eps_grav = .true.
use_Type2_opacities = .true.
newton_itermin = 2
mixing_length_alpha = 1.5
MLT_option = 'Henyey'
allow_semiconvective_mixing = .true.
alpha_semiconvection = 0.01
mesh_delta_coeff = 1.
varcontrol_target = 5d-4
max_allowed_nz = 10000
mesh_dlog_pp_dlogP_extra = 0.4
mesh_dlog_cno_dlogP_extra = 0.4
mesh_dlog_burn_n_dlogP_extra = 0.4
mesh_dlog_3alf_dlogP_extra = 0.4
mesh_dlog_burn_c_dlogP_extra = 0.10
mesh_dlog_cc_dlogP_extra = 0.10
mesh_dlog_co_dlogP_extra = 0.10
mesh_dlog_oo_dlogP_extra = 0.10
velocity_logT_lower_bound=9
dX_nuc_drop_limit=5d-3
dX_nuc_drop_limit_at_high_T = 5d-3 !
screening_mode = 'extended'
max_iter_for_resid_tol1 = 3
tol_residual_norm1 = 1d-5
tol_max_residual1 = 1d-2
max_iter_for_resid_tol2 = 12
tol_residual_norm2 = 1d99
tol_max_residual2 = 1d99
min_timestep_limit = 1d-12 ! (seconds)
delta_lgL_He_limit = 0.1 !
dX_nuc_drop_max_A_limit = 52
dX_nuc_drop_min_X_limit = 1d-4
dX_nuc_drop_hard_limit = 1d99
delta_lgTeff_limit = 0.5
delta_lgL_limit = 0.5
delta_lgRho_cntr_limit = 0.02
T_mix_limit = 0
/ ! end of controls namelist
The most important feature of this model is that it contains significant convective overshoot, especially above convective zones. It is non-rotating, thus there is no rotational mixing.
Just before core O-burning, we change to a 201-isotope reaction network:
change_net = .true.
new_net_name = 'mesa_201.net'
Although our choices affect details of the model (e.g., He core mass), the general features of our model are robust. It always explodes as a red supergiant. It always undergoes convective core C-burning, followed by shell C-burning, core O/Ne-burning, shell O-burning, core Si-burning, shell-Si burning, and then CC. The approximate convective properties (as described by MLT) are not strongly affected by model parameters. Since these properties are most important for IGW AM transport, we argue that the general features of IGWs described in this work are fairly robust against uncertain parameters in our massive star models.
\label{wavestar}
Here we describe our methods for estimating the wave frequency \(\omega_*(r)\) of waves that dominate the AM flux \(J_*(r)\) as IGW propagate through a star and are attenuated. For a more detailed discussion, see F14, whose methods are very similar to ours.
Consider a train of IGWs generated via shell convection which propagate inwards through underlying stable stratification. Upon generation, the IGW carry an energy flux \(\dot{E}_0\) and AM flux \(\dot{J}_0\). We assume the IGW have a frequency spectrum which is initially peaked around \(\omega_c\), with a power law fall off at higher frequencies such that \[\label{eqn:spectrum} \frac{d \dot{E}_0(\omega)}{d\omega} \sim \frac{\dot{E}_0}{\omega_c} \left(\frac{\omega}{\omega_c}\right)^{-a},\] where \(a\) is the slope of the frequency spectrum, which is somewhat uncertain. As in F14, we expect a spectrum slope in the range \(3 \lesssim a \lesssim 7\), and we use a value of \(a=4.5\) in our calculations. Lower frequency waves have shorter wavelengths and slower group velocities, making them more prone to both radiative and non-linear damping. Thus, as waves propagate inwards, low frequency waves may damp out, and increasing the wave frequency \(\omega_*\) of waves that dominate AM transport.
F14 show that radiative damping leads to \[\label{eqn:omstar2} \omega_*(r) = {\rm max} \bigg[ \omega_c \ , \ \bigg(\frac{4}{a} \int^{r_c}_r dr \frac{\lambda^{3/2} N_T^2 N K }{r^3} \bigg)^{1/4} \bigg] .\] and the corresponding AM flux \[\label{eqn:Jstar} \dot{J}_*(r) \sim \bigg[\frac{\omega_*(r)}{\omega_c}\bigg]^{-a} \dot{J}_0.\] In equation \ref{eqn:omstar2}, \(r_c\) is the radius of the inner edge of the convective zone, \(\lambda = l(l+1)\), \(l\) is the angular index of the wave (which corresponds to its spherical harmonic dependence, \(Y_{lm}\)), \(N_T\) is the thermal part of the Brunt-Väisälä frequency, and \(K\) is the thermal diffusivity. In what follows, we focus on \(l=1\) waves because they have the longest damping lengths and may dominate the AM flux when the waves are heavily damped. Moreover, focusing on \(l=1\) waves allows us to estimate maximum spin frequencies, although slower spin frequencies can be obtained when higher values of \(l\) and \(m\) contribute to the AM flux.
IGW will overturn and break, leading to local energy/AM deposition, if they obtain sufficiently non-linear amplitudes. Here we estimate those amplitudes and the AM flux that can be carried toward the center of the star as waves are non-linearly attenuated. For traveling waves in the WKB limit, it is well known that the radial wave number is \[\label{eqn:kr} k_r = \frac{\sqrt{\lambda} N}{r \omega} \, .\] It is straightforward to show that the radial displacement \(\xi_r\) associated with IGW of frequency \(\omega\) carrying an energy flux \(\dot{E}\) is \[\label{eqn:eflux} | \xi_r | = \bigg[ \frac{\sqrt{\lambda} \dot{E} }{\rho N r^3 \omega^2} \bigg]^{1/2} \, .\] The waves become non-linear and break when \[\label{eqn:nl} | k_r \xi_r | = \bigg[ \frac{\lambda^{3/2} N \dot{E} }{\rho r^5 \omega^4} \bigg]^{1/2} \sim 1 \, .\] In the absence of damping, \(\dot{E}\) is a conserved quantity. Therefore, waves become more non-linear as they propagate into regions with larger \(N\), lower density, or smaller radius. In our problem, the geometrical focusing (i.e., the \(r\)-depdendence) is the most important feature of equation \ref{eqn:nl}, and causes waves to non-linearly break as they propagate inward. Note also the \(\omega^{-2}\) dependence of equation \ref{eqn:nl}, which causes low frequency waves to preferentially damp.
Equation \ref{eqn:nl} entails there is a maximum energy flux that can be carried by waves of frequency \(\omega\), \[\label{eqn:emax} \dot{E}_{\rm max} = \frac{ A^2 \rho r^5 \omega^4}{\lambda^{3/2} N} \, ,\] for waves that non-linearly break when \(|k_r \xi_r| = A \sim 1\). When the waves are highly non-linear, the waves which dominate the energy flux are those which are on the verge on breaking. To determine their frequency, we use the frequency spectrum \ref{eqn:spectrum} to find \[\label{eqn:efreq} \dot{E}_0 \bigg( \frac{\omega}{\omega_c} \bigg)^{1-a} \sim \dot{E}_{\rm max}.\] Solving equation \ref{eqn:efreq} yields the wave frequency which dominates energy transport, \[\label{eqn:omstarnl} \omega_* \sim {\rm max} \bigg[ \omega_c \ , \ \omega_c \bigg( \frac{ A^2 \rho r^5 \omega_c^4}{\lambda^{3/2} N \dot{E}_0} \bigg)^{-1/(a+3)} \bigg] \, .\] We expect frequency spectra with slopes somewhere near \(3 \lesssim a \lesssim 7\). Therefore the exponent in equation \ref{eqn:omstarnl} is quite small, and in most cases, \(\omega_*\) does not increase to values much larger than \(\omega_c\).
Substituting equation \ref{eqn:omstarnl} back into equation \ref{eqn:emax} allows us to solve for the energy and AM flux as a function of radius due to non-linear attenuation. The result is \[\label{eqn:jstarnl} \dot{J}_* \sim \bigg[ \frac{\omega_*(r)}{\omega_c} \bigg]^{-a} \dot{J}_0 \sim {\rm min}\bigg[ \dot{J}_0 \ , \ \bigg(\frac{A^2 \rho r^5 \omega_c^4}{\lambda^{3/2} N \dot{E}_0} \bigg)^{a/(a+3)} \dot{J}_0 \bigg] \, .\]
During C-shell burning, we find that radiative diffusion damps waves near the shell burning convective zone, while non-linear breaking damps waves near the center of the star. In this case, we first calculate \(\omega_*\) and its corresponding energy flux \(\dot{E}_*\) via equations \ref{eqn:omstar2} and \ref{eqn:spectrum}. We then substitute the value of \(\dot{E}_*\) for \(\dot{E}_0\) in equation \ref{eqn:omstarnl}. The appropriate value of \(\omega_*\) is then \(\omega_* = {\rm max} \big[ {\rm Eqn.}\) \ref{eqn:omstar2}\(, {\rm Eqn.}\)\ref{eqn:omstarnl}\( \big]\). The corresponding AM flux is \(\dot{J}_* \sim \Big[ \frac{\omega_*(r)}{\omega_c} \Big]^{-a} \dot{J}_0\).
The cores of massive stars nearing death cool primarily through neutrino emission, so it is not unreasonable to think that waves may be damped via neutrino emission. We calculate neutrino energy loss rates in the same manner as \cite{murphy:04}. We find that neutrino damping time scales are always longer than the wave crossing timescale \[\label{eqn:tcross} t_{\rm cross} = \int^{r_c}_0 \frac{dr}{v_g}\, ,\] where the IGW radial group velocity is \(v_g = r \omega^2/(\sqrt{\lambda} N)\). This is not surprising, as \cite{murphy:04} found neutrino growth/damping rates were slower than stellar evolution time scales. Neutrino damping can therefore be safely ignored.