Progenitor Dependence of Neutrino-driven Supernova Explosions
with the Aid of Heavy Axion-like Particles

In this work, we consider a photophilic ALP model where the ALPs are generated through two photon interaction processes; the Primakoff process $(\gamma+p\to a+p)$ and the photon coalescence $(\gamma+\gamma\to a)$.

The Primakoff rate is given as

where $n_{a}$ is the number density of ALPs, $\omega$ is photon energy and $f(\omega)$ is the Bose-Einstein distribution of photons, $\omega_{\rm{pl}}$ is plasma frequency. The ALP energy $E$ is equal to $\omega$ because of the energy conservation, and $\Gamma_{\gamma\to a}$ is given by Di Lella et al. (2000)

where $T$ is the temperature, $\kappa$ is the Debye-Hückel scale, $p$ is the ALP momentum, and $k$ is the wave number of photons in plasma.

The photon coalescence rate is given as Di Lella et al. (2000)

The energy loss rates via these two processes $Q_{\rm{cool}}$ are given by

The photon coalescence contributes to the ALP production only when ALPs are heavier than $2\omega_{\rm{pl}}$ where the plasma frequency $\omega_{\rm{pl}}$ is the “effective photon mass”.

The ALPs produced by these processes propagate through the SN matter and affect the energy transfer. If they decay into photons within the SN matter and the photons are absorbed by the matter in the post-shock region, they could assist the shock revival. We consider the inverse Primakoff process $(a\to\gamma)$ and radiative decay $(a\to\gamma\gamma)$ as the ALP decay processes.

The inverse Primakoff rate is given as Lucente et al. (2020)

where $\beta_{E}$ = $v_{a}/c$ and $v_{a}$ is the velocity of ALPs. The radiative decay rate is given as

The previous study Mori et al. (2022) incorporated the ALP transport into SN simulations as follows. The evolution of the ALP energy per unit volume $\mathcal{E}$ is described by

where $\mathcal{F}$ is the ALP energy flux and $Q_{\textrm{heat}}$ is the heating rate per unit volume due to ALPs. Assuming stationarity and spherical symmetry, Eq. (7) is simplified to

here $L$ is ALP luminosity. We solve Eq. (8) at each time step to obtain $Q_{\mathrm{heat}}$. We then update the internal energy at the $n$-th step as

where $e_{\mathrm{int}}$ is the internal energy and $\Delta t$ is the time step.

The numerical setup in this study is essentially the same as in Mori et al. (2022). We employ the 3DnSNe code Takiwaki et al. (2016), which is a multi-dimensional neutrino radiation hydrodynamics code developed to study core-collapse SNe. The neutrino transport is solved by the three-flavor isotropic diffusion source approximation (IDSA) scheme Liebendoerfer et al. (2009). We use the state-of-the-art neutrino opacity Kotake et al. (2018) and the neutrino energy spectrum is discretized with 20 energy bins for $0<\varepsilon_{\nu}\leq 300\,\rm{MeV}$. We take account of the effective general relativistic effect Marek et al. (2006) for the gravitational potential and the gravitational redshift for the neutrino transport.

We use the 11.2$\,M_{\odot}$ progenitor model from Woosley et al. (2002), and the 20$\,M_{\odot}$ and 25$\,M_{\odot}$ models from Woosley and Heger (2007), to investigate the progenitor dependence of ALP effects on core-collapse SN explosions. There are two reasons for the choice of these progenitors. First, these three progenitor models are red supergiants and massive enough for their core to undergo gravitational collapse during the final stage of their evolution Smartt et al. (2009, 2004). Their explosions will appear as type II SNe, the most commonly observed type among core-collapse SNe. Second, these progenitor models have a different mass-accretion history, due to their differences in compactness, as shown in Fig. 1. Since the gravitational potential of the accreting gas is the dominant energy reservoir for neutrino-driven explosions, the difference in the mass accretion rate leads to different dynamical evolution. Moreover, the thermal evolution of their core (for example, the maximum temperature achieved in the core) strongly affects the ALP production rate. The models examined in this study are suitable for the purpose to inspect the dependence of the ALP emissions on the progenitor structure.

