A Transparent Window into Early-type Stellar Variability

Subsurface convection zones (CZs) are ubiquitous in early-type stars (Cantiello et al. 2009; Cantiello & Braithwaite 2019). Driven by narrow opacity peaks (Stothers & Chin 1993), these thin convective regions transport little heat but play an important role in setting the magnetic fields (Cantiello & Braithwaite 2011; Jermyn & Cantiello 2020) and chemical mixing (Michaud 1970).

Subsurface CZs may also play a role in setting the photometric variability of early-type stars, though the precise connection between these is not well-understood. In particular, one of the most exciting potential observational probes of these subsurface CZs is stochastic low-frequency variability (SLF; Blomme et al. 2011; Bowman et al. 2019), but there is active debate into the origin of the SLF phenomenon (Lecoanet et al. 2019). Subsurface convective motions form one possible explanation (Cantiello et al. 2021), as do internal gravity waves excited by core convection (Bowman et al. 2019) and wind mass loss (Krtička and Feldmeier 2021).

Here we demonstrate that these CZs are not present in as wide a range of stars as previously believed. In particular, there are regions which 1D stellar evolution models report to be convectively unstable but which fall below the critical Rayleigh number for onset of convection. This means that the radiative conductivity is large enough to turn off the convective instability in these regions.

We begin in Section 2 with an overview of the theory of stellar convection and the critical Rayleigh number. In Section 3 we compute stellar evolutionary tracks and use these to demonstrate that over certain mass and age windows some subsurface CZs are subcritical. Importantly, this leaves a metallicity-dependent stability window where there are no subsurface CZs. This region opens at around 8M_⊙ and, for the metallicity of the Small Magellanic Cloud (SMC, Z = 0.002), extends up to 35M_⊙. The stability window is then an excellent target for understanding how much of the SLF phenomenon is due to subsurface convection, as it should be possible to compare the SLF properties of stars inside and outside the window to see if there is any abrupt change, or if the properties better track the more gradual evolution in core properties. We conclude with a discussion of this and other observational implications in Section 4.

In 1D stellar evolution calculations, convection is assumed to be active whenever the Ledoux condition for instability holds, namely

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{rad}}\gt {{\rm{\nabla }}}_{\mathrm{ad}}+{{\rm{\nabla }}}_{\mu }{\left.\displaystyle \frac{\partial \mathrm{ln}T}{\partial \mathrm{ln}\mu }\right|}_{P,\rho }.\end{eqnarray} \tag{ 1 }$$

Here μ is the mean molecular weight,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\equiv \displaystyle \frac{d\mathrm{ln}\mu }{d\mathrm{ln}P}\end{eqnarray} \tag{ 2 }$$

is the composition gradient,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{ad}}\equiv {\left.\displaystyle \frac{\partial \mathrm{ln}T}{\partial \mathrm{ln}P}\right|}_{s}\end{eqnarray} \tag{ 3 }$$

is the adiabatic temperature gradient,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{rad}}\equiv \displaystyle \frac{3\kappa {LP}}{64\pi {GM}\sigma {T}^{4}}\end{eqnarray} \tag{ 4 }$$

is the radiative temperature gradient, κ is the opacity, ρ is the density, L is the luminosity, P is the pressure, T is the temperature, G is the gravitational constant, M is the mass below the point of interest, and σ is the Stefan–Boltzmann constant.

However, there is an additional requirement for convection, namely that the Rayleigh number

$$\begin{eqnarray}&&\mathrm{Ra}\equiv \displaystyle \frac{g({{\rm{\nabla }}}_{\mathrm{rad}}-{{\rm{\nabla }}}_{\mathrm{ad}})\delta {r}^{3}}{\nu \alpha }\left(\displaystyle \frac{\delta r}{h}\right)\left(\displaystyle \frac{4-3\beta }{\beta }\right)\end{eqnarray} \tag{ 5 }$$

exceed the critical value Ra_crit (e.g., Chandrasekhar 1961; Shore 1992). The Rayleigh number expresses the importance of buoyancy compared to diffusive processes; when Ra < Ra_crit, buoyant forces are not strong enough to overcome diffusive processes, and there is no convection. Here

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{16\sigma {T}^{3}}{3\kappa {c}_{p}{\rho }^{2}}\end{eqnarray} \tag{ 6 }$$

is the thermal diffusivity, c_p is the specific heat at constant pressure, g is the acceleration of gravity, δ r is the thickness of the Ledoux-unstable layer, ν is the kinematic viscosity, $h\equiv -{dr}/d\mathrm{ln}P$ is the pressure scale height, and β = P_gas/P is the radiation parameter. The factor depending on β arises from the thermal expansion coefficient, and is equal to the ratio of the density susceptibility to the temperature susceptibility (χ_ρ /χ_T ). See Appendix B for details on how we compute ν.

