INT-PUB-21-024

# Light scalars in neutron star mergers

P. S. Bhupal Dev  Jean-François Fortin  Steven P. Harris  Kuver Sinha
Yongchao Zhang

###### Abstract

Due to their unique set of multimessenger signals, neutron star mergers have
emerged as a novel environment for studies of new physics beyond the Standard
Model (SM). As a case study, we consider the simplest extension of the SM
scalar sector involving a light CP-even scalar singlet $S$ mixing with the SM
Higgs boson. These $S$ particles can be produced abundantly in neutron star
mergers via the nucleon bremsstrahlung process. We show that the $S$ particles
may either be trapped in or stream freely out of the merger remnant, depending
on the $S$ mass, its mixing with the SM Higgs boson, and the temperature and
baryon density in the merger. In the free-streaming region, the scalar $S$
will provide an extra channel to cool down the merger remnant, with cooling
timescales as small as ${\cal O}$(ms). On the other hand, in the trapped
region, the Bose gas of $S$ particles could contribute a larger thermal
conductivity than the trapped neutrinos in some parts of the parameter space,
thus leading to faster thermal equilibration than expected. Therefore, future
observations of the early postmerger phase of a neutron star merger could
effectively probe a unique range of the $S$ parameter space, largely
complementary to the existing and future laboratory and supernova limits. In
view of these results, we hope the merger simulation community will be
motivated to implement the effects of light CP-even scalars into their
simulations in both the free-streaming and trapped regimes.

##  1 Introduction

Dense astrophysical environments have proved extraordinarily useful in
constraining light beyond the Standard Model (BSM) particles (such as axions
and dark photons) [1, 2, 3]. If the BSM particles are sufficiently light and
couple to nucleons, leptons, or photons, they can be readily produced in the
hot and dense nuclear matter in the cores of stars, supernovae and neutron
stars. If the BSM particles interact very weakly with the SM particles, they
will free-stream out of the astrophysical bodies, thus providing an additional
energy loss mechanism that would affect their evolution. Limits can then be
derived on the couplings of these BSM particles to nucleons, leptons and
photons from the requirement that stellar lifetimes and energy-loss rates
should not conflict with observations. There exists a suite of astrophysical
limits from the Sun, from evolved low-mass stars—such as red giants and
horizontal-branch stars in globular clusters, or white dwarfs—as well as from
neutron stars and supernova 1987A [4].

Since the 2017 detection of gravitational waves from the inspiral phase of the
neutron star merger GW170817 [5] along with an electromagnetic counterpart [6,
7], it has become clear that neutron star mergers provide an interestingly new
astrophysical environment to study light BSM physics [8], apart from being a
powerful tool for exploring other fundamental physics topics [9, 10], such as
constraining the nuclear equation of state (EoS) [11, 12], understanding the
$r$-process [13], measuring the Hubble constant [14], studying dark matter
capture [15], placing bounds on the violation of Lorentz invariance [7] and on
gravitational parity violation [16], and testing alternative theories of
gravity [17, 18, 19, 20, 21, 22]. Recently, the possible role of axions [23,
24, 25, 26, 8, 27] and dark gauge bosons [28, 29, 30] in the merger
environment has been investigated. In this work we will consider the impact of
a light CP-even scalar on mergers. We will not attempt to set limits on the
properties of scalars from mergers, but instead make a case for implementing
the effects of scalar particles in numerical simulations of mergers, which can
be used in combination with future observations to constrain or detect the
presence of scalars in mergers.

During the merger, the two constituent neutron stars collide and form a hot,
dense remnant containing nuclear matter where BSM particles can be copiously
produced, similar to the environment that exists in a supernova core, though
perhaps reaching higher densities and temperatures [31] (see the discussion in
Section 3). Neutron star mergers also complement supernovae in terms of their
multimessenger signals. Supernovae produce an electromagnetic signal which
provides information about the surface of last scattering for photons, as well
as a neutrino signal which is sensitive to the temperature of the supernovae
neutrinosphere. Core-collapse supernovae with aspherical mass-energy dynamics
are also believed to generate gravitational waves, though a signal has not yet
been detected [32] and is expected to be beyond the current detection
capabilities unless the supernova occurs close by in our galaxy [33]. Neutron
star mergers, on the other hand, produce a stronger gravitational wave signal
from both the inspiral (measured by LIGO and VIRGO [5, 34]) and the
(unmeasured, thus far) postmerger phase [35]. In addition, they produce
electromagnetic signals in the form of a gamma-ray burst [6] and the kilonova
[36]. Neutron star mergers produce neutrinos as well, but since the number of
neutrinos per event is less than for a supernova, and the merger rate is much
less than the supernova rate, the expected number of detections of merger-
originated neutrinos is, even with a network of gravitational wave detectors
and a megaton neutrino detector operating in tandem, a few per century [37,
38, 39].

A program of using neutron star mergers to discover or constrain BSM physics
must identify the impact a BSM particle will have on either the gravitational
wave or the electromagnetic signal coming from the merger. The case of QCD
axions in the context of mergers was investigated by some of us in Ref. [26]
and it was found that the possibility of axion trapping in the merger remnant
is ruled out by the SN1987A constraint [40, 41, 42, 43, 44, 45]. This is also
expected to be the case for very light axion-like particles (ALPs) [46], with
no allowed region left between the supernova and laboratory constraints [47].
They can still be produced copiously in the merger remnant and then escape,
cooling the hottest regions of the remnant in millisecond timescales. However,
simulations including cooling due to radiation of axions showed little effect
on the gravitational signal or the amount of ejected material [25]. The effect
of ultralight axions on neutron star inspiral phase has been studied in Refs.
[23, 24].

In this work, we investigate the impact of light CP-even scalars on the
cooling and thermal transport properties of a neutron star merger. Such CP-
even scalar particles contrast in several ways with CP-odd axions (and ALPs in
general):

  1. i.

The pseudo-scalar nature of the axion leads to qualitatively different
interactions with SM particles. In particular, in the nucleon bremsstrahlung
process, the axions can only be emitted from either initial or final state
nucleons, but not from the mediator pions, because the axion-pion coupling is
forbidden at leading order. On the other hand, a CP-even scalar can be emitted
from any of the nucleons, as well as from the pion mediator, thus making its
production rate very different from the axion case.

  2. ii.

Light axions can oscillate into SM photons when propagating through background
electromagnetic fields. This property leads to modified axion and SM photon
luminosities at the observation point. On the other hand, the CP-even scalars
in our scenario dominantly decay into $e^{+}e^{-}$ pairs (above the 1 MeV
threshold); their coupling to photons and any associated scalar-photon
oscillation effect (see e.g. Ref. [48]) will be loop-suppressed.

  3. iii.

There exist strong astrophysical [49, 50, 51] and laboratory [52, 53] bounds
on axion couplings, which preclude the possibility of QCD axion (or light ALP)
trapping in neutron star mergers. On the contrary, as we will show in this
work, current constraints on the CP-even scalar, mainly from their decay into
SM particles, still allow for the interesting possibility of trapped scalars
in the merger.

Depending on the thermodynamic conditions in the merger, we find that the CP-
even scalar can have a wide range of mean free path (MFP), for instance—for a
particular choice of mixing angle $\sin{\theta}$—the MFP can range from
hundreds of kilometers to below one meter (see Figure 3). The dependence of
the MFP on the scalar mass $m_{S}$ and its mixing angle $\sin\theta$ with the
SM Higgs is exemplified in Figure 4 for some benchmark baryon densities and
temperatures. Even after taking into account all the laboratory and supernova
limits on the scalar, there is still some parameter space left where the
scalar can either be trapped or free-stream out of the merger. In the free-
streaming region, a given scalar will very likely escape the merger remnant
and cool it down. The emissivity is illustrated in Figure 5 and the
corresponding cooling timescales are shown in Figure 6. It is remarkable that
some regions of the parameter space with the cooling timescale of ${\cal
O}(1\,{\rm ms})$ to ${\cal O}(10\,{\rm ms})$ can be directly probed at future
laboratory experiments, such as DUNE (cf. the lower left panel of Figure 6).
For sufficiently small MFP, the scalar $S$ will form a Bose gas and contribute
to thermal conductivity of the neutron star merger remnants. As shown in
Figures 7, 8, and 9, the thermal equilibration time due to the scalar $S$ can
reach ${\cal O}({\rm ms})$ or even smaller in some regions of the parameter
space that are allowed by current laboratory constraints. The thermal
conductivities due to the trapped scalars and trapped neutrinos are compared
in Figure 10. It turns out that in a large region of parameters, in particular
when the temperature is hotter than a few tens of MeV, the thermal
equilibration will be dominated by the scalars, which could potentially lead
to some detectable effects in future neutron star merger observations.

The rest of the paper is organized as follows: In Section 2, we briefly review
the CP-even scalar scenario and its couplings to the SM from a model-
independent perspective. We give an overview of the conditions inside neutron
star mergers in Section 3 and calculate the MFP of the scalar in these
conditions in Section 4. If free-streaming, the scalar will radiatively cool
the hot regions of the merger. In Section 5 we calculate the timescale of this
cooling. If trapped, the scalar could contribute to the thermal equilibration
of the interior of the merger remnant, which we discuss and calculate the
associated timescale in Section 6. Our conclusions and outlook are given in
Section 7. Details on various aspects of the MFP, emissivity, and thermal
conductivity calculations are presented in several appendices. Throughout this
paper, we work in natural units, where $\hbar=c=k_{B}=1$.

##  2 The CP-even scalar

After the discovery of a Higgs-like particle at the LHC [54, 55], the SM is
now complete. However, there are compelling theoretical arguments and
experimental observations suggesting the existence of BSM physics, although
the energy scale at which it could manifest is still unknown. In a bottom-up
phenomenological approach, i.e. without referring to any specific ultraviolet
(UV) completion, the effects of BSM physics could be captured simply by
extending the SM scalar, fermion and/or gauge sectors. The minimal extension
of the SM scalar sector consists of adding a real SM-gauge-singlet scalar
field (denoted hereafter by $S$) that mixes with the SM Higgs boson [56, 57,
58, 59, 60, 61]. This $S$ particle does not have to be a fundamental degree of
freedom and might actually be the remnant of a more complicated scalar sector
at higher scales. Depending on the model parameters, singlet scalars can be
useful for improving the vacuum stability of the SM [62, 63, 64, 65, 66, 67,
68], first order electroweak phase transition and electroweak baryogenesis
[69, 70, 71, 72, 73, 74, 75, 76, 77, 78], addressing the hierarchy problem in
relaxion models [79, 80, 81, 82, 83], addressing the cosmological constant
problem through radiative breaking of classical scale invariance [84, 85, 86],
or even as a mediator which links the dark and visible sectors [87, 88, 89,
90, 91, 92, 93].

When the mass of the new scalar becomes much larger than the electroweak
scale, the theory can be mapped onto an effective field theory [94, 95, 96,
97]. If the scalar mass is around the electroweak scale, collider constraints
on extended Higgs sector apply [98, 99, 100]. However, for a light scalar well
below the electroweak scale, the constraints from low-energy high-intensity
experiments, such as those from meson decays and beam-dump experiments [92,
101, 102, 103, 104, 105], as well as the astrophysical [92, 106, 107] and
cosmological [92, 108] limits offer a powerful probe. Such a light singlet
scalar may couple to the SM particles via the following ways:

  * •

Through mixing with the SM Higgs $h$. In this case, all the couplings to the
SM particles are proportional to the corresponding SM couplings, rescaled by
the mixing angle $\sin\theta$ of $S$ with the SM Higgs $h$ [56, 57, 58, 59,
60, 61].

  * •

With couplings exclusively or predominantly to some specific SM flavors, such
as the hadronic and leptonic scalars [109, 110, 111, 112].

  * •

With couplings to some SM particles through some heavy BSM particle loops,
such as the $SU(2)_{R}$-breaking scalar in the left-right symmetric model
(LRSM) based on the gauge group $SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$
[113, 114, 115], which couples to photons via the heavy $W_{R}$ and charged
scalars in the LRSM [116, 117, 118, 119].

For simplicity we will consider only the first case above, where all the
phenomenologies of $S$ are determined by only two parameters, i.e. its mass
$m_{S}$ and the mixing angle $\sin\theta$ of $S$ with $h$.

In neutron stars, the dominant production channel for the light $S$ is the
nucleon bremsstrahlung process (just like in the supernova case [106])

$\displaystyle N+N^{\prime}\to N+N^{\prime}+S\,,$ (2.1)

where the scalar $S$ can couple to any of the four external nucleon legs
($N,N^{\prime}$ can be either a neutron or a proton) or to the intermediate
pion mediator, as shown in Figure 1. The $S$ coupling strengths are given by
the Lagrangian

$\displaystyle{\cal
L}\supset\sin\theta\left[y_{hNN}S\bar{N}N+A_{\pi}(S\pi^{0}\pi^{0}+S\pi^{+}\pi^{-})\right]\,,$
(2.2)

with $y_{hNN}\simeq 10^{-3}$ being the coupling of $h$ with nucleons [120,
121]. The coupling of $S$ with pions can be calculated from the chiral
perturbation theory [122, 123]:

$\displaystyle A_{\pi}=\frac{m_{S}^{2}}{v_{\rm
EW}}\left(\frac{2}{9}+\frac{11}{9}\frac{m_{\pi}^{2}}{m_{S}^{2}}\right)\,,$
(2.3)

with $v_{\rm EW}=(\sqrt{2}G_{F})^{-1/2}$ being the electroweak scale ($G_{F}$
is the Fermi constant).

Figure 1: Feynman diagrams for the production of light scalar $S$ in the
nucleon bremsstrahlung process $N+N^{\prime}\to N+N^{\prime}+S$ (with
$N,N^{\prime}=p,\,n$) in the neutron star merger. The scalar $S$ can be
attached to any of the nucleon lines $(a)$, $(b)$, $(c)$, $(d)$,
$(a^{\prime})$, $(b^{\prime})$, $(c^{\prime})$, $(d^{\prime})$, as denoted by
the crosses ($\times$), or to the pion mediator $(e)$, $(e^{\prime})$, as
denoted by the blobs ($\bullet$). The left diagram is for the $t$-channel
process, whereas the right diagram is for the $u$-channel, with the final-
state nucleon 4-momenta interchanged.

After being produced, the light scalar $S$ can decay into the SM particles via
its mixing with the SM Higgs $h$. For a sub-GeV scale mass, the scalar $S$ can
decay into charged lepton pairs $\ell^{+}\ell^{-}=e^{+}e^{-}$,
$\mu^{+}\mu^{-}$ and pion pairs $\pi^{0}\pi^{0}$, $\pi^{+}\pi^{-}$ at the tree
level, and into the photon pairs $\gamma\gamma$ at one-loop level. The
corresponding MFPs can be found in Appendix A.

Let us also summarize the existing (and future) laboratory, astrophysical and
cosmological limits on the light scalar $S$, which will be used later for
comparison with our results:

  * •

The scalar $S$ can obtain flavor-changing neutral current (FCNC) couplings to
the SM quarks at one-loop level via mixing with the SM Higgs, and contribute
to rare FCNC decays of mesons, such as $K\to\pi+X$, $B\to K+X$ and $B\to\pi+X$
with $X=e^{+}e^{-},\,\mu^{+}\mu^{-},\,\gamma\gamma$ or missing energy if $X$
decays outside the detectors. There are stringent limits from
$K^{+}\to\pi^{+}e^{+}e^{-},\;\pi^{+}\mu^{+}\mu^{-}$ in NA48/2 [124, 125],
$K^{+}\to\pi^{+}\gamma\gamma,\;\pi^{+}\nu\bar{\nu}$ at NA62 [126, 127] and
E949 [128], $K_{L}\to\pi^{0}\chi\chi$ (with
$\chi\chi=e^{+}e^{-},\;\mu^{+}\mu^{-},\;\gamma\gamma$) in KTeV [129, 130, 131,
132], $K_{L}\to\pi^{0}\nu\bar{\nu},\;\pi^{0}X$ (with $X$ a long-lived
particle) at KOTO [133], $B\to K+X$ with
$X=e^{+}e^{-},\,\mu^{+}\mu^{-},\,\nu\bar{\nu}$ at LHCb [134], BaBar [135,
136], Belle [137]. The MicroBooNE collaboration has recently performed the
searches of light monoenergetic scalars from kaon decay at rest
$K^{+}\to\pi^{+}+S$ with $S\to e^{+}e^{-}$ [138]. The calculation details of
these FCNC decays can be found e.g. in Refs. [119, 103, 102]. The most
stringent limits from NA62, E949, KOTO, MicroBooNE and LHCb are shown in
Figure 4 respectively as the shaded blue, light green, orange, dark green and
pink regions. The remaining limits are relatively weaker and are not shown in
these figures.

  * •

The light scalar $S$ can also be produced in the high-intensity beam-dump
experiments. The current most stringent limits are from kaon decays
$K\to\pi+S$ in the CHARM experiment [139], which is presented as the magenta
shaded regions in Figure 4. At the LSND experiment, the scalar $S$ is
predominately produced by the proton bremsstrahlung process, and the current
LSND electron and muon data [140, 141] have excluded the brown shaded regions
in Figure 4 [104]. The data from the fixed target experiment PS191 have also
be reinterpreted to set limits on the light scalar from kaon decays
$K\to\pi+S$, and the resultant limits are shown as the red shaded regions in
Figure 4.

  * •

On the astrophysical side, the scalar $S$ can be produced abundantly in the
supernova cores with temperature $T\sim 30$ MeV via the nucleon bremsstrahlung
process (2.1). The observation of neutrino luminosity ${\cal L}_{\nu}$ of
SN1987A can be used to set limits on the mixing angle $\sin\theta$ of $S$ with
the SM Higgs [142, 143, 92, 144, 145, 106]. The most recent supernova limits
can be found in Ref. [106]. Setting conservatively ${\cal L}_{\nu}=3\times
10^{53}$ erg/sec and $5\times 10^{53}$ erg/sec, the corresponding supernova
limits on $m_{S}$ and $\sin\theta$ are shown in Figure 4 respectively as the
lighter and darker black shaded regions, which exclude the mixing angle
$\sin\theta$ from roughly $5.9\times 10^{-7}$ to $7.0\times 10^{-6}$ with
scalar mass up to 148 MeV. With a less conservative luminosity limit, say
$3\times 10^{52}$ erg/sec, the SN1987A observations could exclude larger
regions of $m_{S}$ and $\sin\theta$. Inside the Sun, white dwarfs, red giants
and horizontal-branch stars where the temperature is ${\cal O}({\rm keV})$,
the production of $S$ is expected to be dominated by the process $e+N_{i}\to
e+N_{i}+S$ with $N_{i}$ here denoting different nuclei in the stars [107].
However, the limits on $S$ from these keV-temperature stars can reach only up
to ${\cal O}(100\,{\rm keV})$, and are thus not shown in Figure 4.

  * •

