The Initial Mass–Final Luminosity Relation of Type II Supernova Progenitors: Hints of New Physics?

The solver we use to integrate the differential equations describing the stellar equilibrium structure and the chemical evolution is based on a classical Henyey method. It was originally derived from the former FRANEC code (Chieffi & Straniero 1989). The last version of the FuNS allows us to choose different degrees of coupling between the equations describing the physical structure and those describing the chemical evolution, as due to nuclear burning and mixing. For the models discussed in this work, we have adopted a scheme in which the stellar structure equations (i.e., hydrostatic equilibrium, mass continuity, energy conservation, and energy transport) and the chemical evolution equations due to the nuclear burning are solved simultaneously. Then, the stellar zones that are thermally unstable are separately mixed. When dealing with the advanced stages of the evolution of a massive star, such a scheme is a good compromise between accuracy, velocity, and stability of the numerical solutions. Our code also allows three different choices of the time-dependent mixing scheme, namely, diffusive, advective, or the nonlocal algorithm described in Straniero et al. (2006). For the purpose of this paper the three schemes are equivalent. In order to calculate the temperature gradient and the average turbulent velocity (or diffusion coefficient) in the convective regions, we use, as usual, the mixing-length theory. In particular, we follow the scheme described in Cox & Giuli (1968). The mixing-length parameter (α_ML) has been calibrated by means of a standard solar model (Piersanti et al. 2007). Note that the value of α_ML depends on the adopted solar composition. Since we adopt the Lodders et al. (2009) composition, we found α_ML = 1.9.

In general, the convective boundaries are fixed according to the Ledoux criterion (see Section 2.2). In addition, during the He-burning phase, the external border of the convective core is found by imposing the condition of marginal stability, i.e., ∇_rad strictly equal to ∇_ad (for more details see Straniero et al. 2003). Except for the models described in Section 2.3, no convective overshoot is applied.

The FuNS code also accounts for stellar rotation (Piersanti et al. 2013). As usual, we assume shellular rotation, i.e., constant angular velocity and composition on the isobaric surfaces (see, e.g., Maeder 2009). As illustrated in Section 2.4, the most important instabilities induced by rotation, which are responsible for transport of angular momentum and mixing, are considered. The potential effects of the magnetic field on rotation are ignored.

The FuNS code has been optimized to handle large nuclear networks, such as those required to follow the neutron capture nucleosynthesis in asymptotic giant branch (AGB) stars (see Straniero et al. 2006). However, a minimal nuclear network that includes all the reactions providing the major contribution to the energy production is enough for the purpose of the present paper and allows us to reduce the computational time consumption. Hence, we have considered 33 isotopes, from ¹H to ⁵⁶Ni, coupled by 19 reactions for the H burning (those describing a full pp-chain and CNO cycle) and an α-chain of 51 reactions, plus the ¹²C+¹²C, ¹²C+¹⁶O, and ¹⁶O+¹⁶O, for the more advanced burnings. The reaction rates are from the STARLIB database (Sallaska et al. 2013). Most of the H- and He-burning reaction rates in this database are based on available experimental data. Rate estimates of reactions for which no experimental information exists, or extrapolation to low temperatures at which no experimental rates exist, are obtained by means of statistical (HauserFeshbach) models of nuclear reactions, computed using the code TALYS (Goriely et al. 2008). Finally, electron screening corrections are computed according to Dewitt et al. (1973), Graboske et al. (1973), and Itoh et al. (1979).

Other details on the FuNS input physics are provided here:

• 1.  
  Neutrino Energy Loss:

  • (a)  
    Plasma—Haft et al. (1994)
  • (b)  
    Photo—Itoh et al. (1996a)
  • (c)  
    Pair—Itoh et al. (1996a)
  • (d)  
    Bremsstrahlung—Dicus et al. (1976)
  • (e)  
    Recombination—Beaudet et al. (1967)
• 2.  
  Radiative Opacity:

  • (a)  
    Alexander & Ferguson (1994)+Iglesias & Rogers (1996)+Magee et al. (1995)
• 3.  
  Electron Conductivity:

  • (a)  
    Potekhin et al. (1999), Potekhin (1999)
• 4.  
  Equation of State:

  • (a)  
    Rogers et al. (1996) + Straniero (1988) (see also Prada Moroni & Straniero 2002)
• 5.  
  Mass Loss:

  • (a)  
    Nieuwenhuijzen & de Jager (1990)
• 6.  
  External Boundary Conditions:

  • (a)  
    scaled solar T(τ) relation (Krishna Swamy 1966).

As is well known, the luminosity of an RSG and its He core mass are closely related. Indeed, the mass of the core determines the physical conditions (T and ρ) of the H-burning shell, which is the source of energy that sustains the shining of these stars. After the He-burning phase, the central temperature rapidly increases, becoming larger than ∼5 × 10⁸ K. At that temperature, the production of neutrinos, as due to the Compton scattering and e⁺e⁻ annihilation, becomes very efficient. As a result, the advanced stages of the evolution of a massive star are controlled by the neutrino energy loss that largely overcomes the photon energy loss. Such an occurrence causes a rapid drop of the evolutionary timescale, which becomes much smaller than the H-burning timescale. Hence, since the exhaustion of the central carbon and until the final collapse, the He core mass and, in turn, the stellar luminosity stop to grow. In the rest of this section we will review the most important physical processes that determine the growth of the He core mass up to this constant value and, hence, their impact on the pre-explosive luminosity.

The new models presented here have masses between 11 and 30 M_⊙. If not differently specified, we have assumed solar composition, i.e., Z = 0.014 and Y = 0.27. All the models, except those with M = 11 M_⊙, have been evolved up to the Si burning, when the central temperature is ∼4 × 10⁹ K. Instead, the calculation of the M = 11 M_⊙ models has been stopped just after the off-center Ne ignition.

Among the principal uncertainties affecting H and He-burning models of massive stars, those related to the treatment of turbulent mixing induced by convection still remain the most debated. In this context, a long-standing problem concerns the criterion adopted to establish whether a nonhomogeneous zone is unstable against convection. In particular, in the case of a negative molecular weight gradient, it may happen that a thermally unstable stratification is stabilized against adiabatic convection by a gradient in composition. This phenomenon, often called semiconvection, is particularly relevant for massive stars (see Langer et al. 1985). In this case, the widely adopted Schwarzschild criterion, which ignores the stabilizing effect of the molecular weight gradient, appears inadequate to model the mixing and the consequent heat transport. On the other hand, the Ledoux criterion, which includes composition effects, hampers mixing in semiconvective regions.

Concerning massive stars, semiconvective zones first appear during the H-burning phase, when the convective core recedes, as a consequence of the conversion of H into He, and a region of negative μ gradient is left outside. If according to the Schwarzschild criterion these layers are efficiently mixed, a prompt He ignition occurs after the exhaustion of the central hydrogen, when the star is still a compact blue supergiant. Then, during most of the He-burning phase, the star remains hot, and only toward the end of this phase does it eventually become an RSG. On the contrary, if according to the Ledoux criterion the mixing is inhibited in the semiconvective layers, after the central H exhaustion the stellar envelope expands and the star directly evolves into a cool RSG. In this case the He ignition takes place on the RSG branch.

