Supernova signals of light dark matter

Dark matter direct detection experiments have poor sensitivity to a galactic population of dark matter with mass below the GeV scale. However, such dark matter can be produced copiously in supernovae. Since this thermally-produced population is much hotter than the galactic dark matter, it can be observed with direct detection experiments. In this paper, we focus on a dark sector with fermion dark matter and a heavy dark photon as a specific example. We first extend existing supernova cooling constraints on this model to the regime of strong coupling where the dark matter becomes diffusively trapped in the supernova. Then, using the fact that even outside these cooling constraints the diffuse galactic flux of these dark sector particles can still be large, we show that this flux is detectable in direct detection experiments such as current and next-generation liquid xenon detectors. As a result, due to supernova production, light dark matter has the potential to be discovered over many orders of magnitude of mass and coupling.


Conclusion 21
A Analytic profile of SN

Introduction
The particle nature of dark matter (DM) remains one of the largest outstanding puzzles in physics. Despite the overwhelming evidence for the existence of dark matter from its gravitational imprints on cosmological and astrophysical scales, there have as of yet been no observations of any non-gravitational interactions [1]. The fact that our current measurements leave an enormous range of possibilities for its mass and interactions with the Standard Model has motivated a very rich experimental program exploring a wide range of dark matter models. A large ongoing experimental effort is searching for dark matter candidates with masses in the GeV−TeV range, largely motivated by the WIMP miracle (see e.g. [2]). These experiments have achieved incredible sensitivity to dark matter by searching for very small energy depositions from dark matter scattering with nuclei in extremely clean environments. Despite their great progress, such experiments quickly lose sensitivity to lighter dark matter candidates because the kinetic energy of such candidates is too small to lead to observable signatures in the detectors. This is true under the assumption that we can only detect dark matter particles which are gravitationally bound to our galaxy, which implies a maximum velocity for the dark matter flux.
The existence of astrophysical sources where dark matter could be produced with larger velocities and at a non-negligible flux allows one to significantly extend the reach of these experiments to models of sub-GeV dark matter. 1 Core-collapse supernovae (SN) can reach core temperatures in excess of 30 MeV for O(10) seconds, allowing them to produce vast thermal fluxes of particles with masses O(100) MeV at relativistic speeds. This makes them an ideal astrophysical source for sub-GeV dark matter. Supernovae have already been used extensively to constrain a plethora of models of new physics. In almost all cases, the criterion applied in order to place a bound is the so-called cooling criterion, which states that if any new particle were able to transport energy out of the protoneutron star formed by the SN more quickly than the neutrinos, the cooling timescale of the core would be less than the ten-second timescale observed in SN1987a. This is equivalent to the statement that any new particle that transports greater than 3 × 10 52 erg/s through the radius at which the neutrinos are no longer diffusively trapped is incompatible with observations [7].
There are two distinct regimes relevant for these cooling bounds. At lower couplings one considers bulk emission from the entire protoneutron star volume, which results in a lower bound on the couplings below which the luminosity is less than 3 × 10 52 erg/s due to insufficient production of the new particles. There is a separate regime at large couplings, in which the coupling is large enough that the new particles are diffusively trapped inside the core and the emission effectively occurs from a radius at which the densities and temperatures are low enough to allow the new particles to escape freely. This is the regime in which the upper bound on the couplings constrained by cooling can be derived. The trapped regime is reasonably well-understood analytically for particles that are singly emitted or absorbed (e.g., axions and dark photons). However, for particles that can only be pair-produced, a detailed understanding of different processes (e.g. production versus scattering) is required, making analytic estimates more challenging. For examples of existing cooling constraints on pair-produced particles, see e.g. Refs. [8][9][10][11][12][13][14][15][16].