In the early universe, if the light $S$ keeps in thermal equilibrium with the
SM particles and has a lifetime $\tau_{S}\gtrsim 1$ sec, it will contribute an
extra degree of freedom and spoil the success of big bang nucleosynthesis
(BBN) [146, 108]. However, it turns out that there is no parameter space of
$m_{S}$ and $\sin\theta$ that satisfies both the conditions above, so the
generic BBN limit is not applicable in our setup [106]. Moreover, the actual
BBN constraint will depend on the particular thermal history of the scalar $S$
in the early universe and can be evaded, without affecting the late-time
phenomenology we are interested in here. Therefore, we do not show the BBN
limits in Figure 4. In addition, a sufficiently light scalar $S$ in thermal
equilibrium contributes also to the light degrees of freedom $N_{\rm eff}$ at
the CMB epoch, and the mixing angle $\sin\theta$ is constrained by the
precision Planck data [147]. However, this limit on the mixing from $\Delta
N_{\rm eff}$ is very weak, of order ${\cal O}(0.01)$ [103], therefore it is
not shown in Figure 4.

  * •

The mixing of light $S$ with the SM Higgs $h$ will contribute to the invisible
decay of the SM Higgs via $h\to SS$; as a result, the mixing angle
$\sin\theta$ is constrained by the precision Higgs data at the LHC. The
current limit of ${\rm BR}(h\to{\rm inv.})<0.11$ [148] leads to the limit of
$\sin\theta<0.079$ [107]. As this limit is comparatively weaker than the other
limits discussed above, it is not shown in Figure 4.

It is expected that the laboratory limits on $S$ will be significantly
improved in future high-intensity experiments. For instance, the Fermilab SBN
detectors, which is a combination of using the SBND and ICARUS experiments
with the BNB and NuMI beams respectively, can probe a mixing angle down to
$1.7\times 10^{-7}$ [149]. The near detector at DUNE can further improve the
prospects by over one order of magnitude [150], assuming at least 10 signal
events. Similarly, the FASER 2 experiment at the LHC can detect light long-
lived particles produced in the high-energy $pp$ collisions. Benefiting from
the high colliding energy, the FASER 2 experiment can probe a relatively
heavier $S$ with mass up to few GeV [151, 152]. The corresponding prospects of
$m_{S}$ and $\sin\theta$ at SNB, DUNE and FASER 2 are denoted respectively by
the purple, green and light blue lines in Figure 4. The light scalar $S$ can
also be searched for in some other future high-intensity experiments, such as
NA62 [153], LHCb [154], SHiP [155], CODEX-b [154] and MATHUSLA [156, 157]. For
the sake of clarity, these prospects are not shown in Figure 4. For a summary
of these constraints, see e.g. [158].

##  3 Nuclear matter in neutron star mergers

The matter inside cold neutron stars is a degenerate Fermi liquid of neutrons
and protons along with a nearly ideal, degenerate Fermi gas of electrons and
muons. We ignore the possibility of exotic phases of matter, such as those
containing quarks, hyperons, quarkyonic matter, or Bose condensates [159, 160,
161, 162]. In the interior of the star, the density can reach several times
the nuclear saturation density $n_{0}\equiv 0.16\text{ fm}^{-3}$. Throughout
the neutron star inspiral, there is little change in the thermodynamic
conditions of the stars. When the two stars collide, numerical simulations
[163, 164, 165] indicate that significant shock heating occurs and the
temperature of the nuclear matter at the collision interface rises to several
tens of MeV, while the nuclear matter in the densest part of the neutron star
cores only heats modestly, to perhaps 5–10 MeV [166, 167, 168, 169].
Therefore, in contrast to supernovae (nondegenerate) or cold neutron stars
(degenerate), the nuclear matter inside the neutron star merger remnant spans
the entire range of degeneracies [26, 170].

In this work, we will focus on the regions of the neutron star merger remnant
where the temperature exceeds 5–10 MeV. In these conditions, the neutrino MFP
is well under one kilometer [171, 172] and therefore neutrinos are trapped
inside the remnant through their charged and neutral current interactions, and
attain a Fermi-Dirac distribution. We will consider the neutron star merger
matter to consist of neutrons $n$, protons $p$, electrons $e^{-}$, muons
$\mu^{-}$, electron neutrinos $\nu_{e}$, and muon neutrinos $\nu_{\mu}$ (and
the corresponding antileptons). The content of this matter at a given
temperature $T$ is determined by choosing a baryon density $n_{B}$, an
electron lepton fraction $Y_{L_{e}}\equiv(n_{e}+n_{\nu_{e}})/n_{B}$, and a
muon lepton fraction $Y_{L_{\mu}}\equiv(n_{\mu}+n_{\nu_{\mu}})/n_{B}$. In this
work we choose $Y_{L_{e}}=Y_{L_{\mu}}=0.1$ [167], though the value of
$Y_{L_{\mu}}$ throughout the neutrino-trapped regions of a neutron star merger
remnant is, as far as we know, unstudied. The choice of conserved lepton
fractions has a small effect on the production rate and the MFP of the $S$
particle, typically up to a factor of 2 to 3. This is about the same level of
variation caused by uncertainty in the nuclear EoS (see Appendix F of Ref.
[51] where this was studied for axion production in degenerate nuclear
matter).

We model the nuclear matter with the IUF EoS [173], which is a relativistic
mean field (RMF) theory [174, 175, 176]. In this type of model, the strong
interaction between nucleons is treated as meson exchange, where the meson
fields are frozen to their vacuum expectation values, which vary with density
and temperature. In this mean field approximation, the nucleons have energy
dispersion relationships

$\displaystyle E_{n}$ $\displaystyle=U_{n}+E^{*}_{n}=U_{n}+\sqrt{{\bf
p}_{n}^{2}+m_{*}^{2}}\,,$ (3.1) $\displaystyle E_{p}$
$\displaystyle=U_{p}+E^{*}_{p}=U_{p}+\sqrt{{\bf p}_{p}^{2}+m_{*}^{2}}\,,$
(3.2)

where $U_{n}$ and $U_{p}$ are the nuclear mean fields experienced by the
neutron and proton, respectively. The mean fields $U_{i}$ depend on the baryon
density and temperature, but not the nucleon momentum. In the IUF RMF theory
that we use here, the Dirac effective mass $m_{*}$ of the neutron and proton
are equal. At $n_{B}=n_{0}$, $m_{*}\approx 560$ MeV and as density rises,
$m_{*}$ drops dramatically, falling to $m_{*}\approx 185$ MeV at
$n_{B}=5n_{0}$. As the neutron Fermi momentum rises from 330 to 540 MeV in
this density range, neutrons cannot be considered nonrelativistic for
densities $n_{B}\gtrsim 3n_{0}$ [177]. The leptons are treated as free Fermi
gases. We will discuss the calculation of rate integrals in the RMF theory in
the following section.

##  4 Scalar mean free path

The MFP of the light CP-even scalar $S$ in hot, dense nuclear matter is
primarily determined by the rate of absorption in inverse nucleon
bremsstrahlung processes $N+N^{\prime}+S\rightarrow N+N^{\prime}$ (with
$N,N^{\prime}=n,p$), and to some extent, by the total decay rate of $S$ into
SM particles, such as $S\rightarrow\gamma\gamma,\ e^{+}e^{-},\
\mu^{+}\mu^{-},\ \pi^{0}\pi^{0},\ \pi^{+}\pi^{-}$. In this work, we consider
scalars with masses $m_{S}<500\text{ MeV}$, as scalars in this mass range are
light enough to be produced in significant quantities in neutron star mergers.
However, scalars with masses $m_{S}>2m_{\pi^{0}}$ ($m_{S}>2m_{\pi^{\pm}}$) can
decay into neutral pions (charged pions), yet our nuclear EoS does not account
for thermal populations of pions [178] nor a pion condensate [179, 180, 174].
Inclusion of a thermal population of pions has been shown to slightly increase
the proton fraction and soften the EoS, especially at high temperatures where
the pion population becomes significant [178]. These effects will not
significantly change the MFP of the scalar due to the nucleon bremsstrahlung
processes. However, as we discuss later in this section, the added possibility
of the $S$ decaying to pions can, in certain conditions, noticeably shorten
the MFP of the scalar particle. Additionally, the matrix element for $S$
production from nucleon bremsstrahlung is only valid for scalar masses that
are small compared to the nucleon effective mass (see the description of the
calculation in Appendix E). Therefore, our results for the production rate and
MFP of the scalar particle are most robust for $m_{S}\lesssim 2m_{\pi}$. Since
the MFP of the scalar $S$ is predominantly dictated by absorption via the
inverse nucleon bremsstrahlung processes, we discuss that process here and
leave the contributions of the various decay processes to Appendix A.

The inverse MFP of the CP-even scalar $S$ due to the absorption process
$S+N+N^{\prime}\rightarrow N+N^{\prime}$ is given by

$\displaystyle\lambda_{NN^{\prime}}^{-1}=\int$
$\displaystyle\frac{\mathop{d^{3}{\bf
p}_{1}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{2}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{3}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{4}}}{(2\pi)^{3}}\frac{S_{NN^{\prime}}\sum_{\text{spins}}|\mathcal{M}_{NN^{\prime}}|^{2}}{2^{5}E_{1}^{*}E_{2}^{*}E_{3}^{*}E_{4}^{*}E_{S}}$
(4.1)
$\displaystyle\times(2\pi)^{4}\delta^{4}(p_{S}+p_{1}+p_{2}-p_{3}-p_{4})f_{1}f_{2}(1-f_{3})(1-f_{4}),$

where $p_{i}=(E_{i},{\bf p}_{i})$ are the four-momenta, $f_{i}$ represents a
Fermi-Dirac factor and $E_{i}^{*}=\sqrt{{\bf p}_{i}^{2}+m_{*}^{2}}$ the energy
of the $i$th nucleon in the reaction (without the nuclear mean field
contribution), $E_{S}=\sqrt{{\bf p}_{S}^{2}+m_{S}^{2}}$ is the energy of $S$,
and $S_{NN^{\prime}}$ is a symmetry factor for identical particles. To
calculate the matrix elements ${\cal M}_{NN^{\prime}}$ for the nucleon
bremsstrahlung process, the nucleon-nucleon scattering via the strong
interaction was modeled with the one-pion-exchange (OPE) approximation. For
the reactions $S+n+n\rightarrow n+n$, $S+p+p\rightarrow p+p$ and
$S+n+p\rightarrow n+p$, the matrix elements are respectively given by [106]

$\displaystyle S_{nn/pp}\sum_{\text{spins}}|\mathcal{M}_{nn/pp}|^{2}$
$\displaystyle=\frac{64\pi^{2}\alpha_{\pi}^{2}f^{4}\sin^{2}{\theta}}{m_{*}^{2}}$
$\displaystyle\times\left(\frac{1}{4}\frac{y_{hnn}^{2}m_{S}^{4}}{E_{S}^{4}}I_{A}^{nn/pp}+\frac{1}{2}\frac{m_{*}^{2}}{81v_{\rm
EW}^{2}}I_{B}^{nn/pp}+\frac{2y_{hnn}m_{S}^{2}m_{*}}{9E_{S}^{2}v_{\rm
EW}}I_{C}^{nn/pp}\right),$ (4.2) $\displaystyle
S_{np}\sum_{\text{spins}}|\mathcal{M}_{np}|^{2}$
$\displaystyle=\frac{256\pi^{2}\alpha_{\pi}^{2}f^{4}\sin^{2}{\theta}}{m_{*}^{2}}$
$\displaystyle\times\left(\frac{1}{4}\frac{y_{hnn}^{2}m_{S}^{4}}{E_{S}^{4}}I_{A}^{np}+\frac{1}{2}\frac{m_{*}^{2}}{81v_{\rm
EW}^{2}}I_{B}^{np}+\frac{2y_{hnn}m_{S}^{2}m_{*}}{9E_{S}^{2}v_{\rm
EW}}I_{C}^{np}\right)\,.$ (4.3)

The symmetry factors for identical particles are $S_{nn/pp}=1/4$ and
$S_{np}=1$. The nucleon-pion coupling $f\approx 1$ and
$\alpha_{\pi}\equiv(2m_{N}/m_{\pi})^{2}/(4\pi)\approx 15$ (the nucleon and
pion masses in $\alpha_{\pi}$ are the values in vacuum). The prefactors of
$1/4$ and $1/2$ in the $I_{A}$ and $I_{B}$ terms respectively were absent in
the earlier calculation [106], and come from relativistic corrections to the
nucleons, as discussed in Appendix E. The dimensionless functions $I_{A}$,
$I_{B}$, and $I_{C}$ can be written in terms of the scalar mass, the pion
mass, and the nucleon momentum transfers
$\mathbf{k}=\mathbf{p}_{2}-\mathbf{p}_{4}$ and
$\mathbf{l}=\mathbf{p}_{2}-\mathbf{p}_{3}$, and are listed in Appendix B. In
the phase space integral Eq. (4.1), the factors of energy in the denominator
(or $E^{*}$ for the nucleons) relate to the normalization of the fermion
spinors [171]. The energies in the four-momentum conserving $\delta$-function
are the true particle energies $E$, not $E^{*}$. The twelve-dimensional phase
space integral in Eq. (4.1) can be reduced to a five-dimensional integral, the
details of which are given in Appendix C.

The total MFP of $S$ is calculated by adding together (in inverse [181]) the
MFP of each absorption and decay channel, i.e.

$\lambda_{S}^{-1}=\lambda_{nn}^{-1}+\lambda_{pp}^{-1}+\lambda_{np}^{-1}+\lambda_{\gamma\gamma}^{-1}+\lambda_{e^{+}e^{-}}^{-1}+\lambda_{\mu^{+}\mu^{-}}^{-1}+\lambda_{\pi^{0}\pi^{0}}^{-1}+\lambda_{\pi^{+}\pi^{-}}^{-1}.$
(4.4)

The MFP is a function of the baryon density $n_{B}$, the temperature $T$, the
mixing angle $\sin{\theta}$, and the mass $m_{S}$ and energy $E_{S}$ of the
scalar particle. To illustrate the “typical” MFP of the $S$ particle at a
given temperature and density, we choose to evaluate the MFP at the average
energy $\langle E_{S}\rangle$ of the scalar particle emitted in those
thermodynamic conditions,

$\langle E_{S}\rangle=Q_{S}/\Gamma_{S},$ (4.5)

where $\Gamma_{S}$ is the production rate of scalars per unit volume, and
$Q_{S}$ is the rate of energy production in scalars per unit volume (the
emissivity, see Eq. (5.3) below). The production rate per unit volume
$\Gamma_{S}$ is given by the same phase space integral as $Q_{S}$ (Eq. 5.4)
but without the factor of $E_{S}$. The average energy of the produced $S$
particle is plotted at nuclear saturation density $n_{0}$ in Figure 2. Low-
mass scalars are emitted with typical energies $\langle E_{S}\rangle\approx
3T$, independent of density (assuming the nuclear matter is still strongly
degenerate). This low-mass limit seems to match the result obtained for axions
produced by nucleon bremsstrahlung [142]. As the scalar mass increases, the
rest mass becomes the dominant contribution to the energy of the scalar, but
the energy still increases with temperature on top of the rest mass
contribution, although with a slope less than in the low-mass limit.

Figure 2: Average energy $\langle E_{S}\rangle$ of $S$ particles produced at
$n_{B}=n_{0}$ as a function of temperature $T$ of the fluid element. This
result is independent of $\sin{\theta}$. The nuclear matter is assumed to be
strongly degenerate, making this result independent of density. In the limit
of $m_{S}\to 0$, the average energy $\langle E_{S}\rangle\approx 3T$.

Figure 3: MFP of the scalar $S$ in the $n_{B}-T$ plane with the mixing angle
fixed at $\sin\theta=10^{-5}$, and the scalar mass $m_{S}=1$ MeV (upper left),
20 MeV (upper right), 100 MeV (lower left) and 270 MeV (lower right). For the
trapped (shaded grey) and free streaming (shaded cyan) regions, the MFP is
assumed to be smaller than 100 m and larger than 1 km respectively, although
this demarcation is somewhat arbitrary.

In Figure 3 we present contour plots of the MFP $\lambda_{S}$ of the scalar in
the $n_{B}$ (in units of $n_{0}$)–$T$ plane, for four different values of the
scalar masses $m_{S}=1$ MeV (upper left panel), 20 MeV (upper right panel),
100 MeV (lower left panel) and 270 MeV (lower right panel). The energy of the
scalar is chosen according to Eq. (4.5). We have chosen $\sin{\theta}=10^{-5}$
in each panel, but since $\lambda^{-1}\propto\sin^{2}{\theta}$, it is easy to
get the MFP for other values of $\sin{\theta}$. In each panel, we label
regions where the MFP is longer than 1 km as “free-streaming”, i.e. the $S$
particles created in a particular region of the merger escapes that region—and
perhaps the merger remnant as a whole—without interacting with the nuclear
matter. Regions where the MFP is shorter than 100 m are labeled as “trapped”,
where the $S$ particles form a thermal Bose-Einstein distribution inside a
region of the merger remnant. The situation in between these two regimes is
difficult to treat—scalar particles can be absorbed or decay, but not
frequently enough where they would establish a Bose-Einstein distribution.
Please note that the demarcation chosen here between the trapped and free-
streaming regions is just for concreteness, and our main results will not be
significantly affected by a slightly different choice.

For light scalars, like the $m_{S}=1\text{ MeV}$ case in the upper left panel
of Figure 3, the MFP contours indicate that as the conditions in the merger
remnant get hotter and denser, the MFP grows shorter. At high densities and
low temperatures, where the nuclear matter is strongly degenerate, the MFP
becomes strongly dependent on temperature and weakly dependent on density,
which is expected in degenerate systems because the number of available states
for nucleons to scatter into near the Fermi surface grows as $T$ for each
degenerate nucleon in the reaction [182, 179]. At low densities and high
temperatures, where the nuclear matter is nondegenerate, the MFP depends
somewhat strongly on density and has a weaker temperature dependence. This
behavior of the MFP in nondegenerate matter can be seen in Eq. (4.3) in Ref.
[106]. In fact, the features of this plot match closely the MFP of light QCD
axions (see Figure 2 in Ref. [26]). As the scalar mass increases, the MFP
shrinks and becomes largely dependent only on the baryon density, even in the
degenerate regime.

