Axion Instability Supernovae

New particles coupled to the Standard Model can equilibrate in stellar cores if they are sufficiently heavy and strongly coupled. In this work, we investigate the astrophysical consequences of such a scenario for massive stars by incorporating new contributions to the equation of state into a state of the art stellar structure code. We focus on axions in the"cosmological triangle", a region of parameter space with $300{\rm\,keV} \lesssim m_a \lesssim 2$ MeV, $g_{a\gamma\gamma}\sim 10^{-5}$ GeV$^{-1}$ that is not presently excluded by other considerations. We find that for axion masses $m_a \sim m_e $, axion production in the core drives a new stellar instability that results in explosive nuclear burning that either drives a series of mass-shedding pulsations or completely disrupts the star resulting in a new type of optical transient -- an \textit{Axion Instability Supernova}. We predict that the upper black hole mass gap would be located at $37{\rm M}_\odot \le M\le 107{\rm M}_\odot$ in these theories, a large shift down from the standard prediction, which is disfavored by the detection of the mass gap in the LIGO/Virgo/KAGRA GWTC-2 gravitational wave catalog beginning at $46_{-6}^{+17}{\rm M}_\odot$. Furthermore, axion-instability supernovae are more common than pair-instability supernovae, making them excellent candidate targets for JWST. The methods presented in this work can be used to investigate the astrophysical consequences of any theory of new physics that contains heavy bosonic particles of arbitrary spin. We provide the tools to facilitate such studies.


I. INTRODUCTION
The observation of gravitational waves from merging binary black holes by the LIGO/Virgo/KAGRA (LVK) collaboration has enabled the study of the population statistics of astrophysical black holes [1]. The black hole mass function (BHMF) is predominantly shaped by the processes governing the structure, evolution, and fate of massive stars, making it a sensitive probe of stellar structure theory [2,3]. It is also a novel and powerful probe of new physics beyond the Standard Model. The evolution of massive stars can be altered by modifications of gravity that arise in leading dark energy theories [4], additional energy losses from new light particles [5][6][7], and energy injections from dark matter self-annihilation [8,9], all of which leave imprints on the BHMF that can be distinguished with gravitational wave data from LVK [10]. New physics process can affect both the shape of the BHMF and the locations of features such as the upper black hole mass gap (BHMG), which allows competing theories to be distinguished.
The BHMG is a feature in the BHMF that refers to the absence of astrophysical black holes with masses in the range 50M M 120M . The physical origin of this feature is the so-called pair-instability, a stellar instability due to electron-positron pair-production in the cores of massive stars, which have densities and temperatures that are sufficient to produce copious amounts of nonrelativistic e − e + pairs. Small contractions of the star result in temperature and density increases; a corresponding increase of pressure would be expected to counteract the contraction, but leads instead to an increase in the rate of e − e + production [11]. This causes a runaway contraction that is only halted by explosive oxygen ignition after the core temperature and density have become sufficiently high. The resulting explosion unbinds the star -a process referred to as a pair-instability supernova (PISN) -leaving no black hole (BH) remnant. A recent analysis [10] of the LVK GWTC-2 catalog has determined the location of the lower edge of the BHMG to be M BHMG = 46 +17 −6 M (M BHMG = 54 +6 −6 M if GW190521, a potential outlier, is excluded from the analysis).
In a recent publication [6], we identified the existence of a new stellar instability triggered by heavy new particles. If the new particle couples to the Standard Model (SM) with sufficient strength then it can be produced in the cores of massive stars through SM interactions. Heavy particles are especially interesting because they cannot be produced in well-studied stars with lower masses. If the particles are subsequently unable to free-stream out of the star, they can come into equilibrium with the stellar plasma. Rather than acting as a new source of energy loss, as is the case with light particles, heavy particles alter the equation of state (EOS) of the thermal plasma [11]. This process is analogous to the e − e + pairinstability. The heavy new particles are produced with non-relativistic velocities, and therefore rob the star of pressure support and act to destabilize the star, resulting in a similar runaway contraction. In [6], we demonstrated the existence of such an instability by calculating the equation of state of a gas of coupled ions, photons, electrons, positrons, and heavy new particles. This is sufficient to determine the existence of the instability, but calculating its effect on the structure, evolution, and fate of massive objects requires detailed stellar evolution simulations.
In this paper, we present the results of incorporating heavy new particles into the EOS of stellar matter on the fates of massive stars. We have modified the stellar structure code MESA [12][13][14][15][16] to include an additional sector of heavy particles in the equation of state. We focus on axions with masses m a ∼ 500 keV coupled to photons with coupling strength g aγγ ∼ 10 −5 GeV −1 . This corresponds to the "cosmological triangle" -a small region in the {g aγγ , m a } parameter space that is not currently excluded by other probes of axions [17,18]. We find that stars with initial mass M 48M undergo explosive oxygen burning in the star's core. The explosion is so violent that it unbinds the entire star, leaving no BH remnant. Thus, we predict the existence of a new optical transient: an Axion Instability Supernova (AISN). In stars with initial masses in the range 38M M 48M the instability drives a series of mass-shedding pulsations similar to the pulsational pair instability supernova (PPISN) process driven by the pair-instability. We therefore refer to this process and the associated novel optical transient as a Pulsational Axion Instability Supernova (PAISN). Stars with mass M 38M avoid the instability and core-collapse to form BHs with masses nearly identical to the initial mass. The upper black hole mass gap would therefore be located at M BHMG ≈ 37M , in stark constant to the SM prediction M BHMG ≈ 50M [3,10]. We also find that the the upper edge of the BHMG is located at 107M , significantly lower than the SM prediction of 133M , making it more likely that LVK can detect BHs on the far side of the gap. The results above are robust to variations in metallicity and the 12 C(α, γ) 16 O rate. The detection of the upper BHMG feature in the LVK GWTC-2 gravitational wave catalog beginning at M BHMG = 46 +17 −6 M [10] therefore disfavors this axion mass. This paper is organized as follows. In section II we discuss the cosmological triangle for axions coupled to photons and demonstrate that particles in this region can equilibriate in stellar cores. In section III we describe the stellar structure code used for our simulations, derive the equation of state for heavy new particles, and explain its implementation into MESA. In section IV we present our results. We conclude in section V.