Summary
In this paper, we move beyond these cooling arguments by considering the direct detection of a hot population of supernova-produced dark matter. Even in parameter space outside the cooling bound, a supernova can still produce a vast flux of light dark matter particles. This flux is also hot (semirelativistic), which allows for the possibility of its detection in current and next generation WIMP experiments. Using a Monte Carlo Boltzmann particle transport simulation, we are able to compute the DM flux in the trapped regime, which allows us to estimate the reach of direct detection experiments that were originally expected to only have sensitivity to dark matter with masses above ∼ GeV. Furthermore, our results improve upon and extend previous cooling constraints in the trapped regime which relied on a number of approximations in order to produce an analytic estimate of the flux.
For concreteness we focus on a simple model of dark matter, in which it is a Dirac fermion which interacts with the SM via the four-fermion operator e g d Λ 2χ γ µ χJ µ em (2.1) where χ is the DM field and J µ em is the electromagnetic current of the SM. Such an interaction can be generated if, for example, DM is charged under a dark gauge boson with mass m A = Λ and DM-DM coupling g d that kinetically mixes with the SM with mixing parameter . 2 While this model is a good first test case, it should be noted that the main point of this paper (that SN can produce a flux of light dark matter that is observable in direct detection experiments) also applies to a much wider variety of models. We leave those for future work.
The high temperatures reached in core-collapse SN allow them to produce large abundances of the sub-GeV dark sector fermions considered in this paper. In the regions of parameter space we are interested in, these fermions have a sufficiently strong coupling to the Standard Model that they become diffusively trapped near the protoneutron star that forms from the SN core. The diffusive trapping is primarily due to scatterings off of the free protons generated by the dissociation of nuclei in the SN shock. These scattering interactions are inefficient at changing the dark fermion's energy because of the large mass ratio between the DM and the nuclei. The dark fermions also scatter off of electrons and positrons, allowing thermal exchange with the SM bath. At a certain radius (which we call the electron sphere or energy sphere) the density of electrons and positrons drops to the 2 If the new particle accounts for all of dark matter, there are stringent CMB bounds from DM late annihilation that place strong constraints on this model [17]. Those can be evaded by considering either asymmetric DM or by introducing a small mass splitting, making it a pseudo-Dirac fermion (see e.g. [18,19]). Both of those scenarios would not affect any of our conclusions (as long as the mass splitting is small compared to the temperature of the SN) and so for simplicity, we choose to focus on the simple Dirac fermion model in this paper. point where a dark fermion no longer remains in thermal contact with the SM, see Figure 1. Additionally, there is a streaming sphere at which the density of protons drops low enough that dark fermions are no longer diffusively trapped and begin to free-stream out of the star. Any given dark fermion produced in the SN will therefore diffuse through a protonrich overburden until it either reaches the streaming sphere and escapes or encounters an antiparticle and annihilates.
The dark fermions that do eventually escape are produced with a distribution of semirelativistic velocities. This results in a time-spreading effect during their propagation to Earth. The difference in arrival time of the high-momentum and low-momentum ends of the spectrum is of order the dark fermion travel time between Earth and the SN, hence the dark fermions produced by a single SN arrive on Earth over a timescale of 10 5 years for an average galactic SN. This is in direct contrast to the neutrinos, which are all produced highly relativistically and therefore arrive as a single pulse over a timescale of ten seconds.
Given the typical rates and distances of galactic supernovae, in addition to the inherent signal spread, the dark fermion fluxes from various SN should overlap in time, producing a diffuse galactic SN flux of dark fermions. This signal is reminiscent of the diffuse flux of SN neutrinos, with the distinction that diffuse neutrinos would arise from the sum of a much greater number of extragalactic sources each of short duration, since the time between galactic SN is much larger than the duration of the neutrino flux.
If a diffuse dark matter SN flux is continuously passing through Earth, we must consider ways to detect it with Earth-based experiments. We find that the diffuse flux of DM is detectable in existing and next-generation liquid xenon (LXe) WIMP detectors. Interestingly, though the idea is similar to the direct detection of the diffuse SN neutrino flux, it is not the large neutrino detectors which are best to search for this flux but rather the WIMP detectors due to their low energy thresholds. Though the WIMP detectors were designed to hunt for dark matter on the GeV scale, we show that they are sensitive to recoils by sub-GeV dark sector fermions over a wide range of masses and couplings above even the newly-computed trapped regime cooling bound. As this idea probes extremely weak couplings, it is complementary to most of the experimental proposals searching for sub-GeV dark matter through direct detection or in accelerators (see [20,21] and references therein for details of some of the other proposals for detecting sub-GeV DM).
Existing LXe experiments such as Xenon 1T [22] are already sensitive to the diffuse galactic SN flux of dark sector fermions, and future experiments such as Xenon nT [23], PandaX-4T [24], LUX-Zeplin [25] and, on a longer time scale, DARWIN [26] will cover an even larger region of parameter space.
In Sections 3 and 4, we describe an analytic treatment of the required computation and explain the details of the Monte Carlo Boltzmann transport simulation. We discuss our computation of new cooling bounds in Sec. 5 and direct detection by LXe detectors in Sec. 6. Our results are presented in Section 7.

