# Supernova constraint on self-interacting dark sector particles

Allan Sung<EMAIL_ADDRESS>Institute of Physics, Academia Sinica, Taipei,
11529, Taiwan  Department of Physics, Princeton University, Princeton, New
Jersey 08544, USA  Gang Guo<EMAIL_ADDRESS>Institute of Physics, Academia
Sinica, Taipei, 11529, Taiwan  Key Laboratory of Dark Matter and Space
Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing
210033, China  Meng-Ru Wu<EMAIL_ADDRESS>Institute of Physics,
Academia Sinica, Taipei, 11529, Taiwan  Institute of Astronomy and
Astrophysics, Academia Sinica, Taipei, 10617, Taiwan

###### Abstract

We examine the constraints on sub-GeV dark sector particles set by the proto-
neutron star cooling associated with the core-collapse supernova event
SN1987a. Considering explicitly a dark photon portal dark sector model, we
compute the relevant interaction rates of dark photon ($A^{\prime}$) and dark
fermion ($\chi$) with the Standard Model particles as well as their self-
interaction inside the dark sector. We find that even with a small dark sector
fine structure constant $\alpha_{D}\ll 1$, dark sector self-interactions can
easily lead to their own self-trapping. This effect strongly limits the energy
luminosity carried away by dark sector particles from the supernova core and
thus drastically affects the parameter space that can be constrained by
SN1987a. We consider specifically two mass ratios $m_{A^{\prime}}=3m_{\chi}$
and $3m_{A^{\prime}}=m_{\chi}$ which represent scenarios where the decay of
$A^{\prime}$ to $\chi\bar{\chi}$ is allowed or not. For
$m_{A^{\prime}}=3m_{\chi}$, we show that this effect can completely evade the
supernova bounds on widely-examined dark photon parameter space for a dark
sector with $\alpha_{D}\gtrsim 10^{-7}$. In particular, for the mass range
$m_{\chi}\lesssim 20$ MeV, supernova bounds can only be applied to weakly
self-interacting dark sector with $\alpha_{D}\lesssim 10^{-15}$. For
$3m_{A^{\prime}}=m_{\chi}$, bounds in regions where $\alpha_{D}\gtrsim
10^{-7}$ for $m_{\chi}\lesssim 20$ MeV can be evaded similarly. Our findings
thus imply that the existing supernova bounds on light dark particles can be
generally eluded by a similar self-trapping mechanism. This also implies that
nonstandard strongly self-interacting neutrino is not consistent with the
SN1987a observation. Same effects can also take place for other known stellar
bounds on dark sector particles.

##  I Introduction

The detection of 20 electron antineutrinos emitted from the core-collapse
supernova (SN) explosion event, SN1987a Hirata  _et al._ (1987); Bionta  _et
al._ (1987); Alekseev  _et al._ (1988), 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 nonstandard neutrino interactions Sato
and Suzuki (1987); Spergel  _et al._ (1987); Bahcall  _et al._ (1987); Burrows
and Lattimer (1987); Schramm and Truran (1990); Loredo and Lamb (2002);
Frieman  _et al._ (1988); Kolb and Turner (1989); Berezhiani and Smirnov
(1989); Farzan (2003). 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 Raffelt
and Seckel (1988); Turner (1988); Mayle  _et al._ (1988); Brinkmann and Turner
(1988); Janka  _et al._ (1996); Keil  _et al._ (1997); Fischer  _et al._
(2016); Carenza  _et al._ (2019); Bollig  _et al._ (2020); Lucente  _et al._
(2020); Raffelt and Zhou (2011); Argüelles  _et al._ (2019); Suliga  _et al._
(2019); Syvolap  _et al._ (2019); Mastrototaro  _et al._ (2020); Suliga  _et
al._ (2020); Dent  _et al._ (2012); Rrapaj and Reddy (2016); Mahoney  _et al._
(2017); Chang  _et al._ (2017); Hardy and Lasenby (2017); Chang  _et al._
(2018); Guha  _et al._ (2017, 2019); Ishizuka and Yoshimura (1990); Arndt and
Fox (2003); Keung  _et al._ (2014); Tu and Ng (2017); Dev  _et al._ (2020);
Croon  _et al._ (2021); Camalich  _et al._ (2020); Farzan (2003); Hanhart  _et
al._ (2001a, b); Hannestad  _et al._ (2007); Freitas and Wyler (2007), 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_{\nu}\simeq 3\times
10^{53}$ erg (see, however, a caution from Ref. Bar  _et al._ (2020)). In
addition to the SN cooling (more precisely, the PNS cooling) constraint,
recent works also proposed new constraints on light dark photon or dark photon
portal light dark matter (DM) based on other SN-related observables, such as
the measured SN explosion energy Sung  _et al._ (2019), the $\gamma$-rays
Kazanas  _et al._ (2014); DeRocco  _et al._ (2019a), or the produced
(semi-)relativistic dark matter flux arriving at the terrestrial detectors
DeRocco  _et al._ (2019b).

One important aspect in deriving the SN constraint on light dark sector (DS)
particles is that their interaction with SM particles cannot be too strong for
them being trapped inside the PNS. We note that previous studies always
ignored 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 $\sim
O(10^{-41})$ cm2 for neutrinos of $\sim\mathcal{O}(10)$ MeV, dark sector
particles can trap themselves inside the PNS. Consequently, SN bounds on self-
interacting dark sector particles can be largely evaded (see also a very
recent study discussing SN bound on axion-like particle portal light DM Darmé
_et al._ (2020)). 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 core-cusp, too-big-to-fail, and the missing
satellites problems; see, e.g., Refs. Spergel and Steinhardt (2000);
Vogelsberger  _et al._ (2012); Rocha  _et al._ (2013); Tulin and Yu (2018).
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. Tulin  _et al._ (2013); Bernal
_et al._ (2016); Kawasaki  _et al._ (2015); Bringmann  _et al._ (2017); Elor
_et al._ (2016); Ackermann  _et al._ (2015); Abdallah  _et al._ (2016); Feng
_et al._ (2016); Leane  _et al._ (2017); Kouvaris  _et al._ (2016); Essig  _et
al._ (2019); Chang  _et al._ (2019); Depta  _et al._ (2019); Bernal  _et al._
(2020); Foot and Vagnozzi (2015).

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 in a systematic way 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 Mezzacappa  _et al._
(2020) and references therein).111Reference DeRocco  _et al._ (2019b) adopted
a Monte-Carlo based particle transport scheme to compute the light dark matter
flux emitted from the PNS, without considering their potential self-
interactions. 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 nondiffuse 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 nondiffuse 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 cross
sections, the decay rates, and the diffusion luminosity of dark sector
particles are given in the Appendices. We adopt natural units with $\hbar=c=1$
throughout the paper unless explicitly specified.

##  II Models

###  II.1 Dark sector model

Mass | Interaction | Type | Particle | Coupling
---|---|---|---|---
$m_{A^{\prime}}<2m_{\chi}$ | $A^{\prime}np\rightarrow np$ | Abs. | SM | $\epsilon^{2}$
$A^{\prime}e^{-}\rightarrow e^{-}\gamma$ | Abs. | SM | $\epsilon^{2}$
$A^{\prime}\rightarrow e^{-}e^{+}$ | Abs. | SM | $\epsilon^{2}$
$A^{\prime}A^{\prime}\rightarrow\chi\bar{\chi}$ | Abs. | DS | $\alpha_{D}^{2}$
$A^{\prime}\chi\rightarrow\chi A^{\prime}$ | Sca. | DS | $\alpha_{D}^{2}$
$m_{A^{\prime}}>2m_{\chi}$ | $A^{\prime}np\rightarrow np$ | Abs. | SM | $\epsilon^{2}$
$A^{\prime}e^{-}\rightarrow e^{-}\gamma$ | Abs. | SM | $\epsilon^{2}$
$A^{\prime}\rightarrow e^{-}e^{+}$ | Abs. | SM | $\epsilon^{2}$
$A^{\prime}\rightarrow\chi\bar{\chi}$ | Abs. | DS | $\alpha_{D}$
Table 1: 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^{\prime}}>2m_{\chi}$, $A^{\prime}\rightarrow\chi\bar{\chi}$ also accounts for contribution from $A^{\prime}\chi\rightarrow A^{\prime}\chi$ [see Eq. (37) and text below for details].  Mass | Interaction | Type | Particle | Coupling
---|---|---|---|---
$m_{A^{\prime}}<2m_{\chi}$ | $\chi\bar{\chi}np\rightarrow np$ | Abs. | SM | $\epsilon^{2}\alpha_{D}$
$\chi\bar{\chi}\rightarrow e^{-}e^{+}$ | Abs. | SM | $\epsilon^{2}\alpha_{D}$
$\chi\bar{\chi}\rightarrow\gamma^{\ast}$ | Abs. | SM | $\epsilon^{2}\alpha_{D}$
$\chi p\rightarrow\chi p$ | Sca. | SM | $\epsilon^{2}\alpha_{D}$
$\chi e^{-}\rightarrow\chi e^{-}$ | Sca. | SM | $\epsilon^{2}\alpha_{D}$
$\chi\bar{\chi}\rightarrow A^{\prime}A^{\prime}$ | Abs. | DS | $\alpha_{D}^{2}$
$\chi\chi\rightarrow\chi\chi$ | Sca. | DS | $\alpha_{D}^{2}$
$\chi\bar{\chi}\rightarrow\chi\bar{\chi}$ | Sca. | DS | $\alpha_{D}^{2}$
$\chi A^{\prime}\rightarrow A^{\prime}\chi$ | Sca. | DS | $\alpha_{D}^{2}$
$m_{A^{\prime}}>2m_{\chi}$ | $\chi p\rightarrow\chi p$ | Sca. | SM | $\epsilon^{2}\alpha_{D}$
$\chi e^{-}\rightarrow\chi e^{-}$ | Sca. | SM | $\epsilon^{2}\alpha_{D}$
$\chi\bar{\chi}\rightarrow A^{\prime}$ | Abs. | DS | $\alpha_{D}$
Table 2: Relevant processes of dark fermion interactions. Notations are the
same as in Table 1. Only leading-order processes are included here. Note that
for $m_{A^{\prime}}>2m_{\chi}$, $\chi\bar{\chi}\rightarrow A^{\prime}$ also
accounts for DS processes of $A^{\prime}\chi\rightarrow A^{\prime}\chi$ and
$\chi\bar{\chi}\rightarrow\chi\bar{\chi}$, as well as SM processes involving a
pair of dark fermions [see Eq. (37) and text below for details].

We consider a dark photon portal DS model wherein the Dirac dark fermion
$\chi$ couples to dark photon $A^{\prime}$ and the dark photon kinetically
mixes with the SM photon Holdom (1986); Okun (1982); Pospelov  _et al._
(2008). The corresponding Lagrangian of the dark sector is given by

$\begin{split}\mathcal{L}&\supset-\frac{1}{4}F^{\prime}_{\mu\nu}F^{\prime\mu\nu}-\frac{\epsilon}{2}F^{\prime}_{\mu\nu}F^{\mu\nu}+\frac{1}{2}m_{A^{\prime}}^{2}A^{\prime}_{\mu}A^{\prime\mu}\\\
&+\bar{\chi}\left(i\not{\partial}-m_{\chi}\right)\chi+g_{D}\bar{\chi}\not{A^{\prime}}\chi,\end{split}$
(1)

where $\epsilon$ is the mixing parameter, $m_{A^{\prime}}$ is the dark photon
mass, $m_{\chi}$ is the mass of the dark fermion, and $g_{D}$ is the DS
coupling constant. We define the DS fine structure constant $\alpha_{D}\equiv
g_{D}^{2}/4\pi$ analogous to the electromagnetic fine structure constant
$\alpha_{e}\equiv e^{2}/4\pi$.

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. Chang  _et al._ (2017, 2018) 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 $\simeq 30$ MeV and a core mass
density $\gtrsim 10^{14}$ g cm-3. For the dark photon, we include the nucleon-
nucleon bremsstrahlung $np\rightarrow npA^{\prime}$, Compton-like interaction
$\gamma e^{-}\rightarrow e^{-}A^{\prime}$, and electron-positron annihilation
$e^{-}e^{+}\rightarrow A^{\prime}$. For the dark fermion, we consider three
pair-production channels including nucleon-nucleon bremsstrahlung
$np\rightarrow np\chi\bar{\chi}$, electron-positron annihilation
$e^{-}e^{+}\rightarrow\chi\bar{\chi}$, and plasmon decay
$\gamma^{\ast}\rightarrow\chi\bar{\chi}$.

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^{\prime}A^{\prime}\leftrightarrow\chi\bar{\chi}$, the DS Compton scattering
$A^{\prime}\chi\leftrightarrow\chi A^{\prime}$, the scattering of dark
fermions $\chi\chi\rightarrow\chi\chi$ and
$\chi\bar{\chi}\rightarrow\chi\bar{\chi}$, as well as the dark fermion
scattering with SM charged fermions $\chi p\rightarrow\chi p$ and $\chi
e^{-}\rightarrow\chi 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^{\prime}}>2m_{\chi}$, the
(inverse) decay of dark photon $A^{\prime}\leftrightarrow\chi\bar{\chi}$ 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.222$\chi\bar{\chi}\rightarrow\chi\bar{\chi}$ has a dark photon
resonance and $A^{\prime}\chi\rightarrow A^{\prime}\chi$ has a dark fermion
resonance.; see Eq. (37) and text below in the Appendix VI for details. Other
DS self-interactions without resonances are suppressed by an extra factor of
$\alpha_{D}$ compared to the (inverse) decay rate and thus can be neglected
for $m_{A^{\prime}}>2m_{\chi}$. We list all the included interactions in this
work in Tables 1 and 2 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 $\sim 30$ MeV, we only consider DS
particles of masses below $\simeq 1$ GeV.

Figure 1: Supernova profiles (density, temperature, electron chemical
potential, and plasma frequency) used in this work from the $18$ $M_{\odot}$
progenitor model in Ref. Fischer  _et al._ (2010), extracted at the time 1 s
post the core bounce. The vertical dash line indicates where the
neutrinosphere is located.

###  II.2 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 $\sim 10$ s. The
total energy carried away by neutrinos is $\simeq 3\times 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 Raffelt (1990),

$L_{D}\leq L_{\nu}\equiv 3\times 10^{52}~{}{\rm erg~{}s}^{-1},$ (2)

where $L_{D}$ denotes the energy luminosity of dark emission, evaluated during
the early PNS cooling phase, and $L_{\nu}$ 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. Fischer  _et al._ (2010) to compute the
emission of dark photons and dark fermions. This PNS profile was widely used
recently for similar purposes (see e.g., Refs. Chang  _et al._ (2017); Sung
_et al._ (2019)). Note that the choice of a particular SN model may introduce
uncertainties of a factor of a few for the derived bounds Chang  _et al._
(2017); Mahoney  _et al._ (2017). Figure 1 shows the radial evolution of the
density $\rho$, temperature $T$, the electron chemical potential $\mu_{e}$,
and the plasma frequency $\omega_{p}$ [see Eq. (45)]. The position of the
spectral-averaged neutrino decoupling sphere, i.e., neutrinosphere, is
indicated by the vertical dash line at $R_{\nu}\simeq 22$ km. The density
profile shows a monotonically decreasing behavior as a function of radius,
while the temperature profile exhibits a peak at $r\simeq 11$ km, due to the
inefficient compression heating at the densest core region. The electrons are
highly degenerate inside the PNS. The plasma frequency $\omega_{p}\simeq 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 An  _et al._
(2013); Hardy and Lasenby (2017); Chang  _et al._ (2017). We have included
this effect throughout this work and give the details in Appendix VII.

##  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.1) and for cases where they can be considered as diffusive
due to self-trapping (Sec. III.2). We formulate the criteria that determine if
a DS particle species is in diffuse regime or not in Sec. III.3. In Sec.
III.4, 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^{\prime}$ by
$k=(\omega,\vec{k})$, that of $\chi$ by $p=(E,\vec{p})$ and that of
$\bar{\chi}$ by $p^{\prime}=(E^{\prime},\vec{p}^{\prime})$, unless noted
otherwise.

###  III.1 Nondiffuse regime

In the nondiffuse regime, we consider the bulk emission rates of DS particles
inside the neutrinosphere and the attenuation due to absorption and decay,
following Refs. Chang  _et al._ (2017); Sung  _et al._ (2019). The luminosity
of the dark photon is given by