II. EQUILIBRATION IN THE COSMOLOGICAL TRIANGLE
The methods presented in this work can be applied to any bosonic particle (of arbitrary multiplicity) which equilibrates in post-helium burning stellar cores such as axions, chameleons, symmetrons, dilatons, or hidden photons. An important application, to which we devote our analysis in the rest of this paper, is axion-like particles (henceforward axions) in the so-called cosmological triangle (see, e.g., [17,18]): an open region in the {g aγγ , m a } parameter space bordered by constraints from beam-dump experiments, supernova cooling constraints, and stellar bounds. The cosmological triangle received its name from its shape, and the fact that cosmologydependent constraints can be deduced given additional assumptions on exotic contributions to ∆N eff [19]. This region of parameter space has recently been criticized based on limiting the energy deposition in the mantle and the outer envelopes of low-energy supernovae [20,21]. These limits are suggestive that new optical phenomena are possible in this range of parameter space; however, until observational signatures of this type of energy deposition are better characterized, it is important to pursue additional probes of this parameter space. In this subsection we show that axions in the cosmological triangle will indeed equilibrate on timescales relevant to massive star evolution, and in subsequent sections we work out additional observational clues to the presence of this new kind of particle.
There are two process relevant for heavy axion production in the mass and coupling range considered in this work -photon conversion (the Primakoff process) and photon coalescence [17,22,23]. The production rate for axions by the process of photon conversion is [18,22] where k = | k| is the photon momentum, ω its energy, p = | p| is the axion momentum, T is the plasma temperature, and k s describes plasma screening effects: where Y e is the number of electrons per baryon, Y j is the number per baryon of nuclear species with charge Z j , and m u is the atomic mass unit. The production rate for axions by the process of photon coalescence is (see e.g. [24][25][26]) where the axion decay rate is given by and where the latter is the plasma frequency (n e is the electron number density).
The equilibration time can be estimated as in [6] as τ eq ( Γ a→γ ) −1 , where the sum runs over all processes that allow the absorption of axions, since the rate of absorptive processes must equal the rate of production when in equilibrium. The rate of core temperature evolution in our simulation can be estimated as Γ evol =Ṫ /T where the dot indicates a derivative with respect to time. We can safely assume an equilibrium population of axions if Γ γ→a > Γ evol before the onset of the PISN or PPISN. Given that the decay rate alone, estimated as Γ decay ∼ g 2 a m 3 a , is of order microseconds for the mass and coupling considered in this work, this condition is satisfied in the cores of the stars that we simulate, except possibly during the pulsations. We therefore assume an equilibrium distribution during quasi-static stellar burning and make the conservative choice to disable our modifications to the EOS at the onset of the pulsational stage of evolution.