A simple analysis of stability (see Kippenhahn & Weigert 1990) demonstrates that even if the normal convective modes are stable, oscillatory modes can exist that may induce some mixing in the semiconvective layers. On the other hand, observations of massive stars in the Milky Way and in the Magellanic Clouds clearly favor a prompt evolution to the RSG phase (Stothers & Chin 1992b, 1994). Figure 1 shows some evolutionary tracks of the models we obtained by varying the damping factor (η) of the semiconvective velocity (v_sc) with respect to the fully convective velocity (v_c). In practice, we have assumed that v_sc = ηv_c and varied the η parameter from 0 to 10⁻³. Note that the case η = 0 corresponds to the bare Ledoux criterion, while for η = 1 the effect of the molecular weight gradient is fully suppressed. Only models with η ≤ 10⁻⁵ become RSGs before the He ignition. Note, however, that the variation of the final luminosity is rather small, namely, ${\rm{\Delta }}\mathrm{log}L\lt 0.03$. In the rest of the paper, we will adopt the Ledoux criterion. The resulting initial mass–final luminosity relation is reported in Figure 2. For comparison, other relations, as obtained from nonrotating models with moderate or no overshoot, are also shown, namely, Kepler (Woosley et al. 2002), GENEVA (Meynet et al. 2015), and MESA (Farmer et al. 2016). The small differences are likely due to the adopted treatments of overshoot, semiconvection, and mass loss. The vertical line marks the observed upper limit for the luminosity of SN II progenitors (Smartt 2015). The corresponding masses are reported in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The so-called convective core overshoot is another phenomenon that may affect the final He core mass and, in turn, the final luminosity of a massive star. Indeed, owing to the inertial motion of the convective cells, a transition zone, rather than a sharp discontinuity, should exist between stable and unstable regions. In principle, the extension of this zone and the mixing efficiency within it depend on the convective velocity at the boundary layer and on the steepness of the entropy gradient. Therefore, a unique prescription to be applied to all the convective boundaries appears unrealistic. In most cases hints on the overshooting extension may be provided by the observations of the stellar properties affected by this phenomenon. Concerning massive stars, Stothers & Chin (1992a; but see also Maeder & Meynet 1989) have shown that the convective core overshoot cannot be too large; otherwise, it would be in conflict with observations. Nevertheless, in order to quantify the effect of this phenomenon on the initial mass–final luminosity relation, we have computed two additional sets of models with convective core overshoot. In practice, according to current hydrodynamical simulations of stellar convection (see, e.g., Freytag et al. 1996), we have extended the mixing outside the fully convective core of the H-burning models by imposing an exponential decline of the convective velocity, namely,

$$\begin{eqnarray}&&{v}_{\mathrm{ov}}={v}_{\mathrm{cct}}\exp \left(-\displaystyle \frac{\delta r}{\beta {H}_{p}}\right),\end{eqnarray} \tag{ 1 }$$