$\begin{split}L_{A^{\prime}}&=\sum_{L,T}\int_{0}^{R_{\nu}}4\pi
r^{2}dr\int\frac{d^{3}k}{(2\pi)^{3}}\\\ &\quad\times
g_{L,T}\omega\Gamma^{L,T}_{A^{\prime},\rm
prod}(\omega,r)e^{-\tau_{L,T}(\omega,r)},\end{split}$ (3)

where $g_{L}=1$, $g_{T}=2$, and $\Gamma^{L,T}_{A^{\prime},\rm prod}$ and
$\tau_{L,T}$ are the production rate and optical depth respectively. The
exponential 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
VII). The production rate $\Gamma^{L,T}_{A^{\prime},\rm prod}$ is determined
by the interactions involving the SM particles listed in Table 1. By detailed
balance, the rates of each production process and its inverse process are
related by $\Gamma^{L,T}_{A^{\prime},\rm
prod}=e^{-\omega/T}\Gamma^{L,T}_{A^{\prime},\rm abs}$, where
$\Gamma^{L,T}_{A^{\prime},\rm abs}$ is the total absorption rate of the
inverse process. Therefore,

$\Gamma^{L,T}_{A^{\prime},\rm
prod}=e^{-\omega/T}(\Gamma^{L,T}_{A^{\prime}np\rightarrow
np}+\Gamma^{L,T}_{A^{\prime}e^{-}\rightarrow
e^{-}\gamma}+\Gamma^{L,T}_{A^{\prime}\rightarrow e^{-}e^{+}}).$ (4)

For the optical depth, we include the absorption and decay of $A^{\prime}$ 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).333The pair-annihilation rate of
$A^{\prime}A^{\prime}\rightarrow\chi\bar{\chi}$ is also ignored since the dark
photon abundance inside the PNS is relatively small compared to that of SM
particles in the nondiffuse regime. This gives

$\tau_{L,T}(\omega,r)=f(r)\int_{r}^{R_{\nu}}\frac{d\tilde{r}}{v}\Gamma_{A^{\prime},\rm
abs}^{L,T}(\omega,\tilde{r}),$ (5)

where $v=|\vec{k}|/\omega$ is the dark photon velocity, and $f(r)$ is a
geometric factor used in Chang  _et al._ (2017) 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 VIII.

Next, we compute the luminosity of the dark fermion ($\chi$) in the nondiffuse
regime as

$L_{\chi}=\int_{0}^{R_{\nu}}4\pi
r^{2}dr\int\frac{d^{3}p}{(2\pi)^{3}}g_{\chi}E\,\Gamma_{\chi,\rm prod}(E,r),$
(6)

where $g_{\chi}=2$ is the physical degrees of freedom of $\chi$, and
$\Gamma_{\chi,\rm prod}$ is the production rate of $\chi$ by the SM medium
(see Table 2). Note that since $\chi$ and $\bar{\chi}$ are symmetric in our
model, the total luminosity in the dark fermion pair is
$L_{\chi}+L_{\bar{\chi}}=2L_{\chi}$. Here we do not include the attenuation
due to pair absorption. This is because in the nondiffuse regime, the dark
fermion abundance below the neutrinosphere is very low, which suppresses the
pair absorption 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 $\Gamma^{\rm
eq}_{\chi,\rm prod}=e^{-E/T}\Gamma^{\rm eq}_{\chi,\rm abs}$. We take the
equilibrium production rates $\Gamma^{\rm eq}_{\chi,\rm prod}$ as an
approximation for the production rates $\Gamma_{\chi,\rm prod}$ used in Eq.
(6). That is,

$\Gamma_{\chi,\rm prod}\simeq e^{-E/T}(\Gamma_{\chi\bar{\chi}np\rightarrow
np}+\Gamma_{\chi\bar{\chi}\rightarrow
e^{-}e^{+}}+\Gamma_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}).$ (7)

This approximation in principle underestimates a bit the production rates of
$\chi$ and $\bar{\chi}$ due to the assumed equilibrium occupation number,
which effectively results in Pauli blocking. However, since the $\chi$ and
$\bar{\chi}$ are always produced in pairs, the effective Pauli-blocking
suppression of the rates is small due to their zero chemical potentials.

We note that when $m_{A^{\prime}}\geq 2m_{\chi}$, 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^{\prime}$ to $\chi$ and
$\bar{\chi}$ when $A^{\prime}$ are in nondiffuse 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 nondiffuse regime and when $m_{A^{\prime}}\geq 2m_{\chi}$. We
also note here that if dark photons are in the diffuse limit but dark fermions
are not, then effectively the non-diffuse 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 VI).444Note that we consider
separately the diffuse condition for the longitudinal and the transverse dark
photons. Thus, we further multiply the nondiffuse 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 IX).

Before we discuss the detailed numerical results of the DS luminosity in the
nondiffuse 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_{A^{\prime}}^{2}/\omega^{2}$ (see Appendix VIII) and by the coupling
constant $\alpha_{D}$, respectively, compared to the production rate of the
transverse dark photons. We can thus approximate the nondiffuse luminosity of
DS particles by considering transverse dark photons only. Second, we consider
for simplicity a homogeneous PNS with radius $R_{c}\simeq 10\mbox{ km}$,
temperature $T\simeq 30\mbox{ MeV}$, density $\rho\simeq 3\times 10^{14}\mbox{
g}/\mbox{cm}^{3}$ and electron fraction $Y_{e}\simeq 0.3$. We also assume that
dark photons are relativistic. Taking the nucleon-nucleon bremsstrahlung
$A^{\prime}np\rightarrow np$ without the plasma effects (see Appendix VIII.1),
the luminosity of DS particles, $L_{D}$, is approximately

$L_{D}\simeq L_{\nu}\times\left(\frac{\epsilon}{4\times
10^{-10}}\right)^{2}\times\exp\left(-\frac{m_{A^{\prime}}}{29\mbox{
MeV}}\right)$ (8)

for $m_{A^{\prime}}\lesssim 1$ GeV.

###  III.2 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 X. The energy flux of a particle species $i$ is
approximately given by

$\begin{split}L_{i}&=-\frac{2g_{i}R_{\nu}^{2}T_{\nu}^{3}}{3\pi}\left.\frac{dT}{dr}\right|_{R_{\nu}}\frac{1}{\langle\lambda^{-1}_{i}(R_{\nu})\rangle}\\\
&\quad\times\int_{m_{i}/T_{\nu}}^{\infty}\xi^{3}\sqrt{\xi^{2}-\left(\frac{m_{i}}{T_{\nu}}\right)^{2}}\frac{e^{\xi}}{\left(e^{\xi}\pm
1\right)^{2}}d\xi,\end{split}$ (9)

where the upper (lower) sign is for fermions (bosons), $g_{i}$ is the physical
degrees of freedom of particle $i$, $T_{\nu}$ is the temperature at the
neutrinosphere, $m_{i}$ is the mass of particle $i$, and
$\langle\lambda^{-1}_{i}(r)\rangle$ is the thermally averaged inverse mean
free path (IMFP)555The IMFP used here is rescaled by the relative abundances
of the DS particles. See Sec. III.3 for the definition. of particle $i$ at
radius $r$ defined by

$\langle\lambda^{-1}_{i}(r)\rangle\equiv\frac{\int
d^{3}p\,f_{i}(E,T(r))\lambda^{-1}_{i}(E,r)}{\int d^{3}p\,f_{i}(E,T(r))},$ (10)

where $f_{i}(E,T(r))$ is the distribution function of particle $i$ at radius
$r$. We distinguish between the absorptive and scattering IMFP,
$\lambda_{i,abs}^{-1}$ and $\lambda_{i,sca}^{-1}$, and define the total IMFP
as

$\begin{split}\lambda_{i}^{-1}(E,r)&\equiv\lambda_{i,abs}^{-1}(E,r)(1\pm
e^{-E/T(r)})\\\ &\quad+\lambda_{i,sca}^{-1}(E,r).\end{split}$ (11)

As in the nondiffuse regime, we compute the energy fluxes of the longitudinal
and transverse dark photon separately, and assume that the energy fluxes of
$\chi$ and $\bar{\chi}$ are equal. Hence, the total energy loss rate in the DS
particles is $L_{tot}=L_{A^{\prime},L}+L_{A^{\prime},T}+2L_{\chi}$. The IMFP
calculations can be found in Appendix VIII.

We now estimate the relevant parameter space region in the limit where the DS
self-interaction is the dominant opacity source. When
$m_{A^{\prime}}>2m_{\chi}$, the dominant DS self-interaction is the (inverse)
decay $A^{\prime}\leftrightarrow\chi\bar{\chi}$. Given $T_{\nu}\simeq
3.9\mbox{ MeV}$ and the temperature gradient $|dT/dr|_{R_{\nu}}\simeq
6.2\times 10^{-4}\mbox{ MeV}/\mbox{m}$, if we fix $m_{A^{\prime}}/m_{\chi}=3$,
the luminosities of the DS particles can be fitted by

$\begin{split}L_{A^{\prime}}&\simeq L_{\nu}\times\left(\frac{6.3\times
10^{-12}}{\alpha_{D}}\right)\left(\frac{\rm MeV}{m_{A^{\prime}}}\right)^{2}\\\
&\quad\times\exp\left[-\sqrt{49.0+\left(\frac{m_{A^{\prime}}}{4.3~{}\rm
MeV}\right)^{2}}\right]\end{split}$ (12)

for $m_{A^{\prime}}\lesssim 25$ MeV, and

$\begin{split}L_{\chi}+L_{\bar{\chi}}&\simeq
L_{\nu}\times\left(\frac{2.1\times 10^{-13}}{\alpha_{D}}\right)\left(\frac{\rm
MeV}{m_{\chi}}\right)^{2}\\\
&\quad\times\exp\left[\left(\frac{m_{\chi}}{1.9~{}\rm
MeV}\right)-\sqrt{49.0+\left(\frac{m_{\chi}}{4.3~{}\rm
MeV}\right)^{2}}\right]\end{split}$ (13)

for $m_{\chi}\lesssim 30$ MeV. If we take $m_{A^{\prime}}=9$ MeV and
$m_{\chi}=3$ MeV, then the total DS luminosity is

$L_{D}\simeq L_{\nu}\times\left(\frac{1.5\times 10^{-16}}{\alpha_{D}}\right).$
(14)

The corresponding $\chi\bar{\chi}\rightarrow A^{\prime}$ cross section
is666The estimated value here is the thermally averaged cross section
evaluated at the center of PNS $\sigma_{\chi\bar{\chi}\rightarrow
A^{\prime}}\equiv\langle\lambda^{-1}_{\chi\bar{\chi}\rightarrow
A^{\prime}}(r=0)\rangle/n_{\chi}(r=0)$.

$\sigma_{\chi\bar{\chi}\rightarrow A^{\prime}}\simeq 4.7\times 10^{-41}\mbox{
cm}^{2}\left(\frac{\alpha_{D}}{1.5\times 10^{-16}}\right).$ (15)

When $m_{A^{\prime}}<2m_{\chi}$, the diffuse luminosities of dark photons and
dark fermions depend on the values of $\epsilon$ and $\alpha_{D}$. For regions
where $\epsilon$ is small enough such that the dominating opacity source for
dark photons is DS self-interaction
$A^{\prime}A^{\prime}\rightarrow\chi\bar{\chi}$, one may estimate the DS
luminosity by considering the dark photon luminosity only. If we fix
$m_{A^{\prime}}/m_{\chi}=1/3$, then the dark photon luminosity can be fitted
by

$\begin{split}L_{A^{\prime}}&\simeq L_{\nu}\times\left(\frac{2.0\times
10^{-7}}{\alpha_{D}}\right)^{2}\\\
&\quad\times\exp\left[\left(\frac{m_{A^{\prime}}}{0.74~{}\rm
MeV}\right)-\sqrt{49.0+\left(\frac{m_{A^{\prime}}}{4.3~{}\rm
MeV}\right)^{2}}\right]\end{split}$ (16)

for $m_{A^{\prime}}\lesssim 25$ MeV. Adopting $m_{A^{\prime}}=3\mbox{ MeV}$
and $m_{\chi}=9\mbox{ MeV}$, this gives

$L_{D}\simeq L_{\nu}\times\left(\frac{4.5\times
10^{-8}}{\alpha_{D}}\right)^{2}.$ (17)

The corresponding $\chi\chi\rightarrow\chi\chi$ cross section is777The cross
section is evaluated in the same way as Eq. (15).

$\sigma_{\chi\chi\rightarrow\chi\chi}\simeq 9.7\times 10^{-37}\mbox{
cm}^{2}\left(\frac{\alpha_{D}}{4.5\times 10^{-8}}\right)^{2},$ (18)

and the corresponding $A^{\prime}A^{\prime}\rightarrow\chi\bar{\chi}$ cross
section is

$\sigma_{A^{\prime}A^{\prime}\rightarrow\chi\bar{\chi}}\simeq 3.1\times
10^{-39}\mbox{ cm}^{2}\left(\frac{\alpha_{D}}{4.5\times 10^{-8}}\right)^{2}.$
(19)

These analytical formulas Eqs. (14)–(19) already illustrate that DS particles
can be self-trapped by interaction cross sections of
$\gtrsim\mathcal{O}(10^{-40})$ cm2 such that the corresponding diffuse DS
luminosity $L_{D}\lesssim L_{\nu}$. Only very weakly interacting DS particle
with $\alpha_{D}$ much smaller than those nominal values in Eqs. (14)–(19) can
result in $L_{D}>L_{\nu}$.

###  III.3 Diffuse criteria

Whether DS particles are in the nondiffuse 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., $\chi\bar{\chi}np\rightarrow np$, $\chi\bar{\chi}\rightarrow
e^{-}e^{+}$, and $\chi\bar{\chi}\rightarrow\gamma^{\ast}$) to the dark fermion
optical depth cannot be neglected. Furthermore, the DS abundances lead to
significant DS self-interactions (e.g. $\chi\chi\rightarrow\chi\chi$,
$\chi\bar{\chi}\rightarrow A^{\prime}A^{\prime}$, $\chi A^{\prime}\rightarrow
A^{\prime}\chi$,…) that can trap themselves in the PNS. These DS self-
interactions 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 self-trapping 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

$\begin{split}N_{i}&=\int_{0}^{R_{\nu}}4\pi
r^{2}dr\int\frac{d^{3}p_{i}}{(2\pi)^{3}}g_{i}\Gamma_{i,\mathrm{prod}}\Delta
t_{i},\\\ \end{split}$ (20)

where the quantity $\Delta 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 $\Delta t_{i}$ as $t_{\rm cool}=1$
s. Second, the shortest possible timescale, i.e., the lower bound of $\Delta
t_{i}$, is estimated by the free-escaping timescale $t_{\rm free}\equiv
R_{\nu}/v_{i}$. Another characteristic timescale here in the limit where a DS
particle is trapped is the diffusion timescale $t_{i,\rm diff}\simeq
R_{\nu}^{2}/D_{i}$ where $D_{i}\simeq v_{i}/\lambda^{-1}_{i}$ is the diffusion
coefficient (see below for details). For $t_{\rm free}<t_{i,\rm diff}<t_{\rm
cool}$, we then take $t_{i,\rm diff}$ as $\Delta t_{i}$. Combining these
criteria, the relevant timescale $\Delta t_{i}$ used in Eq. (20) is given by

$\Delta t_{i}=\begin{cases}t_{\rm free},\,{\rm if}~{}t_{i,\rm diff}\leq t_{\rm
free},\\\ t_{i,\rm diff},\,{\rm if}~{}t_{\rm free}<t_{i,\rm diff}<t_{\rm
cool},\\\ t_{\rm cool},\,{\rm if}~{}t_{\rm cool}\leq t_{i,\rm
diff}.\end{cases}$ (21)