Figure 4: MFP of the scalar $S$ in the $m_{S}-\sin{\theta}$ plane, with baryon
density $n_{B}=n_{0}$ (left) or $5n_{0}$ (right), and temperature $T=10$ MeV
(upper), 40 MeV (middle) or 100 MeV (lower). For the trapped (shaded grey) and
free streaming (shaded cyan) regions, the MFP is smaller than 100 m and larger
than 1 km, respectively. The current laboratory and astrophysical constraints
are overlaid as the other shaded regions, and the future prospects at SBN,
DUNE and FASER 2 are shown as the solid purple, green and light blue lines
(see more details in Section 2).

In Figure 4 we plot the MFP of $S$ in the $m_{S}-\sin{\theta}$ plane,
evaluated at the average $S$ energy given in Eq. (4.5). We show the results at
several different densities and temperatures possibly encountered in neutron
star mergers: the three left panels are for the baryon density $n_{B}=n_{0}$,
while the three right panels are at $n_{B}=5n_{0}$. The upper, middle and
lower panels are with the temperatures $T=10$ MeV, 40 MeV and 100 MeV,
respectively. The MFP of the scalar is dominated by absorption from inverse
nucleon bremsstrahlung processes $S+N+N^{\prime}\rightarrow N+N^{\prime}$. All
contributions to the MFP go as $\lambda^{-1}\propto\sin^{2}{\theta}$, making
the mixing angle the dominant factor in the $m_{S}-\sin{\theta}$ plane. In all
panels, as the mixing angle $\sin{\theta}$ grows larger, the MFP of the scalar
shrinks. Likewise, though much less dramatically, as the scalar mass grows,
the MFP shrinks. This behavior was seen in nondegenerate nuclear matter (see
Figure 4 in Ref. [106]). The contribution to the MFP from decay processes is
negligible, with the exception of the decay to pions, which is only possible
when the scalar mass exceeds twice the pion mass, roughly 270 MeV. This decay
becomes significant at high temperatures, because at these temperatures the
thermal pion population is large and the decay of $S$ is subject to Bose
enhancement [45, 1, 183]. A small but sudden decrease in the MFP can be seen
in the bottom left plot ($n_{B}=1n_{0}$ and $T=100\text{ MeV}$) in Figure 4 at
twice the pion mass. At lower temperatures, the thermal population of pions is
not large enough to stimulate sufficient decay to pions to affect the MFP. At
higher densities, the $\pi^{-}$ would likely form a condensate as
$\mu_{\pi^{-}}=\mu_{n}-\mu_{p}$ rises above $m_{\pi^{-}}$, complicating the
calculation of the scalar decay rate to pions, so we neglect this decay
channel at high densities.

Current constraints on the scalar mass $m_{S}$ and the mixing angle
$\sin\theta$ from the laboratory and supernova observations discussed in
Section 2 are displayed on top of the MFP contours in Figure 4 as the shaded
regions. The representative prospects at future experiments SBN, DUNE and
FASER 2 are also indicated by the solid purple, green and light blue lines,
respectively. It is clear that there is unconstrained parameter space in both
the trapped and free-streaming regions of a neutron star merger, where mergers
can in principle provide complementary information. In particular, we find
that mergers have the potential to probe regions of the parameter space of
this model that no other experiment can, including the future experiments
DUNE, SBN, and FASER 2. This motivates a more detailed study of the role of
scalars in the regions of the merger where they free-stream and the regions
where they are trapped. In the rest of the paper, we pursue this analysis,
following the study in the case of axions [26].

##  5 Scalar contribution to cooling

If the scalar is produced in a region of the merger and the MFP of the scalar
is longer than the size of that region, the scalar has a high likelihood of
escaping, cooling down that part of the merger. If the scalar continues to
encounter regions where it has a long MFP, it is likely that the scalar will
free-stream through the merger entirely. The changing conditions the scalar
encounters throughout its trajectory are encapsulated in the optical depth
[184]

$\displaystyle\tau(r,E_{S})=\int_{r}^{\infty}\mathop{dr^{\prime}}\lambda_{S}^{-1}(r^{\prime},E_{S})\,,$
(5.1)

but in this work we will consider trapping or free streaming on the level of
individual fluid elements in the merger, and in this section we will consider
the question of how quickly a fluid element cools when it emits scalar
particles. The decrease in temperature of a fluid element emitting scalars is
related to the emissivity $Q_{S}$ through

$\displaystyle\frac{dT}{dt}\ =\ -\frac{Q_{S}(n_{B},T)}{c_{V}(n_{B},T)}\,,$
(5.2)

where

$\displaystyle Q_{S}=Q_{nn}+Q_{np}+Q_{pp}$ (5.3)

incorporates all the contributions from the $pp$, $nn$ and $np$ bremsstrahlung
processes, with ($NN^{\prime}=pp,\,nn,\,np$)

$\displaystyle Q_{NN^{\prime}}$ $\displaystyle=$
$\displaystyle\int\frac{\mathop{d^{3}{\bf
p}_{1}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{2}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{3}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}}_{4}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{S}}}{(2\pi)^{3}}\frac{S_{NN^{\prime}}\sum_{\text{spins}}|\mathcal{M}_{NN^{\prime}}|^{2}}{2^{5}E_{1}^{*}E_{2}^{*}E_{3}^{*}E_{4}^{*}E_{S}}$
(5.4)
$\displaystyle\times(2\pi)^{4}\delta^{4}(p_{1}+p_{2}-p_{3}-p_{4}-p_{S})f_{1}f_{2}(1-f_{3})(1-f_{4})E_{S}\,,$

with the matrix elements ${\cal M}_{NN^{\prime}}$ same as in Eqs. (4.2) and
(4.3). The specific heat $c_{V}$ per unit volume of degenerate nuclear matter
in Eq. (5.2) is dominated by the neutrons, as they have the largest Fermi
sphere, and is given by [185]

$\displaystyle c_{V}(n_{B},T)\approx\frac{1}{3}m_{n}^{L}p_{Fn}T,$ (5.5)

where $m_{n}^{L}=\sqrt{p_{Fn}^{2}+m_{*}^{2}}$ is the Landau effective mass of
the neutron [186], and $p_{Fn}$ is the Fermi momentum of the neutron. The
specific heat does not deviate much from this expression even at the highest
temperatures encountered in mergers (see Ref. [26] for an expression valid at
arbitrary degeneracy). From Eq. (5.2), one can integrate over the temperature
evolution and estimate the timescale $\tau_{1/2}$ for the fluid element
temperature to drop by a factor of two via

$\displaystyle\tau_{1/2}=\int_{0}^{\tau_{1/2}}dt=\int_{T_{0}/2}^{T_{0}}dT\frac{c_{V}(n_{B},T)}{Q_{S}(n_{B},T)}\,.$
(5.6)

In this integral over temperature, we make the (reasonable) approximation that
the neutron Fermi momentum and Dirac effective mass do not change
significantly with temperature.

Figure 5: Emissivity of the scalar $S$ (left) and resultant radiative cooling
time of fluid elements (right) in the neutron star merger remnant, as a
function of the scalar mass $m_{S}$. In both panels, the solid, short-dashed
and long-dashed curves are respectively for $n_{B}=n_{0}$, $3n_{0}$ and
$5n_{0}$, whereas the black, blue, yellow and red curves represent
respectively $T=20$ MeV, 40 MeV, 75 MeV and 100 MeV. The mixing is fixed at
$\sin{\theta}=10^{-6}$ to ensure that the scalars free-stream in all
thermodynamic conditions for all scalar masses shown in the plot (cf. Figure
4).

In the left panel of Figure 5 we plot the emissivity $Q_{S}$ of the CP-even
scalar as a function of its mass $m_{S}$. The black, blue, yellow and red
curves are respectively for the temperatures $T=20$ MeV, 40 MeV, 75 MeV and
100 MeV, and the solid, short-dashed and long-dashed curves correspond
respectively to the baryon densities of $n_{B}=n_{0}$, $3n_{0}$ and $5n_{0}$.
We see that for all thermodynamic conditions shown in this panel, the
emissivity rises with the scalar mass, reaches a maximum, and then falls off.
This trend is seen in the nondegenerate nuclear matter limit as well (cf.
Figure 2 of Ref. [106]). This mass dependence is more complicated than the
monotonic behavior seen in the emissivity of massive axions [46]. The reason
stems from the CP-even nature of the scalar, thus the matrix element for
$N+N^{\prime}\to N+N^{\prime}+S$ has a nontrivial dependence on the scalar
mass $m_{S}$ [106]. For $m_{S}\gtrsim 10$ MeV, the production of $S$ will be
dominated by the $a^{(^{\prime})}$, $b^{(^{\prime})}$, $c^{(^{\prime})}$ and
$d^{(^{\prime})}$ diagrams in Figure 1 (see Figure 2 of Ref. [106]); as a
result of the partial cancellation of these diagrams, the emissivity $Q_{S}$
will have the mass-dependence in the form of $(m_{S}/E_{S})^{4}$ for $10\text{
MeV}\lesssim m_{S}\lesssim 3T$ (cf. Ref. [106] and Appendix E). This is clear
for all the lines in the left panel of Figure 5. The emissivity increases
strongly with temperature; however, as the mass increases and becomes large
compared to the temperature, the production rate of the scalar becomes
Boltzmann suppressed and falls off as $e^{-m_{S}/T}$. This behavior is evident
for the $T=20\text{ MeV}$ and $T=40\text{ MeV}$ curves, as shown in the left
panel of Figure 5. For the hottest temperatures encountered in mergers, say
$75$ MeV or even $100\text{ MeV}$ in the left panel of Figure 5, the
production rate seems to drop more slowly as the mass $m_{S}$ increases. This
is because at these high temperatures, the Boltzmann suppression is less
dramatic for the range of $m_{S}$ we consider. However, as the mass of the
scalar increases beyond 500 MeV, the emissivity will fall off as
$e^{-m_{S}/T}$. For all the temperatures in Figure 5, the largest emissivity
can be found when the scalar mass $m_{S}\sim 3T$, which is a direct result of
the CP-even nature of the scalar and the resultant $m_{S}$-dependence of the
matrix element of the nucleon bremsstrahlung processes, and is very different
from the axion case. In the right panel of Figure 5, we plot the corresponding
radiative cooling time $\tau_{1/2}$ as a function of the scalar mass $m_{S}$,
for the benchmark temperatures and baryon densities in the left panel of
Figure 5. The nonmonotonicity of the emissivity as a function of $m_{S}$ is
imprinted (in inverse) on the cooling time.

Figure 6: Timescale for fluid elements to cool due to radiating scalar
particles, in the $m_{S}-\sin{\theta}$ plane. Each panel corresponds to a
different density and temperature possibly encountered in a neutron star
merger. The radiative cooling times are only valid for parameters (shown in
cyan) where the scalars free-stream through the matter in the merger. The
parameter space for trapped scalars is also shown, in grey. Constraints from
SN1987a and expected constraints from future collider experiments are
overlaid.

To get a holistic view of the behavior of the radiative cooling time, we plot
the contours of cooling time $\tau_{1/2}$ in the $m_{S}-\sin{\theta}$ plane in
Figure 6. The left and right panels are respectively for the baryon densities
$n_{B}=n_{0}$ and $5n_{0}$. The temperatures chosen reflect the hottest
sections of the merger: in the two upper panels we take $T=40$ MeV, which can
essentially be reached in all numerical simulations for tens of milliseconds
(or until the remnant begins to collapse) [168, 169, 166, 187, 188, 189]; in
the two lower panels the temperature is $T=100$ MeV, which can be present in
some regions of matter in a few simulations for a few milliseconds, and also
as the remnant begins to collapse to a black hole [167, 31]. At these high
temperatures, neutrinos are trapped (and exit the region on the much slower
diffusion timescale [184]) so it would be left to exotic particle species with
long MFPs to radiatively cool the hot fluid elements. Radiative cooling is
only a well-defined concept if the scalars free-stream through the remnant. In
each panel of Figure 6 we color in cyan the parameter space where the scalars
are free-streaming, and the trapped regions are shown in grey. The cooling
time is proportional to $\sin^{2}{\theta}$ (coming from the expression for the
scalar emissivity in Eq. (5.4)), so the cooling timescale dramatically
shortens as the mixing $\sin{\theta}$ increases, up until the MFP shortens to
the point where the scalars cannot free-stream anymore.

Merger remnants can survive for tens, hundreds, or thousands of milliseconds
[190, 191] (or much longer if the remnant mass is below the limit for a
nonrotating star [192, 193]), so radiative cooling that occurs on the
timescales shown in Figure 6 is potentially relevant to the postmerger
evolution. We see from the two upper panels of Figure 6 that areas of the
merger at $T=40\text{ MeV}$ can cool to 20 MeV on timescales shorter than 1
second only if the scalar is relatively light, say $m_{S}\lesssim 250\text{
MeV}$. The mixing angle would need to be roughly
$10^{-7}\lesssim\sin{\theta}\lesssim 10^{-5}$, although a fraction of the
relevant parameter space is ruled out by the SN1987A constraint [106]. If the
remnant is short-lived, say lasting for 10 ms, only scalars with
$m_{S}\lesssim 100\text{ MeV}$ can significantly cool the 40 MeV regions of
the merger before it collapses, and only in a very limited range of
$\sin\theta$. As expected, cooling of a fluid element from $T=100\text{ MeV}$
down to $50\text{ MeV}$ would occur much more rapidly for a given choice of
scalar mass and coupling, as shown in the two lower panels of Figure 6. Rapid
cooling of these extremely hot fluid elements can occur due to radiating
scalars of the full range of masses we discuss in this work. As shown in the
lower left panel of Figure 6, a sizable region of parameter space with ${\cal
O}(1\,{\rm ms})\lesssim\tau_{1/2}\lesssim{\cal O}(10\,{\rm ms})$ can also be
probed at the DUNE experiment.

Cooling of fluid elements due to axion emission was implemented in a neutron
star merger simulation [25]. Even though axion radiation can quickly cool the
merger remnant (see semianalytic calculations in Ref. [26] and the simulation
results in Ref. [25]), it does not lead to significant changes in the dynamics
of the merger remnant. In particular, while the cooling leads to a small
sphericalization of the merger remnant, the gravitational wave signal, remnant
lifetime, and ejecta mass were close to unchanged (compared to the case
without radiative cooling from axion emission), even if the rate of axion
emission was chosen to be unrealistically high [25, 8]. Therefore, we do not
expect radiative cooling to have a direct impact on the dynamics of the
merger, unless it involves physics that was not included in the merger
simulation in Ref. [25]. For example, cooling due to the emission of scalars
could trigger a phase transition in nuclear matter or change the temperature-
dependent transport properties, which could then lead to an observable
signature of the presence of scalars. It is also possible that emitted scalars
can decay once they have escaped the merger remnant. Similar to the discussion
in the case of dark photon emission [30], scalars that escape the merger
remnant can decay into photons directly, or into leptons that scatter and
create photons. These decays may result in distinct $\gamma$-ray signals, the
properties of which deserve further investigation, but beyond the scope of the
current work.

##  6 Scalar contribution to thermal conduction

While experimental data forbid the QCD axion and light ALPs from being trapped
in neutron star merger remnants [26], the CP-even scalar discussed in this
paper has, in certain conditions, a short enough MFP to reach thermal
equilibrium and form a Bose-Einstein distribution. For a given scalar mass and
mixing angle (or, equivalently, coupling strengths), the temperatures and
densities where scalars are trapped are mapped out in Figure 3, and for fixed
densities and temperatures, the values of $m_{S}$ and $\sin{\theta}$ for which
scalars are trapped are mapped out in Figure 4. It is clear that current
laboratory and astrophysical constraints do not rule out the possibility of
trapped scalars.

In regions of the merger remnant where the scalar MFP is sufficiently short,
the scalars will form a free Bose gas with zero chemical potential set by
nucleon bremsstrahlung reactions. Throughout this work, chemical potentials
are defined relativistically, that is, they include the rest mass of the
particle.

In this section, we will calculate the impact of a trapped scalar gas on
transport properties in the merger remnant, and will need a variety of
thermodynamic properties of a quantum gas of scalar particles, which we detail
below.

The number density of a (noninteracting) massive Bose gas with zero chemical
potential is

$n_{S}=\int\dfrac{\mathop{d^{3}{\bf
k}}}{(2\pi)^{3}}\frac{1}{e^{E_{S}/T}-1}=\frac{1}{2\pi^{2}}\int_{m_{S}}^{\infty}\mathop{dE_{S}}\frac{E_{S}\sqrt{E_{S}^{2}-m_{S}^{2}}}{e^{E_{S}/T}-1}\,,$
(6.1)

and the energy density is

$\varepsilon=\int\dfrac{\mathop{d^{3}{\bf
k}}}{(2\pi)^{3}}\frac{E_{S}}{e^{E_{S}/T}-1}=\frac{1}{2\pi^{2}}\int_{m_{S}}^{\infty}\mathop{dE_{S}}\frac{E_{S}^{2}\sqrt{E_{S}^{2}-m_{S}^{2}}}{e^{E_{S}/T}-1}\,.$
(6.2)

The heat capacity at constant volume (or at constant baryon number density, as
baryon number is conserved) is given by

$c_{V,S}=\dfrac{\partial\varepsilon}{\partial
T}=\frac{1}{8\pi^{2}T^{2}}\int_{m_{S}}^{\infty}\mathop{dE_{S}}\dfrac{E_{S}^{3}\sqrt{E_{S}^{2}-m_{S}^{2}}}{\sinh^{2}{[E_{S}/(2T)]}}\,.$
(6.3)

The typical velocity of a scalar particle with mass $m_{S}$ and energy $E_{S}$
is

$v_{S}=\frac{\partial E_{S}(p)}{\partial p}=\sqrt{1-(m_{S}/E_{S})^{2}}\,.$
(6.4)

Trapped species of particles in the merger remnant can participate in various
transport processes, including thermal conductivity, shear viscosity, and
electrical conductivity [194]. Typically particle species with long MFPs (but
still shorter than the system size) will dominate transport, though other
factors are also important for each transport process [194, 195]. In this
section, we will consider the impact of the trapped CP-even scalar on the
thermal conductivity $\kappa$ of the matter in the merger remnant. In the
standard consideration of neutron star mergers, where only neutrons, protons,
electrons, and neutrinos are involved, the thermal conductivity is dominated
by either electrons (in regions of the star where the neutrinos are free-
streaming) or the neutrinos (if they are trapped) [194]. It is useful to
quantify the effect of thermal conductivity by calculating the timescale for
two adjacent fluid elements at slightly different temperatures to reach
thermal equilibrium with each other. A hot fluid element of size $z$ (or
volume $z^{3}$) has excess thermal energy $E_{\text{th}}=c_{V}z^{3}\Delta T$
compared to its neighbors. This heat diffuses through the boundary of the
fluid element (with surface area $6z^{2}$) at rate
$W_{\text{th}}=\kappa(\Delta T/z)6z^{2}$. This thermal equilibration timescale
is given by the ratio of these quantities