Analytic approximation
While the final results used in this paper were computed using a numerical Monte Carlo simulation of particle transport within the supernova and cooling protoneutron star, we first provide a simple physical picture of the expected behavior of the DM flux in the diffusive regime. We then apply the intuitive description to demonstrate how to make rough analytical estimates of the spectrum of the SN-produced dark fermion flux.
Our basic premise is to generalize the idea of a "neutrino sphere" to the case of dark fermions. The term "neutrino sphere" is typically used to describe the radius at which the density of nucleons has dropped such that the neutrinos are no longer diffusively trapped. It is common to approximate the SN neutrino flux as simply the emission of a blackbody sphere with radius and temperature given at the neutrino sphere.
While neutrino decoupling is often assumed to occur at a single radius, this is not a very good approximation for what happens inside a supernovae, specially for 2nd and 3rd generation neutrinos. The number-changing interactions for those neutrinos will freeze out at a different radius than do their scattering interactions with nucleons. Because neutrino interaction cross-sections have a sharp energy dependence, it turns out that the number-changing and energy-changing neutrino interactions freeze out at roughly the same radius [27], but significantly earlier than when they can start free-streaming. This allows their flux to be estimated with reasonable accuracy with a simple combination of a blackbody approximation plus a transmission calculation as described in Ref. [27].
In the case of the dark fermions, the interactions that set the overall number density, the energy spectrum, and the free-streaming radius are all different, so the single sphere approximation is not valid. Instead, we can break the protoneutron star into three radii we have termed characteristic spheres at which different interactions freeze out and cease to affect the dark fermion flux. They are as follows: 1. Annihilation sphere (r N ): This is the radius at which χχ → e + e − freezes out, or, in other words, the DM density has dropped sufficiently that the dark fermions are no longer annihilating with their antiparticles. There are effectively no numberchanging reactions outside this sphere, hence it is this radius that sets the number flux of escaping DM.
2. Electron sphere/energy sphere (r E ): This is the radius at which χe → χe freezes out. Beyond this radius, scattering events of dark fermions with electrons and positrons are no longer sufficient to keep the DM in thermal contact with the SM bath. When r E > r N , this sets the temperature of the escaping DM flux.
3. Streaming sphere (r S ): This is the radius at which χp → χp freezes out. The proton density drops to a point that the DM is no longer diffusively trapped and the DM free-streams out of the star. Note that because the protons are significantly heavier than the DM, they cannot efficiently transfer energy to the DM, hence scattering interactions with protons do not change the energy of the DM appreciably.  Figure 1. The characteristic spheres of the protoneutron star. Outside the annihilation sphere, number-changing processes for the DM (pair production and bremsstrahlung) freeze out, setting the number flux. Outside the energy sphere, the DM thermally decouples and its spectrum is set. Outside the streaming sphere, the DM is no longer diffusively trapped by proton scattering and free-streams out of the star.
In the parameter space of interest, the streaming sphere always lies well outside of the annihilation and electron spheres, hence the number flux and energy distribution are set while the DM is still diffusing. We also find that the electron sphere is always outside of the number sphere (r N < r E ). As we have already discussed, even though the DM continues to scatter off of protons once outside of the electron sphere, the large discrepancy in mass between protons and the dark fermions means that the energy of the dark fermions is not largely affected during these scatterings. As a result, the energy spectrum of the DM flux is set by the temperature at r E . Due to this, we will use the terms electron sphere and energy sphere interchangeably. The characteristic spheres are depicted in Fig. 1.
We can analytically compute these characteristic spheres by finding the radius at which the optical depth associated with a particular interaction becomes O(1) [27]. The optical depth for a given process at some radius r 0 is given by ∞ r 0 λ −1 (r) dr with λ(r) the interaction length of the process as a function of radius. The interaction lengths for χχ annihilation, χe ± scattering, and χp scattering are as follows: with n X the number density for a species X, σ Y the cross-section for a process Y , and T the temperature of the SN at the given radius. The explicit forms of the cross-sections are provided in Appendix B. Note that there is an additional factor of the averaged velocity v χ in the mean free path for interactions with electrons. This comes from the fact that the velocity factor in the cross-section definition is dominated by the relativistic electrons for this reaction, but the mean free path of the DM particle depends only on its own velocity. For the other reactions the velocity factors cancel out since the velocity term in the cross-section is dominated by the DM velocity. 3 Because the mean free-path for scattering with protons is much shorter than that of the other interactions, there will be many scatterings with protons in between scatterings with electrons or annihilations. This must be taken into account by using an effective optical depth for the scattering with electrons and annihilations as discussed in Ref. [27]. The characteristic radii r N , r E , and r S are computed using the following criteria: Having computed these radii, one can use a similar argument to that of the neutrino sphere to make a simple estimate of the outgoing DM spectrum. The logic behind the following methodology is simply that, by definition, r N sets the total number of dark fermions that are produced (since number-changing reactions are insignificant beyond it) and r E sets the energy spectrum (since the DM is thermally decoupled beyond it). Under the assumption that the dark matter flux has reached a temporary steady-state so the DM density profile is not changing in time, the flux that escapes at infinity will be set solely by the total number flux produced and the temperature at the energy sphere.
The analytic estimate proceeds as follows: 1. Treat the protoneutron star as a blackbody of radius r N with a diffusive envelope. The number flux at the blackbody surface is given by where g χ = 4 is the number of degrees of freedom in DM. To obtain an energy flux one can just multiply the integrand by the DM energy.
2. Multiply this total flux by a normalized Fermi-Dirac distribution f F D (T r E ) for the temperature at the energy sphere: 3. Even though the number changing reactions are frozen out at r > r N , some particles emitted from that radius can bounce back and return to the region r < r N as they are trying to diffuse out of the streaming sphere. Therefore one must include a transmission factor to account for the losses due to this effect (see Ref. [27] for details): where τ r N is given by This approach involves a number of simplifying approximations, such as the notion of a sharp radial freeze-out for different processes, but despite its limitations it provides a simple physical picture of how the various interactions affect the DM flux. Furthermore, it serves as a cross check on the results of the full Monte Carlo simulation. We used it as such and found that the analytic estimates agreed with simulation results to within an order of magnitude. To provide a point of comparison, we include at the end of Sec. 4.3 a comparison of the DM profile generated by our analytic estimate to the output of the Monte Carlo transport simulation.

Flux computations from simulation
While the methodology outlined in the previous section provides a convenient way to understand the physics of the trapped regime, it becomes less accurate as the region over which the decoupling of a certain process occurs becomes larger (i.e. the decoupling radii become smeared into decoupling regions). To produce more robust estimates of the DM flux, one must perform a full Boltzmann particle transport simulation. To calculate the transport of dark fermions we use a Monte Carlo (MC) method that can be broken into four main steps: 1. Initial Conditions: To describe the underlying matter distribution of a protoneutron star we refer to detailed multi-physics dynamical simulations of core-collapse supernovae. From these we define fiducial analytic profiles that capture the temperature, density, and electron fraction structure of the cooling protoneutron star.
2. MC flux computation: An iterative MC simulation is used to determine the steadystate DM distribution within the star. The resulting DM profile is used to determine the outgoing number flux of dark fermions.

Energy spectrum:
The energy distribution is set in the same manner as the analytic treatment described in the previous section. The energy sphere is computed and used to set a temperature for the outgoing flux.

Gravitational redshift:
This energy spectrum is subsequently adjusted to take into account the effects of gravitational redshift on the escaping DM.
In the following subsections, we will address each of these steps in turn.