Since the PNS is not homogeneous, the diffusion of DS particles cannot be
described by a constant $D_{i}$. As shown in Tables 1 and 2, 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 timescale.888Overestimating the IMFP makes Eq.
(24) easier to satisfy. Because the luminosity in the diffuse regime is
generally lower than that in the nondiffuse regime for the same parameters,
our estimation leads to a conservative bound. For the SM-type interactions, we
estimate the IMFP at the center ($r_{\rm SM}=0~{}\mathrm{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_{\rm DS}=11\mathrm{km}$ 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 timescale

$t_{i,\rm
diff}\equiv\frac{R^{2}_{\nu}}{v_{i}}\left[\langle\tilde{\lambda}^{-1}_{i,SM}(r_{SM})\rangle+\langle\tilde{\lambda}^{-1}_{i,DS}(r_{DS})\rangle\right]$
(22)

where $\tilde{\lambda}^{-1}_{i,SM}$ and $\tilde{\lambda}^{-1}_{i,DS}$ are SM-
type and DS-type contribution to the IMFP of particle $i$. When computing the
IMFP contribution from DS particles, $\tilde{\lambda}^{-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)–(22) then allows us to compute the DS
abundance $N_{i}$ for particle $i$.

After deriving $N_{i}$, we then define a scaling factor $\eta_{i}\equiv
N_{i}/N_{i}^{\rm eq}$ for each dark particle species $i$ for $N_{i}<N_{i}^{\rm
eq}$ ($\eta_{i}=1$ if $N_{i}\geq N_{\rm eq}$), where

$\begin{split}N_{i}^{\rm eq}\equiv\int_{0}^{R_{\nu}}4\pi
r^{2}dr\int\frac{d^{3}p_{i}}{(2\pi)^{3}}g_{i}f_{i}(E_{i},T(r))\end{split}$
(23)

is the corresponding total equilibrium abundance of $i$. These scaling factors
are then used to better estimate the IMFP, $\lambda^{-1}_{i,DS}$ due to the DS
particles self-interaction, used in Eq. (10) for computing the DS diffuse
luminosities. We also use $\eta_{i}$ to calculate the thermally averaged total
IMFP of $i$ at $R_{\nu}$, $\langle{\lambda}^{-1}_{i}(R_{\nu})\rangle$.

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:

$N_{i}>N_{i}^{\rm eq}{,~{}\rm
and~{}}\langle{\lambda}^{-1}_{i}(R_{\nu})\rangle>R_{\nu}^{-1}.$ (24)

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 $\Delta
t_{i}$. The second condition requires that the IMFP of $i$ at $R_{\nu}$, 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).999Note that in computing the
second criterion, we exclude the contribution from the dark photon decay via
$A^{\prime}\rightarrow 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 $\simeq 30$ MeV, as the electron
chemical potential $\mu_{e}\simeq 15$ MeV at $R_{\nu}$. However, this process
should be strongly Pauli-blocked for $m_{A^{\prime}}\lesssim 200$ MeV in most
region inside the PNS (see Fig. 1) where dark photons are most abundant in the
diffuse limit. When Eq. (24) is satisfied, we use Eq. (9) to compute the DS
diffuse luminosity for $i$. Otherwise, we take the nondiffuse luminosities,
Eqs. (3) and (6) for dark photons and dark fermions, respectively. Note that
in the diffuse regime and when $\Delta t_{i}=t_{i,\rm diff}$, the first
condition in Eq. (24) is approximately equivalent to having the characteristic
thermalization IMFP Keil  _et al._ (2003); Shapiro and Teukolsky (1983) of DS
particles with SM medium, $\sqrt{\lambda^{-1}_{i,SM,\rm
abs}\times\lambda^{-1}_{i,\rm total}}$ averaged over the phase space and
spatial volume inside $R_{\nu}$, being larger than $R_{\nu}^{-1}$. This means
that the energy exchange between the DS particles and the SM medium, enhanced
by the DS self-interactions, is efficient enough to keep themselves in thermal
contact with SM medium.

We provide in Appendix IX a detailed work flow to additionally describe how we
use the results in these sections to compute the dark sector luminosity for
interested readers.

###  III.4 Numerical calculations

Figure 2: Luminosity of different dark sector particles as a function of
$\alpha_{D}$ for a fixed value of $\epsilon=10^{-8}$. In the left (right)
panel, $m_{A^{\prime}}=9\mbox{ MeV}$ and $m_{\chi}=3\mbox{ MeV}$
($m_{A^{\prime}}=3\mbox{ MeV}$ and $m_{\chi}=9\mbox{ MeV}$). The horizontal
dotted lines label the neutrino luminosity $L_{\nu}=3\times 10^{52}\mbox{
erg/s}$ used to set the supernova bound [Eq. (2)]. The vertical dashed lines
indicate the transition from the nondiffuse (smaller $\alpha_{D})$ to the
diffuse (larger $\alpha_{D})$ regime (see text for details). Note that $\chi$
is in the nondiffuse limit for the range of $\alpha_{D}$ shown in the right
panel.

Figure 3: Luminosity of different dark sector particles as a function of
$\epsilon$ for fixed values of $\alpha_{D}=10^{-17}$ (left), $10^{-9}$
(right). In the left (right) panel, $m_{A^{\prime}}=9\mbox{ MeV}$ and
$m_{\chi}=3\mbox{ MeV}$ ($m_{A^{\prime}}=3\mbox{ MeV}$ and $m_{\chi}=9\mbox{
MeV}$). The horizontal dotted lines label the neutrino luminosity
$L_{\nu}=3\times 10^{52}\mbox{ erg/s}$ used to set the supernova bound [Eq.
(2)]. The vertical dashed lines indicate the transition from the nondiffuse
(smaller $\epsilon$) to the diffuse (larger $\epsilon$) regime (see text for
details). Note that $\chi$ is in the nondiffuse limit for the range of
$\epsilon$ shown in the left panel.

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 $\epsilon$ and $\alpha_{D}$. Figure 2 shows the luminosities of
different DS particles as functions of $\alpha_{D}$ with $\epsilon=10^{-8}$
for two different choices of DS masses. In the left panel, the masses are
chosen such that the decay process $A^{\prime}\rightarrow\chi\bar{\chi}$ is
allowed. In this scenario, dark photons of different polarization modes are
trapped diffusively when $\alpha_{D}\gtrsim 5\times 10^{-18}$, while the dark
fermions are diffusive when $\alpha_{D}\gtrsim 10^{-17}$. The nondiffuse
luminosity of the dark photon for small $\alpha_{D}$ mainly depends on the
interaction with the SM particles and thus is independent of
$\alpha_{D}$.101010The difference of $\sim 10^{2.5}$ between the luminosities
of the longitudinal and transverse dark photon is due to the plasma effects..
The dark fermion luminosity is proportional to $\alpha_{D}$ in the diffuse
limit, as analyzed in Appendix VI However, the dark fermion luminosity for
$\alpha_{D}\lesssim 5\times 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.1). In the diffuse regime, the
luminosities of the DS particles are proportional to $\alpha_{D}^{-1}$ due to
the self-trapping interactions effectively dominated by the (inverse) decay
process $A^{\prime}\leftrightarrow\chi\bar{\chi}$.

For the right panel in Fig. 2, we choose DS masses such that the decay process
$A^{\prime}\rightarrow\chi\bar{\chi}$ is not allowed. In this case, dark
photons of different polarization modes are in diffuse regime when
$\alpha_{D}\gtrsim 8\times 10^{-9}$, while dark fermions are always in the
nondiffuse regime for the range of $\alpha_{D}$ shown in the plot. The main
reason leading to several orders of magnitude larger difference in
$\alpha_{D}$ here than the previous case is due to the extra $\alpha_{D}$
dependence in the IMFP of DS self-interactions [see, e.g., Table 1 or Eqs.
(14) and (17)]. Similar to the previous scenario, the nondiffuse luminosity of
the dark photon is independent of $\alpha_{D}$. However, the luminosities in
the diffuse regime scale as $\alpha_{D}^{-2}$. Once again, this is because the
dominant interaction is responsible for trapping the dark photon being
$A^{\prime}A^{\prime}\rightarrow\chi\bar{\chi}$, whose interaction rate is
proportional to $\alpha_{D}^{2}$. The dark fermion luminosity is significantly
smaller than the dark photon luminosity for the range of $\alpha_{D}$ in the
plot, because the dark fermion production rate is proportional to
$\epsilon^{2}\alpha_{D}$, suppressed by an extra factor of $\alpha_{D}$ when
compared with the dark photon production rate.

In Fig. 3, we show the luminosities of the DS particles as functions of
$\epsilon$ with two choices of $\alpha_{D}$ and same DS masses as those in
Fig. 2. For the case with $\alpha_{D}=10^{-17}$ in the left panel where the
decay process $A^{\prime}\rightarrow\chi\bar{\chi}$ is again allowed, the
longitudinal (transverse) dark photons are trapped when $\epsilon\gtrsim
3\times 10^{-9}$ ($\epsilon\gtrsim 2\times 10^{-9}$). For dark fermions, they
just free stream independent of the value of $\epsilon$, because of the small
value of $\alpha_{D}$ (cf., Fig. 2). The nondiffuse luminosities of the dark
photons of both modes are proportional to $\epsilon^{2}$ as determined by
their production rates. The dark fermion luminosity for $\epsilon\lesssim
2\times 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 $\epsilon\gtrsim
10^{-7}$, dark photon–SM interactions dominate over DS self-interactions, so
the dark photon luminosities are proportional to $\epsilon^{-2}$ as determined
by the mean free path of the dark photon–SM interactions. Interestingly, the
dark fermion luminosity becomes independent of $\epsilon$ and can thus remain
larger than $L_{\nu}$ for large $\epsilon$. 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 VI).

For the right panel in Fig. 3, we choose $\alpha_{D}=10^{-9}$ for the case
where $A^{\prime}\rightarrow\chi\bar{\chi}$ is not allowed. There,
longitudinal (transverse) dark photons are diffusively trapped when
$\epsilon\gtrsim 3\times 10^{-7}$ ($\epsilon\gtrsim 2\times 10^{-7}$), while
dark fermions are in the diffuse limit when $\epsilon\gtrsim 2\times 10^{-5}$.
The dark photon luminosities are proportional to $\epsilon^{2}$ for
$\epsilon\lesssim 10^{-9}$, where the optical depths are negligible. The
diffuse luminosities of dark photons scale as $\epsilon^{-2}$ because the
absorption processes by the SM medium are dominant. The dark fermion
luminosity is much smaller than that of the dark photons until
$\epsilon\gtrsim 10^{-5}$, because the dark fermion production rate is
suppressed by an additional factor of $\alpha_{D}$ when compared with the dark
photon production rate. The diffuse luminosity of the dark fermion is
proportional to $\epsilon^{-2}$ for $\epsilon\gtrsim 10^{-4}$, where $\chi
e^{-}\rightarrow\chi e^{-}$ dominates the IMFP at $R_{\nu}$.

Figure 4: Dark sector luminosity $L_{D}$ contours in the
$\alpha_{D}$-$\epsilon$ plane for $m_{A^{\prime}}=9\mbox{ MeV}$,
$m_{\chi}=3\mbox{ MeV}$ (left) and $m_{A^{\prime}}=3\mbox{ MeV}$,
$m_{\chi}=9\mbox{ MeV}$ (right). Solid black (brown) curves indicate where
$L_{D}=L_{\nu}$ ($L_{D}=0.1L_{\nu}$). Regions inside the $L_{D}=L_{\nu}$ are
excluded by the supernova bound. Thin dash lines show where the transition
from the nondiffuse to diffuse regime occur. See text for discussions on
regions labeled (I), (II), and (A), (B).

