Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle

We consider the constraints from Supernova 1987A on particles with small couplings to the Standard Model. We discuss a model with a fermion coupled to a dark photon, with various mass relations in the dark sector; millicharged particles; dark-sector fermions with inelastic transitions; the hadronic QCD axion; and an axion-like particle that couples to Standard Model fermions with couplings proportional to their mass. In the fermion cases, we develop a new diagnostic for assessing when such a particle is trapped at large mixing angles. Our bounds for a fermion coupled to a dark photon constrain small couplings and masses<200 MeV, and do not decouple for low fermion masses. They exclude parameter space that is otherwise unconstrained by existing accelerator-based and direct-detection searches. In addition, our bounds are complementary to proposed laboratory searches for sub-GeV dark matter, and do not constrain several"thermal"benchmark-model targets. For a millicharged particle, we exclude charges between 10^(-9) to a few times 10^(-6) in units of the electron charge; this excludes parameter space to higher millicharges and masses than previous bounds. For the QCD axion and an axion-like particle, we apply several updated nuclear physics calculations and include the energy dependence of the optical depth to accurately account for energy loss at large couplings. We rule out a hadronic axion of mass between 0.1 and a few hundred eV, or equivalently bound the PQ scale between a few times 10^4 and 10^8 GeV, closing the hadronic axion window. For an axion-like particle, our bounds disfavor decay constants between a few times 10^5 GeV up to a few times 10^8 GeV. In all cases, our bounds differ from previous work by more than an order of magnitude across the entire parameter space. We also provide estimated systematic errors due to the uncertainties of the progenitor.


Introduction
In 1987, a core-collapse supernova known as Supernova 1987A (SN1987A) was observed in the Large Magellanic Cloud. SN1987A provided a wealth of information on the supernova explosion itself and also sets unique constraints on the existence of new, low-mass particles that are weakly-coupled to the Standard Model (SM) [1,2]. The existence of new, weakly-coupled particles could provide novel channels to "cool" the protoneutron star and change the neutrino emission from SN1987A. Constraints in this context are derived from the "Raffelt criterion": the luminosity carried by the new particles from the interior of the protoneutron star environment to the outside of the neutrinosphere must be smaller than the luminosity carried by neutrinos [2]. The observed cooling time of the supernova agrees within uncertainties with the SM prediction [3,4]. However, if there were an additional efficient channel for energy flow that could compete with neutrinos, the cooling time of the supernova would have been shorter than observed.
SN1987A provides a hot and dense stellar environment, so even very weakly-coupled particles could have been produced. Since the supernova core temperature, T c , is about 30 MeV, particles with mass less than about a few hundred MeV can be constrained when taking into account the Boltzmann tail. In this paper, we will derive constraints from SN1987A on several possible low-mass particles: various dark sectors consisting of dark matter (DM) and dark photons (including millicharged particles), the QCD axion, and axion-like particles with Yukawa couplings.
A popular and prototypical model for sub-GeV DM is given by a dark sector consisting of a DM particle, χ, interacting with a dark photon, A [5][6][7][8][9][10][11]. Kinetic mixing between the dark photon and the SM photon leads to an interaction between the DM and electrically charged particles of the SM. Various hierarchies between the DM mass, m χ , and the dark photon mass, m , lead to a diverse range of phenomenology. In several cases, sharp benchmark targets can be identified in parameter space for which the DM interacting with a dark photon can obtain the correct DM relic density [11][12][13][14][15] (for a review, see [16,17]). In this paper, we will consider "heavy" DM with m < 2m χ , "light" DM with m > 2m χ , and a dark-sector with a millicharged particle (the latter bounds apply equally well to DM that interacts with a massive, but ultralight, A ). In addition, we will consider "light", inelastic DM, in which the dark-sector consists of two "DM" particles, χ 1 and χ 2 , which have a small mass splitting, and for which the interaction with the dark photon is off-diagonal, i.e. ∼ A χ 1 χ 2 .
In deriving the SN1987A constraints on these particles, we include the thermal effects on Aphoton mixing [18,19], which are very important for the SN1987A constraints on a dark sector consisting of only dark photons [20,21] (for previous bounds see [22][23][24][25][26]). As we will see, the presence of the DM particles changes the constraints in a significant way from the A -only case, even when the A is kinematically forbidden to decay to the DM directly. We also use a novel criterion to calculate the couplings for which the DM is "trapped" inside the proto-neutron star (and thus does not contribute to the cooling): we require them to take a random walk until in their velocity vector is turned by 90 • from their initial direction of motion. The SN1987A constraints on the light DM scenario had been considered previously in [27][28][29][30], but our analysis goes well beyond these references. Our analysis also significantly updates previous work on millicharged particles [31].
In addition to the dark-sector models mentioned above, we also revisit the constraints on another popular and important particle, the QCD axion [32][33][34]. Previous bounds have been extracted with a range of simplifying assumptions, which we attempt to rectify by including additional estimates of known nuclear physics as well as particle physics effects. We find significant differences with the constraints in the literature [35]. Finally, we also revisit constraints on axion-like particles with couplings proportional to the SM Yukawa couplings, updating bounds from [36].
The remainder of the paper is organized as follows. In Sec. 2 we discuss the production of DM and dark photons and describe in some detail how we calculate the luminosity of these dark-sector particles. We discuss the case of "heavy" and "light" DM, defined by whether the dark photon is lighter or heavier than twice the DM mass, respectively. We also discuss variants of the basic DMcoupled-to-a-dark-photon model. This includes inelastic DM, in which case additional mass terms in the Lagrangian allow the fermion states to have different masses and off-diagonal couplings to the dark photon; and the case where the dark photon is essentially massless, so that the DM appears to have a "millicharge." Sec. 3 describes the results for these various models. In Sec. 4 we change gears entirely and address the QCD axion, while Sec. 5 discusses axion-like particles. In Sec. 6 we conclude. We leave many details of our calculations to various Appendices.

Model Description and Preliminary Comments
We consider a variation of the model examined in [20], in which the only new light particle was a dark photon of a new U (1) gauge group that kinetically mixes with the SM hypercharge gauge boson. This dark photon is a massive vector boson with a small coupling to electrically charged particles. Here, we assume that the dark sector also includes a Dirac fermion, χ, charged under the U (1) , that is light enough to be produced in nucleon-nucleon collisions in the proto-neutron star. Thus, pair production of χχ states provides an additional channel through which energy can flow. The low-energy Lagrangian describing such a dark sector is where is the kinetic-mixing parameter, g D is the dark gauge coupling, m (m χ ) is the dark photon (χ) mass, F µν the usual electromagnetic field-strength tensor, and F µν the field-strength tensor of the U (1) gauge boson. We define α D ≡ g 2 D /4π as the "dark fine-structure constant." Conservation of χ-fermion number guarantees that these particles are stable even below the scale of U (1) symmetry breaking, which is why they are a DM candidate. For this reason, we refer to χ as "DM" in this paper, although we will not address its early-universe production nor its cosmological effects. Moreover, as long as χ is stable on the time it takes to escape the proto-neutron star, the SN1987A constraints derived below are applicable even if χ is only a fraction of the DM. Of particular interest for understanding their behavior in SN1987A, conservation of χ-number means that the dark fermions can scatter off SM particles as they make their way out of the star. This affects the energy spectrum of the χ particles and the A compared to the scenario with a solitary A , and qualitatively changes the notion of a trapping limit at large mixing angle. Furthermore, χχ-pairs may be produced through off-shell dark photons, which can avoid the suppression from thermal effects on A -photon mixing described in [20]. As a result, the lower bounds have a different low-mass limiting behavior and become stronger compared to the A -only case.
Our constraints are derived explicitly for a dark sector with a Dirac fermion coupled to the dark photon. However, the constraints should be very similar for a dark sector consisting of a complex scalar coupled to the dark photon. The production and (relativistic) scattering cross sections are slightly modified, but the particle number is still conserved, so the kinematics of scattering are similar. We expect differences due to degree of freedom counting and details of the cross sections to give O(1) corrections to the fermionic bounds. These differences are below the systematic uncertainties due to imperfect knowledge of the supernova temperature and density profiles, discussed in more detail below. Thus, our bounds can serve as rough guidelines on the parameter space of a dark photon coupled to a dark charged scalar.
In what follows, we ignore the presence of a dark Higgs boson, which can affect the phenomenology if the U (1) is broken through a dark Higgs mechanism and the Higgs boson remains very light. The dark Higgs mass is determined by its self-coupling parameter, which we will henceforth take to be large enough that the dark Higgs is heavier than the dark photon and DM and kinematically inaccessible during the supernova. However, while there are additional dark-sector production and decay modes to consider, we claim that these should not lead to a significant change in the supernova cooling rate compared to the rates we derive here. Processes involving dark Higgs production do not suffer suppression from the well-known plasma effects in on-shell A -photon mixing, but such suppression is also absent for DM-pair production. Thus, the A -only case explored in [20] is in some sense unique, and the results in this work should be qualitatively similar even with a light dark Higgs boson.
As we discuss in more detail below, there are several important processes by which DM particles are created in the proto-neutron star environment. The dominant processes at low DM mass are bremsstrahlung of a DM pair during nucleon collisions through an on-shell or off-shell A and SM photon decay to a DM pair in the plasma if m χ ω p,0 , where ω p,0 ≡ ω p (r = 0) ∼ 15 MeV is the plasma frequency at the center of the supernova (see App. A). After being produced, DM particles scatter against nucleons and electrons on their way out of the proto-neutron star environment. Similar to the A -only case, where the rate of A production as well as decay or absorption are proportional to the same parameter ( 2 ), the rate at which the DM particles are produced and the rate at which they scatter are proportional to the same parameter combination (α D 2 ). At increasingly large mixing angle, the DM particles will scatter multiple times during egress, potentially becoming trapped and even returning to chemical equilibrium inside the proto-neutron star. We find that there is some parameter space where α D 2 is large enough that a sufficient number of dark fermions are produced to alter the evolution of the supernova explosion but small enough that the dark-sector particles scatter infrequently on their way out of the star. Thus, for a large range of DM masses, there is both a lower and upper bound on the coupling to the SM. We emphasize here that if the elastic scattering of DM particles is extremely frequent, they may be unable to exit the supernova. If χ andχ particles proliferate throughout the star, DM annihila-tion or pairwise inverse bremsstrahlung will equilibrate everywhere, attaining a thermal abundance on a timescale set by their production rate [37]. 1 Equipartition of degrees of freedom in the thermal bath then implies that the overall temperature as well as the transport properties of the proto-neutron star will change. Although introducing this many new degrees of freedom may have unacceptable consequences for the behavior of the supernova explosion, such effects are hard to resolve analytically and potentially beyond the present day understanding of the proto-neutron star interior, and thus are beyond the scope of this work. Here we will calculate the mixing angle that gives a decoupling radius close to the neutrinosphere using conservative analytic requirements, and we assume a benign thermal population at higher mixing angles. We will show these as plausible upper bounds for the kinetic mixing parameter above which the effects of the dark sector are best investigated with other, laboratory-based techniques, as discussed in e.g. [17]. It would certainly be interesting to use simulations to investigate more comprehensively the effects of a thermally-equilibrated dark sector population on supernova explosions.