Initial Conditions
The collapse of a massive star in a core-collapse supernova leads to the compression of the inner iron core into a dense proto-neutron star (PNS) with a mass M ≈ 1.4 M and a total thermal energy (derived from the gravitational collapse) of order E ≈ 10 53 erg. At its formation, the PNS has a radius of several tens of kilometers, although over a timescale of tens of seconds it will cool and and condense into a cool neutron star of radius ∼ 14 km. The mechanisms of core-collapse SNe and PNS formation are complex multi-physics phenomena that involve a dense matter equation of state and multi-dimensional effects such as turbulence and convection. While constructing high-fidelity simulations of these events remains a work in progress, the essential structure of the remnant PNS can be specified at a level appropriate for making the DM flux estimates of this paper.
The gray lines in Figure 2 show an example PNS structure about 1 second after bounce. In the bulk interior of the PNS, the energy density is dominated by baryons, leading to temperatures of order T ≈ 30 MeV. Deleptonization via weak interaction results in a low electron fraction Y e ≈ 0.1. The very central core (r 5 km) of the PNS typically has a slightly higher Y e ≈ 0.2 and is factor of ∼ 2 colder, as this region is adiabatically compressed in the collapse without experiencing strong heating from shocks.
Above the PNS is a steep, hydrostatic atmosphere where the density falls off exponentially with a scale height 1 km. The temperature structure in the atmosphere is roughly set by neutrino diffusion, which implies a scaling T (r) ∝ ρ(r) 1/4 . The neutrino sphere generally sits somewhere within this steep PNS atmosphere. The electron fraction rises, approaching symmetry (Y e ≈ 0.5) outside the PNS.
Finally, above the PNS atmosphere the densities drop such that radiation pressure dominates and the profile changes. The layers above the PNS are initially convective and hence nearly isentropic, which results in the density and pressure having approximately power-law profiles (ρ(r) ∝ r −3 , T (r) ∝ r −1 ). Matter driven by neutrino winds from the PNS surface may also influence the structure above the atmosphere, but this still results in a similar power-law profile.
To capture these essential features of PNS structure without restricting ourselves to any specific core collapse simulation, we constructed an analytic profile that resembles the results of full simulations. Details of the analytic mapping are given in Appendix A and a comparison of simulation data to the chosen fiducial profile is shown in Figure 2.
We estimated the profile dependence of our results by varying the parameters of our analytic profile. The resulting flux is most sensitive to the overall scale of the temperature since the production terms depend strongly on temperature (see Appendix B.5). Rescaling the profile such that the peak temperature changes from ∼ 30 MeV to ∼ 50 MeV results in an increase in the flux by a small O(1) factor for masses below ∼ 40 MeV and an order of magnitude for large masses. This is unsurprising given that the larger masses are already being produced on a Boltzmann tail, so the production is exponentially-sensitive to the temperature. However, even order-of-magnitude changes in flux make no appreciable change to the sensitivity bounds displayed in this paper due to the fact that the flux changes very rapidly with y. It is true that with a higher temperature, the bounds may extend out  Figure 2. The analytic profile used in this analysis (colored) is displayed alongside the results of one run of the supernova core-collapse simulation (dashed). Note the strong agreement between the analytic profiles and the simulation results. The temperature in the analytic profile is uniformly lower than simulation because it has been adjusted such that it reaches a maximum temperature of 30 MeV, which is a conservative and theoretically-motivated peak core temperature (see discussion in text).
to slightly larger masses, but we have selected the coolest physically-motivated profile and therefore treat our bounds as conservative. Using our analytic profiles for temperature, density, and electron fraction, it is straightforward to compute the resulting abundances of all SM particle species. To compute the proton number density, we assume that the electron and proton fractions are comparable in the protoneutron star (i.e. Y (r) ∼ Y p (r)) and that the protons are the dominant contribution to the total mass density. These assumptions immediately yield n p (r) = Y (r)ρ(r) mp as the proton number density.
To compute the thermal densities of the electrons and positrons, we make the assump-tion of thermal equilibrium and use the associated thermal abundances (see [28]).
The chemical potential µ e can be determined by enforcing charge neutrality, which requires the number density of electrons be equal to the sum of the proton and positron number densities. This yields the following the condition: Critically, these profiles are assumed to be unchanged on the timescale of emission (∼ 10 seconds), hence they are maintained as a fixed background in the following step of the analysis: the MC simulation of the dark fermions.

MC flux computation
Having now found the radial profiles of the SM species, we must determine the DM profile. This necessitates the use of a Monte Carlo simulation of dark fermion diffusion within the protoneutron star.
We begin by computing source and annihilation terms to be used as inputs to the simulation that dictate the DM emissivity and annihilation length, the details of which appear in Appendix B. The two primary interactions that source dark fermions are electronpositron annihilation to DM (e + e − ↔χχ) and proton-neutron bremsstrahlung (np ↔ npχχ). We include both contributions, but we find that for all DM masses considered, the production from electron-positron annihilation dominates over bremsstrahlung except in the innermost region of the core (r 5 km), which coincides with the rise in temperature in the profile we considered.
Within the protoneutron star we represent the dark matter fermion field by a set of N discrete tracer "packets," each of which represents a number of fermions. The initial location and energy of these DM packets is sampled randomly so as to match the total thermal DM emissivity at each location in the protoneutron star. The DM packet are propagated a distance d before experiencing a matter interaction event, where d is determined in standard MC fashion by d = −λ ln(z) where λ is the total mean-free path and z is a uniform random number between (0, 1]. If the interaction event is a scattering, the direction of the DM packet is resampled from an isotropic distribution. If the interaction event is an annihilation, the DM packet is removed from the calculation. DM particles that leave the edge of the domain are tallied as escaped. The inclusion of self-annihilation induces a non-linearity in the transport problem due to the fact that the DM annihilation depends on the background density field of other DM particles. We address this with an iterative approach. Initially, the DM density in each zone is assumed to be thermal at the local temperature. We then run the MC transport procedure and construct an improved DM density profile by counting the DM packets passing through each zone. The entire transport step is repeated using the newly we also display a purely thermal profile in red. Our analytic estimate of the profile is fixed to be thermal up to some decoupling radius, at which point, it free-streams with r −2 . This free-streaming is shown as a blue line. There is an O(1) discrepancy between this analytic decoupled profile and the simulation results due to the fact that the analytic profile assumes instantaneous decoupling, but the scaling behavior at large radii of the two profiles is the same.
constructed DM density profile, and this process is iterated until the density structure and emergent DM flux converges. Typically we find that ≈ 20 iterations is sufficient to converge to a self-consistent DM distribution. For simplicity, the annihilation crosssections are assumed to be angle-and energy-independent, although such effects would be straightforward to include. We include in Figure 3 a comparison between the analytic estimate of the DM profile described in the previous section and the results of the simulation. The simulation results are displayed in green, a purely thermal DM profile is shown in red, and a blue line shows the free-streaming behavior of our analytic profile beyond r S . The profile described in the previous section is fixed to be thermal up until r S and then falls with r −2 beyond it. It is clear that while there is an O(1) difference between this analytic estimate of decoupling and the simulation mainly due to the treatment of the decoupling as occurring at a single radius instead of over an entire region, the general features and scaling behavior at large radii are the same. Part of the discrepancy between the analytical estimate and the Monte Carlo result can be associated with the transmission factor in Eq. 3.10, which was not taken into account in the figure since it only applies to the asymptotic flux.