##  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^{\prime}}$, $m_{\chi}$, $\epsilon$, and
$\alpha_{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^{\prime}\rightarrow\chi\bar{\chi}$ is allowed affects
the DS luminosity significantly. Here we choose two benchmark mass ratios
$m_{A^{\prime}}/m_{\chi}=3$ and $1/3$ to investigate the SN bounds for these
two scenarios.

We first show the contours of the DS luminosity $L_{D}=L_{\nu}$ and
$L_{D}=0.1L_{\nu}$ on the $\alpha_{D}$–$\epsilon$ plane with fixed
$m_{A^{\prime}}$ and $m_{\chi}$ in Fig. 4. The chosen masses are
$m_{A^{\prime}}=9$ MeV and $m_{\chi}=3$ MeV ($m_{A^{\prime}}=3$ MeV and
$m_{\chi}=9$ MeV) for the left (right) panel in the figure, which allows
(forbids) $A^{\prime}\rightarrow\chi\bar{\chi}$. Regions inside the
$L_{D}=L_{\nu}$ contours indicate that the cooling bound [Eq. (2)] is
violated. In the left panel where $A^{\prime}\rightarrow\chi\bar{\chi}$ is
allowed, the excluded region exhibits an L shape, covering the parameter
ranges of (I) $\alpha_{D}\lesssim 10^{-16}$ for
$10^{-10}\lesssim\epsilon\lesssim 10^{-6}$, as well as (II)
$10^{-19}\lesssim\alpha_{D}\lesssim 10^{-16}$ for $\epsilon\gtrsim 10^{-6}$.
For the majority of the parameter space in (I) ($\alpha_{D}<10^{-16}$ and
$10^{-10}<\epsilon<10^{-6}$), the DS luminosity is mainly contributed by the
nondiffuse dark photons (see also Figs. 2 and 3). In (II)
($10^{-19}\lesssim\alpha_{D}\lesssim 10^{-16}$ for $\epsilon\gtrsim 10^{-6}$),
$L_{D}$ is mainly contributed by the nondiffuse dark fermions (see also Fig.
4). Close to the upper edge of the $L_{D}=L_{\nu}$ contour in both regions,
self-trapping of DS particle takes effect such that
$L_{D}\propto\alpha_{D}^{-1}$ decreases with increasing $\alpha_{D}$ (see Fig.
2), leading to the horizontal edge at $\alpha_{D}\sim 10^{-16}$. The lower
bound of $\alpha_{D}\sim 10^{-19}$ for region (II) is due to the inefficient
production of dark fermions. For the right panel in Fig. 4 where
$A^{\prime}\rightarrow\chi\bar{\chi}$ is not allowed, two regions similarly
exist. Region (A) with $10^{-9}\lesssim\epsilon\lesssim 10^{-6}$ and
$\alpha_{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_{\nu}$ at
$\alpha_{D}\sim 10^{-8}$. For the narrow diagonal-shape region (B) at
$\epsilon\gtrsim 10^{-5}$, $L_{D}$ is dominated by dark fermions. This shape
is related to the fact that both the production and pair-absorption rates of
dark fermions from the SM medium are proportional to $\epsilon^{2}\alpha_{D}$,
as discussed also in Sec. III.4.

Figure 5: The $L_{D}=L_{\nu}$ contour on the $\alpha_{D}$-$m_{\chi}$ plane for
a fixed value of $\epsilon=10^{-8}$ and mass ratios
$m_{A^{\prime}}/m_{\chi}=3$ (left), $1/3$ (right). Regions _below_ the
contours are excluded by the SN bound. Above the contours, the self-trapping
of dark sector particles limit their luminosities and cannot be excluded by
the supernova bound.

For both scenarios shown in Fig. 4, there are specific values of $\alpha_{D}$
above which the cooling criterion gives no constraint due to DS self-trapping.
We now investigate the dependence of these critical values of $\alpha_{D}$ on
DS particle masses. Figure 5 shows the excluded regions on the
$\alpha_{D}$-$m_{\chi}$ plane with $\epsilon=10^{-8}$ for fixed mass ratios
$m_{A^{\prime}}=3m_{\chi}$ (left panel) and $3m_{A^{\prime}}=m_{\chi}$ (right
panel). Note that regions below the contours are excluded. In the left panel
where the decay $A^{\prime}\rightarrow\chi\bar{\chi}$ is allowed, the SN bound
only weakly depends on $\alpha_{D}$ for $m_{\chi}\lesssim 10$ MeV. On the
other hand, $\alpha_{D}$ increase sharply with $m_{\chi}$ for $m_{\chi}\gtrsim
10$ MeV. The main reason is the DS particle abundances are insensitive to the
mass for $m_{\chi}\ll 30$ MeV which is the typical temperature inside the PNS.
Larger $m_{\chi}$ (and $m_{A^{\prime}})$ leads to smaller DS abundances, which
in turn gives rise to a smaller IMFP of DS self-trapping for a given
$\alpha_{D}$. Thus, the suppression of $L_{D}$ due to DS self-trapping occurs
at a larger $\alpha_{D}$ for larger $m_{\chi}$ (and $m_{A^{\prime}})$. For
$m_{\chi}\gtrsim 70$ MeV ($m_{A^{\prime}}\gtrsim 210$ MeV), no SN bounds can
be placed due to the inefficient production of DS particles. In the right
panel where $3m_{A^{\prime}}=m_{\chi}$, the bound is completely determined by
the dark photon only because the dark fermion production is further suppressed
by an extra factor of $\alpha_{D}$. Here, the critical $\alpha_{D}$ also
increases with $m_{\chi}$ sharply for $m_{\chi}\gtrsim 10$ MeV as in the left
panel. Once again, this is resulting from smaller dark photon abundance for
larger $m_{\chi}$ (thus $m_{A}^{\prime}$) inside the PNS. It thus requires a
larger $\alpha_{D}$ for $A^{\prime}A^{\prime}\rightarrow\chi\bar{\chi}$ to
self-trap the dark photons. Note that the maximal $m_{\chi}\simeq 800$ MeV
below which SN bound exists corresponds to a maximal $m_{A}^{\prime}\simeq
250$ MeV, consistent with the maximal $m_{A}^{\prime}$ in the left panel.

Figure 6: The excluded regions (enclosed by the solid curves) on the
$\epsilon$-$m_{A^{\prime}}$ plane for different $\alpha_{D}$ values, with
fixed mass ratios $m_{A^{\prime}}/m_{\chi}=3$ (left), $1/3$ (right). Note that
we include bounds from considering the degree of freedom of dark photon only
(dotted black curves) for comparison. The thin gray curve in the right panel
shows regions excluded by the terrestrial experiments extracted from Andreas
_et al._ (2012). These curves show that even with tiny values of $\alpha_{D}$,
supernova bounds on dark photon parameter space can be largely affected.

The above examples clearly demonstrated how efficiently a small $\alpha_{D}$
within the self-interacting DS can affect the SN bound. In Fig. 6, we further
show the excluded regions in the $\epsilon$-$m_{A^{\prime}}$ plane with two
different values of $\alpha_{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. Andreas  _et al._ (2012)). We refrain from showing other
astrophysical or cosmological bounds on dark photon of similar masses; see
e.g., Refs. Jaeckel  _et al._ (2008); Mirizzi  _et al._ (2009); Redondo and
Raffelt (2013); An  _et al._ (2013); Fradette  _et al._ (2014); Sung  _et al._
(2019); DeRocco  _et al._ (2019a); An  _et al._ (2020); Li  _et al._ (2020);
Sieverding  _et al._ (2021). These plots show once again that even for small
values of $\alpha_{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 $\alpha_{D}$, e.g.,
$10^{-17}$ and $10^{-9}$ shown in the left and right panels, respectively, the
cooling bounds can be extended to larger $\epsilon$. This is due to the
contribution of dark fermions through the decay of trapped dark photons, as
discussed in earlier sections. Note that for $\alpha_{D}=10^{-12}$ in the left
panel, the excluded region slightly extends to larger values of
$m_{A^{\prime}}$ for $\epsilon\gtrsim 10^{-7}$. This is because the dark
fermion luminosity depends on both the decay rate of
$A^{\prime}\rightarrow\chi\bar{\chi}$ and the blackbody energy density of dark
photons (see Appendix VI). The larger decay rate can compensate the smaller
dark photon energy density with increasing $m_{A^{\prime}}$, 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 self-interacting 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 1 and 2. We then used these to compute the energy luminosities of dark
sector particles in the nondiffuse 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 suppress 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^{\prime}}/m_{\chi}=3$ and $1/3$, which
represent scenarios where $A^{\prime}\rightarrow\chi\bar{\chi}$ is allowed or
not. For the former (later) case with $m_{A^{\prime}}/m_{\chi}=3$
($m_{A^{\prime}}/m_{\chi}=1/3$), SN bounds only apply to weakly interacting
dark sectors whose dark fine structure constant $\alpha_{D}\lesssim 10^{-15}$
($\lesssim 10^{-7}$), for $m_{\chi}\lesssim\mathcal{O}(20)$ MeV (see Figs. 4
and 5). The dominating dark sector cross sections for these $\alpha_{D}$
values correspond to $\simeq 10^{-40}$ cm2. In particular, there is no SN
bound for the former case for $\alpha_{D}\gtrsim 10^{-7}$ (Fig. 5). Our
results differ from the previous analysis Chang  _et al._ (2018) considering
similar models. This is because Ref. Chang  _et al._ (2018) assumed that the
DS particles decouple from the SM medium at the surface where the $\chi$-p
scattering becomes inefficient, and ignored the DS self-interactions which can
trap themselves and help thermalize the DS particles with the SM medium.

Although the exact excluded regions in the DS parameter space should also
depend on the chosen value of $m_{A^{\prime}}/m_{\chi}$, 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 nonstandard
strongly self-interacting neutrinos proposed to resolve the Hubble tension
Kreisch  _et al._ (2020), the needed strong self-interaction, $\sim 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 Ellis  _et al._ (2018); Nelson  _et al._
(2019); Bauswein  _et al._ (2020). 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 cannot be excluded by the SN cooling bound, as the
required $\chi$-$\chi$ self-interaction cross section
$\sigma_{\chi\chi}/m_{\chi}\sim\mathcal{O}(1)$ cm2 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 nondiffuse 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 nondiffuse
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 multidimensional 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.

###### Acknowledgements.

We thank Jae Hyeok Chang, Yen-Hsun Lin, Huitzu Tu, and Tse-Chun Wang for
useful discussions. We also thank Luc Darme, Sam McDermott, Sunny Vagnozzi,
and an anonymous referee for providing feedback during the reviewing process,
which help improve this work. This work was supported in part by the Ministry
of Science and Technology, Taiwan under Grant No. 108-2112-M-001-010, No.
109-2112-M-001-004, the Academia Sinica under Project No. AS-CDA-109-M11, and
the Physics Division of National Center for Theoretical Sciences.

##  VI Narrow Width Approximation

Consider a process $\chi\bar{\chi}\rightarrow A^{\prime\ast}\rightarrow
f_{1}f_{2}...$ when $m_{A^{\prime}}>2m_{\chi}$ so that the intermediate
$A^{\prime}$ can be on shell. Let $k=(\omega,\vec{k})$ be the momentum of the
intermediate $A^{\prime}$ and $p=(E,\vec{p})$,
$p^{\prime}=(E^{\prime},\vec{p}^{\prime})$ be the momentum of the initial
$\chi$ and $\bar{\chi}$ respectively. The spin-averaged amplitude squared is
given by

$\begin{split}\overline{|\mathcal{M}|^{2}}=g_{D}^{2}\cdot\frac{k^{2}}{2}&\left\\{\frac{|\epsilon_{L}\cdot
J|^{2}}{(k^{2}-m_{A^{\prime}}^{2})^{2}+(\omega\Gamma_{L})^{2}}\left[1-\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right]\right.\\\
&\left.+\frac{\sum_{\lambda}|\epsilon_{T\lambda}\cdot
J|^{2}}{(k^{2}-m_{A^{\prime}}^{2})^{2}+(\omega\Gamma_{T})^{2}}\left[\frac{1}{2}+\frac{2m_{\chi}^{2}}{k^{2}}+\frac{1}{2}\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right]\right\\},\\\
\end{split}$ (25)

where $J^{\mu}$ is the final state current that couples to the dark photon,
$\epsilon_{L}$, $\epsilon_{T\lambda}$ are the polarization vectors of the
longitudinal and transverse dark photon with helicity index $\lambda=1,2$, and
the thermal width of the dark photon $\Gamma_{L,T}$ is given by Weldon (1983)

$\Gamma_{L,T}=(1-e^{-\omega/T})(\Gamma^{L,T}_{A^{\prime},\rm
abs}+\Gamma_{A^{\prime}\rightarrow\chi\bar{\chi}}).$ (26)

In the limit $\Gamma_{L,T}\ll m_{A^{\prime}}$, we can approximate the Breit-
Wigner distribution by a $\delta$ function

$\frac{1}{(k^{2}-m_{A^{\prime}}^{2})^{2}+(\omega\Gamma_{L,T})^{2}}\rightarrow\frac{\pi}{\omega\Gamma_{L,T}}\delta(k^{2}-m_{A^{\prime}}^{2}).$
(27)

With the above approximation, the absorption rate of $\chi$ becomes

$\begin{split}\Gamma_{\chi\bar{\chi}\rightarrow A^{\prime\ast}\rightarrow
f}&\equiv\frac{1}{2E}\int\frac{d^{3}p^{\prime}}{(2\pi)^{3}}\frac{g_{\bar{\chi}}f_{\bar{\chi}}}{2E^{\prime}}d\Pi_{f}(2\pi)^{4}\delta^{4}(p+p^{\prime}-p_{f})\overline{|\mathcal{M}|^{2}}\\\
&\simeq\frac{1}{16\pi
E|\vec{p}|}\int_{\omega_{-}}^{\omega_{+}}d\omega(1+f_{A^{\prime}})g_{\bar{\chi}}f_{\bar{\chi}}\left[\mbox{Br}(A^{\prime}_{L}\rightarrow
f)\overline{|\mathcal{M}|^{2}}_{\chi\bar{\chi}\rightarrow
A^{\prime}_{L}}+\mbox{Br}(A^{\prime}_{T}\rightarrow
f)\overline{|\mathcal{M}|^{2}}_{\chi\bar{\chi}\rightarrow
A^{\prime}_{T}}\right],\end{split}$ (28)

where

$\displaystyle\omega_{\pm}=\frac{m_{A^{\prime}}^{2}}{2m_{\chi}^{2}}\left(E\pm|\vec{p}|\sqrt{1-\frac{4m_{\chi}^{2}}{m_{A^{\prime}}^{2}}}\right),$
(29) $\displaystyle\mbox{Br}(A^{\prime}_{L,T}\rightarrow
f)=\frac{\Gamma_{A^{\prime}\rightarrow f}}{\Gamma^{L,T}_{A^{\prime},\rm
abs}+\Gamma_{A^{\prime}\rightarrow\chi\bar{\chi}}},$ (30)
$\displaystyle\Gamma_{A^{\prime}_{L}\rightarrow f}=\frac{1}{2\omega}\int
d\Pi_{f}(2\pi)^{4}\delta^{4}(k-p_{f})|\epsilon_{L}\cdot J|^{2},$ (31)
$\displaystyle\Gamma_{A^{\prime}_{T}\rightarrow f}=\frac{1}{2\omega}\int
d\Pi_{f}(2\pi)^{4}\delta^{4}(k-p_{f})\frac{1}{2}\sum_{\lambda}|\epsilon_{T\lambda}\cdot
J|^{2},$ (32)
$\displaystyle\overline{|\mathcal{M}|^{2}}_{\chi\bar{\chi}\rightarrow
A^{\prime}_{L}}=\frac{1}{2}g_{D}^{2}m_{A^{\prime}}^{2}\left[1-\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right],$
(33) $\displaystyle\overline{|\mathcal{M}|^{2}}_{\chi\bar{\chi}\rightarrow
A^{\prime}_{T}}=g_{D}^{2}m_{A^{\prime}}^{2}\left[\frac{1}{2}+\frac{2m_{\chi}^{2}}{m_{A^{\prime}}^{2}}+\frac{1}{2}\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right].$
(34)

Note that the inverse decay rates into dark photons of different polarizations
are

$\Gamma_{\chi\bar{\chi}\rightarrow A^{\prime}_{L,T}}=\frac{1}{16\pi
E|\vec{p}|}\int_{\omega_{-}}^{\omega_{+}}d\omega(1+f_{A^{\prime}})g_{\bar{\chi}}f_{\bar{\chi}}\overline{|\mathcal{M}|^{2}}_{\chi\bar{\chi}\rightarrow
A^{\prime}_{L,T}}.$ (35)

Thus, we can relate the differential absorption rate and inverse decay rates
by

$\frac{d\Gamma_{\chi\bar{\chi}\rightarrow A^{\prime\ast}\rightarrow
f}}{d\omega}\simeq\mbox{Br}(A^{\prime}_{L}\rightarrow
f)\frac{d\Gamma_{\chi\bar{\chi}\rightarrow
A^{\prime}_{L}}}{d\omega}+\mbox{Br}(A^{\prime}_{T}\rightarrow
f)\frac{d\Gamma_{\chi\bar{\chi}\rightarrow A^{\prime}_{T}}}{d\omega}.$ (36)

With the above equation, the total pair-absorption rate of $\chi$ is given by

$\begin{split}\Gamma_{\chi,\rm
pair}&=\sum_{f}\Gamma_{\chi\bar{\chi}\rightarrow A^{\prime\ast}\rightarrow
f}\\\
&=\int_{\omega_{-}}^{\omega_{+}}d\omega\,\left[\frac{d\Gamma_{\chi\bar{\chi}\rightarrow
A^{\prime}_{L}}}{d\omega}\sum_{f}\mbox{Br}(A^{\prime}_{L}\rightarrow
f)+\frac{d\Gamma_{\chi\bar{\chi}\rightarrow
A^{\prime}_{T}}}{d\omega}\sum_{f}\mbox{Br}(A^{\prime}_{T}\rightarrow
f)\right]\\\
&=\int_{\omega_{-}}^{\omega_{+}}d\omega\,\left[\frac{d\Gamma_{\chi\bar{\chi}\rightarrow
A^{\prime}_{L}}}{d\omega}+\frac{d\Gamma_{\chi\bar{\chi}\rightarrow
A^{\prime}_{T}}}{d\omega}\right]\\\ &=\Gamma_{\chi\bar{\chi}\rightarrow
A^{\prime}}.\end{split}$ (37)

The total pair-absorption rate is equal to that of the inverse decay process
$\chi\bar{\chi}\rightarrow A^{\prime}$. Therefore, the inverse decay process
$\chi\bar{\chi}\rightarrow A^{\prime}$ accounts for all the pair-absorption
processes with dark photon resonances when $m_{A^{\prime}}>2m_{\chi}$. These
processes include $\chi\bar{\chi}np\rightarrow np$, $\chi\bar{\chi}\rightarrow
e^{-}e^{+}$, $\chi\bar{\chi}\rightarrow\gamma^{\ast}$, and
$\chi\bar{\chi}\rightarrow\chi\bar{\chi}$. Note that the dark Compton
scattering $\chi A^{\prime}\rightarrow A^{\prime}\chi$ also admits an on-shell
dark fermion $\bar{\chi}$. With the approximation analogous to Eq. (27), one
can show that the scattering rate of $\chi$ via $\chi A^{\prime}\rightarrow
A^{\prime}\chi$ is equivalent to the inverse decay rate of
$\chi\bar{\chi}\rightarrow A^{\prime}$, and the scattering rate of
$A^{\prime}$ via the same process is equivalent to the decay rate of
$A^{\prime}\rightarrow\chi\bar{\chi}$. Thus, we do not include the dark
Compton process in our IMFP calculations (see Tables 1 and 2).

With Eq. (36), we can derive the luminosity of $\chi$ and $\bar{\chi}$ in the
nondiffuse regime

$L_{\chi}+L_{\bar{\chi}}\simeq\int dV\int
d\omega\frac{du_{A^{\prime}}}{d\omega}\Gamma_{A^{\prime}\rightarrow\chi\bar{\chi}}\left[\frac{1}{3}\mbox{Br}(A^{\prime}_{L}\rightarrow
f)+\frac{2}{3}\mbox{Br}(A^{\prime}_{T}\rightarrow f)\right],$ (38)

where $u_{A^{\prime}}$ is the blackbody energy density of the dark photon. We
shall make several comments on the above equation.

  1. 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 $\epsilon$
and small $\alpha_{D}$.) In this scenario,
$\mbox{Br}(A^{\prime}_{L,T}\rightarrow SM)\approx 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 $\chi$ according to Eq. (97) for different scenarios.

  2. 2.

In the limit of large $\epsilon$ and small $\alpha_{D}$, the dark sector
luminosity is mostly contributed by the dark fermion since the dark photon is
trapped. The luminosity only depends on $\alpha_{D}$ and does not depend on
$\epsilon$ in this limit. It is possible that the luminosity can exceed
$L_{\nu}$ for some values of $\alpha_{D}$ even if $\epsilon$ is large, as
shown in the left panels of Figs. 3, 4, and 6.

  3. 3.

The luminosity in this regime also depends on $m_{A^{\prime}}$ through the
blackbody energy density $u_{A^{\prime}}$. Therefore, for a larger value of
$\alpha_{D}$, the constraint can extends to larger $m_{A^{\prime}}$ values, as
shown in the left panel of Fig. 6.

##  VII Plasma Effects

In the PNS, the hot ($T\sim 30$ MeV) and dense ($\rho\sim 10^{14}\mbox{
g}/\mbox{cm}^{3}$) plasma consisting of electrons and nucleons can modify the
dispersion relations of the SM photons and give rise to a longitudinally
polarized propagation mode called plasmon. This plasma effect changes the SM
photon propagator, so we must take it into account when calculating the
interaction rates between SM and DS particles through photon–dark photon
mixing. The effective mixing $\epsilon_{k|L,T}$ and the plasma factor
$\beta_{k|L,T}$ for the SM current–dark photon interactions are defined by

$\beta_{k|L,T}^{2}\equiv\frac{\epsilon_{k|L,T}^{2}}{\epsilon^{2}}\equiv\frac{(k^{2})^{2}}{(k^{2}-\mbox{Re}\Pi_{L,T})^{2}+(\mbox{Im}\Pi_{L,T})^{2}},$
(39)

where $k=(\omega,\vec{k})$ is the 4-momentum of the dark photon with
$k^{2}\equiv\omega^{2}-|\vec{k}|^{2}$, and $\Pi_{L,T}$ is the longitudinal (L)
and transverse (T) polarization functions and is related to the SM photon
polarization tensor $\Pi^{\mu\nu}=e^{2}\langle A^{\mu},A^{\nu}\rangle$ by

$\displaystyle\Pi_{L}$ $\displaystyle=\frac{k^{2}}{|\vec{k}|^{2}}\Pi^{00},$
(40) $\displaystyle\Pi_{T}$
$\displaystyle=\frac{1}{2}\left(\delta^{ij}-\frac{k^{i}k^{j}}{|\vec{k}|^{2}}\right)\Pi^{ij}.$
(41)

We would like to point out that our definition for $\Pi_{L}$ aligns with An
_et al._ (2013) but differs from Braaten and Segel (1993). For the on-shell
dark photon, we replace the subscript “$k$” with “$m$” for the effective
mixing, i.e.,
$\epsilon_{m|L,T}^{2}=(\epsilon_{k|L,T}^{2})|_{k^{2}=m_{A^{\prime}}^{2}}$, and
likewise for the plasma factor. On the other hand, the DS current–SM photon
counterpart of Eq. (39) is

$\frac{\bar{\epsilon}_{k|L,T}^{2}}{\epsilon^{2}}=\frac{(k^{2})^{2}}{(k^{2}-m_{A^{\prime}}^{2})^{2}+(\omega\Gamma_{L,T})^{2}}\equiv\bar{\beta}_{k|L,T}^{2},$
(42)

where $\Gamma_{L,T}=(1-e^{-\omega/T})\Gamma^{L,T}_{A^{\prime},\rm abs}$ is the
thermal absorptive width of the dark photon Weldon (1983). 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}=\mathrm{Re}\Pi_{L,T}$ is a transcendental
equation that can only be solved numerically. The derivation of Eqs. (39) and
(42) in diagonalized mass basis can be found in Appendix B of Lin (2019).