where v_cct is the convective velocity at the top of the fully convective core,⁹ δr is the radial distance from the convective border, and H_p is the pressure scale height. β is a free parameter. The results are shown in Figure 3 for β = 0.01 and 0.02 (short- and long-dashed lines, respectively). For a given stellar mass, the final He core mass is larger in the case of convective overshoot, and, in turn, the pre-explosive models are brighter than those obtained by neglecting the overshoot mixing (solid line). For comparison, we have also reported the initial mass–final luminosity relations from two extant sets of stellar models with overshoot, namely, ORFEO (Limongi & Chieffi 2018) and STARS (Eldridge & Tout (2004).¹⁰ Note that a moderate overshoot (of the order of 0.1H_p) has been adopted in the ORFEO models, while the extra-mixing zone is rather extended in STAR models (∼0.5H_p). The corresponding upper mass bounds for the observed progenitors of SNe II are also listed in Table 1. In the case of large overshoot, this maximum mass would be particularly small, thus increasing the tension with the masses estimated by means of LC best fitting or from the oxygen yields measured in the nebular phase (Davies & Beasor 2018; Morozova et al. 2018). In addition, due to the higher luminosity, the mass loss is larger in models with convective core overshoot. As a result, also the maximum mass of collapsing RSGs is smaller than that predicted by models without overshoot.

Zoom In Zoom Out Reset image size

Rotational velocities of the order of 100–200 km s⁻¹ are commonly observed in MS massive stars (Hunter et al. 2008). As is well known, 1D hydrostatic codes may easily account for shellular rotation (Endal & Sofia 1976). There are two major consequences of rotation, namely, (i) the lifting due to the centrifugal force, and (ii) the meridional circulation that arises because of the deviations from thermal equilibrium occurring in rotating structures (Von Zeipel's paradox). The lifting effect of rotation implies more expanded and cooler stellar structures. In addition, the mixing induced by meridional circulation modifies the mean molecular weight and the opacity of the envelope. As a result, larger convective cores are expected in H- and He-burning stars with rotation. However, a quantitative estimation of these effects requires a reliable description of the angular momentum redistribution as due to both dynamical and secular instabilities. Moreover, angular momentum loss driven by stellar wind, which is particularly intense when the star is an RSG, must be considered.

In the present version of the FuNS code the transport of angular momentum is treated as a diffusive process (see details in Piersanti et al. 2013). The most important instabilities induced by rotation are considered. The Eddington–Sweet (ES) and, to a less extent, the Goffrei–Shubert–Fricke (GSF) instabilities may produce significant effects on the evolution of massive stars (Heger & Langer 2000). According to Kippenhahn & Weigert (1990), the ES circulation velocity is given by

$$\begin{eqnarray}&&{v}_{\mathrm{ES}}=\displaystyle \frac{{{\rm{\nabla }}}_{\mathrm{ad}}}{\delta ({{\rm{\nabla }}}_{\mathrm{ad}}-{\rm{\nabla }})}\displaystyle \frac{{\omega }^{2}{r}^{2}L}{{\left({Gm}\right)}^{2}}\left(\displaystyle \frac{2\epsilon {r}^{2}}{L}-\displaystyle \frac{2{r}^{2}}{m}-\displaystyle \frac{3}{4\pi {r}^{2}\rho }\right),\end{eqnarray} \tag{ 2 }$$

where ω is the angular velocity, is the rate of energy loss per unit mass, $\delta =-\tfrac{\partial \mathrm{ln}\rho }{\partial \mathrm{ln}T}$, and ∇ and ∇_ad are the temperature gradient and the adiabatic gradient, respectively. L, m, r, G, and ρ have the usual meaning. As for convection, a negative molecular weight gradient reduces the efficiency of the meridional circulation. According to Kippenhahn (1974), we model the effect of the μ-gradient by defining an equivalent μ current, which works against the ES circulation:

$$\begin{eqnarray}&&{v}_{\mu }={f}_{\mu }\displaystyle \frac{{H}_{P}}{\tau }\displaystyle \frac{\phi {{\rm{\nabla }}}_{\mu }}{{\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}},\end{eqnarray} \tag{ 3 }$$

where H_P is the pressure scale height, $\phi =\tfrac{\partial \mathrm{ln}\rho }{\partial \mathrm{ln}\mu }$, ${{\rm{\nabla }}}_{\mu }=\tfrac{\partial \mathrm{ln}\mu }{\partial \mathrm{ln}P}$ is the μ gradient, and τ is the thermal relaxation timescale (also called the Kelvin–Helmholtz timescale). f_μ is a (free) parameter used to tune the strength of the μ gradient barrier. Hence, the effective circulation velocity is $\left|{v}_{\mathrm{MC}}\right|=\left|{v}_{\mathrm{ES}}\right|-\left|{v}_{\mu }\right|$. Such a velocity is used to calculate the angular momentum diffusion coefficient, while a fraction of it, i.e., v_rot = f_c v_MC, is used to calculate the mixing induced by rotation. In principle, both f_μ and f_c should vary between 0 and 1. These parameters could be calibrated by requiring that models reproduce some observables, such as the modification of the surface composition due to mixing induced by the meridional circulation (see, e.g., Heger & Langer 2000; Maeder 2009; Limongi & Chieffi 2018). Alternatively, asteroseismology may provide more stringent constraints to the efficiency of the angular momentum transport. Indeed, it allows us to measure core rotation rates of stars in different evolutionary phases (see Eggenberger et al. 2019, and references therein). The first results of the application of this method, mainly concerning low- and intermediate-mass stars, seem to indicate that the angular momentum redistribution is very efficient. In general, this occurrence should reduce the effects of rotation. Unfortunately, a clear and unique prescription valid for massive stars is still lacking. This problem represents a further uncertainty affecting massive-star models. In Figure 4 we report the evolutionary tracks of a 20 M_⊙ star as obtained under different assumptions about f_c, f_μ, and the initial (uniform) rotational velocity (v_ini). In general, the final He core mass is larger than that of nonrotating models. As a consequence, a rotating SN progenitor is brighter than a nonrotating one. However, the extent of this effect depends on the assumed set of rotation parameters.

Zoom In Zoom Out Reset image size

In Figure 5, the variation of the initial mass–final luminosity relation for models with rotation is illustrated. The FuNS relation for v_ini = 200 km s⁻¹, f_c = 0.04, and f_μ = 0 is shown, together with that from the ORFEO database (Limongi & Chieffi 2018, v_ini = 150 km s⁻¹, f_c = 1, and f_μ = 0.05) and that from the GENEVA code (Meynet et al. 2015, v_ini = 0.4v_crit). Note that a spread of v_ini implies a spread of the final luminosity. In other words, an observed pre-explosive luminosity does not correspond to a unique value of the initial mass. Instead, a range of masses is possible, depending on the (unknown) initial rotational velocity.

Zoom In Zoom Out Reset image size

Another important consequence of rotation concerns the mass loss (Meynet et al. 2015). Indeed, owing to the higher luminosity, the mass-loss rate is larger in rotating models. As a result, we find that rotating models (v_ini = 200 km s⁻¹) whose initial mass is M ≥ 15 M_⊙ leave the RSG branch before the final core collapse, and for this reason they cannot be progenitors of SNe IIP/L.

The effect of the mass loss on the evolution of massive stars has been deeply investigated by several authors (for a detailed analysis of the mass-loss effects see Meynet et al. 2015). The total amount of mass lost before the final collapse depends on the initial mass, the initial rotation velocity, and the extension of the convective core overshoot. For instance, based on nonrotating FuNS models, a 13 M_⊙ star loses less than 2 M_⊙, while a 25 M_⊙ star is expected to lose about half of its initial mass. In general, this theoretical prediction is in good agreement with the observed mass ejecta of SN IIP progenitors (Morozova et al. 2018). As mentioned in the previous section, the higher RSG luminosity developed in the case of rotation implies a stronger mass-loss rate. For the same reason, also models with convective overshoot develop higher mass-loss rates. In general, the larger the mass-loss rate, the lower the final luminosity. We find, however, that the final luminosity does not depend on the mass lost during the RSG phase, but mostly on the amount of mass lost during the MS (see also Section 2.9). Therefore, since the mass-loss rate is much lower in MS than in the RSG branch, the sensitivity of the initial mass–final luminosity on the adopted mass-loss rate is rather small. For instance, we find that in a nonrotating model of 18 M_⊙ with mass loss, the final luminosity is only 5% lower than that of a model without mass loss.

On the other hand, if the RSG mass loss is strong enough to erode the H-rich envelope, the model leaves the RSG branch prior to the final core collapse, and the resulting SN cannot be of Type IIP/IIL. This occurrence implies a maximum RSG mass-loss rate for a progenitor of these SN types.

In this context, recent evidences of late-time enhanced mass loss in progenitors of some bright SNe IIP/IIL (Moriya et al. 2011; Morozova et al. 2017; Yaron et al. 2017) are not expected to modify the initial mass–final luminosity relation for these stars. However, to explain pre-explosive outbursts commonly observed in SNe IIn, Fuller (2017) has recently suggested that gravity waves generated by vigorous convection during late-stage nuclear burning could transport energy upward, causing a heating of the H-rich envelope and, in turn, an increase of the final luminosity. Although firm conclusions about the occurrence and, eventually, the strength of such a heating process in normal SN IIP progenitors have not been established yet, we may exclude that it could solve the discrepancy between the masses estimated from the SN LC and those derived from the pre-explosive luminosity. On the contrary, this phenomenon would increase the tension. Indeed, the masses from LCs are usually larger than those estimated from the initial mass–final luminosity relation (see the discussion in Section 4).

The rate of neutrino energy loss determines the time at which the evolutionary timescale becomes so short that the He core mass stops to grow and a luminosity freeze-out occurs. In particular, a rate larger than that commonly estimated would imply fainter SN II progenitors. The current rates are computed on the basis of the Weinberg–Salam electro-weak theory (see, e.g., Raffelt 1996). In the past 30 yr, precision experiments have tested this quantum field theory at the level of 1% or better. Note, however, that the neutrino rates depend on the Weinberg angle, whose value is not predicted by the electro-weak theory. In the present calculations we have used the value suggested by Itoh et al. (1996a), namely, ${\sin }^{2}{\theta }_{W}=0.2319$, which is very similar to that reported in the latest compilation of the Particle Data Group¹¹ (0.23155). We have also verified that negligible variations of the final luminosity are obtained when the CODATA 2014 value (0.2223; Mohr et al. 2016) is adopted. As a whole, the accuracy of the Compton and pair neutrino energy-loss rates, which are the dominant neutrino processes in the present stellar model calculations, is better than 5% (Itoh et al. 1996a, 1996b). Such an uncertainty has little effect on the pre-explosive luminosity.

On the other hand, some deviations from the standard electro-weak theory are not completely excluded. In particular, the existence of a nonzero neutrino magnetic moment would enhance the stellar energy loss. Experimental and astrophysical constraints provide upper bounds for this quantity. The more stringent experimental constraint is <2.9 × 10⁻¹¹μ_B. It has been obtained by means of reactor neutrinos (90% CL; Beda et al. 2013). A more stringent astrophysical constraint, as based on the luminosity of the red giant branch (RGB) tip of galactic globular clusters, gives <2.6 × 10⁻¹²μ_B (68% CL; Viaux et al. 2013a). The influence of a nonzero neutrino magnetic moment on the evolution of massive stars has been investigated by Heger et al. (2009). They found that models of initial mass above 15 M_⊙ are practically insensitive to a neutrino magnetic moment lower than 10⁻¹¹ μ_B. However, they do not consider the sensitivity of the final luminosity on the enhanced neutrino rate. Therefore, we have computed a model of 20 M_⊙ by assuming μ_ν = 5 × 10⁻¹¹μ_B, and we find that the resulting final luminosity is just 2% lower than that of the corresponding μ_ν = 0 model.

Similarly to neutrinos, also a variation of the C-burning temperature could anticipate or delay the time of the luminosity freeze-out. Indeed, owing to the strong sensitivity of the neutrino production rates on the temperature, a high C-burning temperature would anticipate the freeze-out, while the opposite occurs in the case of a low C-burning temperature. This temperature depends on the rate of the ¹²C+¹²C reaction and on the amount of carbon left in the core after the He burning.

A large uncertainty affects the ¹²C+¹²C reaction rate for T ≤ 1 GK (for a recent reanalysis, see Zickefoose et al. 2018, and references therein). At present, direct measurements of the cross section of this reaction are only available for energies >2.1 MeV. Owing to a molecular-like structure of the ²⁴Mg compound nucleus, unknown resonances are possible at lower energy. The existence of these molecular states would substantially enhance the low-energy cross section, thus reducing the C-burning ignition temperature. In contrast, a hindrance model (Jiang et al. 2007) predicts a steep drop of the low-energy S(E) factor¹² and, hence, would imply a larger C-burning ignition temperature. In this context, some hints may be obtained with indirect measurements. Tumino et al. (2018) have recently presented new results based on the so-called Trojan Horse Method (THM). The THM is an indirect method, which allows us to measure the S(E) factors of charge particle reactions down to astrophysical relevant energies, where data from direct methods are not available because of the strong Coulomb barrier and the consequent very small cross sections. In the range of energy of more interest to the C-burning stellar rate, Tumino et al. found a complex resonant pattern but no evidence of hindrance. However, the interpretation of these indirect measurements requires reliable theoretical models and accurate calibrations. On this basis, Mukhamedzhanov et al. (2018) have argued that the THM overestimates the actual S(E) factor of the ¹²C+¹²C reaction. Therefore, more accurate experiments and theoretical models are needed to fix this important open issue of massive-star evolution. Meanwhile, we have checked the potential impact of the THM result by computing a 20 M_⊙ model with the Tumino et al. rate. We find that the final luminosity is only 2% higher than that obtained with the widely adopted Caughlan & Fowler (1988) rate. An illustration of the impact of the THM rate on stellar models with lower mass, those that ignite C and Ne in degenerate conditions, can be found in Straniero et al. (2019).

The C-burning rate also scales with the square of the carbon abundance. This carbon has been previously produced during the He-burning phase. The reactions involved in the C production are the triple-α, which represents the production channel, and the ¹²C(α, γ)¹⁶O, which represents the destruction channel. While the rate of the triple-α reaction is known within 10% at the temperature of the He burning in massive stars (Fynbo et al. 2005), the ¹²C(α, γ)¹⁶O reaction rate is still rather uncertain. The rate in the STARLIB compilation is based on the Kunz et al. (2002) study. Recently, a complete reanalysis has been reported by deBoer et al. (2017). Based on an R-matrix fit of all the available experimental data, they find a reaction rate slightly smaller than that of Kunz et al. As a consequence, a larger amount of C is expected in the center of a star at the end of the He-burning phase. In Table 2, we compare the results we obtain for a set of 20 M_⊙ models computed under different assumptions for the ¹²C(α, γ)¹⁶O reaction rate. Within the uncertainty quoted by deBoer et al. (2017), we find a 22% variation of the central carbon mass fraction at the end of the He-burning phase. Such a variation of the C abundance produces sizable effects on the C burning (Weaver & Woosley 1993; Imbriani et al. 2001). Indeed, a lower C abundance implies higher temperatures, smaller convective cores, and, in turn, shorter C-burning lifetimes. It affects the compactness of the pre-explosive stellar structure, as well as the hydrostatic and explosive nucleosynthesis. In principle, these occurrences could also affect the final luminosity. However, as noted by Straniero et al. (2016), a larger ¹²C(α, γ)¹⁶O rate also implies a larger He-burning lifetime. It occurs because the ¹²C(α, γ)¹⁶O reaction consumes one-third of He fuel compared to the triple-α but releases a similar amount of nuclear energy. As a consequence, the larger the ¹²C(α, γ)¹⁶O reaction rate, the slower the He consumption. In practice, this variation of the He-burning lifetime compensates the effect of the variation of the central C abundance, so that a small variation of the final He core mass is found. As reported in the third column of Table 2, when the¹²C(α, γ)¹⁶O is varied within the uncertainty quoted by deBoer et al., such an occurrence implies negligible variations (<0.2%) of the final luminosity.

Table 2.  Central C Mass Fraction at the End of the He Burning (X_c) and Pre-explosive Luminosity for 20 M_⊙ Models Computed with Different ¹²C(α, γ)¹⁶O Reaction Rates

Rate	Xc	$\mathrm{log}L/{L}_{\odot }$
Kunz et al. (2002)	0.285	5.1739
deBoer et al. (2017)	0.309	5.1739
High	0.275	5.1740
Low	0.342	5.1746

Note. Low and high refer to the lower and upper bounds quoted by deBoer et al. (2017), respectively.

Download table as:  ASCIITypeset image

SNe II are found preferentially on the disks of spiral galaxies. Moreover, the host galaxies where the progenitors can be resolved belong to the local universe. For these reasons, the parent stellar populations should be rather similar. Nevertheless, a certain spread of the initial composition cannot be excluded. Usually a solar composition is adopted. An increase of the initial He mass fraction and a reduction of the initial metallicity both imply smaller envelope opacity and less efficient H burning. The resulting stellar structures are more compact and brighter. In Figure 6 we compare the initial mass–final luminosity relations obtained by increasing the initial He mass fraction, from Y = 0.27 to Y = 0.32, and by decreasing the initial metallicity, from Z = 0.014 down to Z = 0.006, to the standard one. The masses corresponding to a pre-explosive luminosity $\mathrm{log}L/{L}_{\odot }=5.1$ are reported in Table 1.

Zoom In Zoom Out Reset image size

Summarizing, variations of Y and Z have opposite effects on the final luminosity. Notice that a positive correlation between the He mass fraction and the metallicity is a natural consequence of the chemical evolution of galaxies (Pagel & Portinari 1998). For this reason a spread of the initial composition should not produce a substantial modification of our findings.

More than 50% of massive stars are in binary (or multiple) stellar systems, and most of these systems are close enough that one or more mass transfer episodes through Roche lobe overflow may occur during the pre-SN evolution (Sana et al. 2012). In this case, not only mass but also angular momentum is exchanged between the components of the binary. Tidal as well as dissipative forces may also induce a synchronization of the stellar rotation rate and the orbital motion. These occurrences are expected to modify the standard paradigm of massive stars, as based on models of single-star evolution. Close-binary evolution may indeed produce a variety of different evolutionary channels (Paczyński 1967; Podsiadlowski et al. 1992; Nomoto et al. 1995; Vanbeveren et al. 1998; Eldridge et al. 2008; Yoon et al. 2010). For this reason, the binarity is often invoked to explain peculiar objects or, more in general, to interpret the diversity of the core-collapse SN inventory. Famous examples are (i) the anomalous SN 1993J, whose progenitor was supposed to be a K-type supergiant hosted in a binary system with a hotter companion, and (ii) SN 1987A, whose progenitor was clearly identified as a blue supergiant, possibly the secondary component of a binary with initial mass ratio close to 1 (for a review see Smartt 2009).

In general, three cases of mass transfer are usually distinguished, depending on the evolutionary status of the Roche lobe filling star. In case A a star fills the Roche lobe during the MS, while in the B and C cases the mass exchange occurs when the donor is becoming or already is an RSG. Then, the mass lost by the donor flows through the internal Lagrangian point (L1) and eventually feeds an accretion disk around the smaller companion. It should be noted that MS massive stars are rather compact objects with radii of just a few R_⊙, so that the separation between the two components must be particularly short for a case A Roche lobe overflow. Such kinds of binaries are probably a small fraction of all the binary systems harboring massive stars. More frequent should be the other two cases, in which the Roche lobe filling star is a red or yellow supergiant.

Among the various close-binary evolutionary channels, here we are interested in those that at the end of the mutual interaction leave stars with an H-rich envelope whose mass is larger than about 2–3 M_⊙, so that they could be still considered as possible progenitors of SNe IIP.¹³ Let us discuss the mass donors first. In order to investigate the consequences of an enhanced mass loss due to a binary interaction on the pre-explosive luminosity, we have calculated a small set of 20 M_⊙ models, in which 5 M_⊙ of the H-rich envelope are artificially removed at a quite fast rate. The fast mass-loss episode is switched on at different epochs. Then, when the total mass of 15 M_⊙ is attained, the mass loss is switched off and the evolution has been continued up to the Si burning. No further mass transfer episodes are considered during the remaining part of the evolution. In addition, we ignore the possibility of a dynamical mass transfer with the consequent common envelope episode, a phenomenon that may possibly occur in cases B and C of Roche lobe overflow, due to the presence of an extended convective envelope in the donor star. Nonetheless, this simple experiment can provide a reasonable estimation of the effects of the binary interaction on the final luminosity. Four cases are here illustrated, namely, case A-1 and A-2, in which the fast mass loss takes place during the MS, when the central H mass fractions are X = 0.49 and X = 0.04, respectively, and the corresponding radii are 7 and 12 R_⊙; case B, in which the mass-loss phase occurs when the star, after the MS, moves toward the RSG and the radius attains 281 R_⊙; and case C, in which the fast mass loss starts after the He-burning phase, when the radius attains 10³ R_⊙. The resulting evolutionary tracks are shown in Figure 7. Note that in cases A-2, B, and C the final luminosity practically coincides with that of a single-star model with 20 M_⊙ (the initial mass value), while in case A-1 it is intermediate between that of the 15 and 20 M_⊙ single-star models. In other words, the initial mass–final luminosity relation for the donor stars is insensitive to the mass transfer, provided that it occurs close to the end or after the MS. In all cases, if most of the mass lost by the donor is subtracted by the accretor (conservative case), the initial mass obtained by means of the initial mass–final luminosity relation for single stars will be systematically larger than that estimated by means of the LC-fitting method. Among the eight SN progenitors whose mass has been derived with both methods, only one, i.e., SN 2012ec, shows a similar behavior, while the opposite is found for the others (see the discussion in Section 4).

Zoom In Zoom Out Reset image size

Concerning the accretors, also in this case the consequences of the mass transfer depend on the evolutionary stage at the beginning of the accretion phase. According to Podsiadlowski et al. (1992), if the accretion occurs when the star is on the MS, the post-accretion evolution will be very similar to that of a single star whose mass is equal to the total mass, i.e., initial mass plus the accreted mass. On the other hand, when the accretion phase starts after the exhaustion of the central H, the pre-explosive structure will be a blue or yellow supergiant (see Figure 10 in Podsiadlowski et al. 1992). Likely, this is a consequence of the decrease of the core mass-total mass ratio (${\xi }_{c}={M}_{c}/M$) caused by the mass accretion process. Summarizing, only MS accretors may produce RSG progenitors and, possibly, generate an SN IIP. In this case the mass estimated by means of the initial mass–final luminosity relation for single-star progenitors should coincide with that obtained with the LC-fitting method.

Finally, the possible gain of angular momentum, as due to the spin–orbit coupling, would have the same consequences illustrated in Section 2.4, namely, larger core masses, brighter progenitors, and, in turn, lower masses estimated by means of the initial mass–final luminosity relation.

The Henyey method employed in the majority of the stellar evolution codes is a first-order implicit method, and the accuracy of the solutions depends on the assumed mass and time resolutions. In the FuNS code, we adopt an adaptive algorithm to select the grid of mass shells. This algorithm controls the variations between adjacent mesh points of luminosity, pressure, temperature, mass, and radius. For instance, if A is one of L, P, T, M, or R, it should be

$$\begin{eqnarray}&&0.8\times \delta \lt \displaystyle \frac{\left|A(N+1)-A(N-1)\right|}{A(N))}\lt \delta ,\end{eqnarray} \tag{ 4 }$$

for each shell N. In the present work we have used δ = 0.05 everywhere, except for the shells located around some critical points, such as the boundaries of the convective zones, for which we use δ = 0.005. This choice implies about 1000 mesh points for MS models and up to 5000 for the most advanced models. Recently, Farmer et al. (2016) argue that a minimum mass resolution of 0.01 M_⊙ is necessary to achieve convergence in the He core mass within 5%. However, as noted by Sukhbold et al. (2018), the effects of an improved resolution on the observable quantities, such as luminosity or effective temperature, are generally small. In particular, they found a weak trend to produce slightly larger He core masses and, in turn, brighter pre-SN models, when the zoning is finer. To check this occurrence, we have calculated a few additional evolutionary tracks with the further constraint that the difference in mass between two adjacent mesh points cannot exceed 0.01 M_⊙. The resulting pre-explosive luminosities are less than 1% larger than that obtained under our standard choice for the grid resolution.

As discussed above, the enormous majority of results on the impact of axions on stellar evolution involves low-mass stars, that is, stars of mass similar to or smaller than the mass of the Sun. Part of the reason is that low-mass stars are the most abundant constituents of the galactic stellar populations, allowing a quite larger statistics. However, the energy-loss rate via axions is very sensitive to the temperature (see the Appendix), and this favors the hotter environment of the core of more massive stars. In spite of that, only a few attempts have been made to study the axion effects on the evolution of stars a few times heavier than the Sun. Dominguez et al. (1999) studied the impact of axions on the AGB phase of low- and intermediate-mass stars (up to 9M_⊙), showing that they produce relevant changes in the evolution and, in particular, cause a deficit of massive C–O WDs. Later, Friedland et al. (2013) considered the impact of axions on the He-burning phase of solar-metallicity stars with mass in the 8–12 M_⊙ range, focusing, in particular, on the axion effects on the blue loops. In the present work, we are considering the impact of the axion emission on the evolution of massive stars and show how axions may reduce the pre-explosive luminosity of SN II progenitors and, in turn, increase the estimated mass.

The axion luminosity is given by

$$\begin{eqnarray}&&{L}_{a}={\int }_{0}^{R}4\pi {r}^{2}\rho {\epsilon }_{a}{dr},\end{eqnarray} \tag{ 6 }$$

where R is the stellar radius, ρ(r) and T(r) are the density and temperature internal profiles, and _a(ρ, T) is the rate of energy drained by axions per unit mass (see the Appendix). In the hot plasma of a star, _a takes contributions mostly from the following processes:

• 1.  
  Primakoff, i.e., photon–axion conversion in the electrostatic field of electrons and ions.
• 2.  
  Compton, i.e., scattering of photons on electrons or ions.
• 3.  
  Electron–positron annihilation e⁺e⁻ → γ + a, also known as pair production.
• 4.  
  Bremsstrahlung e + Ze → e + Ze + a on electrons and nucleons.

The numerical recipes we have used to estimate these rates are summarized in the Appendix, so the results can be reproduced and tested in future analyses.

The expected axion luminosity along the whole evolution of a 20 M_⊙ star is shown in Figure 9. In the right panel, we show the contributions to the total axion luminosity due to the four production processes described above. It is evident that the Primakoff and Compton contributions are the major sources of axions in these massive stars. In the left panel, we compare the total axion luminosity to the thermal neutrino luminosity. Axions clearly represent a significant energy sink in these stars, which is much larger than that induced by neutrinos for a considerable portion of the evolution. In particular, major effects are produced starting from the late part of the He-burning phase. As shown in Figure 10, when the axion energy loss is considered, the C burning and the Ne burning take place at higher temperature and density. Moreover, the envelope freeze-out occurs earlier than in standard models. Consequently, the final luminosity is substantially lower (see Figure 11). In addition, axions may also affect the compactness of pre-explosive models. We recall that the evolutionary tracks presented in this paper have been stopped when the maximum temperature attains ∼4 GK. The final mass distributions within the core of the two 20 M_⊙ evolutionary tracks, with and without axion energy loss, are compared in Figure 12. As expected, the axion model shows a higher mass concentration in the 0.006–0.03 R_⊙ region. However, the innermost portion of the core (M < 1.6 M_⊙) is not affected by axion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Then, we have computed two additional sets of stellar models. In the first one we have included the energy loss induced by the axion–photon coupling alone, namely, g_aγ = 0.6 × 10⁻¹⁰ GeV⁻¹ and g_ae = 0. In this case, only the Primakoff process is activated. This choice would be appropriate for, e.g., a pure hadronic axion model, such as the KSVZ, or an ALP coupled with photons only, with a coupling close to the astrophysical and experimental upper bound. In the second set of stellar models we have switched on the contribution from the axion–electron coupling by fixing g_aγ = 0.6 × 10⁻¹⁰ GeV⁻¹ and g_ae = 4 × 10⁻¹³, corresponding roughly to the limit from globular cluster stars. The initial mass–final luminosity relations for nonrotating and rotating massive-star models are reported in Figure 13.

Zoom In Zoom Out Reset image size

The initial masses corresponding to a final luminosity $\mathrm{log}L/{L}_{\odot }=5.1$ are 19.7 and 20.9 for the models with the axion–photon coupling only and the models with both the axion–photon and the axion–electron couplings, respectively.

Worthy of notice is the peculiar low final luminosity we found in the case of the 14 M_⊙ model, namely, $\mathrm{log}L/{L}_{\odot }=4.37$, of the second set of axion models. The plot of the corresponding evolutionary track illustrates the origin of this result (see Figure 14). Indeed, this model experiences an extended blue loop during the He-burning phase, but at variance with the standard evolutionary track (no axions), the model with both the axion–photon and axion–electron couplings does not complete this loop, returning on the RSG branch at the end of the He burning. Instead, the luminosity freeze-out occurs before the star closes its blue loop, so that the final point is substantially fainter and bluer than that obtained in the standard case. We found that this behavior affects the evolutionary tracks of models with mass in a restricted range around 14 M_⊙, while the 13 M_⊙ model, which also experiences a blue loop, shows normal behavior.

Zoom In Zoom Out Reset image size

We have shown that axions would substantially modify the initial mass–final luminosity relation of massive stars, thus moving to higher values the mass range of SN II progenitors. Axions also represent an appealing solution for various stellar anomalies, such as those related to the WD cooling rate. The considerable larger amount of WDs and Type II progenitors expected from the LSST deep photometric survey will definitely shed light on these problems (Bechtol et al. 2019; Drlica-Wagner et al. 2019). Moreover, the axion experimental proposals of the next generation could possibly test this hypothesis. In order to make our results clearer from an experimental point of view and to guide the experimental effort, here we make some general considerations and study the potential of proposed axion experiments.

As seen, the emission of axions shifts up the initial mass that corresponds to $\mathrm{log}L/{L}_{\odot }=5.1$, where L is the final stellar luminosity. The exact shift depends on several physical parameters. However, we have found a reasonable parameterization as the sum of terms proportional to the square of the coupling constants:¹⁴

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}M}{{M}_{\odot }}\approx 2.2{\left(\displaystyle \frac{{g}_{a\gamma }}{{10}^{-10}}\right)}^{2}+\displaystyle \frac{1.5}{16}{\left(\displaystyle \frac{{g}_{{ae}}}{{10}^{-13}}\right)}^{2}.\end{eqnarray} \tag{ 7 }$$