Energy spectrum
With the number flux computed from the MC simulation, we must set the energy spectrum for the escaping DM. While it would in principle be possible to extract a complete spectrum from the MC simulation itself, we find that only the dark fermions living in the high-momentum tail of the spectrum will be observable in liquid xenon detectors. Comput-ing this tail with any precision is computationally prohibitive in that it would require the simulation to track a vast number of dark fermions such that the tail would not be dominated by statistical noise. Therefore, we instead choose to employ the analytic method detailed here to compute the spectrum because it allows for a robust prediction of the quantity of the escaping flux living in the high-momentum tail of the spectrum.
We compute the spectrum in the same manner as in the analytic methodology outlined in the previous section. Namely, we compute r N and r E using Eqs. 3.4 -3.6, with number densities for protons, electrons, and positrons set by the abundances computed in Section 4.1. Note that the cross-sections that appear in the interaction lengths are momentum-dependent. For these computations, the momentum is taken to be the average center-of-mass momentum at a given radius. This is simply p CM = 3T (r) for DM scattering off of electrons/positrons and p CM = 6m χ T (r) for DM scattering off of protons.
As before, we take the temperature at thermal decoupling to be T (r E ). We then enforce that the DM energy spectrum take the form of a Fermi-Dirac distribution at this temperature, but with normalization set by the number flux determined via the MC simulation. Hence, we have the following differential flux: where Φ χ denotes the total DM flux in number per second and Φ MC χ denotes the number flux computed with the simulation.

Gravitational redshift
Finally, we must take into account the effect of gravitational redshift on the spectrum computed in the previous step. The redshifted momentum of a DM particle emitted with p 0 at r E is given by with ∆Φ the change in potential between r E and r = ∞, defined as where m enc (r) is the mass enclosed within r.
In the region of parameter space we are interested in, this effect does not decrease the momenta of escaping dark fermions by more than an O(1) factor. However, the effect does introduce a sharp cutoff in the spectrum corresponding to where the DM no longer has sufficient initial momentum to escape the gravitational well. This cutoff momentum is given by We find that including these effects decreases the DM flux above detector threshold by ∼ 30 − 40%.

Cooling
As mentioned in the Introduction, supernovae have been used for decades to constrain models of new physics by way of a cooling argument. Our observations of the neutrino emission from SN1987a suggest a cooling timescale for the protoneutron star of ∼ 10 seconds. For new degrees of freedom to be compatible with this cooling timescale, they must transport energy out of the star at a rate less than the neutrinos. This simply means that new degrees of freedom must transport energy out of the neutrino sphere at a rate less than 3 × 10 52 erg/s [7]. The cooling bound is usually computed more carefully in the free-streaming regime, where analytic computations can produce robust estimates of the escaping flux. However for the trapped regime it usually relies on many approximations and many important aspects have not been taken into account in previous analysis. In this paper, we both extend this bound to the trapped regime using the results of our MC simulation and recompute the bound in the free-streaming regime with gravitational redshift folded in, an effect that was not included by previous papers. The upper bound and lower bound are placed in two different manners due to the fact that the upper bound (stronger couplings to the SM) will be in the trapped regime, while the lower bound (weaker couplings to the SM) will be in the free-streaming regime.
The upper bound is computed straightforwardly using the results of the simulation. The DM profiles produced by the simulation are taken to be steady-state solutions, hence the total flux going through any given radius must be constant throughout the profile. Though the cooling constraint refers to energy transport through the neutrino sphere (∼ 20 km), the flux of dark fermions through this radius will be equal to the flux at infinity. In all regions of parameter space that can be constrained by cooling, the energy sphere for the DM lies well within the neutrino sphere (r E < r ν ), hence we can compute the energy transfer simply by computing the fraction of the non-redshifted spectrum above p min and multiplying by the flux at infinity. The cooling constraint can therefore be expressed as ∞ p min ∂Φ χ ∂E E= √ p 2 +m 2 χ p dp < 3 × 10 52 erg/s (5.1) with ∂Φχ ∂E defined by Eq. 4.3 and p min defined by Eq. 4.6. For the lower bound we can assume that all DM particles produced in the core will free-stream and if their velocity is above the escape velocity, they will carry energy out of the neutrino sphere. The luminosity can be calculated by the volume integral where dL brem /dV and dL e + e − /dV are respectively the local luminosities due to np → npχχ and e + e − →χχ in an infinitesimal volume dV around a point r, and R ν is the radius of the neutrino sphere. This functions are described in Appendix D, and only include particles produced with velocities above the escape velocity at a point r.

Detection
As described in the Summary, liquid xenon (LXe) WIMP detectors are well-suited to observing the high-energy dark fermion flux emitted by supernovae. It may seem at first surprising that a detector designed to detect weak-scale WIMPs would be sensitive to MeV-scale particles. Recall, however, that LXe detectors hunt for WIMPs as a constituent of the ambient galactic dark matter density. As such, the WIMPs are generally fairly cold, traveling with the galactic virial velocity of 10 −3 . In contrast, the dark fermions produced by SN are boosted to semi-relativistic velocities, hence have v ∼ 1. The maximum recoil energy that an impinging DM particle with momentum p could possibly deliver to a xenon nucleus is given by ∼ 2p 2 m Xe . With a WIMP of O(10) GeV (approximately the lower design limit for most LXe experiments) and v = 10 −3 , this is a recoil energy of O(1) keV. Similarly, with a dark fermion of mass O(10) MeV and v ∼ 1, we find a maximum recoil energy of O(1) keV. Unsurprisingly, given the values we chose, LXe detectors typically have thresholds on this order [22]. Since LXe detectors are already searching for WIMPs at the zero-background limit, they make for ideal targets for hunting for sub-GeV DM produced in SN.