Model Variation: Inelastic Dark Matter
The model above describes a Dirac fermion coupled to a dark photon. An economical, UV-complete model of dark-sector masses is provided by the introduction of a doubly-charged dark Higgs field h D that experiences spontaneous symmetry breaking with a nonzero vacuum expectation value h D . If the U (1) -charged fermions have Yukawa couplings to the dark Higgs as well as a U (1) -invariant Dirac mass, the mass eigenstates can undergo a small splitting after symmetry breaking. The Lagrangian for this scenario is where χ = (λ ξ † ); λ, ξ are two-component Weyl fermions; and m ξ , m λ ∝ h D arise after symmetry breaking. A field redefinition allows us to work in terms of mass eigenstate fields χ 1 , χ 2 . These fermions have different masses and in principle can couple to the dark photon inelastically as well as elastically: with We define the mass splitting ∆ ≡ m 2 − m 1 = m ξ + m λ . The lighter of the two Majorana fermions, χ 1 , could be a DM candidate. We refer to this model as "inelastic DM" [38], since the scattering of χ 1 off SM particles could be dominated by a transition from χ 1 to the heavier state χ 2 . This model was studied in [39,40] in the sub-GeV mass range of interest in this paper. If m ξ = m λ exactly, the elastic Figure 1. The Feynman diagram for the interaction between χ and Standard Model particles f , which have a charge qe, where e is the electron's charge.
coupling vanishes and only inelastic scattering is possible at tree-level. We will calculate below the SN1987A constraints for this simple variant of the inelastic DM model. We note that even in this case χ 1 could scatter elastically at one-loop, although this is highly suppressed, as we will discuss further in Sec. 3.3.

Model Variation: Millicharged Particles
A millicharged particle is a dark-sector particle with a small electric charge. One well-motivated way to attain such small charge is by introducing a massless dark photon, which can be removed by a field redefinition of the SM photon. Under such a field redefinition, dark-sector particles coupled to the A acquire a small coupling to the SM photon, g D , where is again the kinetic-mixing parameter between SM photon and dark photon. We define Q ≡ g D /e so that dark-sector particles charged under the dark photon have an electric charge Qe. We will derive below the constraints from SN1987A on millicharged particles, which update bounds presented in [31]. These bounds are equally applicable for DM particles that couple to an "ultralight", massive dark-photon mediator, with a mass below the typical momentum transfer in DM production and scattering processes, i.e. so long as the mediator mass can be neglected in the calculations [15].

Dark-Sector Particle Production in the Proto-Neutron Star
Dark photons and particles charged under U (1) can be produced in the proto-neutron star through the kinetic mixing between the SM photon and the dark photon, see Fig. 1. The total luminosity in dark-sector particles is L dark = L χ + L A (note that L χ denotes the total luminosity in χ andχ), and the criterion L dark = L ν determines the boundary of constraints, where L ν = 3 × 10 52 erg / s is the neutrino luminosity at one second [2]. The dominant production mechanisms for DM particles are through bremsstrahlung (via an off-or on-shell A ) during neutron-proton collisions and through SM photon decay in the plasma, shown in the left and middle panels, respectively, of Fig. 2. For on-shell A production, bremsstrahlung dominates, while for χ production both bremsstrahlung and SM photon decay are important, with the latter dominating by a factor of a few. Before discussing the calculation of L dark , it is worth making some general comments on the DM production in the supernova. We then describe the calculation of L A and L χ in Sec. 2.3 and Sec. 2.4, respectively, leaving detailed formulae to the Appendices.

Sam McDermott
Due to plasma effects in the proto-neutron star interior, we find it convenient to calculate scattering amplitudes in the gauge boson interaction basis, as in [18]. With this choice, all Feynman diagrams describing the interaction of DM with electrically charged particles such as the proton implicitly contain the diagram of Fig. 1, for which the amplitude is where K µ = (ω, k) is the momentum four-vector of the intermediate state (carried by both the SM photon and the dark photon), Π D is the self-energy of the dark photon in a plasma of DM particles, P T µν and P Lµν are the transverse and longitudinal projection operators of the SM polarization states, respectively, and Π T and Π L are polarization tensors of the SM photon from thermal effects. The dark photon absorptive width Γ is dominated by its decay width to DM when this decay is on shell, In principle, the presence of the dark photon self-energy Π D suggests that we should include separate longitudinal and transverse projection operators for the dark photon like we do for the SM photon. However, Π D is negligible on the lower boundary of the excluded parameter space where we calculate dark-sector production rates, since the dark-sector particles free stream. The effect of Π D is only important near the upper boundary where the dark-sector number densities can be very high. However, as we discuss in more detail below, in this part of parameter space it is safe to assume that number densities are simply given by a thermal distribution inside some radius, so their exact production rate will not be important. Thus, we will not use Π D in any explicit calculation.
The dark sector luminosity admits two kinds of resonances, as can be seen in Eq. (2.5): the on-shell peak from the dark photon propagator, attained for K 2 = m 2 , and the "thermal peak," at K 2 = Π L,T from the SM photon propagator. The on-shell peak dominates if m ω p,0 , the thermal peak dominates if m ω p , and both peaks can be attained for m ∼ ω p . Thus for m ω p , off-shell DM production is suppressed, and the lower bounds are same as the dark photon only case. However, off-shell DM production dominates for m ω p,0 , so the low-bounds at small m are stronger than the bound from the A -only case, which decouples like m 2 due to suppression by thermal effects [18,20]. With the inclusion of dark sector fermions, the production rate for darksector particles becomes independent of their mass, so the lower bound is flat.