In the medium consisting of relativistic and degenerate electron, the real
parts of the scalar polarization functions are Braaten and Segel (1993)

$\displaystyle\mbox{Re}\Pi_{L}$
$\displaystyle=\frac{3\omega_{p}^{2}(1-v^{2})}{v^{2}}\left[\frac{1}{2v}\ln{\left(\frac{1+v}{1-v}\right)}-1\right],$
(43) $\displaystyle\mbox{Re}\Pi_{T}$
$\displaystyle=\frac{3\omega_{p}^{2}}{2v^{2}}\left[1-\frac{1-v^{2}}{2v}\ln{\left(\frac{1+v}{1-v}\right)}\right],$
(44)

where $v\equiv|\vec{k}|/\omega$ and the plasma frequency in this limit is

$\omega_{p}^{2}=\frac{4\alpha_{e}}{3\pi}\left(\mu_{e}^{2}+\frac{1}{3}\pi^{2}T^{2}\right).$
(45)

With detailed balance $\Gamma^{L,T}_{\rm prod}=e^{-\omega/T}\Gamma^{L,T}_{\rm
abs}$, the imaginary part of the polarization functions is

$\mbox{Im}\Pi_{L,T}=-\omega(\Gamma^{L,T}_{\rm prod}-\Gamma^{L,T}_{\rm
abs})=-\omega\Gamma^{L,T}_{\rm abs}(1-e^{-\omega/T}),$ (46)

where $\Gamma^{L,T}_{\rm prod}$ ($\Gamma^{L,T}_{\rm 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,

$\tilde{\epsilon}_{L,T}^{\mu}=\sqrt{Z_{L,T}}\epsilon_{L,T}^{\mu},$ (47)

where the renormalization factor is given by

$Z_{L,T}^{-1}=1-\left.\frac{\partial\mbox{Re}\Pi_{L,T}}{\partial\omega^{2}}\right|_{\rm
pole},$ (48)

which can be calculated from Eqs. (43) and (44),

$\displaystyle Z_{L}^{-1}$
$\displaystyle=\frac{3\omega_{p}^{2}}{2\omega^{2}v^{2}}\left[1-\frac{1-v^{2}}{2v}\ln{\left(\frac{1+v}{1-v}\right)}\right]_{\rm
pole},$ (49) $\displaystyle Z_{T}^{-1}$
$\displaystyle=1-\frac{3\omega_{p}^{2}}{2\omega^{2}v^{2}}\left[\frac{3}{2}-\frac{3-v^{2}}{4v}\ln{\left(\frac{1+v}{1-v}\right)}\right]_{\rm
pole}.$ (50)

##  VIII Interaction Rates and Inverse Mean Free Path

The interaction rate of particle $1$ via a process $1+...+n\rightarrow
1^{\prime}+...+m^{\prime}$ is given by

$\Gamma=\frac{1}{2E_{1}}\int\prod_{i=1}^{n}\frac{d^{3}p_{i}}{(2\pi)^{3}}\frac{g_{i}f_{i}}{2E_{i}}\prod_{j=1^{\prime}}^{m^{\prime}}\frac{d^{3}p_{j}}{(2\pi)^{3}}\frac{1\pm
f_{j}}{2E_{j}}\overline{|\mathcal{M}|^{2}}(2\pi)^{4}\delta^{4}(\sum_{i=1}^{n}p_{i}-\sum_{j=1^{\prime}}^{m^{\prime}}p_{j}),$
(51)

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 by

$\tilde{\lambda}^{-1}=\frac{\Gamma}{v_{1}}=\frac{E_{1}}{|\vec{p}_{1}|}\Gamma.$
(52)

In the following, we denote the 4-momentum of $A^{\prime}$ by
$k=(\omega,\vec{k})$, that of $\chi$ by $p=(E,\vec{p})$, and that of
$\bar{\chi}$ by $p^{\prime}=(E^{\prime},\vec{p}^{\prime})$, unless noted
otherwise.

###  VIII.1 $A^{\prime}np\rightarrow np$

The absorption rate via inverse nucleon bremsstrahlung
$A^{\prime}np\rightarrow np$ is given by Chang  _et al._ (2017)

$\Gamma^{L,T}_{A^{\prime}np\rightarrow
np}=\frac{32}{3\pi}\frac{\alpha_{e}\epsilon^{2}_{m|L,T}n_{n}n_{p}}{\omega^{3}}\left(\frac{\pi
T}{m_{N}}\right)^{3/2}\left\langle\sigma_{np}^{(2)}(T)\right\rangle h_{L,T},$
(53)

where $n_{n,p}$ is the number density of neutrons and protons, respectively,
$m_{N}$ is the nucleon mass, $\left\langle\sigma_{np}^{(2)}(T)\right\rangle$
is the thermally averaged $n$-$p$ cross section defined as Rrapaj and Reddy
(2016)

$\left\langle\sigma^{(2)}_{np}(T)\right\rangle=\frac{1}{2}\int_{0}^{\infty}dxx^{2}e^{-x}\int_{-1}^{1}d\cos{\theta_{cm}}(1-\cos{\theta_{cm}})\frac{d\sigma_{np}}{d\cos{\theta_{cm}}}\left(T_{cm}=xT,\cos{\theta_{cm}}\right),$
(54)

and $h_{L,T}$ is defined as

$\displaystyle h_{L}$ $\displaystyle=\frac{m_{A^{\prime}}^{2}}{\omega^{2}},$
(55) $\displaystyle h_{T}$ $\displaystyle=1.$ (56)

In the derivation of Eq. (53), 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.

###  VIII.2 $A^{\prime}e^{-}\rightarrow e^{-}\gamma$

The absorption rate via the Compton-like process $A^{\prime}e^{-}\rightarrow
e^{-}\gamma$ is approximately Chang  _et al._ (2017)

$\Gamma^{L,T}_{A^{\prime}e^{-}\rightarrow
e^{-}\gamma}=\frac{8\pi\alpha_{e}^{2}\epsilon^{2}_{m|L,T}n_{e}}{3E_{F}^{2}}\sqrt{\frac{\omega_{p}}{\omega}}h_{L,T},$
(57)

where $n_{e}$ is the electron number density,
$E_{F}=\sqrt{(3\pi^{2}n_{e})^{2/3}+m_{e}^{2}}$ is the electron Fermi energy,
and $h_{L,T}$ is defined in Eqs. (55) and (56). The approximated formula is
valid for $\omega\lesssim 200$ MeV and $r\leq R_{\nu}$ Chang  _et al._ (2017).

###  VIII.3 $A^{\prime}\rightarrow e^{-}e^{+}$

The decay rates of dark photon via $A^{\prime}\rightarrow e^{-}e^{+}$, given
that $m_{A^{\prime}}>2m_{e}$, are given by

$\displaystyle\Gamma_{A^{\prime}\rightarrow e^{-}e^{+}}^{L}$
$\displaystyle=\frac{\epsilon_{m|L}^{2}\alpha_{e}m_{A^{\prime}}^{2}}{4|\vec{k}|}\int_{-\xi_{0}}^{\xi_{0}}d\xi\left[1+\exp{\left(\frac{\mu_{e}-\omega/2}{T}-\frac{\omega}{2T}\xi\right)}\right]^{-1}\left(1-\frac{\omega^{2}}{|\vec{k}|^{2}}\xi^{2}\right),$
(58) $\displaystyle\Gamma_{A^{\prime}\rightarrow e^{-}e^{+}}^{T}$
$\displaystyle=\frac{\epsilon_{m|T}^{2}\alpha_{e}m_{A^{\prime}}^{2}}{4|\vec{k}|}\int_{-\xi_{0}}^{\xi_{0}}d\xi\left[1+\exp{\left(\frac{\mu_{e}-\omega/2}{T}-\frac{\omega}{2T}\xi\right)}\right]^{-1}\left(\frac{1}{2}+\frac{2m_{e}^{2}}{m_{A^{\prime}}^{2}}+\frac{\omega^{2}}{2|\vec{k}|^{2}}\xi^{2}\right),$
(59)

where the integral limit is

$\xi_{0}=\frac{|\vec{k}|}{\omega}\sqrt{1-\frac{4m_{e}^{2}}{m_{A^{\prime}}^{2}}}.$
(60)

###  VIII.4 $\chi\bar{\chi}np\rightarrow np$

The pair-annihilation rate via inverse nucleon bremsstrahlung
$\chi\bar{\chi}np\rightarrow np$ is

$\Gamma_{\chi\bar{\chi}np\rightarrow
np}=\frac{16\epsilon^{2}\alpha_{D}\alpha_{e}n_{n}n_{p}}{3\pi^{2}E}\left(\frac{\pi
T}{m_{N}}\right)^{3/2}\left\langle\sigma_{np}^{(2)}(T)\right\rangle\int_{0}^{\infty}d|\vec{k}|\int_{-1}^{1}d\cos{\theta}\,\frac{|\vec{k}|^{2}(\mathcal{F}^{L}_{\chi\bar{\chi}np\rightarrow
np}+2\mathcal{F}^{T}_{\chi\bar{\chi}np\rightarrow
np})}{E^{\prime}\omega^{2}k^{2}(1+e^{E^{\prime}/T})},$ (61)

where the momentum transfer is $k=p+p^{\prime}$ and
$\cos{\theta}=\vec{k}\cdot\vec{p}/|\vec{k}||\vec{p}|$. Note that $E,p$
($E^{\prime},p^{\prime}$) are for $\chi$ ($\bar{\chi}$). Thus
$E^{\prime}=\sqrt{E^{2}+|\vec{k}|^{2}-2|\vec{k}||\vec{p}|\cos{\theta}}$ and
$\omega=E+E^{\prime}$. The two terms
$\mathcal{F}^{L}_{\chi\bar{\chi}np\rightarrow np}$,
$\mathcal{F}^{T}_{\chi\bar{\chi}np\rightarrow np}$ in the integrand are

$\displaystyle\mathcal{F}^{L}_{\chi\bar{\chi}np\rightarrow np}$
$\displaystyle=\beta_{k|L}^{2}\bar{\beta}_{k|L}^{2}\cdot\frac{k^{2}}{\omega^{2}}\left[1-\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right],$
(62) $\displaystyle\mathcal{F}^{T}_{\chi\bar{\chi}np\rightarrow np}$
$\displaystyle=\beta_{k|T}^{2}\bar{\beta}_{k|T}^{2}\cdot\left[\frac{1}{2}+\frac{2m_{\chi}^{2}}{k^{2}}+\frac{1}{2}\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right].$
(63)

Similar to the calculation of $A^{\prime}np\rightarrow np$, we use soft
radiation approximation for the dark fermion pair absorption and ignore the
Pauli-blocking effect of the nucleons.

When $m_{A^{\prime}}>2m_{\chi}$, we use the narrow width approximation (NWA)
detailed in Appendix VI to approximate the integral. The resulting formula is

$\Gamma^{\rm NWA}_{\chi\bar{\chi}np\rightarrow
np}=\frac{8\epsilon^{2}\alpha_{D}\alpha_{e}m_{A^{\prime}}^{2}n_{n}n_{p}}{3\pi
E|\vec{p}|}\left(\frac{\pi
T}{m_{N}}\right)^{3/2}\left\langle\sigma_{np}^{(2)}(T)\right\rangle\int_{\omega_{-}}^{\omega_{+}}d\omega\,\frac{\mathcal{F}^{L\ast}_{\chi\bar{\chi}np\rightarrow
np}+2\mathcal{F}^{T\ast}_{\chi\bar{\chi}np\rightarrow
np}}{\omega^{3}[1+e^{(\omega-E)/T}]},$ (64)

where the two terms $\mathcal{F}^{L\ast}_{\chi\bar{\chi}np\rightarrow np}$,
$\mathcal{F}^{T\ast}_{\chi\bar{\chi}np\rightarrow np}$ in the integrand are

$\displaystyle\mathcal{F}^{L\ast}_{\chi\bar{\chi}np\rightarrow np}$
$\displaystyle=\frac{\beta_{m|L}^{2}}{\Gamma_{L}}\cdot\frac{m_{A^{\prime}}^{2}}{\omega^{2}}\left[1-\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right],$
(65) $\displaystyle\mathcal{F}^{T\ast}_{\chi\bar{\chi}np\rightarrow np}$
$\displaystyle=\frac{\beta_{m|T}^{2}}{\Gamma_{T}}\cdot\left[\frac{1}{2}+\frac{2m_{\chi}^{2}}{m_{A^{\prime}}^{2}}+\frac{1}{2}\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right],$
(66)

and we now have $|\vec{k}|=\sqrt{\omega^{2}-m_{A^{\prime}}^{2}}$. The dark
photon widths must include the contribution from the dark photon decay via
$A^{\prime}\rightarrow\chi\bar{\chi}$. Thus we have
$\Gamma_{L,T}=(1-e^{-\omega/T})(\Gamma^{L,T}_{A^{\prime},\rm
abs}+\Gamma_{A^{\prime}\rightarrow\chi\bar{\chi}})$. Momentum conservation
demands that $\omega$ is bounded by

$\omega_{\pm}=\frac{m_{A^{\prime}}^{2}}{2m_{\chi}^{2}}\left(E\pm|\vec{p}|\sqrt{1-\frac{4m_{\chi}^{2}}{m_{A^{\prime}}^{2}}}\right).$
(67)

###  VIII.5 $\chi\bar{\chi}\rightarrow e^{-}e^{+}$

The pair-annihilation rate via $\chi(p)+\bar{\chi}(p^{\prime})\rightarrow
e^{-}(p_{3})+e^{+}(p_{4})$ is

$\Gamma_{\chi\bar{\chi}\rightarrow
e^{-}e^{+}}=\frac{\epsilon^{2}\alpha_{D}\alpha_{e}}{4\pi
E}\int_{0}^{\infty}d|\vec{k}|\int_{-1}^{1}d\cos{\theta}\,\int_{E_{3-}}^{E_{3+}}dE_{3}\,\frac{|\vec{k}|}{E^{\prime}}f^{\prime}(1-f_{3})(1-f_{4})\,(\mathcal{F}^{L}_{\chi\bar{\chi}\rightarrow
e^{-}e^{+}}+2\mathcal{F}^{T}_{\chi\bar{\chi}\rightarrow e^{-}e^{+}}),$ (68)

where we have the same definitions for $\cos{\theta}$, $E^{\prime}$, and
$\omega$ as in Eq. (61). Moreover, $f^{\prime}$, $f_{3}$, and $f_{4}$ are the
distribution functions of $\bar{\chi}$, $e^{-}$, and $e^{+}$, respectively.
That is,

$\displaystyle f^{\prime}$ $\displaystyle=\frac{1}{e^{E^{\prime}/T}+1},$ (69)
$\displaystyle f_{3}$ $\displaystyle=\frac{1}{e^{(E_{3}-\mu_{e})/T}+1},$ (70)
$\displaystyle f_{4}$ $\displaystyle=\frac{1}{e^{(E_{4}+\mu_{e})/T}+1},$ (71)

where the positron energy is expressed as $E_{4}=\omega-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 $\mathcal{F}^{L}_{\chi\bar{\chi}\rightarrow
e^{-}e^{+}}$, $\mathcal{F}^{T}_{\chi\bar{\chi}\rightarrow e^{-}e^{+}}$ in the
integrand are

$\displaystyle\mathcal{F}^{L}_{\chi\bar{\chi}\rightarrow e^{-}e^{+}}$
$\displaystyle=\beta_{k|L}^{2}\bar{\beta}_{k|L}^{2}\left[1-\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right]\left[1-\left(\frac{\omega-2E_{3}}{|\vec{k}|}\right)^{2}\right],$
(72) $\displaystyle\mathcal{F}^{T}_{\chi\bar{\chi}\rightarrow e^{-}e^{+}}$
$\displaystyle=\beta_{k|T}^{2}\bar{\beta}_{k|T}^{2}\left[\frac{1}{2}+\frac{2m_{\chi}^{2}}{k^{2}}+\frac{1}{2}\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right]\left[\frac{1}{2}+\frac{2m_{e}^{2}}{k^{2}}+\frac{1}{2}\left(\frac{\omega-2E_{3}}{|\vec{k}|}\right)^{2}\right].$
(73)

Momentum conservation demands that $E_{3}$ is bounded by

$E_{3\pm}=\frac{1}{2}\left[\omega\pm|\vec{k}|\sqrt{1-\frac{4m_{e}^{2}}{k^{2}}}\right].$
(74)

When $m_{A^{\prime}}>2m_{\chi}$, we use NWA (see Appendix VI) to approximate
the annihilation rate. The resulting formula is

$\Gamma^{\rm NWA}_{\chi\bar{\chi}\rightarrow
e^{-}e^{+}}=\frac{\epsilon^{2}\alpha_{e}\alpha_{D}m_{A^{\prime}}^{4}}{8E|\vec{p}|}\int_{\omega_{-}}^{\omega_{+}}d\omega\int_{E_{3-}}^{E_{3+}}dE_{3}\frac{f_{2}(1-f_{3})(1-f_{4})}{\omega|\vec{k}|}(\mathcal{F}^{L\ast}_{\chi\bar{\chi}\rightarrow
e^{-}e^{+}}+2\mathcal{F}^{T\ast}_{\chi\bar{\chi}\rightarrow e^{-}e^{+}}),$
(75)

where we now have $|\vec{k}|=\sqrt{\omega^{2}-m_{A^{\prime}}^{2}}$. The
energies of $\bar{\chi}$ and $e^{+}$ on which their distribution functions
depend are expressed as $E^{\prime}=\omega-E$ and $E_{4}=\omega-E_{3}$. The
two terms $\mathcal{F}^{L\ast}_{\chi\bar{\chi}\rightarrow e^{-}e^{+}}$,
$\mathcal{F}^{T\ast}_{\chi\bar{\chi}\rightarrow e^{-}e^{+}}$ in the integrand
now are

$\displaystyle\mathcal{F}^{L\ast}_{\chi\bar{\chi}\rightarrow e^{-}e^{+}}$
$\displaystyle=\frac{\beta_{m|L}^{2}}{\Gamma_{L}}\left[1-\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right]\left[1-\left(\frac{\omega-2E_{3}}{|\vec{k}|}\right)^{2}\right],$
(76) $\displaystyle\mathcal{F}^{T\ast}_{\chi\bar{\chi}\rightarrow e^{-}e^{+}}$
$\displaystyle=\frac{\beta_{m|T}^{2}}{\Gamma_{T}}\left[\frac{1}{2}+\frac{2m_{\chi}^{2}}{m_{A^{\prime}}^{2}}+\frac{1}{2}\left(\frac{\omega-2E}{|\vec{k}|}\right)^{2}\right]\left[\frac{1}{2}+\frac{2m_{e}^{2}}{m_{A^{\prime}}^{2}}+\frac{1}{2}\left(\frac{\omega-2E_{3}}{|\vec{k}|}\right)^{2}\right],$
(77)

where $\Gamma_{L,T}$ is the same as in Eqs. (65) and (66). The limit of
$E_{3}$ now becomes

$E_{3\pm}=\frac{1}{2}\left[\omega\pm|\vec{k}|\sqrt{1-\frac{4m_{e}^{2}}{m_{A^{\prime}}^{2}}}\right],$
(78)

while the limit of $\omega$ is the same as Eq. (67).

###  VIII.6 $\chi\bar{\chi}\rightarrow\gamma^{\ast}$

Let $\omega_{L,T}=\omega_{L,T}(|\vec{k}|)$ be the energy of the longitudinal
and transverse SM photon as functions of $|\vec{k}|$, and $m_{L,T}$ be the
effective mass satisfying
$m_{L,T}^{2}=\omega_{L,T}^{2}-|\vec{k}|^{2}=\mathrm{Re}\Pi_{L,T}(\omega_{L,T},|\vec{k}|)$.
Then the pair-annihilation rate can be expressed as

$\Gamma_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}=\frac{\epsilon^{2}\alpha_{D}}{4E|\vec{p}|}\sum_{L,T}\int_{0}^{\infty}\frac{|\vec{k}|d|\vec{k}|}{\omega_{L,T}}\frac{g_{L,T}\mathcal{F}^{L,T}_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}\Theta(1-\cos^{2}{\theta_{L,T}})}{[1+e^{(\omega_{L,T}-E)/T}](1-e^{-\omega_{L,T}/T})},$
(79)