Diffuse galactic flux
It is an interesting physical consequence of the semi-relativistic velocities with which the dark fermions are emitted that they will form a diffuse galactic flux of energetic DM. This flux is similar to the diffuse supernova neutrino background (DSNB) (see, e.g., [29] for a review), but with the significant difference that it is due to the overlapping emissions of galactic supernovae, while the DSNB is due to extragalactic supernovae. The reason for this is that, in contrast to the neutrinos, the dark fermions are emitted traveling with an O(1) spread in velocity. This distribution of velocities at emission means that the DM arrives at Earth over a long period of time (comparable to the light travel time to the SN). For galactic SN, this timescale is of order 10 5 years. With an estimated galactic SN rate of roughly 1 per century [30], we see immediately that the dark fermion emissions from up to 10 4 galactic SN can overlap simultaneously at Earth, resulting in a diffuse galactic flux of SN-produced dark fermions. 4 (Note that since SN neutrinos are produced at c, they arrive in a ten-second window. It is clear that the galactic SN rate is insufficient for neutrino emissions from different SN to ever overlap, however the extragalactic rate is suitably large enough for overlap, leading to the existence of the DSNB.) To compute this diffuse DM flux from galactic SN, we take the double-exponential profile of Adams et al. [30] for the core-collapse SN density rate in our galaxy: Since the flux from a given SN falls off with 1/r 2 with r the distance from the SN, we can integrate over this distribution, weighting by the 1/( r − R E ) 2 . The integral therefore takes the form where N χ denotes the total number of DM particles produced in a single SN. Computing this at the Earth's location in the galaxy gives a flux on Earth of Φ diffuse = (2.69 × 10 −54 cm −2 s −1 )N χ . In Fig. 4, we use the N χ produced by our MC simulation to display the magnitude of this diffuse galactic flux on Earth as a function of y, a convenient variable that encapsulates the strength of the DM-SM coupling. It is defined as with the small parameter controlling the kinetic mixing of the SM photon with the dark photon, α D the fine-structure constant of the dark U(1) sector, and m A the mass of the dark photon [19]. The free-streaming and trapped regimes are both apparent in the figure.
At low couplings, the DM free-streams from the PNS and the production scales linearly with y, hence the flux on Earth scales linearly with y as well. For larger couplings, we enter the trapped regime, where the DM is emitted from some approximately blackbody surface. As the coupling increases, this surface moves out to larger radii where the PNS is cooler, hence the DM flux decreases. This diffuse source can be compared to the flux from a hypothetical nearby point source. We find that in order for a single SN to produce a comparable flux of DM on Earth, it would have to sit within roughly 1 parsec of Earth and would have had to have occurred recently enough that the DM flux would still be passing through us. There are no observed SN that unambiguously satisfy these criteria, hence our sensitivity limits are placed using exclusively the galactic diffuse flux. However, if future observations detect such a SN, this would potentially enhance experimental sensitivity to DM flux from supernovae. Point sources and their associated recoil spectra are further discussed in Appendix C.

Count rates in liquid xenon detector
The final necessary piece of this analysis is to determine the detection rate of the diffuse flux in liquid xenon detectors. This is given by the following expression: with dσ dErec the differential DM-Xe cross-section defined in Appendix B.7, E rec the recoil energy of the xenon nucleus, [E thresh , E max ] the recoil energies measured by the detector, (6.5) the differential diffuse galactic flux of dark fermions. Note that the outer integral is taken over p ∞ , the dark fermion momentum at infinity (given by Eq. 4.4), since the scattering is occurring on Earth, however the factors corresponding to the energy spectrum of the DM are in terms of p 0 , the momentum at production, since the distribution is defined at T r E . It is trivial to find p 0 by inverting Eq. 4.4. The limits of integration derive from requiring that the recoil energy be above threshold and less than the maximum recoil energy probed by the detector. Note that since the DM is usually very near the lower threshold for energy deposition and typical values of E max are usually several tens of keV [22], E max plays little role in determining the event rate.
In Figure 5, we show three recoil spectra for a liquid xenon detector. We have set log y = −15.3 and plot a variety of masses. All of these points lie within the interesting region of parameter space for direct detection. It is clear from the figure that lower masses result in lower average recoil energies while the tail of recoil energies can be fairly large for heavier DM owing to its larger kinetic energy. Integrating these distributions allows us to

Results
Our results are summarized in Figures 6 and 7. We have chosen to display the sensitivity limits of the following detectors: 1. Xenon1T: Xenon1T has already completed a one ton-year exposure with no observation of a signal above background [22]. As such, we choose to display the sensitivity region for this exposure. The Xenon1T sensitivity region is shown in red.
2. LUX-Zeplin: LUX-Zeplin is a LXe WIMP experiment currently under construction. When completed, it is projected to be the most sensitive LXe detector to date. It is expected to run for a total integrated exposure of 15 ton-years [25], which is the value we have used in computing our limits. Its reach is shown in yellow.
3. DARWIN: DARWIN is a future LXe experiment designed to be the ultimate LXe WIMP detector, with sensitivity down to the neutrino floor [26]. If constructed, it will have an integrated exposure of 200 ton-years. Its reach is shown in red.
Existing LXe detectors generally have nuclear recoil thresholds of 5 keV [22] but future improvements aim to lower this to 2.5 keV, where solar neutrinos begin to become a large The detector threshold has been taken to be 2.5 keV and the emission timescale from the SN to be log 10 seconds. We compute these curves using the diffuse galactic flux. The region bounded by our cooling bound is overlaid in blue. background [26]. As a result, we have chosen to display the sensitivity limits for both values. Our emission timescale has been chosen conservatively to be log(10) seconds as we do not have a precise notion of the time dependence of the profile at the radii of interest and thus assumed that the χ luminosity will decrease approximately as 1/t in the first 10 seconds, in analogy with the neutrino case. The vertical axis is defined in terms of the convenient variable y, which is an oft-used variable in discussions of these models that serves as a measure of the coupling of the DM to the SM. Recall that y is defined as y = 2 α D mχ m A 4 with the small parameter controlling the kinetic mixing of the SM photon with the dark photon, α D the fine-structure constant of the dark U(1) sector, and m A the mass of the dark photon [19]. There is clearly a degeneracy between the parameters of the dark sector for a given value of y. It should be noted that all of the detection curves presented here are sensitivity regions, not exclusion limits. In other words, at any given point within the reach, the detector is sensitive to some choice of parameters that yields a given y, but is not necessarily sensitive to all choices of parameters. This is an important distinction given that for certain values of α D , the scattering of the dark fermions within the protoneutron star will be dominated by self-scattering, rather than scattering off of protons, an effect neglected in this analysis. We will treat these self-interactions in upcoming work, as well as considering models with extra structure, including a lighter dark photon and cannibalistic interactions [31,32].

