# Supernova Limits on Muonic Dark Forces

Claudio Andrea Manzari  0000-0001-8114-3078<EMAIL_ADDRESS>Berkeley Center
for Theoretical Physics, Department of Physics, University of California,
Berkeley, CA 94720, USA  Theoretical Physics Group, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  Jorge Martin Camalich
0000-0002-5673-4500<EMAIL_ADDRESS>Instituto de Astrofísica de Canarias, C/
Vía Láctea, s/n E38205 - La Laguna, Tenerife, Spain  Universidad de La Laguna,
Departamento de Astrofísica, La Laguna, Tenerife, Spain  CERN, Theoretical
Physics Department, CH-1211 Geneva 23, Switzerland  Jonas Spinner
0000-0003-4046-948X<EMAIL_ADDRESS>Institut für
Theoretische Teilchenphysik, Karlsruhe Institute of Technology, Karlsruhe,
Germany  Institut für Theoretische Physik, Universität Heidelberg, Germany
Robert Ziegler  0000-0003-2653-7327<EMAIL_ADDRESS>Institut für
Theoretische Teilchenphysik, Karlsruhe Institute of Technology, Karlsruhe,
Germany  Physikalisches Institut, Albert-Ludwigs-Universität Freiburg, 79104
Freiburg, Germany

(August 28, 2024)

###### Abstract

Proto-neutron stars formed during core-collapse supernovae are hot and dense
environments that contain a sizable population of muons. If these interact
with new long-lived particles with masses up to roughly 100 MeV, the latter
can be produced and escape from the stellar plasma, causing an excessive
energy loss constrained by observations of SN 1987A. In this article we
calculate the emission of light dark fermions that are coupled to leptons via
a new massive vector boson, and determine the resulting constraints on the
general parameter space. We apply these limits to the gauged
$L_{\mu}-L_{\tau}$ model with dark fermions, and show that the SN 1987A
constraints exclude a significant portion of the parameter space targeted by
future experiments. We also extend our analysis to generic effective four-
fermion operators that couple dark fermions to muons, electrons, or neutrinos.
We find that SN 1987A cooling probes a new-physics scale up to $\sim 7$ TeV,
which is an order of magnitude larger than current bounds from laboratory
experiments.

††preprint: CERN-TH-2023-134,TTP23-023, FR-PHENO-2023-06, P3H-23-042

##  I Introduction

Understanding the fundamental nature of Dark Matter (DM), which comprises
$\sim 84\%$ of the matter of the Universe Planck:2018vyg , has become one of
the most pressing problems in contemporary physics (see Ref. Bertone:2010zza
for a review). A wide class of theoretical models describe DM as new light
(sub-GeV) particles, which couple only very weakly to the particles of the
Standard Model (SM). The simplest possible interactions between these two
sectors have been systematically classified within the so-called portal
framework Wilczek:1977pj ; Weinberg:1977ma ; Shifman:1979if ;
Zhitnitsky:1980tq ; Dine:1981rt ; Holdom:1985ag ; Dodelson:1993je ;
Boehm:2003hm ; Asaka:2005an ; Asaka:2005pn ; OConnell:2006rsp ; Patt:2006fw ;
Pospelov:2007mp ; Gorbunov:2007ak ; Batell:2009di ; Krnjaic:2015mbs ;
Arina:2021nqi ; Manzari:2022iyn , giving rise to several benchmark scenarios
for SM interactions with a dark sector that can be tested experimentally (see
Refs. Alekhin:2015byh ; Battaglieri:2017aum ; Beacham:2019nyx ;
Lanfranchi:2020crw ; Krnjaic:2022ozp  for reviews).

Interestingly, if the dark particles are sufficiently light to be produced in
stellar plasmas, then their emission modifies the standard picture of stellar
evolution and stringent constraints on SM interactions with the dark sector
can be obtained also from astrophysical observations Sato:1975vy ;
Dicus:1978fp ; Vysotsky:1978dc ; Dicus:1979ch ; Raffelt:1996wa ;
Balaji:2022noj ; Lella:2022uwi ; Turner:1991ax ; Raffelt:1993ix ;
Caputo:2021rux ; Caputo:2022rca . Ordinary stars yield strong constraints on
DM coupled to electrons and photons Raffelt:1996wa , while the extreme
temperatures and densities reached in the proto-neutron stars (PNS) formed
during core-collapse supernovae (SN) allow one to probe also DM couplings to
nucleons Raffelt:1987yt ; Turner:1987by ; Mayle:1987as ; Burrows:1988ah ;
Burrows:1990pk ; Dent:2012mx ; Rrapaj:2015wgs ; Chang:2018rso ;
Carenza:2019pxu ; Sung:2021swd ; Calore:2021klc ; Choi:2021ign ; Shin:2021bvz
; Calore:2021lih ; Lucente:2022vuo ; Balaji:2022noj ; Lella:2022uwi , pions
Turner:1991ax ; Raffelt:1993ix ; Keil:1996ju ; Carenza:2020cis ;
Fischer:2021jfm ; Shin:2022ulh , hyperons MartinCamalich:2020dfe ;
Camalich:2020wac ; Vonk:2021sit  and muons Bollig:2020xdr ; Calibbi:2020jvd ;
Croon:2020lrf ; Caputo:2021rux . Most of these analyses focus on direct
production and emission of light dark bosons (such as axions or dark photons)
from the stellar medium.

In this article we aim to study instead the case where these bosons merely
serve as massive mediators between dark fermions and SM leptons; i.e. they
serve only as a portal to the dark sector, which allows the production of
sufficiently light dark fermions in stellar plasmas. Prominent examples of
this scenario are gauged lepton flavor models such as
$U(1)_{L_{\mu}-L_{\tau}}$, which contain a light massive gauge boson He:1990pn
; He:1991qd ; Ma:2001md ; Baek:2001kca ; Salvioni:2009jp ; Heeck:2011wj ;
Feng:2012jn ; Harigaya:2013twa ; Bauer:2018onh ; Amaral:2021rzw ;
Greljo:2021npi ; Greljo:2022dwn ; Langacker:2008yv  and a dark sector charged
under the corresponding group Foldenauer:2018zrz ; Asai:2020qlp ;
Holst:2021lzm ; Drees:2021rsg ; Hapitas:2021ilr ; Carpio:2021jhu ;
Heeck:2022znj ; Qi:2022kgs ; Deka:2022ltk ; Baek:2022ozm ; Escudero:2022gez ;
Hooper:2023fqn . This type of scenarios has attracted much attention as they
can simultaneously address the $(g-2)_{\mu}$ anomaly Aoyama:2020ynm ;
Muong-2:2006rrc ; Muong-2:2021ojo ; Borsanyi:2020mff ; Alexandrou:2022amy ;
FermilabLattice:2022izv ; Ce:2022kxy , provide a DM candidate with the right
abundance, and contribute to the effective number of neutrino species,
alleviating the Hubble constant tension Planck:2018vyg ; Brout:2022vxf ;
Bernal:2016gxb ; Escudero:2019gzq ; DiValentino:2021izs . Although we will
present novel constraints on this scenario from SN 1987A later on, we perform
our analysis within a more general setup.

For definiteness, let us consider a vector mediator $Z^{\prime}$ with mass
$m_{Z^{\prime}}$ and couplings to SM leptons and dark fermions $\chi$ with
mass $m_{\chi}$,

$\displaystyle{\cal L}$ $\displaystyle\supset
Z_{\nu}^{\prime}\left(g_{\ell}\overline{\ell}\gamma^{\nu}\ell+g_{\nu_{\ell}}\overline{\nu_{\ell}}\gamma^{\nu}P_{L}\nu_{\ell}+g_{\chi}\overline{\chi}\gamma^{\nu}\chi\right)\,,$
(1)

where $\ell=e,\mu,\tau$, and $g_{\ell}$, $g_{\nu_{\ell}}$ and $g_{\chi}$ are
generic couplings. Assuming that the dark fermions are sufficiently light
($m_{\chi}\lesssim 150\,{\rm MeV}$), we will show in the following that their
production from the PNS in SN 1987A leads to stringent constraints on their
couplings to SM leptons. For the benchmark $U(1)_{L_{\mu}-L_{\tau}}$ model,
this excludes large regions of the parameter space targeted by future
experiments Alekhin:2015byh ; Krnjaic:2019rsv ; Kahn:2018cqs ; Galon:2019owl ;
Asai:2021wzx ; Sieber:2021fue ; Bandyopadhyay:2022klg ; Cesarotti:2022ttv ;
Moroi:2022qwz ; NA64:2022rme ; Rella:2022len ; Ariga:2023fjg .

While our analysis is valid for any mass of the $Z^{\prime}$, we can integrate
it out and describe its contribution with an effective four-fermion operator,
if $m_{Z^{\prime}}$ is much larger than the PNS’s temperature and chemical
potentials. This leads to significant simplifications in the analysis, and
allows us to extend it to derive SN 1987A bounds on completely generic
interactions of the dark sector with leptons through heavy portal mediators.
In the heavy $Z^{\prime}$ limit, our calculations are in fact analogous to the
ones necessary to study the SM production of neutrinos from leptons in stellar
plasmas, which have received much attention over the past half century, since
the seminal works in the 1960’s Chiu:1961zza ; ritus_1961 ;
1963PhRv..129.1383A ; Petrosian:1967alk ; 1967ApJ…150..979B ; Dicus:1972yr ;
1976PhRvD..13.2700S ; 1985ApJ…296..197M ; 1987ApJ…313..531S ;
1989ApJ…339..354I ; 1992ApJ…392…70B ; 1993PhRvD..48.1478B ; 1994ApJ…425..222H
; 1996ApJS..102..411I ; 2003NuPhB.658..217E ; 2004PhRvD..69b3005D ;
2007MNRAS.381.1702K  (see also Refs. 1976PhRvD..13.2700S ; 1977PhRvD..15..977D
; Barbieri:1988av ; Guha:2018mli  for production of light dark fermions from
heavy new physics). On the other hand, in case the $Z^{\prime}$ is a light and
narrow state, the calculation is similar to the on-shell production of massive
vector bosons coupled to leptons in the plasma Croon:2020lrf ;
Kachelriess:2000qc ; Lai:2022csw ; Akita:2022etk ; 2023arXiv230412907L  (see
also Refs. Bollig:2020xdr ; Croon:2020lrf ; Caputo:2021rux ; Lucente:2021hbp
for massive axions).

In this article we focus on the portal interactions with muons, but we also
study neutrinos, which could be naturally linked to muons via $SU(2)_{L}$.
Electrons however form a highly degenerate and ultra-relativistic plasma in
the PNS, which might lead to important medium effects in the electron and
photon dispersion relations, requiring the inclusion of other production
mechanisms not relevant for muons and neutrinos. With this caveat in mind, our
calculations are easily extensible to electrons, generalizing and updating the
pioneering work of Ref. Barbieri:1988av  and improving the results presented
in Ref. Guha:2018mli .

The rest of the paper is organized as follows. In Sec. II we outline the
classical SN argument to constrain new exotic cooling agents using the
neutrino flux observed from SN 1987A. Besides describing the general
theoretical framework, we specify the SN simulations that we employ in our
numerical analysis. In Sec. III we describe and compute the rates of the main
emission mechanisms induced by the model in Eq. (1). We focus on extracting
the main physical features of the rates using different approximations in the
various regimes of the $Z^{\prime}$ mass, and on deriving analytical
estimates. However, our final results rely on exact numerical computations
whose details are deferred to Appendices. In Sec. IV, then, we implement these
calculations of the rates in the SN simulations and derive the constraints on
the parameter space of the $Z^{\prime}$ model in Eq. (1). We also generalize
this analysis in terms of effective operators in the heavy $Z^{\prime}$ limit,
and to one particular realization of the model arising from a
$L_{\mu}-L_{\tau}$ gauged symmetry. Finally, in Sec. V we summarize the
results of our paper and close with a brief outlook.

##  II Supernova Cooling

In the dense and hot environment within proto-neutron stars Burrows:1986me ;
Bethe:1990mw ; 1997PhDT……..18R ; Janka:2012wk  neutrinos become trapped and a
thermal population of muons is predicted to arise Bollig:2017lki ;
Fischer:2020vie . New light dark particles that couple to leptons, e.g. via
the interactions in Eq. (1), can be produced efficiently in the stellar
plasma, leading to a significant loss of energy if they can escape from the
PNS. The corresponding dark luminosity $L_{\chi}$ is then subject to the
classical bound $L_{\chi}\lesssim L_{\nu}$ at $1\,$s post-bounce, where
$L_{\nu}$ is the neutrino luminosity Raffelt:1996wa ; Burrows:1988ah ;
Burrows:1990pk . This limit is obtained from the observation of a neutrino
pulse over $\sim 10\,$s Loredo:2001rx ; Vissani:2014doa  during SN 1987A
Bionta:1987qt ; Hirata:1987hu ; Alekseev:1988gp , which is in accordance with
the predictions of the standard theory of core-collapse SN (see Refs.
Bar:2019ifz ; Page:2020gsx ; Bollig:2020xdr  for a critical reappraisal of
this limit111It has also been noted in Ref. Li:2023ulf  that there is a
disagreement between the results of simulations and the observed neutrino
signal of SN 1987A during the first second. However, this conclusion has been
disputed in Ref. Fiorillo:2023frv , where a more extensive analysis with a
wider time window has been performed.).

Here, we apply this argument to scenarios where light dark fermions couple to
leptons with interactions such as those in Eq. (1). One needs to distinguish
two regimes based on the mean free path (MFP) of the dark fermions in the
plasma or, equivalently, the strength of the portal interactions. If the dark
particles are very weakly coupled (or the MFP is much larger than the radius
of the PNS) then they free stream out from the SN once produced, whereas for
large couplings (or MFP much shorter than the radius of the PNS) they
thermalize with the medium and get trapped inside of the PNS.

In the free-streaming regime the general expression for the total energy-loss
rate per unit volume, $Q$, for a given emission process is

$\displaystyle Q=$
$\displaystyle\int\Big{[}\prod_{\text{init},\,i}\frac{d^{3}p_{i}}{(2\pi)^{3}2E_{i}}f_{i}\Big{]}\Big{[}\prod_{\text{final},\,j}\frac{d^{3}p_{j}}{(2\pi)^{3}2E_{j}}(1\pm
f_{j})\Big{]}$
$\displaystyle\times(2\pi)^{4}\delta^{4}(\sum_{i}p_{i}-\sum_{j}p_{j})\sum_{\rm
spins}|\mathcal{M}|^{2}E_{\chi}\,.$ (2)

These are thermal integrals over the phase space of all the initial- and
final-state particles weighted by their number density distributions $f_{i}$
and the Pauli blocking or Bose enhancement factors $(1\mp f_{j})$,
respectively. Furthermore, $|\mathcal{M}|^{2}$ is the the squared matrix
element of the given production process and $E_{\chi}$ is the total energy
carried away by the dark particles. In the calculations of the free-streaming
regime, one conventionally uses $f_{\chi}=0$ for the new particles, because
their occupation numbers inside the PNS are very low and not thermalized by
assumption.

In the deep trapping regime, on the other hand, the dark-sector particles are
in thermal equilibrium with the plasma and they are emitted from a surface
with radius $r_{\chi}$ (dark sphere) following a law analogous to the one of
the black body radiation,

$\displaystyle L_{\chi}^{\rm
trap}=\frac{\mathfrak{g}_{\chi}}{\pi}r_{\chi}^{2}T_{\chi}^{4}\int_{x_{m}}^{\infty}dx\frac{x^{2}\sqrt{x^{2}-x_{m}^{2}}}{e^{x}+1}\,,$
(3)

where $\mathfrak{g}_{\chi}$ is the number of degrees of freedom of the $\chi$
particle ($\mathfrak{g}_{\chi}=2$ for massive dark fermions),
$x_{m}=m_{\chi}/T_{\chi}$ and $T_{\chi}=T(r_{\chi})$ is the temperature of the
dark sphere. The radius $r_{\chi}$ is defined, as is conventional in
astrophysics 1967pswh.book…..Z ; 1967aits.book…..C ; 2019leas.book…..W ,
through the optical depth $\tau_{\chi}(r)$, by requiring

$\displaystyle\tau_{\chi}(r_{\chi})=\int_{r_{\chi}}^{\infty}\frac{dr}{\lambda(r)}=\frac{2}{3}\,,$
(4)

where $\lambda(r)$ is a suitable spectral average of the dark fermion’s MFP at
a radius $r$. In this work, we use a “naive” thermal average

$\displaystyle\lambda(r)=\langle\lambda(r,p_{\chi})\rangle_{\chi}\equiv\frac{\mathfrak{g}_{\chi}}{n_{\chi}(r)}\int\frac{d^{3}p_{\chi}}{(2\pi)^{3}}\frac{\lambda(r,p_{\chi})}{e^{E_{\chi}/T(r)}+1}\,,$
(5)

for computational simplicity222We have checked that other averages, such as
the conventional Rosseland MFP 1967pswh.book…..Z ; 1967aits.book…..C ;
2019leas.book…..W , give very similar results. In the recent literature also
other forms of taking the spectral average have been discussed, for example by
including an energy-dependent opacity in the free-streaming spectral
luminosity, integrating over the energy only in the very end Croon:2020lrf .
For a discussion and comparison of the various approaches see Ref.
Caputo:2021rux . (see Appendix A for our definitions of thermal averages). The
energy-dependent MFP $\lambda(r,p_{\chi})$ is related to the total rate of
interaction of a dark-sector particle in the medium,
$\Gamma_{\chi}=v_{\chi}/\lambda(r,p_{\chi})$, through its velocity
$v_{\chi}=p_{\chi}/E_{\chi}$.

The contribution to $\Gamma_{\chi}$ of a given process with a bunch of target
particles $b$ colliding with $\chi$ in the initial state is defined through

$\mathcal{C}_{\rm
abs}^{b}=\mathfrak{g}_{\chi}\int\frac{d^{3}p_{\chi}}{(2\pi)^{3}}f_{\chi}\Gamma_{\chi}^{b}\,.$
(6)

The quantity $\mathcal{C}_{\rm abs}$ is the collision operator describing the
absorption rate per unit volume of the medium

$\displaystyle\mathcal{C}_{\rm abs}^{b}=$
$\displaystyle\int\Big{[}\prod_{\text{init},\,i}\frac{d^{3}p_{i}}{(2\pi)^{3}2E_{i}}f_{i}\Big{]}\Big{[}\prod_{\text{final},\,j}\frac{d^{3}p_{j}}{(2\pi)^{3}2E_{j}}(1\pm
f_{j})\Big{]}$
$\displaystyle\times(2\pi)^{4}\delta^{4}(\sum_{i}p_{i}-\sum_{j}p_{j})\sum_{\rm
spins}|\mathcal{M}|^{2}]\,,$ (7)

which uses the same definitions as in Eq. (2), except for $|\mathcal{M}|^{2}$
which is now the squared matrix element of the given absorption
process333Denoting $\Gamma_{\text{ours}}$ for our definition, the standard
approach in the literature Caputo:2021rux  is to work in terms of _emission_
rates $\Gamma_{E}=\mathfrak{g}_{\chi}f_{\chi}\Gamma_{\text{ours}}$. However,
one then has to use Boltzmann-equation arguments to argue which rate has to be
used for absorptive processes, leading to the _reduced absorption_ rate
$\Gamma_{A}$. The definition $\Gamma_{\text{ours}}$ is chosen such that
$\Gamma_{\text{ours}}=\Gamma_{\chi}^{b}=n_{b}\langle\sigma v\rangle_{b}$,
which naturally leads to $\Gamma_{\text{ours}}=\Gamma_{A}$ (see Appendix A)..

For the numerical analyses of this paper we use SN simulations including muons
presented in Ref. Bollig:2020xdr  and whose radial profiles for the relevant
quantities are reported in Garching . Our fiducial results are obtained using
the simulation labeled as SFHo-18.80, which reaches the lowest temperatures
and, therefore, will lead to the most conservative limits on the dark
luminosity (at 1 s post-bounce). The upper bound is set by the neutrino
luminosity calculated within the same simulation, which for SFHo-18.80 is
given by444We use the total co-moving neutrino luminosities reported in the
simulations at the radius of the neutrino-sphere $\sim 16$ km. Note that
including relativistic corrections would reduce both the neutrino and the dark
fermion flux as seen by a distant observer by roughly 30% Caputo:2021rux . We
thank R. Bollig and H-. T. Janka for facilitating us the necessary data to
make these estimates. See also Ref. Caputo:2021rux .