$\tau_{\kappa}=\frac{E_{\text{th}}}{W_{\text{th}}}=\frac{c_{V}z^{2}}{6\kappa}\,.$
(6.5)

The specific heat of nuclear matter $c_{V}$ is given by Eq. (5.5). The thermal
conductivity is the sum of the thermal conductivity generated by each particle
species in the remnant

$\kappa_{i}=\frac{1}{3}c_{V,i}v_{i}\lambda_{i}\,,$ (6.6)

where $c_{V,i}$, $v_{i}$ and $\lambda_{i}$ are respectively the specific heat,
speed and MFP of the $i$th particle species. Particle species with the optimal
combination of large density and long MFP (but still shorter than the system
size) dominate the thermal conductivity. In this paper, we focus on the
neutrino-trapped regime with temperature $T>5$ to 10 MeV, so out of the SM
particles, neutrinos dominate the thermal conductivity. The details of the
neutrino thermal conductivity calculation can be found in Appendix F.

Figure 7: Thermal equilibration timescale due to the scalar $S$ over a
distance of 1 km as a function of the scalar mass $m_{S}$. The solid, short-
dashed and long-dashed curves are respectively for $n_{B}=n_{0}$, $3n_{0}$ and
$5n_{0}$, the black, blue, yellow and red curves represent respectively $T=20$
MeV, 40 MeV, 75 MeV and 100 MeV. The mixing $\sin{\theta}=10^{-4}$ is chosen
so that the scalar particles are trapped for all displayed thermodynamic
conditions and scalar masses.

Under certain conditions, trapped CP-even scalars would give rise to a larger
thermal conductivity than the trapped neutrinos. The thermal conductivity of
scalars can be calculated using Eq. (6.6), with input from Eqs. (6.3) and
(6.4), as well as the scalar MFP calculated in Section 4. To obtain an
estimate for the thermal equilibration timescale due to diffusion of a gas of
CP-even scalars in a neutron star merger remnant, we consider a temperature
gradient occurring over a distance scale of $z=1\text{ km}$. For example,
numerical simulations indicate the presence of a ring-like hot region in the
merger remnant, occurring 4 to 7 km out from the remnant’s center [196, 197,
167, 166, 198, 199]. This ring has a thickness of 1 to 3 km. We will show
results only for $z=1\text{ km}$, but since $\tau_{\kappa}\propto z^{2}$, the
thermal equilibration times for other distance scales can be easily obtained.
In particular, it is clear that temperature gradients on shorter distance
scales, for example from turbulence [200], will equilibrate more quickly.

In Figure 7 we plot the timescale for thermal equilibration of matter in the
merger due to diffusion of trapped CP-even scalars (cf. Eq. (6.5)), over a
distance scale of 1 km, as a function of the scalar mass $m_{S}$. As in Figure
5, we have chosen the benchmark temperatures $T=20$ MeV (black), 40 MeV
(blue), 75 MeV (yellow) and 100 MeV (red), and the densities $n_{B}=n_{0}$
(solid), $3n_{0}$ (short-dashed) and $5n_{0}$ (long-dashed). The mixing angle
of $\sin\theta=10^{-4}$ has been chosen such that the scalar $S$ is trapped in
the neutron star merger remnant in the whole mass range we are considering,
cf. Figure 4. It should be noted, however, that for $\sin\theta=10^{-4}$ the
scalar mass range of $220\,{\rm MeV}\lesssim m_{S}\lesssim 310$ MeV has been
excluded by the limits from PS191, CHARM, and LSND (see Figure 4). It is clear
that higher temperature regions in the merger remnant have a higher density of
scalars (cf. Eq. (6.1)) and thus thermally equilibrate faster. The density of
scalars drops as their mass increases, causing the thermal equilibration time
to grow longer as the scalar mass $m_{S}$ increases.

Figure 8: Thermal equilibration timescale due to the scalar $S$ over a
distance of 1 km in the $m_{S}-T$ plane, for the baryon density of
$n_{B}=n_{0}$ (left) and $5n_{0}$ (right). Here $\sin{\theta}=10^{-4}$ is
chosen to ensure that the scalar is trapped for all displayed thermodynamic
conditions and scalar masses.

In Figure 8, we show contour plots of the thermal equilibration time in the
$m_{S}-T$ plane, for two choices of density $n_{B}=n_{0}$ (left panel) and
$5n_{0}$ (right panel), again with $\sin\theta=10^{-4}$ to ensure the scalars
are trapped for all conditions shown in the plots. This presentation of the
thermal equilibration time makes clear that hotter regions of the merger reach
thermal equilibrium more quickly, for the reasons mentioned above, especially
when the scalars are light. The kink in the contours in the left panel at
$m_{S}=2m_{\pi^{\pm}}\simeq 270$ MeV is due to the appearance of the decay
channel $S\rightarrow\pi^{+}\pi^{-}$, which shortens the MFP of the scalar,
slowing its ability to attenuate temperature gradients.

Figure 9: Thermal equilibration timescale due to the scalar $S$ over a
distance of 1 km in the $m_{S}-\sin\theta$ plane, for the baryon density of
$n_{B}=n_{0}$ (left) and $5n_{0}$ (right), and temperatures $T=10$ MeV
(upper), 40 MeV (middle) and 100 MeV (bottom). The thermal equilibration times
are only valid for parameters where the scalars are trapped in the matter in
the merger (shown in grey). The other labels are the same as in Figure 4.

The thermal equilibration time due to a Bose gas of scalars in the
$m_{S}-\sin{\theta}$ plane is shown in Figure 9. As in Figure 4, the left and
right panels are respectively for the baryon density $n_{B}=n_{0}$ and
$5n_{0}$, and the upper, middle and lower panels correspond to the
temperatures of $T=10$ MeV, 40 MeV and 100 MeV, respectively. The labels for
all the current laboratory constraints and future prospects, the supernova
limits and the trapped and free streaming regions are the same as in Figure 4.
The thermal equilibration time depends on the mixing angle via
$\tau_{\kappa}\propto\sin^{2}{\theta}$, meaning that as the coupling between
nucleons and scalars grows weaker, the MFP of the scalar grows longer,
shortening the timescale needed to even out temperature gradients in the
remnant (provided that the scalar MFP remains much shorter than the system
size).

Figure 10: Thermal conductivity ratio $\kappa_{S}/\kappa_{\nu}$ in the
$m_{S}-\sin\theta$ plane, for the baryon density of $n_{B}=n_{0}$ (left) and
$5n_{0}$ (right), and temperatures $T=10$ MeV (upper), 40 MeV (middle) and 100
MeV (bottom). The labels are the same as in Figures 4 and 9 with the current
laboratory limits combined together as the black-shaded region.

While a trapped Bose gas of scalars can quickly even out temperature gradients
in the merger remnant as shown in Figure 9, a clear signature of BSM physics
is only present if that rapid thermal conduction was unable to be attributed
to SM particles, namely trapped neutrinos [194]. In Figure 10 we plot the
ratio of the thermal conductivity $\kappa_{S}$ due to trapped scalars to that
of trapped neutrinos $\kappa_{\nu}$, with baryon density $n_{B}=n_{0}$ and
$5n_{0}$ respectively in the left and right panels. For the sake of clarity,
we have combined all the current laboratory limits, and show them in black in
Figure 10. All other labels are the same as in Figures 4 and 9. In places
where the ratio $\kappa_{S}/\kappa_{\nu}$ is larger than one, fast thermal
equilibration is a signature of BSM physics, as no SM mechanism yields such a
high thermal conductivity. Figure 10 shows that at moderate temperatures like
$T=10\text{ MeV}$ (upper panels), the neutrino thermal conductivity is fairly
high, since the neutrino MFP is long (neutrinos become free-streaming at
temperatures just a few MeV lower) and so any contribution to the thermal
conductivity from the scalars is hidden by the neutrino contribution. One has
to go to higher temperatures where the neutrino MFP shrinks to find situations
where the scalars provide the dominant contribution to thermal equilibration
inside the merger remnant. In regions of the merger that reach $T=40\text{
MeV}$ (middle panels), the scalars can quickly equilibrate temperature
gradients, sometimes even one or two orders of magnitude faster than
neutrinos. This is only possible for scalars with masses below about 400 MeV,
and couplings $10^{-6}\lesssim\sin{\theta}\lesssim 10^{-4}$. If regions of the
merger remnant were to reach temperatures of 100 MeV (lower panels), trapped
scalars could quickly bring those regions to thermal equilibrium with nearby
regions for a wide range of potential scalar masses. It is remarkable that
some of the regions with $\kappa_{S}/\kappa_{\nu}>1$ can be directly probed at
the future laboratory experiments such as SBN, DUNE and FASER 2.

##  7 Discussion and conclusions

In this work, we have studied the effect of a CP-even scalar on the thermal
transport properties of a neutron star merger remnant. These scalars are
produced in the hot, dense nuclear matter in the merger by nucleon
bremsstrahlung processes. The effect of the scalar depends on its MFP in the
merger environment, which we plot in Figures 3 and 4. The parameter space of
$m_{S}-\sin{\theta}$ plane excluded by current laboratory limits and the
supernova limits from SN1987A are also depicted in Figure 4. It is clear that
taking into account the current constraints on the model, the MFP of the
scalars could be as short as hundreds of micrometers. On the other hand, as
the coupling between the scalar and nucleons becomes weaker, the MFP is free
to grow arbitrarily large.

If the MFP is longer than the size of a particular macroscopic region of the
merger remnant, the scalar will free-stream out of the region, cooling it. If
the scalar fails to encounter any region of the star where it has a short MFP
throughout its trajectory, it will escape the remnant entirely, leading to a
net loss of energy from the merger remnant as a whole. In Figures 5 and 6, we
plot the time for a fluid element to cool to half of its original temperature
by radiating scalar particles. As long as the scalar is not too heavy, it can
provide a mechanism to cool the hottest parts of the merger remnant on
timescales shorter than the remnant lifetime. While it does not seem that this
cooling affects the gravitational wave signal or the amount of ejecta
significantly [25], the emission of scalars and their subsequent decay might
be able to produce an additional electromagnetic signal, or the cooling to
trigger a phase transition or affect the transport properties in the nuclear
matter, which requires further investigation.

Perhaps more intriguingly, if the scalar is trapped in the merger remnant due
to a short MFP, forming a Bose-Einstein distribution, it could contribute to
transport processes in its own right. To this effect, we studied the
contribution to the thermal conductivity coming from the Bose gas of CP-even
scalar particles. The timescale for the scalar to smooth temperature gradients
in the merger remnant is plotted in Figures 7, 8, and 9. We find that there is
a significant amount of unconstrained space in the $m_{S}-\sin{\theta}$ plane
where the scalars can cause fast thermal equilibration. However, Figure 10
shows that only some of that parameter space represents a case where the
scalar thermal conductivity dominates over that coming from trapped neutrinos.
In particular, only at temperatures above 10 MeV can the scalar dominate the
thermal equilibration of the merger remnant. This could potentially lead to
some detectable signatures in future postmerger data. Furthermore, a part of
this parameter space can also be directly probed at the future laboratory
experiments such as SBN, DUNE and FASER 2 (cf. the middle and lower panels of
Figure 10), thus providing a complementary and independent probe of the model.

Before closing, we would like to comment on potential uncertainties and
possible future improvements. First of all, in the nucleon bremsstrahlung
process that produces $S$ particles, the nucleon scattering occurs via the
strong interaction. In our calculation of the matrix element for this process,
following Ref. [106] (and Refs. [201, 202, 203, 204, 205] for bremsstrahlung
processes producing other BSM particles), we used the OPE approximation. It is
well known that the OPE approximation underestimates the nucleon scattering
cross section at low energies and overestimates it at high energies [206]. In
the future, this calculation of the $S$ production rate should be conducted in
the soft-radiation approximation (SRA), which allows the $S$ production rate
to be written in terms of the on-shell nucleon-nucleon scattering amplitude
[207, 206, 8]. For other bremsstrahlung processes, the SRA predicts a
reduction of the production rate by a factor of a few compared to the rate
calculated in the OPE approximation. An alternative approach is to consider
the exchange of heavier mesons or multiple pions, which also reduces the
production rate relative to the OPE result [208]. If an improved treatment of
the nuclear interaction were to decrease the emissivity of scalars compared to
the OPE result we use here by, say, an order of magnitude (it would also
increase the MFP by a similar factor due to the similarity of the phase space
integrals for $Q$ and $\lambda$), our qualitative results would be unchanged,
though the specific parameter space where, for example, trapped scalars would
quickly smooth out temperature gradients in the system, would be different.

Secondly, scalars with masses $m_{S}>2m_{\pi^{0}}$ ($m_{S}>2m_{\pi^{\pm}}$)
can decay into neutral pions (charged pions). For even higher mass scalars,
decays into more than two pions or to heavier mesons are possible. Nuclear
matter in neutron star mergers likely also contains pions, either in a thermal
population or in a condensate, either of which would alter the decay
properties or even production mechanisms of a scalar particle. In this work we
allowed for a thermal population of pions at low density, which leads to Bose
enhancement of the decay process $S\rightarrow\pi^{+}\pi^{-}$. We treated the
pions as free particles, although the inclusion of pion-nucleon interactions
has been studied [178]. At high densities, the pions likely condense, but no
calculation has yet been done of the decay rate of the scalar to pions in the
presence of a pion condensate, nor is there a consensus on the density at
which condensation occurs. A unified, finite-temperature framework that can
describe pions in and out of the condensed phase is needed.

Ultimately, to use neutron star mergers to find signatures of scalar
particles, or to constrain them, will require numerical simulations which
incorporate scalar particles. In the free-streaming regime, this can be done
by a leakage scheme where each fluid element has an energy loss term
corresponding to the emissivity of scalars at the density and temperature of
the fluid element. This procedure was used to study cooling due to axion
radiation [25]. In the regime where scalars are trapped, a thermal conduction
term could be added to the relativistic hydrodynamic equations. However, it is
likely that for a scalar of a particular mass $m_{S}$ and coupling
$\sin{\theta}$, there will be thermodynamic conditions in certain parts of the
remnant where it will be trapped, but other parts of the remnant where it will
free-stream. Therefore, it would be necessary to develop a transport scheme
for scalars similar to those developed for neutrinos in merger or supernova
simulations [209, 210, 211, 164]. We hope the results presented here will
motivate the merger simulation community to implement the effects of light BSM
particles like the one studied here into their simulations.

## Acknowledgments

We thank Mark Alford for discussions on thermal conductivity; Albino Perego
for discussions on lepton number in neutron star merger simulations; Katy
Clough and Tim Dietrich for discussions on the possibility of including
scalars in merger simulations; and Shyam Balaji, Bhaskar Dutta, Kevin Kelly,
Rabi Mohapatra and Joe Silk for discussions on laboratory and astrophysical
constraints on light scalars. The work of BD is supported in part by the U.S.
Department of Energy under Grant No. DE-SC0017987 and by a Fermilab Intensity
Frontier Fellowship. This work was partly performed at the Aspen Center for
Physics, which is supported by National Science Foundation grant PHY-1607611.
The work of JFF is supported by NSERC. The work of SPH is supported by the U.
S. Department of Energy grant DE-FG02-00ER41132 as well as the National
Science Foundation grant No. PHY-1430152 (JINA Center for the Evolution of the
Elements). The work of KS is supported in part by DOE Grant DE-SC0009956. The
work of YZ is supported by the National Natural Science Foundation of China
under Grant No. 12175039, the 2021 Jiangsu Shuangchuang (Mass Innovation and
Entrepreneurship) Talent Program No. JSSCBS20210144, and the “Fundamental
Research Funds for the Central Universities”.

##  Appendix A Mean free path for $S$ decay processes

In this appendix we list the MFPs for the $S$ particle through its decays with
respect to the relevant SM particles. For convenience, we define for particle
species $i$

$\beta_{i}\equiv\left(1-\frac{m_{S}^{2}}{E_{S}^{2}}\right)\left(1-\frac{1}{\tau_{i}}\right),$
(A.1)

where $\tau_{i}\equiv m_{S}^{2}/4m_{i}^{2}$. The matrix element for
$S\rightarrow e^{+}e^{-}$ is

$\sum_{\text{spins}}|\mathcal{M}|_{e^{+}e^{-}}^{2}=2\sqrt{2}G_{F}\sin^{2}{\theta}m_{e}^{2}(m_{S}^{2}-4m_{e}^{2}),$
(A.2)

and the MFP of the scalar from this decay is

$\displaystyle\lambda^{-1}_{e^{+}e^{-}}$
$\displaystyle=\int\frac{\mathop{d^{3}{\bf
p}_{+}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{-}}}{(2\pi)^{3}}(2\pi)^{4}\delta^{4}(p_{S}-p_{+}-p_{-})\frac{\sum_{\text{spins}}|\mathcal{M}|_{e^{+}e^{-}}^{2}}{2^{3}E_{s}E_{+}E_{-}}$
$\displaystyle\quad\times(1-f(E_{+}+\mu_{e}))(1-f(E_{-}-\mu_{e}))$
$\displaystyle=\frac{\sin^{2}{\theta}}{8\pi v_{\rm
EW}^{2}}\frac{m_{e}^{2}(m_{S}^{2}-4m_{e}^{2})}{E_{S}\sqrt{E_{S}^{2}-m_{S}^{2}}}\frac{T}{1-e^{-E_{S}/T}}\Theta(m_{S}-2m_{e})$
$\displaystyle\quad\times\ln{\left\\{\frac{\cosh{(\mu_{e}/T)}+\cosh{[\frac{E_{S}}{2T}(1+\sqrt{\beta_{e}})]}}{\cosh{(\mu_{e}/T)}+\cosh{[\frac{E_{S}}{2T}(1-\sqrt{\beta_{e}})]}}\right\\}},$
(A.3)

where $p_{\pm}$ and $E_{\pm}$ are respectively the 4-momenta and energies of
$e^{\pm}$, and $f$ denotes a Fermi-Dirac factor. The MFP due to the decay
$S\rightarrow\mu^{+}\mu^{-}$ is the same, but with the electron mass and
chemical potential replaced by those of the muon. The matrix element of the
$S\rightarrow\gamma\gamma$ process is

$\sum_{\text{spins}}|\mathcal{M}|_{\gamma\gamma}^{2}=\frac{\sqrt{2}G_{F}\alpha^{2}m_{S}^{4}\sin^{2}{\theta}}{16\pi^{2}}\left|\sum_{f}N_{C}^{f}Q_{f}^{2}A_{1/2}(\tau_{f})+A_{1}(\tau_{W})\right|^{2}\,,$
(A.4)

where the one-loop functions are