Re lic De ns ity
The cooling region is shown in blue. The upper region is calculated in the trapped regime, and is valid under our assumption that the self-interactions can be neglected. The bottom of the exclusion region is obtained from the free streaming regime and should be valid even when considering large self-interactions. Our bounds are stronger than those obtained in Ref. [15] for two main reasons: (1) their analysis only included production through nucleon-nucleon bremsstrahlung, which is subdominant in all of the parameter space we considered to the production from e + e − and (2) their treatment of the trapped regime is more conservative in that they only consider the equivalent of the free-streaming sphere and approximate the dark matter flux as a black-body at that radius. 5 The relic density line is reproduced from Ref. [19] and corresponds to where the relic abundance of dark fermions produced by freeze-out matches the observed dark matter density. It is included for reference. The parameter space constrained by our analysis lies beneath this, meaning that for a standard cosmological history, the dark fermions would not have sufficient cross-section to be depleted down to the measured dark matter density and thus would be overabundant. However these constraints can be avoided by considering non-standard cosmologies with e.g. late entropy injections or by including extra interactions in the dark sector.

Conclusion
The extreme temperatures and densities that are reached during supernovae would create vast abundances of any sub-GeV degrees of freedom in a dark sector. In regions of parameter space where the coupling of the dark sector to the Standard Model is too large to allow the produced dark matter to free-stream out of the cooling protoneutron star, the DM becomes diffusively trapped. In this paper, we focus on a model of an additional U(1) dark sector populated by O(1 − 100) MeV fermions and a heavy dark photon that mixes kinetically with the Standard Model photon. As the dark fermions diffuse out of the star, the flux and spectrum are set by the freeze-out of various interactions. Here, we have used this to calculate the DM flux by employing a dedicated Monte Carlo simulation of particle transport within the protoneutron star. The results allow us to extend the well-known cooling bound into the diffusive regime.
In addition, the fluxes can also be sufficiently large to be detectable in existing liquid xenon WIMP detectors. Due the semirelativistic velocities with which the fermions escape from the star, the arrival time on Earth of the flux from a single SN overlaps with 10 4 other SN, leading to a diffuse galactic flux of dark fermions permeating the Earth. We show that existing and proposed liquid xenon detectors are sensitive to this flux over a large region of parameter space. Future LXe experiments may provide the first direct detection of dark matter at the MeV-scale.
Although we have focused on a particular model of such light dark matter, the same idea applies broadly to many models of DM with mass below ∼ GeV. Existing direct detection results along with SN cooling in the trapped regime may already set important limits on these other models. Perhaps most excitingly, future direct detection experiments could very well discover a wide variety of light dark matter through supernova production.

A Analytic profile of SN
To provide an analytic fit to the results of the full multi-physics supernova simulation described in Section 4.1, we defined a fiducial profile in the following way: with the following fiducial parameters: See Figure 3 for a comparison of this profile to the output of the simulation.

B Cross-sections
In this Appendix we list the cross-sections and rates relevant for the DM dynamics in the supernova.

B.1 χe → χe
This cross-section is relevant for the energy decoupling of dark matter. Since for the cases of interest this is dominated at radii 15 km, we ignore the effect of Pauli-blocking (which would decrease the cross-section and thus lead to a colder spectrum, which is conservative for our estimates). With this approximation the cross-section in the center of momentum (COM) frame is given by where p is the COM momentum and we neglected the electron mass.

B.2 χp → χp
The cross-section in the COM frame is given by where we neglected terms that were suppressed by the proton mass.
For the DM annihilation into electron-positron pairs we take Fermi blocking of the electrons into account since this is a large effect in the the core, where the electron chemical potential is large. Because the cross-section now depends on the electron distribution function we work in the frame of the proton-neutron star and the cross-section will be in terms of the two incoming dark matter momenta p and k. First let us define the following auxiliary functions which appear frequently in the cross-section due to the Pauli-blocking term In order to simplify the expression we will also use the following definitions where cos θ is the cosine of the angle between p and k.
With those definitions, the cross-section is (B.5) Where E and Q where defined in Eq. B.4, and all B i are to be interpreted as B i (E, Q, T, µ) as defined in Eq. B.3.