$\displaystyle L_{\chi}\leq L_{\nu}=5.7\times 10^{52}\text{ erg s}^{-1}\,.$
(8)

For a rough estimate of the systematic uncertainties related to SN modeling,
we will also show the more stringent limits obtained from using the hotter
SFHo-20.0 simulation, which gives $L_{\chi}\leq 1.0\times 10^{53}\text{ erg
s}^{-1}$. In the free-streaming regime, the dark luminosity is obtained as a
volume integral of Eq. (2), $L_{\chi}=\int QdV$, while in the trapping regime
we use Eq. (3).

Finally, it will be useful to estimate the contributions to the dark
luminosity of the different processes to understand their relative importance.
For this, we define “typical PNS conditions” as those at 1 s post-bounce and
at a radius $\approx 10$ km. This region dominates the volume emission in
$L_{\chi}$ and is representative of the bounds in the free-streaming regime.
Using the simulation SFHo-18.80 Garching , this approximately corresponds to:

$\boxed{\begin{array}[]{c}\text{\bf Typical PNS conditions}\\\
\;\;\;T=30\,{\rm MeV},\;\rho=2\times 10^{14}$ g cm$^{-3},\\\
\;\;\;\mu_{\mu}=100\,{\rm MeV},\;\mu_{\nu_{e}}=20\,{\rm MeV},\\\
\;\;\;\mu_{\nu_{\mu}}=-10\,{\rm MeV},\;\mu_{e}=130\,{\rm MeV},\\\
Y_{\mu}=0.026,\;\;\;Y_{e}=0.12.\end{array}}$ (9)

Here $T$ denotes the temperature, $\rho$ the density, $\mu_{l}$ the chemical
potential of the lepton $l$ and $Y_{\ell}$ is the number density fraction of
the charged lepton $\ell$ relative to the one of baryons. For the $Y_{\ell}$
we quote the results derived from the rounded temperature and chemical
potentials in Eq. (9) and, therefore, they are slightly different to those
reported in Garching . Let us stress again that Eq. (9) will be only used for
numerical estimates, while our final results and constraints on the models
will be obtained using the full radial profiles of all relevant
thermodynamical quantities.

##  III Production and absorption rates

There are two main production mechanisms of $\chi\bar{\chi}$ pairs from muons
and neutrinos in SN (see top and bottom panel of Fig. 1): (i) Annihilation
$\mu^{-}\mu^{+}\to\chi\bar{\chi}$ and
$\nu_{\ell}\bar{\nu}_{\ell}\to\chi\bar{\chi}$; (ii) photoproduction
$\gamma\mu^{-}\to\mu^{-}\chi\bar{\chi}$. We do not consider bremsstrahlung
processes, $\mu^{-}p\to\mu^{-}p\chi\bar{\chi}$, because they are suppressed
with respect to photo-production (or semi-Compton production) of
(pseudo)scalars and $Z^{\prime}$ Croon:2020lrf ; Bollig:2020xdr ;
Caputo:2021rux .

Note also that this is not generally true for production from electrons
because they are ultra-relativistic and form a highly degenerate system that
suppresses photoproduction compared to bremsstrahlung and annihilation
Raffelt:1996wa . Moreover, there are important plasma effects which, for
example, dress the electron with an effective mass $m_{e}^{*}\sim 10\text{
MeV}$ and give rise to pseudo-particle excitations that need to be taken into
account in a realistic analysis (see e.g. Ref. Lucente:2021hbp  for the
emission of massive axions from electrons in SN). Nonetheless, the
calculations we present in this work can be easily extended to electrons and
compared with previous literature where all these effects have been neglected
Barbieri:1988av ; Guha:2018mli . We will estimate some of them below in Sec.
III.3.

In case of absorption, there are the inverse processes
$\mu^{-}\chi\bar{\chi}\to\gamma\mu^{-}$ and $\chi\bar{\chi}\to\mu^{-}\mu^{+}$,
whose rates are related by detailed balance to those of the photoproduction
and annihilation production, respectively, provided that the $\chi$ and
$\bar{\chi}$ particles reach thermal equilibrium. In addition, other
scattering processes may contribute to the diffusion and energy transport in
the trapping regime, such as $\chi\mu^{-}\to\chi\mu^{-}$ and
$\chi\nu_{\ell}\to\chi\nu_{\ell}$ (see middle panel of Fig. 1), or processes
in the dark sector such as $\chi\bar{\chi}\to\chi\bar{\chi}$.

In Appendices B and C we provide the cross sections for all relevant $2\to 2$
and $2\to 3$ processes needed to calculate the energy-loss and absorption
rates. In the following, we discuss in detail the contributions of the
annihilation and photoproduction topologies.

###  III.1 Annihilation

Figure 1: Most relevant $\ell-\chi$ processes for supernova cooling where
$\ell$ is a charged lepton. Top-to-bottom panels represent annihilation,
scattering and photoproduction processes, respectively. Annihilation and
scattering diagrams will also contribute when $\ell$ is replaced by a
neutrino.

The energy-loss rate per unit volume in Eq. (2) for
$\mu^{-}\mu^{+}\to\chi\bar{\chi}$ annihilation can be simplified to

$\displaystyle Q$
$\displaystyle=\frac{\mathfrak{g}_{\mu}^{2}}{16\pi^{4}}\int_{m_{\mu}}^{\infty}dE_{+}\int_{m_{\mu}}^{\infty}dE_{-}\left(E_{+}+E_{-}\right)\overline{p}_{+}\overline{p}_{-}f_{+}f_{-}I_{s}\,.$
(10)

In this equation, $\mathfrak{g}_{\mu}$ are the muon’s spin degrees of freedom,
$E_{-(+)}$ denote the muon (anti-muon) energy in the PNS rest frame,
$\overline{p}_{\pm}$ are the absolute values of their 3-momenta,
$\overline{p}_{\pm}\equiv\sqrt{E_{\pm}^{2}-m_{\mu}^{2}}$, and $f_{{\pm}}\equiv
f_{{\pm}}(E_{\pm})=1/(e^{(E_{\pm}\pm\mu)/T}+1)$ are their Fermi-Dirac
distributions in the medium. The function $I_{s}$ is a (dimensionless) angular
integral over the annihilation cross section
$\sigma(s)=\sigma({\mu^{+}\mu^{-}\to\chi\overline{\chi}})$

$\displaystyle I_{s}=\int d\cos\theta\,s\sqrt{1-4m_{\mu}^{2}/s}\,\sigma(s)\,,$
(11)

which depends on the colliding angle $\theta$ between the muon and anti-muon
through the Mandelstam variable $s$ in the PNS frame,
$s=2\left(m_{\mu}^{2}+E_{+}E_{-}-\overline{p}_{+}\overline{p}_{-}\cos\theta\right)$.
The cross section is physical only above the 2-particle threshold which
imposes the kinematic constraint $s\geq 4\max({m_{\chi}^{2},m_{\mu}^{2}})$ in
the angular integral.

The annihilation cross section in Eq. (11) for the $Z^{\prime}$ model in Eq.
(1) is

$\displaystyle\sigma(s)=\frac{g_{\mu}^{2}g_{\chi}^{2}}{3\mathfrak{g}_{\mu}^{2}\pi}\frac{s}{(s-m_{Z^{\prime}}^{2})^{2}+m_{Z^{\prime}}^{2}\Gamma_{Z^{\prime}}^{2}}\frac{\beta_{\chi}(s)}{\beta_{\mu}(s)}\kappa_{\mu}(s)\kappa_{\chi}(s)\,,$
(12)

where we have introduced $\beta_{i}(s)=\sqrt{1-4m_{i}^{2}/s}$,
$\kappa_{i}(s)=1+2m_{i}^{2}/s$ and the total $Z^{\prime}$ width

$\displaystyle\Gamma_{Z^{\prime}}=\frac{m_{Z^{\prime}}}{12\pi}\sum_{i}g_{i}^{2}\kappa_{i}(m_{Z^{\prime}}^{2})\beta_{i}(m_{Z^{\prime}}^{2})\theta(m_{Z^{\prime}}-2m_{i})\,,$
(13)

and where $\theta(x)$ is the Heavyside step function.

Note that the average energy carried away by the dark fermions in the
annihilation process is equal to the thermally averaged center-of-mass (CM)
energy of the leptons,

$\displaystyle\mathcal{E}_{l}\equiv\sqrt{\langle s\rangle_{l\bar{l}}}.$ (14)

For typical conditions in the PNS, Eq. (9), $\mathcal{E}_{\mu}\sim 280\,{\rm
MeV}$, which sets the scale for $2m_{\chi}$ above which the production of
$\chi$’s in the plasma will become exponentially (“Boltzmann”) suppressed by
the distribution functions $f_{\pm}$. In addition, this energy scale allows
one to define three regimes of $m_{Z^{\prime}}$ in Eq. (12) depending on which
term dominates the denominator: (i) A “heavy regime” in which
$m_{Z^{\prime}}\gtrsim 1\,{\rm GeV}\gg\mathcal{E}_{\mu}$, so that the
$Z^{\prime}$ is too heavy to be produced on-shell; (ii) the “resonant regime”
where the $Z^{\prime}$ can be produced on-shell, $2m_{\mu}\leq
m_{Z^{\prime}}\lesssim\mathcal{E}_{\mu}$; and (iii) the “light regime” with
$Z^{\prime}$ masses below the two-muon threshold, $m_{Z^{\prime}}<2m_{\mu}$,
so that the $Z^{\prime}$ is too light to be produced on-shell.

Analogous expressions can be defined for electrons and neutrinos by replacing
the couplings and masses accordingly (notice that for neutrinos
$\mathfrak{g}_{\nu_{\ell}}=1$). Moreover, the average CM energies in these
cases are $\mathcal{E}_{e}\sim 160\,{\rm MeV}$ and $\mathcal{E}_{\nu}\sim
130\,{\rm MeV}$, and analogous regimes to those for the muons can be
formulated for neutrino-antineutrino and electron-positron annihilation. The
light regime in these cases is restricted to extremely small $Z^{\prime}$
masses, making it irrelevant for the range of vector boson masses we are
considering here.

The demarcation of these regimes is useful because one can use approximations
to derive analytic results and isolate the main physical factors in control.
In the following, we discuss these approximations and describe their
contributions to the absorption rate $\Gamma_{\chi}$.

####  III.1.1 Heavy regime

In this case, the denominator in the propagator (see Eq. (12)) is dominated by
the $Z^{\prime}$ mass. Expanding in powers of $s/m_{Z^{\prime}}^{2}$ up to
leading order, the cross section can be easily integrated analytically, giving
a function $I_{s}(E_{+},E_{-},m_{\chi},m_{\mu})$ proportional to the effective
coupling $g_{\chi}^{2}g_{\mu}^{2}/m_{Z^{\prime}}^{4}$. We can further
approximate this expression by taking the high-energy limit $m_{\chi}\to 0$
and $m_{\mu}\to 0$, obtaining555In the following discussion we fix
$\mathfrak{g}_{\mu}=2$ with the understanding that some intermediate formulas
change by factors of 2 for the neutrino case.

$\displaystyle I_{s}^{\rm heavy}(E_{+},E_{-},0,0)$
$\displaystyle=\frac{8g_{\chi}^{2}g_{\mu}^{2}}{9\pi}\frac{E_{+}^{2}E_{-}^{2}}{m_{Z^{\prime}}^{4}}\,.$
(15)

Also setting $m_{\mu}\to 0$ in the integrals in Eq. (10), the integrations can
be carried out analytically, giving

$\displaystyle Q^{\rm heavy}$
$\displaystyle=\frac{2g_{\chi}^{2}g_{\mu}^{2}}{9\pi^{5}}\frac{T^{9}}{m_{Z^{\prime}}^{4}}\left[H_{4}(y)H_{3}(-y)+(y\to-y)\right]\,.$
(16)

Here we have re-scaled the chemical potential $y=\mu/T$ and introduced the
functions

$\displaystyle H_{n}(y)=\int_{0}^{\infty}dx\frac{x^{n}}{e^{x-y}+1}=-n!\,{\rm
Li}_{n+1}(-e^{y})\,,$ (17)

where ${\rm Li}_{n+1}(z)$ is the polylogarithm of order $n+1$. If we also take
vanishing chemical potentials, we recover the results in Ref. Barbieri:1988av

$\displaystyle Q^{\rm
heavy}_{0}=\frac{4g_{\chi}^{2}g_{\mu}^{2}}{9\pi^{5}}\frac{T^{9}}{m_{Z^{\prime}}^{4}}F_{4}\,F_{3}\,,$
(18)

in terms of the Riemann $\zeta$-function

$\displaystyle F_{n}=H_{n}(0)=n!(1-2^{-n})\zeta(n+1)\,,$ (19)

with $F_{4}F_{3}\approx 133$. We have also included a subindex in $Q$ to
indicate that this is a zeroth-order approximation neglecting masses and
chemical potentials of the leptons.

In order to assess the accuracy of the above approximations, we compare $Q$ in
Eq. (10) for massless dark fermions and the cross section in the heavy
$Z^{\prime}$ limit with Eq. (18) for different SM particles at the typical
conditions of PNS in Eq. (9). For muons one finds $Q_{\mu}/Q_{0}\approx 0.33$
while for neutrinos and electrons one finds $Q_{\nu}/Q_{0}\approx 0.99$ and
$Q_{e}/Q_{0}\approx 0.54$ (using the physical electron mass in vacuum),
respectively.

The thermal suppression of the muon population is mild for these conditions in
the PNS, $Y_{e}\simeq 4Y_{\mu}$. The positron abundance is also suppressed by
the large electronic chemical potential and, hence, for the same couplings to
electrons and muons one obtains similar rates.

With these approximations one can estimate the parametric dependence of the
energy loss rate per unit mass (i.e. the emissivity) produced by lepton
annihilation in the heavy regime as

$\displaystyle\epsilon_{\rm ann}^{\rm heavy}=\epsilon_{\rm
max}\left(\frac{T}{30\,{\rm
MeV}}\right)^{9}\left(\sqrt{g_{\chi}g_{l}}\;\frac{4.1\,{\rm
TeV}}{m_{Z^{\prime}}}\right)^{4}\,,$ (20)

where we have divided Eq. (18) by the density in Eq. (9), and $\epsilon_{\rm
max}=2.1\times 10^{19}$ erg s-1 g-1 has been estimated dividing $L_{\chi}$ in
(8) by the total mass of the PNS in this simulation $M_{\rm
PNS}=1.351\leavevmode\nobreak\ M_{\odot}$.

####  III.1.2 Resonant regime

If the $Z^{\prime}$ can be produced on-shell, then the denominator in Eq. (12)
is dominated by the $Z^{\prime}$ width $\Gamma_{Z^{\prime}}$, and it can be
replaced by
$\pi/(m_{Z^{\prime}}\Gamma_{Z^{\prime}})\delta(s-m_{Z^{\prime}}^{2})$ in the
narrow width approximation. The $\delta$-function can be used to perform the
angular integration in Eq. (11) and for $m_{\chi}\to 0$ this gives (neglecting
terms of relative size $2m_{\mu}^{2}/m_{Z^{\prime}}^{2}$)

$\displaystyle I_{s}^{\rm res}$
$\displaystyle=\frac{g_{\chi}^{2}g_{\mu}^{2}}{24}\frac{m_{Z^{\prime}}^{3}}{\Gamma_{Z^{\prime}}E_{+}E_{-}}\,.$
(21)

The energy integrations of Eq. (10) can be well-approximated by neglecting the
chemical potentials, but keeping a non-zero muon mass in the integrand, giving

$\displaystyle Q^{\rm res}_{\mu=0}$ $\displaystyle=\frac{g_{\mu}^{2}{\rm
BR}_{\chi}}{4\pi^{3}}m_{Z^{\prime}}^{2}T^{2}m_{\mu}e^{-2m_{\mu}/T}\,,$ (22)

where ${\rm BR}_{\chi}\equiv{\rm BR}(Z^{\prime}\to\chi\overline{\chi})$
denotes the invisible $Z^{\prime}$ branching ratio, and numerically
$m_{\mu}/Te^{-2m_{\mu}/T}\approx 0.004$. For muon annihilation this indeed
yields a good approximation, with $Q_{\mu}/Q_{\mu=0}\approx 0.95$, while for
electrons one can set $m_{e}\to 0$ in the integrals, giving a result similar
to Eq. (16), which can be further approximated by666For $y\gg 1$, one has
$H_{n}(y)\approx y^{n+1}/(n+1),H_{n}(-y)\approx e^{-y}n!$.

$\displaystyle Q^{\rm res}_{m=0}$ $\displaystyle=\frac{g_{\mu}^{2}{\rm
BR}_{\chi}}{16\pi^{3}}m_{Z^{\prime}}^{2}T\mu_{e}^{2}e^{-\mu_{e}/T}\,,$ (23)

and numerically $\mu_{e}^{2}/4T^{2}e^{-\mu_{e}/T}\approx 0.06$, resulting in
$Q_{e}/Q_{m=0}\approx 1.6$. Finally neglecting lepton masses and chemical
potential simplifies to

$\displaystyle Q^{\rm res}_{0}=\frac{g_{\mu}^{2}{\rm
BR}_{\chi}}{4\pi^{3}}m_{Z^{\prime}}^{2}T^{3}\,F_{1}F_{0}\,,$ (24)

where numerically $F_{1}F_{0}\approx 0.57$. This is a good approximation for
neutrino annihilation with $Q_{\nu}/Q_{0}\approx 0.95$, while energy-loss
rates for electrons (muons) are smaller by a factor 10 (100).

Importantly, the contribution to the energy-loss rates in the resonant regime
scales perturbatively with the couplings as $\sim\mathcal{O}(g^{2})$ instead
of $\sim\mathcal{O}(g^{4})$ in the heavy or light regimes. In fact, for ${\rm
BR}_{\chi}=1$, one should recover the results obtained for the coalescence
$Z^{\prime}$ production mechanisms in Ref. Croon:2020lrf 777Indeed, using the
approximation $f_{+}(E_{+})f_{-}(E_{-})\approx f_{Z^{\prime}}(E_{+}+E_{-})$,
we reproduce their expression for the $Z^{\prime}$ production rate in the
light $Z^{\prime}$ and for massless leptons and $\chi$, up to a factor $2/3$.
See also footnote 3 in Ref. Caputo:2021rux  for a discussion of the possible
origin of these discrepancies.. Also, the various rates in Eqs. (22), (23) and
(24) all scale quadratically with the $Z^{\prime}$ mass.

From these expressions one can readily obtain the emissivities, which for
neutrino annihilation read

$\displaystyle\epsilon_{\rm ann,\nu}^{\rm res}=\epsilon_{\rm
max}\left(\frac{T}{30\,{\rm
MeV}}\right)^{3}\left(\frac{g_{\nu_{\ell}}}{10^{-9}}\;\frac{m_{Z^{\prime}}}{10\;{\rm
MeV}}\right)^{2}{\rm BR}_{\chi}\,,$ (25)

where we have used the same approximations as in Sec III.1.1 which are valid
up to $m_{Z^{\prime}}\sim 200$ MeV. As discussed above, emissivities for
electrons and muons are expected to be smaller by a factor 10 and 100,
respectively.

####  III.1.3 Light regime

In this case the denominator of the propagator is dominated by $s$, and the
cross section can be integrated analytically. For massless $\chi$ and muons
one obtains

$\displaystyle I_{s}^{\rm light}$
$\displaystyle=\frac{g_{\chi}^{2}g_{\mu}^{2}}{6\pi}\,,$ (26)

which is independent of $E_{+},E_{-}$. Similarly as in the resonant regime,
the energy integrations of Eq. (10) can be approximated by neglecting the
chemical potentials, but keeping a non-zero lepton mass in the integrand

$\displaystyle Q^{\rm light}_{\mu=0}$
$\displaystyle=\frac{g_{\chi}^{2}g_{\mu}^{2}}{12\pi^{5}}T^{2}m_{\mu}^{3}e^{-2m_{\mu}/T}\,,$
(27)

where $m_{\mu}^{3}/T^{3}e^{-2m_{\mu}/T}\approx 0.04$. For muon annihilation
this indeed yields a good approximation, with $Q_{\mu}/Q_{\mu=0}\approx 0.85$,
while for electrons one can set $m_{e}\to 0$ in the integrals (but keep non-
zero chemical potentials), giving a result similar to Eq. (23)

$\displaystyle Q^{\rm light}_{m_{\ell}=0}$
$\displaystyle=\frac{g_{\chi}^{2}g_{\mu}^{2}}{72\pi^{5}}T^{2}\mu_{e}^{3}e^{-\mu_{e}/T}\,,$
(28)