where $g_{L}=1$, $g_{T}=2$, and
$\mathcal{F}^{L}_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}$,
$\mathcal{F}^{T}_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}$ are defined as

$\displaystyle\mathcal{F}^{L}_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}$
$\displaystyle=\bar{\beta}^{2}_{m|L}Z_{L}m_{L}^{2}\left[1-\left(\frac{\omega_{L}-2E}{|\vec{k}|}\right)^{2}\right],$
(80) $\displaystyle\mathcal{F}^{T}_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}$
$\displaystyle=\bar{\beta}^{2}_{m|T}Z_{T}m_{T}^{2}\left[\frac{1}{2}+\frac{2m_{\chi}^{2}}{m_{T}^{2}}+\frac{1}{2}\left(\frac{\omega_{T}-2E}{|\vec{k}|}\right)^{2}\right],$
(81)

where $Z_{L,T}$ is the renormalization factor for the longitudinal plasmon and
transverse photon detailed in Appendix VII. The bosonic enhancement for the
plasmon is included in Eq. (79). The step function in Eq. (79) is to ensure
that the momenta are kinematically allowed. The definition of
$\cos{\theta_{L,T}}$ is

$\cos{\theta_{L,T}}=\frac{2\omega_{L,T}E-m_{L,T}^{2}}{2|\vec{k}||\vec{p}|}.$
(82)

We use NWA to approximate the annihilation rate when
$m_{A^{\prime}}>2m_{\chi}$. The resulting expression is

$\Gamma_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}=\frac{\pi\epsilon^{2}\alpha_{D}m_{A^{\prime}}^{6}}{4E|\vec{p}|}\sum_{L,T}\frac{g_{L,T}\mathcal{F}^{L,T\ast}_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}\Theta(\omega_{+}-\omega^{\ast}_{L,T})\Theta(\omega^{\ast}_{L,T}-\omega_{-})}{[1+e^{(\omega^{\ast}_{L}-E)/T}](1-e^{-\omega^{\ast}_{L}/T})},$
(83)

where $\omega^{\ast}_{L,T}$ is the solution of the equation
$\mbox{Re}\Pi_{L,T}(\omega,\sqrt{\omega^{2}-m_{A^{\prime}}^{2}})=m_{A^{\prime}}^{2}$.
$\mathcal{F}^{L,T\ast}_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}$ is now
defined as

$\displaystyle\mathcal{F}^{L\ast}_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}$
$\displaystyle=\frac{1}{\omega^{\ast}_{L}\Gamma_{L}}\left[1-\frac{(\omega^{\ast}_{L}-2E)^{2}}{\omega^{\ast
2}_{L}-m_{A^{\prime}}^{2}}\right]\left|\left.\frac{d\mbox{Re}\Pi_{L}(\omega,\sqrt{\omega^{2}-m_{A^{\prime}}^{2}})}{d\omega}\right|_{\omega^{\ast}_{L}}\right|^{-1},$
(84)
$\displaystyle\mathcal{F}^{T\ast}_{\chi\bar{\chi}\rightarrow\gamma^{\ast}}$
$\displaystyle=\frac{1}{\omega^{\ast}_{T}\Gamma_{T}}\left[\frac{1}{2}+\frac{2m_{\chi}^{2}}{m_{A^{\prime}}^{2}}+\frac{1}{2}\frac{(\omega^{\ast}_{T}-2E)^{2}}{\omega^{\ast
2}_{T}-m_{A^{\prime}}^{2}}\right]\left|\left.\frac{d\mbox{Re}\Pi_{T}(\omega,\sqrt{\omega^{2}-m_{A^{\prime}}^{2}})}{d\omega}\right|_{\omega^{\ast}_{T}}\right|^{-1}.$
(85)

$\omega_{\pm}$ is the same as Eq. (67) and $\Gamma_{L,T}$ is the same as in
Eqs. (65) and (66).

###  VIII.7 $\chi e^{-}\rightarrow\chi e^{-}$ and $\chi p\rightarrow\chi p$

The scattering rate of
$\chi(p)+e^{-}(p_{2})\rightarrow\chi(p^{\prime})+e^{-}(p_{4})$ is given by

$\Gamma_{\chi e^{-}\rightarrow\chi
e^{-}}=\frac{\epsilon^{2}\alpha_{D}\alpha_{e}}{4\pi E}\int
d|\vec{k}|\,d\cos{\theta}\,\int_{E_{2+}}^{\infty}dE_{2}\,\frac{|\vec{k}|}{E^{\prime}}f_{2}(1-f^{\prime})(1-f_{4})\,(\mathcal{F}^{L}_{\chi
e^{-}\rightarrow\chi e^{-}}+2\mathcal{F}^{T}_{\chi e^{-}\rightarrow\chi
e^{-}}),$ (86)

where $k=p-p^{\prime}$ is the momentum transfer, $\cos{\theta}$ and
$E^{\prime}$ are the same in Eq. (61) while $\omega=E-E^{\prime}$ instead.
$f_{2}$, $f^{\prime}$, and $f_{4}$ are the distribution functions of the
initial state $e^{-}$, final state $\chi$, and final state $e^{-}$,
respectively. The final state electron energy is expressed as
$E_{4}=\omega+E_{2}$ by energy conservation. The two terms
$\mathcal{F}^{L}_{\chi e^{-}\rightarrow\chi e^{-}}$, $\mathcal{F}^{T}_{\chi
e^{-}\rightarrow\chi e^{-}}$ are the same as Eqs. (72) and (73) except that
$-E_{3}$ is replaced with $+E_{2}$, and $E_{2}$ is bounded below by

$E_{2+}=\frac{1}{2}\left(|\vec{k}|\sqrt{1-\frac{4m_{e}^{2}}{k^{2}}}-\omega\right).$
(87)

Note that the momentum transfer is spacelike ($k^{2}<0$) for the
$\chi$-$e^{-}$ scattering. In this scenario, we take the dark photon width to
be zero, while $\mbox{Im }\Pi_{L,T}$ is determined by the imaginary parts of
Eqs. (43) and (44) with $\omega\rightarrow\omega+i\varepsilon$ for an
infinitesimal $\varepsilon>0$ Kuznetsov and Mikheev (2013).

The scattering rate via $\chi p\rightarrow\chi p$ is the same as $\chi
e^{-}\rightarrow\chi e^{-}$ if the differences in masses and distribution
functions of the proton and the electron are taken into account. However, we
can approximate the $\chi$-$p$ scattering as an elastic scattering due to the
fact that $m_{p}$ is much larger than $m_{\chi}$. We also ignore the Pauli-
blocking effect of proton. Thus the scattering rate of $\chi$-$p$ scattering
can be approximated by

$\Gamma_{\chi p\rightarrow\chi
p}=\frac{4\pi\epsilon^{2}\alpha_{e}\alpha_{D}n_{p}|\vec{p}|}{E(1+e^{-E/T})}\int_{-1}^{1}d\cos{\theta}\beta_{k|L}^{2}\cdot\frac{2E^{2}-|\vec{p}|^{2}(1-\cos{\theta})}{(k^{2}-m_{A^{\prime}}^{2})^{2}},$
(88)

where the momentum transfer is
$k^{2}\simeq-|\vec{k}|^{2}\simeq-2|\vec{p}|^{2}(1-\cos{\theta})$. Note that in
this static limit ($\omega\rightarrow 0$), the contribution from the
transverse mode, corresponding to the classical magnetic field, is suppressed.
And the plasma factor $\beta_{k|L}^{2}$ of the longitudinal mode,
corresponding to the static electric field, accounts for the screening effect
in medium Kapusta and Gale (2006).

###  VIII.8 $A^{\prime}\chi\rightarrow\chi A^{\prime}$ and
$A^{\prime}A^{\prime}\leftrightarrow\chi\bar{\chi}$

The spin-averaged matrix element squared of
$A^{\prime}(p_{1})+\chi(p_{2})\rightarrow\chi(p_{3})+A^{\prime}(p_{4})$ is

$\begin{split}\overline{|\mathcal{M}|^{2}}_{A^{\prime}\chi\rightarrow\chi
A^{\prime}}&=\frac{64\pi^{2}\alpha_{D}^{2}}{3}\left\\{(m_{A^{\prime}}^{2}+2m_{\chi}^{2})^{2}\left[\frac{1}{(s-m_{\chi}^{2})^{2}}+\frac{1}{(t-m_{\chi}^{2})^{2}}\right]+\frac{8(m_{\chi}^{4}-m_{A^{\prime}}^{4})}{(s-m_{\chi}^{2})(t-m_{\chi}^{2})}\right.\\\
&\left.+\,4(m_{\chi}^{2}+m_{A^{\prime}}^{2})\left(\frac{1}{s-m_{\chi}^{2}}+\frac{1}{t-m_{\chi}^{2}}\right)-\left(\frac{s-m_{\chi}^{2}}{t-m_{\chi}^{2}}+\frac{t-m_{\chi}^{2}}{s-m_{\chi}^{2}}\right)\right\\},\end{split}$
(89)

where $s=(p_{1}+p_{2})^{2}$ and $t=(p_{1}-p_{3})^{2}$ are the Mandelstam
variables. By crossing symmetry, the matrix element of
$A^{\prime}(p_{1})+A^{\prime}(p_{2})\rightarrow\chi(p_{3})+\bar{\chi}(p_{4})$
is the same as Eq. (89) except that $s$ is replaced by $u=(p_{1}-p_{4})^{2}$
and the matrix element is multiplied by a factor of $-2/3$. (The minus sign is
due to the crossing of one fermion state.) For
$\chi(p_{1})+\bar{\chi}(p_{2})\rightarrow
A^{\prime}(p_{3})+A^{\prime}(p_{4})$, we can again reuse Eq. (89), replace $s$
by $u$, and multiply by a factor of $-3/2$.