The critical Rayleigh number depends on the precise setup of the CZ, including geometry, the shape of the opacity profile, rotation, and the diffusivities (Net et al. 2008; Calkins et al. 2015; Garcia et al. 2018). For simplicity, we adopt Chandrasekhar's theory (Chandrasekhar 1961, Ch. 29), and assume that convection occurs by either direct onset or overstability, whichever has a lower Ra_crit. For non-rotating systems this produces Ra_crit ∼ 10³, while for rotating systems it is somewhat greater because rotation stabilizes against convection. See Appendix D for details of how we compute Ra_crit.

We calculated stellar evolutionary tracks for stars ranging from 1.1 to 60M_⊙ using revision 15140 of the Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019) software instrument. Details on the MESA microphysics inputs are provided in Appendix A. Our models use convective premixing (Paxton et al. 2019) and determine the convective boundary using the Ledoux criterion.

We use Equation (5) to calculate the Rayleigh numbers in subsurface convective layers caused by hydrogen, helium, and iron recombination (H i, He i, He ii, and the FeCZ). As they trace different ionization stages of H, He and Fe, these convective regions occur in specific temperature ranges (Cantiello & Braithwaite 2019). Moreover the FeCZ is only found in stars above a certain luminosity threshold, which in turns depends on the stellar metallicity (Cantiello et al. 2009).

Figure 1 shows the Rayleigh number normalized to the critical Ra_crit ( ∼ 10³; see Appendix D) on a Hertzsprung–Russel diagram (HRD) for the four different subsurface CZs with Milky Way metallicity Z = 0.014. A few features are striking. First, the He i CZ is almost always subcritical. This is a missing CZ: it is present in stellar models, and is still a superadiabatic layer in stars, but it likely does not undergo convective motions in most stars and has ∇ = ∇_rad because convection transports no heat there.

Zoom In Zoom Out Reset image size

Similarly the H i zone is subcritical for masses 2–3M_⊙ and the He ii zone is subcritical above 5M_⊙. The FeCZ by contrast is always supercritical where it exists.

The differences in Ra between these different layers can be explained by a few factors. First, the four CZs exist at very different densities. At lower density photons have a longer mean free path, and since photons carry both heat and momentum this results in a higher thermal and kinetic diffusivity. Examples of these trends are shown in the lower two panels of Figure 2. Moreover, the size of the convective region δ r is comparable to the local scale height, and this gets smaller for regions closer to the surface (Cantiello & Braithwaite 2019; see also the middle panels of Figure 2).

Zoom In Zoom Out Reset image size

A careful examination of Figure 1 reveals that there is always a subsurface CZ. At Milky Way metallicity (Z = 0.014) at least one of the four CZs is active for every mass we have examined and across the HRD. This is not the case for all Z, however.

Toward smaller Z the FeCZ moves up in $\mathrm{log}{T}_{\mathrm{eff}}$ and $\mathrm{log}L$ and a window opens where there are no subsurface CZs. This is shown in Figure 3, where we have plotted the window with no subsurface CZs for four different metallicities ranging from the SMC (Z = 0.002) up to the Milky Way (where the window vanishes). We call this the stability window, because it is a region where subsurface convection is stabilized by microscopic diffusivities.

Zoom In Zoom Out Reset image size

At SMC metallicity the stability window is wide, from 8 to 35M_⊙ and covering the whole main-sequence for most of those masses. As Z rises the upper end of the window comes down, tracing the FeCZ, first to 16M_⊙ at Large Magellanic Cloud (LMC) metallicity (Z = 0.006) then 10M_⊙ at Z = 0.01. Finally by Milky Way metallicity of Z = 0.014 the start of the FeCZ pulls down to 8M_⊙, matching the upper mass at which He ii convection occurs, and the window disappears. For convenience the same quantities are shown in Figure 4 versus spectroscopic luminosity ${\mathscr{L}}\equiv {T}_{\mathrm{eff}}^{4}/g$, where g is the surface gravity.

Zoom In Zoom Out Reset image size

Rotation modifies our scenario in two ways. First, it slightly changes the Rayleigh number in these subsurface superadiabatic zones by changing the hydrostatic structure of the star, but this is quite a small effect. The reason for this is that these zones are generated by opacity bumps which occur at nearly fixed T. Taking the density to be unchanged, and treating these zones as of order the pressure scale height h in radial extent, we find

$$\begin{eqnarray}&&{\mathrm{Ra}}_{\mathrm{rotating}}\sim \displaystyle \frac{{h}^{3}{g}_{\mathrm{eff}}}{\alpha \nu }\propto \displaystyle \frac{1}{{g}_{\mathrm{eff}}^{2}}\propto \displaystyle \frac{1}{{\left(1-{{\rm{\Omega }}}^{2}/{{\rm{\Omega }}}_{{\rm{c}}}^{2}\right)}^{2}},\end{eqnarray} \tag{ 7 }$$

