Supernova Constraint on Self-Interacting Dark Sector Particles

Allan Sung, 2, ∗ Gang Guo, 3, † and Meng-Ru Wu 4, 5, ‡ Institute of Physics, Academia Sinica, Taipei, 11529, Taiwan Department of Physics, Princeton University, Princeton, NJ 08544, USA Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210033, China Institute of Astronomy and Astrophysics, Academia Sinica, Taipei, 10617, Taiwan National Center for Theoretical Sciences, Physics Division, Taipei, 10617, Taiwan (Dated: February 10, 2021)


I. INTRODUCTION
The detection of 20 electron antineutrinos emitted from the core-collapse supernova (SN) explosion event, SN1987a [1][2][3], not only broadly confirmed the prevalent SN theory, but also led to several important consequences to fundamental physics, including e.g., bounds on the neutrino decay lifetime, the absolute masses of neutrinos, and non-standard neutrino interactions [4][5][6][7][8][9][10][11][12][13]. In particular, important constraints on a variety of particles beyond the Standard Model (SM) including the axion, sterile neutrino, dark photon, etc., whose masses are sub-GeV, were derived , which complement ongoing experimental searches for those particles. These constraints were based on the requirement that the exotic particles should not carry away an amount of energy from the cooling proto-neutron star (PNS) more than the inferred total energy carried by neutrinos, E ν 3×10 53 erg (see, however, a caution from Ref. [49]). In addition to the SN cooling (more precisely, the PNS cooling) constraint, recent work also proposed new constraints on light dark photon or dark photon portal light dark matter based on other SN-related observables, such as the measured SN explosion energy [50], the γ-rays [51,52], or the produced (semi-)relativistic dark matter flux arriving at the terrestrial detectors [53]. * as111@princeton.edu † gangg23@gmail.com ‡ mwu@gate.sinica.edu.tw One important aspect in deriving the SN constraint on light dark sector particles is that their interaction with SM particles cannot be too strong for them being trapped inside the PNS. We note that previous studies ignored completely the self-interactions between dark sector particles when deriving the SN bounds. However, if the abundance of dark sector particles inside the SN core can be as large as SM particles, and if the self-interaction cross-section can be as large as the neutrino-nucleon scattering cross-section ∼ O(10 −41 ) cm 2 for neutrinos of ∼ O(10) MeV, dark sector particles can trap themselves inside the PNS. Consequently, SN bounds on selfinteracting dark sector particles can be largely evaded. Note that self-interacting dark matter has been considered to be a viable option to resolve a number of tensions in the scale of galaxies or galaxy clusters, e.g., the corecusp, too-big-to-fail, and the missing satellites problems; see e.g., Refs. [54][55][56][57]. Extensive efforts investigating consequences of self-interacting or annihilating dark matter on various cosmological and astrophysical signatures have been pursued in recent years, e.g., Refs. [58][59][60][61][62][63][64][65][66][67][68][69][70][71].
In this work, we aim to address the issue of SN constraints on self-interacting dark sector particles by taking into account their self-trapping effect for the first time. We use a widely-examined dark photon portal dark sector model explicitly to compute all the relevant interaction cross-sections and the decay rates. In principle, to determine precisely the dark sector particle fluxes emerging from the PNS requires solving full Boltzmann transport equations in a way similar to the neutrino transport problem in SNe (see e.g., a recent review [72]   I: Relevant processes of dark photon interactions considered in this work. "Abs." refers to a process which absorbs dark photon(s) or decay of dark photon. "Sca." refers to a scattering process of a dark photon with a Standard Model (SM) particle or a dark sector (DS) particle. Only leading-order processes are included here. Note that for m A > 2m χ , A → χχ also accounts for contribution from A χ → A χ (see Eq. (A13) and text below for details).
therein) 1 . Such approach demands intensive computational power to fully incorporate the scattering kernels and particle annihilation. Instead of directly pursuing full numerical simulations, we adopt an approximated approach to estimate the energy fluxes carried by dark photons and dark fermions evaluated in the non-diffuse regime and diffuse regime separately, and formulate a physically motivated criterion to switch from one regime to another. This approach allows us to estimate the effect of dark sector self-trapping on SN bounds for a wide range of parameter space, which turns out to be very important even for small couplings in the dark sector. The rest of the paper is organized as follows. In Sec. II, we describe the underlying dark photon portal dark sector model, the considered SN model, and list the relevant interactions and decay processes that we included in this work. In Sec. III, we compute the energy luminosity of dark sector particles leaving the PNS in the non-diffuse and the diffuse regime, respectively, and formulate the criterion to switch from one regime to another. We apply this method to derive SN bounds on self-interacting dark sector particles in Sec. IV. Our conclusion and discussions of potential caveats as well as other implications are given in Sec. V. All detailed derivations of the crosssections, the decay rates, and the diffusion luminosity of dark sector particles are given in the Appendices. We adopt natural units withh = c = 1 throughout the paper unless explicitly specified.   Table I. Only leading-order processes are included here. Note that for m A > 2m χ , χχ → A also accounts for DS processes of A χ → A χ and χχ → χχ, as well as SM processes involving a pair of dark fermions (see Eq. (A13) and text below for details).