For these progenitor models, we performed one-dimensional core-collapse SN simulations to investigate the impact of ALPs on the core-collapse SN dynamics. The ALP parameter space has been constrained by various studies (e.g. Lucente et al., 2020, 2022; Jaeckel et al., 2018b; NA64 Collaboration et al., 2020; Carenza et al., 2020; Döbrich et al., 2020; Dolan et al., 2017; Diamond et al., 2024; Dev et al., 2024). ALPs with $m_{a}=\mathcal{O}(100)$ MeV are reported to have significant effects on core-collapse SN dynamics Caputo et al. (2022b); Mori et al. (2022). Considering these results, we have chosen the parameter space of ALP mass, $m_{a}=100-800\,\rm{MeV}$, and ALP-photon coupling constant, $g_{\rm{10}}=g_{a\gamma}/(10^{-10}\,\rm{GeV^{-1}})=4-10$. Note that there is also discussion that excludes some parts of this range Caputo et al. (2022a); Fiorillo et al. (2025); Diamond et al. (2024).

The models are labeled as follows. For example, a model with the $20.0\,M_{\odot}$ progenitor and ALPs with $m_{a}=200\,\rm{MeV}$ and $g_{\rm{10}}=8$ is called ‘s20_m2_g08’. In addition, the model without ALPs is labeled as ‘s20_m0_g00’ and referred to as the no-ALP model.

When ALPs are not taken into account, shock revival does not occur in one-dimensional core-collapse SN models (e.g. Liebendörfer et al., 2001; O’Connor et al., 2018), except for models with a low-mass progenitor Kitaura et al. (2006); Fischer et al. (2010); Hüdepohl et al. (2010). However, in some of our models considering ALPs, the shock wave is successfully revived and leads to a supernova explosion.

Figure 2 shows the time evolution of the shock radius for the s20_m2 model series with different coupling constants $g_{10}=4$, 6, 8 and 10, along with the no-ALP model denoted by s20_m0_g00. For the no-ALP model and the $g_{10}=4$ model (solid black and dash-dot blue lines), the shock wave stalls at $r\sim 150\,\rm{km}$ and fails to revive. This behavior is similar to conventional one-dimensional simulations without ALPs. However, in the cases of the $g_{10}=8$ and 10 models (dotted magenta and solid red lines), the shock wave successfully revives and turns to a runaway expansion around $t_{\rm{pb}}=0.3$ s. For the $g_{10}=6$ model (dashed green line), the shock wave gradually expands after $t_{\rm pb}\sim 0.3$ s, but the shock radius remains at 100 km at the end of the simulation and it is not clear if the shock expansion continues afterward. These models indicate the trend that the shock wave is more likely to revive for the models with higher coupling constants.

This trend can be explained by the dependence of the ALP heating rate $Q_{\mathrm{\rm{heat}}}$ induced by ALP decay photons on the coupling constant. In Fig. 3, the radial profiles of $Q_{\rm{heat}}$ are shown for the s20_m2 model series at the time $t_{\rm{pb}}=0.17$ s. We can see that a higher coupling constant induces higher ALP heating rates, which can lead to the successful shock revival.

Figure 3 shows that the ALP heating rates peak at $r\sim 10$ km regardless of the coupling constant. This is because the temperature distribution (drawn in the solid gray line) peaks at this location, where the ALP production rate is also maximized due to its sensitivity to temperature. Furthermore, because of the ray-by-ray approximation used in solving ALP transport, ALPs cannot propagate inward from the radius where they are produced. As a result, there is a precipitous cutoff in $Q_{\rm{heat}}$ at $r\sim 10$ km. Outside of the peak, $Q_{\rm{heat}}$ decreases following a $1/r^{2}$ dependence shown in Eq. (15).

We can also see that the neutrino heating rate (drawn in the solid black line) is dominant between the gain radius $R_{\rm{gain}}\sim 60$ km and the the shock radius $R_{\rm{shock}}\sim 85$ km. The neutrino cooling dominates over the heating inside the gain radius, and they balance out at $r=R_{\rm{gain}}$. The neutrino heating rate drops at $r=R_{\rm{shock}}$ because nucleons, which are the primary targets of the neutrino absorption process, are scarcely present outside the shock radius. The excess of the neutrino heating rate over the heating rate of ALPs implies that supernova explosions are primarily driven by the neutrino heating and ALPs play a secondary role.

In Fig. 4, we present the neutrino luminosity and mean energy for the s20_m2 model series, along with the residuals relative to the no-ALP model. Before shock revival ($t_{\rm{pb}}=0.2-0.3$ s), the luminosity and the mean energy for all flavors are essentially unchanged by considering the effects of ALPs. However, the models with higher coupling constants show decrease in the $\nu_{e}$ and $\bar{\nu}_{e}$ luminosities and their mean energies, because these models experience early shock revival, leading to a decrease in the mass accretion rate. Only a slight difference is observed in the $\nu_{x}$ mean energy among these models because $\nu_{x}$ is predominantly produced in the dense core via neutral current pair production (e.g. Betranhandy and O’Connor, 2020) and less affected by the mass accretion history. Although ALP emission can rapidly extract energy from the PNS, it takes a time on the order of the neutrino diffusion time (several seconds) for the information to reach the neutrino sphere Bar et al. (2020).