where g_eff is the effective acceleration of gravity accounting for rotation, Ω is the angular velocity, and Ω_c is the critical rotation rate. Except very near critical rotation this is a small effect that does not compete with the many orders of magnitude variation in Ra across the HRD.

Second, and more importantly, rotation increases the critical Rayleigh number, stabilizing more stars against convection and growing the stability window. In the calculations above we assumed a very slow rotation rate of 10³ days. We do this because Ra_crit has a different value when no rotation is assumed than when an extremely slow rotation is assumed. Using slow but nonzero rotation rates slightly expands the stability windows. The effects of increasing rotation further are small, however, and we see no difference in the stability windows calculated with rotation periods varying from P = 100 days down to P = 2 days.

How much do these missing CZs matter? For stellar structure, not much: subsurface CZs carry very little heat, the temperature gradient ∇ ≈ ∇_rad, so in the absence of convective heat transport the near-surface temperature profiles of these stars are likely not that different. However, for more subtle measurements such as magnetism and photometric variability the properties of subsurface CZs are critically important.

For instance in Jermyn & Cantiello (2020, 2021) we made predictions of the interaction between surface magnetism and subsurface CZs. At Milky Way metallicities and above these are unaffected, but at lower metallicities (e.g., SMC and LMC) those predictions are significantly altered because there is a window in mass where no subsurface convection occurs (Figure 3), and we now expect:

• 1.  
  there should be no magnetic desert for stars in the stability window because there is no subsurface convection that can erase a magnetic field;
• 2.  
  a corollary is that stars should typically have stronger surface magnetic fields in the window than those just below the window in mass, because there can be no surface convective dynamo fields in these stars, only fossil fields, which tend to be stronger;
• 3.  
  if discrete absorption components (DACs; Kaper et al. 1996) are generated by magnetic spots then DACs should be steady rather than time-varying in the stability window, because there is no subsurface convection present to generate a disordered magnetic field;
• 4.  
  if wind clumping is seeded close to the stellar surface via convectively-driven fluctuations in density and velocity (e.g., Puls et al. 2006, 2008; Cantiello et al. 2009) then the clumping properties of stars in the window could be very different.

Similarly, the stability window provides an excellent opportunity to test theories of photometric variability in early-type stars. In particular there is active debate over the origin of surface low-frequency variability (Lecoanet et al. 2019). We are aware of three proposed explanations: subsurface convective motions (Cantiello et al. 2021), internal gravity waves excited by core convection (Bowman et al. 2019), and wind mass loss (Krtička and Feldmeier 2021). Stars inside the stability window, however, have no subsurface convection. Any photometric variability observed in the window, then, cannot come from convection and must come from either one of the two other proposals or else some process not yet considered. If the SLF phenomenon has multiple causes, we should expect a discontinuity on the edge of the window, and if it is monocausal we should either expect it to be indifferent to the window (ruling out a subsurface convective origin) or to be excluded from the window (requiring a subsurface convective origin).

The Flatiron Institute is supported by the Simons Foundation. This research was supported in part by the National Science Foundation under grant No. PHY-1748958. We thank Dominic Bowman for suggesting the addition of Figure 4 and Daniel Lecoanet for advice on computing the critical Rayleigh number in rotating convection.

Software: MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019,http://mesa.sourceforge.net), MESASDK 20190830 (Townsend 2019a, 2019b), matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011).

The MESA EOS is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), FreeEOS (Irwin 2004), HELM (Timmes & Swesty 2000), and PC (Potekhin & Chabrier 2010) EOSes.

Radiative opacities are primarily from OPAL (Iglesias & Rogers 1993, 1996), with low-temperature data from Ferguson et al. (2005) and the high-temperature, Compton-scattering dominated regime by Buchler & Yueh (1976). Electron conduction opacities are from Cassisi et al. (2007).

Nuclear reaction rates are from JINA REACLIB (Cyburt et al. 2010) plus additional tabulated weak reaction rates (Fuller et al. 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2000). Screening is included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from Itoh et al. (1996).

Models were constructed on the pre-main sequence with Y = 0.024 + 2Z and X = 1 − Y − Z and evolved from there. We neglect rotation and associated chemical mixing.

Computing the viscosity of a plasma is complicated. For simplicity we use the viscosity of pure hydrogen plus radiation, so that

$$\begin{eqnarray}&&\nu \approx {\nu }_{{\rm{H}}}+{\nu }_{\mathrm{rad}}.\end{eqnarray} \tag{ B1 }$$