A. Dark Sector Model
We consider a dark photon portal dark sector (DS) model wherein the Dirac dark fermion χ couples to dark photon A and the dark photon kinetically mixes with the SM photon [73][74][75]. The corresponding Lagrangian of the dark sector is given by where is the mixing parameter, m A is the dark photon mass, m χ is the mass of the dark fermion, and g D is the DS coupling constant. We define the DS fine structure constant α D ≡ g 2 D /4π analogous to the electromagnetic fine structure constant α e ≡ e 2 /4π.
Through the mixing of dark photon with the SM photon, dark photons and dark fermions can be produced via processes analogous to the SM electromagnetic ones. We follow Refs. [33,35] to consider the following production channels for light dark photons and dark fermions inside the hot and dense interior of a PNS with a temperature 30 MeV and a core mass density > ∼ 10 14 g cm −3 . For the dark photon, we include the nucleon-nucleon bremsstrahlung np → npA , Comptonlike interaction γe − → e − A , and electron-positron annihilation e − e + → A . For the dark fermion, we consider three pair-production channels including nucleon-nucleon bremsstrahlung np → npχχ, electron-positron annihilation e − e + → χχ, and plasmon decay γ * → χχ.
When the PNS interior is optically thin to dark photons or dark fermions, the rates of the above production channels directly determine the energy luminosity carried away by the dark particles. However, when the interactions between dark particles and the SM medium, as well as the self-interactions in the DS, are strong enough, dark photons and fermions can be trapped in the PNS. These interactions include the inverse processes of the above production channels, the dark photon (fermion) pair-annihilation A A ↔ χχ, the DS Compton scattering A χ ↔ χA , the scattering of dark fermions χχ → χχ and χχ → χχ, as well as the dark fermion scattering with SM charged fermions χp → χp and χe − → χe − . These interactions determine the mean free path of the DS particles and thus the energy loss rate through the neutrinosphere in the diffuse regime. Moreover, when m A > 2m χ , the (inverse) decay of dark photon A ↔ χχ needs to be considered. In this scenario, the inverse decay process dominates all the dark fermion pair-absorption processes and the DS self-interactions due to the resonances 2 . Other DS self-interactions without resonances are suppressed by an extra factor of α D compared to the (inverse) decay rate and thus can be neglected for m A > 2m χ . We list all the included interactions in this work in Table I and Table II for dark photon and dark  fermion, respectively. Throughout this work, we will use perturbative calculation for the interaction rates. Since the temperature in the PNS is ∼ 30 MeV, we only consider DS particles of masses below 1 GeV.