and numerically $\mu_{e}^{3}/6T^{3}e^{-\mu_{e}/T}\approx 0.2$, resulting in
$Q_{e}/Q_{m=0}\approx 1.2$. Finally neglecting lepton masses and chemical
potential simplifies to

$\displaystyle Q^{\rm
light}_{0}=\frac{g_{\chi}^{2}g_{\mu}^{2}}{12\pi^{3}}T^{5}F_{2}F_{1}\,,$ (29)

where numerically $F_{2}F_{1}\approx 1.5$. This is a good approximation for
neutrino annihilation with $Q_{\nu}/Q_{0}\approx 0.97$, while for electrons
(muons) the energy-loss rates and emissivities are smaller by a factor 10
(100) than predicted by this formula. Nevertheless, using Eq. (29) we obtain
the emissivity

$\displaystyle\epsilon_{\rm ann}^{\rm light}=\epsilon_{\rm
max}\left(\frac{T}{30\,{\rm
MeV}}\right)^{5}\left(\frac{\sqrt{g_{\chi}g_{\mu}}}{3\times
10^{-5}}\right)^{4}\,.$ (30)

####  III.1.4 Annihilation and scattering contributions to trapping

Given the scattering of a $\chi$ with another particle $b$, the absorption
rate can be approximated by (see Appendix A)

$\displaystyle\Gamma_{\chi}^{b}\approx\left(\prod_{i}F_{{\rm
deg},i}\right)\mathfrak{g}_{b}\int\frac{d^{3}p_{b}}{2\pi^{3}}f_{b}\sigma(s)v\,,$
(31)

where $\mathfrak{g}_{b}$ are the number of degrees of freedom of the particle
$b$ ($\mathfrak{g}_{b}=2$ for $b=\bar{\chi},\mu,e$ and $\mathfrak{g}_{b}=1$
for $b=\nu_{\ell}$), $v=\sqrt{(p_{\chi}\cdot
p_{b})^{2}-m_{\chi}^{2}m_{b}^{2}}/E_{\chi}E_{b}$ is the Møller velocity and
$\sigma(s)$ is the scattering cross section for $\chi+b\to X_{1}+\ldots+X_{n}$
Gondolo:1990dk . Moreover, the index $i$ runs over the final-state particles
and we have approximated the effect of the Pauli blocking by its thermal
average or degeneracy factors Raffelt:1996wa

$F_{\text{deg},i}=\langle
1-f_{i}\rangle_{i}=\frac{\mathfrak{g}_{i}}{n_{i}}\int\frac{d^{3}p_{i}}{(2\pi)^{3}}f_{i}(1-f_{i})\,,$
(32)

where $\mathfrak{g}_{i}$ and $n_{i}$ denote their degrees of freedom and
number densities, respectively.

There are two types of processes related by crossing to the annihilation
diagram that are relevant for the trapping regime: Inverse annihilation,
$\chi\,\overline{\chi}\to\mu^{-}\,\mu^{+}$ and scattering,
$\chi\,\mu^{-}\to\chi\,\mu^{-}$ and
$\bar{\chi}\,\mu^{-}\to\bar{\chi}\,\mu^{-}$. Scattering processes are
kinematically more involved as they exchange a $Z^{\prime}$ in the
$t$-channel. For very light $Z^{\prime}$’s ($m_{Z^{\prime}}\ll T$), they have
a differential cross section with a Coulombian enhancement in the forward
direction which involves a small momentum transfer and, therefore, little
contribution to the thermalization rate between the dark and SM sectors.

The interplay between the contributions of inverse annihilation and scattering
to the absorption rate is similar to the case of heavy-lepton neutrinos in SN
Raffelt:2001kv , as recently emphasized in Ref. DeRocco:2019jti . In the
absence of self-interactions between the dark fermions, these two processes
really define two surfaces that determine different properties of the dark
luminosity in the trapping regime. The freeze-out of inverse annihilation
first fixes the number flux of $\chi$’s. The outgoing flow then thermalizes
via scattering processes with the leptons until they decouple at a larger
radius. For $Z^{\prime}$ masses in the resonant regime, absorption rates will
be dominated by the inverse annihilation and the two surfaces coalesce into a
single dark sphere that determines the thermal emission of the $\chi$’s. In
the heavy regime both mechanisms can be important and this distinction must be
kept in mind.

This is also reminiscent of models with new neutrino self-interactions and
their effect in the dynamical evolution of the SN Manohar:1987ec ;
Dicus:1988jh ; Bustamante:2020mep ; Chang:2022aas ; Cerdeno:2023kqo ;
Fiorillo:2023cas ; Fiorillo:2023ytr . Notably, when considering the model in
Eq. (1) with neutrino interactions, processes such as
$\nu_{\ell}\bar{\nu}_{\ell}\to\nu_{\ell}\bar{\nu}_{\ell}$ also occur. These
effects may lead to a fundamentally different incarnation of the SN 1987A
cooling limit, valid in the trapping regime, which is however beyond the scope
of the classical SN 1987A cooling bound that we apply in our analysis.

With all this in mind, our fiducial analysis in the trapping regime includes
both processes, inverse annihilation and scattering, in the calculation of the
MFP to obtain one single dark sphere that determines $L_{\chi}$. However, we
repeat the calculations for the case where we do not include the scattering
processes. We take the variation of our results as an indicator of the
potential systematic uncertainties involved in our treatment of the trapping
regime.

In addition, dark elastic scattering, $\chi\bar{\chi}\to\chi\bar{\chi}$, may
become relevant in the trapping regime, as recently discussed in Refs.
Darme:2020sjf ; Sung:2021swd . In our case, however, this will not play an
important role for the calculations of $L_{\chi}$ in the trapping regime. The
reason is that dark elastic scattering does not directly contribute to
maintaining the population of $\chi$’s in thermal equilibrium with the SM
plasma after they freeze out (after inverse annihilation turns off). Moreover,
for low $m_{Z^{\prime}}$, the rate of dark elastic scattering is resonant and
overwhelmingly larger than scattering or inverse photoproduction (discussed
below in Sec. III.2.3), that do tend to maintain the thermal equilibrium
between the dark and SM sectors. This situation is particularly relevant for
muons, whose population drops significantly at the outer layers of the PNS
where freeze out of the $\chi$’s occurs. Therefore, they would effectively
decouple immediately after, except for the contributions to the MFP induced by
their interactions with the leptons, which are already accounted for in our
fiducial analysis888A proper treatment of the thermalization rates of dark
matter with self interactions is, actually, an important issue when studying
the structure predicted by these models at the center of galactic halos (see
e.g. Ref. Tulin:2017ara ).. Nevertheless, we have studied the contributions of
dark elastic scattering processes for completeness, and discuss in Sec. IV the
consequences for our results if these are included at face value in the
calculation of $\Gamma_{\chi}$ and the MFP.

###  III.2 Photoproduction

For muons we use two approximations: (1) Describe the Pauli blocking of the
final muon by $(1-f_{i})\to F_{{\rm deg},i}$ in Eq. (32) and, (2) in the phase
space integrals of the initial particles we take the extreme non-relativistic
limit where the muon is static and recoilless. Thus, the kinematics is
evaluated in the muon’s rest frame, $s=m_{\mu}^{2}+2m_{\mu}\omega$, with
$\omega$ the photon’s energy, and $E_{\chi}+E_{\bar{\chi}}=\omega$. One
arrives at the conventional formula Raffelt:1996wa

$\displaystyle Q_{\gamma}^{\rm nr}=\frac{n_{\mu}\,F_{{\rm
deg},\mu}}{\pi^{2}}\int_{\omega_{0}}^{\infty}d\omega\omega^{3}f_{\gamma}\,\sigma(s)\,,$
(33)

where $\omega_{0}=(s_{0}-m_{\mu}^{2})/2m_{\mu}$ and $s_{0}$ is the kinematic
threshold of the process.

While Eq. (33) is approximately valid for muons, electrons are instead ultra-
relativistic in the PNS. If we neglect the electron’s mass in the phase space
integrals999Since the photoproduction cross sections have a kinematic
singularity in the limit $m_{e}\to 0$, we keep the physical value of the mass
in numerical estimates, neglecting also plasma effects. and assume that
$E_{\chi}+E_{\bar{\chi}}\approx(\omega+\omega^{\prime})/2$, then

$\displaystyle Q_{\gamma}^{\rm r}$ $\displaystyle=\frac{F_{{\rm
deg},e}}{8\pi^{4}}\int^{\infty}_{0}d\omega\omega
f_{\gamma}\int^{\infty}_{0}d\omega^{\prime}\omega^{\prime}(\omega+\omega^{\prime})f_{e}$
$\displaystyle\times\int^{+1}_{-1}d(\cos\theta)s\sigma(s)\,,$ (34)

where $\omega^{\prime}$ is the electron energy and
$s=2\omega\omega^{\prime}(1-\cos\theta)$.

Similarly to the annihilation topologies, one can estimate the average energy
available for the production of dark particles in the photoproduction process.
For this we use the average CM energy at threshold,

$\displaystyle\mathcal{E}_{\ell}^{\gamma}=\sqrt{\langle
s\rangle_{\ell\bar{\ell}}}-m_{\ell},$ (35)

and one can also define different regimes of $m_{Z^{\prime}}$ for
photoproduction: (1) the “heavy regime” where $m_{Z^{\prime}}\gtrsim 1\;{\rm
GeV}\gg\mathcal{E}_{\ell}^{\gamma}$; (2) the “resonant regime” with
$m_{Z^{\prime}}\lesssim\mathcal{E}_{\ell}^{\gamma}$, where the $Z^{\prime}$
can be produced on-shell. For the conditions specified in Eq. (9) we find
$\mathcal{E}_{\mu}^{\gamma}\sim 90\,{\rm MeV}$ and
$\mathcal{E}_{e}^{\gamma}\sim 150\,{\rm MeV}$.

In Appendix C we provide the results for the photoproduction cross sections
along with some details of the calculations. In addition, we have checked the
validity of the approximations in Eqs. (33) and (III.2) by calculating the
full 5-body phase space thermal integrals exactly, as described in the seminal
work of Ref. Petrosian:1967alk , where this was done for effective
(axial)vector operators. For completeness, we outline this calculation in
Appendix D. We find that the approximate formulas provide a very accurate
description, within $\sim$20% (muons) and $\sim 40\%$ (electrons) of the full
photoproduction rates in the resonant regime and neglecting degeneracy of the
final state lepton. In the heavy regime, Eq. (33) underestimates the full
result for muons by a factor $\sim 8$. Nevertheless, as shown below in Sec.
III.2.1, the emission rate for heavy $m_{Z^{\prime}}$ is dominated by
annihilation by orders of magnitude. On the other hand, for the highly
degenerate electrons in the PNS, we find that approximating the Pauli blocking
effects by $F_{\rm deg,e}$ overestimates the full energy loss rate by a factor
$\sim 3$. In the case of muons, the degeneracy factor $F_{\rm deg,\mu}$
instead agree with the exact result within $\sim 5\%$ accuracy.

In the following, we discuss in more detail the two regimes of photoproduction
and the contribution of inverse photoproduction to the interaction rate
$\Gamma_{\chi}$.

####  III.2.1 Heavy regime

The cross section of the process $\mu^{-}\gamma\to\mu^{-}\chi\bar{\chi}$,
induced by the exchange of a heavy $Z^{\prime}$, is

$\displaystyle\sigma(s)$ $\displaystyle=\frac{\alpha
g_{\mu}^{2}g_{\chi}^{2}m_{\mu}^{2}}{1728\pi^{2}m_{Z^{\prime}}^{4}}\frac{1}{\hat{s}^{2}(\hat{s}-1)^{3}}$
(36)
$\displaystyle\times\Big{(}-55\hat{s}^{6}+682\hat{s}^{5}+483\hat{s}^{4}-968\hat{s}^{3}-169\hat{s}^{2}+30\hat{s}$
$\displaystyle-3+12\hat{s}^{2}\big{(}2\hat{s}^{4}-14\hat{s}^{3}-87\hat{s}^{2}-52\hat{s}+1\big{)}\log\hat{s}\Big{)}\,,$

with $\hat{s}=s/m_{\mu}^{2}$, $\alpha$ denotes the fine-structure constant and
we have taken for simplicity the massless limit $m_{\chi}\to 0$. This cross
section is equivalent to the expression obtained by Dicus in Dicus:1972yr  for
the photoproduction of a neutrino pair using only vectorial couplings101010In
fact, our results agree with Ref. Dicus:1972yr  but disagrees with Ref.
Raffelt:1996wa  which cites Ref. Dicus:1972yr  with the opposite sign in the
term $\propto 120\hat{s}/(\hat{s}-1)^{2}$ in Eq. (3.12). We thank Georg
Raffelt for confirming this typo. See App. C for more details.. One might
attempt to obtain a more simplified expression of the rate by taking the non-
relativistic limit in Eq. (36). However, the cross section in this case,

$\displaystyle\sigma(\omega)=\frac{2\alpha
g_{\mu}^{2}g_{\chi}^{2}}{105\pi^{2}m_{\mu}^{2}}\frac{\omega^{4}}{m_{Z^{\prime}}^{4}}\,,$
(37)

grows rapidly with energy and leads to a gross overestimation of the integral
in Eq. (33) (see also Ref. Petrosian:1967alk ). For instance, for the typical
PNS conditions used above, one obtains a rate that is larger than the one
obtained using the relativistic expression of the cross section by a factor
$\sim 25$.

Nevertheless we use the non-relativistic approximation together with Eq. (33)
to get a rough estimate of the emissivity produced by the photoproduction
process in this regime, giving

$\displaystyle\epsilon_{\gamma}^{\rm heavy}=\epsilon_{\rm
max}\left(\frac{Y_{\mu}}{0.025}\right)\left(\frac{T}{30\;{\rm
MeV}}\right)^{8}\left(\sqrt{g_{\chi}g_{\mu}}\,\frac{{\rm
TeV}}{m_{Z^{\prime}}}\right)^{4}\,.$ (38)

Comparing with the equivalent contribution from annihilation in Eq. (20), we
conclude that photoproduction gives an emissivity rate that is smaller by a
factor $\sim 250$ and thus can be neglected in the heavy regime.

####  III.2.2 Resonant regime

The resonant $\mu^{-}\gamma\to\mu^{-}Z^{\prime}(\to\chi\bar{\chi})$ cross
section is

$\displaystyle\sigma(s)$
$\displaystyle=\frac{\pi\alpha\alpha_{\chi}}{m_{\mu}^{2}}\frac{\text{BR}_{\chi}}{\hat{s}^{2}(\hat{s}-1)^{3}}$
(39)
$\displaystyle\times\Big{(}\big{(}\hat{s}(\hat{s}(\hat{s}+7x^{\prime}+15)+2x^{\prime}-1)-x^{\prime}+1\big{)}R(x^{\prime},\hat{s})$
$\displaystyle+4\hat{s}^{2}\big{(}\hat{s}^{2}-2\hat{s}(x^{\prime}+3)+2x^{\prime}(x^{\prime}+1)-3\big{)}$
$\displaystyle\times\tanh^{-1}\frac{R(x^{\prime},\hat{s})}{\hat{s}-x^{\prime}+1}\Big{)}\,,$

where we have introduced the notation $\alpha_{\chi}=g_{\mu}^{2}/4\pi$, and
where $\hat{s}=s/m_{\mu}^{2}$, $x^{\prime}=m_{Z^{\prime}}^{2}/m_{\mu}^{2}$ and
$R(x,\hat{s})=(\hat{s}^{2}-2\hat{s}(x+1)+(x-1)^{2})^{1/2}$. For
$\text{BR}_{\chi}\to 1$ we recover the cross section for semi-Compton
production of massive vector bosons, $\gamma\mu^{-}\to Z^{\prime}\mu^{-}$. In
the $m_{Z^{\prime}}\to 0$ limit, one obtains

$\displaystyle\sigma(s)$
$\displaystyle=\frac{\pi\alpha\alpha_{\chi}}{m_{\mu}^{2}}\frac{\text{BR}_{\chi}}{\hat{s}^{2}(\hat{s}-1)^{3}}\;\Big{(}\hat{s}^{4}+14\hat{s}^{3}-16\hat{s}^{2}+2\hat{s}-1$
$\displaystyle+\left(2\hat{s}^{4}-12\hat{s}^{3}-6\hat{s}^{2}\right)\log\hat{s}\Big{)}\,,$
(40)

and if we now perform the non-relativistic expansion,

$\displaystyle\sigma(s)\approx\frac{8\pi\alpha\alpha_{\chi}\text{BR}_{\chi}}{3m_{\mu}^{2}}\,,$
(41)

we recover the Thomson cross section for $\alpha_{\chi}\to\alpha$ and
$\text{BR}_{\chi}\to 1$. This is the expression commonly used for the semi-
Compton production of vector particles in stellar plasmas Raffelt:1996wa ;
Croon:2020lrf , but it is less appropriate for the PNS where leptons are
relativistic. In fact, in case of muons for $m_{Z^{\prime}}=0$ and in the
typical conditions we have been using for the PNS, we find that the Thomson
cross section overestimates the relativistic one by a factor $\sim 2$.

On the other hand, the energy-loss rate of the full resonant photoproduction
cross section is insensitive to $m_{Z^{\prime}}$ up to
$m_{Z^{\prime}}\gtrsim\leavevmode\nobreak\ T$, at which point it starts
dropping due to increased Boltzmann suppression and defines the onset of the
heavy regime. Taking as a reference the Thomson cross section, we can estimate
the emissivity of the photoproduction in the resonant regime as

$\displaystyle\epsilon_{\gamma}^{\rm res}=\epsilon_{\rm
max}\left(\frac{Y_{\mu}}{0.025}\right)\left(\frac{T}{30\;{\rm
MeV}}\right)^{4}\left(\frac{g_{\mu}}{5\times 10^{-10}}\right)^{2}\,,$ (42)

using our typical PNS conditions in Eq. (9). Comparing this to the emissivity
from $\mu^{+}\mu^{-}$ annihilation, Eq. (30), we observe that the rate of
$\chi\bar{\chi}$ production from muons for light $Z^{\prime}$ will be
dominated by photoproduction for many orders of magnitude. For
$m_{Z^{\prime}}\lesssim 10$ MeV this process is even more important than
resonant neutrino-antineutrino annihilation for the case
$g_{\nu_{\ell}}=g_{\mu}$, as demonstrated by comparing to Eq. (25).

####  III.2.3 Contribution of inverse photoproduction to trapping

Assuming thermal equilibrium, which is adequate in the trapping regime, the
contribution of inverse photoproduction
$\ell^{-}\chi\bar{\chi}\to\ell^{-}\gamma$ to $\Gamma_{\chi}$ can be related to
the production rates by means of detailed balance (see Appendix A). In case of
muons, the photoproduction rate of $\chi$ per unit volume can be calculated
using the same approximations as in Eq. (33),

$\displaystyle\mathcal{C}_{\rm prod,\mu}^{\gamma}=\frac{n_{\mu}\,F_{{\rm
deg},\mu}}{\pi^{2}}\int_{\omega_{0}}^{\infty}d\omega\omega^{2}f_{\gamma}\,\sigma(s)\,,$
(43)

while for electrons we use

$\displaystyle\mathcal{C}_{\rm prod,e}^{\gamma}$ $\displaystyle=\frac{F_{{\rm
deg},e}}{8\pi^{4}}\int^{\infty}_{0}d\omega\omega
f_{\gamma}\int^{\infty}_{0}d\omega^{\prime}\omega^{\prime}f_{e}$
$\displaystyle\times\int^{+1}_{-1}d(\cos\theta)s\sigma(s)\,.$ (44)

Then, we estimate the contribution of photoproduction to the MFP by applying
detailed balance, using the inverse of $\langle\Gamma_{\chi}^{\gamma}\rangle$
(see Appendix A)

$\displaystyle\langle\Gamma_{\chi}^{\gamma}\rangle_{\chi}=\frac{\mathcal{C}_{\rm
prod}^{\gamma}}{n_{\chi}}\,.$ (45)

Note that this approximation differs from the direct calculation of the
contribution of inverse annihilation and scattering in Eq. (31). In order to
combine the two contributions in our estimate of the MFP we use the
approximate formula