The SN dynamics is also dependent on the ALP mass. In Fig. 5, we plot the shock evolution for the s20_g08 model series with $m_{a}=100$, 200, 400, and 600 MeV. For the model with $m_{a}=100$ MeV (dash-dot blue line), the shock wave shows almost no change compared to the no-ALP model and fails to revive, whereas the shock waves of the models with $m_{a}=$ 200 and 400 MeV models (dashed green and dotted magenta lines) turn to expand at $\sim 300$ ms and 250 ms after bounce. Therefore, within the range of $m_{a}=100-400\,\rm{MeV}$, higher ALP masses lead to an earlier shock revival. However, in the case of the model with $m_{a}=600$ MeV (solid red line), the propagation of the bounce shock is slower than in the models with lighter ALPs, and in the $m_{a}=800$ MeV model, the shock wave does not initiate a runaway expansion, similar to the no-ALP model (solid black line). These results indicate that a higher ALP mass does not necessarily facilitate the shock revival and that the relationship between the ALP mass and the shock wave behavior exhibits a non-monotonic trend.

The non-monotonic dependence on $m_{a}$ can be explained by the $m_{a}$ dependence of the ALP production rate. Fig. 6 shows the radial profile of $Q_{\mathrm{cool}}$ at $t_{\rm{pb}}=0.17$ s. The cooling rates reach a peak at $r\sim 10\,\rm{km}$ regardless of the ALP masses, because the ALP production processes are sensitive to the matter temperature, which also peaks in this region. In the case of the s20_m6_g08 model (and the s20_m8_g08 model as well), the peak value is about an order of magnitude lower than the value in the other models. This is due to the Boltzmann suppression, as the temperature is not high enough to produce such heavy ALPs.

Figure 7 shows the time evolution of the ALP heating rate in the gain region defined as

For the s20_m1_g08 model (dash-dot blue line), the ALP heating rate is as low as $L_{\rm{heat}}\sim 10^{48}$ erg s^-1, which is insufficient for the revival of the shock wave. The s20_m2_g08 and s20_m4_g08 models (dashed green dotted magenta lines) show higher heating rates of $L_{\rm{heat}}\sim 10^{50}$ erg s^-1 at $t_{\rm{pb}}=0.1$ s compared to the other models, thus photons from decaying ALPs effectively heating the region behind the shock, promoting shock revival. The differences in the heating rates among these models can be explained by the fact that heavier ALP masses contribute to the higher heating rate and more efficient heating in the gain region because of their shorter mean free path. However, in the s20_m6_g08 model (solid red line), despite the higher mass, the heating rate $L_{\rm{heat}}$ is lower than that of the s20_m2_g08 and s20_m4_g08 models due to a reduction in the production rate. Furthermore, the s20_m6_g08 model shows a delayed increase in the heating rate compared to the other models, because it takes longer to reach a temperature to sufficiently produce 600 MeV ALPs. These delays in the increases in heating rate and temperature lead to the delayed shock expansion.

In this section, we discuss the progenitor dependence of the shock evolution and its underlying causes, focusing on the s11_g08, s20_g08 and s25_g08 model series with ALP masses $m_{a}=200$ and $600\,\rm{MeV}$.

Figure 8 shows the shock evolution for each model. For the $11.2M_{\odot}$ progenitor model, a bump in the shock radius is observed at $t_{\rm{pb}}\sim 0.1$ s, which is induced by the infall of the oxygen-silicon layer. For the models with $m_{a}=200$ MeV (the left panel of Figure 8), the $11.2M_{\odot}$ model (solid green line) shows slow outward propagation of the shock wave, although it is not clear whether the model successfully explodes. The $20M_{\odot}$ and $25M_{\odot}$ models (dotted magenta and dash-dot red lines) show a successful shock revival at $t_{\rm{pb}}\sim 0.3$ s.

Figure 8 also shows that the shock expansion is suppressed for the models with heavy ALPs ($m_{a}=600$ MeV, the right panel). This is because the production rate of heavy ALPs is relatively low. In particular, the s20_m6_g08 model (loosely dashed line) exhibits a pronounced delay in shock expansion, with the shock failing to reach $200\,\rm{km}$ by $t_{\rm{pb}}\sim 0.52$ s. These results show that the heavier progenitor models facilitated by ALPs are more likely to undergo shock expansion. This trend becomes more pronounced for the models with heavier ALPs such as $600$ MeV.