B. PNS cooling constraint and SN model
Without the presence of dark particles, long-term SN simulations predicted that the PNS cools by emitting neutrinos of all flavors in ∼ 10 s. The total energy carried away by neutrinos is 3 × 10 53 erg, which is fixed by the available gravitational energy released to form a compact neutron star. The potential emission of exotic particles that may be produced inside the PNS will thus reduce the duration of the neutrino emission and can be constrained by the observed neutrino events from SN1987a. Based on comparisons with simulations including the emission of axions from the PNS, a well-known criterion (Raffelt's criterion) was formulated to constrain the maximal energy luminosity carried by any exotic particles [77]: where L D denotes the energy luminosity of dark emission, evaluated during the early PNS cooling phase, and L ν is approximately the time-averaged neutrino energy luminosity in SM.
In this work, we use the PNS density and temperature profile at 1 s post the SN core-bounce obtained in Ref. [76] to compute the emission of dark photons and dark fermions. This PNS profile was widely used recently for similar purposes (see e.g., Refs. [33,50]). Note that the choice of a particular SN model may introduce uncertainties of a factor of a few for the derived bounds [32,33]. Fig. 1 shows the radial evolution of the density ρ, temperature T , the electron chemical potential µ e , and the plasma frequency ω p (see Eq. B7). The position of the spectral-averaged neutrino decoupling sphere, i.e., neutrinosphere, is indicated by the vertical dash line at R ν 22 km. The density profile shows a monotonically decreasing behavior as a function of radius, while the temperature profile exhibits a peak at r 11 km, due to the inefficient compression heating at the densest core region. The electrons are highly degenerate inside the PNS. The plasma frequency ω p 14 MeV at the PNS center and decreases at larger radii. The plasma effect effectively alters the mixing between the dark photon and the SM photon differently for the transverse and the longitudinal polarizations [33,34,78]. We have included this effect throughout this work and give the details in the Appendix B.

III. LUMINOSITY OF DARK SECTOR PARTICLES
In this section, we first describe how we compute the energy luminosity of DS particles leaving the PNS for the scenario where the DS self-trapping can be ignored (Sec. III A), and for cases where they can be considered as diffusive due to self-trapping (Sec. III B). We formulate the criteria that determine if a DS particle species is in diffuse regime or not in Sec. III C. In Sec. III D, we then apply our formalism to the adopted PNS profile to compute the total luminosity carried away by DS particles from the PNS interior.
In the rest of the paper, we denote the 4-momentum of A by k = (ω, k), that of χ by p = (E, p) and that of χ by p = (E , p ), unless noted otherwise.

A. Non-diffuse regime
In the non-diffuse regime, we consider the bulk emission rates of DS particles inside the neutrinosphere and the attenuation due to absorption and decay, following Refs. [33,50]. The luminosity of the dark photon is given by where g L = 1, g T = 2, and Γ L,T A ,prod and τ L,T are the production rate and optical depth respectively. The expo- nential factor accounts for the absorption of dark photons by the medium and their decay. We separate the longitudinal (L) and transverse (T) modes because the medium effect modifies their dispersion relations and leads to different effective mixings of the two modes with the SM photon (see Appendix B). The production rate Γ L,T A ,prod is determined by the interactions involving the SM particles listed in Table I. By detailed balance, the rates of each production process and its inverse process are related by Γ L,T A ,prod = e −ω/T Γ L,T A ,abs , where Γ L,T A ,abs is the total absorption rate of the inverse process. Therefore, For the optical depth, we include the absorption and decay of A by the same processes with SM particle as above and ignore those involving DS particles, to avoid double counting the total DS luminosity (see also later discussions in this subsection) 3 . This gives where v = | k|/ω is the dark photon velocity, and f (r) is a geometric factor used in [33] that effectively takes into account different path-lengths of dark photons emitting locations to the neutrinosphere. The explicit forms of these absorption rates are given in Appendix C. Next, we compute the luminosity of the dark fermion (χ) in the non-diffuse regime as where g χ = 2 is the physical degrees of freedom of χ, and Γ χ,prod is the production rate of χ by the SM medium (see Table II). Note that since χ andχ are symmetric in our model, the total luminosity in the dark fermion pair is L χ + Lχ = 2L χ . Here we do not include the attenuation due to pair absorption. This is because in the non-diffuse regime, the dark fermion abundance below the neutrinosphere is very low, which suppresses the pairabsorption of dark fermions. If the dark fermions were in equilibrium with the medium, detailed balance could relate the production rates to the absorption rates by Γ eq χ,prod = e −E/T Γ eq χ,abs . We take the equilibrium production rates Γ eq χ,prod as an approximation for the production rates Γ χ,prod used in Eq. (6). That is, Γ χ,prod e −E/T (Γ χχnp→np + Γ χχ→e − e + + Γ χχ→γ * ). (7) This approximation in principle underestimates a bit the production rates of χ andχ due to the assumed equilibrium occupation number, which effectively results in Pauli-blocking. However, since the χ andχ are always produced in pairs, the effective Pauli-blocking suppression of the rates are small due to their zero chemical potentials.
We note that when m A ≥ 2m χ , the dark fermion pair production is in fact dominated by the decay of on-shell dark photons. Since in Eq. (3) we do not include the decay of A to χ andχ when A are in non-diffuse regime, it leads to double-counting of the total dark luminosity if we also include the dark fermion production. Thus, we do not consider the contribution of Eq. (6) when dark photons are in non-diffuse regime and when m A ≥ 2m χ . We also note here that if dark photons are in the diffuse limit but dark fermions are not, then effectively the nondiffuse dark fermion luminosity can be determined by the decay rate of the trapped dark photons that are in thermal equilibrium with the SM medium (see Appendix A) 4 .
Before we discuss the detailed numerical results of the DS luminosity in the non-diffuse regime, let us provide an analytic estimation for the relevant region of the parameter space. First, for most cases, the production rates of longitudinal dark photons and dark fermions are suppressed by a factor of m 2 A /ω 2 (see Appendix C) and by the coupling constant α D , respectively, compared to the production rate of the transverse dark photons. We can thus approximate the non-diffuse luminosity of DS particles by considering transverse dark photons only. Second, we consider for simplicity a homogeneous PNS with radius R c 10 km, temperature T 30 MeV, density ρ 3 × 10 14 g/cm 3 and electron fraction Y e 0.3. We also assume that dark photons are relativistic. Taking the nucleon-nucleon bremsstrahlung A np → np without 4 Note that we consider separately the diffuse condition for the longitudinal and the transverse dark photons. Thus, we further multiply the non-diffuse dark fermion luminosity by a factor of 1/3 (2/3) if only the longitudinal (transverse) mode of the dark photons are trapped (see Appendix D).
the plasma effects (See Appendix C 1), the luminosity of DS particles, L D , is approximately for m A < ∼ 1 GeV.

B. Diffuse regime
In the diffuse regime, we assume that DS particles are in good thermal contact with the SM medium. Due to the temperature gradient, the DS phase space distributions are slightly anisotropic, which induces an outward energy flux. We use the radiative transfer equation for the DS particle energy flux through the neutrinosphere in Appendix E. The energy flux of a particle species i is approximately given by where the upper (lower) sign is for fermions (bosons), g i is the physical degrees of freedom of particle i, T ν is the temperature at the neutrinosphere, m i is the mass of particle i, and λ −1 i (r) is the thermally averaged inverse mean free path (IMFP) 5 of particle i at radius r defined by where f i (E, T (r)) is the distribution function of particle i at radius r. We distinguish between the absorptive and scattering IMFP, λ −1 i,abs and λ −1 i,sca , and define the total IMFP as As in the non-diffuse regime, we compute the energy fluxes of the longitudinal and transverse dark photon separately, and assume that the energy fluxes of χ andχ are equal. Hence, the total energy loss rate in the DS particles is L tot = L A ,L +L A ,T +2L χ . The IMFP calculations can be found in Appendix C. We now estimate the relevant parameter space region in the limit where the DS self-interaction is the dominant opacity source. When m A > 2m χ , the dominant DS selfinteraction is the (inverse) decay A ↔ χχ. Given T ν 3.9 MeV and the temperature gradient |dT /dr| Rν 6.2× 10 −4 MeV/m, if we fix m A /m χ = 3, the luminosities of the DS particles can be fitted by for m A < ∼ 25 MeV, and for m χ < ∼ 30 MeV. If we take m A = 9 MeV and m χ = 3 MeV, then the total DS luminosity is The corresponding χχ → A cross-section is 6 When m A < 2m χ , the diffuse luminosities of dark photons and dark fermions depend on the values of and α D . For regions where is small enough such that the dominating opacity source for dark photons is DS selfinteraction A A → χχ, one may estimate the DS luminosity by considering the dark photon luminosity only. If we fix m A /m χ = 1/3, then the dark photon luminosity can be fitted by for m A < ∼ 25 MeV. Adopting m A = 3 MeV and m χ = 9 MeV, this gives (17) 6 The estimated value here is the thermally averaged cross section evaluated at the center of PNS σ χχ→A ≡ λ −1 χχ→A (r = 0) /nχ(r = 0).
The corresponding χχ → χχ cross-section is 7 and the corresponding A A → χχ cross section is These analytical formulas Eqs. (14)- (19) already illustrate that DS particles can be self-trapped by interaction cross-sections of > ∼ O(10 −40 ) cm 2 such that the corresponding diffuse DS luminosity L D < ∼ L ν . Only very weakly interacting DS particle with α D much smaller than those nominal values in Eqs. (14)- (19) can result in L D > L ν .

C. Diffuse Criteria
Whether DS particles are in the non-diffuse or in the diffuse regime sensitively depends on their abundances in the PNS. For example, if the dark fermion abundance is significant, the contribution of their pair-absorption processes (e.g. χχnp → np, χχ → e − e + and χχ → γ * ) to the dark fermion optical depth cannot be neglected. Furthermore, the DS abundances lead to significant DS selfinteractions (e.g. χχ → χχ, χχ → A A , χA → A χ,...) that can trap themselves in the PNS. These DS selfinteractions delay the escape times for both the dark photon and the dark fermion, which enhances their abundances and optical depths. Therefore, the abundances of DS particles is critical in determining whether selftrapping is important. Below, we formulate the diffuse criteria in terms of the DS abundances and their IMFP at the neutrinosphere.
To obtain the exact DS abundances, detailed transport incorporating their production, scattering, and absorption needs be solved. Here we roughly estimate the abundance of each particle species by their production rate and a relevant timescale. The abundance of a DS particle species i is approximated as where the quantity ∆t i is subject to the following considerations. First, the longest possible time for DS particles to accumulate is the cooling time of the PNS. We set this upper bound of ∆t i as t cool = 1 s. Second, the shortest possible timescale, i.e., the lower bound of ∆t i , is estimated by the free-escaping timescale t free ≡ R ν /v i . Another characteristic timescale here in the limit where 7 The cross section is evaluated in the same way as Eq. (15).
is the diffusion coefficient (see below for details). For t free < t i,diff < t cool , we then take t i,diff as ∆t i . Combining these criteria, the relevant timescale ∆t i used in Eq. (20) is given by Since the PNS is not homogeneous, the diffusion of DS particles cannot be described by a constant D i . As shown in Tables I and II, the interactions relevant to the IMFP of DS particles can be categorized as SM-type or DS-type by particles involved in the interactions. We choose particular locations in the PNS where the IMFPs of each type are maximal to roughly estimate their contribution to the diffusion time scale 8 . For the SM-type interactions, we estimate the IMFP at the center (r SM = 0 km) of the PNS where the density of SM particles is the largest. As for the DS-type interactions, we estimate the IMFP at r DS = 11km where the temperature is the highest. This location has the largest DS particle density, if DS particles are in thermal equilibrium with the SM medium. Knowing the IMFPs of both types at these respective locations, we can then compute the diffusion time scale whereλ −1 i,SM andλ −1 i,DS are SM-type and DS-type contribution to the IMFP of particle i. When computing the IMFP contribution from DS particles,λ −1 i,DS in Eq. (22), we take the assumption that DS particles are in full thermal equilibrium with the SM medium. Combining Eqs. (20), (21), and (22) then allows us to compute the DS abundance N i for particle i.
After deriving N i , we then define a scaling factor η i ≡ N i /N eq i for each dark particle species i for N i < N eq is the corresponding total equilibrium abundance of i. These scaling factors are then used to better estimate the IMFP, λ −1 i,DS due to the DS particles self-interaction, used in Eq. (10) for computing the DS diffuse luminosities. We also use η i to calculate the thermally averaged total We are now finally ready to write down our diffuse criteria. We say that a DS particle species i can be treated as in diffuse limit if the following conditions are satisfied.
The first condition ensures that a DS particle species i is only considered to be diffusive when their production is efficient enough such that their amount can exceed the equilibrium number within the relevant timescale ∆t i . The second condition requires that the IMFP of i at R ν , where both the density and temperature are the lowest in the PNS, is large enough to trap particle i inside the PNS (see Fig. 1) 9 . When Eq. (24) is satisfied, we use Eq. (9) to compute the DS diffuse luminosity for i. Otherwise, we take the non-diffuse luminosities, Eqs. (3) and (6) for dark photons and dark fermions, respectively. We provide in Appendix D a detailed workflow to additionally describe how we use the results in these sections to compute the dark sector luminosity for interested readers.

D. Numerical Calculations
We now apply formulas derived in previous sections to compute the dark photon and dark fermion luminosities and examine how they depend on the coupling constants and α D . Fig. 2 shows the luminosities of different DS particles as functions of α D with = 10 −8 for two different choices of DS masses. In the left panel, the masses are chosen such that the decay process A → χχ is allowed. In this scenario, dark photons of different polarization modes are trapped diffusively when α D > ∼ 5 × 10 −18 , while the dark fermions are diffusive when α D > ∼ 10 −17 . The non-diffuse luminosity of the dark photon for small α D mainly depends on the interaction with the SM particles, and thus is independent of α D 10 . The dark fermion luminosity is proportional to α D in the diffuse limit, as analyzed in Appendix A. However, the dark fermion luminosity for α D < ∼ 5 × 10 −18 is set to zero to avoid double counting because they are predominantly produced through the decay of on-shell dark photons (see discussions in Sec. III A). In the diffuse regime, the luminosities of the DS particles are proportional to α −1 D due to the selftrapping interactions effectively dominated by the (inverse) decay process A ↔ χχ. 9 Note that in computing the second criterion, we exclude the contribution from the dark photon decay via A → e − e + This is because including this decay process leads to an artificially enhancement of the IMFP by several orders of magnitudes for dark photon heavier than 30 MeV, as the electron chemical potential µe 15 MeV at Rν . However, this process should be strongly Pauli-blocked for m A < ∼ 200 MeV in most region inside the PNS (see Fig. 1) where dark photons are most abundant in the diffuse limit. 10 The difference of ∼ 10 2.5 between the luminosities of the longitudinal and transverse dark photon is due to the plasma effects. Note that χ is in the non-diffuse limit for the range of α D shown in the right panel. For the right panel in Fig. 2, we choose DS masses such that the decay process A → χχ is not allowed. In this case, dark photons of different polarization modes are in diffuse regime when α D > ∼ 8 × 10 −9 , while dark fermions are always in the non-diffuse regime for the range of α D shown in the plot. The main reason leading to several orders of magnitude larger difference in α D here than the previous case is due to the extra α D dependence in the IMFP of DS self-interactions (see e.g., Table I or Eqs. (14) and (17)). Similar to the previous scenario, the non-diffuse luminosity of the dark photon is independent of α D . However, the luminosities in the diffuse regime scale as α −2 D . Once again, this is because the dominant interaction responsible for trapping the dark photon being A A → χχ, whose interaction rate is proportional to α 2 D . The dark fermion luminosity is significantly smaller than the dark photon luminosity for the range of α D in the plot, because the dark fermion production rate is proportional to 2 α D , suppressed by an extra factor of α D when compared with the dark photon production rate.
In Fig. 3 we show the luminosities of the DS particles as functions of with two choices of α D and same DS masses as those in Fig. 2. For the case with α D = 10 −17 in the left panel where the decay process A → χχ is again allowed, the longitudinal (transverse) dark photons are trapped when > ∼ 3 × 10 −9 ( > ∼ 2 × 10 −9 ). For dark fermions, they just free stream independent of the value of , because of the small value of α D (c.f., Fig. 2). The non-diffuse luminosities of the dark photons of both modes are proportional to 2 as determined by their production rates. The dark fermion luminosity for < ∼ 2 × 10 −9 is set to zero for the same reason discussed above. The diffuse luminosities of dark photons are affected by both the interactions with the SM particles and the DS self-interactions. For > ∼ 10 −7 , dark photon-SM interactions dominate over DS self-interactions, so the dark photon luminosities are proportional to −2 as determined by the mean free path of the dark photon-SM interactions. Interestingly, the dark fermion luminosity becomes independent of and can thus remain larger than L ν for large . This is because the branching ratio of the dark photon absorption processes approaches to 1, such that dark fermions can be produced via the decay of trapped dark photons (see Appendix A).
For the right panel in Fig. 3, we choose α D = 10 −9 for the case where A → χχ is not allowed. There, longitudinal (transverse) dark photons are diffusively trapped , while dark fermions are in the diffuse limit when > ∼ 2 × 10 −5 . The dark photon luminosities are proportional to 2 for < ∼ 10 −9 , where the optical depths are negligible. The diffuse luminosities of dark photons scale as −2 because the absorption processes by the SM medium is dominant. The dark fermion luminosity is much smaller than that of the dark photons until > ∼ 10 −5 , because the dark fermion production rate is suppressed by an additional factor of α D when compared with the dark photon production rate. The diffuse luminosity of the dark fermion is proportional to −2 for > ∼ 10 −4 , where χe − → χe − dominates the IMFP at R ν .

IV. COOLING BOUNDS ON SELF-INTERACTING DARK SECTOR PARTICLES
In this section, we examine the excluded parameter space of the DS. Since the model has four free parameters, m A , m χ , , and α D , we choose to project the exclusion contours on various combinations of these parameters. As discussed in the previous sections, whether the decay of A → χχ is allowed affects the DS luminosity significantly. Here we choose two benchmark mass ratios m A /m χ = 3 and 1/3 to investigate the SN bounds for these two scenarios.
We  2)) is violated. In the left panel where A → χχ is allowed, the excluded region exhibits an L-shape, covering the parameter ranges of (I) α D < ∼ 10 −16 for 10 −10 < ∼ < ∼ 10 −6 , as well as (II) 10 −19 < ∼ α D < ∼ 10 −16 for > ∼ 10 −6 . For the majority of the parameter space in (I) (α D < 10 −16 and 10 −10 < < 10 −6 ), the DS luminosity is mainly contributed by the non-diffuse dark photons (see also Figs. 2 and 3). In (II) (10 −19 < ∼ α D < ∼ 10 −16 for > ∼ 10 −6 ), L D is mainly contributed by the non-diffuse dark fermions (see also Fig. 4). Close to the upper edge of the L D = L ν contour in both regions, self-trapping of DS particle takes effect such that L D ∝ α −1 D decreases with increasing α D (see Figs. 2), leading to the horizontal edge at α D ∼ 10 −16 . The lower bound of α D ∼ 10 −19 for region (II) is due to the inefficient production of dark fermions. For the right panel in Fig. 4 where A → χχ is not allowed, two regions similarly exist. Region (A) with 10 −9 < ∼ < ∼ 10 −6 and α D < 10 −8 receives dominant contribution from nondiffuse dark photons as region (I) in the left panel (see also Figs. 2 and 3). Similarly, the self-trapping of DS particles defines the upper edge of L D = L ν at α D ∼ 10 −8 . For the narrow diagonal-shape region (B) at > ∼ 10 −5 , L D is dominated by dark fermions. This shape is related to the fact that both the production and pairabsorption rates of dark fermions from the SM medium are proportional to 2 α D , as discussed also in Sec. III D.
For both scenarios shown in Fig. 4, there are specific values of α D above which the cooling criterion gives no constraint due to DS self-trapping. We now investigate the dependence of these critical values of α D on DS particle masses. Fig. 5 shows the excluded regions on the α Dm χ plane with = 10 −8 for fixed mass ratios m A = 3m χ (left panel) and 3m A = m χ (right panel). Note that regions below the contours are excluded. In the left panel where the decay A → χχ is allowed, the SN bound only weakly depends on α D for m χ < ∼ 10 MeV. On the other hand, α D increase sharply with m χ for m χ > ∼ 10 MeV. The main reason is the DS particle abundances are insensitive to the mass for m χ 30 MeV which is the typical temperature inside the PNS. Larger m χ (and m A ) leads to smaller DS abundances, which in turn gives rise to a smaller IMFP of DS self-trapping for a given α D . Thus, the suppression of L D due to DS self-trapping occurs at a larger α D for larger m χ (and m A ). For m χ > ∼ 70 MeV (m A > ∼ 210 MeV), no SN bounds can be placed due to the inefficient production of DS particles. In the right panel where 3m A = m χ , the bound is completely determined by the dark photon only beacuse the dark fermion production is further suppressed by an extra factor of α D . Here, the critical α D also increases with m χ sharply for m χ > ∼ 10 MeV as in the left panel. Once again, this is resulting from smaller dark photon abundance for larger m χ (thus m A ) inside the PNS. It thus requires a larger α D for A A → χχ to self-trap the dark photons. Note that the maximal m χ 800 MeV below which SN bound exists corresponds to a maximal m A 250 MeV, consistent with the maximal m A in the left panel.    The above examples clearly demonstrated how efficiently a small α D within the self-interacting DS can affect the SN bound. In Fig. 6, we further show the excluded regions in the -m A plane with two different values of α D for each choice of the DS mass ratio. Also shown are the bounds derived by considering the dark photon degree of freedom only, as well as the existing experimental constraints on dark photon (extracted from Ref. [79]). We refrain from showing other astrophysical or cosmological bounds on dark photon of similar masses; see e.g., Refs. [50,52,78,[80][81][82][83][84][85][86]. These plots show once again that even for small values of α D = 10 −12 (left panel) and 10 −5 (right panel), the SN exclusion regions shrink significantly due to the self-trapping effects when compared to results derived by considering only dark photons without self-trapping. However, for very small values of α D , e.g., 10 −17 and 10 −9 shown in the left and right panels, respectively, the cooling bounds can be extended to larger . This is due to the contribution of dark fermions through the decay of trapped dark photons, as discussed in earlier sections. Note that for α D = 10 −12 in the left panel, the excluded region slightly extends to larger values of m A for > ∼ 10 −7 . This is because the dark fermion luminosity depends on both the decay rate of A → χχ and the blackbody energy density of dark photons (see Appendix A). The larger decay rate can compensate the smaller dark photon energy density with increasing m A , provided that the branching ratio   of the dark photon absorption processes is close to 1.

V. CONCLUSION
In this work we examined the SN bounds on selfinteracting dark sector particles. Adopting a dark photon portal dark sector, we derived the relevant interaction cross-sections and (inverse) decay rates for reactions listed in Tables I and II. We then used these to compute the energy luminosities of dark sector particles in the non-diffuse and diffuse regimes separately, and formulated a simple criterion to connect these two regimes. The self-interaction of dark sector particles can efficiently trap themselves inside the proto-neutron star and thus suppresses their energy luminosities.
Comparing the dark sector luminosity with the neutrino luminosity inferred from the SN1987a event, we derived SN bounds for two assumed dark photon to dark fermion mass ratios, m A /m χ = 3 and 1/3, which represent scenarios where A → χχ is allowed or not. For the former (later) case with m A /m χ = 3 (m A /m χ = 1/3), SN bounds only apply to weakly interacting dark sectors whose dark fine structure constant α D < ∼ 10 −15 ( < ∼ 10 −7 ), for m χ < ∼ O (20) MeV (see Figs. 4 and 5). The dominating dark sector cross-sections for these α D values correspond to 10 −40 cm 2 . In particular, there is no SN bound for the former case for α D > ∼ 10 −7 (Fig. 5). Although the exact excluded regions in the DS parameter space should also depend on the chosen value of m A /m χ , which are unexplored in this work, our results demonstrated that when applying the supernova bounds to dark sector particles, their self-interactions, which can evade the bounds, must be taken into considerations.
Our results here also imply that other stellar bounds, e.g., from the horizontal branch stars, tips of red giants, or white dwarfs, on dark sector particles may also be sensitive to the structure of the dark sector. Similarly, our results also indicate that for non-standard strongly selfinteracting neutrinos proposed to resolve the Hubble tension [87], the needed strong self-interaction, ∼ 10 6 or 10 9 times stronger than the SM weak interaction, will likely lead to self-trapping of neutrinos and results in inconsistency with the SN1987a observation (cf. Eqs. (15) and (19)). Moreover, we would like to point out that although the self-interaction of dark sector can completely evade the SN bound, new constraints may be further derived by considering their potential late-time heating to the remnant NS via decays or annihilations in a longer timescale. Furthermore, such self-trapping effect might provide an efficient mechanism to produce DM-admixed NS, which might have implication for the GW detection of binary NS merger events [88][89][90]. All these aspects are beyond the scope of this paper and can be further investigated in future work. We also note that our results indicate that self-interacting DM that can help solve the small-scale issues in galaxies can not be excluded by the SN cooling bound, as the required χ-χ self-interaction cross-section σ χχ /m χ ∼ O(1) cm 2 g −1 largely exceeds values that can be constrained by SN cooling (cf. Eq. (18)).
Finally, we comment on potential caveats in this work. First, our criterion of switching from the non-diffuse to diffuse regime is rather abrupt and can sometimes create non-negligible discontinuity in dark sector luminosities as shown in e.g., Figs. 2 and 3. In reality, the transition should be smooth and this can possibly introduce errors of a factor of a few in all our derived bounds. For example, while evaluating the luminosities in the non-diffuse regime, we only included absorption and decay processes in the opacity. In principle, scattering can somewhat reduce the energy luminosity of dark sector particles leaving the PNS, before the condition for diffusion is fully satisfied. Also, we used a sharp neutrinosphere as a boundary to estimate the luminosities in the diffuse regime. This may lead to some errors when the dark sector particles are not fully in the diffuse limit. Second, when we evaluated the diffusion timescale used to determine the diffuse criterion, we selected two specific locations where the IMFPs are largest for simplicity. This approximation may overestimate the diffusion time a bit. All these sources of uncertainties can only be addressed by performing a full numerical calculation of multi-dimensional Boltzmann transport and can be pursued in future. However, the main conclusion derived in this work -self-interactions inside the dark sector can crucially affect stellar bounds -should remain relatively solid and needs to be considered in all relevant studies.
where J µ is the final state current that couples to the dark photon, L , T λ are the polarization vectors of the longitudinal and transverse dark photon with helicity index λ = 1, 2, and the thermal width of the dark photon Γ L,T is given by [91] Γ L, In the limit Γ L,T m A , we can approximate the Breit-Wigner distribution by a δ-function With the above approximation, the absorption rate of χ becomes where Note that the inverse decay rates into dark photons of different polarizations are (A11) Thus, we can relate the differential absorption rate and inverse decay rates by With the above equation, the total pair-absorption rate of χ is given by The total pair-absorption rate is equal to that of the inverse decay process χχ → A . Therefore, the inverse decay process χχ → A accounts for all the pair-absorption processes with dark photon resonances when m A > 2m χ . These processes include χχnp → np, χχ → e − e + , χχ → γ * , and χχ → χχ. Note that the dark Compton scattering χA → A χ also admits an on-shell dark fermionχ. With the approximation analogous to Eq. (A3), one can show that the scattering rate of χ via χA → A χ is equivalent to the inverse decay rate of χχ → A , and the scattering rate of A via the same process is equivalent to the decay rate of A → χχ. Thus, we do not include the dark Compton process in our IMFP calculations (See Tables I, II). With Eq. (A12), we can derive the luminosity of χ andχ in the non-diffuse regime where u A is the blackbody energy density of the dark photon. We shall make several comments on the above equation.
1. It is assumed that dark photons are in thermal equilibrium with the SM medium while the dark fermion streams freely. (This corresponds to large and small α D .) In this scenario, Br(A L,T → SM ) ≈ 1. So 1/3 of the dark fermion is produced by the longitudinal dark photon, while 2/3 of the dark fermion is produced by the transverse dark photon. However, if the dark photon escapes the PNS freely, this calculation is invalid, and the dark photon usually cannot decay into the dark fermions before escaping the PNS. Thus the dark fermion luminosity would be negligible. Therefore, we compute the luminosity of χ according to Eq. (D4) for different scenarios.
where Γ L,T = (1 − e −ω/T )Γ L,T A ,abs is the thermal absorptive width of the dark photon [91]. Similarly, for the on-shell transverse photon or plasmon, we replace the subscript "k" with "m" for the effective mixing and the plasma factor. Note that the dispersion relation of the SM photon k 2 = ReΠ L,T is a transcendental equation that can only be solved numerically. The derivation of Eq. (B1) and Eq. (B4) in diagonalized mass basis can be found in Appendix B of [93].
In the medium consisting of relativistic and degenerate electron, the real part of the scalar polarization functions are [92] where v ≡ | k|/ω and the plasma frequency in this limit is With detailed balance Γ L,T prod = e −ω/T Γ L,T abs , the imaginary part of the polarization functions are where Γ L,T prod (Γ L,T abs ) are the production (absorption) rates of the SM photon. When calculating the plasmon decay, we should include the renormalization factor for the polarization vectors of the external photons where the renormalization factor is given by which can be calculated from Eqs. (B5) and (B6) Appendix C: Interaction Rates and Inverse Mean Free Path The interaction rate of particle 1 via a process 1 + ... + n → 1 + ... + m is given by where g i is the degrees of freedom of the initial state particles and f i and f j are the distribution functions of the initial and final state particles. The upper (lower) sign in front of f j is for bosonic (fermionic) final states, which takes into account the bosonic enhancement and Pauli blocking effects. The above expression already assumes that all the particles are in thermal equilibrium. The IMFP of the process is related to the interaction rate bỹ In the following, we denote the 4-momentum of A by k = (ω, k), that of χ by p = (E, p) and that ofχ by p = (E , p ), unless noted otherwise.

A np → np
The absorption rate via inverse nucleon bremsstrahlung A np → np is given by [33] Γ L,T A np→np = 32 3π α e 2 m|L,T n n n p ω 3 where n n,p is the number density of neutrons and protons respectively, m N is the nucleon mass, σ np (T ) is the thermally averaged n-p cross section defined as [31] σ (2) np (T ) = and h L,T is defined as In the derivation of Eq. (C3), the soft radiation approximation is used to connect the bremsstrahlung rate with the experimental n-p scattering cross section. Additionally, the nucleons are assumed to follow Maxwell-Boltzmann distribution, and the Pauli-blocking effect is ignored.
The absorption rate via the Compton-like process A e − → e − γ is approximately [33] where n e is the electron number density, E F = (3π 2 n e ) 2/3 + m 2 e is the electron Fermi energy, and h L,T is defined in Eqs. (C5) and (C6). The approximated formula is valid for ω < ∼ 200 MeV and r ≤ R ν [33].
The decay rates of dark photon via A → e − e + , given that m A > 2m e , are given by where the integral limit is The pair-annihilation rate via inverse nucleon bremsstrahlung χχnp → np is where the momentum transfer is k = p + p and cos θ = k · p/| k|| p|. Note that E, p (E , p ) are for χ (χ). Thus E = E 2 + | k| 2 − 2| k|| p| cos θ and ω = E + E . The two terms F L χχnp→np , F T χχnp→np in the integrand are Similar to the calculation of A np → np, we use soft radiation approximation for the dark fermion pair absorption, and ignore the Pauli-blocking effect of the nucleons. When m A > 2m χ , we use the narrow width approximation (NWA) detailed in Appendix A to approximate the integral. The resulting formula is where the two terms F L * χχnp→np , F T * χχnp→np in the integrand are and we now have | k| = ω 2 − m 2 A . The dark photon widths must include the contribution from the dark photon decay via A → χχ. Thus we have Γ L,T = (1 − e −ω/T )(Γ L,T A ,abs + Γ A →χχ ). Momentum conservation demands that ω is bounded by The pair-annihilation rate via χ(p) +χ(p ) → e − (p 3 ) + e + (p 4 ) is where we have the same definitions for cos θ, E and ω as in Eq. (C11). Moreover, f , f 3 and f 4 are the distribution functions ofχ, e − and e + , respectively. That is, where the positron energy is expressed as E 4 = ω − E 3 by energy conservation. With the electron and positron distribution function included, we take the Pauli-blocking effect into account, which is significant near the center of the PNS. The two terms F L χχ→e − e + , F T χχ→e − e + in the integrand are Momentum conservation demands that E 3 is bounded by When m A > 2m χ , we use NWA (See Appendix A) to approximate the annihilation rate. The resulting formula is where we now have | k| = ω 2 − m 2 A . The energies ofχ and e + on which their distribution functions depend are expressed as E = ω − E and E 4 = ω − E 3 . The two terms F L * χχ→e − e + , F T * χχ→e − e + in the integrand now are where Γ L,T is the same as in Eqs. (C15) and (C16). The limit of E 3 now becomes while the limit of ω is the same as Eq. (C17). 6. χχ → γ * Let ω L,T = ω L,T (| k|) be the energy of the longitudinal and transverse SM photon as functions of | k|, and m L,T be the effective mass satisfying m 2 L,T = ω 2 L,T − | k| 2 = ReΠ L,T (ω L,T , | k|). Then the pair-annihilation rate can be expressed as where g L = 1, g T = 2, and F L χχ→γ * , F T χχ→γ * are defined as where Z L,T is the renormalization factor for the longitudinal plasmon and transverse photon detailed in Appendix B. The bosonic enhancement for the plasmon is included in Eq. (C29). The step function in Eq. (C29) is to ensure that the momenta are kinematically allowed. The definition of cos θ L,T is We use NWA to approximate the annihilation rate when m A > 2m χ . The resulting expression is where ω * L,T is the solution of the equation ReΠ L,T (ω, (C35) ω ± is the same as Eq. (C17) and Γ L,T is the same as in Eqs. (C15) and (C16). 7. χe − → χe − and χp → χp The scattering rate of χ(p) + e − (p 2 ) → χ(p ) + e − (p 4 ) is given by where k = p − p is the momentum transfer, cos θ and E are the same in Eq. (C11) while ω = E − E instead. f 2 , f and f 4 are the distribution functions of the initial state e − , final state χ, and final state e − respectively. The final state electron energy is expressed as E 4 = ω + E 2 by energy conservation. The two terms F L χe − →χe − , F T χe − →χe − are the same as Eqs. (C22) and (C23) except that −E 3 is replaced with +E 2 , and E 2 is bounded below by Note that the momentum transfer is spacelike (k 2 < 0) for the χ-e − scattering. In this scenario, we take the dark photon width to be zero, while Im Π L,T is determined by the imaginary parts of Eqs. (B5) and (B6) with ω → ω + iε for an infinitesimal ε > 0 [94]. The scattering rate via χp → χp is the same as χe − → χe − if the differences in masses and distribution functions of the proton and the electron are taken into account. However, we can approximate the χ-p scattering as an elastic scattering due to the fact that m p is much larger than m χ . We also ignore the Pauli-blocking effect of proton. Thus the scattering rate of χ-p scattering can be approximated by where the momentum transfer is k 2 −| k| 2 −2| p| 2 (1−cos θ). Note that in this static limit (ω → 0), the contribution from the transverse mode, corresponding to the classical magnetic field, is suppressed. And the plasma factor β 2 k|L of the longitudinal mode, corresponding to the static electric field, accounts for the screening effect in medium [95].

A ↔ χχ
For the decay process A → χχ, we take into account the Pauli-blocking effect of the dark fermions. The decay rate is given by where ξ 0 is similar to Eq. (C10) (C42) Note that in the limit T → 0, Eq. (C41) reduces to the decay rate in vacuum given by Eq. (2.6) in [35]. The inverse decay rate of χχ → A , taking into account the bosonic enhancement of the dark photon, is given by where ω ± is the same as Eq. (C17).
Appendix D: Detailed Workflow for Computing DS Luminosity 1. We take the following steps to determine if a DS particle species is in the diffuse regime or not.
(a) Computeλ −1 i (E i , r), the IMFP of each species i assuming all DS particle species are in thermal equilibrium with the SM medium at temperature T (r).
(b) Compute ∆t i , the escape time scale for each particle species i with Eqs. (21) and (22).
(c) Compute N i , the estimated abundance of each particle species with Eq. (20), and obtain the relative abundance η i ≡ N i /N eq i , where N eq i is the abundance of DS particle species i assuming they are in thermal equilibrium with the SM medium. If η i > 1, we simply set η i = 1.
where j runs over all possible DS particle species.
(e) Check the diffuse criteria for the dark photons: If η A L,T = 1 and λ −1 A L,T (R ν ) − λ −1 A L,T →e − e + (R ν ) > R −1 ν , then A L,T is in the diffuse limit. Otherwise, A L,T is treated as non-diffuse particles.
2. Depending on whether the DS particles are in the diffuse limit or not, we compute the luminosity of each particle species as follows.
(a) Regardless of the DS masses, the dark photon luminosity is given by L A L,T = L A L,T ,diff. , for diffuse A L,T (use Eq. (9)), L A L,T ,non−diff. , for non-diffuse A L,T (use Eq. (3)).
(c) When m A > 2m χ , the dark fermion luminosity is given by , if χ is in the diffuse limit, L χ,non−diff. , if only A L , A T are in the diffuse limit, 2 3 L χ,non−diff. , if only A T is in the diffuse limit, 1 3 L χ,non−diff. , if only A L is in the diffuse limit, 0, if no DS particle species is in the diffuse limit.
(D4) (d) The total DS luminosity is L X = L A L + L A T + 2L χ .

Appendix E: Radiative Transfer
We now derive the energy flux carried by the DS particles through a surface of radius r near the neutrinosphere, assuming they are in the diffuse limit. The derivation follows from Appendix I of [96]. The equation of radiative transfer is given by where ρ is the matter density, I is the intensity of radiation per unit solid angle per unit frequency, θ is the angle between the direction of radiation and the radial direction, κ is the opacity, j is the total radiation power emitted per unit mass per unit frequency. We distinguish between the opacity contributions from scattering processes and absorption processes κ = κ s +κ a , and between the radiation by scattering and radiation emitted by matter j = j s +j em . In equilibrium and isotropic environment, where g is the degree of freedom, and the upper (lower) sign is for fermions (bosons). However, in an anisotropic environment, the relation between j em and I is given by j em = κ a 1 ± e −E/T B ∓ κ a e −E/T I.
In the second term, the minus sign for fermion is due to Pauli blocking, while the plus sign for boson is due to stimulated emission. Substitute the above equation back into Eq. (E1) 1 ρ ∂I ∂r cos θ + κ s I − j s + κ * a (I − B) = 0, where κ * a ≡ κ a (1 ± e −E/T ). We assume that the radiation intensity I is very close to I iso,eq . Therefore, we can expand I in terms of Legendre polynomials P n (cos θ) and substitute it back into Eq. (E5). Keeping the terms up to n = 1, we obtain It follows that the energy flux through the spherical surface of radius r is L(r) = 4πr 2 I cos θdEdΩ where m is the mass of the particle carrying the radiation, and λ −1 is the effective IMFP defined as λ −1 = ρ(κ * a + κ s ).

(E9)
We note that since the integral in Eq. (E7) is dominated by the energies with small λ −1 (E, r), it is possible that Eq. (E7) could overestimate the energy flux if λ −1 (E, r) < ∼ R ν for some energies E. We found that it occurs near the switching of diffuse and non-diffuse regimes, and thus it leads to orders of magnitude jumps of energy luminosity. To avoid this caveat, we approximate the IMFP λ −1 (E, r) by the thermally averaged IMFP λ −1 (r) defined as λ −1 (r) = d 3 p f (E, T (r))λ −1 (E, r) d 3 p f (E, T (r)) .
Thus, the energy flux is approximately given by