$\displaystyle A_{1/2}(\tau)$ $\displaystyle\equiv$ $\displaystyle
2\left[\tau+(\tau-1)f(\tau)\right]\tau^{-2}\,,$ (A.5) $\displaystyle
A_{1}(\tau)$ $\displaystyle\equiv$
$\displaystyle-\left[2\tau^{2}+3\tau+3(2\tau-1)f(\tau)\right]\tau^{-2}\,,$
(A.6) $\displaystyle{\rm with}\quad f(\tau)$ $\displaystyle\equiv$
$\displaystyle\left\\{\begin{array}[]{cl}{\rm\arcsin}^{2}\sqrt{\tau}&({\rm
for}\leavevmode\nobreak\ \tau\leq 1)\\\
-{\displaystyle\frac{1}{4}}\left[\log\left(\frac{1+\sqrt{1-1/\tau}}{1-\sqrt{1-1/\tau}}\right)-i\pi\right]^{2}&({\rm
for}\leavevmode\nobreak\ \tau>1)\end{array}\right..$ (A.9)

The corresponding MFP is given by

$\displaystyle\lambda^{-1}_{\gamma\gamma}=$
$\displaystyle\int\frac{\mathop{d^{3}{\bf
p}_{1}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{2}}}{(2\pi)^{3}}(2\pi)^{4}\delta^{4}(p_{S}-p_{1}-p_{2})\frac{\sum_{\text{spins}}|\mathcal{M}|_{\gamma\gamma}^{2}}{2^{3}E_{s}E_{1}E_{2}}(1+g(E_{1}))(1+g(E_{2}))$
$\displaystyle=$
$\displaystyle\frac{\alpha^{2}\sin^{2}{\theta}}{128\pi^{3}v_{\rm
EW}^{2}}\frac{m_{S}^{4}}{E_{S}\sqrt{E_{S}^{2}-m_{S}^{2}}}\left|\sum_{f}N_{C}^{f}Q_{f}^{2}A_{1/2}(\tau_{f})+A_{1}(\tau_{W})\right|^{2}$
$\displaystyle\times\frac{T}{1-e^{-E_{S}/T}}\ln{\left\\{\frac{\sinh{\left[\frac{E_{S}}{4T}(1+\sqrt{\beta_{\gamma}})\right]}}{\sinh{\left[\frac{E_{S}}{4T}(1-\sqrt{\beta_{\gamma}})\right]}}\right\\}},$
(A.10)

where $p_{1,2}$ and $E_{1,2}$ are respectively the 4-momenta and energies of
the two photons, and $g$ is the Bose distribution factor.

We also consider the possibility of $S$ decaying into pions. We will treat the
pions as a free Bose gas, neglecting for simplicity modifications to their
dispersion relation due to the strong interaction (see Ref. [178]).
Furthermore, the strong interaction would likely alter the pion mass [212],
but we do not consider this effect, as pions have not yet been consistently
included in a RMF theory. The neutral pion has zero chemical potential, while
the charged pions have chemical potential
$\mu_{\pi^{\mp}}=\pm(\mu_{n}-\mu_{p})$. The matrix element for
$S\rightarrow\pi^{0}\pi^{0}$ is

$\sum_{\text{spins}}|\mathcal{M}|_{\pi^{0}\pi^{0}}^{2}=\frac{2}{81}\frac{\sin^{2}{\theta}}{v_{\rm
EW}^{2}}\left(m_{S}^{2}+\frac{11}{2}m_{\pi^{0}}^{2}\right)^{2},$ (A.11)

and the corresponding MFP is

$\displaystyle\lambda_{\pi^{0}\pi^{0}}^{-1}=$
$\displaystyle\int\frac{\mathop{d^{3}{\bf
p}_{1}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{2}}}{(2\pi)^{3}}(2\pi)^{4}\delta^{4}(p_{S}-p_{1}-p_{2})\frac{\sum_{\text{spins}}|\mathcal{M}|_{\pi^{0}\pi^{0}}^{2}}{2^{3}E_{S}E_{1}E_{2}}(1+g(E_{1}))(1+g(E_{2}))$
$\displaystyle=$ $\displaystyle\frac{\sin^{2}{\theta}}{324\pi v_{\rm
EW}^{2}}\frac{\left(m_{S}^{2}+\frac{11}{2}m_{\pi^{0}}^{2}\right)^{2}}{E_{S}\sqrt{E_{S}^{2}-m_{S}^{2}}}\frac{T}{1-e^{-E_{S}/T}}\ln{\left\\{\frac{\sinh{\left[\frac{E_{S}}{4T}(1+\sqrt{\beta_{\pi^{0}}})\right]}}{\sinh{\left[\frac{E_{S}}{4T}(1-\sqrt{\beta_{\pi^{0}}})\right]}}\right\\}}\Theta(m_{S}-2m_{\pi^{0}})\,,$
(A.12)

where $p_{1,2}$ and $E_{1,2}$ respectively the 4-momenta and energies for the
two neutral pions. The matrix element for $S\rightarrow\pi^{+}\pi^{-}$ is

$\sum_{\text{spins}}|\mathcal{M}|_{\pi^{+}\pi^{-}}^{2}=\frac{4}{81}\frac{\sin^{2}{\theta}}{v_{\rm
EW}^{2}}\left(m_{S}^{2}+\frac{11}{2}m_{\pi^{\pm}}^{2}\right)^{2},$ (A.13)

and the corresponding MFP is

$\displaystyle\lambda_{\pi^{+}\pi^{-}}^{-1}=$
$\displaystyle\int\frac{\mathop{d^{3}{\bf
p}_{+}}}{(2\pi)^{3}}\frac{\mathop{d^{3}{\bf
p}_{-}}}{(2\pi)^{3}}(2\pi)^{4}\delta^{4}(p_{S}-p_{+}-p_{-})\frac{\sum_{\text{spins}}|\mathcal{M}|_{\pi^{+}\pi^{-}}^{2}}{2^{3}E_{s}E_{+}E_{-}}$
$\displaystyle\quad\times(1+g(E_{+}+\mu_{\pi^{-}}))(1+g(E_{-}-\mu_{\pi^{-}}))$
$\displaystyle=$ $\displaystyle\frac{\sin^{2}{\theta}}{324\pi v_{\rm
EW}^{2}}\frac{\left(m_{S}^{2}+\frac{11}{2}m_{\pi^{\pm}}^{2}\right)^{2}}{E_{S}\sqrt{E_{S}^{2}-m_{S}^{2}}}\frac{T}{1-e^{-E_{S}/T}}\Theta(m_{S}-2m_{\pi^{\pm}})$
$\displaystyle\times\ln{\left\\{\frac{\cosh{(\mu_{\pi^{-}}/T)}-\cosh{\left[\frac{E_{S}}{2T}(1+\sqrt{\beta_{\pi^{\pm}}})\right]}}{\cosh{(\mu_{\pi^{-}}/T)}-\cosh{\left[\frac{E_{S}}{2T}(1-\sqrt{\beta_{\pi^{\pm}}})\right]}}\right\\}}\,,$
(A.14)

where $p_{\pm}$ and $E_{\pm}$ respectively the 4-momenta and energies for
$\pi^{\pm}$. This formula for scalar decaying into charged pions is only valid
when the pions have not formed a Bose-Einstein condensate.

##  Appendix B Dimensionless functions in nucleon bremsstrahlung matrix
element

The dimensionless functions in the nucleon bremsstrahlung matrix elements were
derived in Ref. [106] and are reproduced here:

$\displaystyle I_{A}^{nn/pp}$
$\displaystyle=\frac{\mathbf{k}^{4}}{(\mathbf{k}^{2}+m_{\pi}^{2})^{2}}+\frac{\mathbf{l}^{4}}{(\mathbf{l}^{2}+m_{\pi}^{2})^{2}}-\frac{\mathbf{k}^{2}\mathbf{l}^{2}-2(\mathbf{k}\cdot\mathbf{l})^{2}}{(\mathbf{k}^{2}+m_{\pi}^{2})(\mathbf{l}^{2}+m_{\pi}^{2})},$
(B.1) $\displaystyle I_{B}^{nn/pp}$
$\displaystyle=\left(m_{S}^{2}+\frac{11}{2}m_{\pi}^{2}\right)^{2}\left[\frac{\mathbf{k}^{4}}{(\mathbf{k}^{2}+m_{\pi}^{2})^{4}}+\frac{\mathbf{l}^{4}}{(\mathbf{l}^{2}+m_{\pi}^{2})^{4}}-\frac{\mathbf{k}^{2}\mathbf{l}^{2}-2(\mathbf{k}\cdot\mathbf{l})^{2}}{(\mathbf{k}^{2}+m_{\pi}^{2})^{2}(\mathbf{l}^{2}+m_{\pi}^{2})^{2}}\right],$
(B.2) $\displaystyle I_{C}^{nn/pp}$
$\displaystyle=\left(m_{S}^{2}+\frac{11}{2}m_{\pi}^{2}\right)\bigg{\\{}\frac{\mathbf{k}^{4}}{(\mathbf{k}^{2}+m_{\pi}^{2})^{3}}+\frac{\mathbf{l}^{4}}{(\mathbf{l}^{2}+m_{\pi}^{2})^{3}}$
$\displaystyle\hskip
142.26378pt-\frac{1}{2}\frac{(\mathbf{k}^{2}+\mathbf{l}^{2}+2m_{\pi}^{2})[\mathbf{k}^{2}\mathbf{l}^{2}-2(\mathbf{k}\cdot\mathbf{l})^{2}]}{(\mathbf{k}^{2}+m_{\pi}^{2})^{2}(\mathbf{l}^{2}+m_{\pi}^{2})^{2}}\bigg{\\}},$
(B.3) $\displaystyle I_{A}^{np}$
$\displaystyle=\frac{\mathbf{k}^{4}}{(\mathbf{k}^{2}+m_{\pi}^{2})^{2}}+\frac{4\mathbf{l}^{4}}{(\mathbf{l}^{2}+m_{\pi}^{2})^{2}}+2\frac{\mathbf{k}^{2}\mathbf{l}^{2}-2(\mathbf{k}\cdot\mathbf{l})^{2}}{(\mathbf{k}^{2}+m_{\pi}^{2})(\mathbf{l}^{2}+m_{\pi}^{2})},$
(B.4) $\displaystyle I_{B}^{np}$
$\displaystyle=\left(m_{S}^{2}+\frac{11}{2}m_{\pi}^{2}\right)^{2}\left[\frac{\mathbf{k}^{4}}{(\mathbf{k}^{2}+m_{\pi}^{2})^{4}}+\frac{4\mathbf{l}^{4}}{(\mathbf{l}^{2}+m_{\pi}^{2})^{4}}+2\frac{\mathbf{k}^{2}\mathbf{l}^{2}-2(\mathbf{k}\cdot\mathbf{l})^{2}}{(\mathbf{k}^{2}+m_{\pi}^{2})^{2}(\mathbf{l}^{2}+m_{\pi}^{2})^{2}}\right],$
(B.5) $\displaystyle I_{C}^{np}$
$\displaystyle=\left(m_{S}^{2}+\frac{11}{2}m_{\pi}^{2}\right)\bigg{\\{}\frac{\mathbf{k}^{4}}{(\mathbf{k}^{2}+m_{\pi}^{2})^{3}}+\frac{4\mathbf{l}^{4}}{(\mathbf{l}^{2}+m_{\pi}^{2})^{3}}$
$\displaystyle\hskip
142.26378pt+\frac{(\mathbf{k}^{2}+\mathbf{l}^{2}+2m_{\pi}^{2})[\mathbf{k}^{2}\mathbf{l}^{2}-2(\mathbf{k}\cdot\mathbf{l})^{2}]}{(\mathbf{k}^{2}+m_{\pi}^{2})^{2}(\mathbf{l}^{2}+m_{\pi}^{2})^{2}}\bigg{\\}},$
(B.6)

where $\mathbf{k}=\mathbf{p}_{2}-\mathbf{p}_{4}$ and
$\mathbf{l}=\mathbf{p}_{2}-\mathbf{p}_{3}$ measure the nucleon 3-momentum
transfers. We do not choose to average over the direction of the momentum of
the emitted BSM particle, as was done for axion emission, for example, in Ref.
[213].

##  Appendix C MFP phase space integral

In this appendix, we detail the calculation of the MFP of the $S$ particle due
to absorption via $n+n+S\rightarrow n+n$, and then comment on the other two
inverse bremsstrahlung absorption processes. The phase space integral in Eq.
(4.1) can be reduced to a five-dimensional integral via the following steps.
First, we neglect the 3-momentum of $S$ in the momentum-conserving
$\delta$-function; however, we keep its energy in the energy
$\delta$-function. We perform the 3-dimensional integral over $\mathbf{p}_{1}$
using the 3-momentum $\delta$-function and then we convert each of the
remaining 3-momentum integrals to spherical coordinates. We choose our
$z$-axis to align with $\mathbf{p}_{2}$, and then choose the $xz$-plane to be
the plane in which $\mathbf{p}_{2}$ and $\mathbf{p}_{3}$ lie, with angle
$r\equiv\cos{\theta}$ between them. This choice of axes allows us to write
$\mathbf{p}_{2}$, $\mathbf{p}_{3}$, and $\mathbf{p}_{4}$ as

$\displaystyle\mathbf{p}_{2}$ $\displaystyle=p_{2}(0,0,1),$
$\displaystyle\mathbf{p}_{3}$ $\displaystyle=p_{3}(\sqrt{1-r^{2}},0,r),$
$\displaystyle\mathbf{p}_{4}$
$\displaystyle=p_{4}(\sqrt{1-s^{2}}\cos{\phi},\sqrt{1-s^{2}}\sin{\phi},s),$
(C.1)

where $s$ is the cosine of the angle between $\mathbf{p}_{2}$ and
$\mathbf{p}_{4}$. This choice of coordinate system renders three of the
angular integrals trivial, yielding a factor of $8\pi^{2}$. The MFP is now
given by

$\displaystyle\lambda_{nn}^{-1}=$
$\displaystyle\frac{\alpha_{\pi}^{2}f^{4}\sin^{2}{\theta}}{16\pi^{4}m_{*}^{2}E_{S}}\int_{0}^{\infty}\mathop{dp_{2}}\mathop{dp_{3}}\mathop{dp_{4}}\int_{-1}^{1}\mathop{dr}\mathop{ds}\int_{0}^{2\pi}\mathop{d\phi}p_{2}^{2}p_{3}^{2}p_{4}^{2}\frac{f_{1}f_{2}(1-f_{3})(1-f_{4})}{E_{1}^{*}E_{2}^{*}E_{3}^{*}E_{4}^{*}}$
$\displaystyle\times\delta(E_{1}^{*}+E_{2}^{*}-E_{3}^{*}-E_{4}^{*}+E_{S})\left(\frac{y_{hnn}^{2}m_{S}^{4}}{4E_{S}^{4}}I_{A}^{nn}+\frac{m_{*}^{2}}{192v_{\rm
EW}^{2}}I_{B}^{nn}+\frac{2y_{hnn}m_{S}^{2}m_{*}}{9E_{S}^{2}v_{\rm
EW}}I_{C}^{nn}\right)\,,$ (C.2)

where
$E_{1}^{*}=\sqrt{(\mathbf{p}_{3}+\mathbf{p}_{4}-\mathbf{p}_{2})^{2}+m_{*}^{2}}$.
Now, we integrate over $\phi$, using the energy $\delta$-function. The
argument of the $\delta$-function can only go through zero if

$E_{3}^{*}+E_{4}^{*}>E_{2}^{*}+E_{S}.$ (C.3)

If

$\displaystyle|m_{*}^{2}+E_{S}^{2}/2+E_{3}^{*}E_{4}^{*}-E_{2}^{*}E_{3}^{*}-$
$\displaystyle
E_{2}^{*}E_{4}^{*}+E_{S}(E_{2}^{*}-E_{3}^{*}-E_{4}^{*})+p_{2}p_{3}r+p_{2}p_{4}s-p_{3}p_{4}rs|$
$\displaystyle\leq p_{3}p_{4}\sqrt{1-r^{2}}\sqrt{1-s^{2}},$ (C.4)

then there are two values of $\phi$ between $0$ and $2\pi$ which are picked
out by the $\delta$-function. Both values $\phi_{0}$ have the same cosine. If
Eq. (C.4) does not hold, then the integral is zero. Condition (C.3) can be
turned into an upper bound on $p_{2}$ along with the further constraint

$E_{3}^{*}+E_{4}^{*}-E_{S}\geq m_{*}.$ (C.5)

Therefore the MFP of the scalar $S$ from the neutron-neutron inverse
bremsstrahlung process can be written as

$\displaystyle\lambda_{nn}^{-1}=$
$\displaystyle\frac{\alpha_{\pi}^{2}f^{4}\sin^{2}{\theta}}{8\pi^{4}m_{*}^{2}E_{S}}\int_{0}^{\infty}\mathop{dp_{3}}\mathop{dp_{4}}\int_{0}^{\sqrt{(E_{3}^{*}+E_{4}^{*}-E_{S})^{2}-m_{*}^{2}}}\mathop{dp_{2}}\int_{-1}^{1}\mathop{dr}\mathop{ds}\frac{p_{2}^{2}p_{3}^{2}p_{4}^{2}}{E_{2}^{*}E_{3}^{*}E_{4}^{*}}f_{1}f_{2}(1-f_{3})(1-f_{4})$
$\displaystyle\times\left(\frac{y_{hnn}^{2}m_{S}^{4}}{4E_{S}^{4}}I_{A}^{nn}+\frac{m_{*}^{2}}{192v_{\rm
EW}^{2}}I_{B}^{nn}+\frac{2y_{hnn}m_{S}^{2}m_{*}}{9E_{S}^{2}v_{\rm
EW}}I_{C}^{nn}\right)$
$\displaystyle\times\bigg{(}p_{3}^{2}p_{4}^{2}(1-r^{2})(1-s^{2})-(m_{*}^{2}+E_{S}^{2}/2+E_{3}^{*}E_{4}^{*}-E_{2}^{*}E_{3}^{*}-E_{2}^{*}E_{4}^{*}$
$\displaystyle+E_{S}(E_{2}^{*}-E_{3}^{*}-E_{4}^{*})+p_{2}p_{3}r+p_{2}p_{4}s-p_{3}p_{4}rs)^{2}\bigg{)}^{-1/2}\,,$
(C.6)

with the restrictions (C.4) and (C.5). The Fermi-Dirac factors are

$\displaystyle f_{1}$
$\displaystyle=\frac{1}{1+e^{(E_{3}^{*}+E_{4}^{*}-E_{2}^{*}-E_{S}-\mu_{n}^{*})/T}}\,,$
$\displaystyle f_{2}$
$\displaystyle=\frac{1}{1+e^{(E_{2}^{*}-\mu_{n}^{*})/T}}\,,$ $\displaystyle
1-f_{3}$ $\displaystyle=\frac{1}{1+e^{-(E_{3}^{*}-\mu_{n}^{*})/T}}\,,$
$\displaystyle 1-f_{4}$
$\displaystyle=\frac{1}{1+e^{-(E_{4}^{*}-\mu_{n}^{*})/T}}\,.$ (C.7)

The integration variables $p_{2},p_{3},p_{4}$ can be nondimensionalized by
$m_{*}$, if desired.

The MFP due to absorption via $p+p+S\rightarrow p+p$ is the same, except that
the neutron chemical potential is replaced by that of the proton. The MFP due
to the absorption $n+p+S\rightarrow n+p$ has the neutron chemical potentials
of particles 2 and 4 replaced by the proton chemical potential, and the result
is divided by four (or, $\lambda_{np}^{-1}$ is multiplied by a factor of four)
due to the difference in matrix elements (cf. Eqs. (4.2) and (4.3)).

##  Appendix D Emissivity phase space integral

To evaluate the phase space integral in Eq. (5.4) for, say, $n+n\rightarrow
n+n+S$, we neglect the 3-momentum of the scalar particle in the 3-momentum
conserving $\delta$-function and then, considering the integral over
$\mathbf{p}_{S}$ in spherical coordinates, integrate over the angular part,
giving a factor of $4\pi$. To simplify the 3-momentum conserving
$\delta$-function, we change the coordinate system from
$\\{\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4}\\}$ to
$\\{\mathbf{p}_{+},\mathbf{p}_{-},\mathbf{a},\mathbf{b}\\}$ where