$\displaystyle\lambda=\langle\frac{v_{\chi}}{\Gamma_{\chi}^{\bar{\chi}}+\Gamma_{\chi}^{l}+\Gamma_{\chi}^{\gamma}}\rangle\approx\frac{1}{\langle
v_{\chi}/(\Gamma_{\chi}^{\bar{\chi}}+\Gamma_{\chi}^{l})\rangle^{-1}+\langle\Gamma_{\chi}^{\gamma}\rangle/\langle
v_{\chi}\rangle}\,,$ (46)

where all the thermal averages are understood to be taken with respect to the
$\chi$ kinematics (see Appendix A).

###  III.3 Other processes and neglected plasma effects

In our analysis, we have selected the processes that are dominant for the
production and absorption of $\chi$’s in muons and have neglected plasma
effects which are expected to be small in this case. Electrons in the PNS, on
the other hand, are highly degenerate. Moreover, in the plasma the electron
and photon dispersion relations are significantly modified. The electron mass
effectively increases while the photon acquires a longitudinal mode and an
effective mass that could enable the decay $\tilde{\gamma}\to\chi\bar{\chi}$,
where the “plasmon” $\tilde{\gamma}$ includes these collective plasma modes
1963PhRv..129.1383A ; 1992ApJ…392…70B ; 1994ApJ…425..222H ; Raffelt:1996wa .

Nonetheless, the production of $\chi$’s in the heavy regime and neglecting the
$\chi$ mass would be analogous to the SM pair-production of heavy-lepton
neutrinos from the electrons in the stellar plasma. Plasmon decay is indeed an
important process in the conditions of high densities predicted in the PNS.
However, at the high temperatures reached in the SN explosions considered in
this work, $e^{+}e^{-}$ annihilation becomes the dominant process
1994ApJ…425..222H . Adding mass to the $\chi$’s will not affect this
conclusion and may, in fact, kinematically close the plasmon decay if
$m_{\chi}\gtrsim 10\leavevmode\nobreak\ {\rm MeV}$, which is the scale of the
plasma frequencies expected in the PNS. Finally, accounting for the increase
of the electron mass in the plasma by a similar amount does not affect
significantly the annihilation rates as discussed in Sec. III.1.1.

On the other hand, for the light $Z^{\prime}$ some of the neglected effects
can become important. For instance, resonant bremsstrahlung production
$e^{-}p\to e^{-}p(Z^{\prime}\to\chi\bar{\chi})$ could become the dominant
process for $m_{Z^{\prime}}\sim 10\leavevmode\nobreak\ {\rm MeV}$, as it is
the case for on-shell production of axion-like particles Lucente:2022vuo . In
addition, for these masses one needs to consider medium-induced
$\tilde{\gamma}-Z^{\prime}$ mixing which may also have an impact An:2013yfc .
For all of these reasons, our results for the emission from electrons in the
light regime, $m_{Z^{\prime}}\lesssim 50$ MeV, should be considered as an
intermediate step towards a more refined and robust calculation, and the
results that we report for this case regarded as a rough approximation.

##  IV Results

In this section, we apply the upper limit on the dark luminosity in Eq. (8) to
obtain the SN 1987A constraint on the parameter space of the different dark
sector models we consider. We start by presenting an analysis in terms of
effective operators, which corresponds to the heavy regime introduced in Sec.
II for the calculation of the relevant processes. This allows us to generalize
our analysis to generic interactions of dark fermions coupled to leptons via
dimension-6 operators, and extract a SN 1987A limit for any portal mediator
with mass much larger than the temperatures and chemical potentials in the
PNS. It is interesting to note that this approach was already applied in the
context of neutrino emission from stellar plasmas in the very early days of
the SM 1976PhRvD..13.2700S ; 1977PhRvD..15..977D . We then focus on the
constraints on the parameter space of the simplified $Z^{\prime}$ model of Eq.
(1), which can be regarded as the continuous extension of the EFT bounds to
low mediator masses for one particular operator ($VV$). Finally, we present
the SN 1987A constraints for a phenomenologically relevant and UV-motivated
version of the $Z^{\prime}$ model, obtained by gauging the $L_{\mu}-L_{\tau}$
symmetry.

###  IV.1 Effective field theory

$XY$ |  $\Lambda_{\mu}^{\rm eff}$ [TeV] |  $\Lambda_{\nu_{\mu}}^{\rm eff}$ [TeV] |  $\Lambda_{e}^{\rm eff}$ [TeV]
---|---|---|---
$SS$ | $0.0017-4.4$ | $0.062-5.2$ | $0.070-5.4$
$PS$ | $0.00044-5.1$ | $0.062-5.2$ | $0.070-5.4$
$VV$ | $0.0017-5.7$ | $0.072-5.6$ | $0.11-5.9$
$AV$ | $0.0022-4.7$ | $0.072-5.6$ | $0.11-5.8$
$LL$ | $0.0015-3.7$ | $0.051-4.0$ | $0.074-4.1$
$LV$ | $0.0018-4.4$ | $0.061-4.7$ | $0.088-4.9$
$TT$ | $0.0033-6.8$ | $0.10-6.7$ | $0.17-7.0$
Table 1: SN 1987A exclusion range of the effective scale $\Lambda_{l}^{\rm
eff}\equiv\Lambda_{l}/\sqrt{C_{XY}^{l}}$ for the EFT interactions defined in
Eq. (47), with $m_{\chi}=0$ and for the simulation SFHo-18.8. The lower limit
of the constraint is set by the trapping regime. The bounds for $XY=SP$, $PP$,
$AA$, $VA$, and $T^{\prime}T$ are equal to those of $SS$, $PS$, $VV$, $AV$,
and $TT$, respectively. The bounds for $RR$, $LR$ and $RL$ are identical to
those of the $LL$ operators. The constraints on $\Lambda^{\rm eff}_{\nu_{e}}$
essentially coincide with those on $\Lambda^{\rm eff}_{\nu_{\mu}}$.  Figure 2:
Comparison between SN 1987A and collider limits on the effective scale
$\Lambda_{e}^{\rm eff}\equiv\Lambda_{e}/\sqrt{C_{XY}^{e}}$ for the EFT
interactions of electrons defined in Eq. (47), with $m_{\chi}=0$. The orange
bars show the exclusion range obtained from SN 1987A, using the simulation
SFHo-18.8. The upper limit is obtained by the free-streaming regime and the
lower limit is set by the trapping regime. The red bars show the lower limits
obtained at LEP Fox:2011fx , while the purple and green hatched bars show the
projections for Belle-II with $50\,{\rm ab}^{-1}$ of integrated luminosity
Liang:2021kgw  and ILC with $\sqrt{s}=1\,{\rm TeV}$ and $L_{\rm
int}=1000\,{\rm fb}^{-1}$ Chae:2012bq . Figure 3:  SN 1987A constraints on the
heavy scale $\Lambda_{l}^{\rm eff}$ as a function of the dark fermion mass
$m_{\chi}$ for the $VV$ operator and different leptons. The solid orange
region is excluded using numerical input from the simulation SFHo-18.8, while
the dotted line encloses the region excluded using SFHo-20.0. The hatched
region gives an estimate of the uncertainty of the trapping regime, obtained
by omitting the contribution to the MFP from scattering processes, as
discussed in Sec. III.1.4. The constraints on $\nu_{e}$ and $\nu_{\mu}$ are
essentially similar, so we show those on $\nu_{\mu}$ and label them as $\nu$.

The analysis for the heavy $Z^{\prime}$ can be generalized in the context of
an EFT with four-fermion operators. Focusing on the couplings of dark-sector
fermions $\chi$ to leptons $l=e$, $\mu$, $\nu_{\ell}$, the most general
effective Lagrangian at leading order is

$\displaystyle{\cal L}_{\rm EFT}$
$\displaystyle=\frac{1}{\Lambda_{l}^{2}}\sum_{X,Y}C_{XY}^{l}\left(\overline{l}\,\Gamma_{X}\,l\right)\cdot\left(\overline{\chi}\,\Gamma_{Y}\,\chi\right)\,,$
(47)

where $X,Y$ run over $V,A,S,P,L,R,T,T^{\prime}$, with
$\Gamma_{V}=\gamma^{\mu}$, $\Gamma_{A}=\gamma^{\mu}\gamma_{5}$,
$\Gamma_{S}=\mathbb{1}$, $\Gamma_{P}=\gamma_{5}$,
$\Gamma_{R,L}=\gamma^{\mu}(\mathbb{1}\pm\gamma_{5})/2$,
$\Gamma_{T}=\sigma^{\mu\nu}=i/2[\gamma^{\mu},\gamma^{\nu}]$ and
$\Gamma_{T^{\prime}}=\sigma^{\mu\nu}\gamma_{5}$, and Lorentz indices properly
contracted111111Only two tensor operators are independent and we will take
those corresponding to $C^{l}_{TT}$ and $C^{l}_{T^{\prime}T}$.. Matching the
effective operators to the $Z^{\prime}$ model in Eq. (1) coupled to charged
leptons yields $\Lambda_{\ell}=m_{Z^{\prime}}$ and
$C_{VV}^{\ell}=g_{\ell}g_{\chi}$. For neutrinos their bilinears in the EFT
Lagrangian are constructed with left-handed fields and contribute, instead, to
$C_{LV}^{\nu_{\ell}}=g_{\nu_{\ell}}g_{\chi}$.

In Table 1, we show the limits obtained on $\Lambda_{l}^{\rm
eff}\equiv(\Lambda_{l}/C^{l}_{XY})^{1/2}$ for these interactions in the limit
$m_{\chi}=0$ for muons, neutrinos and electrons. The upper limits correspond
to the free-streaming regime, the lower limits to the trapping regime. Notice
that the excluded regions of $\Lambda_{l}^{\rm eff}$ are in the EFT range of
validity as long as $\Lambda_{l}\gtrsim 1\,{\rm GeV}$, which is much larger
than all other energy scales relevant in the PNS. This is the case for
essentially all operators.

The SN 1987A bounds on the EFT operators are very strong, reaching up to $4-7$
TeV. This sensitivity to the mediator mass scale is approximately one order of
magnitude better than the one achieved by laboratory experiments for similar
leptonic interactions Fox:2011fx ; Alcaide:2019pnf ; Escrihuela:2021mud ;
Liang:2021kgw  121212See also Fernandez-Martinez:2023phj  for an estimate of
the SN 1987A bounds obtained by recasting the results derived in
DeRocco:2019jti  for an operator coupled like the electromagnetic current..
For instance, monophoton searches at LEP have been used to set the following
lower bounds on $\Lambda^{\rm eff}_{e}$ Fox:2011fx ,

$\displaystyle\begin{split}VV,\leavevmode\nobreak\ AA:&\hskip
17.07182pt0.48\text{ TeV}\,,\\\ SS:&\hskip 17.07182pt0.44\text{
TeV}\,.\end{split}$ (48)

In Fig. 2, we show the SN 1987A limits on these EFT operators compared with
those obtained at LEP Fox:2011fx , and with the projections of the sensitivity
that could be achieved at Belle II Liang:2021kgw  or at a future $e^{+}e^{-}$
linear collider Chae:2012bq . Remarkably, the SN 1987A bounds are stronger
than LEP by roughly one order of magnitude, and will even dominate over future
collider limits. Note however that SN 1987A constraints apply only to
sufficiently light dark fermions, while collider bounds typically extend to
larger $\chi$ masses, such as LEP, which provides constraints for
$m_{\chi}\lesssim 100$ GeV.

As discussed in Sec. II, the dominant effect in the heavy regime of the
$Z^{\prime}$ model is annihilation (for typical SN conditions), and
photoproduction can be neglected. We have checked that this is indeed true for
any EFT operator by calculating their contributions to photoproduction
explicitly (see Appendix C). In fact, we note that, for the EFT limit of the
$Z^{\prime}$ model, the approximate expression in the high-energy limit
($m_{\ell}\to 0$) in Eq. (18) leads to a bound on the effective scale that is
off only by $\sim 30\%$ with respect to the full calculation. Therefore, since
the high-energy limit gives quite accurate results, the numerical bounds for
different Lorentz structures and different leptons are of the same order, as
in this limit the corresponding cross sections differ at most by
${\mathcal{O}}(1)$ numbers.

In Fig. 3 we show the dependence of the limits on the heavy scale as a
function of the dark fermion mass $m_{\chi}$ for the effective $VV$
interaction and $LV$ for neutrinos. As expected, the excluded regions shrink
with increasing $m_{\chi}$ and they are limited to masses $m_{\chi}\lesssim
300$ MeV. We also analyze the variations of the constrained region produced by
using the simulation SFHo-20.0 or different prescriptions for the processes
included in the trapping regime. The upper limits of $\Lambda_{l}^{\rm eff}$
obtained from this hotter simulation are a factor $\sim 2$ stronger in the
free-streaming regime, because in this case emission is dominated by the
hottest region inside the supernova. On the other hand, the uncertainties of
the boundary with the trapping regime for electrons and neutrinos are
relatively small and our fiducial calculation is, again, on the conservative
side. The reason for this behavior is that the dark sphere is not located in
the region of highest temperature. This makes trapping sensitive to the shape
of the temperature profile, which is similar for both simulations.

Moreover, the scale marking the onset of the trapping regime in these cases is
$\Lambda_{l}^{\rm eff}\approx 100$ GeV, which is of the order of the
electroweak scale. This is consistent with the fact that the boundary of the
trapping regime is set by the luminosity of the trapped neutrinos via Eq. (8),
which interact with SM leptons precisely through dimension-6 operators
suppressed by the Fermi scale.

For muons, however, the bound extends to scales that are a few orders of
magnitude lower than for electrons and neutrinos. In addition, there is a
large variation in the location of the boundary depending on the selection of
processes contributing to trapping. The reason is that muons are relatively
heavy and there is a maximal radius where they can be produced by thermal
fluctuations in the plasma. Putting it differently, there are no muons to
scatter with in case of inverse photoproduction and inverse annihilation
becomes ineffective because of the strong phase-space suppression. In case
that only $\mu-\chi$ interactions are included in the calculation of the MFP,
the radius of the dark sphere is typically smaller than the one of the
neutrino sphere and the bounds on $\Lambda_{\mu}^{\rm eff}$ extend down to
$\sim 1\leavevmode\nobreak\ {\rm GeV}$. However, if one were to also include
$\chi\bar{\chi}$ elastic scattering in the calculation of the MFP, and the
$\chi$ is light, $m_{\chi}\lesssim T$, then they would produce contributions
analogous in size to those given by the neutrino and electron interactions and
the boundary of the trapping regime would be increased again to $\sim
100\leavevmode\nobreak\ {\rm GeV}$, cf. the discussion in Section III.1.4.

###  IV.2 Simplified $Z^{\prime}$ models

Figure 4:  Excluded parameter space of the $Z^{\prime}$ model in Eq. (1) using
constraints from SN 1987A with the simulation labeled as SFHo-18.8 in Ref.
Bollig:2020xdr . The dotted line encloses the region excluded using SFHo-20.0.
The hatched region gives an estimate of the uncertainty of the trapping
regime, obtained by omitting the contribution to the MFP from elastic
scattering processes, as discussed in Sec. III.1.4.

We now present the SN 1987A constraints in the parameter space of the
simplified $Z^{\prime}$ model in Eq. (1). In the derivation of our results for
the annihilation contributions to $Q$ we numerically solve the integral in Eq.
(10) (and the equivalent ones for neutrinos and electrons), as described in
Appendix E. This procedure allows us to track the contribution of the
annihilation rates across the whole $m_{Z^{\prime}}$ range, including the
transition regions between the three regimes defined in Sec. III.1. Instead,
for photoproduction we only use the expressions in the resonant regime,
because photoproduction is only relevant for the low-$m_{Z^{\prime}}$ regime
of the charged leptons, see Eq. (38).

Figure 5: SN 1987A cooling constraint (orange region) in the mass-coupling
plane of the $L_{\mu}-L_{\tau}$ model, assuming $m_{Z^{\prime}}=3m_{\chi}$ and
$g_{\chi}=g_{\mu-\tau}$. The dotted orange line indicates the bound obtained
from using the SN simulation SFHo-20.0, instead of SFHo-18.8 which gives the
weakest constraints. Also shown are the preferred regions to explain the
$(g-2)_{\mu}$ anomaly (green band) and the $H_{0}$ tension (yellow band),
along with current constraints (gray region) and the forecasts for future
experiments. See main text for details.

In Fig. 4 we show the SN 1987A limits for the model in Eq. (1), for the case
when only one of the leptonic couplings is present and with $g_{\chi}=g_{l}$.
In the lower (upper) part of the plots we identify the dominant process for
production (absorption) in the indicated mass range. We also show the
variation of the constrained region obtained by using the SFHo-20.0 simulation
(dotted curve) and, independently, by omitting the scattering contribution to
the MFP in the trapping regime (hatched region). In all cases we observe the
onset at high masses of the power-law behavior of the constraints $\propto
1/m_{Z^{\prime}}^{4}$, characteristic of the EFT. Interestingly, this occurs
for $m_{Z^{\prime}}\gtrsim 1$ GeV in the free-streaming domain but already for
$m_{Z^{\prime}}\gtrsim 100$ MeV in the trapping regime. This is because in the
former case emission is governed by the conditions in the hottest region of
the PNS, while in the latter it corresponds to the cooler outmost layer of the
dark sphere, where all relevant energy scales are smaller. The boundary of the
trapping regime changes very little with respect to adding or not scattering
contributions in the absorption rate. Adding dark elastic scattering
($\chi\bar{\chi}\to\chi\bar{\chi}$) in the muonic case (for $m_{\chi}\lesssim
T$) would instead make the trapping region similar to the neutrino case.

The shape of the constrained region below $m_{Z^{\prime}}\approx 1$ GeV
depends on the lepton considered. For muons the low $m_{Z^{\prime}}$ region is
dominated by (resonant) photoproduction, which gives a flat bound up to
$m_{Z^{\prime}}\gtrsim T$, where the on-shell production of the $Z^{\prime}$
starts decreasing due to Boltzmann suppression. Nevertheless, it remains more
important than $\mu^{+}\mu^{-}$-annihilation (in the light regime), until it
becomes resonant at the $\mu^{+}\mu^{-}$ threshold, $m_{Z^{\prime}}\geq
2m_{\mu}$. However, this occurs already at large energies and the production
quickly suffers from the Boltzmannian suppression, converging to the EFT
scaling at higher $m_{Z^{\prime}}$.

For electrons, photoproduction dominates again for low $m_{Z^{\prime}}$, even
though the $e^{+}e^{-}$ threshold is much lower than for muons. Resonant
annihilation of electrons quickly replaces photoproduction as the dominant
process above $m_{Z^{\prime}}\gtrsim T$ in the free-streaming domain, until
the resonance starts suffering from Boltzmann suppression and the EFT takes
over. For the trapping regime, there is no range of $m_{Z^{\prime}}$ where
resonant annihilation dominates and the EFT directly replaces inverse-
photoproduction at $m_{Z^{\prime}}\gtrsim 100$ MeV.

For neutrinos there is no photoproduction and the production and absorption
rates are given by annihilation in the resonant and heavy regimes. In the
free-streaming regime we observe a strengthening of the bound up to $T\lesssim
100$ MeV. This is due to the $m_{Z^{\prime}}^{2}$ scaling of the emission rate
in the resonant regime, see Eq. (24), which is quickly overcome by Boltzmann
suppression until the heavy regime takes over.

###  IV.3 Gauged $L_{\mu}-L_{\tau}$ models

Finally, we study SN 1987A constraints on a UV-motivated realization of the
simplified $Z^{\prime}$ model in Eq. (1). This is the gauged
$L_{\mu}-L_{\tau}$ model coupled to DM fermions, described by the interaction
Lagrangian

$\displaystyle{\cal L}_{\rm int}$
$\displaystyle=Z_{\mu}^{\prime}\left(g_{\mu-\tau}j^{\mu}_{\rm
SM}+g_{\chi}\overline{\chi}\gamma^{\mu}\chi\right)\,,$ (49)

where $\chi$ is the dark fermion and $j^{\mu}_{\rm SM}$ is the SM part of the
$L_{\mu}-L_{\tau}$ current. These interactions induce an irreducible
contribution to the kinetic mixing of the $Z^{\prime}$ with the photon,
through $\mu$ and $\tau$ loops, giving a mixing parameter of the order
$\epsilon\sim g_{\mu-\tau}/70$ Kamada:2015era . This would presumably give
only small corrections to our analysis (see e.g. Croon:2020lrf ) and thus will
be neglected.