The parameter band corresponding to the mass shift 0.5 M_⊙ ≤ ΔM ≤ 3 M_⊙, as derived from the above equation, is shown as a gray hatched area in Figure 15. The figure also shows the experimental potential to explore the region with proposed experiments such as BabyIAXO and IAXO (Armengaud et al. 2014), ALPS II (Bähre et al. 2013), and DARWIN (Aalbers et al. 2016). In the case of CAST, BabyIAXO, IAXO, and ALPS II, the experimental potential assumes low-mass ALPS. The region hatched in purple in the figure is the one inferred from observations of low-mass stars (Giannotti et al. 2016b). The area overlaps, albeit in a small region, with the parameters required for a mass shift ΔM ≲ 1 M_⊙. Larger mass shifts require a larger electron coupling (larger photon couplings are already excluded by CAST, at least at mass below ∼20 meV). In any case, the parameter range is accessible to the next generation of axion helioscope experiments, as well as to ALPS II. On the other hand, as evident from the figure, WIMP detectors such as LUX, or even the future DARWIN detector, do not have the capability to probe ALP parameters that give a mass shift up to a few M_⊙.

Zoom In Zoom Out Reset image size

QCD axions deserve a separate analysis. Because of the mass-coupling relation, the axion couplings of interest for stellar evolution imply masses of ∼10–100 meV and could not be detected by the helioscope experiments, except for IAXO. Figure 16 shows the parameter space for the DFSZ axion, which may interact with both electrons and photons, together with the regions excluded by stellar evolution considerations. The main result is that an upward shift of about 2 M_⊙ of the initial masses of SN IIP progenitors is allowed when the RGB and HB bounds to the axion couplings are assumed. Moreover, IAXO would definitely be able to test the possibility that DFSZ axions are indeed responsible for the missing energy required in our problem.