$\displaystyle\mathbf{p}_{+}$
$\displaystyle=\frac{1}{2}(\mathbf{p}_{1}+\mathbf{p}_{2})\,,$
$\displaystyle\mathbf{a}$ $\displaystyle=\mathbf{p}_{3}-\mathbf{p}_{+}\,,$
$\displaystyle\mathbf{p}_{-}$
$\displaystyle=\frac{1}{2}(\mathbf{p}_{1}-\mathbf{p}_{2})\,,$
$\displaystyle\mathbf{b}$ $\displaystyle=\mathbf{p}_{4}-\mathbf{p}_{+}\,.$
(D.1)

The Jacobian gives a factor of 8, and the emissivity becomes

$\displaystyle Q_{nn}=$
$\displaystyle\frac{\alpha_{\pi}^{2}f^{4}\sin^{2}{\theta}}{32\pi^{8}m_{*}^{2}}\int\mathop{d^{3}{\bf
p}_{+}}\mathop{d^{3}{\bf p}_{-}}\mathop{d^{3}{\bf a}}\mathop{d^{3}{\bf
b}}\int_{0}^{\infty}\mathop{dp_{S}}p_{S}^{2}\delta(E_{1}^{*}+E_{2}^{*}-E_{3}^{*}-E_{4}^{*}-E_{S})\delta^{(3)}(\mathbf{a}+\mathbf{b})$
$\displaystyle\times\frac{f_{1}f_{2}(1-f_{3})(1-f_{4})}{E_{1}^{*}E_{2}^{*}E_{3}^{*}E_{4}^{*}}\left(\frac{y_{hnn}^{2}m_{S}^{4}}{E_{S}^{4}}I_{A}^{nn}+\frac{m_{*}^{2}}{81v_{\rm
EW}^{2}}I_{B}^{nn}+\frac{2y_{hnn}m_{S}^{2}m_{*}}{9E_{S}^{2}v_{\rm
EW}}I_{C}^{nn}\right).$ (D.2)

Now we do the integral over $\mathbf{b}$ using the 3-momentum
$\delta$-function and then the integral over $p_{S}$, making use of the energy
$\delta$-function. This yields

$\displaystyle Q_{nn}=$
$\displaystyle\frac{\alpha_{\pi}^{2}f^{4}\sin^{2}{\theta}}{32\pi^{8}m_{*}^{2}}\int\mathop{d^{3}{\bf
p}_{+}}\mathop{d^{3}{\bf p}_{-}}\mathop{d^{3}{\bf
a}}(E_{1}^{*}+E_{2}^{*}-E_{3}^{*}-E_{4}^{*})\sqrt{(E_{1}^{*}+E_{2}^{*}-E_{3}^{*}-E_{4}^{*})^{2}-m_{S}^{2}}$
$\displaystyle\times\frac{f_{1}f_{2}(1-f_{3})(1-f_{4})}{E_{1}^{*}E_{2}^{*}E_{3}^{*}E_{4}^{*}}\left(\frac{y_{hnn}^{2}m_{S}^{4}}{4E_{S}^{4}}I_{A}^{nn}+\frac{m_{*}^{2}}{192v_{\rm
EW}^{2}}I_{B}^{nn}+\frac{2y_{hnn}m_{S}^{2}m_{*}}{9E_{S}^{2}v_{\rm
EW}}I_{C}^{nn}\right),$ (D.3)

with the constraint

$E_{1}^{*}+E_{2}^{*}-E_{3}^{*}-E_{4}^{*}\geq m_{S}.$ (D.4)

We choose a coordinate system such that

$\displaystyle\mathbf{a}$ $\displaystyle=a(0,0,1),$
$\displaystyle\mathbf{p_{-}}$ $\displaystyle=p_{-}(\sqrt{1-r^{2}},0,r),$
$\displaystyle\mathbf{p_{+}}$
$\displaystyle=p_{+}(\sqrt{1-s^{2}}\cos{\phi},\sqrt{1-s^{2}}\sin{\phi},s),$
(D.5)

with $r$ and $s$ here respectively the cosines of the angles between ${\bf
p}_{\mp}$ and ${\bf a}$, and then do the integrals over three trivial angles,
yielding

$\displaystyle Q_{nn}=$
$\displaystyle\frac{\alpha_{\pi}^{2}f^{4}\sin^{2}{\theta}}{4\pi^{6}m_{*}^{2}}\int_{0}^{\infty}\mathop{dp_{+}}\mathop{dp_{-}}\mathop{da}\int_{-1}^{1}\mathop{dr}\mathop{ds}\int_{0}^{2\pi}\mathop{d\phi}p_{+}^{2}p_{-}^{2}a^{2}(E_{1}^{*}+E_{2}^{*}-E_{3}^{*}-E_{4}^{*})$
$\displaystyle\times\sqrt{(E_{1}^{*}+E_{2}^{*}-E_{3}^{*}-E_{4}^{*})^{2}-m_{S}^{2}}\frac{f_{1}f_{2}(1-f_{3})(1-f_{4})}{E_{1}^{*}E_{2}^{*}E_{3}^{*}E_{4}^{*}}$
$\displaystyle\times\left(\frac{y_{hnn}^{2}m_{S}^{4}}{4E_{S}^{4}}I_{A}^{nn}+\frac{m_{*}^{2}}{192v_{\rm
EW}^{2}}I_{B}^{nn}+\frac{2y_{hnn}m_{S}^{2}m_{*}}{9E_{S}^{2}v_{\rm
EW}}I_{C}^{nn}\right),$ (D.6)

where the condition (D.4) holds, and
$E_{S}=E_{1}^{*}+E_{2}^{*}-E_{3}^{*}-E_{4}^{*}$. In the functions $I_{A,B,C}$
(cf. Appendix B), the momentum transfer variables are

$\displaystyle\mathbf{k}^{2}=a^{2}+p_{-}^{2}-2p_{-}ar,\quad\mathbf{l}^{2}=a^{2}+p_{-}^{2}+2p_{-}ar,\quad\mathbf{k}\cdot\mathbf{l}=p_{-}^{2}-a^{2}.$
(D.7)

In addition, the former integration variables are now expressed as functions
of the new integration variables

$\displaystyle{\bf p}_{1}^{2}$
$\displaystyle=p_{-}^{2}+p_{-}^{2}+2p_{+}p_{-}rs+2p_{+}p_{-}\sqrt{1-r^{2}}\sqrt{1-s^{2}}\cos{\phi}\,,$
$\displaystyle{\bf p}_{2}^{2}$
$\displaystyle=p_{+}^{2}p_{-}^{2}-2p_{+}p_{-}rs-2p_{+}p_{-}\sqrt{1-r^{2}}\sqrt{1-s^{2}}\cos{\phi}\,,$
$\displaystyle{\bf p}_{3}^{2}$ $\displaystyle=p_{+}^{2}+a^{2}+2p_{+}as\,,$
$\displaystyle{\bf p}_{4}^{2}$ $\displaystyle=p_{+}^{2}+a^{2}-2p_{+}as\,.$
(D.8)

##  Appendix E Higher-order corrections for scalar production amplitude

We perform in this appendix a more complete calculation for the nucleon
bremsstrahlung process Eq. (2.1), following the procedure in Ref. [106]. In
this section, we use $m_{N}$ to label the nucleon mass, which in a RMF theory
is replaced with the Dirac effective mass $m_{*}$ (see Section 3). For the
diagrams (a) through (d) in Figure 1, the amplitudes are

$\displaystyle{\cal M}_{a,b,c,d}\propto\frac{1}{(p_{i}\mp
k_{s})^{2}-m_{N}^{2}}\,,$ (E.1)

where the “$-$” and “$+$” signs are respectively for $i=1,\,2$ and $i=3,\,4$,
and

$\displaystyle(p_{i}+k_{s})^{2}-m_{N}^{2}$ $\displaystyle=$ $\displaystyle\mp
2(p_{i}\cdot k_{s})+m_{S}^{2}$ (E.2) $\displaystyle=$ $\displaystyle\mp
2m_{N}E_{S}\left[1+\frac{E_{i}-m_{N}}{m_{N}}-\frac{{\bf p}_{i}\cdot{\bf
k}_{S}}{m_{N}E_{S}}\mp\frac{m_{S}^{2}}{2m_{N}E_{S}}\right]\,.$

For convenience, we define the small dimensionless parameters

$\displaystyle a_{i}\equiv\frac{E_{i}-m_{N}}{m_{N}}\,,\quad
b_{i}\equiv-\frac{{\bf p}_{i}\cdot{\bf k}_{S}}{m_{N}E_{S}}\,,\quad
c_{i}\equiv\mp c\,,$ (E.3)

where

$c\equiv\frac{m_{S}^{2}}{2m_{N}E_{S}}\,.$ (E.4)

Then the propagators become

$\displaystyle\frac{1}{(p_{i}\mp k_{s})^{2}-m_{N}^{2}}$ $\displaystyle=$
$\displaystyle\mp\frac{1}{2m_{N}E_{S}}\frac{1}{1+a_{i}+b_{i}+c_{i}}\,,$ (E.5)
$\displaystyle{\rm with}\quad\frac{1}{1+a_{i}+b_{i}+c_{i}}$ $\displaystyle=$
$\displaystyle
1+\sum_{k=1}^{\infty}(-1)^{k}\left(a_{i}^{k}+b_{i}^{k}+c_{i}^{k}\right)$ (E.6)
$\displaystyle+\sum_{k=2}^{\infty}\sum_{j=1}^{k-1}(-1)^{k}\frac{k!}{j!(k-j)!}\left(a_{i}^{j}b_{i}^{k-j}+a_{i}^{j}c_{i}^{k-j}+b_{i}^{j}c_{i}^{k-j}\right)$
$\displaystyle+\sum_{k=3}^{\infty}\sum_{j,m=1}^{k-1}(-1)^{k}\frac{k!}{j!m!(k-j-m)!}a_{i}^{j}b_{i}^{m}c_{i}^{k-j-m}\,.$

When we sum up the (a) through (d) diagrams, the “1” term in Eq. (E.6) cancels
out, which is the cancellation mentioned in Ref. [106]. For the $a_{i}$ terms,
let us first consider the leading order terms with $k=1$, which produce

$\displaystyle-\sum_{i}\mp a_{i}$ $\displaystyle=$
$\displaystyle\frac{E_{1}-m_{N}}{m_{N}}+\frac{E_{2}-m_{N}}{m_{N}}-\frac{E_{3}-m_{N}}{m_{N}}-\frac{E_{4}-m_{N}}{m_{N}}=a\,,$
(E.7) $\displaystyle-\sum_{i}\mp b_{i}$ $\displaystyle=$
$\displaystyle-\frac{{\bf p}_{1}\cdot{\bf k}_{S}}{m_{N}E_{S}}-\frac{{\bf
p}_{2}\cdot{\bf k}_{S}}{m_{N}E_{S}}+\frac{{\bf p}_{3}\cdot{\bf
k}_{S}}{m_{N}E_{S}}+\frac{{\bf p}_{4}\cdot{\bf k}_{S}}{m_{N}E_{S}}=-b\,,$
(E.8) $\displaystyle(-1)^{k}\sum_{i}\mp c_{i}^{k}$ $\displaystyle=$
$\displaystyle-2\left(1+(-1)^{k-1}\right)c^{k}\,,$ (E.9)

where we introduced

$\displaystyle a\equiv\frac{E_{S}}{m_{N}}\,,\quad
b\equiv\frac{E_{S}^{2}-m_{S}^{2}}{m_{N}E_{S}}\,.$ (E.10)

For the high-order terms $a_{i}^{k}$ with $k=2,3,4...$ in Eq. (E.6), in the
center-of-mass frame and in the limit of $m_{S}$, $E_{S}\to 0$, $\Delta
E_{i}=E_{i}-m_{N}$ are all the same; as a result the $(E_{i}-m_{N})^{k}$ terms
with $k\geq 2$ cancel out. In other words, these higher-order terms are highly
suppressed by powers of $(E_{S}/m_{N})^{k}$. For the higher-order terms
$b_{i}^{k}$ (with $k=2,3,4...$) in the center-of-mass frame and in the limit
of $m_{S},\,E_{S}\to 0$, the scalar products

$\displaystyle{\bf p}_{i}\cdot{\bf k}_{S}=|{\bf p}_{i}||{\bf
k}_{S}|\cos\theta_{i}\,,$ (E.11)

where $\theta_{i}$ are the angles between ${\bf p}_{i}$ and ${\bf k}_{S}$,
then

$\displaystyle\sum_{i}\mp b_{i}^{k}$ $\displaystyle\propto$
$\displaystyle\left(\frac{|{\bf p}_{i}||{\bf
k}_{S}|}{m_{N}E_{S}}\right)^{k}\left(\cos^{k}\theta_{1}+\cos^{k}\theta_{2}-\cos^{k}\theta_{3}-\cos^{k}\theta_{4}\right)\,.$
(E.12)

Integrating over all the parameter space, the angles $\theta_{1,2,3,4}$ are
averaged, thus the factor

$\displaystyle\left(\cos^{n}\theta_{1}+\cos^{n}\theta_{2}-\cos^{n}\theta_{3}-\cos^{n}\theta_{4}\right)\to
0\,.$ (E.13)

In other words, the higher-order terms of $b_{i}^{k}$ (with $k\geq 2$) are
also highly suppressed by $(E_{S}/m_{N})^{k}$.

For the cross terms $a_{i}^{j}b_{i}^{k-j}$ in Eq. (E.6), let us first consider
the first term with $j=1$ and $k=2$, which turns out to be

$\displaystyle\sum_{i}\mp 2a_{i}b_{i}$ $\displaystyle=$
$\displaystyle\frac{2}{m_{N}^{2}E_{S}}\left[E_{1}{\bf p}_{1}\cdot{\bf
k}_{S}+E_{2}{\bf p}_{2}\cdot{\bf k}_{S}-E_{3}{\bf p}_{3}\cdot{\bf
k}_{S}-E_{4}{\bf p}_{4}\cdot{\bf k}_{S}\right]-\frac{2{\bf
k}_{S}^{2}}{m_{N}E_{S}}\,.$ (E.14)

To evaluate the $E_{i}{\bf p}_{i}\cdot{\bf k}_{S}$ terms, we take the limit of
$E_{i}-m_{N}\to\Delta E\simeq{\bf p}_{N}^{\,2}/2m_{N}^{2}$ in the center-of-
mass frame with the limit of $m_{S},\,E_{S}\to 0$. In this limit, Eq. (E.14)
can be greatly simplified to

$\displaystyle\sum_{i}\mp 2a_{i}b_{i}$ $\displaystyle=$ $\displaystyle 2bd\,,$
(E.15)

where we have defined $d\equiv\Delta E/m_{N}$. Generalizing to higher-order
terms of the form $a_{i}^{k-1}b_{i}$ (with $k\geq 2$) in Eq. (E.6), we have

$\displaystyle(-1)^{k}\sum_{i}\mp ka_{i}^{k-1}b_{i}$ $\displaystyle=$
$\displaystyle(-1)^{k}kbd^{k-1}\,.$ (E.16)

The $a_{i}^{k-j}b_{i}^{j}$ terms (with $j\geq 2$) are vanishing as

$\displaystyle a_{i}^{m}b_{i}^{n}\propto\left(\frac{\Delta
E}{m_{N}}\right)^{m}\left(\cos^{n}\theta_{1}+\cos^{n}\theta_{2}-\cos^{n}\theta_{3}-\cos^{n}\theta_{4}\right)\,.$
(E.17)

Similarly it is straightforward to obtain that

$\displaystyle\sum_{i}\mp a_{i}^{j}c_{i}^{k-j}$ $\displaystyle=$
$\displaystyle\frac{1}{2}c^{k-j}\left[4d^{j}\left(1+(-1)^{k-j-1}\right)-\delta_{j1}a\left(1+(-1)^{k-j}\right)\right]\,,$
(E.18) $\displaystyle\sum_{i}\mp b_{i}^{j}c_{i}^{k-j}$ $\displaystyle=$
$\displaystyle\frac{1}{2}bc^{k-j}\delta_{j1}\left(1+(-1)^{k-j}\right)\,,$
(E.19) $\displaystyle\sum_{i}\mp a_{i}^{j}b_{i}^{m}c_{i}^{k-j-m}$
$\displaystyle=$
$\displaystyle\frac{1}{2}bd^{j}c^{k-j-m}\delta_{m1}\left(1+(-1)^{k-j-m}\right)\,.$
(E.20)

Summing up all the contributions, we get

$\displaystyle
a-b-4\sum_{n=1}^{\infty}c^{2n-1}+\sum_{n=1}^{\infty}(-1)^{n+1}(n+1)d^{n}b$
(E.21)
$\displaystyle+4\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}(-1)^{m+2n-1}C_{m+2n-1}^{m}d^{m}c^{2n-1}+\sum_{n=1}^{\infty}(1+2n)ac^{2n}$
$\displaystyle-\sum_{n=1}^{\infty}(1+2n)bc^{2n}+\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}(-1)^{m+2n+1}\frac{(m+2n+1)!}{m!(2n)!}d^{m}bc^{2n}$
$\displaystyle=$
$\displaystyle\frac{(1+c^{2})a}{(1-c^{2})^{2}}-\frac{[(1+d)^{2}+c^{2}]b}{[(1+d)^{2}-c^{2}]^{2}}-\frac{4c}{(1+d)^{2}-c^{2}}\,.$