It is well known that such models can accommodate the present $(g-2)_{\mu}$
anomaly with couplings of order $g_{\mu-\tau}\sim 10^{-4}$ for a light
$Z^{\prime}$ gauge boson, $m_{Z^{\prime}}\ll m_{\mu}$ Fayet:2007ua ;
Pospelov:2008zw ; Chen:2017awl . They also allow to reproduce the DM relic
abundance through resonant $s$-channel annihilation, when muon and DM fermion
have similar couplings to the $Z^{\prime}$ gauge boson and the latter is
heavier than the DM fermion by a factor $2-3$ Holst:2021lzm ; Drees:2021rsg .
For such light masses, $Z^{\prime}$ decays and DM annihilation can heat the SM
bath after neutrino decoupling, thereby increasing the effective number of
relativistic degrees of freedom, usually expressed as an effective number of
neutrino species $N_{\rm eff}$. Such a contribution could help to reduce the
long-standing tension between local and cosmological determinations of the
Hubble constant Bernal:2016gxb ; DiValentino:2021izs , if the new
contributions is of order $\Delta N_{\rm eff}\sim 0.1-0.4$ Planck:2018vyg .

At present, the most relevant laboratory constraints (see Ref. Bauer:2018onh
for an overview) stem from BaBar searches for $Z^{\prime}$-bosons above 212
MeV decaying to muons BaBar:2016sci , neutrino trident production
Altmannshofer:2014pba  at CCFR CCFR:1991lpl , bounds on coherent elastic
neutrino nucleus scattering from the COHERENT collaboration Cadeddu:2020nbr
and constraints on neutrino-electron scattering Harnik:2012ni ; Kaneta:2016uyt
at BOREXINO Bellini:2011rx . A variety of accelerator-based searches have been
proposed to explore the unconstrained parameter space, such as NA62
Krnjaic:2019rsv , which looks for final state radiation of $Z^{\prime}$-bosons
in $K^{+}\to\mu^{+}\nu_{\mu}$, and dedicated searches using muon beam
facilities such as the NA64$\mu$ experiment Gninenko:2020hbd  at CERN and
$M^{3}$ Kahn:2018cqs  at Fermilab.

Astrophysical limits from white dwarfs have been studied already in Refs.
Dreiner:2013tja ; Bauer:2018onh . Here we show that also constraints from SN
1987A disfavor a large region of parameter space with significant overlap with
the expected reach of planned experiments and with the region that could
address the $H_{0}$ tension. In Fig. 5, we show the SN 1987A limits on the
muon coupling as a function of the $Z^{\prime}$ mass, in the scenario where
muons and dark matter couple equally to the $Z^{\prime}$,
$g_{\mu-\tau}=g_{\chi}$, and $m_{Z^{\prime}}=3m_{\chi}$. This is compared to
the current bounds discussed above, shown in gray, and the regions preferred
at $95\%$ C.L. by $(g-2)_{\mu}$ and $H_{0}$.

In the free-streaming regime we use the same method as for the simplified
model described in Appendix E, including the contributions from
$\mu^{+}\mu^{-}$, $\nu_{\mu}\overline{\nu}_{\mu}$ and
$\nu_{\tau}\overline{\nu}_{\tau}$ annihilations. We also add the resonant
photoproduction off muons which dominates the rate up to
$m_{Z^{\prime}}\lesssim 10$ MeV, where neutrino annihilation starts to give
the largest contribution to the rate, producing the characteristic
strengthening of the bound with $m_{Z^{\prime}}$. At
$m_{Z^{\prime}}=2\,m_{\mu}$ a small feature signalizes the onset of resonant
$\mu^{+}\mu^{-}$ annihilation. Instead, the boundary of the constrained region
in the deep trapping regime is dominated by inverse resonant annihilation into
neutrinos. Note that the $\chi$ mass scales with the $Z^{\prime}$ mass, which
has the effect of slightly suppressing the rate as compared to the
massless-$\chi$ case, cf. Fig. 4. A rough estimate of the uncertainty of the
excluded region is indicated with a dashed orange line, which shows the limits
obtained from employing the hottest SN simulation SFHo-20.0, as described in
Sec. II.

##  V Summary and Conclusions

In this paper, we have studied the SN 1987A cooling constraints on dark-sector
models induced by the emission of new light dark fermions $\chi$ coupled to
leptons. To provide a concrete framework, we consider general vector portal
interactions arising from the exchange of a massive $Z^{\prime}$ vector
mediator. We focus primarily on the couplings to muons, which are predicted to
have sizable number densities within the hot and dense environments of the
proto-neutron star formed during core-collapse supernovae. However, we also
extend our analysis to couplings to neutrinos and electrons.

We have considered various mechanisms for the production and absorption of the
$\chi$’s and different regimes that depend on the ranges of the parameters of
the model. Firstly, the constraints depend on the mass of $\chi$, as their
pair production becomes Boltzmann suppressed for $m_{\chi}\gg T$. Secondly,
different regions of $m_{Z^{\prime}}$ can be identified based on whether the
dark fermions are resonantly produced or generated from the tail of the
$Z^{\prime}$ resonance. This distinction arises, for instance, when the
$Z^{\prime}$ is heavy and cannot be produced on-shell by the thermal
fluctuations in the medium. Finally, there exist two distinct regimes of
coupling values, depending on whether the dark fermions free-stream or become
trapped within the PNS. Consequently, by analyzing the $\chi\bar{\chi}$
production within these two regimes, for given masses of the $Z^{\prime}$ and
the $\chi$, we can determine the range of couplings that is excluded by the
observations from SN 1987A.

For $Z^{\prime}$ particles with masses $m_{Z^{\prime}}\lesssim T\sim 10$ MeV
and massless $\chi$, the observations from SN 1987A place constraints on the
couplings between $\sim 10^{-1}$ and $\sim 10^{-9}$, for equal couplings of
the vector mediator to leptons and dark fermions, cf. Fig. 4. However, the
range of these bounds strongly depends on the $Z^{\prime}$ mass and the
specific lepton to which the $Z^{\prime}$ couples.

These calculations can be readily extended to explicit $Z^{\prime}$-models,
for example motivated by a gauged lepton flavor symmetry with a dark sector
charged under it. We specifically investigate the case of $L_{\mu}-L_{\tau}$,
which has been proposed in the literature as a combined solution to the
$(g-2)_{\mu}$ anomaly, the Hubble tension and the dark matter puzzle. The SN
1987A limit covers a large region of parameter space that overlaps with the
forecasts of future experiments and with part of the region that could address
some of the tensions, see Fig. 5.

On the other hand, when the $Z^{\prime}$ mass is larger than the temperature
and chemical potentials in the PNS, the interactions mediated by the
$Z^{\prime}$ can be accurately described by effective operators. This allows
us to generalize the analysis to completely generic heavy portal interactions
between dark fermions and SM leptons, summarized in Fig. 3. We find that SN
1987A cooling can probe new-physics scales up to $4-7$ TeV (cf. Table 1),
which surpasses current bounds from laboratory experiments by an order of
magnitude (see Fig. 2).

We emphasize that our analysis is not complete when a light $Z^{\prime}$ is
coupled to electrons. In this case, bremsstrahlung processes are expected to
provide the dominant contributions to the emission of dark fermions, and
plasma effects can have a significant impact on the analysis. Nevertheless, we
consider our results as an important step into this direction, which
significantly extends previous studies in the literature.

Finally, other aspects of SN physics could lead to constraints complementary
to the ones obtained in this work (see e.g. Refs. Caputo:2022mah ;
Fiorillo:2022cdq ). In particular, there have been recent efforts towards a
better understanding of the effect of neutrino self-interactions Chang:2022aas
; Fiorillo:2023cas ; Fiorillo:2023ytr  in the trapping regime. Moreover, it
would be interesting to check whether bounds obtained from energy transport by
dark particles exceed the bounds derived from energy loss alone, following the
analysis in Ref. Caputo:2022mah .

## Code availability

We provide a minimal code example to test different parameter points of the
$L_{\mu}-L_{\tau}$ model at https://github.com/spinjo/SNforMuTau.git.

## Acknowledgments

We thank Andrea Caputo, Pilar Coloma, Miguel Escudero, Patrick Foldenauer,
Felix Kahlhöfer, Enrico Nardi, Uli Nierste, Filippo Sala and Stefan Vogl for
useful discussions. We want to thank Robert Bollig and Hans Thomas Janka for
sharing data of the simulations with us. This work is partially supported by
project C3b of the DFG-funded Collaborative Research Center TRR257, “Particle
Physics Phenomenology after the Higgs Discovery” and has received support from
the European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska -Curie grant agreement No 860881-HIDDeN. The work of C.A.M.
is supported by the Office of High Energy Physics of the U.S. Department of
Energy under contract DE- AC02-05CH11231. Work by JMC is supported by
PGC2018-102016-A-I00, and the “Ramón y Cajal” program RYC-2016-20672.

##  Appendix A Generalities of absorption rates

In the case of generic $2\to n$ processes, where a dark particle $\chi$
interacts with a particle $b$ in the initial state, the resulting
contributions $\Gamma_{\chi}^{b}$ to the absorption rate of $\chi$ can be
approximated by

$\displaystyle\Gamma_{\chi}^{b}$
$\displaystyle\approx\frac{1}{2\mathfrak{g}_{\chi}E_{\chi}}\prod_{\text{final
}i}F_{\text{deg,}i}\int\frac{d^{3}p_{b}}{(2\pi)^{3}2E_{b}}f_{b}\frac{d^{3}p_{i}}{(2\pi)^{3}2E_{i}}$
$\displaystyle\times(2\pi)^{4}\delta^{4}(\sum_{i}p_{i})\sum_{\rm
spins}|\mathcal{M}|^{2}$ $\displaystyle=\prod_{\text{final
}i}F_{\text{deg,}i}\times n_{b}\langle\sigma v\rangle_{b}\,,$ (50)

where $\langle\cdot\rangle_{b}$ denotes the thermal average taken over the
particle $b$, defined by

$\displaystyle\langle X\rangle_{b}$
$\displaystyle=\frac{\mathfrak{g_{b}}}{n_{b}}\int\frac{d^{3}p_{b}}{(2\pi)^{3}}f_{b}X\,.$
(51)

Furthermore, $\sigma$ denotes the cross section of the process,
$v=\sqrt{(p_{b}\cdot p_{\chi})^{2}-m_{\chi}^{2}m_{b}^{2}}/E_{b}E_{\chi}$ is
the Møller velocity and we have approximated the Pauli-blocking effects by
introducing the degeneracy factors defined in Eq. (32).

By performing the $d^{3}p_{\chi}$ integral in Eq. (6), one arrives at

$\mathcal{C}_{\rm
abs}^{b}=n_{\chi}\langle\Gamma_{\chi}^{b}\rangle_{\chi}=n_{\chi}n_{b}\langle\sigma
v\rangle_{\chi\,b}\,,$ (52)

where $\langle\sigma v\rangle_{\chi\,b}$ denotes the thermal average over the
complete initial state kinematics, i.e.

$\displaystyle\langle\sigma
v\rangle_{\chi\,b}\equiv\frac{\mathfrak{g}_{\chi}\mathfrak{g}_{b}}{n_{\chi}n_{b}}\int\frac{d^{3}p_{\chi}}{(2\pi)^{3}}f_{\chi}\frac{d^{3}p_{b}}{(2\pi)^{3}}f_{b}\sigma(p_{\chi},p_{b})v\,.$
(53)

The inverse process defines an analogous collision operator for production of
$\chi$’s, $\mathcal{C}_{\rm prod}^{b}$, which in conditions of thermal and
chemical equilibrium reads

$\displaystyle\mathcal{C}_{\rm prod}^{b}=\mathcal{C}_{\rm abs}^{b}\,.$ (54)

This detailed balance relation can then be used to estimate the MFP in the PNS

$\displaystyle\langle\lambda(p_{\chi})\rangle_{\chi}=\langle\frac{v_{\chi}}{\Gamma_{\chi}^{b}}\rangle_{\chi}\approx\frac{\langle
v_{\chi}\rangle_{\chi}}{\langle\Gamma_{\chi}^{b}\rangle_{\chi}}=\frac{n_{\chi}}{\mathcal{C}_{\rm
prod}^{b}}\langle v_{\chi}\rangle_{\chi}\,.$ (55)

##  Appendix B Cross sections for annihilation ($2\to 2$)

In this appendix, we list cross sections for the $2\to 2$ processes
$\ell\ell\to\chi\chi,\chi\chi\to\ell\ell$ ($s$-channel) and
$\ell\chi\to\ell\chi$ ($t$-channel), where $\ell$ generically refer to any
lepton, including neutrinos. For the effective interactions in Eq. (47), a
linear independent basis is given by the operators $O_{SS}$, $O_{PP}$,
$O_{SP}$, $O_{PS}$, $O_{VV}$, $O_{AA}$, $O_{AV}$, $O_{VA}$, $O_{TT}$,
$O_{T^{\prime}T}$ and we find131313We disagree with Ref. Guha:2018mli  on the
$t$-channel recovering their results only in the limit $m_{\ell},m_{\chi}\to
0$.

$\displaystyle\begin{split}\sigma_{\ell\ell\to\chi\chi}=&\frac{\sqrt{s-4m_{\chi}^{2}}}{48\pi
s\Lambda_{l}^{4}\sqrt{s-4m_{\ell}^{2}}}\bigg{[}3\,C_{SS}^{2}(s-4m_{\ell}^{2})(s-4m_{\chi}^{2})+3\,C_{PP}^{2}s^{2}+3\,C_{PS}^{2}s(s-4m_{\chi}^{2})+3\,C_{SP}^{2}s(s-4m_{\ell}^{2})\\\
&+4\,C_{VV}^{2}(s+2m_{\ell}^{2})(s+2m_{\chi}^{2})+4\,C_{AA}^{2}(s^{2}-4s(m_{\ell}^{2}+m_{\chi}^{2})+28m_{\ell}^{2}m_{\chi}^{2})+4\,C_{VA}^{2}(s+2m_{\ell}^{2})(s-4m_{\chi}^{2})\\\
&+4\,C_{VA}^{2}(s-4m_{\ell}^{2})(s+2m_{\chi}^{2})+8\,C_{TT}^{2}(s^{2}+2s(m_{\ell}^{2}+m_{\chi}^{2})+40m_{\ell}^{2}m_{\chi}^{2})+8\,C_{T^{\prime}T}^{2}(s^{2}+2s(m_{\ell}^{2}+m_{\chi}^{2})-32m_{\ell}^{2}m_{\chi}^{2})\\\
&-24\,C_{AA}\,C_{PP}\,s\,m_{\ell}m_{\chi}+144\,C_{VV}\,C_{TT}\,s\,m_{\ell}m_{\chi}\bigg{]}\,,\end{split}$
(56)

$\displaystyle\begin{split}\sigma_{\chi\ell\to\chi\ell}=&\frac{1}{48\pi\Lambda_{l}^{4}s^{3}}\bigg{[}C_{SS}^{2}\Big{(}s^{4}+2s^{3}(m_{\ell}^{2}+m_{\chi}^{2})-2s^{2}(3m_{\ell}^{4}-14m_{\ell}^{2}m_{\chi}^{2}+3m_{\chi}^{4})+2s(m_{\ell}^{2}+m_{\chi}^{2})(m_{\ell}^{2}-m_{\chi}^{2})^{2}+(m_{\ell}^{2}-m_{\chi}^{2})^{4}\Big{)}\\\
&+C_{SP}^{2}\Big{(}m_{\ell}^{4}-2m_{\ell}^{2}(m_{\chi}^{2}-2s)+(s-m_{\chi}^{2})^{2}\Big{)}\Big{(}m_{\ell}^{4}-2m_{\ell}^{2}(m_{\chi}^{2}+s)+(s-m_{\chi}^{2})^{2}\Big{)}\\\
&+C_{PS}^{2}\Big{(}m_{\ell}^{4}-2m_{\ell}^{2}(s+m_{\chi}^{2})+(s-m_{\chi}^{2})^{2}\Big{)}\Big{(}m_{\ell}^{4}+m_{\chi}^{4}-2m_{\ell}^{2}(m_{\chi}^{2}+s)+s(s+4m_{\chi}^{2})\Big{)}\\\
&+C_{PP}^{2}\Big{(}m_{\ell}^{2}-2m_{\ell}^{2}(s+m_{\chi}^{2})+(s-m_{\chi}^{2})^{2}\Big{)}^{2}\\\
&+2C_{VV}^{2}\Big{(}4s^{4}-10s^{3}(m_{\ell}^{2}+m_{\chi}^{2})+s^{2}(9m_{\ell}^{2}+22m_{\ell}^{2}m_{\chi}^{2}+9m_{\chi}^{4}\Big{)}-4s(m_{\ell}^{2}+m_{\chi}^{2})(m_{\ell}^{2}-m_{\chi}^{2})^{2}+(m_{\ell}^{2}-m_{\chi}^{2})^{4}\Big{)}\\\
&+2\,C_{VA}^{2}\Big{(}s-(m_{\ell}+m_{\chi})^{2}\Big{)}\Big{(}s-(m_{\ell}-m_{\chi})^{2}\Big{)}\Big{(}m_{\ell}^{4}-2m_{\ell}^{2}(m_{\chi}^{2}+s)+(2s+m_{\chi}^{2})^{2}\Big{)}\\\
&+2\,C_{AV}^{2}\Big{(}s-(m_{\ell}+m_{\chi})^{2}\Big{)}\Big{(}s-(m_{\ell}-m_{\chi})^{2}\Big{)}\Big{(}4s^{2}+2s(2m_{\ell}^{2}-m_{\chi}^{2})+(m_{\ell}^{2}-m_{\chi}^{2})^{2}\Big{)}\\\
&+2\,C_{AA}^{2}\Big{(}4s^{4}-4s^{3}(m_{\ell}^{2}+m_{\chi}^{2})-s^{2}(3m_{\ell}^{4}-46m_{\ell}^{2}m_{\chi}^{2}+3m_{\chi}^{4})+2s(m_{\ell}^{2}+m_{\chi}^{2})(m_{\ell}^{2}-m_{\chi}^{2})^{2}+(m_{\ell}^{2}-m_{\chi}^{2})^{4}\Big{)}\\\
&+8\,C_{TT}^{2}\Big{(}7s^{4}-13s^{3}(m_{\ell}^{2}+m_{\chi}^{2})+2s^{2}(3m_{\ell}^{4}+26m_{\ell}^{2}m_{\chi}^{2}+3m_{\chi}^{4})-s(m_{\ell}^{2}+m_{\chi}^{2})(m_{\ell}^{2}-m_{\chi}^{2})^{2}+(m_{\ell}^{2}-m_{\chi}^{2})^{4}\Big{)}\\\
&+8\,C_{T^{\prime}T}^{2}\Big{(}s-(m_{\ell}+m_{\chi})^{2}\Big{)}\Big{(}s-(m_{\ell}-m_{\chi})^{2}\Big{)}\Big{(}7s^{2}+s(m_{\ell}^{2}+m_{\chi}^{2})+(m_{\ell}^{2}-m_{\chi}^{2})^{2}\Big{)}\\\
&+4\,C_{SS}\Big{(}C_{TT}(m_{\ell}^{8}-m_{\ell}^{6}(4m_{\chi}^{2}+s)+m_{\ell}^{4}m_{\chi}^{2}(6m_{\chi}^{2}+s)+m_{\ell}^{2}(-4m_{\chi}^{6}+m_{\chi}^{4}s-8m_{\chi}^{2}s^{2}+5s^{3})\\\
&+m_{\chi}^{8}-m_{\chi}^{6}s+5m_{\chi}^{2}s^{3}-2s^{4})-3\,C_{VV}m_{\ell}m_{\chi}s(m_{\ell}^{4}-2m_{\ell}^{2}(m_{\chi}^{2}-s)+m_{\chi}^{4}+2m_{\chi}^{2}s-3s^{2})\Big{)}\\\
&+4\,C_{PP}\Big{(}3\,C_{AA}m_{\chi}s(m_{\ell}^{4}-2m_{\ell}^{2}(m_{\chi}^{2}+s)+(m_{\chi}^{2}-s)^{2})+C_{TT}(m_{\ell}^{8}-m_{\ell}^{6}(4m_{\chi}^{2}+s)+m_{\ell}^{4}m_{\chi}^{2}(6m_{\chi}^{2}+s)\\\
&+m_{\ell}^{2}(-4m_{\chi}^{6}+m_{\chi}^{4}s-8m_{\chi}^{2}s^{2}+5s^{3})+m_{\chi}^{8}-m_{\chi}^{6}s+5m_{\chi}^{2}s^{3}-2s^{4})\Big{)}\\\
&+4(C_{SP}+C_{PS})\Big{(}C_{T^{\prime}T}(m_{\ell}^{8}-m_{\ell}^{6}(4m_{\chi}^{2}+s)+m_{\ell}^{4}m_{\chi}^{2}(6m_{\chi}^{2}+s)+m_{\ell}^{2}(-4m_{\chi}^{6}+m_{\chi}^{4}s-8m_{\chi}^{2}s^{2}+5s^{3})+m_{\chi}^{8}\\\
&-m_{\chi}^{6}s+5m_{\chi}^{2}s^{3}-2s^{4})+6C_{VA}m_{\ell}m_{\chi}s^{2}(m_{\chi}^{2}-m_{\ell}^{2})\Big{)}\\\
&-4\,C_{VV}\Big{(}C_{AA}(m_{\ell}^{8}-m_{\ell}^{6}(4m_{\chi}^{2}+s)+m_{\ell}^{4}m_{\chi}^{2}\left(6m_{\chi}^{2}+s\right)+m_{\ell}^{2}(-4m_{\chi}^{6}+m_{\chi}^{4}s-8m_{\chi}^{2}s^{2}+5s^{3})+m_{\chi}^{8}\\\
&-m_{\chi}^{6}s+5m_{\chi}^{2}s^{3}-2s^{4})+18C_{TT}m_{\ell}m_{\chi}s(m_{\ell}^{4}-2m_{\ell}^{2}(m_{\chi}^{2}+s)+(m_{\chi}^{2}-s)^{2})\Big{)}\\\
&+72\,C_{AA}\,C_{TT}m_{\ell}m_{\chi}s(m_{\ell}^{4}-2m_{\ell}^{2}\left(m_{\chi}^{2}-s\right)+m_{\chi}^{4}+2m_{\chi}^{2}s-3s^{2})\\\
&+C_{VA}\Big{(}144C_{T^{\prime}T}m_{\ell}m_{\chi}s^{2}(m_{\ell}^{2}-m_{\chi}^{2})-4C_{AV}(m_{\ell}^{8}-m_{\ell}^{6}(4m_{\chi}^{2}+s)+m_{\ell}^{4}m_{\chi}^{2}(6m_{\chi}^{2}+s)\\\
&+m_{\ell}^{2}(-4m_{\chi}^{6}+m_{\chi}^{4}s-8m_{\chi}^{2}s^{2}+5s^{3})+m_{\chi}^{8}-m_{\chi}^{6}s+5m_{\chi}^{2}s^{3}-2s^{4})\Big{)}\\\
&+144\,C_{AV}\,C_{T^{\prime}T}\,s^{2}m_{\ell}m_{\chi}(m_{\ell}^{2}-m_{\chi}^{2})\bigg{]}\,.\end{split}$
(57)