Zoom In Zoom Out Reset image size

The Primakoff production of axions in a stellar core, that is, the production of axions from thermal photons conversion in the electrostatic field of the nuclei or of the electrons γ + Ze → γ + a, is widely discussed in the literature (e.g., Raffelt 1986; Raffelt & Dearborn 1987). Our method is based on Raffelt & Dearborn (1987) but includes a new parameterization of the degeneracy effects.

The axion emission rate is calculated summing the contributions from the scattering of photons on electrons and ions. We assume that ions are nondegenerate while electrons may have a certain degree of degeneracy. The Coulomb interaction between charges and thermal photons is screened. The typical screening scale is related to the longitudinal component of the polarization tensor, which, in the static and nonrelativistic limit, is (Raffelt 1996)

$$\begin{eqnarray}&&{\pi }_{L}=\displaystyle \frac{4{Z}^{2}\alpha \,m}{\pi }{\int }_{0}^{\infty }\displaystyle \frac{1}{\exp [\beta (E-\mu )]+1}\,{dp},\end{eqnarray} \tag{ 11 }$$

where β = 1/k_BT, with k_B the Boltzmann constant, μ is the chemical potential, and E = m + p²/2m is the fermion energy. In the limit of a nondegenerate plasma, π_L converges to

$$\begin{eqnarray}&&{\pi }_{L}\to {\kappa }_{{\rm{D}}}^{2}=\displaystyle \frac{4\pi {Z}^{2}\alpha \,n}{T},\end{eqnarray} \tag{ 12 }$$