Dark Photon Luminosity (L A ) for Small Couplings
We now discuss the calculation of the luminosity L A for small couplings and α D . This is similar to the pure A case discussed in [20] only if the DM particles free stream out of the proto-neutron star. If, however, the couplings are large, then the dark matter abundance is also large and the A will experience a large optical depth. We will discuss our treatment of the bounds at large coupling values in Sec. 2.5; Eq. (2.12) gives the A luminosity for large couplings. In this section, we will ignore dark-sector interactions.
Bremsstrahlung production of the A dominates over purely electromagnetic-like processes such as semi-Compton scattering, because the QCD coupling α s ∼ O(4π) α EM 1/137 and because nucleons are highly abundant but not Pauli blocked like electrons. SM fermion annihilation contributes negligibly because the chemical potential of all such particles is very high and their antiparticles are very scarce. Mixed nucleon scattering dominates over proton-proton scattering because the former emits dipole radiation while the latter only emits like a quadrupole [26,41]. The "direct luminosity" in A particles is where R ν is the neutrinosphere radius, τ is the optical depth, Γ br is the dark photon production rate via bremsstrahlung, Γ ibr is the inverse bremsstrahlung rate, Γ e (Γ χ ) is the A decay width to electromagnetic (DM) particles, and Γ dC (ω, r ) is the rate for "dark Compton" scattering (e.g. A χ → A χ) that only contributes when the χ are trapped (see Sec. 2.5). The far radius R far is taken to be the neutrino gain radius R g 100 km [20,42]. App. A contains the definitions of all these rates.
In principle, DM annihilation as well as semi-Compton scattering involving a single SM photon and a single dark photon can also contribute to the power, but these are negligible unless the DM is trapped, which we explicitly ignore for the time being. The widths Γ ibr and Γ e are suppressed by 2 , and Γ χ is nonzero only if m > 2m χ . If m > 2m χ dark photons decay to DM particles on very short distances compared to the size of the proto-neutron star (unless α D is very small), which sends e −τ (ω,r) → 0 such that L A 0 and L dark L χ .

Dark Matter Luminosity (L χ ) for Small Couplings
There are two main contributions to the DM production rate: (i) bremsstrahlung of DM pairs in proton-neutron collisions, which we call L b χ , and (ii) SM photon decays in the thermal plasma, which we call L d χ ; see left and middle panel of Fig. 2. Assuming that the DM particles do not scatter on their way out of the star (valid for small values of α D 2 ), it is straightforward to calculate the resulting dark-fermion luminosity, L χ = L b χ + L d χ . We will give the corresponding expressions in Sec. 2.4.1 and Sec. 2.4.2.
For large values of α D and , the DM particles may scatter and thermalize with SM material, rendering their escape energy different from their energy at production. Since DM number is conserved, the dark fermion luminosity then does not have a description analogous to Eq. (2.7). Moreover, as mentioned in Sec. 2.3, the dark-photon luminosity L A needs to be modified from Eq. (2.7) in the presence of DM particles and large values of α D and , since in this case the A can scatter off DM particles, which conserves A number. For large couplings, the calculations of L χ and L A then require defining a "trapping criterion". Various choices for such a definition are in principle possible; we describe our criterion in Sec. 2.5.

Bremsstrahlung of Dark Matter Pairs
DM pairs can be produced by bremsstrahlung in proton-neutron collisions. If the DM particles do not scatter on their way out of the star, the differential luminosity per unit volume is 3 is the Lorentz-invariant phase space of the particle with four-momentum P µ = (E i , p i ), the nucleon phase space densities are f i , and the incoming (outgoing) nucleons have momenta p 1 , p 2 (p 3 , p 4 ). As in [20] we employ the soft radiation approximation to calculate the matrix element, which is valid for the mass range in which E χ + Eχ | p N | 2 /2m N [26]. Details of the calculation of Eq. (2.8) are in App. B. The luminosity due to emission of DM in this limit is the dV , assuming that dark matter particles free stream out of the star.

Standard-Model Photon Decay in the Thermal Plasma
In vacuum, dark-sector particles charged under U (1) do not couple to on-shell SM photons because the mixing term is K 2 and K 2 = 0 for on-shell photons. Since the dispersion relation for photons is altered in a thermal plasma, however, this coupling does arise. The SM photon dispersion relation picks up a real part Re(Π L,T ) ≤ 3/2 ω p (where ω p is a function of the distance from the center of the proto-neutron star, r) and a nonzero imaginary part. DM with mass less than 3/2 ω p /2 can therefore be produced from the decay of SM photons, shown in Fig. 2, much like the plasmon process that leads to neutrino production [43].
We can write the differential luminosity of the DM from SM photon decay as where K µ T,L = (ω T,L , k) is the photon four-momentum and Γ d T,L is its decay rate to a DM pair (2.10) k k is the longitudinal polarization vector, and θ χ is the angle between the incoming photon and the outgoing DM. The amplitudes in Eq. (2.10) are not Lorentz invariant because the plasma frame breaks Lorentz invariance, and the transverse and longitudinal modes have different dispersion relations. In Eq. (2.9), the momentum integral of the longitudinal mode is cut off at k max , the largest three-momentum that the longitudinal photon can have, while transverse photons can have any value of k [43].
The luminosity of DM from Standard Model photon decay is a factor of a few higher than from bremsstrahlung for the mass range in which photon decay is possible; it is also independent of m χ for m χ ω p . The luminosity due to emission of DM in this limit is the volume integral of Eq. (2.9), , assuming that dark matter particles free stream out of the star.

Dark-Sector Luminosity and Trapping Criterion for Large Couplings
If the DM scatters many times on its way out of the star, it can equilibrate with the SM particles. If the DM is in equilibrium, dark photons will be trapped as well. It is therefore important to define a decoupling radius R d for which DM is in equilibrium inside and free-streaming outside. 2 For a given m χ and m A , this radius is obviously a function of and α D . We shall assume that dark sector energy emission from R d is free-streaming and thermal and the luminosity is an effective blackbody. We set the upper bound of to be where blackbody emission from R d (α D , ) equals L ν ; we specify an algorithm for computing R d (α D , ) below. For larger values of α D or , the decoupling radius increases and the energy emission decreases because temperature falls sharply with radius. This is reminiscent of the calculations done for sterile neutrino emission from SN1987A [44]. The blackbody luminosity of a fermion from a radius R d is 2 In fact, there are two "decoupling radii," one each for chemical and kinetic equilibrium. The chemical decoupling radius r cd determines the number density of DM particles emitted from the supernova, and the kinetic decoupling radius r kd determines the energy of each DM particle being emitted. The chemical decoupling radius should be smaller than the kinetic decoupling radius because kinetic equilibrium is necessary for chemical equilibrium in this model. However, these radii are difficult to find exactly without simulations. As a conservative assumption, appropriate considering all the other uncertainties of the problem, we solve only for the radius of kinetic equilibrium. This is conservative because it provides a lower limit on the number density of DM particles, since r cd ≤ r kd in reality.
where E χ = p 2 χ + m 2 χ , v χ = p χ /E χ , g χ = 4 counts degrees of freedom, and the final approximation assumes a massless fermion. The blackbody luminosity of a dark photon is To derive a constraint on (given some choice for α D and other model parameters), we (i) calculate the radius R * d at which the dark sector blackbody luminosity equals the neutrino luminosity, L dark (R * d )| therm. = L ν , and then (ii) find the value of that gives thermal decoupling at this radius.
Step (i) is computationally straightforward given that Eqs. (2.11) and (2.12) are simple to compute for a given mass. However, step (ii) is more involved, and we discuss our approach next.
Since the "decoupling radius" is only an approximate concept and finding the zone of decoupling is impossible without simulations, we propose a simple criterion: the kinetic decoupling radius R * d is defined to be where the expected angular deflection of a thermal DM particle starting at R * d and ending at R f is π/2. At smaller radii r < R * d most DM particles are redirected and find antiparticles to annihilate with, while at larger radii most DM particles escape and drain energy from the supernova explosion. To calculate the expected angular deflection, we first define the total number of scatters and the maximum angular deflection that a DM particle would experience if it was scattered in the same plane and in the same direction every time: Here, Γ s is the event rate for χ + p → χ + p elastic scattering, ∆θ is the average angular deflection per collision, and E(R d ) is the thermally averaged energy at the radius R d . (We ignore the energy change of the DM after each collision, since it is small. We also ignore the difference in total path length due to the angle change, since this is a higher order effect.) Most particles will not move in the same angular direction upon each scattering, of course, but instead will take a random walk in solid angle. The expected displacement due to a random walk is the mean of the chi distribution in d dimensions, which differs from the root-mean-square deviation θ max / N steps by a factor 2/d × The expected angular deflection for a typical particle is thus (2.14) , so the expected angular deflection and the upper bound given in Eq. (2.14) rise linearly in √ α D ; as a result, if we choose a different critical angular deflection, for example |θ(α D , u , R * d )| = π instead of π/2, our bounds would be twice as restrictive. Details of the calculation of Eq. (2.13) are given in App. C.
We emphasize that our approach to calculating a constraint for large values of is quite different when DM is present compared to the A -only case. In [20] we calculated the dark photon energy emission by weighting the differential power from all radii with the probability of escape, e −τ , regardless of the value of . Since both the power and τ have nontrivial energy dependence, we found that the luminosity in dark photons is dominated by higher energies at higher mixing angle, though for large enough the Boltzmann suppression becomes important and the total luminosity decreases. This calculation is valid because dark photons do not survive scattering with SM particles. However, DM particles can elastically scatter many times and still escape, as can a dark photon that scatters off of DM. At large mixing angles, the dark sector energy distribution at escape may therefore be different from the energy distribution at production, unlike in [20]. We also note that the DM elastic scattering cross section can be forward peaked if the mediator is light compared to the typical momentum transfer (a few MeV at the supernova core), so calculating the mixing angle at which Rν R d drΓ s 1 is misleading, as it is reasonable to expect that the DM can scatter at least once on its way out of the proto-neutron star without returning to chemical equilibrium. Our more involved calculation is necessary to obtain accurate limits.
When calculating the upper boundary for large DM masses for the inelastic DM scenario, we will revert to a much simpler criterion than the one discussed above: we will require, very conservatively, that the DM scatters only once. We will further explain and justify this approach in Sec. 3.3.
Finally, we note that for values of above our upper boundaries, DM inside of the proto-neutron star attains a thermal abundance out to radii larger than R * d . It is possible that some DM particles are able to escape and travel to the detectors that registered the SN1987A neutrinos. Upon arriving, they may be observed through elastic scattering with the water in the neutrino detectors. However, by assumption the thermal energy of DM particles at these masses is ∼ O(T (R d )) 3 MeV, which is below the threshold of the Kamiokande detector [45], so that only a few of the emitted particles far along the Boltzmann tail are detectable even in principle. In addition, unless the DM mass is very small, their arrival at the detector will be significantly delayed compared to the neutrino signal, and there will be a large spread in arrival times due to a large velocity dispersion.