We can also generalize the cross sections for the $Z^{\prime}$ model in Eq.
(1) by including generic vector ($V_{i}$) and axial ($A_{i}$) couplings to the
leptons (such that $V_{\ell}=1$, $A_{\ell}=0$ gives back the model in Eq.
(1)):

$\displaystyle\sigma_{\ell\ell\to\chi\chi}^{V}=$
$\displaystyle\frac{g_{\ell}^{2}g_{\chi}^{2}}{12\pi
s}\sqrt{\frac{s-4m_{\chi}^{2}}{s-4m_{\ell}^{2}}}\frac{s+2m_{\chi}^{2}}{(s-m_{Z^{\prime}}^{2})^{2}+m_{Z^{\prime}}^{2}\Gamma_{Z^{\prime}}^{2}}$
$\displaystyle\times\Big{[}V_{\ell}^{2}(s+2m_{\ell}^{2})+A_{\ell}^{2}(s-4m_{\ell}^{2})\Big{]}\,.$

The $t$-channel for small $Z^{\prime}$ masses is more involved because the
propagator depends on the Mandelstam variable $t$, over which one integrates
to obtain the total cross section. Neglecting the $Z^{\prime}$ decay width,
the resulting expression reads

$\displaystyle\sigma_{\chi\ell\to\chi\ell}^{V}$
$\displaystyle=\frac{g_{\ell}^{2}g_{\chi}^{2}}{8\pi}\Big{[}(V_{\ell}^{2}+A_{\ell}^{2})\frac{2s+m_{Z^{\prime}}^{2}}{sm_{Z^{\prime}}^{2}}$
(58)
$\displaystyle+\frac{1}{m_{Z^{\prime}}^{2}}\frac{V_{\ell}^{2}(m_{Z^{\prime}}^{4}+8m_{\ell}^{2}m_{\chi}^{2})+A_{\ell}^{2}m_{Z^{\prime}}^{2}(m_{Z^{\prime}}^{2}-4m_{\ell}^{2})}{m_{\ell}^{4}-2m_{\ell}^{2}(s+m_{\chi}^{2})+(s-m_{\chi}^{2})^{2}+sm_{Z^{\prime}}^{2}}$
$\displaystyle-2\frac{V_{\ell}^{2}(s+m_{Z^{\prime}}^{2})+A_{\ell}^{2}(s+m_{Z^{\prime}}^{2}-2m_{\ell}^{2})}{m_{\ell}^{4}-2m_{\ell}^{2}(s+m_{\chi}^{2})+(s-m_{\chi}^{2})^{2}}$
$\displaystyle\times\log\frac{m_{\ell}^{4}-2m_{\ell}^{2}(s+m_{\chi}^{2})+(s-m_{\chi}^{2})^{2}+sm_{Z^{\prime}}^{2}}{sm_{Z^{\prime}}^{2}}\Big{]}\,.$

The cross sections for the inverse process $\chi\chi\to\ell\ell$ can be
obtained using

$\sigma_{\chi\chi\to\ell\ell}=\frac{s-4m_{\ell}^{2}}{s-4m_{\chi}^{2}}\sigma_{\ell\ell\to\chi\chi}\,.$
(59)

##  Appendix C Cross sections for photoproduction ($2\to 3$) in the EFT limit

In this appendix, we derive the cross sections for the photoproduction
processes $\ell^{-}\gamma\to\ell^{-}\bar{\chi}\chi$ with the effective
operators given in Eq. (47). For simplicity, we rewrite Eq. (47) factorizing
the Wilson Coefficients in terms of a leptonic and dark current,
$C_{XY}^{l}=X_{l}Y_{\chi}$. We start with the most simple case of scalar
interactions. The amplitude reads

$\displaystyle i\mathcal{M}=$
$\displaystyle\frac{e}{\Lambda_{\ell}^{2}}\epsilon_{\mu}(p_{b})\bar{u}_{\chi}(p_{1})(S_{\chi}+iP_{\chi}\gamma_{5})v_{\chi}(p_{2})$
$\displaystyle\times\bar{u}_{\ell}(p_{3})\Big{[}(S_{\ell}+iP_{\ell}\gamma_{5})\frac{\not{p}_{a}+\not{p}_{b}+m_{\ell}}{(p_{a}+p_{b})^{2}-m_{\ell}^{2}}\gamma^{\mu}$
$\displaystyle+\gamma^{\mu}\frac{\not{p}_{3}-\not{p}_{b}+m_{\ell}}{(p_{3}-p_{b})^{2}-m_{\ell}^{2}}(S_{\ell}+iP_{\ell}\gamma_{5})\Big{]}u_{\ell}(p_{a})\,.$
(60)

The squared and spin-averaged amplitude factorizes in two contributions

$\displaystyle\overline{|\mathcal{M}|^{2}}=X(p_{1},p_{2})L(p_{a},p_{b},p_{3})\,,$
(61)

where the $X$ and $L$ denote the traces over dark and SM particles,
respectively. The phase space can be factorized

$\displaystyle d\Phi_{3}$ $\displaystyle(p_{a}+p_{b};p_{1},p_{2},p_{3})=$
$\displaystyle\frac{dm_{12}^{2}}{2\pi}d\Phi_{2}(p_{a}+p_{b};p_{12},p_{3})d\Phi_{2}(p_{12};p_{1},p_{2})\,,$
(62)

where we introduced $s=(p_{a}+p_{b})^{2}$ and the momentum of the effective
two-body system of dark particles $p_{12}=p_{1}+p_{2}$ with invariant mass
$m_{12}^{2}=p_{12}^{2}$.

We start with the dark system, i.e. $d\Phi_{2}(p_{12};p_{1},p_{2})$. The
function $X(p_{1},p_{2})$ is a scalar and can therefore only depend on the
scalar product $p_{1}p_{2}$, which can be rewritten in terms of $m_{12}^{2}$
using $p_{1}p_{2}=(m_{12}^{2}-2m_{\chi}^{2})/2$, leading to
$X(p_{1},p_{2})=\tilde{X}(m_{12})$. Thus, we obtain

$\displaystyle\int
d\Phi_{2}(p_{12};p_{1},p_{2})X(p_{1},p_{2})=\frac{\tilde{X}(m_{12})}{8\pi}\sqrt{1-4\frac{m_{\chi}^{2}}{m_{12}^{2}}}\,.$

The second phase-space integral can be simplified as

$\displaystyle\int d\Phi_{2}(p_{a}+p_{b};$ $\displaystyle
p_{12},p_{3})=\frac{1}{(4\pi)^{2}\sqrt{s}}$
$\displaystyle\times\int\bar{p}_{3}d\bar{p}_{3}\delta\Big{(}\bar{p}_{3}-\frac{\sqrt{s}}{2}\beta\Big{)}d\cos\theta
d\phi\,,$ (63)

with

$\displaystyle\beta$
$\displaystyle=\sqrt{1-2\frac{m_{12}^{2}+m_{\ell}^{2}}{s}+\frac{(m_{12}^{2}-m_{\ell}^{2})^{2}}{s^{2}}}\,,$
(64)

where $\bar{p}_{3}$ is the spatial component of the 4-vector $p_{3}$, i.e.
$p_{3}^{2}=E_{3}^{2}-\bar{p}_{3}^{2}=m_{\ell}^{2}$.

The $p_{3}$-dependence in the function $L(p_{a},p_{b},p_{3})$ can be rewritten
in terms of $s,\bar{p}_{3}$ and $\cos\theta$ in the center-of-mass (CM) frame
using

$\displaystyle\begin{split}p_{a}p_{3}&=\frac{\sqrt{s}}{2}\Big{[}E_{3}\Big{(}1+\frac{m_{\ell}^{2}}{s}\Big{)}-\bar{p}_{3}\cos\theta\Big{(}1-\frac{m_{\ell}^{2}}{s}\Big{)}\Big{]}\,,\\\
p_{b}p_{3}&=\frac{\sqrt{s}}{2}\Big{(}1-\frac{m_{\ell}^{2}}{s}\Big{)}(E_{3}+\bar{p}_{3}\cos\theta)\,.\end{split}$
(65)

One can now evaluate the $d\Phi_{2}(p_{a}+p_{b};p_{12},p_{3})$ integral. After
the trivial integrals in $d\phi$ and $d\bar{p}_{3}$ one can perform the
$d\cos\theta$ integration analytically. At this point, we are only left with
the $m_{12}^{2}$ integration, which can not be done analytically in the
general case. Introducing the dimensionless variable

$\displaystyle\beta_{\chi}=\sqrt{1-4m_{\chi}^{2}/m_{12}^{2}}\,,$ (66)

and
$x=m_{12}^{2}/m_{\ell}^{2},x_{\chi}=m_{\chi}^{2}/m_{\ell}^{2},\hat{s}=s/m_{\ell}^{2}$,
we arrive at the result

$\displaystyle\begin{split}\sigma=&\frac{\alpha
m_{\ell}^{2}}{64\pi^{2}\Lambda_{\ell}^{4}}\frac{1}{(\hat{s}-1)^{3}\hat{s}^{2}}\int_{4x_{\chi}}^{(\sqrt{\hat{s}}-1)^{2}}dx\,x\,\beta_{\chi}\\\
&\times\Big{(}\beta_{\chi}^{2}S_{\chi}^{2}+P_{\chi}^{2}\Big{)}\Big{[}S_{\ell}^{2}f_{S}(x,\hat{s})+P_{\ell}^{2}f_{P}(x,\hat{s})\Big{]}\,,\end{split}$
(67)

where the functions $f_{S},f_{P}$ are defined by

$\displaystyle f_{S}=$ $\displaystyle
4\hat{s}^{2}\Big{(}\hat{s}^{2}-2\hat{s}(x-3)+2x^{2}-10x+9\Big{)}T(x,\hat{s}),$
$\displaystyle-\Big{(}3\hat{s}^{3}+\hat{s}^{2}(25-7x)+\hat{s}(5-2x)+x-1\Big{)}R(x,\hat{s})\,,$
$\displaystyle f_{P}=$ $\displaystyle
4\hat{s}^{2}\Big{(}\hat{s}^{2}-2\hat{s}(x+1)+2x^{2}-2x+1\Big{)}T(x,\hat{s}),$
$\displaystyle-\Big{(}3\hat{s}^{3}-7\hat{s}^{2}(x+1)+\hat{s}(5-2x)+x-1\Big{)}R(x,\hat{s})\,,$
(68)

with,

$\displaystyle R(x,\hat{s})=$
$\displaystyle\sqrt{\hat{s}^{2}-2\hat{s}(x+1)+(x-1)^{2}}\,,$ $\displaystyle
T(x,\hat{s})=$ $\displaystyle\tanh^{-1}\frac{R(x,\hat{s})}{\hat{s}-x+1}\,.$
(69)

In the limit $m_{\chi}\to 0$, the final integral can be evaluated
analytically,

$\displaystyle\sigma$ $\displaystyle=\frac{\alpha
m_{\ell}^{2}}{2304\pi^{2}\Lambda_{\ell}^{4}}\frac{1}{(\hat{s}-1)^{3}\hat{s}^{2}}$
$\displaystyle\times\Big{(}S_{\chi}^{2}+P_{\chi}^{2}\Big{)}\Big{(}S_{\ell}^{2}g_{S}(\hat{s})+P_{\ell}^{2}g_{P}(\hat{s})\Big{)}\,,$
(70)

where

$\displaystyle\begin{split}g_{S}=&-79\hat{s}^{6}-14\hat{s}^{5}-189\hat{s}^{4}-296\hat{s}^{3}+527\hat{s}^{2}+54\hat{s}-3\\\
&+12\hat{s}^{2}\Big{(}2\hat{s}^{4}+10\hat{s}^{3}+9\hat{s}^{2}+44\hat{s}+25\Big{)}\log\hat{s}\,,\\\
g_{P}=&-79\hat{s}^{6}+338\hat{s}^{5}+675\hat{s}^{4}-1160\hat{s}^{3}+175\hat{s}^{2}+54\hat{s}\\\
&-3+12\hat{s}^{2}\Big{(}2\hat{s}^{4}+2\hat{s}^{3}-63\hat{s}^{2}-28\hat{s}+17\Big{)}\log\hat{s}\,.\end{split}$
(71)

The matrix element for vector-like interactions reads

$\displaystyle i\mathcal{M}=$
$\displaystyle\frac{e}{\Lambda_{\ell}^{2}}\epsilon_{\mu}(p_{b})\bar{u}_{\chi}(p_{1})\gamma_{\nu}(V_{\chi}+iA_{\chi}\gamma_{5})v_{\chi}(p_{2})$
(72)
$\displaystyle\times\bar{u}_{\ell}(p_{3})\Big{[}\gamma^{\nu}(V_{\ell}+iA_{\ell}\gamma_{5})\frac{\not{p}_{a}+\not{p}_{b}+m_{\ell}}{(p_{a}+p_{b})^{2}-m_{\ell}^{2}}\gamma^{\mu}$
$\displaystyle+\gamma^{\mu}\frac{\not{p}_{3}-\not{p}_{b}+m_{\ell}}{(p_{3}-p_{b})^{2}-m_{\ell}^{2}}\gamma^{\nu}(V_{\ell}+iA_{\ell}\gamma_{5})\Big{]}u_{\ell}(p_{a})\,.$

After squaring and spin-averaging, one again obtains two traces for the SM and
dark part of the amplitude, respectively. Due to the vector-like nature of the
interaction, these traces are contracted with two Lorentz indices, one arising
from $\mathcal{M}$ and one from $\mathcal{M}^{\dagger}$

$\displaystyle\overline{|\mathcal{M}|^{2}}=X_{\mu\nu}(p_{1},p_{2})L^{\mu\nu}(p_{a},p_{b},p_{3})\,.$
(73)

Due to the fact that $X_{\mu\nu}$ can not only depend on $m_{12}^{2}$, but
also on $p_{1}^{\mu},p_{2}^{\mu}$, we can not simply factor $X_{\mu\nu}$ out
of the $d\Phi(p_{12};p_{1},p_{2})$ integral as in the scalar case. However,
Lorentz covariance implies that the integral can only depend on the vector
$p_{12}^{\mu}=p_{1}^{\mu}+p_{2}^{\mu}$

$\displaystyle\int
d\Phi_{2}(p_{12};p_{1},p_{2})X^{\mu\nu}(p_{1},p_{2})=A_{1}m_{12}^{2}g^{\mu\nu}+A_{2}p_{12}^{\mu}p_{12}^{\nu}\,,$
(74)

where $A_{1},A_{2}$ are functions of $m_{12}^{2}$

$\displaystyle A_{1}$
$\displaystyle=-\Big{(}(V_{\chi}^{2}+A_{\chi}^{2})+2\frac{m_{\chi}^{2}}{m_{12}^{2}}(V_{\chi}^{2}-2A_{\chi}^{2})\Big{)}\,,$
$\displaystyle A_{2}$
$\displaystyle=\left(1+\frac{2m_{\chi}^{2}}{m_{12}^{2}}\right)(V_{\chi}^{2}+A_{\chi}^{2})\,.$
(75)

From this point on, the calculation is equivalent to the scalar case. We
arrive at the result

$\displaystyle\sigma=$ $\displaystyle\frac{\alpha
m_{\ell}^{2}}{48\pi^{2}\Lambda_{\ell}^{4}}\frac{1}{(\hat{s}-1)^{3}\hat{s}^{2}}$
(76)
$\displaystyle\times\int_{4x_{\chi}}^{(\sqrt{\hat{s}}-1)^{2}}dx\,\beta_{\chi}\Big{(}V_{\ell}^{2}f_{V}(x,x_{\chi},\hat{s})+A_{\ell}^{2}f_{A}(x,x_{\chi},\hat{s})\Big{)}\,,$

with

$\displaystyle f_{V}$
$\displaystyle=-4A_{1}\hat{s}^{2}x\Big{(}\hat{s}^{2}-2\hat{s}(x+3)+2x^{2}+2x-3\Big{)}T(x,\hat{s})$
$\displaystyle-
A_{1}x\Big{(}\hat{s}^{3}+\hat{s}^{2}(7x+15)+\hat{s}(2x-1)-x+1\Big{)}R(x,\hat{s})\,,$
$\displaystyle f_{A}$
$\displaystyle=-4A_{1}\hat{s}^{2}x\Big{(}\hat{s}^{2}-2\hat{s}(x-5)+2x^{2}-14x+13\Big{)}T(x,\hat{s})$
$\displaystyle-
A_{1}x\Big{(}\hat{s}^{3}+7\hat{s}^{2}(x-7)+\hat{s}(2x-1)-x+1\Big{)}R(x,\hat{s})$
$\displaystyle+8A_{2}\hat{s}^{2}\Big{(}\hat{s}^{2}-2\hat{s}(x+1)+2x^{2}-2x+1\Big{)}T(x,\hat{s})$
$\displaystyle-2A_{2}\Big{(}3\hat{s}^{3}-7\hat{s}^{2}(1+x)+\hat{s}(5-2x)+x-1\Big{)}R(x,\hat{s})\,.$
(77)

In the limit $m_{\chi}\to 0$, one can again perform the $dx$ integration
analytically and arrive at

$\displaystyle\sigma$ $\displaystyle=\frac{\alpha
m_{\ell}^{2}}{1728\pi^{2}\Lambda_{\ell}^{4}}\frac{1}{\hat{s}^{2}(\hat{s}-1)^{3}}$
$\displaystyle\times(V_{\chi}^{2}+A_{\chi}^{2})\Big{(}V_{\ell}^{2}g_{V}(\hat{s})+A_{\ell}^{2}g_{A}(\hat{s})\Big{)}\,,$
(78)

with