###  VIII.9 $\chi\bar{\chi}\rightarrow\chi\bar{\chi}$ and
$\chi\chi\rightarrow\chi\chi$

The spin-averaged matrix element squared of
$\chi(p_{1})+\bar{\chi}(p_{2})\rightarrow\chi(p_{3})+\bar{\chi}(p_{4})$ is

$\begin{split}&\overline{|\mathcal{M}|^{2}}_{\chi\bar{\chi}\rightarrow\chi\bar{\chi}}=64\pi^{2}\alpha_{D}^{2}\left\\{\left(\frac{5}{2}m_{A^{\prime}}^{4}-4m_{\chi}^{2}m_{A^{\prime}}^{2}+4m_{\chi}^{4}\right)\left[\frac{1}{(t-m_{A^{\prime}}^{2})^{2}}+\frac{1}{(s-m_{A^{\prime}}^{2})^{2}}\right]+\frac{4(m_{A^{\prime}}^{4}-m_{\chi}^{4})}{(t-m_{A^{\prime}}^{2})(s-m_{A^{\prime}}^{2})}\right.\\\
&\left.+(3m_{A^{\prime}}^{2}-4m_{\chi}^{2})\left[\frac{s-m_{A^{\prime}}^{2}}{(t-m_{A^{\prime}}^{2})^{2}}+\frac{t-m_{A^{\prime}}^{2}}{(s-m_{A^{\prime}}^{2})^{2}}\right]+6m_{A^{\prime}}^{2}\left(\frac{1}{t-m_{A^{\prime}}^{2}}+\frac{1}{s-m_{A^{\prime}}^{2}}\right)+\left(\frac{s-m_{A^{\prime}}^{2}}{t-m_{A^{\prime}}^{2}}+1+\frac{t-m_{A^{\prime}}^{2}}{s-m_{A^{\prime}}^{2}}\right)^{2}\right\\}.\end{split}$
(90)

Due to crossing symmetry, we can reuse Eq. (LABEL:eq:xx'xx') for the
scattering process $\chi(p_{1})+\chi(p_{2})\rightarrow\chi(p_{3})+\chi(p_{4})$
with $s$ replaced by $u$.

###  VIII.10 $A^{\prime}\leftrightarrow\chi\bar{\chi}$

For the decay process $A^{\prime}\rightarrow\chi\bar{\chi}$, we take into
account the Pauli-blocking effect of the dark fermions. The decay rate is
given by

$\Gamma_{A^{\prime}\rightarrow\chi\bar{\chi}}=\frac{\alpha_{D}(m_{A^{\prime}}^{2}+2m_{\chi}^{2})}{3\omega|\vec{k}|}\frac{2T}{1-e^{-\omega/T}}\ln{\left\\{\frac{\cosh{\left[(1+\xi_{0})\omega/4T\right]}}{\cosh{\left[(1-\xi_{0})\omega/4T\right]}}\right\\}},$
(91)

where $\xi_{0}$ is similar to Eq. (60),

$\xi_{0}=\frac{|\vec{k}|}{\omega}\sqrt{1-\frac{4m_{\chi}^{2}}{m_{A^{\prime}}^{2}}}.$
(92)

Note that in the limit $T\rightarrow 0$, Eq. (91) reduces to the decay rate in
vacuum given by Eq. (2.6) in Chang  _et al._ (2018).

The inverse decay rate of $\chi\bar{\chi}\rightarrow A^{\prime}$, taking into
account the bosonic enhancement of the dark photon, is given by

$\Gamma_{\chi\bar{\chi}\rightarrow
A^{\prime}}=\frac{\alpha_{D}(m_{A^{\prime}}^{2}+2m_{\chi}^{2})}{2E|\vec{p}|}\frac{T}{1+e^{-E/T}}\ln{\left\\{\frac{\sinh{(\omega_{+}/2T)}\cosh{[(\omega_{-}-E)/2T]}}{\sinh{(\omega_{-}/2T)}\cosh{[(\omega_{+}-E)/2T]}}\right\\}},$
(93)

where $\omega_{\pm}$ is the same as Eq. (67).

##  IX Detailed Work Flow for Computing DS Luminosity

  1. 1.

We take the following steps to determine if a DS particle species is in the
diffuse regime or not:

    1. (a)

Compute $\tilde{\lambda}^{-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)$.

    2. (b)

Compute $\Delta t_{i}$, the escape timescale for each particle species $i$
with Eqs. (21) and (22).

    3. (c)

Compute $N_{i}$, the estimated abundance of each particle species with Eq.
(20), and obtain the relative abundance $\eta_{i}\equiv N_{i}/N_{i}^{eq}$,
where $N_{i}^{\rm eq}$ is the abundance of DS particle species $i$ assuming
they are in thermal equilibrium with the SM medium. If $\eta_{i}>1$, we simply
set $\eta_{i}=1$.

    4. (d)

Compute $\lambda^{-1}_{i}$, the IMFP of each species rescaled by the relative
abundance as

$\lambda^{-1}_{i}=\sum_{\rm
possible\,final\,states}(\tilde{\lambda}^{-1}_{i\rightarrow\rm
final\,states}+\sum_{j}\eta_{j}\tilde{\lambda}^{-1}_{i+j\rightarrow\rm
final\,states}),$ (94)

where $j$ runs over all possible DS particle species.

    5. (e)

Check the diffuse criteria for the dark photons: if
$\eta_{A^{\prime}_{L,T}}=1$ and
$\langle\lambda^{-1}_{A^{\prime}_{L,T}}(R_{\nu})\rangle-\langle\tilde{\lambda}^{-1}_{A^{\prime}_{L,T}\rightarrow
e^{-}e^{+}}(R_{\nu})\rangle>R_{\nu}^{-1}$, then $A^{\prime}_{L,T}$ is in the
diffuse limit. Otherwise, $A^{\prime}_{L,T}$ is treated as nondiffuse
particles.

    6. (f)

Check the diffuse criteria for the dark fermion: if $\eta_{\chi}=1$ and
$\langle\lambda^{-1}_{\chi}(R_{\nu})\rangle>R_{\nu}^{-1}$, then $\chi$ is in
the diffuse regime. Otherwise, $\chi$ is treated as nondiffuse particles.

  2. 2.

Depending on whether the DS particles are in the diffuse limit or not, we
compute the luminosity of each particle species as follows:

    1. (a)

Regardless of the DS masses, the dark photon luminosity is given by