As all the parameters $a$, $b$, $c$ and $d$ are very small, the summation can
be approximated to be

$\displaystyle
a-b-4c=\frac{E_{S}}{m_{N}}-\frac{E_{S}^{2}-m_{S}^{2}}{m_{N}E_{S}}-\frac{2m_{S}^{2}}{m_{N}E_{S}}=-\frac{m_{S}^{2}}{m_{N}E_{S}}\,,$
(E.22)

which corresponds to the lowest order terms $a_{i}$, $b_{i}$ and $c_{i}$ in
Eq. (E.6). This means that we have to scale the amplitudes for the (a) through
(d) diagrams in Ref. [106] by (roughly) a factor of 1/2.

##  Appendix F Neutrino thermal conductivity

The thermal conductivity of a neutrino gas has been calculated for one
neutrino species in the degenerate limit [214] and at arbitrary neutrino
degeneracy (though assuming nondegenerate nuclear matter) [215]. In mergers,
there exist neutrinos and antineutrinos of electron and muon flavors, so we
perform an estimate of the neutrino thermal conductivity using the simple
expression given by Eq. (6.6) that includes all four particle types and does
not make assumptions about the degeneracy of neutrinos. However, Eq. (6.6)
assumes that the MFP is independent of the neutrino energy, which is certainly
not true for neutrinos in dense matter (see, for example, Figure 1 in Ref.
[216] or Figures 5 and 6 in Ref. [217]), so the following calculation is just
a rough estimate.

The neutrino contribution to the thermal conductivity is the sum of the
individual thermal conductivities of electron and muon neutrinos and their
antineutrino counterparts. The specific heat of arbitrary degeneracy gases of
$\nu_{e}$ and $\bar{\nu}_{e}$ are

$\displaystyle c_{V,{\nu_{e}}}$
$\displaystyle=\frac{3T^{3}}{\pi^{2}}\left[\frac{\mu_{\nu_{e}}}{T}\phi_{3}(-e^{\mu_{\nu_{e}}/T})-4\phi_{4}(-e^{\mu_{\nu_{e}}/T})\right],$
(F.1) $\displaystyle c_{V,{\bar{\nu}_{e}}}$
$\displaystyle=-\frac{3T^{3}}{\pi^{2}}\left[\frac{\mu_{\nu_{e}}}{T}\phi_{3}(-e^{-\mu_{\nu_{e}}/T})+4\phi_{4}(-e^{-\mu_{\nu_{e}}/T})\right],$
(F.2)

where $\phi_{n}(x)$ is the polylogarithm function of order $n$. The sum of the
specific heat of electron neutrinos and antineutrinos is

$\displaystyle
c_{V,{\nu_{e}}}+c_{V,\bar{\nu}_{e}}=\frac{7}{30}\pi^{2}T^{3}+\frac{1}{2}\mu_{\nu_{e}}^{2}T\,,$
(F.3)

which is consistent with the expression given for the energy density of a
fermion-antifermion gas in Ref. [218]. But, we will treat neutrinos and
antineutrinos as separate species for this calculation because they have
different MFPs. The specific heat expressions for muon neutrinos and
antineutrinos are analogous to Eqs. (F.1) and (F.2).

The neutrino thermal conductivity is then given by

$\kappa_{\text{all
neutrinos}}=\frac{1}{3}\left(c_{V,{\nu_{e}}}\lambda_{\nu_{e}}+c_{V,\bar{\nu}_{e}}\lambda_{\bar{\nu}_{e}}+c_{V,{\nu_{\mu}}}\lambda_{\nu_{\mu}}+c_{V,\bar{\nu}_{\mu}}\lambda_{\bar{\nu}_{\mu}}\right),$
(F.4)

where the velocity term in the thermal conductivity expression (see Eq. (6.6))
has been set to one because the neutrinos are essentially massless. The
neutrino/antineutrino MFPs were calculated in the same manner as in Refs.
[216, 217] (that is, using nucleon dispersion relations obtained from a RMF
theory and conducting the full integration over phase space), but neglecting
for simplicity the weak magnetism and pseudoscalar parts of the hadronic
current. The following processes were included in the calculation:

$\displaystyle\nu_{e}+n\rightarrow p+e^{-},$
$\displaystyle\nu_{e}+n\rightarrow\nu_{e}+n,$
$\displaystyle\bar{\nu}_{e}+p\rightarrow e^{+}+n,$
$\displaystyle\bar{\nu}_{e}+p+e^{-}\rightarrow n,$
$\displaystyle\bar{\nu}_{e}+n\rightarrow\bar{\nu}_{e}+n,$
$\displaystyle\nu_{\mu}+n\rightarrow p+\mu^{-},$
$\displaystyle\nu_{\mu}+n\rightarrow\nu_{\mu}+n,$
$\displaystyle\bar{\nu}_{\mu}+p\rightarrow\mu^{+}+n,$
$\displaystyle\bar{\nu}_{\mu}+p+\mu^{-}\rightarrow n,$
$\displaystyle\bar{\nu}_{\mu}+n\rightarrow\bar{\nu}_{\mu}+n.$ (F.5)

We evaluate the neutrino MFP at $E_{\nu}\approx 3T$, which is its average
energy when the chemical potential is much less than the temperature [219]. A
proper calculation of the neutrino thermal conductivity would involve solving
the Boltzmann equation, properly accounting for the full distribution of
neutrinos and the energy-dependence of the neutrino MFP.