$\displaystyle g_{V}=$
$\displaystyle-55\hat{s}^{6}+682\hat{s}^{5}+483\hat{s}^{4}-968\hat{s}^{3}-169\hat{s}^{2}+30\hat{s}$
$\displaystyle-3+12\hat{s}^{2}\Big{(}2\hat{s}^{4}-14\hat{s}^{3}-87\hat{s}^{2}-52\hat{s}+1\Big{)}\log\hat{s}\,,$
$\displaystyle g_{A}=$
$\displaystyle-55\hat{s}^{6}-254\hat{s}^{5}+219\hat{s}^{4}-296\hat{s}^{3}+119\hat{s}^{2}+294\hat{s}$
$\displaystyle-27+12\hat{s}^{2}\Big{(}2\hat{s}^{4}+10\hat{s}^{3}+33\hat{s}^{2}-4\hat{s}+49\Big{)}\log\hat{s}\,.$

Finally, the matrix element for tensor-like interactions reads

$\displaystyle i\mathcal{M}=$
$\displaystyle\frac{e}{\Lambda_{\ell}^{2}}\epsilon_{\mu}(p_{b})\bar{u}_{\chi}(p_{1})\sigma_{\nu\rho}(T_{\chi}+T^{\prime}_{\chi}\gamma_{5})v_{\chi}(p_{2})$
(79)
$\displaystyle\times\bar{u}_{\ell}(p_{3})\Big{[}\sigma^{\nu\rho}\frac{\not{p}_{a}+\not{p}_{b}+m_{\ell}}{(p_{a}+p_{b})^{2}-m_{\ell}^{2}}\gamma^{\mu}$
(80)
$\displaystyle+\gamma^{\mu}\frac{\not{p}_{3}-\not{p}_{b}+m_{\ell}}{(p_{3}-p_{b})^{2}-m_{\ell}^{2}}\sigma^{\nu\rho}\Big{]}u_{\ell}(p_{a})\,.$
(81)

The squared and spin-averaged matrix element can be factorized as

$\displaystyle\overline{|\mathcal{M}|^{2}}=X_{\mu\nu\rho\sigma}(p_{1},p_{2})L^{\mu\nu\rho\sigma}(p_{a},p_{b},p_{3})\,.$
(82)

We define the index ordering in $X^{\mu\nu\rho\sigma}$ such that the first and
second pair of indices each corresponds to one $\sigma^{\mu\nu}$ in the trace
that constitutes $X^{\mu\nu\rho\sigma}$, i.e. $X^{\mu\nu\rho\sigma}$ is
antisymmetric under $\mu\leftrightarrow\nu$ and $\rho\leftrightarrow\sigma$.
After proper antisymmetrization, there are only two Lorentz structures that
satisfy this condition, leading to

$\displaystyle\int$ $\displaystyle
d\Phi_{2}(p_{12};p_{1},p_{2})X^{\mu\nu\rho\sigma}(p_{1},p_{2})=m_{12}^{2}B_{1}\Big{(}g^{\mu\rho}g^{\nu\sigma}-g^{\nu\rho}g^{\mu\sigma}\Big{)}$
$\displaystyle+B_{2}\Big{(}p_{12}^{\mu}(p_{12}^{\rho}g^{\nu\sigma}-p_{12}^{\sigma}g^{\nu\rho})-p_{12}^{\nu}(p_{12}^{\rho}g^{\mu\sigma}-p_{12}^{\sigma}g^{\mu\rho})\Big{)}\,,$

where $B_{1},B_{2}$ are functions of $m_{12}^{2}$

$\displaystyle B_{1}=$
$\displaystyle\frac{1}{12\pi}\beta_{\chi}\Big{(}(T_{\chi}^{2}+T^{\prime
2}_{\chi})+4\frac{m_{\chi}^{2}}{m_{12}^{2}}(2T^{\prime
2}_{\chi}-T_{\chi}^{2})\Big{)}\,,$ $\displaystyle B_{2}=$
$\displaystyle-\frac{1}{6\pi}\beta_{\chi}\frac{m_{12}^{2}+2m_{\chi}^{2}}{m_{12}^{2}}(T_{\chi}^{2}+T^{\prime
2}_{\chi})\,.$

After following the standard procedure as outlined for the scalar case, one
arrives at the result

$\displaystyle\sigma=$ $\displaystyle\frac{\alpha
m_{\ell}^{2}}{4\pi\Lambda_{\ell}^{4}}\frac{1}{(\hat{s}-1)^{3}\hat{s}^{2}}$
(83)
$\displaystyle\int_{4x_{\chi}}^{(\sqrt{\hat{s}}-1)^{2}}dx\Big{(}B_{1}f_{1}(x,x_{\chi},\hat{s})+B_{2}f_{2}(x,x_{\chi},\hat{s})\Big{)}\,,$

where

$\displaystyle f_{1}=$ $\displaystyle
96\hat{s}^{2}\Big{(}(\hat{s}-x+1)T(x,\hat{s})-R(x,\hat{s})\Big{)}\,,$
$\displaystyle f_{2}=$
$\displaystyle-\Big{(}4\hat{s}^{4}+\hat{s}^{3}(x-16)+\hat{s}^{2}(7x^{2}+59x+24)$
$\displaystyle+\hat{s}(2x^{2}+7x-16)-x^{2}-3x+4\Big{)}R(x,\hat{s})$
$\displaystyle-4\hat{s}^{2}x\Big{(}\hat{s}^{2}-2\hat{s}(x+9)+2x^{2}+14x-15\Big{)}T(x,\hat{s})\,.$
(84)

In the limit $m_{\chi}\to 0$, the final $dx$ integration can be done
analytically and we arrive at

$\displaystyle\sigma=$ $\displaystyle\frac{\alpha
m_{\ell}^{2}}{864\pi^{2}\Lambda_{\ell}^{4}}\frac{1}{(\hat{s}-1)^{3}\hat{s}^{2}}(T_{\chi}^{2}+T^{\prime
2}_{\chi})$ (85)
$\displaystyle\times\Big{[}17\hat{s}^{6}+370\hat{s}^{5}+675\hat{s}^{4}-344\hat{s}^{3}-1153\hat{s}^{2}+486\hat{s}-51$
$\displaystyle+12\hat{s}^{2}\Big{(}2\hat{s}^{4}-26\hat{s}^{3}-27\hat{s}^{2}-136\hat{s}+37\Big{)}\log\hat{s}\Big{]}\,.$

##  Appendix D Full calculation of photoproduction

The energy-loss rate per unit volume for photoproduction
$\ell^{-}(p_{a})\gamma(p_{b})\to\chi(p_{1})\overline{\chi}(p_{2})\ell^{-}(p_{3})$
can be written as

$\displaystyle Q$
$\displaystyle=\frac{1}{32\pi^{4}}\int_{m_{\ell}}^{\infty}dE_{a}\bar{p}_{a}f_{a}\int_{0}^{\infty}dE_{b}E_{b}f_{b}$
$\displaystyle\times\int_{-1}^{1}dc_{\theta}J_{s}(s,E_{a},E_{b})\,,$ (86)
$\displaystyle J_{s}$ $\displaystyle=\int
d\Phi_{3}(p_{a}+p_{b};p_{1},p_{2},p_{3})|{\cal M}|^{2}$
$\displaystyle\times(E_{a}+E_{b}-E_{3})(1-f_{3})\,,$ (87)

where $d\Phi_{3}$ is the Lorentz-invariant 3-body phase space volume,
$\bar{p}_{a}\equiv|\vec{p}_{a}|=\sqrt{E_{a}^{2}-m_{\ell}^{2}}$,
$f_{a}^{-1}=e^{(E_{a}-\mu_{\ell})/T}+1$, $f_{b}^{-1}=e^{E_{b}/T}-1$,
$f_{3}^{-1}=e^{(E_{3}-\mu_{\ell})/T}+1$,
$s=m_{\ell}^{2}+2E_{b}(E_{a}-\bar{p}_{a}c_{\theta})$,
$c_{\theta}\equiv\cos\theta$ and all energies refer to the PNS frame.

Using the decomposition of the 3-body phase space in Eq. (62), the 2-body
phase space integral over the dark fermion momenta $p_{1}$ and $p_{2}$ can be
easily performed, as the spin-summed, squared matrix element $|{\cal M}|^{2}$
factorizes into two contributions, cf. Eq. (73)

$\displaystyle M(p_{a},p_{b},p_{3},m_{12})$ $\displaystyle\equiv\int
d\Phi_{2}(p_{12};p_{1},p_{2})\left|{\cal M}\right|^{2}\,.$ (88)

The remaining phase space integral $d\Phi_{2}(p_{a}+p_{b};p_{12},p_{3})$ is
suitably evaluated in the CM frame, giving

$\displaystyle\int d\Phi_{2}$
$\displaystyle=\frac{2\pi}{\sqrt{s}}\int\frac{\bar{p}_{3}^{\prime}d\bar{p}_{3}^{\prime}dc_{\delta}d\phi}{4(2\pi)^{3}}\delta(\bar{p}_{3}^{\prime}-\beta\sqrt{s}/2)\,,$
(89)

where $c_{\delta}\equiv\cos\delta$,
$\bar{p}_{3}^{\prime}=\sqrt{(E_{3}^{\prime})^{2}-m_{\ell}^{2}}$ with
$E_{3}^{\prime}$ the final lepton energy in the CM frame, $\phi$ and $\delta$
denote the polar and azimuthal angle of ${\vec{p}_{3}}^{\,\prime}$ and
$\beta(s,m_{12})$ is given in Eq. (64).

To evaluate the Lorentz-invariant quantity $M(p_{a},p_{b},p_{3},m_{12}^{2})$
in this frame, we express the relevant Lorentz scalars in the CM frame, where
$p_{a}^{\prime}=(E_{a}^{\prime},0,0,E_{b}^{\prime})$,
$p_{b}^{\prime}=(E_{b}^{\prime},0,0,-E_{b}^{\prime})$

$\displaystyle p_{a}p_{b}$ $\displaystyle=\frac{s-m_{\ell}^{2}}{2}\,,$
$\displaystyle p_{a}p_{3}$
$\displaystyle=E_{a}^{\prime}E_{3}^{\prime}+E_{b}^{\prime}\bar{p}_{3}^{\prime}c_{\delta}\,,$
$\displaystyle p_{b}p_{3}$
$\displaystyle=E_{b}^{\prime}E_{3}^{\prime}-E_{b}^{\prime}\bar{p}_{3}^{\prime}c_{\delta}\,,$
(90)

where the energies in the CM frame can be written in terms of $s$ only using

$\displaystyle E_{a}^{\prime}$
$\displaystyle=\frac{s+m_{\ell}^{2}}{2\sqrt{s}}\,,$ $\displaystyle
E_{b}^{\prime}$ $\displaystyle=\frac{s-m_{\ell}^{2}}{2\sqrt{s}}\,.$ (91)

These expressions allow to write $M(p_{a},p_{b},p_{3},m_{12})$ as a function
of $s,\cos\delta,E_{3}^{\prime}$ and $m_{12}$, and it remains to evaluate
$(E_{a}+E_{b}-E_{3})(1-f_{3})$ in the CM frame. $E_{3}$ is related to
$E_{3}^{\prime}$ by a Lorentz boost with velocity
$(\vec{p}_{a}+\vec{p}_{b})/(E_{a}+E_{b})\equiv\vec{p}_{ab}/E_{ab}=\bar{p}_{ab}/E_{ab}(0,s_{\eta},c_{\eta})$
with $s_{\eta}\equiv\sin\eta$, $c_{\eta}\equiv\cos\eta$ and boost factor
$\gamma=E_{ab}/\sqrt{s}$ (since $s=E_{ab}^{2}-\bar{p}_{ab}^{2}$). The same
boost relates $E_{b}$ to $E_{b}^{\prime}$, which can be used to obtain an
expression for $\cos\eta$ as a function of $s$ and the initial energies in the
PNS frame

$\displaystyle c_{\eta}$
$\displaystyle=\frac{s(E_{b}-E_{a})+m_{\ell}^{2}(E_{b}+E_{a})}{(s-m_{\ell}^{2})\sqrt{(E_{a}+E_{b})^{2}-s}}\,.$
(92)

This angle determines the desired relation between $E_{3}$ and
$E_{3}^{\prime}$ as

$\displaystyle E_{3}$
$\displaystyle=\frac{1}{\sqrt{s}}\left(E_{ab}E_{3}^{\prime}+\bar{p}_{ab}\bar{p}_{3}^{\,\prime}(s_{\eta}s_{\phi}s_{\delta}+c_{\eta}c_{\delta})\right)\,.$
(93)

Having expressed all final state energies in the CM frame, one can finally
perform the $d\bar{p}_{3}^{\prime}$ using the $\delta$-function in Eq. (89),
giving

$\displaystyle\bar{p}_{3}^{\prime}$ $\displaystyle=\beta\frac{\sqrt{s}}{2}\,,$
$\displaystyle E_{3}^{\prime}$
$\displaystyle=\frac{\sqrt{s}}{2}\left(1-\frac{m_{12}^{2}-m_{\ell}^{2}}{s}\right)\,.$
(94)

Putting everything together, we finally obtain for Eq. (87)

$\displaystyle J_{s}$
$\displaystyle=\frac{1}{64\pi^{3}}\int_{4m_{\chi}^{2}}^{(\sqrt{s}-m_{l})^{2}}dm^{2}_{12}\,\beta(s,m_{12})\int_{-1}^{1}dc_{\delta}M(s,m_{12}^{2},c_{\delta})$
$\displaystyle\times\int_{0}^{2\pi}d\phi(E_{a}+E_{b}-E_{3})(1-f_{3})\,,$ (95)

where $\beta(s,m_{12})$ is given in Eq. (64), $M(s,m_{12}^{2},c_{\delta})$ is
obtained from Eqs. (88),(90),(91),(94), and
$E_{3}=E_{3}(s,m_{12},E_{a},E_{b},\phi,\delta)$ from Eqs. (92),(93),(94). When
$1-f_{3}$ is to good approximation independent of $E_{3}$, the $d\phi$
integration can be trivially carried out, giving an overall factor of $2\pi$
and $\left(s_{\eta}s_{\phi}s_{\delta}+c_{\eta}c_{\delta}\right)\to
c_{\eta}c_{\delta}$ in Eq (93).

##  Appendix E Numerical resolution of the kinematic integrals for the
annihilation rates

Throughout this work, we use Monte Carlo techniques Lepage:2020tgj  to
evaluate the thermal phase space integrals numerically, e.g. Eq. (10) for
annihilation and Eq. (31) for trapping. For processes where the $Z^{\prime}$
can be produced resonantly, the Breit-Wigner propagator yields a peak in the
Mandelstam variable $s$ located at $m_{Z^{\prime}}$ with width
$\sqrt{m_{Z^{\prime}}\Gamma_{Z^{\prime}}}\sim gm_{Z^{\prime}}$. We want to
calculate these integrals for very small couplings $g\sim 10^{-10}$,
translating into very narrow peaks. When using only several thousands of
points for the Monte Carlo estimate, it is unlikely that this estimate can
resolve the peak structure. In particular, iterative algorithms like VEGAS are
unable to adapt to the peak if they can not extract information about the peak
from samples in early iterations.

Fortunately, we have knowledge about the shape of the peak, allowing us to
provide this information to the Monte Carlo sampler. To achieve this, we
include the Mandelstam variable $s$, corresponding to the peak location, as
one of the integration variables. Subsequently, we divide the integral into
three regions. One of these regions includes only the peak

$m_{Z^{\prime}}^{2}-\sigma m_{Z^{\prime}}\Gamma_{Z^{\prime}}\leq s\leq
m_{Z^{\prime}}^{2}+\sigma m_{Z^{\prime}}\Gamma_{Z^{\prime}}\,,$ (96)

with $\sigma\sim 5$. The other two regions cover the regions to the left and
to the right of the peak, respectively. This approach ensures that the sum of
the three integrals effectively captures and resolves the peak.

These integrals scale as $\propto g^{2}$ when the contribution from the
resonance dominates the integral, otherwise they scale as $\propto g^{4}$. We
can employ this observation to validate our approach of calculating the
integrals, as shown in Fig. 6.

|
---|---
Figure 6: Free-streaming luminosity at R=9.2 km for a $Z^{\prime}$ coupled to
electrons (left panel) and muons (right panel), rescaled by a factor $g^{4}$
to highlight the regions where the resonant contribution dominates.