where n is the number density of the charged particles, setting the Debye length as the typical screening scale. For a very degenerate plasma, on the other hand, the relevant screening parameter becomes the Thomas–Fermi scale

$$\begin{eqnarray}&&{\pi }_{L}\to {\kappa }_{\mathrm{TF}}^{2}=\displaystyle \frac{4{Z}^{2}\alpha \,m\,{p}_{F}}{\pi }.\end{eqnarray} \tag{ 13 }$$

In general, to account for any degree of degeneracy, it is convenient to introduce the degeneracy parameter R_deg, defined such that

$$\begin{eqnarray}&&{\kappa }^{2}={R}_{\deg }{\kappa }_{{\rm{D}}}^{2}.\end{eqnarray} \tag{ 14 }$$

We have derived a good numerical fit for the degeneracy parameter as a function of the temperature and density of the electron gas. Introducing the dimensionless parameter

$$\begin{eqnarray}&&z=\displaystyle \frac{\rho }{{T}^{3/2}{\mu }_{e}},\end{eqnarray} \tag{ 15 }$$

where T is in Kelvin and ρ in g cm⁻³, we find

$$\begin{eqnarray}&&{R}_{\deg }=1,\end{eqnarray} \tag{ 16 }$$

for z ≤ 5.45 × 10⁻¹¹,