## References

  * [1] E. W. Kolb and M. S. Turner, The Early Universe.  CRC Press, 1990.
  * [2] G. G. Raffelt, Astrophysical methods to constrain axions and other novel particle phenomena, Phys. Rept. 198 (1990) 1–113.
  * [3] G. G. Raffelt, Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles.  5, 1996.
  * [4] R. A. Batista et al., EuCAPT White Paper: Opportunities and Challenges for Theoretical Astroparticle Physics in the Next Decade, 2110.10074.
  * [5] LIGO Scientific, Virgo collaboration, B. Abbott et al., GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119 (2017) 161101, [1710.05832].
  * [6] LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT collaboration, B. P. Abbott et al., Multi-messenger Observations of a Binary Neutron Star Merger, Astrophys. J. Lett. 848 (2017) L12, [1710.05833].
  * [7] LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL collaboration, B. P. Abbott et al., Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A, Astrophys. J. Lett. 848 (2017) L13, [1710.05834].
  * [8] J.-F. Fortin, H.-K. Guo, S. P. Harris, D. Kim, K. Sinha and C. Sun, Axions: From magnetars and neutron star mergers to beam dumps and BECs, Int. J. Mod. Phys. D 30 (2021) 2130002, [2102.12503].
  * [9] E. Burns, Neutron Star Mergers and How to Study Them, Living Rev. Rel. 23 (2020) 4, [1909.06085].
  * [10] M. Evans et al., A Horizon Study for Cosmic Explorer: Science, Observatories, and Community, 2109.09882.
  * [11] S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger and C. M. Biwer, Tidal Deformabilities and Radii of Neutron Stars from the Observation of GW170817, Phys. Rev. Lett. 121 (2018) 091102, [1804.08583]. [Erratum: Phys.Rev.Lett. 121, 259902 (2018)].
  * [12] LIGO Scientific, Virgo collaboration, B. P. Abbott et al., GW170817: Measurements of neutron star radii and equation of state, Phys. Rev. Lett. 121 (2018) 161101, [1805.11581].
  * [13] D. Kasen, B. Metzger, J. Barnes, E. Quataert and E. Ramirez-Ruiz, Origin of the heavy elements in binary neutron-star mergers from a gravitational wave event, Nature 551 (2017) 80, [1710.05463].
  * [14] LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER collaboration, B. P. Abbott et al., A gravitational-wave standard siren measurement of the Hubble constant, Nature 551 (2017) 85–88, [1710.05835].
  * [15] M. Bezares, D. Viganò and C. Palenzuela, Gravitational wave signatures of dark matter cores in binary neutron star mergers by using numerical simulations, Phys. Rev. D 100 (2019) 044049, [1905.08551].
  * [16] S. H. Alexander and N. Yunes, Gravitational wave probes of parity violation in compact binary coalescences, Phys. Rev. D 97 (2018) 064033, [1712.01853].
  * [17] P. Creminelli and F. Vernizzi, Dark Energy after GW170817 and GRB170817A, Phys. Rev. Lett. 119 (2017) 251302, [1710.05877].
  * [18] J. Sakstein and B. Jain, Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories, Phys. Rev. Lett. 119 (2017) 251303, [1710.05893].
  * [19] J. M. Ezquiaga and M. Zumalacárregui, Dark Energy After GW170817: Dead Ends and the Road Ahead, Phys. Rev. Lett. 119 (2017) 251304, [1710.05901].
  * [20] S. Boran, S. Desai, E. O. Kahya and R. P. Woodard, GW170817 Falsifies Dark Matter Emulators, Phys. Rev. D 97 (2018) 041501, [1710.06168].
  * [21] T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller and I. Sawicki, Strong constraints on cosmological gravity from GW170817 and GRB 170817A, Phys. Rev. Lett. 119 (2017) 251301, [1710.06394].
  * [22] N. Sennett, R. Brito, A. Buonanno, V. Gorbenko and L. Senatore, Gravitational-Wave Constraints on an Effective Field-Theory Extension of General Relativity, Phys. Rev. D 102 (2020) 044056, [1912.09917].
  * [23] A. Hook and J. Huang, Probing axions with neutron star inspirals and other stellar processes, JHEP 06 (2018) 036, [1708.08464].
  * [24] J. Huang, M. C. Johnson, L. Sagunski, M. Sakellariadou and J. Zhang, Prospects for axion searches with Advanced LIGO through binary mergers, Phys. Rev. D 99 (2019) 063013, [1807.02133].
  * [25] T. Dietrich and K. Clough, Cooling binary neutron star remnants via nucleon-nucleon-axion bremsstrahlung, Phys. Rev. D 100 (2019) 083005, [1909.01278].
  * [26] S. P. Harris, J.-F. Fortin, K. Sinha and M. G. Alford, Axions in neutron star mergers, JCAP 07 (2020) 023, [2003.09768].
  * [27] D. F. G. Fiorillo and F. Iocco, Axion-like particle production in neutron star mergers: gamma-ray signal detectability, 2109.10364.
  * [28] D. Croon, A. E. Nelson, C. Sun, D. G. E. Walker and Z.-Z. Xianyu, Hidden-Sector Spectroscopy with Gravitational Waves from Binary Neutron Stars, Astrophys. J. Lett. 858 (2018) L2, [1711.02096].
  * [29] J. A. Dror, R. Laha and T. Opferkuch, Probing muonic forces with neutron star binaries, Phys. Rev. D 102 (2020) 023005, [1909.12845].
  * [30] M. Diamond and G. Marques-Tavares, Gamma-rays flashes from dark photons in neutron star mergers, 2106.03879.
  * [31] S. Bernuzzi, Neutron Star Merger Remnants, Gen. Rel. Grav. 52 (2020) 108, [2004.06419].
  * [32] LIGO Scientific, Virgo collaboration, B. P. Abbott et al., A First Targeted Search for Gravitational-Wave Bursts from Core-Collapse Supernovae in Data of First-Generation Laser Interferometer Detectors, Phys. Rev. D 94 (2016) 102001, [1605.01785].
  * [33] E. Abdikamalov, G. Pagliaroli and D. Radice, Gravitational Waves from Core-Collapse Supernovae, 2010.04356.
  * [34] LIGO Scientific, Virgo collaboration, B. P. Abbott et al., GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\odot}$, Astrophys. J. Lett. 892 (2020) L3, [2001.01761].
  * [35] L. Baiotti, Gravitational waves from neutron star mergers and their relation to the nuclear equation of state, Prog. Part. Nucl. Phys. 109 (2019) 103714, [1907.08534].
  * [36] B. D. Metzger, Kilonovae, Living Rev. Rel. 23 (2020) 1, [1910.01617].
  * [37] K. Kyutoku and K. Kashiyama, Detectability of thermal neutrinos from binary-neutron-star mergers and implication to neutrino physics, Phys. Rev. D 97 (2018) 103001, [1710.05922].
  * [38] T. S. H. Schilbach, O. L. Caballero and G. C. McLaughlin, Black Hole Accretion Disk Diffuse Neutrino Background, Phys. Rev. D 100 (2019) 043008, [1808.03627].
  * [39] Z. Lin and C. Lunardini, Observing cosmological binary mergers with next generation neutrino and gravitational wave detectors, Phys. Rev. D 101 (2020) 023016, [1907.00034].
  * [40] M. S. Turner, Axions from SN 1987a, Phys. Rev. Lett. 60 (1988) 1797.
  * [41] G. Raffelt and D. Seckel, Bounds on Exotic Particle Interactions from SN 1987a, Phys. Rev. Lett. 60 (1988) 1793.
  * [42] 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–196.
  * [43] A. Burrows, M. S. Turner and R. P. Brinkmann, Axions and SN 1987a, Phys. Rev. D39 (1989) 1020.
  * [44] R. Mayle, J. R. Wilson, J. R. Ellis, K. A. Olive, D. N. Schramm and G. Steigman, Updated Constraints on Axions from SN 1987a, Phys. Lett. B 219 (1989) 515.
  * [45] A. Burrows, M. T. Ressell and M. S. Turner, Axions and SN1987A: Axion trapping, Phys. Rev. D42 (1990) 3297–3309.
  * [46] M. Giannotti and F. Nesti, Nucleon-nucleon Bremsstrahlung emission of massive axions, Phys. Rev. D72 (2005) 063005, [hep-ph/0505090].
  * [47] 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].
  * [48] M. K. Jaiswal, D. Singh, V. Singh and A. K. Ganguly, Differentiating Dilatons from Axions by their mixing with photons, 2107.11594.
  * [49] G. G. Raffelt, Astrophysical axion bounds, Lect. Notes Phys. 741 (2008) 51–71, [hep-ph/0611350].
  * [50] J.-F. Fortin and K. Sinha, Constraining Axion-Like-Particles with Hard X-ray Emission from Magnetars, JHEP 06 (2018) 048, [1804.01992].
  * [51] J.-F. Fortin, H.-K. Guo, S. P. Harris, E. Sheridan and K. Sinha, Magnetars and axion-like particles: probes with the hard X-ray spectrum, JCAP 06 (2021) 036, [2101.05302].
  * [52] I. G. Irastorza and J. Redondo, New experimental approaches in the search for axion-like particles, Prog. Part. Nucl. Phys. 102 (2018) 89–159, [1801.08127].
  * [53] K. Choi, S. H. Im and C. S. Shin, Recent Progress in the Physics of Axions and Axion-Like Particles, Ann. Rev. Nucl. Part. Sci. 71 (2021) 225–252, [2012.05029].
  * [54] ATLAS collaboration, G. Aad et al., Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC, Phys. Lett. B 716 (2012) 1–29, [1207.7214].
  * [55] CMS collaboration, S. Chatrchyan et al., Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC, Phys. Lett. B 716 (2012) 30–61, [1207.7235].
  * [56] V. Silveira and A. Zee, Scalar Phantoms, Phys. Lett. 161B (1985) 136–140.
  * [57] H. Davoudiasl, R. Kitano, T. Li and H. Murayama, The New minimal standard model, Phys. Lett. B 609 (2005) 117–123, [hep-ph/0405097].
  * [58] R. M. Schabinger and J. D. Wells, A Minimal spontaneously broken hidden sector and its impact on Higgs boson physics at the large hadron collider, Phys. Rev. D 72 (2005) 093007, [hep-ph/0509209].
  * [59] 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].
  * [60] V. Barger, P. Langacker, M. McCaskey, M. J. Ramsey-Musolf and G. Shaughnessy, LHC Phenomenology of an Extended Standard Model with a Real Scalar Singlet, Phys. Rev. D 77 (2008) 035005, [0706.4311].
  * [61] T. Robens and T. Stefaniak, Status of the Higgs Singlet Extension of the Standard Model after LHC Run 1, Eur. Phys. J. C 75 (2015) 104, [1501.02234].
  * [62] M. Gonderinger, Y. Li, H. Patel and M. J. Ramsey-Musolf, Vacuum Stability, Perturbativity, and Scalar Singlet Dark Matter, JHEP 01 (2010) 053, [0910.3167].
  * [63] M. Gonderinger, H. Lim and M. J. Ramsey-Musolf, Complex Scalar Singlet Dark Matter: Vacuum Stability and Phenomenology, Phys. Rev. D 86 (2012) 043511, [1202.1316].
  * [64] O. Lebedev, On Stability of the Electroweak Vacuum and the Higgs Portal, Eur. Phys. J. C 72 (2012) 2058, [1203.0156].
  * [65] J. Elias-Miro, J. R. Espinosa, G. F. Giudice, H. M. Lee and A. Strumia, Stabilization of the Electroweak Vacuum by a Scalar Threshold Effect, JHEP 06 (2012) 031, [1203.0237].
  * [66] N. Khan and S. Rakshit, Study of electroweak vacuum metastability with a singlet scalar dark matter, Phys. Rev. D 90 (2014) 113008, [1407.6015].
  * [67] A. Falkowski, C. Gross and O. Lebedev, A second Higgs from the Higgs portal, JHEP 05 (2015) 057, [1502.01361].
  * [68] P. Ghorbani, Vacuum stability vs. positivity in real singlet scalar extension of the standard model, Nucl. Phys. B 971 (2021) 115533, [2104.09542].
  * [69] J. R. Espinosa and M. Quiros, The Electroweak phase transition with a singlet, Phys. Lett. B305 (1993) 98–105, [hep-ph/9301285].
  * [70] J. Choi and R. R. Volkas, Real Higgs singlet and the electroweak phase transition in the Standard Model, Phys. Lett. B 317 (1993) 385–391, [hep-ph/9308234].
  * [71] S. W. Ham, Y. S. Jeong and S. K. Oh, Electroweak phase transition in an extension of the standard model with a real Higgs singlet, J. Phys. G 31 (2005) 857–871, [hep-ph/0411352].
  * [72] S. Profumo, M. J. Ramsey-Musolf and G. Shaughnessy, Singlet Higgs phenomenology and the electroweak phase transition, JHEP 08 (2007) 010, [0705.2425].
  * [73] J. R. Espinosa, T. Konstandin and F. Riva, Strong Electroweak Phase Transitions in the Standard Model with a Singlet, Nucl. Phys. B 854 (2012) 592–630, [1107.5441].
  * [74] V. Barger, D. J. H. Chung, A. J. Long and L.-T. Wang, Strongly First Order Phase Transitions Near an Enhanced Discrete Symmetry Point, Phys. Lett. B 710 (2012) 1–7, [1112.5460].
  * [75] S. Profumo, M. J. Ramsey-Musolf, C. L. Wainwright and P. Winslow, Singlet-catalyzed electroweak phase transitions and precision Higgs boson studies, Phys. Rev. D 91 (2015) 035018, [1407.5342].
  * [76] D. Curtin, P. Meade and C.-T. Yu, Testing Electroweak Baryogenesis with Future Colliders, JHEP 11 (2014) 127, [1409.0005].
  * [77] A. V. Kotwal, M. J. Ramsey-Musolf, J. M. No and P. Winslow, Singlet-catalyzed electroweak phase transitions in the 100 TeV frontier, Phys. Rev. D 94 (2016) 035022, [1605.06123].
  * [78] C.-Y. Chen, J. Kozaczuk and I. M. Lewis, Non-resonant Collider Signatures of a Singlet-Driven Electroweak Phase Transition, JHEP 08 (2017) 096, [1704.05844].
  * [79] P. W. Graham, D. E. Kaplan and S. Rajendran, Cosmological Relaxation of the Electroweak Scale, Phys. Rev. Lett. 115 (2015) 221801, [1504.07551].
  * [80] T. Flacke, C. Frugiuele, E. Fuchs, R. S. Gupta and G. Perez, Phenomenology of relaxion-Higgs mixing, JHEP 06 (2017) 050, [1610.02025].
  * [81] C. Frugiuele, E. Fuchs, G. Perez and M. Schlaffer, Relaxion and light (pseudo)scalars at the HL-LHC and lepton colliders, JHEP 10 (2018) 151, [1807.10842].
  * [82] A. Banerjee, H. Kim, O. Matsedonskyi, G. Perez and M. S. Safronova, Probing the Relaxed Relaxion at the Luminosity and Precision Frontiers, JHEP 07 (2020) 153, [2004.02899].
  * [83] P. Brax and C. Burrage, Screening the Higgs portal, Phys. Rev. D 104 (2021) 015011, [2101.10693].
  * [84] R. Foot and A. Kobakhidze, Electroweak Scale Invariant Models with Small Cosmological Constant, Int. J. Mod. Phys. A 30 (2015) 1550126, [1112.0607].
  * [85] M. Heikinheimo, A. Racioppi, M. Raidal, C. Spethmann and K. Tuominen, Physical Naturalness and Dynamical Breaking of Classical Scale Invariance, Mod. Phys. Lett. A 29 (2014) 1450077, [1304.7006].
  * [86] Z.-W. Wang, T. G. Steele, T. Hanif and R. B. Mann, Conformal Complex Singlet Extension of the Standard Model: Scenario for Dark Matter and a Second Higgs Boson, JHEP 08 (2016) 065, [1510.04321].
  * [87] M. Pospelov, A. Ritz and M. B. Voloshin, Secluded WIMP Dark Matter, Phys. Lett. B662 (2008) 53–61, [0711.4866].
  * [88] S. Baek, P. Ko and W.-I. Park, Search for the Higgs portal to a singlet fermionic dark matter at the LHC, JHEP 02 (2012) 047, [1112.1847].
  * [89] S. Baek, P. Ko, W.-I. Park and E. Senaha, Vacuum structure and stability of a singlet fermion dark matter model with a singlet scalar messenger, JHEP 11 (2012) 116, [1209.4163].
  * [90] S. Baek, P. Ko, W.-I. Park and E. Senaha, Higgs Portal Vector Dark Matter : Revisited, JHEP 05 (2013) 036, [1212.2131].
  * [91] K. Schmidt-Hoberg, F. Staub and M. W. Winkler, Constraints on light mediators: confronting dark matter searches with B physics, Phys. Lett. B 727 (2013) 506–510, [1310.6752].
  * [92] G. Krnjaic, Probing Light Thermal Dark-Matter With a Higgs Portal Mediator, Phys. Rev. D94 (2016) 073009, [1512.04119].
  * [93] A. Beniwal, F. Rajec, C. Savage, P. Scott, C. Weniger, M. White et al., Combined analysis of effective Higgs portal dark matter models, Phys. Rev. D 93 (2016) 115016, [1512.06458].
  * [94] R. Contino, M. Ghezzi, C. Grojean, M. Muhlleitner and M. Spira, Effective Lagrangian for a light Higgs-like scalar, JHEP 07 (2013) 035, [1303.3876].
  * [95] C.-W. Chiang and R. Huo, Standard Model Effective Field Theory: Integrating out a Generic Scalar, JHEP 09 (2015) 152, [1505.06334].
  * [96] S. Dawson and C. W. Murphy, Standard Model EFT and Extended Scalar Sectors, Phys. Rev. D 96 (2017) 015041, [1704.07851].
  * [97] M. Jiang, N. Craig, Y.-Y. Li and D. Sutherland, Complete one-loop matching for a singlet scalar in the Standard Model EFT, JHEP 02 (2019) 031, [1811.08878]. [Erratum: JHEP 01, 135 (2021)].
  * [98] LHC Higgs Cross Section Working Group collaboration, S. Dittmaier et al., Handbook of LHC Higgs Cross Sections: 1. Inclusive Observables, 1101.0593.
  * [99] S. Dawson et al., Working Group Report: Higgs Boson, in Community Summer Study 2013: Snowmass on the Mississippi, 10, 2013, 1310.8361.
  * [100] G. Brooijmans et al., Les Houches 2017: Physics at TeV Colliders New Physics Working Group Report, in 10th Les Houches Workshop on Physics at TeV Colliders, 3, 2018, 1803.10379.
  * [101] M. W. Winkler, Decay and detection of a light scalar boson mixing with the Higgs boson, Phys. Rev. D 99 (2019) 015018, [1809.01876].
  * [102] D. Egana-Ugrinovic, S. Homiller and P. Meade, Light Scalars and the Koto Anomaly, Phys. Rev. Lett. 124 (2020) 191801, [1911.10203].
  * [103] P. S. B. Dev, R. N. Mohapatra and Y. Zhang, Constraints on long-lived light scalars with flavor-changing couplings and the KOTO anomaly, Phys. Rev. D 101 (2020) 075014, [1911.12334].
  * [104] S. Foroughi-Abari and A. Ritz, LSND Constraints on the Higgs Portal, Phys. Rev. D 102 (2020) 035015, [2004.14515].
  * [105] G. Lanfranchi, M. Pospelov and P. Schuster, The Search for Feebly-Interacting Particles, 2011.02157.
  * [106] P. S. B. Dev, R. N. Mohapatra and Y. Zhang, Revisiting supernova constraints on a light CP-even scalar, JCAP 2008 (2020) 003, [2005.00490]. [Erratum: JCAP2011,E01(2020)].
  * [107] P. S. B. Dev, R. N. Mohapatra and Y. Zhang, Stellar limits on light CP-even scalar, JCAP 05 (2021) 014, [2010.01124].
  * [108] A. Fradette and M. Pospelov, BBN for the LHC: constraints on lifetimes of the Higgs portal scalars, Phys. Rev. D 96 (2017) 075033, [1706.01920].
  * [109] P. S. B. Dev, R. N. Mohapatra and Y. Zhang, Lepton Flavor Violation Induced by a Neutral Scalar at Future Lepton Colliders, Phys. Rev. Lett. 120 (2018) 221804, [1711.08430].
  * [110] B. Batell, A. Freitas, A. Ismail and D. Mckeen, Flavor-specific scalar mediators, Phys. Rev. D 98 (2018) 055026, [1712.10022].
  * [111] A. de Gouvêa, P. S. B. Dev, B. Dutta, T. Ghosh, T. Han and Y. Zhang, Leptonic Scalars at the LHC, JHEP 07 (2020) 142, [1910.01132].
  * [112] P. S. B. Dev, B. Dutta, T. Ghosh, T. Han, H. Qin and Y. Zhang, Leptonic Scalars and Collider Signatures in a UV-complete Model, 2109.04490.
  * [113] J. C. Pati and A. Salam, Lepton Number as the Fourth Color, Phys. Rev. D10 (1974) 275–289. [Erratum: Phys. Rev.D11,703(1975)].
  * [114] R. N. Mohapatra and J. C. Pati, A Natural Left-Right Symmetry, Phys. Rev. D11 (1975) 2558.
  * [115] G. Senjanović and R. N. Mohapatra, Exact Left-Right Symmetry and Spontaneous Violation of Parity, Phys. Rev. D12 (1975) 1502.
  * [116] M. Nemevšek, G. Senjanović and Y. Zhang, Warm Dark Matter in Low Scale Left-Right Theory, JCAP 1207 (2012) 006, [1205.0844].
  * [117] A. Maiezza, G. Senjanović and J. C. Vasquez, Higgs sector of the minimal left-right symmetric theory, Phys. Rev. D 95 (2017) 095004, [1612.09146].
  * [118] P. S. B. Dev, R. N. Mohapatra and Y. Zhang, Displaced photon signal from a possible light scalar in minimal left-right seesaw model, Phys. Rev. D 95 (2017) 115001, [1612.09587].
  * [119] P. S. B. Dev, R. N. Mohapatra and Y. Zhang, Long Lived Light Scalars as Probe of Low Scale Seesaw Models, Nucl. Phys. B923 (2017) 179–221, [1703.02471].
  * [120] M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Remarks on Higgs Boson Interactions with Nucleons, Phys. Lett. 78B (1978) 443–446.
  * [121] H.-Y. Cheng, Low-energy Interactions of Scalar and Pseudoscalar Higgs Bosons With Baryons, Phys. Lett. B219 (1989) 347–353.
  * [122] M. B. Voloshin, Once Again About the Role of Gluonic Mechanism in Interaction of Light Higgs Boson with Hadrons, Sov. J. Nucl. Phys. 44 (1986) 478. [Yad. Fiz.44,738(1986)].
  * [123] J. F. Donoghue, J. Gasser and H. Leutwyler, The Decay of a Light Higgs Boson, Nucl. Phys. B343 (1990) 341–368.
  * [124] NA48/2 collaboration, J. R. Batley et al., Precise measurement of the $K^{\pm}\to\pi^{\pm}e^{+}e^{-}$ decay, Phys. Lett. B677 (2009) 246–254, [0903.3130].
  * [125] NA48/2 collaboration, J. R. Batley et al., New measurement of the $K^{\pm}\to\pi^{\pm}\mu^{+}\mu^{-}$ decay, Phys. Lett. B697 (2011) 107–115, [1011.4817].
  * [126] NA62 collaboration, C. Lazzeroni et al., Study of the $K^{\pm}\to\pi^{\pm}\gamma\gamma$ decay by the NA62 experiment, Phys. Lett. B732 (2014) 65–74, [1402.4334].
  * [127] NA62 collaboration, E. Cortina Gil et al., An investigation of the very rare ${K}^{+}\to{\pi}^{+}\nu\overline{\nu}$ decay, JHEP 11 (2020) 042, [2007.08218].
  * [128] BNL-E949 collaboration, A. V. Artamonov et al., Study of the decay $K^{+}\to\pi^{+}\nu\bar{\nu}$ in the momentum region $140<P_{\pi}<199$ MeV/c, Phys. Rev. D79 (2009) 092004, [0903.0030].
  * [129] KTeV collaboration, A. Alavi-Harati et al., Search for the rare decay $K_{L}\to\pi^{0}e^{+}e^{-}$, Phys. Rev. Lett. 93 (2004) 021805, [hep-ex/0309072].
  * [130] KTeV collaboration, A. Alavi-Harati et al., Search for the Decay $K_{L}\to\pi^{0}\mu^{+}\mu^{-}$, Phys. Rev. Lett. 84 (2000) 5279–5282, [hep-ex/0001006].
  * [131] KTeV collaboration, T. Alexopoulos et al., Measurements of $K_{L}$ branching fractions and the CP violation parameter $|\eta^{\pm}|$, Phys. Rev. D70 (2004) 092006, [hep-ex/0406002].
  * [132] KTeV collaboration, E. Abouzaid et al., Final Results from the KTeV Experiment on the Decay $K_{L}\to\pi^{0}\gamma\gamma$, Phys. Rev. D77 (2008) 112004, [0805.0031].
  * [133] KOTO collaboration, J. K. Ahn et al., Search for the $K_{L}\\!\to\\!\pi^{0}\nu\overline{\nu}$ and $K_{L}\\!\to\\!\pi^{0}X^{0}$ decays at the J-PARC KOTO experiment, Phys. Rev. Lett. 122 (2019) 021802, [1810.09655].
  * [134] LHCb collaboration, R. Aaij et al., Differential branching fraction and angular analysis of the $B^{+}\rightarrow K^{+}\mu^{+}\mu^{-}$ decay, JHEP 02 (2013) 105, [1209.4284].
  * [135] BaBar collaboration, B. Aubert et al., Evidence for the rare decay $B\to K^{*}\ell^{+}\ell^{-}$ and measurement of the $B\to K\ell^{+}\ell^{-}$ branching fraction, Phys. Rev. Lett. 91 (2003) 221802, [hep-ex/0308042].
  * [136] BaBar collaboration, J. P. Lees et al., Search for $B\to K^{(*)}\nu\overline{\nu}$ and invisible quarkonium decays, Phys. Rev. D87 (2013) 112005, [1303.7465].
  * [137] Belle collaboration, J. T. Wei et al., Measurement of the Differential Branching Fraction and Forward-Backword Asymmetry for $B\to K^{(*)}\ell^{+}\ell^{-}$, Phys. Rev. Lett. 103 (2009) 171801, [0904.0770].
  * [138] MicroBooNE collaboration, P. Abratenko et al., Search for a Higgs portal scalar decaying to electron-positron pairs in the MicroBooNE detector, Phys. Rev. Lett. 127 (2021) 151803, [2106.00568].
  * [139] CHARM collaboration, F. Bergsma et al., Search for Axion Like Particle Production in 400-GeV Proton - Copper Interactions, Phys. Lett. B 157 (1985) 458–462.
  * [140] LSND collaboration, A. Aguilar-Arevalo et al., Evidence for neutrino oscillations from the observation of $\bar{\nu}_{e}$ appearance in a $\bar{\nu}_{\mu}$ beam, Phys. Rev. D 64 (2001) 112007, [hep-ex/0104049].
  * [141] LSND collaboration, L. B. Auerbach et al., Search for pi0 —$>$ nu(mu) anti-nu(mu) decay in LSND, Phys. Rev. Lett. 92 (2004) 091801, [hep-ex/0310060].
  * [142] N. Ishizuka and M. Yoshimura, Axion and Dilaton Emissivity From Nascent Neutron Stars, Prog. Theor. Phys. 84 (1990) 233–250.
  * [143] R. Diener and C. P. Burgess, Bulk Stabilization, the Extra-Dimensional Higgs Portal and Missing Energy in Higgs Events, JHEP 05 (2013) 078, [1302.6486].
  * [144] J. S. Lee, Revisiting Supernova 1987A Limits on Axion-Like-Particles, 1808.10136.
  * [145] D. Arndt and P. J. Fox, Saxion emission from SN1987A, JHEP 02 (2003) 036, [hep-ph/0207098].
  * [146] K. Kainulainen, K. Tuominen and V. Vaskonen, Self-interacting dark matter and cosmology of a light scalar mediator, Phys. Rev. D 93 (2016) 015016, [1507.04931]. [Erratum: Phys.Rev.D 95, 079901 (2017)].
  * [147] Planck collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, [1807.06209]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [148] ATLAS collaboration, Combination of searches for invisible Higgs boson decays with the ATLAS experiment, . ATLAS-CONF-2020-052.
  * [149] B. Batell, J. Berger and A. Ismail, Probing the Higgs Portal at the Fermilab Short-Baseline Neutrino Experiments, Phys. Rev. D 100 (2019) 115039, [1909.11670].
  * [150] J. M. Berryman, A. de Gouvea, P. J. Fox, B. J. Kayser, K. J. Kelly and J. L. Raaf, Searches for Decays of New Particles in the DUNE Multi-Purpose Near Detector, JHEP 02 (2020) 174, [1912.07622].
  * [151] J. L. Feng, I. Galon, F. Kling and S. Trojanowski, Dark Higgs bosons at the ForwArd Search ExpeRiment, Phys. Rev. D 97 (2018) 055034, [1710.09387].
  * [152] FASER collaboration, A. Ariga et al., FASER’s physics reach for long-lived particles, Phys. Rev. D 99 (2019) 095011, [1811.12522].
  * [153] NA62 collaboration, B. Döbrich, Dark Sectors at fixed targets: The example of NA62, Frascati Phys. Ser. 66 (2018) 312–327, [1807.10170].
  * [154] V. V. Gligorov, S. Knapen, M. Papucci and D. J. Robinson, Searching for Long-lived Particles: A Compact Detector for Exotics at LHCb, Phys. Rev. D 97 (2018) 015023, [1708.09395].
  * [155] 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].
  * [156] D. Curtin et al., Long-Lived Particles at the Energy Frontier: The MATHUSLA Physics Case, Rept. Prog. Phys. 82 (2019) 116201, [1806.07396].
  * [157] J. P. Chou, D. Curtin and H. J. Lubatti, New Detectors to Explore the Lifetime Frontier, Phys. Lett. B 767 (2017) 29–36, [1606.06298].
  * [158] L. A. Anchordoqui et al., The Forward Physics Facility: Sites, Experiments, and Physics Potential, 2109.10905.
  * [159] L. McLerran and S. Reddy, Quarkyonic Matter and Neutron Stars, Phys. Rev. Lett. 122 (2019) 122701, [1811.12503].
  * [160] E. Annala, T. Gorda, A. Kurkela, J. Nättilä and A. Vuorinen, Evidence for quark-matter cores in massive neutron stars, Nature Phys. 16 (2020) 907–910, [1903.09121].
  * [161] M. G. Alford, S. Han and K. Schwenzer, Signatures for quark matter from multi-messenger observations, J. Phys. G 46 (2019) 114001, [1904.05471].
  * [162] J. Schaffner-Bielich, Compact Star Physics.  Cambridge University Press, 8, 2020, 10.1017/9781316848357.
  * [163] L. Rezzolla and O. Zanotti, Relativistic hydrodynamics.  Oxford University Press, 2013.
  * [164] L. Baiotti and L. Rezzolla, Binary neutron star mergers: a review of Einstein’s richest laboratory, Rept. Prog. Phys. 80 (2017) 096901, [1607.03540].
  * [165] M. D. Duez and Y. Zlochower, Numerical Relativity of Compact Binaries in the 21st Century, Rept. Prog. Phys. 82 (2019) 016902, [1808.06011].