## References

  * (1) Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * (2) J. Silk et al., _Particle Dark Matter: Observations, Models and Searches_ , Cambridge Univ. Press, Cambridge (2010), 10.1017/CBO9780511770739.
  * (3) F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279.
  * (4) S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * (5) M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, _Can Confinement Ensure Natural CP Invariance of Strong Interactions?_ , _Nucl. Phys. B_ 166 (1980) 493.
  * (6) A.R. Zhitnitsky, _On Possible Suppression of the Axion Hadron Interactions. (In Russian)_ , _Sov. J. Nucl. Phys._ 31 (1980) 260.
  * (7) M. Dine, W. Fischler and M. Srednicki, _A Simple Solution to the Strong CP Problem with a Harmless Axion_ , _Phys. Lett. B_ 104 (1981) 199.
  * (8) B. Holdom, _Two U(1)’s and Epsilon Charge Shifts_ , _Phys. Lett. B_ 166 (1986) 196.
  * (9) S. Dodelson and L.M. Widrow, _Sterile-neutrinos as dark matter_ , _Phys. Rev. Lett._ 72 (1994) 17 [hep-ph/9303287].
  * (10) C. Boehm and P. Fayet, _Scalar dark matter candidates_ , _Nucl. Phys. B_ 683 (2004) 219 [hep-ph/0305261].
  * (11) T. Asaka, S. Blanchet and M. Shaposhnikov, _The nuMSM, dark matter and neutrino masses_ , _Phys. Lett. B_ 631 (2005) 151 [hep-ph/0503065].
  * (12) T. Asaka and M. Shaposhnikov, _The $\nu$MSM, dark matter and baryon asymmetry of the universe_, _Phys. Lett. B_ 620 (2005) 17 [hep-ph/0505013].
  * (13) D. O’Connell, M.J. Ramsey-Musolf and M.B. Wise, _Minimal Extension of the Standard Model Scalar Sector_ , _Phys. Rev. D_ 75 (2007) 037701 [hep-ph/0611014].
  * (14) B. Patt and F. Wilczek, _Higgs-field portal into hidden sectors_ , hep-ph/0605188.
  * (15) M. Pospelov, A. Ritz and M.B. Voloshin, _Secluded WIMP Dark Matter_ , _Phys. Lett. B_ 662 (2008) 53 [0711.4866].
  * (16) D. Gorbunov and M. Shaposhnikov, _How to find neutral leptons of the $\nu$MSM?_, _JHEP_ 10 (2007) 015 [0705.1729].
  * (17) B. Batell, M. Pospelov and A. Ritz, _Exploring Portals to a Hidden Sector Through Fixed Targets_ , _Phys. Rev. D_ 80 (2009) 095024 [0906.5614].
  * (18) G. Krnjaic, _Probing Light Thermal Dark-Matter With a Higgs Portal Mediator_ , _Phys. Rev. D_ 94 (2016) 073009 [1512.04119].
  * (19) C. Arina, J. Hajer and P. Klose, _Portal Effective Theories. A framework for the model independent description of light hidden sector interactions_ , _JHEP_ 09 (2021) 063 [2105.06477].
  * (20) C.A. Manzari and S. Profumo, _Flavor inspired model for dark matter_ , _Phys. Rev. D_ 106 (2022) 075025 [2206.06768].
  * (21) S. Alekhin et al., _A facility to Search for Hidden Particles at the CERN SPS: the SHiP physics case_ , _Rept. Prog. Phys._ 79 (2016) 124201 [1504.04855].
  * (22) M. Battaglieri et al., _US Cosmic Visions: New Ideas in Dark Matter 2017: Community Report_ , in _U.S. Cosmic Visions: New Ideas in Dark Matter_ , 7, 2017 [1707.04591].
  * (23) J. Beacham et al., _Physics Beyond Colliders at CERN: Beyond the Standard Model Working Group Report_ , _J. Phys. G_ 47 (2020) 010501 [1901.09966].
  * (24) G. Lanfranchi, M. Pospelov and P. Schuster, _The Search for Feebly-Interacting Particles_ , 2011.02157.
  * (25) G. Krnjaic et al., _A Snowmass Whitepaper: Dark Matter Production at Intensity-Frontier Experiments_ , 2207.00597.
  * (26) K. Sato and H. Sato, _Higgs Meson Emission from a Star and a Constraint on Its Mass_ , _Prog. Theor. Phys._ 54 (1975) 1564.
  * (27) D.A. Dicus, E.W. Kolb, V.L. Teplitz and R.V. Wagoner, _Astrophysical Bounds on the Masses of Axions and Higgs Particles_ , _Phys. Rev. D_ 18 (1978) 1829.
  * (28) M. Vysotsky, Y. Zeldovich, M. Khlopov and V. Chechetkin, _Some Astrophysical Limitations on Axion Mass_ , _Pisma Zh. Eksp. Teor. Fiz._ 27 (1978) 533.
  * (29) D.A. Dicus, E.W. Kolb, V.L. Teplitz and R.V. Wagoner, _Astrophysical Bounds on Very Low Mass Axions_ , _Phys. Rev. D_ 22 (1980) 839.
  * (30) G.G. Raffelt, _Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles_ (5, 1996).
  * (31) S. Balaji, P.S.B. Dev, J. Silk and Y. Zhang, _Improved stellar limits on a light CP-even scalar_ , 2205.01669.
  * (32) A. Lella, P. Carenza, G. Lucente, M. Giannotti and A. Mirizzi, _Proto-neutron stars as cosmic factories for massive axion-like-particles_ , 2211.13760.
  * (33) M.S. Turner, _Dirac neutrinos and SN1987A_ , _Phys. Rev. D_ 45 (1992) 1066.
  * (34) G. Raffelt and D. Seckel, _A selfconsistent approach to neutral current processes in supernova cores_ , _Phys. Rev. D_ 52 (1995) 1780 [astro-ph/9312019].
  * (35) A. Caputo, G. Raffelt and E. Vitagliano, _Muonic boson limits: Supernova redux_ , _Phys. Rev. D_ 105 (2022) 035022 [2109.03244].
  * (36) A. Caputo, G. Raffelt and E. Vitagliano, _Radiative transfer in stars by feebly interacting bosons_ , _JCAP_ 08 (2022) 045 [2204.11862].
  * (37) G. Raffelt and D. Seckel, _Bounds on Exotic Particle Interactions from SN 1987a_ , _Phys. Rev. Lett._ 60 (1988) 1793.
  * (38) M.S. Turner, _Axions from SN 1987a_ , _Phys. Rev. Lett._ 60 (1988) 1797.
  * (39) R. Mayle, J.R. Wilson, J.R. Ellis, K.A. Olive, D.N. Schramm and G. Steigman, _Constraints on Axions from SN 1987a_ , _Phys. Lett. B_ 203 (1988) 188.
  * (40) A. Burrows, M.S. Turner and R. Brinkmann, _Axions and SN 1987a_ , _Phys. Rev. D_ 39 (1989) 1020.
  * (41) A. Burrows, M. Ressell and M.S. Turner, _Axions and SN1987A: Axion trapping_ , _Phys. Rev. D_ 42 (1990) 3297.
  * (42) J.B. Dent, F. Ferrer and L.M. Krauss, _Constraints on Light Hidden Sector Gauge Bosons from Supernova Cooling_ , 1201.2683.
  * (43) E. Rrapaj and S. Reddy, _Nucleon-nucleon bremsstrahlung of dark gauge bosons and revised supernova constraints_ , _Phys. Rev. C_ 94 (2016) 045805 [1511.09136].
  * (44) J.H. Chang, R. Essig and S.D. McDermott, _Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle_ , _JHEP_ 09 (2018) 051 [1803.00993].
  * (45) P. Carenza, T. Fischer, M. Giannotti, G. Guo, G. Martínez-Pinedo and A. Mirizzi, _Improved axion emissivity from a supernova via nucleon-nucleon bremsstrahlung_ , _JCAP_ 10 (2019) 016 [1906.11844].
  * (46) A. Sung, G. Guo and M.-R. Wu, _Supernova Constraint on Self-Interacting Dark Sector Particles_ , _Phys. Rev. D_ 103 (2021) 103005 [2102.04601].
  * (47) F. Calore, P. Carenza, M. Giannotti, J. Jaeckel, G. Lucente and A. Mirizzi, _Supernova bounds on axionlike particles coupled with nucleons and electrons_ , _Phys. Rev. D_ 104 (2021) 043016 [2107.02186].
  * (48) K. Choi, H.J. Kim, H. Seong and C.S. Shin, _Axion emission from supernova with axion-pion-nucleon contact interaction_ , _JHEP_ 02 (2022) 143 [2110.01972].
  * (49) C.S. Shin and S. Yun, _Dark gauge boson production from neutron stars via nucleon-nucleon bremsstrahlung_ , _JHEP_ 02 (2022) 133 [2110.03362].
  * (50) F. Calore, P. Carenza, M. Giannotti, J. Jaeckel, G. Lucente, L. Mastrototaro et al., _511 keV line constraints on feebly interacting particles from supernovae_ , _Phys. Rev. D_ 105 (2022) 063026 [2112.08382].
  * (51) G. Lucente, L. Mastrototaro, P. Carenza, L. Di Luzio, M. Giannotti and A. Mirizzi, _Axion signatures from supernova explosions through the nucleon electric-dipole portal_ , _Phys. Rev. D_ 105 (2022) 123020 [2203.15812].
  * (52) W. Keil, H.-T. Janka, D.N. Schramm, G. Sigl, M.S. Turner and J.R. Ellis, _A Fresh look at axions and SN-1987A_ , _Phys. Rev. D_ 56 (1997) 2419 [astro-ph/9612222].
  * (53) P. Carenza, B. Fore, M. Giannotti, A. Mirizzi and S. Reddy, _Enhanced Supernova Axion Emission and its Implications_ , _Phys. Rev. Lett._ 126 (2021) 071102 [2010.02943].
  * (54) T. Fischer, P. Carenza, B. Fore, M. Giannotti, A. Mirizzi and S. Reddy, _Observable signatures of enhanced axion emission from protoneutron stars_ , _Phys. Rev. D_ 104 (2021) 103012 [2108.13726].
  * (55) C.S. Shin and S. Yun, _Dark gauge boson emission from supernova pions_ , 2211.15677.
  * (56) J. Martin Camalich, M. Pospelov, P.N.H. Vuong, R. Ziegler and J. Zupan, _Quark Flavor Phenomenology of the QCD Axion_ , _Phys. Rev. D_ 102 (2020) 015023 [2002.04623].
  * (57) J.M. Camalich, J. Terol-Calvo, L. Tolos and R. Ziegler, _Supernova Constraints on Dark Flavored Sectors_ , _Phys. Rev. D_ 103 (2021) L121301 [2012.11632].
  * (58) T. Vonk, F.-K. Guo and U.-G. Meißner, _The axion-baryon coupling in SU(3) heavy baryon chiral perturbation theory_ , _JHEP_ 08 (2021) 024 [2104.10413].
  * (59) R. Bollig, W. DeRocco, P.W. Graham and H.-T. Janka, _Muons in supernovae: implications for the axion-muon coupling_ , _Phys. Rev. Lett._ 125 (2020) 051104 [2005.07141].
  * (60) L. Calibbi, D. Redigolo, R. Ziegler and J. Zupan, _Looking forward to lepton-flavor-violating ALPs_ , _JHEP_ 09 (2021) 173 [2006.04795].
  * (61) D. Croon, G. Elor, R.K. Leane and S.D. McDermott, _Supernova Muons: New Constraints on $Z$’ Bosons, Axions and ALPs_, _JHEP_ 01 (2021) 107 [2006.13942].
  * (62) X.G. He, G.C. Joshi, H. Lew and R.R. Volkas, _NEW Z-prime PHENOMENOLOGY_ , _Phys. Rev. D_ 43 (1991) 22.
  * (63) X.-G. He, G.C. Joshi, H. Lew and R.R. Volkas, _Simplest Z-prime model_ , _Phys. Rev. D_ 44 (1991) 2118.
  * (64) E. Ma, D.P. Roy and S. Roy, _Gauged L(mu) - L(tau) with large muon anomalous magnetic moment and the bimaximal mixing of neutrinos_ , _Phys. Lett. B_ 525 (2002) 101 [hep-ph/0110146].
  * (65) S. Baek, N.G. Deshpande, X.G. He and P. Ko, _Muon anomalous g-2 and gauged L(muon) - L(tau) models_ , _Phys. Rev. D_ 64 (2001) 055006 [hep-ph/0104141].
  * (66) E. Salvioni, A. Strumia, G. Villadoro and F. Zwirner, _Non-universal minimal Z’ models: present bounds and early LHC reach_ , _JHEP_ 03 (2010) 010 [0911.1450].
  * (67) J. Heeck and W. Rodejohann, _Gauged $L_{\mu}-L_{\tau}$ Symmetry at the Electroweak Scale_, _Phys. Rev. D_ 84 (2011) 075007 [1107.5238].
  * (68) W.-Z. Feng, P. Nath and G. Peim, _Cosmic Coincidence and Asymmetric Dark Matter in a Stueckelberg Extension_ , _Phys. Rev. D_ 85 (2012) 115016 [1204.5752].
  * (69) K. Harigaya, T. Igari, M.M. Nojiri, M. Takeuchi and K. Tobe, _Muon g-2 and LHC phenomenology in the $L_{\mu}-L_{\tau}$ gauge symmetric model_, _JHEP_ 03 (2014) 105 [1311.0870].
  * (70) M. Bauer, P. Foldenauer and J. Jaeckel, _Hunting All the Hidden Photons_ , _JHEP_ 07 (2018) 094 [1803.05466].
  * (71) D.W.P. Amaral, D.G. Cerdeno, A. Cheek and P. Foldenauer, _Confirming $U(1)_{L_{\mu}-L_{\tau}}$ as a solution for $(g-2)_{\mu}$ with neutrinos_, _Eur. Phys. J. C_ 81 (2021) 861 [2104.03297].
  * (72) A. Greljo, Y. Soreq, P. Stangl, A.E. Thomsen and J. Zupan, _Muonic force behind flavor anomalies_ , _JHEP_ 04 (2022) 151 [2107.07518].
  * (73) A. Greljo, P. Stangl, A.E. Thomsen and J. Zupan, _On $(g-2)_{\mu}$ from gauged $U(1)_{X}$_, _JHEP_ 07 (2022) 098 [2203.13731].
  * (74) P. Langacker, _The Physics of Heavy $Z^{\prime}$ Gauge Bosons_, _Rev. Mod. Phys._ 81 (2009) 1199 [0801.1345].
  * (75) P. Foldenauer, _Light dark matter in a gauged $U(1)_{L_{\mu}-L_{\tau}}$ model_, _Phys. Rev. D_ 99 (2019) 035007 [1808.03647].
  * (76) K. Asai, S. Okawa and K. Tsumura, _Search for $\mathrm{U}{(1)}_{L_{\mu}-{L}_{\tau}}$ charged dark matter with neutrino telescope_, _JHEP_ 03 (2021) 047 [2011.03165].
  * (77) I. Holst, D. Hooper and G. Krnjaic, _Simplest and Most Predictive Model of Muon g-2 and Thermal Dark Matter_ , _Phys. Rev. Lett._ 128 (2022) 141802 [2107.09067].
  * (78) M. Drees and W. Zhao, _U(1)(L(mu)-L(tau)) for light dark matter, g-2, the 511 keV excess and the Hubble tension_ , _Phys. Lett. B_ 827 (2022) 136948 [2107.14528].
  * (79) T. Hapitas, D. Tuckler and Y. Zhang, _General kinetic mixing in gauged U(1) $L_{\mu}-L_{\tau}$ model for muon g-2 and dark matter_, _Phys. Rev. D_ 105 (2022) 016014 [2108.12440].
  * (80) J.A. Carpio, K. Murase, I.M. Shoemaker and Z. Tabrizi, _High-energy cosmic neutrinos as a probe of the vector mediator scenario in light of the muon g-2 anomaly and Hubble tension_ , _Phys. Rev. D_ 107 (2023) 103057 [2104.15136].
  * (81) J. Heeck and A. Thapa, _Explaining lepton-flavor non-universality and self-interacting dark matter with $L_{\mu}-L_{\tau}$_, _Eur. Phys. J. C_ 82 (2022) 480 [2202.08854].
  * (82) X. Qi, W. Liu and H. Sun, _Leptogenesis and light scalar dark matter in a $L_{\mu}-L_{\tau}$ model_, 2204.01086.
  * (83) K. Deka, S. Sadhukhan and M.P. Singh, _GeV Dark Matter in $U(1)_{L_{\mu}-L_{\tau}}$: Role of $(g-2)_{\mu}$ Anomaly_, 2203.17122.
  * (84) S. Baek, J. Kim and P. Ko, _Muon $(g-2)$ and Thermal WIMP DM in ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ Models_, 2204.04889.
  * (85) M. Escudero, T. Schwetz and J. Terol-Calvo, _A seesaw model for large neutrino masses in concordance with cosmology_ , _JHEP_ 02 (2023) 142 [2211.01729].
  * (86) D. Hooper, J. Iguaz Juan and P.D. Serpico, _Signals of a New Gauge Boson from IceCube and Muon $g-2$_, 2302.03571.
  * (87) T. Aoyama et al., _The anomalous magnetic moment of the muon in the Standard Model_ , _Phys. Rept._ 887 (2020) 1 [2006.04822].
  * (88) Muon g-2 collaboration, _Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL_ , _Phys. Rev. D_ 73 (2006) 072003 [hep-ex/0602035].
  * (89) Muon g-2 collaboration, _Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm_ , _Phys. Rev. Lett._ 126 (2021) 141801 [2104.03281].
  * (90) S. Borsanyi et al., _Leading hadronic contribution to the muon magnetic moment from lattice QCD_ , _Nature_ 593 (2021) 51 [2002.12347].
  * (91) C. Alexandrou et al., _Lattice calculation of the short and intermediate time-distance hadronic vacuum polarization contributions to the muon magnetic moment using twisted-mass fermions_ , 2206.15084.
  * (92) Fermilab Lattice, MILC, HPQCD collaboration, _Windows on the hadronic vacuum polarization contribution to the muon anomalous magnetic moment_ , _Phys. Rev. D_ 106 (2022) 074509 [2207.04765].
  * (93) M. Cè et al., _Window observable for the hadronic vacuum polarization contribution to the muon g-2 from lattice QCD_ , _Phys. Rev. D_ 106 (2022) 114502 [2206.06582].
  * (94) D. Brout et al., _The Pantheon+ Analysis: Cosmological Constraints_ , _Astrophys. J._ 938 (2022) 110 [2202.04077].
  * (95) J.L. Bernal, L. Verde and A.G. Riess, _The trouble with $H_{0}$_, _JCAP_ 10 (2016) 019 [1607.05617].
  * (96) M. Escudero, D. Hooper, G. Krnjaic and M. Pierre, _Cosmology with A Very Light L μ $-$ Lτ Gauge Boson_, _JHEP_ 03 (2019) 071 [1901.02010].
  * (97) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri et al., _In the realm of the Hubble tension—a review of solutions_ , _Class. Quant. Grav._ 38 (2021) 153001 [2103.01183].
  * (98) G. Krnjaic, G. Marques-Tavares, D. Redigolo and K. Tobioka, _Probing Muonphilic Force Carriers and Dark Matter at Kaon Factories_ , _Phys. Rev. Lett._ 124 (2020) 041802 [1902.07715].
  * (99) Y. Kahn, G. Krnjaic, N. Tran and A. Whitbeck, _M 3: a new muon missing momentum experiment to probe (g $-$ 2)? and dark matter at Fermilab_, _JHEP_ 09 (2018) 153 [1804.03144].
  * (100) I. Galon, E. Kajamovitz, D. Shih, Y. Soreq and S. Tarem, _Searching for muonic forces with the ATLAS detector_ , _Phys. Rev. D_ 101 (2020) 011701 [1906.09272].
  * (101) K. Asai, K. Hamaguchi, N. Nagata, S.-Y. Tseng and J. Wada, _Probing the L $\mu$-L$\tau$ gauge boson at the MUonE experiment_, _Phys. Rev. D_ 106 (2022) L051702 [2109.10093].
  * (102) H. Sieber, D. Banerjee, P. Crivelli, E. Depero, S.N. Gninenko, D.V. Kirpichnikov et al., _Prospects in the search for a new light Z’ boson with the NA64 $\mu$ experiment at the CERN SPS_, _Phys. Rev. D_ 105 (2022) 052006 [2110.15111].
  * (103) T. Bandyopadhyay, S. Chakraborty and S. Trifinopoulos, _Displaced searches for light vector bosons at Belle II_ , _JHEP_ 05 (2022) 141 [2203.03280].
  * (104) C. Cesarotti, S. Homiller, R.K. Mishra and M. Reece, _Probing New Gauge Forces with a High-Energy Muon Beam Dump_ , 2202.12302.
  * (105) T. Moroi and A. Niki, _Leptophilic Gauge Bosons at Lepton Beam Dump Experiments_ , 2205.11766.
  * (106) NA64 collaboration, _Search for a light Z’ in the L $\mu$-L$\tau$ scenario with the NA64-e experiment at CERN_, _Phys. Rev. D_ 106 (2022) 032015 [2206.03101].
  * (107) C. Rella, B. Döbrich and T.-T. Yu, _Searching for muonphilic dark sectors with proton beams_ , _Phys. Rev. D_ 106 (2022) 035023 [2205.09870].
  * (108) A. Ariga, R. Balkin, I. Galon, E. Kajomovitz and Y. Soreq, _Hunting muonic forces at emulsion detectors_ , 2305.03102.
  * (109) H.-Y. Chiu and R.C. Stabler, _Emission of Photoneutrinos and Pair Annihilation Neutrinos from Stars_ , _Phys. Rev._ 122 (1961) 1317.
  * (110) V. Ritus, _Photoproduction of neutrinos on electrons and neutrino radiation from stars_ , _Journal of Experimental and Theoretical Physics_ 14 (1961) 915–920.
  * (111) J.B. Adams, M.A. Ruderman and C.-H. Woo, _Neutrino Pair Emission by a Stellar Plasma_ , _Physical Review_ 129 (1963) 1383.
  * (112) V. Petrosian, G. Beaudet and E.E. Salpeter, _Photoneutrino Energy Loss Rates_ , _Phys. Rev._ 154 (1967) 1445.
  * (113) G. Beaudet, V. Petrosian and E.E. Salpeter, _Energy Losses due to Neutrino Processes_ , _Astrophys. J._ 150 (1967) 979.
  * (114) D.A. Dicus, _Stellar energy-loss rates in a convergent theory of weak and electromagnetic interactions_ , _Phys. Rev. D_ 6 (1972) 941.
  * (115) P. Sutherland, J.N. Ng, E. Flowers, M. Ruderman and C. Inman, _Astrophysical Limitations on Possible Tensor Contributions to Weak Neutral Current Interactions_ , _Phys. Rev. D_ 13 (1976) 2700.
  * (116) H. Munakata, Y. Kohyama and N. Itoh, _Neutrino Energy Loss in Stellar Interiors_ , _Astrophys. J._ 296 (1985) 197.
  * (117) P.J. Schinder, D.N. Schramm, P.J. Wiita, S.H. Margolis and D.L. Tubbs, _Neutrino Emission by the Pair, Plasma, and Photo Processes in the Weinberg-Salam Model_ , _Astrophys. J._ 313 (1987) 531.
  * (118) N. Itoh, T. Adachi, M. Nakagawa, Y. Kohyama and H. Munakata, _Neutrino Energy Loss in Stellar Interiors. III. Pair, Photo-, Plasma, and Bremsstrahlung Processes_ , _Astrophys. J._ 339 (1989) 354.
  * (119) E. Braaten, _Neutrino Emissivity of an Ultrarelativistic Plasma from Positron and Plasmino Annihilation_ , _Astrophys. J._ 392 (1992) 70.
  * (120) E. Braaten and D. Segel, _Neutrino energy loss from the plasma process at all temperatures and densities_ , _Phys. Rev. D_ 48 (1993) 1478 [hep-ph/9302213].
  * (121) M. Haft, G. Raffelt and A. Weiss, _Standard and Nonstandard Plasma Neutrino Emission Revisited_ , _Astrophys. J._ 425 (1994) 222 [astro-ph/9309014].
  * (122) N. Itoh, H. Hayashi, A. Nishikawa and Y. Kohyama, _Neutrino energy loss in stellar interiors. vii. pair, photo-, plasma, bremsstrahlung, and recombination neutrino processes_ , _Astrophysical Journal Supplement Series_ (1996) .
  * (123) S. Esposito, G. Mangano, G. Miele, I. Picardi and O. Pisanti, _Neutrino energy loss rate in a stellar plasma_ , _Nuclear Physics B_ 658 (2003) 217 [astro-ph/0301438].