Supernova Temperature and Density Profiles
The calculations of dark-sector particle production and their luminosity require knowledge of the temperature and density profiles of the proto-neutron star. There are large uncertainties in these profiles, and we thus use four different profiles in order to estimate the systematic uncertainties in our resulting constraints. We choose the same profiles as in [20], and will refer to these as the "fiducial" [2], "Fischer, 11.8M " [46], "Fischer, 18M " [46], and "Nakazato, 13M " [47,48]. The profiles from [46] use the AGILE-BOLTZTRAN code [49][50][51][52][53], while the profile from [47,48] is based on a solution to a neutrino radiative hydrodynamical code before shock revival and a solution to the flux-limited diffusion equation after cooling has commenced. See [20] for further details and comparisons of these simulations. (See also [54][55][56] for a qualitatively different explanation of the observed neutrino burst.)

Dark Matter Coupled to a Dark Photon: Results
The phenomenology of DM particles interacting with dark photons in the supernova depends on whether or not the dark photon can decay to DM, so we have to consider the two possible mass hierarchies between the A and χ separately. We thus study two scenarios: (i) "heavy dark matter", where we choose the specific mass relation m χ = 3m for illustration, and (ii) "light dark matter", where we choose the specific mass relation m = 3m χ for illustration. In the former case, the dark photon can be stable against decays on the supernova timescales, since the decay is suppressed by 2 ; in the latter case, all dark photons promptly decay to DM (we will assume that α D is large enough to allow this). The constrained parameter space is very similar in both cases regardless of the mass hierarchy, since the luminosity is approximately a blackbody at large couplings, see Eqs. (2.11) and (2.12). Even if the A is stable against decay to χχ, the increased optical depth from the abundant DM particles, manifested as Γ dC in Eq. (2.7), dramatically increases the A optical depth and reduces the energy released in dark photons.
We will also study the two model variations mentioned in Sec. 2.1, namely (iii) "inelastic dark matter" for several choices of the mass splitting ∆, and (iv) "millicharged particles".
We will not show any parameter space for m χ ≤ 100 keV, but, as discussed at length above, these bounds do not decouple in the small m χ limit. Note also that we can safely ignore the suppression due to the Landau-Pomeranchuk-Migdal effect, since this only suppresses production of particles with energy less than the quasiparticle width γ, with γ 5 MeV for the densities of interest here [57,58]. In all cases, we will compare the SN1987A bounds to laboratory bounds and projections for proposed experiments, as appropriate.