$L_{A^{\prime}_{L,T}}=\begin{cases}L_{A^{\prime}_{L,T},\rm diff.},\,\mbox{for
diffuse }A^{\prime}_{L,T}\mbox{ [use Eq.~{}\eqref{eq:flux}]},\\\
L_{A^{\prime}_{L,T},\rm nondiff.},\,\mbox{for nondiffuse
}A^{\prime}_{L,T}\mbox{ [use Eq.~{}\eqref{eq:LA'emis}]}.\end{cases}$ (95)

    2. (b)

When $m_{A^{\prime}}<2m_{\chi}$, the dark fermion luminosity is given by

$L_{\chi}=\begin{cases}L_{\chi,\rm diff.},\,\mbox{for diffuse }\chi\mbox{ [use
Eq.~{}\eqref{eq:flux}]},\\\ L_{\chi,\rm nondiff.},\,\mbox{for nondiffuse
}\chi\mbox{ [use Eq.~{}\eqref{eq:Lchiemis}]}.\end{cases}$ (96)

    3. (c)

When $m_{A^{\prime}}>2m_{\chi}$, the dark fermion luminosity is given by

$L_{\chi}=\begin{cases}L_{\chi,\rm diff.},\,\mbox{if }\chi\mbox{ is in the
diffuse limit},\\\ L_{\chi,\rm nondiff.},\,\mbox{if only
}A^{\prime}_{L},A^{\prime}_{T}\mbox{ are in the diffuse limit},\\\
\frac{2}{3}L_{\chi,\rm nondiff.},\,\mbox{if only }A^{\prime}_{T}\mbox{ is in
the diffuse limit},\\\ \frac{1}{3}L_{\chi,\rm nondiff.},\,\mbox{if only
}A^{\prime}_{L}\mbox{ is in the diffuse limit},\\\ 0,\,\mbox{if no DS particle
species is in the diffuse limit}.\end{cases}$ (97)

    4. (d)

The total DS luminosity is $L_{\rm
X}=L_{A^{\prime}_{L}}+L_{A^{\prime}_{T}}+2L_{\chi}$.

##  X 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 Prialnik (2000). The equation of
radiative transfer is given by

$\frac{1}{\rho}\frac{\partial I}{\partial r}\cos{\theta}+\kappa I-j=0,$ (98)

where $\rho$ is the matter density, $I$ is the intensity of radiation per unit
solid angle per unit frequency, $\theta$ is the angle between the direction of
radiation and the radial direction, $\kappa$ is the opacity, and $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 $\kappa=\kappa_{s}+\kappa_{a}$, and between the radiation by
scattering and radiation emitted by matter $j=j_{s}+j_{em}$. In equilibrium
and isotropic environment,

$\displaystyle I^{\rm iso,eq}$ $\displaystyle\equiv
B=\frac{g}{(2\pi)^{3}}\frac{E^{2}p}{e^{E/T}\pm 1},$ (99) $\displaystyle
j_{em}^{\rm iso,eq}$ $\displaystyle=\kappa_{a}B,$ (100)

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}=\kappa_{a}\left(1\pm e^{-E/T}\right)B\mp\kappa_{a}e^{-E/T}I.$ (101)

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. (98),

$\frac{1}{\rho}\frac{\partial I}{\partial
r}\cos{\theta}+\kappa_{s}I-j_{s}+\kappa_{a}^{\ast}(I-B)=0,$ (102)

where $\kappa_{a}^{\ast}\equiv\kappa_{a}(1\pm e^{-E/T})$. We assume that the
radiation intensity $I$ is very close to $I^{\rm iso,eq}$. Therefore, we can
expand $I$ in terms of Legendre polynomials $P_{n}(\cos{\theta})$ and
substitute it back into Eq. (102). Keeping the terms up to $n=1$, we obtain

$I\simeq B-\frac{1}{\rho(\kappa_{a}^{\ast}+\kappa_{s})}\frac{\partial
B}{\partial r}\cos{\theta}.$ (103)

It follows that the energy flux through the spherical surface of radius $r$ is

$\begin{split}L(r)&=4\pi r^{2}\int I\cos{\theta}dEd\Omega\\\
&=-\frac{16}{3}\pi^{2}r^{2}\frac{dT}{dr}\int\frac{1}{\rho(\kappa_{a}^{\ast}+\kappa_{s})}\frac{\partial
B}{\partial T}dE\\\
&=-\frac{2gr^{2}}{3\pi}\frac{T^{3}dT}{dr}\int_{m/T}^{\infty}\frac{\xi^{3}\sqrt{\xi^{2}-\left(\frac{m}{T}\right)^{2}}}{\lambda^{-1}(E=\xi
T,r)}\frac{e^{\xi}}{\left(e^{\xi}\pm 1\right)^{2}}d\xi,\end{split}$ (104)

where $m$ is the mass of the particle carrying the radiation, and
$\lambda^{-1}$ is the effective IMFP defined as

$\lambda^{-1}=\rho(\kappa_{a}^{\ast}+\kappa_{s}).$ (105)

We can identify the absorptive IMFP $\lambda_{\rm abs}^{-1}$ with
$\rho\kappa_{a}$ and the scattering IMFP $\lambda_{\rm sca}^{-1}$ with
$\rho\kappa_{s}$. Thus, we can express the effective IMFP as

$\lambda^{-1}(E,r)\equiv\lambda_{\rm abs}^{-1}(E,r)(1\pm
e^{-E/T})+\lambda_{\rm sca}^{-1}(E,r).$ (106)

We note that since the integral in Eq. (104) is dominated by the energies with
small $\lambda^{-1}(E,r)$, it is possible that Eq. (104) could overestimate
the energy flux if $\lambda^{-1}(E,r)\lesssim R_{\nu}$ for some energies $E$.
We found that it occurs near the switching of diffuse and nondiffuse regimes,
and thus it leads to orders of magnitude jumps of energy luminosity. To avoid
this caveat, we approximate the IMFP $\lambda^{-1}(E,r)$ by the thermally
averaged IMFP $\langle\lambda^{-1}(r)\rangle$ defined as

$\langle\lambda^{-1}(r)\rangle=\frac{\int
d^{3}p\,f(E,T(r))\lambda^{-1}(E,r)}{\int d^{3}p\,f(E,T(r))}.$ (107)

Thus, the energy flux is approximately given by

$L(r)\simeq-\frac{2gr^{2}}{3\pi}\frac{T^{3}dT}{dr}\frac{1}{\langle\lambda^{-1}(r)\rangle}\int_{m/T}^{\infty}\xi^{3}\sqrt{\xi^{2}-\left(\frac{m}{T}\right)^{2}}\frac{e^{\xi}}{\left(e^{\xi}\pm
1\right)^{2}}d\xi.$ (108)

## References

  * Hirata  _et al._ (1987) K. Hirata _et al._ (Kamiokande-II Collaboration), Phys. Rev. Lett. 58, 1490 (1987).
  * Bionta  _et al._ (1987) R. M. Bionta _et al._ , Phys. Rev. Lett. 58, 1494 (1987).
  * Alekseev  _et al._ (1988) E. N. Alekseev, L. N. Alekseeva, I. V. Krivosheina, and V. I. Volchenko, Phys. Lett. B205, 209 (1988).
  * Sato and Suzuki (1987) K. Sato and H. Suzuki, Phys. Rev. Lett. 58, 2722 (1987).
  * Spergel  _et al._ (1987) D. N. Spergel, T. Piran, A. Loeb, J. Goodman, and J. N. Bahcall, Science 237, 1471 (1987).
  * Bahcall  _et al._ (1987) J. N. Bahcall, T. Piran, W. H. Press, and D. N. Spergel, Nature (London) 327, 682 (1987).
  * Burrows and Lattimer (1987) A. Burrows and J. M. Lattimer, Astrophys. J. 318, L63 (1987).
  * Schramm and Truran (1990) D. Schramm and J. Truran, Phys. Rep. 189, 89 (1990).
  * Loredo and Lamb (2002) T. J. Loredo and D. Q. Lamb, Phys. Rev. D 65, 063002 (2002), arXiv:astro-ph/0107260 [astro-ph] .
  * Frieman  _et al._ (1988) J. A. Frieman, H. E. Haber, and K. Freese, Phys. Lett. B 200, 115 (1988).
  * Kolb and Turner (1989) E. W. Kolb and M. S. Turner, Phys. Rev. Lett. 62, 509 (1989).
  * Berezhiani and Smirnov (1989) Z. Berezhiani and A. Smirnov, Phys. Lett. B 220, 279 (1989).
  * Farzan (2003) Y. Farzan, Phys. Rev. D 67, 073015 (2003), arXiv:hep-ph/0211375 [hep-ph] .
  * Raffelt and Seckel (1988) G. Raffelt and D. Seckel, Phys. Rev. Lett. 60, 1793 (1988).
  * Turner (1988) M. S. Turner, Phys. Rev. Lett. 60, 1797 (1988).
  * Mayle  _et al._ (1988) R. Mayle, J. R. Wilson, J. R. Ellis, K. A. Olive, D. N. Schramm, and G. Steigman, Phys. Lett. B 203, 188 (1988).
  * Brinkmann and Turner (1988) R. P. Brinkmann and M. S. Turner, Phys. Rev. D 38, 2338 (1988).
  * Janka  _et al._ (1996) H.-T. Janka, W. Keil, G. Raffelt, and D. Seckel, Phys. Rev. Lett. 76, 2621 (1996), arXiv:astro-ph/9507023 [astro-ph] .
  * Keil  _et al._ (1997) W. Keil, H.-T. Janka, D. N. Schramm, G. Sigl, M. S. Turner, and J. R. Ellis, Phys. Rev. D 56, 2419 (1997), arXiv:astro-ph/9612222 [astro-ph] .
  * Fischer  _et al._ (2016) T. Fischer, S. Chakraborty, M. Giannotti, A. Mirizzi, A. Payez, and A. Ringwald, Phys. Rev. D 94, 085012 (2016), arXiv:1605.08780 [astro-ph.HE] .
  * Carenza  _et al._ (2019) P. Carenza, T. Fischer, M. Giannotti, G. Guo, G. Martínez-Pinedo, and A. Mirizzi, J. Cosmol. Astropart. Phys. 10, 016 (2019), [Erratum: JCAP 05, E01 (2020)], arXiv:1906.11844 [hep-ph] .
  * Bollig  _et al._ (2020) R. Bollig, W. DeRocco, P. W. Graham, and H.-T. Janka, Phys. Rev. Lett. 125, 051104 (2020), arXiv:2005.07141 [hep-ph] .
  * Lucente  _et al._ (2020) G. Lucente, P. Carenza, T. Fischer, M. Giannotti, and A. Mirizzi, J. Cosmol. Astropart. Phys. 12, 008 (2020), arXiv:2008.04918 [hep-ph] .
  * Raffelt and Zhou (2011) G. G. Raffelt and S. Zhou, Phys. Rev. D 83, 093014 (2011), arXiv:1102.5124 [hep-ph] .
  * Argüelles  _et al._ (2019) C. A. Argüelles, V. Brdar, and J. Kopp, Phys. Rev. D 99, 043012 (2019), arXiv:1605.00654 [hep-ph] .
  * Suliga  _et al._ (2019) A. M. Suliga, I. Tamborra, and M.-R. Wu, J. Cosmol. Astropart. Phys. 12, 019 (2019), arXiv:1908.11382 [astro-ph.HE] .
  * Syvolap  _et al._ (2019) V. Syvolap, O. Ruchayskiy, and A. Boyarsky, (2019), arXiv:1909.06320 [hep-ph] .
  * Mastrototaro  _et al._ (2020) L. Mastrototaro, A. Mirizzi, P. D. Serpico, and A. Esmaili, J. Cosmol. Astropart. Phys. 01, 010 (2020), arXiv:1910.10249 [hep-ph] .
  * Suliga  _et al._ (2020) A. M. Suliga, I. Tamborra, and M.-R. Wu, J. Cosmol. Astropart. Phys. 08, 018 (2020), arXiv:2004.11389 [astro-ph.HE] .
  * Dent  _et al._ (2012) J. B. Dent, F. Ferrer, and L. M. Krauss, (2012), arXiv:1201.2683 [astro-ph.CO] .
  * Rrapaj and Reddy (2016) E. Rrapaj and S. Reddy, Phys. Rev. C 94, 045805 (2016), arXiv:1511.09136 [nucl-th] .
  * Mahoney  _et al._ (2017) C. Mahoney, A. K. Leibovich, and A. R. Zentner, Phys. Rev. D 96, 043018 (2017).
  * Chang  _et al._ (2017) J. H. Chang, R. Essig, and S. D. McDermott, J. High Energy Phys. 01, 107 (2017), arXiv:1611.03864 [hep-ph] .
  * Hardy and Lasenby (2017) E. Hardy and R. Lasenby, J. High Energy Phys. 02, 033 (2017), arXiv:1611.05852 [hep-ph] .
  * Chang  _et al._ (2018) J. H. Chang, R. Essig, and S. D. McDermott, J. High Energy Phys. 09, 051 (2018), arXiv:1803.00993 [hep-ph] .
  * Guha  _et al._ (2017) A. Guha, J. Selvaganapathy, and P. K. Das, Phys. Rev. D 95, 015001 (2017), arXiv:1509.05901 [hep-ph] .
  * Guha  _et al._ (2019) A. Guha, P. S. B. Dev, and P. K. Das, J. Cosmol. Astropart. Phys. 02, 032 (2019), arXiv:1810.00399 [hep-ph] .
  * Ishizuka and Yoshimura (1990) N. Ishizuka and M. Yoshimura, Prog. Theor. Phys. 84, 233 (1990).
  * Arndt and Fox (2003) D. Arndt and P. J. Fox, J. High Energy Phys. 02, 036 (2003), arXiv:hep-ph/0207098 [hep-ph] .
  * Keung  _et al._ (2014) W.-Y. Keung, K.-W. Ng, H. Tu, and T.-C. Yuan, Phys. Rev. D 90, 075014 (2014), arXiv:1312.3488 [hep-ph] .
  * Tu and Ng (2017) H. Tu and K.-W. Ng, J. High Energy Phys. 07, 108 (2017), arXiv:1706.08340 [hep-ph] .
  * Dev  _et al._ (2020) P. B. Dev, R. N. Mohapatra, and Y. Zhang, J. Cosmol. Astropart. Phys. 08, 003 (2020), [Erratum: JCAP 11, E01 (2020)], arXiv:2005.00490 [hep-ph] .
  * Croon  _et al._ (2021) D. Croon, G. Elor, R. K. Leane, and S. D. McDermott, J. High Energy Phys. 01, 107 (2021), arXiv:2006.13942 [hep-ph] .
  * Camalich  _et al._ (2020) J. M. Camalich, J. Terol-Calvo, L. Tolos, and R. Ziegler, (2020), arXiv:2012.11632 [hep-ph] .
  * Hanhart  _et al._ (2001a) C. Hanhart, D. R. Phillips, S. Reddy, and M. J. Savage, Nucl. Phys. B595, 335 (2001a), arXiv:nucl-th/0007016 [nucl-th] .
  * Hanhart  _et al._ (2001b) C. Hanhart, J. A. Pons, D. R. Phillips, and S. Reddy, Phys. Lett. B509, 1 (2001b), arXiv:astro-ph/0102063 [astro-ph] .
  * Hannestad  _et al._ (2007) S. Hannestad, G. Raffelt, and Y. Y. Y. Wong, Phys. Rev. D 76, 121701 (2007), arXiv:0708.1404 [hep-ph] .
  * Freitas and Wyler (2007) A. Freitas and D. Wyler, J. High Energy Phys. 12, 033 (2007), arXiv:0708.4339 [hep-ph] .
  * Bar  _et al._ (2020) N. Bar, K. Blum, and G. D’Amico, Phys. Rev. D 101, 123025 (2020), arXiv:1907.05020 [hep-ph] .
  * Sung  _et al._ (2019) A. Sung, H. Tu, and M.-R. Wu, Phys. Rev. D 99, 121305 (2019), arXiv:1903.07923 [hep-ph] .
  * Kazanas  _et al._ (2014) D. Kazanas, R. N. Mohapatra, S. Nussinov, V. L. Teplitz, and Y. Zhang, Nucl. Phys. B890, 17 (2014), arXiv:1410.0221 [hep-ph] .
  * DeRocco  _et al._ (2019a) W. DeRocco, P. W. Graham, D. Kasen, G. Marques-Tavares, and S. Rajendran, J. High Energy Phys. 02, 171 (2019a), arXiv:1901.08596 [hep-ph] .
  * DeRocco  _et al._ (2019b) W. DeRocco, P. W. Graham, D. Kasen, G. Marques-Tavares, and S. Rajendran, Phys. Rev. D 100, 075018 (2019b), arXiv:1905.09284 [hep-ph] .
  * Darmé  _et al._ (2020) L. Darmé, F. Giacchino, E. Nardi, and M. Raggi, (2020), arXiv:2012.07894 [hep-ph] .
  * Spergel and Steinhardt (2000) D. N. Spergel and P. J. Steinhardt, Phys. Rev. Lett. 84, 3760 (2000), arXiv:astro-ph/9909386 .
  * Vogelsberger  _et al._ (2012) M. Vogelsberger, J. Zavala, and A. Loeb, Mon. Not. R. Astron. Soc. 423, 3740 (2012), arXiv:1201.5892 [astro-ph.CO] .
  * Rocha  _et al._ (2013) M. Rocha, A. H. G. Peter, J. S. Bullock, M. Kaplinghat, S. Garrison-Kimmel, J. Onorbe, and L. A. Moustakas, Mon. Not. R. Astron. Soc. 430, 81 (2013), arXiv:1208.3025 [astro-ph.CO] .
  * Tulin and Yu (2018) S. Tulin and H.-B. Yu, Phys. Rep. 730, 1 (2018), arXiv:1705.02358 [hep-ph] .
  * Tulin  _et al._ (2013) S. Tulin, H.-B. Yu, and K. M. Zurek, Phys. Rev. D 87, 115007 (2013), arXiv:1302.3898 [hep-ph] .
  * Bernal  _et al._ (2016) N. Bernal, X. Chu, C. Garcia-Cely, T. Hambye, and B. Zaldivar, J. Cosmol. Astropart. Phys. 03, 018 (2016), arXiv:1510.08063 [hep-ph] .
  * Kawasaki  _et al._ (2015) M. Kawasaki, K. Kohri, T. Moroi, and Y. Takaesu, Phys. Lett. B 751, 246 (2015), arXiv:1509.03665 [hep-ph] .
  * Bringmann  _et al._ (2017) T. Bringmann, F. Kahlhoefer, K. Schmidt-Hoberg, and P. Walia, Phys. Rev. Lett. 118, 141802 (2017), arXiv:1612.00845 [hep-ph] .
  * Elor  _et al._ (2016) G. Elor, N. L. Rodd, T. R. Slatyer, and W. Xue, J. Cosmol. Astropart. Phys. 06, 024 (2016), arXiv:1511.08787 [hep-ph] .
  * Ackermann  _et al._ (2015) M. Ackermann _et al._ (Fermi-LAT Collaboration), J. Cosmol. Astropart. Phys. 09, 008 (2015), arXiv:1501.05464 [astro-ph.CO] .
  * Abdallah  _et al._ (2016) H. Abdallah _et al._ (H.E.S.S. Collaboration), Phys. Rev. Lett. 117, 111301 (2016), arXiv:1607.08142 [astro-ph.HE] .
  * Feng  _et al._ (2016) J. L. Feng, J. Smolinsky, and P. Tanedo, Phys. Rev. D 93, 115036 (2016), [Erratum: Phys. Rev.D96,no.9,099903(2017)], arXiv:1602.01465 [hep-ph] .
  * Leane  _et al._ (2017) R. K. Leane, K. C. Y. Ng, and J. F. Beacom, Phys. Rev. D 95, 123016 (2017), arXiv:1703.04629 [astro-ph.HE] .
  * Kouvaris  _et al._ (2016) C. Kouvaris, K. Langæble, and N. G. Nielsen, J. Cosmol. Astropart. Phys. 10, 012 (2016), arXiv:1607.00374 [hep-ph] .
  * Essig  _et al._ (2019) R. Essig, S. D. McDermott, H.-B. Yu, and Y.-M. Zhong, Phys. Rev. Lett. 123, 121102 (2019), arXiv:1809.01144 [hep-ph] .
  * Chang  _et al._ (2019) J. H. Chang, D. Egana-Ugrinovic, R. Essig, and C. Kouvaris, J. Cosmol. Astropart. Phys. 03, 036 (2019), arXiv:1812.07000 [hep-ph] .
  * Depta  _et al._ (2019) P. F. Depta, M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, J. Cosmol. Astropart. Phys. 04, 029 (2019), arXiv:1901.06944 [hep-ph] .
  * Bernal  _et al._ (2020) N. Bernal, X. Chu, S. Kulkarni, and J. Pradler, Phys. Rev. D 101, 055044 (2020), arXiv:1912.06681 [hep-ph] .
  * Foot and Vagnozzi (2015) R. Foot and S. Vagnozzi, Phys. Rev. D 91, 023512 (2015), arXiv:1409.7174 [hep-ph] .
  * Mezzacappa  _et al._ (2020) A. Mezzacappa, E. Endeve, O. B. Messer, and S. W. Bruenn, (2020), arXiv:2010.09013 [astro-ph.HE] .
  * Holdom (1986) B. Holdom, Phys. Lett. B 166, 196 (1986).
  * Okun (1982) L. Okun, Sov. Phys. JETP 56, 502 (1982).
  * Pospelov  _et al._ (2008) M. Pospelov, A. Ritz, and M. B. Voloshin, Phys. Lett. B 662, 53 (2008), arXiv:0711.4866 [hep-ph] .
  * Fischer  _et al._ (2010) T. Fischer, S. C. Whitehouse, A. Mezzacappa, F. K. Thielemann, and M. Liebendorfer, Astron. Astrophys. 517, A80 (2010), arXiv:0908.1871 [astro-ph.HE] .
  * Raffelt (1990) G. G. Raffelt, Phys. Rep. 198, 1 (1990).
  * An  _et al._ (2013) H. An, M. Pospelov, and J. Pradler, Phys. Lett. B 725, 190 (2013), arXiv:1302.3884 [hep-ph] .
  * Keil  _et al._ (2003) M. T. Keil, G. G. Raffelt, and H.-T. Janka, Astrophys. J. 590, 971 (2003), arXiv:astro-ph/0208035 .
  * Shapiro and Teukolsky (1983) S. L. Shapiro and S. A. Teukolsky,  _Black Holes, White Dwarfs, and Neutron Stars : The Physics of Compact Objects_ (Wiley, New York, 1983).
  * Andreas  _et al._ (2012) S. Andreas, C. Niebuhr, and A. Ringwald, Phys. Rev. D 86, 095019 (2012), arXiv:1209.6083 [hep-ph] .
  * Jaeckel  _et al._ (2008) J. Jaeckel, J. Redondo, and A. Ringwald, Phys. Rev. Lett. 101, 131801 (2008), arXiv:0804.4157 [astro-ph] .
  * Mirizzi  _et al._ (2009) A. Mirizzi, J. Redondo, and G. Sigl, J. Cosmol. Astropart. Phys. 03, 026 (2009), arXiv:0901.0014 [hep-ph] .
  * Redondo and Raffelt (2013) J. Redondo and G. Raffelt, J. Cosmol. Astropart. Phys. 08, 034 (2013), arXiv:1305.2920 [hep-ph] .
  * Fradette  _et al._ (2014) A. Fradette, M. Pospelov, J. Pradler, and A. Ritz, Phys. Rev. D 90, 035022 (2014), arXiv:1407.0993 [hep-ph] .
  * An  _et al._ (2020) H. An, M. Pospelov, J. Pradler, and A. Ritz, Phys. Rev. D 102, 115022 (2020), arXiv:2006.13929 [hep-ph] .
  * Li  _et al._ (2020) J.-T. Li, G. M. Fuller, and E. Grohs, J. Cosmol. Astropart. Phys. 12, 049 (2020), arXiv:2009.14325 [astro-ph.CO] .
  * Sieverding  _et al._ (2021) A. Sieverding, E. Rrapaj, G. Guo, and Y. Z. Qian, (2021), arXiv:2101.08672 [astro-ph.SR] .
  * Kreisch  _et al._ (2020) C. D. Kreisch, F.-Y. Cyr-Racine, and O. Doré, Phys. Rev. D 101, 123505 (2020), arXiv:1902.00534 [astro-ph.CO] .
  * Ellis  _et al._ (2018) J. Ellis, A. Hektor, G. Hütsi, K. Kannike, L. Marzola, M. Raidal, and V. Vaskonen, Phys. Lett. B 781, 607 (2018), arXiv:1710.05540 [astro-ph.CO] .
  * Nelson  _et al._ (2019) A. Nelson, S. Reddy, and D. Zhou, J. Cosmol. Astropart. Phys. 07, 012 (2019), arXiv:1803.03266 [hep-ph] .
  * Bauswein  _et al._ (2020) A. Bauswein, G. Guo, J.-H. Lien, Y.-H. Lin, and M.-R. Wu, (2020), arXiv:2012.11908 [astro-ph.HE] .
  * Weldon (1983) H. Weldon, Phys. Rev. D 28, 2007 (1983).
  * Braaten and Segel (1993) E. Braaten and D. Segel, Phys. Rev. D 48, 1478 (1993), arXiv:hep-ph/9302213 [hep-ph] .
  * Lin (2019) T. Lin, Proc. Sci. 333, 009 (2019), arXiv:1904.07915 [hep-ph] .
  * Kuznetsov and Mikheev (2013) A. Kuznetsov and N. Mikheev,  _Electroweak Processes in External Active Media_, Springer Tracts in Modern Physics, Vol. 252 (Springer, New York, 2013).
  * Kapusta and Gale (2006) J. I. Kapusta and C. Gale,  _Finite-Temperature Field Theory: Principles and Applications_, 2nd ed., Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, United Kingdom, 2006).
  * Prialnik (2000) D. Prialnik,  _An Introduction to the Theory of Stellar Structure and Evolution_ (Cambridge University Press, Cambridge, England, 2000).