$$\begin{eqnarray}&&{R}_{\deg }=0.63+0.3\arctan \left(0.65-9316\,{z}^{0.48}+\displaystyle \frac{0.019}{{z}^{0.212}}\right),\end{eqnarray} \tag{ 17 }$$

for 5.45 × 10⁻¹¹ < z ≤ 7.2 × 10⁻⁸, and

$$\begin{eqnarray}&&{R}_{\deg }=4.78\times {10}^{-6}{z}^{-0.667},\end{eqnarray} \tag{ 18 }$$

for z > 7.2 × 10⁻⁸. The precision of the fit is always better than 10%, in the whole region of interest.¹⁵

A further effect of the partial degeneracy is the reduction of the effective number of targets n → n_eff, defined, for example, in Raffelt & Dearborn (1987). As shown in Payez et al. (2015), this correction factor is again R_deg. Therefore, we find n_eff = R_deg n.

We are now ready to provide a numerical recipe for the Primakoff axion production rate. Let us first define the quantities

$$\begin{eqnarray}&&{y}_{\mathrm{pl}}=\displaystyle \frac{{\omega }_{\mathrm{pl}}}{T}=\displaystyle \frac{{\left(\rho /{\mu }_{e}\right)}^{1/2}}{{\left(1+{\left(1.02\times {10}^{-6}\rho /{\mu }_{e}\right)}^{2/3}\right)}^{1/4}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{y}_{\mathrm{ions}}=2.57\times {10}^{10}{\left(\displaystyle \frac{\rho }{{T}^{3}}\displaystyle \sum _{\mathrm{ions}}\displaystyle \frac{{Z}_{j}^{2}{X}_{j}}{{A}_{j}}\right)}^{1/2},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{y}_{\mathrm{el}}=2.57\times {10}^{10}{R}_{\deg }{\left(\displaystyle \frac{\rho }{{T}^{3}}\displaystyle \sum _{\mathrm{ions}}\displaystyle \frac{{Z}_{j}{X}_{j}}{{A}_{j}}\right)}^{1/2},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{y}_{s}={\left({y}_{\mathrm{ions}}^{2}+{y}_{\mathrm{el}}^{2}\right)}^{1/2},\end{eqnarray} \tag{ 22 }$$

where T is in K and ρ in g cm⁻³, and the function

$$\begin{eqnarray}&&f({y}_{\mathrm{pl}},{y}_{s})=\displaystyle \frac{1}{4\pi }{\int }_{{y}_{\mathrm{pl}}}^{\infty }\displaystyle \frac{{y}^{2}\sqrt{{y}^{2}-{y}_{\mathrm{pl}}^{2}}}{{e}^{y}-1}\,I({y}_{\mathrm{pl}},{y}_{s}){dy},\end{eqnarray} \tag{ 23 }$$

with

$$\begin{eqnarray}&&\begin{array}{l}I=\displaystyle \frac{{r}^{2}-1}{s}\mathrm{ln}\left(\displaystyle \frac{r-1}{r+1}\right)+\displaystyle \frac{{\left(r+s\right)}^{2}-1}{s}\\ \mathrm{ln}\left(\displaystyle \frac{s+r+1}{s+r-1}\right)-2,\end{array}\end{eqnarray} \tag{ 24 }$$

and

$$\begin{eqnarray}&&r=\displaystyle \frac{2{y}^{2}-{y}_{\mathrm{pl}}^{2}}{2y\sqrt{{y}^{2}-{y}_{\mathrm{pl}}^{2}}},\qquad s=\displaystyle \frac{{y}_{s}^{2}}{2y\sqrt{{y}^{2}-{y}_{\mathrm{pl}}^{2}}}.\end{eqnarray} \tag{ 25 }$$

Therefore, the axion Primakoff emission rate reads

$$\begin{eqnarray}&&\varepsilon =4.71\times {10}^{-31}{g}_{10}^{2}{T}^{4}\left[\displaystyle \sum _{\mathrm{ions}}\left({Z}_{j}^{2}+{R}_{\deg }{Z}_{j}\right)\displaystyle \frac{{X}_{j}}{{A}_{j}}\right]f({y}_{\mathrm{pl}},{y}_{s})\displaystyle \frac{\mathrm{erg}}{{\rm{g}}\cdot {\rm{s}}},\end{eqnarray} \tag{ 26 }$$

with g₁₀ = g_aγ/10⁻¹⁰ GeV⁻¹ and T in K.

The Compton axion production, γ + e → γ + a, is the production of axions from the scattering of thermal photons on electrons, and it is driven by the axion–electron coupling. The nonrelativistic cross section is (Raffelt 1996)

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{1}{3}\alpha {\left(\displaystyle \frac{{g}_{{ae}}}{{m}_{e}}\right)}^{2}{\left(\displaystyle \frac{\omega }{{m}_{e}}\right)}^{2},\end{eqnarray} \tag{ 27 }$$

where ω is the photon energy, assumed to be ω ≪ m_e. Here we have neglected the plasma frequency. This simplification, however, does not substantially modify the result. Additionally, the Compton process is relevant only in the nondegenerate regime since in a degenerate plasma bremsstrahlung dominates.

The induced energy loss is

$$\begin{eqnarray}&&\varepsilon ={R}_{\deg }\displaystyle \frac{{n}_{e}}{\rho }\int \displaystyle \frac{2\,{d}^{3}k}{{\left(2\pi \right)}^{3}}\displaystyle \frac{\sigma \omega }{{e}^{\omega /T}-1},\end{eqnarray} \tag{ 28 }$$

where R_deg accounts for the reduction in number of effective electron targets when the plasma is degenerate.¹⁶ We already have a fit for θ_deg. Extracting the dependence on ω and setting x = ω/T, we find

$$\begin{eqnarray}\begin{array}{rcl}\varepsilon & = & {R}_{\deg }\displaystyle \frac{{n}_{e}}{\rho }\displaystyle \frac{8\pi }{{\left(2\pi \right)}^{3}}\left(\displaystyle \frac{\sigma }{{\omega }^{2}}\right){T}^{6}{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{x}^{5}}{{e}^{x}-1}\,{dx}\\ & \simeq & {R}_{\deg }\displaystyle \frac{122.1}{3{\pi }^{2}}\displaystyle \frac{1}{{m}_{e}^{4}\,{m}_{u}\,{\mu }_{e}}\,\alpha \,{g}_{{ae}}^{2}\,{T}^{6}.\end{array}\end{eqnarray} \tag{ 29 }$$

Numerically,

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{compton}}={\theta }_{\deg }2.66\times {10}^{-48}\,{g}_{13}^{2}\,\displaystyle \frac{{T}^{6}}{{\mu }_{e}},\end{eqnarray} \tag{ 30 }$$

where g₁₃ = g_ae/10⁻¹³.