Heavy Dark Matter
For m < 2m χ , the energy is carried by both A and χ particles, since the A is stable at leading order in . Production of χ particles is independent of m χ as long as m χ < 3/2 ω p /2, but is Boltzmann-suppressed at high masses. In contrast, A production is suppressed at small m . For a given α D not too small, the lower boundary on the -constraint is thus determined from χ production when m χ is small and from A production when m χ is large. For large , both A and χ particles are abundant, but the A can experience a large optical depth if a dense "cloud" of DM particles is created in the explosion, 3 written as Γ dC (ω, r ) in Eq. (2.7). As discussed in Sec. 2.3, to determine the upper boundary on the -constraint, we use the more conservative bound (i.e. the bound that excludes less parameter space) between the criterion from Sec. 2.5 and from [20].
We show the resulting SN1987A bound in the left plot of Fig. 3. The black solid and dashed lines are constraints on the heavy DM model for two values of α D , while the blue dotted line is the constraint for the A -only case presented in [20]. The lower bounds are stronger for small m χ than in the A -only case because they are not lifted as m → 0, and they return to the same value as the A -only case for m ω p . For small m χ , as discussed above, the upper boundary lies below the A -only case due to the large contribution of Γ dC (ω, r ) to the optical depth.
For additional insight, we display the contour along which a typical χ scatters off a proton either once or ten times on the way out of the supernova with the brown dashed and dot-dashed lines, respectively. This diagnostic clearly gives us much less sensitivity than asking where the χ is expected to satisfy Eq. (2.14). This reflects a real physical effect: in order for a light DM particle coupled Figure 3. SN1987A constraints on "heavy" dark matter coupled to a dark photon, for the specific mass relation m χ = 3m . Left: Solid (dashed) black line shows the constraint for α D = 0.5 (α D = 0.005). Along the brown dashed and dot-dashed lines, the dark matter scatters once and 10 times, respectively, on its way out of the star, for α D = 0.5. The blue dotted line is the constraint on a dark sector that contains only a dark photon and no dark matter. We assume the fiducial temperature and density profile for the supernova. Right: Black lines are the same as in the left plot for the fiducial temperature and density profile, while colored lines are the constraints for the other profiles with α D = 0.5. through a light mediator to become trapped and return to chemical equilibrium, it must scatter much more than once on its way out of the proto-neutron star.
Varying α D changes the asymptotically flat part of the upper boundary in such that α D 2 is kept fixed, since this boundary is determined by the dark-matter-proton scattering cross section. In addition, the flat part of the lower boundary in (i.e. for small m χ ) also changes such that α D 2 is kept fixed, since that region is dominated by dark-matter pair production from bremsstrahlung (as opposed to A production with the A decaying to DM). In addition, while we do not show this explicitly, changing the mass ratio of m χ /m affects the value of m χ below which the lower bound becomes independent of m χ . In the top-right panel of Fig. 3, we show constraints for different temperature and density profiles as reviewed in Sec. 2.6 and given in [20]. The variation with different profiles can be taken as a systematic uncertainty on the bound. The upper boundary is higher for profiles that have a lower density beyond R ν (see Fig. 3 in [20]), since it is easier in this case for the χ and A to leave the proto-neutron star.
It is interesting to show the constrained region in relation to laboratory searches for this dark sector model. The left plot in Fig. 4 shows the SN1987A constraints together with the latest laboratorybased searches, including colliders, beam-dump and fixed-target experiments, and precision measure- ments [17,23,. The SN1987A bounds are complementary and constrain lower values of than these laboratory bounds. The plot also shows bounds from direct-detection experiments. The bounds assume that the χ make up all the DM. There are two types of direct-detection bounds: from electronrecoil searches and from nuclear-recoil searches. For the former, we use the constraints from [81][82][83], which are based on XENON10 [84], XENON100 [85], and DarkSide-50 data [83], while for the latter, we use the combined bounds from the CRESST [86,87], SuperCDMS [88], and LUX [89] collaborations. In order to put these bounds onto the versus m parameter space, we follow the definitions in [11,15] and define the direct-detection cross section for DM scattering off electrons (protons) as The electron-recoil searches constrainσ e , while the nuclear-recoil searches constrainσ p . Given a specific mass relation (we choose m χ = 3m ) and value for α D (we choose 0.5), we can display these constraints on the versus m parameter space. Note that for smaller values of α D , the direct-detection constraints will weaken as 1/ √ α D , as indicated on the plot.
We also show a combined projection from future collider and beam-dump searches (black dotted line) [17], as well as from SuperCDMS SNOLAB [90], SENSEI [15,91], and a possible search with an experiment using a silicon target sensitive to single electrons with a 1 kg-year exposure [15]; for other projections, we refer the reader to [17]. These projections are largely complementary to the SN1987A bounds.
One additional bound on this dark-sector model that we do not show on the plots can be relevant for large couplings and m χ 10 MeV, assuming the χ are all of the DM and in thermal equilibrium with the SM sector in the early Universe. This bound comes from constraints on the effective number of degrees of freedom, N eff , from Big Bang Nucleosynthesis and the Cosmic Microwave Background [92,93], since a light DM particle could affect the relation between the photon and neutrino temperatures after neutrino decoupling. Our SN1987A bounds are largely complimentary to this, since they apply to dark matter with small couplings.

Light Dark Matter
If the DM is "light", i.e. for 2m χ < m , the dark photons will decay to DM pairs quickly (assuming α D is not too small), so that all of the energy in the dark sector is in χχ pairs. For this reason, the Figure 6. Thick solid black, red, green, and red lines show the SN1987A constraints for various temperature and density profiles on "light" dark matter coupled to a dark photon, assuming the specific mass relation m = 3m χ , and α D = 0.5. Dashed black line shows the constraint for α D = 0.005 using the fiducial profile. Thick orange lines show several benchmark model "targets", along (or above) which the DM can obtain the correct relic abundance in various scenarios. Note that the SN1987A have been evaluated specifically for a DM particle that is a Dirac fermion (and not a scalar), while we show targets for both scalar and fermionic DM; however, the SN1987A constraints are expected to be similar in both cases (see discussion in Sec. 2.1). Existing laboratory-based searches are shown in gray, including colliders, beam-dump and fixed-target experiments that search for A → χχ decay. Left: Under the assumption that the χ is all of the dark matter, we show constraints from dark-matter electron scattering from XENON10, XENON100, and DarkSide-50, and constraints on dark-matter-nucleus scattering from the CRESST, SuperCDMS, and LUX collaborations. Dotted lines show projections SuperCDMS SNOLAB (green), as well as SENSEI and a possible experiment using a silicon target sensitive to single electrons with a 1 kg-year exposure (both blue). Right: Same as left plot, but in the y versus m χ parameter space. We again show the same accelerator-based constraints as in the left plot, but now show projections in dotted lines from Belle-2 (cyan) as well as the proposed experiments BDX (blue), LDMX (magenta), and MiniBooNE (dark green). See text for references and details.
lower boundary of the SN1987A constraint is determined from Eqs. (2.8) and (2.9) and the upper boundary is determined from the criterion in Eq. (2.14). We show the resulting SN1987A constraint in Fig. 5. We use the same profiles and parameters as in Fig. 3, except we now choose m = 3m χ instead of m = m χ /3. The qualitative features in both scenarios are largely the same, except that the upper boundary now increases at large m relative to Fig. 3. This is because it is harder to trap the DM particles compared to the A alone, so that the decoupling radius R d decreases for m χ T c , where T c 30 MeV is the supernova core temperature.
It is again instructive to show the SN1987A constraints for this particular dark-sector model together with other, laboratory-based constraints and projected sensitivities from selected future directdetection and accelerator-based experiments. We show this in Fig. 6, in the σ e versus m χ plane (left) and the y versus m χ plane (right), where The accelerator constraints are based on LSND [94], E137 [95], BaBar [28,96], and MiniBooNE [97], and are discussed in detail in e.g. [14,15]. Projections for accelerator-based searches are shown for Belle-2 [28], BDX [98,99], LDMX [30], and MiniBooNE [17]. The direct-detection constraints and projections are the same as in Sec. 3.1. For other projections see [17]. The bounds and projections from accelerator-based searches strengthen for smaller α D : for experiments in which the DM is produced in a beam dump and then scatters in a downstream detector (LSND, E137, BDX, and MiniBooNE), the bound scales as √ α D , while for experiments searching for missing-energy signals (BaBar, Belle-2, LDMX), the bound scales as α D . The direct-detection constraints and projections do not change when varying α D . We see that a significant amount of parameter space is unconstrained by the SN1987A bound, and ripe for exploration by these future searches. The thick orange lines show specific experimental "targets" corresponding to several benchmark models (they do not change when varying α D ). These are based on work from [5,14,15,[100][101][102][103], and we refer the reader to [17] for a summary. Note that we derived the SN1987A constraints for a dark sector consisting only of a Dirac fermion that is coupled to a dark photon (solid lines). Some of the targets shown (in dashed) assume a scalar DM particle, additional interactions within the dark sector, and/or a resonance in the process χ + χ → A * → SM + SM. A scalar φ, even with strong self-interactions (as a SIMP), will have a similar production rate as fermions, and the upper bound will only change by an equivalent number of effective blackbody degrees of freedom, ∼ g φ /(g χ × 7/8). Likewise, a resonance in the DM annihilation cross section, parameterized by R ≡ (m 2 −4m 2 χ )/4m 2 χ , lowers the required couplings to achieve the correct relic abundance, but this resonance does not impact the dark-sector production rate in the proto-neutron star to an appreciable extent. None of these changes to the dark sector content will thus drastically affect the SN1987A constraint, and we find it instructive to show all the "targets" on the same plot. We see that most of the targets are unconstrained by the SN1987A bound; only the resonant thermal targets with R 0.1 are partly constrained.

Inelastic Dark Matter
We now discuss the SN1987A constraints on an inelastic DM model consisting of two states, χ 1 and χ 2 . We will only consider the "light" DM scenario, where the dark photon is heavy and allows for the decay A → χ 1 χ 2 . As discussed in Sec. 2.1.1, we will focus on the case where the elastic, tree-level coupling χ i χ i (i = 1, 2) vanishes. If such a coupling is present, it is velocity suppressed. We thus expect the bounds at small DM masses to be similar to the elastic case discussed in previous sections, but at large DM masses, T c , when the DM does not have much kinetic energy, the bound will likely be similar to the inelastic case discuss in this subsection. Defining ∆ ≡ m 2 − m 1 , there are two cases of interest: (i) ∆ m 1 m 2 and (ii) ∆ m 1 . For case (i), the bounds are essentially the same as the elastic cases discussed in the previous sections. However, for larger ∆, i.e. case (ii), the SN1987A bounds become significantly stronger at large couplings (along the upper boundary) for χ 1 with masses above T c , since it is harder for the DM particles to scatter and become trapped.
Let us discuss case (ii) in more detail. Here χ 1 and χ 2 are produced from (on-shell) A decay and via bremsstrahlung in proton-neutron collisions. However, for sizable ∆, any χ 2 that is produced will quickly decay to χ 1 e + e − through an on-or off-shell A , so that the proto-neutron star essentially contains only χ 1 . In order for the χ 1 to become trapped, they must scatter off protons into the heavier particle χ 2 . This is only possible for those χ 1 that find a proton with energy ∆; the population of such protons is exponentially suppressed if ∆ T c . Therefore, if ∆ 15 MeV, even very large couplings will be excluded by the SN1987A data, since the χ 1 can freely escape.
A simulation is required to calculate the upper boundary accurately for large ∆: the χ 1 scatter into χ 2 , which in turn decay to χ 1 e + e − , with the resulting χ 1 typically having less energy than the original χ 1 . This process can then repeat multiple times as the χ 1 attempt to escape the proto-neutron star. It is computationally challenging to calculate the upper boundary using our trapping criterion, Eq. (2.14), as we did for the elastic case. Instead, we will use a simpler and very conservative criterion: we calculate the couplings needed for which a typical χ 1 scatters off a proton once on its way out of the proto-neutron star. This criterion is appropriate given the other uncertainties and also because after a single scatter a good fraction of the energy is immediately reprocessed into the SM sector via the e + e − produced in the χ 2 decay.
We present our results in Fig. 7. The left plot shows the constraint in the versus m χ 1 plane for various ∆, for the fiducial temperature and density profiles and α D = 0.1. The upper boundary for the ∆ = 0 constraint uses our trapping criterion, Eq. (2.14), while the upper boundaries for ∆ = 0 are derived by requiring the χ 1 to scatter once as discussed above. As expected, the upper boundary of the bounds strengthens dramatically for m 1 m 2 . The SN1987A data constrains very large couplings for large ∆. However, for very large couplings, a two-dark-photon-exchange process at one-loop allows for an elastic scatter of χ 1 to χ 1 , which can dominate over the kinematically suppressed χ 1 → χ 2 transition. We do not calculate in detail this one-loop diagram, but give a simple estimate above which the bounds shown in solid lines in Fig. 7 are not applicable. The cross section for the one-loop diagram is proportional to α 2 α 2 D 4 /16π 2 , while for the tree-level A -exchange process, the cross section scales as αα D 2 . In order to estimate when the one-loop elastic process is important in trapping the χ 1 , we simply set The left-hand side of this equation is set by the value of α D 2 calculated for ∆ = 0 (the elastic case) with our trapping criterion, Eq. (7); setting this equal to α D 2 on the right-hand side then determines when the elastic one-loop scattering process contributes at a similar level. We find that 7 × 10 −3 for α D = 0.1, which is indicated by the dotted line in Fig. 7 (left).
The right plot in Fig. 7, shows the constraints on the y versus m χ 1 parameter space for ∆ = 0.4m χ 1 for our four temperature and density profiles. Here the upper boundary of the SN1987A bound is derived by requiring either the trapping criterion for ∆ = 0 or a single scatter, whichever is stronger. We see that for ∆ T c , the SN1987A data constrains larger couplings, while for smaller ∆, the upper boundary is essentially the same as in the elastic case. We also show current constraints from accelerator-based searches (in gray) and projections from proposed experiments (dotted lines), including Belle-2, MiniBooNE, BDX, and LDMX [40]. We again see that existing constraints, projected searches, and the SN1987A constraints are all complimentary and probe different regions in parameter space.

Millicharged Particles
In this subsection, we consider millicharged particles, as discussed in Sec. 2.1.2. Our SN1987A constraints improve on prior work [31,104] by considering the plasma effects on the SM photon, an improved trapping criterion, and an improved treatment of the high-mass region. We also consider several detailed temperature and density profiles.
Our main results are shown in Fig. 8 (left) in the Q versus m χ parameter space. The solid colored lines show the constraint from using different temperature and density profiles for the proto-neutron star. Note that plasma effects self-consistently cut off the potential divergence at low-momentum transfers in our calculation. The dotted line shows the constraint from [31]. While our lower boundary is slightly higher, our upper boundary is stronger by more than an order of magnitude. We also Other constraints on millicharged particles are taken from [105][106][107]. Right: The constraints on millicharged particles are also applicable to dark matter coupled to an ultralight dark-photon mediator, which we show here in theσ e versus m χ parameter space. We also show a bound on dark-matter-electron scattering using XENON10 data (shaded light blue region) [82], and projections from the upcoming direct-detection experiment SENSEI and a possible experiment using a silicon target sensitive to single electrons with a 1 kg-year exposure (both dotted blue) [15,91].
show constraints from the SLAC millicharge experiment [105], as well as white-dwarf, red-giant, and horizontal-branch stars, all of which are independent of whether the χ is present in the early Universe [106]. In addition, we show constraints from N eff considerations at the time of BBN and the CMB, assuming no dark sector population after reheating [106], and we also show a region in which the DM has not decoupled from the SM at the time of the formation of the CMB [107,108]. The constraints on millicharged particles can also be applied to DM interacting with a massive, but ultralight, mediator. Such a mediator can mediate DM-electron scattering, leading to a cross section that scales as 1/q 4 , where q is the momentum transfer. We show the SN1987A constraints on theσ e versus m χ plane in Fig. 8 (right), together again with the other constraints also shown in the left plot. We also now include a constraint on DM-electron scattering using XENON10 data from [82]. This plot updates the bounds presented in [15]. Projections from selected future direct-detection experiments are shown with dotted lines [15,91]; for other projections see [17].
Note that since our upper boundary is more than an order of magnitude stronger than prior bounds, it could have ramifications for the recently claimed detection of an anomalous absorption strength in the 21cm line from the epoch of first star formation [109,110]; for example, our bounds disfavor some of the parameter space shown to be open in [111].

The Hadronic QCD Axion
We now shift to discuss a different DM candidate, the QCD axion [32][33][34]. Unlike the fermionic DM discussed in Sec. 2, the axion has no conserved quantum numbers. This enables us to make a straightforward calculation of the luminosity as in Eq. (2.7), where we must of course apply suitable substitutions for the bremsstrahlung production rate and the optical depth. In this section, we evaluate this luminosity for the KSVZ or "hadronic" axion. The KSVZ axion couples to the CP-odd combination of gluon and hypercharge field strengths, and also to nucleons with On the equations of motion, the fermion mass may be substituted for the derivative, where m N is the mass of a nucleon. Calculations for the axion bremsstrahlung rate have been obtained previously under a variety of simplifying assumptions. Limits on the axion coupling and mass have been extracted in these contexts starting immediately after the observation of SN1987A. We provide a chronological summary of related prior work in App. E. Here, we evaluate the limit on the axion beyond the diagrammatic calculation of the nuclear scattering cross section and without approximating the luminosity as a blackbody spectrum at large coupling, where absorption is important. Instead, we use results for the spin-flip current obtained at N 3 LO order in chiral perturbation theory [112][113][114] to consistently "correct" the diagrammatic rate. The higher order contributions are stable, but due to a large, wellunderstood destructive interference at NLO they display qualitative differences from the leading order result. In order to make the comparison with previous bounds as clear as possible, we will phrase the N 3 LO results in terms of multiplicative corrections to the leading order result. In practice, we multiply the tree-level result by suitable density-and energy-dependent factors to reproduce the N 3 LO result.
These corrections change existing limits in important ways. At low (high) coupling, our constraints point to a bound on the axion mass that is a factor of about five (one to two orders of magnitude) higher than previously extracted [35]. Equivalently, this implies a bound on the Peccei-Quinn breaking scale that is lower by a factor of about five (one to two orders of magnitude). We discuss the nature of these corrections now.

Corrections to the Axion Bremsstrahlung Rate
Our results incorporate three classes of corrections to the tree-level, massless pion calculation: a cutoff for scattering at arbitrarily low energies, a factor for the nucleon phase space that accounts for the finite pion mass, and a factor that introduces higher orders in the nucleon scattering. These effects have been known in some cases for many years, but they have not been consistently applied to the scattering rate of the axion.
Collecting all effects and setting the notation, we write an amended form of the canonical expression for the axion absorptive width (found, e.g., in [115]) as The factors that appear in Eq. (4.2) are: C i is the coupling of the axion to nucleon i = n, p; Y i is the mass fraction of nucleon i; f a is the axion "decay constant," the scale of breaking of the global U (1) symmetry of which the axion is the pseudo-Nambu-Goldstone boson; σ npπ is the nucleon-nucleon scattering cross section from exchange of a single pion with vanishing pion mass, with canonical value σ npπ = 4α 2 π πT /m 5 N [115][116][117], where α π 15 and T is the temperature of the SM matter in the proto-neutron star; γ f is introduced to cut off the low-energy divergence of Eq. (4.2) [118,119]. We use the form 1 + (n B σ npπ /2ω) 2 −1 [118], which mimics plasma effects that cut off small-angle scattering; γ p accounts for the finite pion mass and nucleon degeneracy, for which we use the dimensionless phase space integral s n B , Y i , ω T , mπ T described in the case of neutron-neutron scattering at arbitrary degeneracy in Eq. (49) of [120]; and γ h is the ratio of the dynamical spin structure function calculated in chiral perturbation theory for nucleons i, j to the value in the one-pion exchange approximation, schematically γ h = S σ / S σ | OPE , for S σ defined in [58,112,113]. The Y i = 0.5 case was originally obtained with a nuclear potential calculation by [119]; the Y i = 0 case was addressed in [58] using the soft radiation approximation and measured nucleon phase shifts; and the extension to arbitrary proton fraction using chiral effective field theory at high densities and measured nucleon phase shift at low densities data was developed in [112][113][114]. For simplicity, we use the fitting function in Eq. 5 and Tab. 1 of [114] and we assume no energy dependence, which is roughly compatible with [58,119] away from the deuteron resonance at relatively low energies.
We reproduce these various correction factors in Fig. 9, fixing ω = T for illustration. Critically, each correction factor individually reduces the original rate by a non-negligible multiplicative factor. Very roughly speaking, we find that the rates are suppressed by a factor between 5 and 100 from the core to the neutrinosphere. We discuss alternate parameterizations of these effects in App. D and find very similar results. Some of these effects, specifically γ f , have been included in calculations of the axion luminosity before, as discussed in App. E, but this is the first attempt to collect all known effects together. Combined with our improved treatment of the energy dependence of the optical depth and our inclusion of different supernova temperature and density profiles, we find that bounds on the axion mass may change significantly from the canonical values. Since we will be interested in the sum of the scattering rate over all nucleon pairs, we define a reduced coupling constant C 2 = Y n C 2 n + Y p C 2 p . We then go on to model-independently bound C 2 along with the axion decay constant f a . Following convention, we will show this as a bound on the axion mass, which is in one-to-one correspondence with the decay constant. The relation between f a . For radii close to the core the suppression is more than two orders of magnitude, so we zoom in on the product of corrections at small radii in the right panel.
and m a is, at leading order, m 2 a f 2 a = m 2 π f 2 π m u m d /(m u + m d ) 2 [121], and including NLO effects the relation is m a = 5.7 eV(10 6 GeV/f a ) [122]. Finally, we have where the density and temperature are of order ρ c = 3 × 10 14 g / cm 3 , T c = 30 MeV, and the reduced coupling in the case of the KSVZ axion, with C n 0, C p −0.47 [122], is C 2 KSVZ 0.066 for Y p = 0.3. We use Eq. (2.7) with the replacements τ = Γ ibr dr → Γ a dr and Γ br → e −ω/T Γ a to get the total instantaneous luminosity in axions. The hadronic axion that we consider is stable against decay and other absorptive processes in the proto-neutron star, so Eq. (4.3) is the only width we need to calculate.
We emphasize that our various correction factors collectively reproduce the N 3 LO calculation in chiral perturbation theory [112][113][114] and together should consistently "correct" the leading order calculation of the axion emission rate. In other words, the product γ f γ p γ h is a self-consistent correction: starting from a simplified calculation for which a closed-form solution is easy to calculate, we wind up with the N 3 LO ChPT result. However, a full calculation should include additional effects and error bars. New nuclear potentials could also be used to expand on our treatment of higher-order corrections, e.g. by including additional energy dependence that we did not model. It is also important to understand more systematically the exact nature of the low-energy cutoff. For these and other reasons, an exact calculation is still desirable. f a · C KSVZ /C [GeV] Figure 10. Left: Luminosity of the QCD axion for a variety of correction factors. The red dashed line labelled "thermal" is the bound one would obtain at large couplings (equivalently, small f a or large m a ) if one assumes that the emission is a blackbody. The black lines instead assume a more accurate calculation as described in the text. The thin black line labelled "uncorr." does not include any correction factors given in Eq. (4.2), while the other black lines include one, two, or all three correction factors, respectively. Right: Luminosity of the QCD axion for a variety of supernova temperature and density profiles.

Results
We plot the luminosity 4 as a function of QCD axion mass times reduced coupling in Fig. 10, and we show the corresponding excluded regions of the axion mass times reduced coupling in Fig. 11. In the left panel of Fig. 10, we show the breakdown of effects arising from the different correction factors γ and also from the novel treatment of the optical depth at high coupling. The improvement in the treatment of the optical depth leads to big effects at large coupling, while the low-energy cutoff and higher-order diagrams have bigger effects at low coupling. In the right panel of Fig. 10, and in Fig. 11, we show the effect of using numerical proto-neutron star temperature and density profiles rather than the "fiducial" profile adapted from [2]. Interestingly, we find that the fiducial profile leads to the most conservative excluded region. In all cases, we are able to close the "hadronic axion window" that had previously existed between the luminosity bounds [115] and the bounds from additional counts in Kamiokande for a more tightly coupled axion [123], labeled "(counts)" in Fig. 11. We also point out  Figure 11. Constraints on the QCD axion mass and axion decay constant for various supernova temperature and density profiles. The "canonical" bound from the PDG [35,115] is shown with a solid gray line, while the bound labelled "counts" comes from [123]. Our bounds close the gap between these constraints, known as the "hadronic axion window". that our revised bound has implications for the claim that stellar cooling anomalies can be explained by weakly coupled, non-hadronic axions [124], and new joint constraints are warranted.
Our results differ from those in canonical references by up to roughly two orders of magnitude at large coupling and a factor of a few at small coupling [35,115]. This comes from several effects, all of which point in the same direction. Our approach to taking into account the energy dependence of the optical depth, following [20], increases the extent of the bounds at high coupling by approximately a factor of five compared to assuming that axions thermalize and are emitted with a blackbody spectrum. The remaining difference between our final bounds and the ones shown in [35] is slightly less than an order of magnitude: the difference is apparent at both high and low coupling, and is attributable to our inclusion of the correction factors γ in Eq. (4.2). The factors γ f and γ p lead to approximately a factor of a few discrepancy with [35], and the corrections to the nucleon scattering rate encapsulated by γ h lead to a similar correction. We illustrate this breakdown in the left panel of Fig. 10. These nuclear corrections have been incorporated for neutrino interactions in various nuclear physics codes that evolve supernova explosions, in particular in [114], but to our knowledge this is the first time these effects have been included in bounds on the interactions of the axion.

Axion-like Particles with Yukawa Couplings
Our analysis of the QCD axion naturally extends to variations on the single-parameter axion model. This allows us to investigate bounds in a general two-parameter space for what is commonly referred Thick solid black, red, green, and red lines show the SN1987A constraints for various temperature and density profiles on axion-like particles (ALPs) with Yukawa couplings, updating the bounds presented in [36]. Other constraints are taken from [36,125], omitting bounds due to Kaon decays from [125,126].
to as an "axion-like particle," or ALP. The ALP mass m A and "decay constant" f A are not related, and we will explore the part of the parameter space for which the finite mass of the ALP becomes non-negligible (our bounds can then be simply extrapolated to lower masses at fixed f A ). As in [36], we consider an axion-like particle for which the mass and coupling are no longer related, m A f A / Λ 2 QCD . The Lagrangian for the ALP scenario is similar to the QCD axion Lagrangian, but the ALP couples to particles other than the nucleons, Since the mass and coupling are now independent parameters, we assume for simplicity that all of the couplings C i are equal, but because of the identity ∂ µf γ µ γ 5 f → 2im ff γ 5 f noted above, the ALP will couple less strongly to lighter SM fermions. For this reason, as well as for the reasons enumerated at the beginning of Sec. 2.3, the ALP production due to scattering or annihilation of electrons in the proto-neutron star is negligible. The luminosity is then given by . For convenience, we have rescaled f A → f A /C i . We explicitly include a contribution to the absorptive width of the axion for its decay to leptons if m a > 2m , although in practice we find that this does not affect the limits at all, since m e T and the Boltzmann suppression effectively depresses the production rate for m a 2m µ . We assume that the coupling to leptons does not affect the early stages of the supernova explosion, but this must be checked for self-consistency. If we omit these couplings the SN1987A bounds are not impacted, but the accelerator and rare-decay bounds from [36,125] are not applicable.
We show our results in Fig. 12, updating the bounds in [36]. In particular, as in the DM case, axions with masses that are kinematically accessible but too weakly coupled to be produced at accelerator-based experiments are potentially probed by SN1987A [36]. Interestingly, a small gap remains between the SN1987A bounds and accelerator-based searches.

Conclusion
In this paper, we have considered constraints derived from the duration of the neutrino cooling phase of SN1987A on two broad classes of DM particles: a dark sector fermion coupled to a kinetically mixed dark photon, as well as the QCD axion and axion-like particles.
For the dark sector fermion, we derive constraints for several different mass hierarchies of the dark-sector particles. We show these constraints for the case of a heavy dark photon that can decay to dark fermions (m = 3m χ ) for elastic DM in Figs. 5 and 6 and for inelastic DM in Fig. 7; for the case of a massive but light dark photon (m m χ T c ) in Figs. 3 and 4; and for the millicharged case (m m χ T c ) in Fig. 8. To derive these constraints, we have suggested a novel criterion for highly mixed dark fermions, wherein they return to chemical equilibrium if they take a random walk in their velocity vector that turns them 90 • from their initial direction of motion. We use this requirement because the scattering cross section for dark sector fermions can be very forward peaked, and scattering an O(1) number of times does not change a light DM trajectory enough to prevent the DM from escaping.
Our bounds have important implications for popular sub-GeV DM models, in which the DM couples to a dark photon of similar mass. They suggest that large regions of otherwise unexplored sub-GeV DM parameter space are now disfavored. However, the SN1987A bounds are complementary to both existing bounds and proposed experimental searches: they lie well below current bounds, and many motivated and concrete benchmark-model "targets" remain unconstrained. This further emphasizes the need for a robust experimental program to search for sub-GeV DM as envisioned in [17], at least down to the SN1987A constraint, if not beyond.
The QCD axion has been studied in some detail previously, but bounds on axion properties from SN1987A have heretofore been extracted with a range of simplifying assumptions that are known to be violated at the order-of-magnitude level. Here we attempted to rectify this situation by including some estimates of known nuclear physics and particle physics effects. In particular, recent progress in chiral effective theories demonstrates that corrections up to N 3 LO can have a substantial impact on the spin fluctuation rate of free nucleons [112][113][114], confirming earlier calculations using nuclear phase shift data [58,119] that have long been applied to the neutrino emissivity. These effects conspire to point in the same direction, resulting in large changes to the expected axion emission rate. Coupled with our improved description of boson luminosity in the high-mixing limit, the axion bounds are changed significantly from the "canonical" range, as shown in Fig. 11. We also re-visited the constraints on axion-like particles with Yukawa couplings, shown in Fig. 12, finding some difference with the previous literature.
The wealth of information that has been gained over the years from the observation of SN1987A is rather remarkable. As simulations of core-collapse supernova keep improving, it will be highly desirable to continue the effort to include new, weakly-coupled particles directly into the simulations. from the YITP when this work commenced. For the second half of this work, SDM was supported by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

A Production and Decay of Dark Photons
In the relativistic, degenerate regime we use Eq. (77) of [43] with the conventions of [18] to define the SM photon polarization tensor: As a result of the structure of the mixing angle, there is a particular energy ω * at which Re Π = m 2 , where the mixing angle hits a resonance. When production is resonant, the magnitude of the differential power exactly compensates the narrow width of the resonance, and the luminosity becomes independent of the production mechanism. The rates for resonant inverse bremsstrahlung and electromagnetic decay of the A particle are [20] Γ L,T ibr. = 32 3π where σ (2) is an angle-and energy-averaged neutronproton scattering cross section extracted from measured nuclear phase shifts [26]; we introduce a kine- (1 − 4m 2 e /m 2 ) (1 − m 2 /ω 2 ) for the outgoing electron-positron pair. In the soft radiation approximation, the detailed balance factor e −ω/T between bremsstrahlung and inverse bremsstrahlung becomes unity, and so we define a bremsstrahlung production rate np (T ) only in the lower limit of the energy integral.

B Dark Matter Bremsstrahlung
One of the dominant production modes for DM in the supernova is via on-or off-shell A bremsstrahlung during nucleon elastic scattering events, as in Fig. 2. We calculate this amplitude of this process in the limit of soft bremsstrahlung following §6.1 of [127].
The matrix element for DM production is where M np is the matrix element for the process with no bremsstrahlung, which is n − p scattering [26]; the incoming p [n] has four momentum P µ 1 = (E 1 , p 1 ) [P µ 2 = (E 2 , p 2 )]; the outgoing p [n] has four momentum P µ 3 = (E 3 , p 3 ) [P µ 4 = (E 4 , p 4 )]; the DM particles have four momenta P µ χ = (E χ , p χ ) and P μ χ = (Eχ, pχ); and the dark photon carries an interior momentum K µ = P µ χ + P μ χ = (ω, k). In what follows, lower-case letters without the vector symbol indicate the magnitude of the three vector, e.g. k = | k|. We include different contributions from the longitudinal and transverse modes, which can contribute differently in the dense environment of the proto-neutron star.
In the low momentum or "soft" limit, Eq. (B.1) becomes

(B.2)
We square the amplitude and take the trace. Current conservation, K µ P µν = 0, implies P µν (P ν χ P β χ + P ν χ P β χ )P αβ = −2P µν P ν χ P β χ P αβ , leading to We find that we may in general rearrange this as where, for the sake of brevity, we separate the contribution due to the n − p dynamics from a "soft factor" S due to the DM emission. This soft factor is a function only of the virtual and radiated particle momenta. Assuming that the DM does not scatter on its way out of the star, we can calculate the local differential luminosity per unit volume, where the f i are distribution functions. In the following, we will assume the particles are nondegenerate such that we may ignore all (1−f ) factors. We also approximate the effect of the soft radiation limit (invoked above to obtain the matrix element) by taking δ 4 (P 1 +P 2 −P 3 −P 4 −P χ −Pχ) δ 4 (P 1 + P 2 − P 3 − P 4 )e −ω/T . Note we use a different approximation from [20], which gives more conservative results. Changing variables Pχ → K for convenience, we have Since the first line only depends on the nucleon scattering and the second line is only sensitive to the DM kinematics, we can calculate them separately. The first line of Eq. (B.6) involves many of the same features as the result in [20] and summarized in App. A, and we follow a similar procedure. In particular, we assume the nucleons are nonrelativistic and invoke the relations where σ which is calculable numerically.

C Dark Matter Elastic Scattering
Once DM particles are produced, they elastically scatter off protons on their way out of the supernova, as in Fig. 2. This can lead to the thermalization of the DM particles, which can allow them to return to chemical equilibrium, as described in Sec. 2.5. The matrix element squared for this process is where P 1 (P 2 ), P 3 (P 4 ) are incoming and outgoing DM (proton) momenta and K = P 1 − P 3 is the momentum transfer. We define a scattering rate and an average angular deflection per scatter by (2π) 4 δ 4 (P 1 + P 2 − P 3 − P 4 )θ 13 |M| 2 s , (C.3) where θ 13 is an angle between incoming and outgoing DM, and we assume the protons and the DM are nondegenerate, such that 1 − f 3 1 − f 4 1. Also, we assume f 2 follows the Maxwell-Boltzmann distribution. With these definitions, we can define the total number of scatters and the maximum angular deflection in Eq. (2.13). Scattering through a light mediator has a t-channel singularity and thus is weighted towards small angles, θ 13 | typical π/ few. This characteristically "forward peaked" scattering indicates that light DM scattering through a light mediator neither loses a significant fraction of its energy nor deviates far from its initial trajectory in a typical scattering event. Thus, measures of DM decoupling that assume stationary initial state nucleons and calculate mixing angles for which the DM scatters an order one number of times are bound to overestimate the tendency of DM to be trapped in the neutron star and thus to predict an overly small value of the kinetic mixing as the trapping line.
We also point out that if the dark photon is massless, as we assume in Sec. 3.4, then the longitudinal mode decouples and the propagators P L should be omitted from Eq. (C.1) because the different polarizations do not mix.
Finally, note that in Eq. (2.13), we numerically evaluate the integrals from R d to R ν , and we approximate the integrand for the region from R ν and R f as Γ(R ν )R ν /5v.

D Alternate Parameterizations of Axion Corrections
In Sec. 4.1 we chose to correct the rate for axion bremsstrahlung to match the energy-averaged rate calculated at N 3 LO order in chiral perturbation theory. Other parameterizations of similar effects exist in the literature, and the exact rate could in principle deviate from the chiral perturbation theory result. Here, we summarize some possible alternative correction factors that account for similar physical effects in different ways: γ mπ accounts for the finite pion mass rather than γ p , which we model as γ mπ = 1 + m 2 π 3m N T −2 , roughly matching [117] (this analytic prescription falls between the numerical work of [128], obtained with non-degenerate nucleons, and the result of [129,130], calculated for degenerate nucleons); γ SRA is the ratio of the dynamical spin structure function for nucleons i, j in the soft radiation approximation to the value in the one-pion exchange approximation with finite pion mass, obtained numerically with the aid of nuclear phase shift measurements. The Y i = 0.5 case was originally discussed by [119] while the Y i = 0 case was addressed in [58], each finding reductions of order a few. For the purpose of illustration, we will neglect the density dependence and simply assume a constant factor of 5 reduction in the rate compared to the uncorrected result, which potentially underestimates the axion luminosity; and γ LPM accounts for the Landau-Pomeranchuk-Migdal effect, for which we use the semi-analytic fit  Together, γ SRA and γ LPM should roughly combine to account for the same physics as γ f and γ h in our main results. Likewise, the factor γ mπ potentially mimics the effect of the pion propagator in place of γ p . We show results for the luminosity in axions using these alternate correction factors in Fig. 13. This figure is the counterpart of Fig. 10. By comparison of these two figures we see that the excluded regions change very little regardless of how we choose to account for these nuclear physics corrections.

E Summary of Previous Work on the Hadronic Axion
The absorption rate of the QCD axion has been obtained to varying degrees of precision since before the explosion of SN1987A. Because the axion is predominantly produced during nucleon-nucleon bremsstrahlung, the exact result requires evaluating the fifteen-dimensional integral of a non-perturbative matrix element with a partially degenerate phase space. This technical challenge has taken quite some time to thoroughly understand.
Here we summarize the evolution of the work that has previously put bounds on the QCD axion, listed in chronological order: [1,132] provided the first calculations for axions emitted from nucleon-nucleon bremsstrahlung, modeling the nuclear interaction with a single (massless) pion exchange and assuming that the squared matrix element is constant in the nucleon momenta; [133] compared measured pion production rates in p−p scattering to those found from a diagrammatic one-pion exchange calculation and found that these agreed to within a factor of a few; [134] conducted supernova explosion simulations including a free-streaming axion energy sink and backreaction on the star for a wide variety of proto-neutron star profiles; [116] computed phase space integrals over the one-pion exchange diagram for arbitrary nucleon degeneracies, justifying the use of non-degenerate phase space; [135] conducted supernova explosion simulations for a tightly coupled axion, confirming prior bounds in the trapping regime; [117] verified the calculation of [116] and discusses ways of cutting off pathological limits, including the first appearance of the 1/(ω 2 + aΓ 2 ) prescription and a discussion of when the pion mass should not be neglected; [118] also advocates the 1/(ω 2 + aΓ 2 ) approach and additionally proposes a "saturation width" that cuts off the rates at some maximum spin fluctuation rate.
All of these authors roughly agree in the free-streaming limit, finding the requirement that the Peccei-Quinn scale must respect f a 10 8−9 GeV, with the uncertainty on this limit primarily arising from the difference in treatment of the low-energy scattering, which can be cut off by the 1/(ω 2 + aΓ 2 ) factor. Many other works have calculated neutrino couplings and luminosities, which are important because neutrinos and axions couple to the same nuclear current. There are too many developments to name here, but we do clarify the origin of the chiral effective theory corrections that we utilize above: [119] gives the spin density structure function for n − p scattering in a variety of ways, indicating a qualitative difference in the magnitude of the scattering rate and in the density dependence, ultimately due to the different contribution to the partition function of n − p scattering, which can be resonant near the formation of a deuteron; [58] gives the ratio of spin density structure function for identical nucleon scattering based on measured phase shifts; [112][113][114]136] use a chiral effective theory approach at high densities and show that this matches to the phase shift analyses at intermediate densities, all of which confirm the high-density suppression and low-density enhancement suggested by [58,119] These corrections are a major ingredient that lead us to the modified limits shown in Fig. 10.