III. STELLAR MODELLING
Our simulations were performed using the publicly available stellar structure code MESA version 12778 [12][13][14][15][16]. Our modifications to MESA as well as a Mathematica code to produce equations of state for user-supplied masses and spins are publicly available at the following URL https://zenodo.org/record/6347632 [27].

A. Input Physics
MESA is a one-dimensional (meaning that it assumes spherical symmetry) hydrostatic code. It is equipped with a Harten-Lax-vanLeer-Contact (HLLC) hydrodynamical solver that can be switched on to simulate departures from hydrostatic equilibrium such as pulsations and shocks [16]. Our simulation of the pulsations and stellar explosions due to new particles is identical to the implementation described in [2,4,6,16,28].
Relevant prescriptions for our simulations are as follows. We treat convection according to the Cox prescription for mixing length theory [29] with efficiency parameter α MLT = 2.0. Semi-convection is modeled according to [30] with efficiency parameter α SC = 1.0. Convective overshooting is described using an exponential profile parameterized by f 0 , which sets the point inside the convective boundary where overshooting begins, and f ov , which determines the scale height of the overshoot. We set f 0 = 0.005 and f ov = 0.01. Our prescription for mass loss due to stellar winds follows that of [31] with the wind efficiency parameter (clumping parameter) η fixed to 0.1. Finally, we use the MESA default nuclear burning rates (these are a mixture of the NACRE [32] and REACLIB [33] tables) with the one exception of the 12 C(α, γ) 16 O reaction. Previous studies have found this to be the most important rate for determining the location of the upper BHMG [2,3,34] so we use the most recent results from [35]. In particular, we use their median values (changing this by ±1σ changes the location of the SM BHMG by ±2.5M [3,10]).
We derive the final black hole mass by simulating helium cores from the zero age helium branch (ZAHB) until either core collapse or PISN/AISN. The black hole mass is calculated as the mass of the layers expanding with velocity less than the escape velocity. The inlists given in the reproduction package for this paper are identical to those used to generate our results.