The radiation component is (Spitzer 1962)

$$\begin{eqnarray}&&{\nu }_{\mathrm{rad}}=\displaystyle \frac{4{{aT}}^{4}}{15c\kappa {\rho }^{2}},\end{eqnarray} \tag{ B2 }$$

where a is the radiation gas constant. We obtain the hydrogen viscosity using the Braginskii–Spitzer formula (Braginskii 1957; Spitzer 1962)

$$\begin{eqnarray}&&{\nu }_{\mathrm{ie}}\approx 2.21\,\times \,{10}^{-15}\displaystyle \frac{{\left(T/{\rm{K}}\right)}^{5/2}}{(\rho /{\rm{g}}\,{\mathrm{cm}}^{-3})\mathrm{ln}{\rm{\Lambda }}}{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ B3 }$$

where $\mathrm{ln}{\rm{\Lambda }}$ is the Coulomb logarithm, given for hydrogen by (Wendell et al. 1987)

$$\begin{eqnarray}&&\mathrm{ln}{\rm{\Lambda }}=-17.9\,+\,1.5\mathrm{ln}\displaystyle \frac{T}{{\rm{K}}}-0.5\mathrm{ln}\displaystyle \frac{\rho }{{\rm{g}}\,{\mathrm{cm}}^{-3}}\end{eqnarray} \tag{ B4 }$$

for temperatures T < 4.2 × 10⁵ K and

$$\begin{eqnarray}&&\mathrm{ln}{\rm{\Lambda }}=-11.5+\mathrm{ln}\displaystyle \frac{T}{{\rm{K}}}-0.5\mathrm{ln}\displaystyle \frac{\rho }{{\rm{g}}\,{\mathrm{cm}}^{-3}}\end{eqnarray} \tag{ B5 }$$

for T > 4.2 × 10⁵ K. Corrections owing to different compositions are generally small relative to the many orders of magnitude we are interested in here. For instance the difference between pure hydrogen and a cosmic mixture of hydrogen and helium is under 30% (Balbus & Henri 2008).

The plotting scripts used in this work, as well as the MESA run directories used in producing Figure 2, may be found in this GitHub repository in the commit with short-sha d7cdf2f. The inlists and run scripts used in producing the HR diagrams in this work are available in this other GitHub repository on the main branch in the commit with short-sha 964e07d. The data those scripts produced are available in Jermyn et al. (2022).

The critical Rayleigh number depends on both the value of the Taylor number (rotation) and the Prandtl number (diffusivity ratio),

$$\begin{eqnarray}&&\mathrm{Ta}\equiv \displaystyle \frac{4\delta {r}^{4}{{\rm{\Omega }}}^{2}}{{\nu }^{2}}\end{eqnarray} \tag{ D1 }$$

and

$$\begin{eqnarray}&&\Pr \equiv \displaystyle \frac{\nu }{\alpha }.\end{eqnarray} \tag{ D2 }$$

Recent work has explored the intricacies of convective onset in rotating, spherical domains (Net et al. 2008; Garcia et al. 2018). For simplicity, we use the derivation of Chandrasekhar (1961) to determine convective onset. For a given value of Ta and Pr, we use his Ch. 3, Equations (222) and (227) to respectively calculate the critical Rayleigh number of an overstable onset and direct onset,

$$\begin{eqnarray}&&{\mathrm{Ra}}_{\mathrm{crit},\mathrm{over}}=2{\pi }^{4}\displaystyle \frac{(1+\Pr )}{x}\left({\left(1+x\right)}^{3}+\displaystyle \frac{{\Pr }^{2}\mathrm{Ta}/{\pi }^{4}}{{\left(1+\Pr \right)}^{2}}\right),\end{eqnarray} \tag{ D3 }$$

$$\begin{eqnarray}&&{\mathrm{Ra}}_{\mathrm{crit},\mathrm{direct}}={\pi }^{4}\left(\displaystyle \frac{{\left(1+x\right)}^{3}+\mathrm{Ta}/{\pi }^{4}}{x}\right),\end{eqnarray} \tag{ D4 }$$

where x is a nondimensional wavenumber of convective onset, which must be compared to a reference value (Equation (225)),

$$\begin{eqnarray}&&{x}^{* }={\left(\displaystyle \frac{\mathrm{Ta}(1-\Pr )}{{\pi }^{4}(1+\Pr )}\right)}^{1/3}-1.\end{eqnarray} \tag{ D5 }$$

If x < x^*, Ra_crit = Ra_crit,over, otherwise Ra_crit = Ra_crit,direct. For a given value of Ta and Pr, we evaluate Ra_crit for x ∈ [10⁻³, 10], and take the minimum value to be the true value of Ra_crit.