B.4 Inverse bremsstrahlung annihilation term
Here we compute the DM absorption rate through inverse bremsstrahlung:χχnp → np. We will use the Soft Radiation Approximation (SRA), which is also used in the neutrino production (and absorption) through (inverse) bremsstrahlung [33] and also for computing dark photon production in the proto-neutron star [34]. This approximation allows one to factorize the nucleon-nucleon interaction from the emission process, and the latter can be directly measured by experiment. This approximation is well-justified when the energy of the emitted dark matter pair is much smaller than the COM kinetic energy of the nucleons, ω χ E CM . For us, this is not satisfied for most of the DM masses, and we are usually in a regime where ω χ ∼ E CM . In Ref. [34] it was argued that even in such regime the SRA approximation only resulted in a factor of 2 error in the case of dark photon production. We expect that this approximation leads to an O(1) error in the rate, but as we will find, this rate is subdominant to the annihilation to e + e − almost everywhere in the proton-neutron star by a significant margin.
The absorption rate for DM via inverse bremsstrahlung is given by is the distribution function for DM (including the number of spin dof), n χ is the number density of DM, f p/n is the distribution function for protons/neutrons and|M| 2 χχnp is the averaged matrix element squared for theχχpn → pn process. Now, using SRA, we can rewrite this as with p 1(3) the momentum of the incoming (outgoing) proton and the sum over spin in the previous equation being over the DM spin. Note that in this approximation we drop the momentum of DM in the energy momentum conservation delta function, since in the SRA these are soft compared to the nucleon energy and momentum. Because of this, we can first perform the d 3 k i integrals. For this is it useful to first compute (B.9) of the DM particle. This gives the expression where p is the momentum of the incoming DM particle and E rec is the recoil energy of the xenon nucleus. The nuclear charge is Z = 54 for xenon.

C Recoil spectra from nearby SN
Though there are no observed supernovae that would produce a singular flux in excess of the diffuse SN background of dark fermions discussed in the body of this paper, it is still interesting to consider the recoil spectrum from a single point source. Since the fermions are produced with an O(1) spread in velocities, the arrival time varies between different parts of the spectrum. Dark fermions living on the high-energy (high-velocity) end of the spectrum will arrive far sooner than those on the low-energy (hence low-velocity) end. The majority of the flux will arrive with a delay behind the neutrinos of order the light-travel time to the SN. As a result of this, the recoil spectrum of xenon in a liquid xenon detector on Earth changes over time. Shortly after the arrival of light from the SN, we expect to see a recoil spectrum that extends to high recoil energies (due to the highly-boosted fermions) but with low event rates (due to the fact that the high-velocity fermions live on a tail of the spectrum). As time passes, event rates will increase but the average recoil energy will decrease as the more abundant, less energetic part of the dark fermion distribution begins to arrive on Earth. This evolution is displayed in Figure 8. For the purposes of computation, we have focused on the case of a 30 MeV dark fermion with log y = −16.3 and an SN occurring 30 kpc from Earth (the distance to the galactic center). The recoil spectra are plotted for three different time delays: 10 3 , 10 4 , and 10 5 years after the arrival of the neutrinos on Earth. As expected, the shortest time delay corresponds to the highest energies of dark fermions, hence we have a relatively low yield, but energetic recoil spectrum. As we move towards longer delays, the average recoil energy decreases, but the event rate increases. At 10 5 years (the light-travel time for 30 kpc), we reach the maximal event rate since this corresponds to the arrival of the peak of the dark fermion spectrum. By 10 6 years (not shown), the dark fermion flux is once again very low since it corresponds to the arrival of the low-energy tail. The average recoil energy is well below detector threshold. We find this change in recoil spectrum a noteworthy feature of the SN flux as it could provide a discriminatory tool for detecting a DM flux from a future nearby SN and have included it for completeness.

D Cooling in the free streaming regime
The lower limits of the cooling bound in Figs. 6 and 7 are obtained by considering the free streaming regime of DM produced in the SN. In this case all DM produced in the core can free stream out of the SN as long as it has enough kinetic energy to escape the gravitational attraction due to the proto-neutron star.
In order to compute the minimum escape energy from a region of radius r we need to compute the metric inside the proto-neutron star. Following Ref. [35], the metric can be written as ds 2 = B(r)dt 2 − A(r)dr 2 − r 2 dΩ 2 . with p the pressure in the star. Given that the pressure term is subdominant, we can approximate the pressure by treating the protons and neutrons as a gas of degenerate fermions to the level of precision we are interested in. The minimum energy required for DM to escape from a radius r is given by As discussed in previous sections there are two important production channels which contribute to the DM production: electron-positron annihilation to DM and DM bremsstrahlung from proton-neutron scattering. For the profile used in our work we found that the latter yields a larger production for all masses of interest, but we include both contributions for completeness.
For the bremsstrahlung case, we use a similar calculation to what was done in Sec. B.4. However, because we are now interested in the energy flux and not the number flux, and because we must impose a minimum energy due to gravitational trapping, we cannot utilize that result which was obtained via detailed balance. The steps to compute the production cross-section are almost identical, except that one must impose a maximum energy cutoff for the DM produced by hand, since due to the SRA the energy of the DM no longer appears in the energy conserving delta function. For that purpose we include an exponential regulator exp[−(ω 1 + ω 2 )/T ], where ω i is the DM energy and T the temperature. 7 Using this, the local DM luminosity from this channel is given by dL brem dV = 64αy 9π n p n n (πm p T ) 3/2 4 4 + x 1 x 2 (3x 2 1 + 4x 1 x 2 + 3x 2 2 ) +(5x 2 1 + 12x 1 x 2 + 5x 2 2 ) , (D.6) where the first integral over dx i correspond when both pair-produced DM have energy above the escape energy m χ / √ B and the second one when only one of them does. For the electron-positron annihilation term the full form of the production above a certain energy threshold is very complicated due to the average over the initial electron 7 Another option is to introduce a hard cutoff on the DM energy such that ω1 + ω2 ≤ | p1 − p2| 2 /(4mp), where p 1(2) are the nuclei initial momentum. This guarantees that the produced energy is smaller than the COM kinetic energy of the nuclei. We found that the exponential regulator gives a smaller (and thus more conservative) rate, and also that it gives an answer that is closer to satisfying detailed balance when compared to the absorption rate. and positron momentum. In order to simplify our treatment we compute the luminosity such that half the COM energy √ s/2 is above m χ / √ B and consider that √ s/2 of energy is carried away (i.e. we only consider the energy carried by the particle which gains energy from the boost from the COM frame to the star frame). Since we do not include the enhancement of the energy due to the boost, and only consider one of the produced DM particles for the luminosity, this leads to a conservative estimate. Using this the luminosity from electron positron annihilation is given by where the Θ ensures that the COM energy is above the escape energy.