B. Equation of State for Heavy Bosonic Particles
We assume that the new particle is bosonic with g = 2s + 1 (with s the spin) degrees of freedom and further assume that it has no gauge quantum numbers, implying that it has zero chemical potential. This assumption can easily be relaxed in more complicated models. The contribution to the EOS can then be calculated by integrating over the distribution. Here we assume a thermal equilibrium Bose-Einstein distribution. We leave a detailed study of the thermalization processes for future work. Our procedure for calculating the EOS for heavy new particles and its implementation into MESA is described in Appendix A.
We assume that all of the axions that are produced by the process described in Eq. 1 remain trapped in the stellar core. This is only an approximation: some of the axions produced from photon conversion will be generated with velocities exceeding the escape velocity from the star. However, for the parameter values we choose, such particles should be able to scatter with SM particles on their way out of the star, and in doing so will on average lose energy to the gravitationally bound stellar material, which will lead them to eventually become trapped and populate the thermal phase space as assumed here. The process by which they scatter and become bound could lead to transport of internal stellar energy (see e.g. [23,[36][37][38][39]). The energy transport and subsequent redistribution in the core reduces the temperature gradient and hence the amount of convection. Convective mixing of 12 C inhibits the explosive burning of Oxygen [3] so additional energy transport by axions is likely to result in stronger pulsations and the onset of instability at lower ZAHB mass. Full numerical studies are needed to confirm this. For this reason, we believe that a full, out-of-equilibrium simulation of the dynamics of the axion population in the star will be an interesting topic to study in greater detail in future work. Nevertheless, our study is of value for the limiting cases of very tightly coupled new particles or modifications of the EOS induced by modified gravity or new couplings between SM particles.

IV. RESULTS
We report here on the impact of an axion with mass m a = m e = 511 keV -near the lower edge of the cosmological triangle. The effects of less-massive axions will be even more dramatically visible than an electron-mass axion. Intriguingly, we have verified that if the axion has a mass m a 2m e then stars will encounter the PISN before the AISN, thus recovering Standard Model-like behavior anticipated with no axion at all. Below, we comment on the significance of these relations for LVK observations.

A. Physics of the Axion Instability
The inclusion of axions with masses equal to electrons exacerbates pair-instability. The temperatures and densities where the axion-instability manifests are shown in Figure 1. This figure demonstrates that adding a bosonic particle with the same mass as electrons implies an extension of the region with Γ 1 = 4/3, where the stellar core is unstable. We therefore expect lower mass objects to experience mass-shedding pulsations (PAISN) and complete disruption (AISN) compared with the SM prediction. We also expect that these events will be more violent than the PPISN and PISN predicted by the SM. For m a 2m e , axion production is efficient under similar conditions as e − e + production, such that the effects of the pair-instability and axion-instability compound, resulting in a smaller value of Γ 1 at fixed temperature and density and therefore a stronger inward contraction when the instability is first encountered. This stronger contraction can only be countered by a more violent oxygen explosion.
A second consequence of the compounding of the axion and pair instabilities is that lighter BHs can be formed at the upper edge of the upper black hole mass gap. The BHs reappear at high masses (M BH = 133M in the SM) because the core temperatures and densities are sufficient for photodisintegration reactions to occur. These result in their own independent instability because the temperature and density increase due to small contractions further photodisintegrates the heavy elements rather than raising the pressure to counteract the contraction. Photodisintegrations compound with the pair-instability to lower Γ 1 further below 4/3 to the point where the oxygen ignition is not sufficient to reverse the contraction and core collapse ensues. One would therefore expect a similar phenomenon resulting from the compounding of the axion and pair instabilities. This compounding happens at smaller temperatures and densities than those required for photodisintegration reactions, implying that the oxygen explosion will not be sufficient to reverse the implosions of some lower mass objects.

B. Location of the Upper Black Hole Mass Gap
The results of our simulations are presented in figure 2 where we plot the BH mass as a function of initial helium core mass. The metallicity was taken to be Z = 10 −5 corresponding to population-III stars. The qualitative effects described above are evident. The pulsations and total disruptions (AISN) indeed begin at lower masses compared with the SM, and stars that would have experienced pulsations are instead instead totally disrupted, leaving no BH remnant. Furthermore, the BHs reappear at smaller masses compared with the SM. Our simulations predict that, when electron-mass axions are included in the EOS, the lower edge of the upper BHMG is located at 37M and the upper edge of the BHMG is It is seen that the value of MBHMG is lower than the value found in the Standard Model (see [10]) by more than 10 M .
located at ∼ 107M (compared to M BH = 133M in the SM). Thus, we predict that the BHMG is narrowed from 86M to 71M . We have run additional grids with Z = Z /10 (corresponding to population-II stars) where the solar metallicity is Z = 0.0142 and where the rate of the 12 C(α, γ) 16 O is varied by ±1σ from its median value. In all cases we find 4M changes in the location of the lower edge of the BHMG. Thus, our predictions are similarly robust as the predictions of the PISN [2].

C. Black Hole Population Signatures
The results of our simulations are shown in the left panel of Fig. 3. To understand the function M BH (M i ), we perform a numerical fit to a continuous, parametric function which can reproduce the shape of the curve: a constant, two power-laws with arbitrary coefficients to model both the BHs unaffected by pair-instability and PAISN BHs, and an exponential fall-off capturing AISN, as introduced in [10]. In the right panel, we perform a fit to the black hole mass function introduced in [10]: It is seen that the proposed mass function approximates the shape of this MESA-fit mass function well with just three parameters. The best fit parameters are also given in Fig. 3. Strikingly, it is seen that the model (5)  (54 +6 −6 M if GW190521 is excluded from the analysis) [10], which seemingly disfavors the possibility of AISN. However, for axion masses m a 2 MeV, we have determined that stars will encounter the PISN instead of the AISN. In this case, the functions M BH (M i ) and dN BH (M BH ) are the same as if there is no axion at all. Thus, we suggest the intriguing possibility that a spatially varying axion mass, such as could be generated by locally varying CP-violation in different host galaxies [40], could lead to two distinct peaks in the black hole mass population. Finding a complete, cosmologically viable model including this physical effect is an important target for further research.
Black holes beyond the upper edge of the upper BHMG are currently outside the detection range of LIGO/Virgo/KAGRA, but planned upgrades to their sensitivities and the addition of planned future interferometers to the network (e.g. LIGO-India) will enable the location of the upper edge to be measured with percentlevel precision [41]. Future detectors such as Cosmic Explorer (CE), Einstein Telescope (ET), and LISA will also be able to measure the location of the upper edge. It is therefore possible that the location of the upper edge will be measured within the next decade, at which point it could be used to place additional bounds on heavy axions.

D. Optical Signatures
Before concluding, we briefly comment on the possibility of directly detecting AISN using optical observations. AISN will be more common than PISN since their progenitors have lower masses, and lower mass stars are more numerous due to the negative slope of the stellar initial Additionally, we find that at fixed progenitor mass AISN are brighter than PISN owing to the more violent oxygen explosion needed to counteract the stronger contraction. We plot the luminosity of a 52M and 70M AISN alongside a 70M PISN in figure 4. Evidently, a 52M AISN is brighter than a 70M PISN, whereas the 70M AISN is brighter by more than an order of magnitude. JWST will therefore see more AISN than PISN, making AISN excellent targets. It would be interesting to use population synthesis and semi-analytic N-body codes to predict the AISN rate, and to use radiation hydrodynamics codes to predict their light curves. Such undertakings are beyond the scope of this work, but would provide additional detection channels.

V. SUMMARY AND OUTLOOK
In this work we have explored the effects of heavy axions in an unconstrained region of parameter space -the cosmological triangle -on the structure, evolution, and fate of massive stars. At these masses (m a ∼ m e ), these particles can remain in the star and equilibrate with the stellar plasma, leading to a new instability akin to the pair-instability. We have explored the ultimate effects of this instability by incorporating heavy axions into the equation of state module of the stellar structure code MESA.
Our simulations predict that the axion-instability results in a contraction that ultimately unbinds the entire star. Thus we predict the existence of a new optical transient, an axion-instability supernova.
We investigated the axion-instability supernovae induced by an axion of mass m a = m e . We explored the impacts on the upper black hole mass gap and found that both the upper and lower edges are shifted to smaller masses. In particular, the lower edge is located at 37M , compared to the SM prediction of 49M , and the upper edge is located at 107M , compared to the SM prediction of 133M , as shown in Fig. 2. Thus, the width of the gap narrows. The location of the lower edge is shifted to lower masses because the axion-instability manifests at higher densities than the pair-instability and can therefore be encountered by lower mass objects. The location of the upper edge is shifted to lower masses because the axion and pair instabilities compound as a result of axions in the cosmological triangle having similar masses to the electron. The resulting contraction is stronger than in the SM and the oxygen explosion is not sufficient to reverse it and prevent core collapse. We have checked that these results are robust to variations in metallicity and the 12 C(α, γ) 16 O rate. We also determined that axion-instability supernovae disappears for axions of mass m a 2 MeV, reverting to the standard astrophysical prediction.
We studied the signatures of AISN on black hole populations. In the case of the lower edge, we demonstrated that an analysis of the population statistics of BHs in the LVK gravitational wave catalogs can distinguish between the two models. This could provide a robust detection of heavy axions in the cosmological triangle by fitting the BH mass catalog to the fitting function in Eq. 5, which enables extraction of the location of the lower edge of the mass gap and the sharpness of the peak of the BH mass spectrum. Although these results appear to be in tension with current data, if there is spatial variation in the axion mass, then this could lead to a "smearing" of the peak in the black hole mass function, or even a mass function with multiple peaks. Evidence in favor of this possibility will be interesting to search for in black hole population catalogs. Planned upgrades to the sensitivities of LVK and the inclusion of additional detectors within the next decade will enable a robust measurement of the location of the upper edge, at which point independent corroborations can be placed.
A central assumption we have made in our analysis is that the new particles fill out the entirety of the local thermal phase space. This is inaccurate if they are produced with velocities greater than the escape velocity and have a mean free path much larger than the core size. In this case, they would contribute to anomalous energy transport in the star. In light of the findings here, a suite of out-of-equilibrium, momentum-dependent simulations of such new particles will be valuable to conduct. We have also restricted to a single, spin-0, CP-odd degree of freedom, motivated by the existence of the "cosmological triangle" for axions, but our methods can be used to investigate the effects of new bosons with arbitrary masses and spins.
The MESA code used to generate our results is included in the reproduction package for this paper [27]. Also included are the codes necessary to generate modifications to the MESA EOS due to bosons with usersupplied masses and spins. We encourage the use of these for further science studies of new heavy particles coupled to the SM.
Finally, we note that it would be interesting to predict the optical and neutrino signatures of axion-instability supernovae. Optical signatures could be predicted by generating synthetic light curves using photometric evolution codes, and neutrino signatures could be predicted using full three-dimensional simulations of axioninstability superovae. We found that AISN are brighter and occur more frequently than PISN (we gave an example in Fig. 4), raising the tantalizing possibility that they may be directly detected by JWST. The MESA EOS module requires 16 thermodynamic variables to be specified. Three of them -the mean molecular weight per gas particle µ, the mean number of free electrons per nucleon µ −1 e , and the ratio of the electron chemical potential to k B T (degeneracy parameter η) -are not modified by the presence of new particles. The remaining 13 are listed in table I. In addition to these, derivatives of each of these quantities with respect to ln T and ln ρ must also be provided. The EOS (and derivatives) are critical for determining the MESA time-step so must be accurate to high precision.
The default MESA EOS at temperatures and densities relevant for post-main-sequence massive star evolution is the HELM EOS [42]. Our procedure for including new particles in the EOS is to first call this EOS and then modify each of the quantities listed in table I. The MESA EOS uses density and temperature as input variables. Consequentially, the contribution of new particles to quantities defined as derivatives at constant ρ and T e.g. c V and χ ρ , as well as those with no derivatives e.g. pressure can be found by calculating the corresponding quantity for new particles and adding this to the MESA default value. Quantities defined as derivatives with other variables held constant e.g. Γ 1 (which involves constant entropy) and c P (which involves constant pressure), cannot be added in this manner. The reason for this is that the MESA default value is calculated at constant default quantity which, with the exception of T and ρ, is altered by the presence of new particles. In these cases, one must calculate the relevant quantity using thermodynamic relations when all species are included simultaneously (visible matter and new particles). We have derived formulas that can be used to calculate these quantities from a combination the MESA default EOS quantities, which we denote using overbars, and the EOS variables that can be updated to include new particles by simply adding their contribution. These formulas are listed in the final column of table I.
Our procedure is then to first calculate all of the EOS variables that can be found by adding the contribution from new particles, and then to use these to calculate the remaining quantities. This mandates that the EOS variables be calculated in a specific order since some nonadditive EOS variables must be specified before others can be calculated. From top to bottom, the order in which the quantities appear in table I reflects the order in which they are calculated in our MESA code. Given any single row, quantities from previous rows may act as inputs to the formula but not quantities from subsequent rows. The starting point for this process is the gas pressure P g , specific internal energy E, and specific entropy s. All other quantities can be calculated from combinations of these three and further quantities derived therefrom.
Assuming that the new particle is a spin-s boson φ with mass m φ and degeneracy g φ = 2s + 1 which is in thermal equilibrium with the stellar material, these three quantities are calculated analytically from first principles by integrating over the Bose-Einstein distribution as follows. We begin by defining where m φ is the mass of the new particle. The relevant thermodynamic quantities are the pressure, density, internal energy, and specific entropy given respectively by where The assumption that the phase space distribution of axions follows its thermal equilibrium form is nontrivial, since the production processes are athermal. Nevertheless, thermal broadening, redshifting, and rescattering will push the distribution towards thermalizing. Exploring the process of thermalization in detail is an interesting topic for future work.
The formulas for the EOS variables that we use in our MESA code are given in table I. The formulas for  the derivatives of these variables with respect to temperature are given in table II and the formula for the  derivatives with respect to density are given in table III. The derivatives are crucial for determining the MESA timesteps and cannot be omitted. In order to enable rapid computation and maintain high precision we have derived analytic fitting functions for P φ , u φ , and s φ , and their derivatives with respect to temperature and density. Specifically, we fit an eighth-order polynomial in T to each formula (the ρ-dependence is simple and can be included analytically) and evaluate the formulas in the final columns of the tables in MESA using these fitting functions. A copy of the specific fitting formulas used in this work (corresponding to a 511 keV scalar) is included with the reproduction package for this paper (see https://zenodo.org/record/6347632 [27]). Our reproduction package includes a Mathematica script that generates these fitting formulas for arbitrary boson masses m φ and degeneracies g φ for immediate use. We have defined the quantities X = P χT and Y = ρT cv. The specific enthalpy is h = E + P/ρ. The gas pressure is defined as the total of all sources of pressure except the radiation pressure P rad = aT 4 /3 with a the radiation constant. The total pressure is then P = Pg + aT 4 /3. Quantities with an overbar are those returned by the MESA default EOS and do not include DM. The formulas in the final column are those used in our modified MESA code.