The bremsstrahlung process, e + Ze → e + Ze + a, is the most important axion production mechanism in a degenerate plasma. Therefore, we will discuss the degenerate limit first. In this limit we neglect the contribution from the scattering on a single electron, since very few electron targets are available in this limit. Additionally, the Debye screening length is a good approximation for the screening length for ions. With these approximations, we find (in units of erg g⁻¹ s⁻¹)

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{BD}}=8.6\,F\times {10}^{-33}\,{g}_{13}^{2}{T}^{4}\left(\sum \displaystyle \frac{{X}_{j}{Z}_{j}^{2}}{{A}_{j}}\right),\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&F=\displaystyle \frac{2}{3}\mathrm{ln}\left(\displaystyle \frac{2+{\kappa }^{2}}{{\kappa }^{2}}\right)+\left[\displaystyle \frac{2+5{\kappa }^{2}}{15}\mathrm{ln}\left(\displaystyle \frac{2+{\kappa }^{2}}{{\kappa }^{2}}\right)-\displaystyle \frac{2}{3}\right]{\beta }_{F}^{2},\end{eqnarray} \tag{ 32 }$$

with

$$\begin{eqnarray}&&{\kappa }^{2}=\displaystyle \frac{{k}_{{\rm{D}}}^{2}}{2{p}_{F}^{2}}=9.27\times {10}^{4}\,\displaystyle \frac{{\rho }^{1/3}}{T}{\left(\sum \displaystyle \frac{{X}_{j}{Z}_{j}^{2}}{{A}_{j}}\right)}^{1/3},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\beta }_{F}=\displaystyle \frac{{p}_{F}}{{E}_{F}}=\displaystyle \frac{{p}_{F}}{\sqrt{{m}_{e}^{2}+{p}_{F}^{2}}},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{p}_{F}={\left(3{\pi }^{2}\displaystyle \frac{\rho }{{m}_{u}}\sum \displaystyle \frac{{X}_{j}{Z}_{j}}{{A}_{j}}\right)}^{1/3}=5155{\left(\rho \sum \displaystyle \frac{{X}_{j}{Z}_{j}}{{A}_{j}}\right)}^{1/3}\mathrm{eV}.\end{eqnarray} \tag{ 35 }$$

In the nondegenerate limit one finds (again, in units of erg g⁻¹ s⁻¹)

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{BND}}=4.7\times {10}^{-25}{g}_{13}^{2}{T}^{2.5}\,\displaystyle \frac{\rho }{{\mu }_{e}}\sum \displaystyle \frac{{X}_{j}\,{Z}_{j}}{{A}_{j}}\left({Z}_{j}+\displaystyle \frac{1}{\sqrt{2}}\right).\end{eqnarray} \tag{ 36 }$$

Equation (36) does not take into account the screening effects, which are always small in the limits of validity of this expression. To take those effects into account, one has to substitute the sum in Equation (36) with

$$\begin{eqnarray}&&\sum \displaystyle \frac{{X}_{j}}{{A}_{j}}\left[{Z}_{j}^{2}\left(1-\displaystyle \frac{5}{8}\displaystyle \frac{{k}_{S}^{2}}{{m}_{e}T}\right)+\displaystyle \frac{{Z}_{j}}{\sqrt{2}}\left(1-\displaystyle \frac{5}{4}\displaystyle \frac{{k}_{S}^{2}}{{m}_{e}T}\right)\right],\end{eqnarray} \tag{ 37 }$$

where k_S accounts for the screening that, in the nondegenerate limit, takes contributions from both electrons and ions. Numerically,

$$\begin{eqnarray}&&\displaystyle \frac{{k}_{S}^{2}}{{m}_{e}T}=1.12\times {10}^{11}\,\displaystyle \frac{\rho }{{T}^{2}}\,\sum \displaystyle \frac{{X}_{j}({Z}_{j}+{Z}_{j}^{2})}{{A}_{j}}.\end{eqnarray} \tag{ 38 }$$

For the intermediate degeneracy case, we use the approach discussed in Raffelt & Weiss (1995):

$$\begin{eqnarray}&&\varepsilon ={\left(\displaystyle \frac{1}{{\varepsilon }_{{\rm{D}}}}+\displaystyle \frac{1}{{\varepsilon }_{\mathrm{ND}}}\right)}^{-1}.\end{eqnarray} \tag{ 39 }$$

Finally, axions can be produced from the annihilation of an electron–positron pair, e⁺e⁻ → γ + a. According to Pantziris & Kang (1986), we define the dimensionless parameters

$$\begin{eqnarray}&&\lambda =T/{m}_{e}\simeq 1.69\times {10}^{-10}{T}_{K}\end{eqnarray} \tag{ 40 }$$

and

$$\begin{eqnarray}&&\nu =\mu /T\end{eqnarray} \tag{ 41 }$$

and distinguish among different plasma conditions.

In the nondegenerate and nonrelativistic regime, we find

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{{g}_{{ae}}^{2}\alpha }{4{\pi }^{3}\rho }{m}_{e}^{2}{T}^{3}{e}^{-2/\lambda }\left(1+\displaystyle \frac{\lambda }{2}\right).\end{eqnarray} \tag{ 42 }$$

The degenerate case is given in Pantziris & Kang (1986) and applies to the case of λ ≪ 1 and 1/λ ≪ ν ≪ 2/λ. The result is cumbersome,

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{{g}_{{ae}}^{2}\alpha }{{\pi }^{4}\rho }{m}_{e}^{2}{T}^{3}{e}^{\eta }{\int }_{0}^{\infty }\displaystyle \frac{{dx}\sqrt{x}}{{e}^{x-a}+1}{\int }_{0}^{\infty }\displaystyle \frac{{dy}\sqrt{y}}{{e}^{y}}\left(1+\displaystyle \frac{x+y}{6}\lambda \right),\end{eqnarray} \tag{ 43 }$$

where η = ν − 1/λ.¹⁷ A simplified expression, accurate to about 15% in the whole range of interest, is

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{{g}_{{ae}}^{2}\alpha }{{\pi }^{4}\rho }{m}_{e}^{2}{T}^{3}{e}^{\nu -1/\lambda }f(\nu -1/\lambda ),\end{eqnarray} \tag{ 44 }$$

where

$$\begin{eqnarray}&&f(x)=0.605\,{e}^{-0.84x}.\end{eqnarray} \tag{ 45 }$$

Finally, in the nonrelativistic and mildly degenerate case (μ ≃ m_e), we find the expression

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{{g}_{{ae}}^{2}\alpha }{4{\pi }^{3}\rho }{m}_{e}^{2}{T}^{3}{e}^{-2/\lambda }0.76\left(1+0.53\lambda \right).\end{eqnarray} \tag{ 46 }$$

A numerical expression for the axion pair production in cgs units that summarizes all the previous results is

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{2.14\times {10}^{16}{g}_{13}^{2}{\lambda }^{3}}{\rho }\,\displaystyle \frac{\mathrm{erg}}{{\rm{g}}\,{\rm{s}}}\,{e}^{-\mathrm{Max}(\eta ,0)-2/\lambda }F,\end{eqnarray} \tag{ 47 }$$

where F is

• 1.  
  nondegenerate:

  $$\begin{eqnarray}&&F=\displaystyle \frac{1}{4}\left(1+\displaystyle \frac{\lambda }{2}\right);\end{eqnarray} \tag{ 48 }$$
• 2.  
  mildly degenerate:

  $$\begin{eqnarray}&&F=0.19\left(1+0.53\lambda \right);\end{eqnarray} \tag{ 49 }$$
• 3.  
  degenerate:

  $$\begin{eqnarray}&&F=\displaystyle \frac{1}{\pi }f(\eta ).\end{eqnarray} \tag{ 50 }$$