The dependence of the shock behavior on the progenitors is induced by the difference in the maximum temperature of the SN models. Fig. 9 shows the maximum temperature in the SN core, which is achieved at $r\sim 10$ km, with a focus on the s20 and s25 model series. For all models, the temperature increases over time and reaches $55-60\,\rm{MeV}$ at $t_{\rm{pb}}=0.3$ s due to the gravitational compression of the accretion layer associated with the growth of the PNS mass. After $t_{\rm{pb}}=0.3$ s, the models taking account of the ALP effects show that the increase in temperature is suppressed because the shock expansion achieved at this time causes a decrease in the accretion rate.

Figure 9 also shows that, for any ALP parameter sets, the peak temperature of the s25 series models (red lines) is higher than that of the s20 series models (magenta lines). The difference in the PNS temperature can be interpreted in terms of the positive correlation between the iron core mass and the compactness parameter O’Connor and Ott (2011). For massive or high-compactness progenitors, the iron core masses tend to be heavier due to the thermal pressure support. This explains a general tendency that a progenitor with higher compactness leads to a hotter PNS.

The ALP production rate $Q_{\mathrm{cool}}$ is higher for a heavier progenitor model because of the higher temperature. As a result, the ALP heating rate $Q_{\mathrm{heat}}$ becomes also higher for a heavier progenitor.

In Fig. 10, we show the evolution of the ALP cooling and heating rates in the gain region. Here the cooling rate is defined as

For the m2 model series, the cooling rates begin to increase early, with little noticeable difference between the $20\,M_{\odot}$ and $25\,M_{\odot}$ models (dotted magenta and dash-dot red lines). This is because the temperature required to produce 200 MeV ALPs is comparatively lower than that for heavier ALPs and realized shortly after post-bounce for both progenitor models. Accordingly, the differences in the ALP heating rates are observed to be negligible before the shock revival ($t_{\rm{pb}}\sim 0.3$ s). After the shock revival, the ALP heating rates show a substantial increase due to the extension of the gain region.

On the other hand, for the m6 model series, the cooling rates increase with a delay and are lower compared to the m2 model series. The reason for this is that it takes a longer time to realize a sufficient temperature to produce 600 MeV ALPs. We can also see that the s25_m6_g08 model (dash-dot-dot red line) exhibits the higher cooling and heating rates than the s20_m6_g08 model (loosely dashed magenta line), because of its higher temperature. From our systematic investigation, it is suggested that since $Q_{\mathrm{cool}}$ is sensitive to the temperature, $L_{\mathrm{cool}}$ and $L_{\mathrm{heat}}$ tend to be higher in the heavier star. As a result, in the s25_m6_g08 model, ALPs heat the gain region more effectively, making shock wave revival more likely than in the s20_m6_g08 model.

In Ref. Lucente et al. (2020), post-processing calculations of the ALP cooling rate were performed using one-dimensional simulations that artificially amplify neutrino heating and trigger the explosion. They showed that, at $t_{\rm{pb}}=1$ s for the $18\,M_{\odot}$ model, the cooling rate induced by ALPs with $m_{a}=600\,\rm{MeV}$ is $2-3$ orders of magnitude lower compared to the case for $m_{a}=200\,\rm{MeV}$. However, in our study, although the progenitor model and time differ — $20\,M_{\odot}$ and $t_{\rm{pb}}=0.5$ s, respectively — the difference in the cooling rates is only $\sim 0.5$ orders of magnitude. One of the reasons for this discrepancy is that, in our model, the delay in shock expansion for $m_{a}=600\,\rm{MeV}$ suppresses the decrease in mass accretion rate. This allows the temperature to continue increasing up to $6-8\,\rm{MeV}$ higher compared to the case for $m_{a}=200\,\rm{MeV}$ at $t_{\rm{pb}}=0.5$ s. As a result, ALPs with $m_{a}=600\,\rm{MeV}$ are more likely to be produced in our self-consistent simulations.

Figure 11 shows the explodability and explosion energies of our SN models with various ALP parameters. The symbols denote the final fate of the models and the color filling each symbol shows the explosion energy $E_{\rm{exp}}$ at the time when the shock wave reaches $r=500\,\rm{km}$. In this plot, the crosses represent the models that fail to explode, while the filled circles represent the models with successful explosion. The triangles represent models in which the shock wave propagates outward but does not reach $r=500\,\rm{km}$ at the end of the simulations. This figure indicates that the heavier progenitor models have a more extended ALP parameter region that can lead to shock revival even with the spherically symmetric geometry. It is also shown that the heavier progenitor models exhibit higher explosion energies.