# How Macroscopic Limits on Neutron Star Baryon Loss Yield Microscopic Limits
on Non-Standard-Model Baryon Decay

Jeffrey M. Berryman<EMAIL_ADDRESS>Department of Physics,
University of California, Berkeley, CA 94720, USA  Institute for Nuclear
Theory, University of Washington, Seattle, WA 98195, USA  Center for Neutrino
Physics, Physics Department, Virginia Tech, Blacksburg, VA 24061, USA  Susan
Gardner<EMAIL_ADDRESS>Department of Physics and Astronomy, University
of Kentucky, Lexington, KY 40506-0055, USA  Mohammadreza Zakeri
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of Kentucky,
Lexington, KY 40506-0055, USA

###### Abstract

We investigate how our baryon-loss limits from anomalous binary-pulsar period
lengthening can be interpreted microscopically to yield specific constraints
on the particle physics of baryon number violation within a neutron star. We
focus on the possibility of anomalous baryon disappearance via dark baryon
processes and on scenarios in which the produced dark-sector particles do not
survive to influence the response of the star to baryon-number-violating
effects. We flesh out the conditions for which this may occur, as well as
other key assumptions. We then turn to the analysis of particle processes in
the dense nuclear medium found at the core of a neutron star, employing the
techniques of relativistic mean-field theory. Using our study of in-medium
effects and limits on macroscopic baryon number violation we extract limits on
in-vacuum baryon-number-violating processes, and we determine them for various
equations of state. We conclude by noting the implications of our results for
models of dark-sector-enabled baryogenesis.

††preprint: FERMILAB-PUB-23-224-T, INT-PUB-23-001, N3AS-23-012

###### Contents

  1. 1 Introduction
  2. 2 Particle Physics Models of Baryon Dark Decays
    1. 2.1 Dark Baryon Removal Rates
  3. 3 Macroscopic Baryon Number Violation in Neutron Stars
    1. 3.1 Assumptions
    2. 3.2 Framework
    3. 3.3 Observables
    4. 3.4 Interpretation
  4. 4 Dense Matter Considerations for Particle Processes
    1. 4.1 Modelling Dense Matter
    2. 4.2 Medium Effects: Effective Masses and Beyond
  5. 5 Baryon Dark Decay Rates in Dense Matter
    1. 5.1 General Considerations
    2. 5.2 Dark Decay Rate Estimates
    3. 5.3 Medium-Enabled Dark Decay Processes
    4. 5.4 Total Rates
  6. 6 Inferred Limits on Baryon Dark Decays
  7. 7 Implications for Models of Baryogenesis and Dark Matter
  8. 8 Summary
  9. A In-Medium Electromagnetic Form Factors
  10. B Nonrelativistic Limit of In-Medium Scattering
  11. C In-Medium Compton Scattering
  12. D Fermion Mixing in Dense Matter
  13. E Baryon Decays to $\chi+\gamma$
    1. E.1 Matrix Element
    2. E.2 Integrated Rates

##  1 Introduction

The cosmic excess of baryons over antibaryons is well established [1], but the
theoretical mechanism by which it is produced is not. The essential
theoretical ingredients are thought to be known: baryon number violation
(BNV), along with C and CP violation, must all be present in a non-equilibrium
environment [2]. Thus BNV would seem to play an essential role, though in the
Standard Model (SM) BNV is thought to occur appreciably only at extremely high
temperature [3, 4] — and the existence of BNV at low energies has as yet to be
established. In this paper we continue our scrutiny of such effects through
observations of neutron stars, which contain enormous reservoirs of baryons.
In earlier work we identified sensitive limits on BNV through the
interpretation of precise observations of energy loss in isolated neutron
stars and in neutron-star binary systems [5]. These studies limit the baryon-
number-violating effects that occur across the entirety of a neutron star. In
this sense they are macroscopic limits. In this paper we interpret these
limits in a microscopic way, in that we develop a framework in which they can
be translated to limits on the parameters of particular particle physics
models that generate baryon-number-violating effects.

The particular models to which our studies are most sensitive are those in
which baryons decay or otherwise transform to dark-sector fermions, of ${\cal
O}(1\,\rm GeV)$ in mass, that carry baryon number. In such cases BNV becomes
an apparent, rather than explicit, effect, because the dark-sector particles
are unobserved, even if baryon number is not broken. Although the existence of
dark matter is certainly established through astrometric observations, both
its nature and origin continue to be open questions. It is possible that the
origins of dark matter and of the cosmic baryon asymmetry are related, so that
the loosely similar value of the cosmic baryon and dark matter energy
densities today may follow from a single underlying model [6]. The possibility
of baryons that connect to hidden-sector baryons of comparable mass figure in
many such explanations. In this paper we constrain this possibility through
the study of neutron and hyperon transitions to final states with dark baryons
in the neutron star. To our knowledge, an in-depth, quantitative study of non-
SM processes within dense nuclear matter has not previously been realized
111Albeit studies of exotic light particle emission in dense matter, which
possesses simplifying aspects, are of long standing [7] and continue to be
investigated [8]., and its execution necessitates much care. The existence of
neutron stars of about $2\,M_{\astrosun}$ in mass speaks to central densities
in excess of three times nuclear matter saturation density, so that in this
paper we employ relativistic mean field theory in baryonic degrees of freedom
for our dense matter description, as its accuracy should improve with
increasing density — and thus it should work best at the core of the star. We
note that a neutron star may become a hybrid star, i.e., one with a quark-
based core predicated by a finite-density quark-hadron phase transition, if it
is sufficiently heavy, and this possibility can also be constructed within
this framework [9]. Transitions to dark baryons could also occur within the
quark-based core, though we will set aside this possibility in this paper —
and revisit it only in offering an assessment of our uncertainties in our
concluding summary.

The broader possibility of dark decays of the neutron has been noted in
explanation [10, 11] of the long-standing neutron lifetime anomaly [12], in
which the lifetime inferred from counting surviving neutrons is significantly
different from that inferred from counting the protons subsequent to ordinary
neutron decay. Although the discrepancy may arise from experimental effects,
the possibility that dark decays contribute to it in some measure is a
continuing possibility [5]. In this paper we provide severe limits on the
flavor structure of possible new-physics models with dark baryonic sectors,
such as Refs. [13, 10, 14, 15, 16], that arise from the interpretation of
neutron-star energy loss constraints we developed in Ref. [5]. In this paper
we flesh out the general assumptions of that earlier analysis and note how the
specific models we consider can satisfy them.

Let us conclude our introduction with a brief outline of the body of our
paper. In Sec. 2 we detail the models of baryon dark decays that we are able
to constrain through our neutron star studies, and we note how they are
distinct from models that we cannot. We also compute baryon dark decay rates
in vacuum, for later reference, as well as dark baryon removal rates, because
our analysis assumes that SM dynamics determine the response of the star in
the presence of BNV. In Sec. 3 we consider macroscopic baryon number violation
in neutron stars, revisiting our earlier work [5] and fleshing out constraints
following from its assumptions in greater detail. In Sec. 4 we develop how to
evaluate particle processes within dense matter, employing RMF, as developed
in Refs. [17, 18, 19], to describe the neutron star medium in
$\beta$-equilibrium [20, 21]. In this context, uncertainties in our
description of the dense medium are captured through variations in the
equation of state (EoS). With these developments in hand, we evaluate particle
processes within our framework for the dense nuclear medium of a neutron star
in Sec. 5 and use our macroscopic limits on BNV from Sec. 3 to report limits
on the parameters of the microscopic models we consider in Sec. 6. In Sec. 7
we consider the implications of our results for models of dark-sector
baryogenesis and dark matter, and we offer a summary and outlook in Sec. 8.

##  2 Particle Physics Models of Baryon Dark Decays

The possibility of hadronic processes with dark-sector particles naturally
emerges in models that explain both the origin of dark matter and the cosmic
baryon asymmetry, particularly if the dark sector candidate carries a baryonic
charge [22, 13, 14]. Although it has long been thought that dark matter could
also be described as a relic asymmetry [23, 24], in these models, rather, the
two problems are solved simultaneously [6]. More recently, highly testable
scenarios [25] have been developed [26, 27, 28, 29, 16, 30], and we probe
their flavor structure through the studies of this paper — and in Sec. 7 we
consider the implications of the constraints that we find. Since the dark-
sector particles are presumably SM gauge singlets, they could be light in
mass, potentially with masses comparable to that of the known hadrons, and yet
have escaped experimental detection thus far.

Our current discussion is loosely inspired by models connected to explanations
of the neutron lifetime anomaly [10, 31, 16], with neutrons decaying to a dark
baryon with a photon or an $e^{+}e^{-}$ pair. Models with similar content have
been considered for broader purposes [32, 33, 34, 15], and alternative
solutions have also been noted [35, 36, 37]. The dark channels in the various
models would impact the determined bottle lifetime, with a mirror neutron
model [35] serving as a rare exception. There, neutron-to-mirror-neutron
conversion occurs in a strong magnetic field, impacting the ability to detect
protons in the beam-lifetime experiment. This last possibility has been
excluded as a complete explanation of the anomaly by a direct experimental
search [38]. We note that models that would explain the anomaly through
neutron disappearance or decay to dark-sector final states can also be
constrained by the close empirical agreement of the neutron lifetime with its
measured $A$ decay correlation as interpreted in the SM [39, 40, 5]. This
agreement limits the branching ratio on such exotic processes to [41]:

${\rm Br}(n\to\rm exotics)<0.16\%\,\,(95\%\,\text{one-sided C.L.})\,,$ (2.1)

where we note that the neutron lifetime anomaly is roughly a 1% effect
[10]222The most precise measurement of the $A$ correlation coefficient yields
the ratio of the axial-vector to vector coupling constants
$|\lambda|=1.27641(56)$ [42], but recent measurements of the $a$ correlation
do not completely fit this picture, yielding $|\lambda|=1.2677(28)$ [43] and
$|\lambda|=1.2796(62)$ [44].. Direct experimental limits on $n\to\chi\gamma$
[45] and $n\to\chi e^{+}e^{-}$ [46] decays also exist, removing ranges of
parameter space as an explanation of the anomaly. We will be able to set much
more severe limits through our studies, where we note the limit on
$\Lambda\to\chi\gamma$ from SN1987 for reference [16]. We regard the neutron
lifetime anomaly as a motivation for further investigation of baryon dark
decays, with new limits constraining the manner in which the co-genesis of
dark matter and the cosmic baryon asymmetry could possibly occur. We now turn
to the development of models of dark baryon decays.

Following Ref. [33], we introduce a Dirac fermion $\chi$ with baryon number
$B=+1$ which interacts with SM quarks via the generic form

${\cal
L_{\chi}}=\bar{\chi}(i/{\\!\\!\\!\partial}-m_{\chi})\chi+\left(\frac{u_{i}d_{j}d_{k}\chi_{L}^{c}}{\Lambda_{ijk}^{2}}+\frac{Q_{i}Q_{j}d_{k}\chi_{L}^{c}}{\tilde{\Lambda}_{ijk}^{2}}+{\rm
h.c.}\right)\,,$ (2.2)

where $i,j,k$ are generational indices, $Q$ and $q$ denote a left-handed quark
doublet and a right-handed quark, respectively — and color and Lorentz indices
are left implicit. Such interactions can generate both decay and scattering
processes involving dark final states, which we consider closely in this
paper. First, though, we address their flavor structure. We could neglect this
possibility altogether, dropping all subscript dependence, but simple,
renormalizable models that produce Eq. (2.2), at energies below the mass scale
of their new physics, show that strong flavor sensitivity can nevertheless
exist. Turning to models with leptoquarks [10, 15], we consider colored
scalars $S_{1}$ and ${\bar{S}}_{1}$ transforming as $(\bar{3},1,1/3)$ 333This
variant was first considered in Ref. [10]. and $(\bar{3},1,-2/3)$,
respectively, under the SM gauge groups and SM invariant scalar-fermion
interactions. Non-trivial flavor structure follows from the choice of
leptoquark in that $S_{1}$ can mediate both $n\to\chi\gamma$ and
$\Lambda\to\chi\gamma$ decay at tree level, whereas ${\bar{S}}_{1}$ can
mediate $\Lambda\to\chi\gamma$ at tree level but to mediate $n\to\chi\gamma$
requires a one-loop process with $W^{\pm}$ exchange as well [15]. Thus in this
paper we strive to probe both $n\to\chi\gamma$ and $\Lambda\to\chi\gamma$
decay processes. These models also readily generate proton decay [33, 10, 15],
noting $p\to\chi\pi^{+}$ or $p\to\chi K^{+}$ decay as examples, so that the
possible range of $\chi$ masses is rather restricted as a result. We note that
the stability of the 9Be nucleus [10], particularly stability against
$\rm{}^{9}Be\to\chi\alpha\alpha$ decay [47], requires

$m_{\chi}\geq 0.937993\,\rm GeV\,,$ (2.3)

slightly in excess of the proton stability constraint $m_{\chi}>m_{p}-m_{e}$,
and that atomic hydrogen is stable if $m_{\chi}>m_{p}+m_{e}=0.93878\,\rm GeV$
[34]. If either constraint were not satisfied, then the empirical limit on the
pertinent lifetime would bound the parameters of the model. Within the SM both
systems are absolutely stable, yet empirical tests of that, with a determined
lifetime as an outcome, should be possible. We note H lifetime estimates, made
finite through a model with a suitably light $\chi$, are made in Ref. [34].
Moreover, the radiative decay $\rm H\to\nu\chi\gamma$, which is subdominant
relative to $\rm H\to\nu\chi$, can be probed through measurements at Borexino
[34, 48]. Similar expectations follow from violating Eq. (2.3) — and a
concrete estimate of the 9Be lifetime can be found in Ref. [49].

In what follows we ignore the possible chiral structure of the quark-$\chi$
couplings and simply consider [16]

${\cal L}\supset\frac{u^{i}d^{j}d^{k}\chi}{\Lambda^{2}}+\rm h.c.\,$ (2.4)

Since the quarks carry electric charge, we have, at the energy scales for
which baryonic degrees of freedom are pertinent,

${\cal L}_{\rm
n}=\bar{n}\left(i/{\\!\\!\\!\partial}-m_{n}+\frac{g_{n}e}{8m_{n}}\sigma^{\alpha\beta}F_{\alpha\beta}\right)n+\bar{\chi}(i/{\\!\\!\\!\partial}-m_{\chi})\chi+\varepsilon_{n\chi}(\bar{n}\chi+\bar{\chi}n)\,,$
(2.5)

noting $g_{n}=-3.826$ is the $g$-factor of the neutron [50]. This form also
holds for the $\Lambda$ upon the replacement $n\to\Lambda$, taking
$g_{\Lambda}=-1.22$. After redefining the fields to remove the mixing term in
Eq. (2.5), then if $\varepsilon\ll m_{n}-m_{\chi}$, with $m_{\chi}<m_{n}$, we
have [10, 15]

${\cal
L}_{n\to\chi\gamma}=\frac{g_{n}e}{8m_{n}}\frac{\varepsilon_{n\chi}}{m_{n}-m_{\chi}}\bar{\chi}\sigma^{\alpha\beta}F_{\alpha\beta}n\,,$
(2.6)

though potentially this operator could also stem from a distinct higher-energy
source. Generally, the interaction of Eq. (2.4) can also generate transitions
to dark baryon states with mesons, such as the decays $n\to\chi+{\rm meson}$
or $\Lambda\to\chi+\rm meson$. Ref. [16] uses chiral effective theory [51] to
relate the possibilities. We eschew this path because chiral effective theory
ceases to be valid if the density of the neutron star medium much exceeds that
of nuclear matter saturation density. Since our particular purpose is to set
limits on microscopic models given BNV limits determined from observables
associated with the entire neutron star, we set aside the study of final
states containing both dark and hadronic degrees of freedom. They are distinct
from the final states we do study, and cancellations cannot occur. We thus
expect that including these additional decays with hadrons can only make our
limits more severe, though the inclusion of hadronic channels would make our
estimates less sure.

Figure 1: (Color Online) Illustration of $n\to\chi\gamma$ decay in the degrees
of freedom of Eq. (2.4). The decay $\Lambda\to\chi\gamma$ follows from the
replacement of one $d$ quark with a $s$ quark.

At low energies, the magnetic interaction of Eq. (2.6), employed in Refs. [10,
15], can be used to compute $n\to\chi\gamma$ or $\Lambda\to\chi\gamma$. A
pertinent Feynman diagram is illustrated in Fig. 1. Denoting ${\cal B}$ as
either $n$ or $\Lambda$, the total decay rate for ${\cal B}\to\chi\gamma$ is
given by

$\displaystyle\Gamma\left({\cal B}\to\chi\gamma\right)=\frac{g_{{\cal
B}}^{2}\,e^{2}\,\varepsilon_{{\cal B}\chi}^{2}}{128\pi}\frac{\left(m_{{\cal
B}}+m_{\chi}\right)^{2}}{m_{{\cal B}}^{5}}\left(m_{{\cal
B}}^{2}-m_{\chi}^{2}\right),$ (2.7)

in agreement with Ref. [15]. The mixing parameter $\varepsilon_{{\cal B}\chi}$
follows once the UV model is given, and it is what we constrain through the
analysis of this paper.

To determine the impact of these microscopic processes on the neutron star
requires further model building. Thus far, at low energies we have a dark
baryon $\chi$, which we take to be a massive Dirac fermion. If it is a stable
particle, then it can also be a dark matter candidate. If so, then it may
already exist within the material that collapsed to form the protoneutron
star, though likely only in small amounts, and through dark decays or
adsorption on the star it may accumulate within the star. If it is able to
give up its kinetic energy, then it may settle in the core of the star,
ultimately impacting its properties and evolution.

Figure 2: Illustration of various $\chi$-nucleon ($N$) processes at low
energies, with the heavy black dot denoting the $n\chi\pi^{0}$ effective
vertex noted in the text, namely, (a) $n-N$ scattering to produce $\chi-N$,
(b) $\chi-N$ elastic scattering, and (c) $N-N$ annihilation to produce
$\chi-\chi$. Processes with $\gamma$ in place of $\pi^{0}$ are also possible.
The reverse of reactions (a) and (c) should be strongly suppressed by Pauli-
blocking effects in the interior of a neutron star.

There are many processes in which $\chi$ could participate, though the
interactions with baryons are severely limited by the cold, degenerate nature
of the interior of the neutron star. In principle, given the $n\chi\gamma$ and
$n\chi\pi^{0}$ effective interactions in the models we have noted, and using
$N$ to denote either a neutron or a proton, $\chi$ could (i) be produced via
$nN\to\chi N$ scattering, (ii) interact elastically with another nucleon via a
$nN$ intermediate state, (iii) be formed via the annihilation $nn\to\chi\chi$
or (iv) it can decay via $\chi\to p+e^{-}+{\bar{\nu}}_{e}$ if it is heavy
enough. The reverses of the reactions in (i) and (iii) could also occur.
Pauli-blocking effects associated with the cold, dense neutron medium strongly
suppress all of the reactions in which nucleons appear in the final state.
Moreover, $\chi-N$ elastic scattering is further suppressed in that it occurs
at ${\cal O}(\varepsilon_{n\chi}^{2})$ at amplitude level. We note Fig. 2 for
an illustration. Given this, and our interest in limiting BNV within the star
in a model-independent way, implying that the response of the star to BNV
ought be controlled by SM dynamics, we think that ensuring $\chi$
disappearance is important. Thus we consider two different pathways to do just
that. In the first, we add $\chi$-lepton interactions [5], which intrinsically
break baryon number and are intrinsically very poorly constrained. We would
also want the rate for $\chi$ decay to be no less of that for $\chi$
production. This path, however, is potentially subject to severe constraints
from proton decay experiments. For example, we could have $\chi\to
e^{+}e^{-}\nu$ or $\chi\to 3\nu$, and these channels could give rise to proton
decay via an off-shell $\chi^{*}$ state as in

$p\to\pi^{+}\chi^{*}\to\pi^{+}e^{+}e^{-}\nu\,.$ (2.8)

(Exotic proton decays of just this ilk also emerge in models with quark and
lepton compositeness [52].) Admittedly, this process, as well as the
$p\to\pi^{+}3\nu$ channel, may evade severe constraints due to the particular
nature of existing $|\Delta B|=1$ searches, both because of the final-states
studied and the cuts on the final-state particle momenta needed to control
backgrounds. Although this path could prove to be viable, we favor an
alternate choice: we will allow $\chi$ to decay to other dark particles. A
simple realization of this is given by [14]

${\cal L}_{\rm dark}\supset y_{d}{\bar{\chi}}\phi_{B}\xi+\rm h.c.\,,$ (2.9)

where $\phi_{B}$ is a complex scalar with $B=+1$ and $\xi$ is a Majorana
fermion — and both are dark-matter candidates. Introducing a $Z_{2}$ symmetry,
so that $\chi$, $n$, and $p$ are all $Z_{2}$ even, but $\phi_{B}$ and $\xi$
are $Z_{2}$ odd, we see Eq. (2.9) is the only surviving interaction that
traces to the visible sector, with $n-\xi$ oscillations, say, forbidden by the
$Z_{2}$ symmetry. One interesting consequence of this new path is that dark
decays can be induced in the scattering with either $\phi_{B}$ or $\xi$ in the
initial state, as developed in Ref. [53] and illustrated in Fig. 3. A similar
mechanism, considered in the context of the neutron lifetime anomaly, has been
studied in Ref. [36].

Figure 3: (Color Online) Feynman diagram contributing to induced neutron decay
via a ${\bar{\chi}}\phi_{B}\xi$ interaction, as per Ref. [54] — a $p$ decay
channel follows from the replacement of $d\to u$ in the spectator quark.

The same process can destabilize the proton, with $|\Delta B|=1$ experimental
studies constraining the model parameters [53]. We note that the
${\bar{\chi}}\phi_{B}\xi$ interaction can also induce $\chi\chi$ annihilation,
as noted and illustrated in Fig. 4. The scope of possibilities can be limited
through judicious choices of the parameters of the dark sector. For example,
if $m_{\xi}>m_{n}$ then dark-matter induced nucleon decay will not occur, and
Fig. 4 depicts the only possible tree-level annihilation channel. This last
effect acts to remove $\chi$ produced through neutron decay from the star, yet
$\phi_{B}$ could potentially accumulate in its core — and impact the survival
of the neutron star [55]. If we suppose, rather, that $\phi_{B}$ is light
enough to escape the star, then that outcome can be avoided. We now turn to
the explicit evaluation of proceses that can remove $\chi$ from the neutron
star.

Figure 4: Feynman diagram contributing to $\chi$-$\chi$ annihilation via $\xi$
exchange to yield $B$-carrying scalars, as per the conventions of Fig. 3.
Alternatively, $\chi\chi$ annihilation via $\phi_{B}$ exchange in $t$-channel
would yield a $\xi\xi$ final state, which could ultimately rematerialize as a
$\bar{\chi}\bar{\chi}$ pair.

###  2.1 Dark Baryon Removal Rates

If the masses of $\xi$ and $\phi_{B}$ sum to less than the mass of $\chi$,
then the decay $\chi\to\xi\phi_{B}$ is operative. Using Eq. (2.9) and Refs.
[54, 56], we calculate the width of this decay to be

$\Gamma_{\chi\to\xi\phi_{B}}=\frac{y_{d}^{2}}{16\pi
m_{\chi}^{3}}\left[(m_{\chi}+m_{\xi})^{2}-m_{\phi_{B}}^{2}\right]^{3/2}\left[(m_{\chi}-m_{\xi})^{2}-m_{\phi_{B}}^{2}\right]^{1/2}.$
(2.10)

However, if this decay is operative and if
$m_{\phi_{B}}+m_{\xi}<m_{p}-m_{\pi}$, then this allows for proton decay via
$p^{+}\to\pi^{+}\xi\phi_{B}$. We avoid potentially running afoul of these
constraints by insisting that this decay not be operative and thus require
$m_{\xi}>m_{\chi}$.

Instead, we focus on possible annihilation processes of $\chi$, where we have
assumed that only $\phi_{B}$ is lighter than $\chi$. Adopting the same tools
to compute $\chi\chi\to\phi_{B}\phi_{B}$ we have:

$\displaystyle\sigma_{\chi\chi\to\phi_{B}\phi_{B}}(s)=\frac{y_{d}^{4}m_{\xi}^{2}}{64\pi
s}\Biggl{[}$
$\displaystyle\frac{2\sqrt{(s-4m_{\chi}^{2})(s-4m_{\phi_{B}}^{2})}}{m_{\xi}^{4}+m_{\xi}^{2}(s-2m_{\chi}^{2}-2m_{\phi_{B}}^{2})+(m_{\chi}^{2}-m_{\phi_{B}}^{2})^{2}}$
(2.11)
$\displaystyle+\frac{4}{s+2m_{\xi}^{2}-2m_{\chi}^{2}-2m_{\phi_{B}}^{2}}$
$\displaystyle\quad\times\ln\left(\frac{s+2m_{\xi}^{2}-2m_{\chi}^{2}-2m_{\phi_{B}}^{2}+\sqrt{(s-4m_{\chi}^{2})(s-4m_{\phi_{B}}^{2})}}{s+2m_{\xi}^{2}-2m_{\chi}^{2}-2m_{\phi_{B}}^{2}-\sqrt{(s-4m_{\chi}^{2})(s-4m_{\phi_{B}}^{2})}}\right)\Biggr{]}.$

We note that this cross section goes to zero as $m_{\xi}\to 0$. This must
occur, so that this outcome serves as a non-trivial check of our procedure.
Our cross section result does not depends on whether the scalar is real or
complex, but its interpretation does. If the scalar is real, it cannot carry
baryon-number, and $\chi\chi$ annihilation to scalars would then break B by
two units. This can only occur if $m_{\xi}$ has a nonzero baryon-number-
violating mass. Thus its rate vanishes if $m_{\xi}$ does.

We would like to understand how these annihilation processes operate within a
neutron star. As we will see in Sec. 3.1, these cross sections would need to
be averaged over the true distribution of $\chi$s produced in baryon decays
within the star. Generically, $\chi$s need not be distributed thermally, and
the process of thermalization would require self-interactions, which do not
appear at tree level in our simple model. The problem of $\chi$ transport in
the neutron star is beyond the scope of this paper, so that we assume that the
thermally averaged cross section is a reasonable estimate of what the true
averaged cross section would be.

We proceed by employing pertinent results from the seminal Ref. [57]. The
thermally averaged cross section $\langle\sigma v\rangle$ is given formally by

$\langle\sigma
v\rangle=\frac{1}{8m_{\chi}^{4}T_{\chi}K_{2}^{2}(m_{\chi}/T_{\chi})}\int_{4m_{\chi}^{2}}^{\infty}ds\,\sigma(s)\times(s-4m_{\chi}^{2})\sqrt{s}K_{1}(\sqrt{s}/T_{\chi}),$
(2.12)

where $T_{\chi}$ is the $\chi$ temperature (which is generically nonzero and
may be different from the temperature of the rest of the neutron star) and
$K_{1,2}$ are modified Bessel functions of the second kind. This expression
assumes that it is appropriate to describe the $\chi$ fluid as abiding by a
Maxwell-Boltzmann distribution; it would be inappropriate to apply this
expression to a cold, degenerate population of $\chi$, but such a population
does not occur in our framework. To perform the thermal averaging, we expand
$\sigma(s)\times v$ in powers of $\epsilon\equiv s/(4m_{\chi}^{2})-1$:

$\sigma v=a^{(0)}+a^{(1)}\epsilon+\frac{1}{2}a^{(2)}\epsilon^{2}+\ldots;$
(2.13)

this requires that $v=2\sqrt{\epsilon(1+\epsilon)}/(1+2\epsilon)$. In the
limit in which the $\chi$ fluid is nonrelativistic, the thermally averaged
cross section can we written in terms of the coefficients $a^{(n)}$ as
follows:

$\langle\sigma
v\rangle=a^{(0)}+\frac{3}{2}a^{(1)}\left(\frac{T_{\chi}}{m_{\chi}}\right)+\frac{15}{8}a^{(2)}\left(\frac{T_{\chi}}{m_{\chi}}\right)^{2}+\ldots$
(2.14)

This prescription is expected to be valid as long as $T_{\chi}\lesssim
3m_{\chi}$ [57]. For $\chi\chi\to{\phi_{B}}{\phi_{B}}$, we find the leading-
order contribution to the thermally averaged cross section in $T_{\chi}$ to be

$\langle\sigma
v\rangle_{\chi\chi\to{\phi_{B}}{\phi_{B}}}=\frac{3}{2}\left(\frac{h^{4}m_{\xi}^{2}\sqrt{m_{\chi}^{2}-m_{\phi_{B}}^{2}}}{8\pi
m_{\chi}(m_{\xi}^{2}+m_{\chi}^{2}-m_{\phi_{B}}^{2})^{2}}\right)\left(\frac{T_{\chi}}{m_{\chi}}\right)+\ldots$
(2.15)

Since the $a^{(0)}$ term vanishes, we conclude that the $s$-wave annihilation
contribution vanishes, resulting in a suppression at low temperatures. We
expect our $\chi$s to have a nonzero average kinetic energy from decays, so we
do not expect to encounter a scenario in which these annihilations are
completely quenched by the low energies of their parents, but it is an
interesting feature to note.

We conclude by noting some relevant qualitative features of this model. Since
$\chi$ self-interactions do arise at the one-loop level as a result of
interactions with ${\phi_{B}}$ and $\xi$, we can expect the $\chi$ population
would thermalize, but that timescale is likely slow relative to that of their
annihilation to scalars. There are many more interesting phenomenological
consequences of this model that one could explore, but for our purposes, it is
enough to assume that the masses and coupling conspire such that $\chi$ can be
removed from neutron stars quickly enough that our formalism is valid.

##  3 Macroscopic Baryon Number Violation in Neutron Stars

We set out this section by elaborating the main assumptions for our analysis,
followed by a description of the resulting formalism, which we flesh out in
greater detail than in Ref. [5]. We then discuss the observable effects
associated with our framework, along with methods of interpreting pulsar
observations to yield limits on BNV in such systems. We use the limits derived
at the end of this section to constrain specific baryon dark decay rates in
Sec. 6, though we develop our description of dense matter, as well as of
particle processes within it, in intervening sections before doing so.

###  3.1 Assumptions

The structure of a neutron star can be approximated by a static and
spherically symmetric metric ($g_{\mu\nu}$) with a line element given by [58]

$d\tau^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=e^{2\nu(r)}\,dt^{2}-e^{2\lambda(r)}\,dr^{2}-r^{2}\,d\theta^{2}-r^{2}\,\sin^{2}\theta\,d\phi^{2},$
(3.1)

in which $\nu(r),\lambda(r)$ are solutions to the Einstein field equations
[59], $G^{\mu\nu}=-8\pi GT^{\mu\nu}$, in which $G^{\mu\nu}$ is Einstein’s
tensor, $G$ is the gravitational constant, and $T^{\mu\nu}$ is the stress-
energy tensor. The rotation effects on the neutron star structure, which are
$\mathcal{O}\left(\Omega^{2}/\left(G\,M/R^{3}\right)\right)$ [60], amount to
less than $3\%$ for the fastest spinning pulsar (J1614$-$2230) that we
consider in this work. Furthermore, the inclusion of quasi-static BNV
processes, which are sourced by the matter in the star, would keep the
spherical symmetry intact and changes to the metric ($g_{\mu\nu}$) very slow
in time, such that the use of Eq. (3.1) is warranted.

We also assume that the medium in neutron star can be described by a perfect
fluid with

$T_{0}^{\,0}={\cal E},\quad T_{i}^{\,i}=-P\qquad(i=1,2,3),$ (3.2)

as the only nonzero components of the stress-energy tensor in which $P$ and
${\cal E}$ are the local pressure and energy density of the fluid respectively
which in general depend on the local baryon number density ($n$) and
temperature ($T$) via the EoS. In the standard picture, neutron stars cool
down to internal temperatures $T\lesssim 10^{11}\,{\rm K}\ll E_{\rm
F}\lesssim{\rm GeV}$ within a minute after formation [61], such that the
thermal contribution to the pressure and energy density can be neglected. The
neutron star fluid can then be described as a cold degenerate Fermi gas at
$\beta$–equilibrium. The existing terrestrial constraints on neutron dark
decay, Eq. (2.1), along with the BNV limits we find in Table 1, show that BNV
rates should be slow with respect to other dynamical processes in the neutron
star. We have also devised a model in which $\chi$, the dark baryon-like
particle, can be removed efficiently from the star. Thus we expect the
deviations from a degenerate state at $\beta$–equilibrium due to BNV should be
negligibly small, and we leave a more detailed study of possible thermal
effects on neutron stars from BNV to future work.

In order to be able to apply our model-independent formalism [5], we are going
to focus on a subset of models in which the dark contributions to the EoS are
negligible relative to the energy density and pressure of the visible sector.
In other words, we demand that the following (local) conditions

$\frac{{\cal E}_{\chi}(r)}{{\cal E}(r)}\ll 1,\qquad{\rm
and}\qquad\frac{P_{\chi}(r)}{P(r)}\ll 1,$ (3.3)

hold throughout the neutron star at all times, which can be equivalently
written as a condition on the local number density of $\chi$: $n_{\chi}(r)\ll
n(r)$. This means that $\chi$ has to decay or annihilate either back to the
visible sector or to some other dark particles that can escape the neutron
star. We assume that $\chi$ participates in self-annihilation to lighter dark
particles that can escape the neutron star (see Sec. 2.1 for more details).

We can express the condition $n_{\chi}(r)\ll n(r)$ in terms of the BNV rate,
$\Gamma_{\rm BNV}$, and the annihilation cross section that is averaged over
$\chi$ distribution, which we denote by $\langle\sigma v\rangle$. We note that
the exact distribution of $\chi$ in the neutron star can in principle be found
by solving the Boltzmann transport equation in the star, but this is not
practical for our estimation purposes. We instead consider two scenarios for
$\chi$: one in which the annihilation rate is much faster than the self-
interactions which help establish a thermal equilibrium, and another in which
self-interactions of $\chi$ are much faster than its annihilation rate.

We first consider the scenario in which dark particles have a non-thermal
distribution at the time of their annihilation. If we ignore the effects due
to radial redistribution of $\chi$ after their production and prior to their
annihilation, their number density ($n_{\chi}$) would approximately satisfy

$\dot{n}_{\chi}(t)=n_{i}(t)\times\Gamma_{\rm BNV}-n_{\chi}^{2}(t)\langle\sigma
v\rangle,$ (3.4)

in which $n_{i}(t)$ is the decaying baryon number density which we take to be
constant on short timescales, and $\langle\sigma v\rangle$ is the annihilation
cross section averaged over the distribution of $\chi$. The asymptotic value
for $\chi$ number density (at times $t\gg 1/\sqrt{n_{i}\,\Gamma_{\rm
BNV}\,\langle\sigma v\rangle}$) is then equal to
$n_{\chi}^{\infty}=\sqrt{n_{i}\,\Gamma_{\rm BNV}/\langle\sigma v\rangle}$,
which relative to the local baryon number density $n(r)$ is given by

$\frac{n_{\chi}^{\infty}(r)}{n(r)}=1.4\times
10^{-15}\sqrt{f_{i}(r)}\left(\frac{n_{\rm
sat}}{n(r)}\right)^{1/2}\left(\frac{\Gamma_{\rm BNV}}{10^{-10}\,{\rm
yr}^{-1}}\right)^{1/2}\left(\frac{10^{-26}\,{\rm cm}^{3}\,{\rm
s}^{-1}}{\langle\sigma v\rangle}\right)^{1/2},$ (3.5)

in which $f_{i}(r)\equiv n_{i}(r)/n(r)<1$ is the fraction of baryon $i$
relative to the total baryon number density, $n_{\rm sat}=0.15\,{\rm fm}^{-3}$
is the nuclear saturation density, and we used the scale of the canonical
weak-scale cross section ($10^{-26}\,{\rm cm}^{3}\,{\rm s}^{-1}$) for
comparison. We can see that this ratio is negligible for the reference values
in this equation if $\langle\sigma v\rangle\gg 10^{-56}\,{\rm cm}^{3}\,{\rm
s}^{-1}$.

We can generalize Eq. (3.4) to scenarios in which the redistribution of
$\chi$’s , after their production and prior to their annihilation, is not
negligible, by noting that the total $\chi$-population satisfies

$\dot{N}_{\chi}(t)=B_{i}(t)\times\Gamma_{\rm BNV}-C_{\rm ann}N_{\chi}^{2}(t),$
(3.6)

in which $B_{i}(t)$ is the number of decaying baryons of type $i$, and $C_{\rm
ann}$ is the annihilation rate per particle, such that the total annihilation
rate is identified as $\Gamma_{\rm ann}\equiv C_{\rm ann}N_{\chi}^{2}/2$. We
are interested in short timescales during which $B_{i}(t)$ can be taken as a
constant ($t\ll\Gamma_{\rm BNV}^{-1}$). In this case, the solutions to Eq.
(3.6), assuming $N_{\chi}(0)=0$, are given by

$N_{\chi}(t)=\sqrt{\frac{B_{i}\,\Gamma_{\rm BNV}}{C_{\rm
ann}}}\tanh\left(\sqrt{B_{i}\,\Gamma_{\rm BNV}\,C_{\rm
ann}}\,t\right),\qquad\qquad t\ll\Gamma_{\rm BNV}^{-1}$ (3.7)

in which the timescale for achieving an equilibrium between the production and
annihilation of $\chi$ ($\dot{N}_{\chi}(\tau_{\infty})\approx 0$) can be
identified as $\tau_{\infty}=1/\sqrt{B_{i}\,\Gamma_{\rm BNV}\,C_{\rm ann}}$,
which can be achieved for $\tau_{\infty}<t\ll\Gamma_{\rm BNV}^{-1}$, if

$\frac{\Gamma_{\rm BNV}}{B_{i}}\ll C_{\rm ann}.$ (3.8)

The total number of $\chi$s can then be approximated by its equilibrium value
given by $N_{\chi}^{\infty}=\sqrt{B_{i}\,\Gamma_{\rm BNV}/C_{\rm ann}}$. We
can see that if the condition in Eq. (3.8) holds, then $N_{\chi}^{\infty}\ll
B_{i}$.

We now calculate $C_{\rm ann}$ in the scenario in which the annihilation rate
of $\chi$ is slower than its self-interaction rate, and the $\chi$’s are
distributed spherically with an average radius of $R_{\chi}$, according to
Boltzmann distribution. Using the virial theorem and assuming a radially
uniform distribution of background neutron star matter (over $R_{\chi}$) with
an average energy density $\bar{\cal E}$, we can write

$n_{\chi}(r)=n_{\chi}(0)\,e^{-r^{2}/R_{\chi}^{2}},\quad\quad
R_{\chi}=\sqrt{\frac{3\,k_{B}\,T_{\chi}}{2\pi\,G\,\bar{\cal E}\,m_{\chi}}},$
(3.9)

in which $k_{B}$ is the Boltzmann constant, and $T_{\chi}$ is the dark sector
temperature. The total annihilation rate ($\Gamma_{\rm ann}$) and $N_{\chi}$
can then be evaluated as

$\displaystyle\Gamma_{\rm ann}=$
$\displaystyle\frac{1}{2}\int_{0}^{R_{\chi}}4\pi
r^{2}\left[n_{\chi}(r)\right]^{2}\,\langle\sigma
v\rangle\,dr=0.24\,\left(\frac{k_{B}\,T_{\chi}}{G\,\bar{\cal
E}\,m_{\chi}}\right)^{3/2}\,\langle\sigma v\rangle\,[n_{\chi}(0)]^{2},$ (3.10)
$\displaystyle N_{\chi}=$ $\displaystyle\int_{0}^{R_{\chi}}4\pi
r^{2}n_{\chi}(r)\,dr=0.78\,\left(\frac{k_{B}\,T_{\chi}}{G\,\bar{\cal
E}\,m_{\chi}}\right)^{3/2}\,n_{\chi}(0),$ (3.11)

in which $\langle\sigma v\rangle$ is the thermally averaged annihilation cross
section. Using the definition of $C_{\rm ann}$ we have

$C_{\rm ann}\equiv\frac{2\,\Gamma_{\rm
ann}}{N_{\chi}^{2}}=0.78\,\left(\frac{G\,\bar{\cal
E}\,m_{\chi}}{k_{B}\,T_{\chi}}\right)^{3/2}\,\langle\sigma v\rangle,$ (3.12)

and an equilibrium between the production and annihilation can be achieved on
timescales $t\ll\Gamma_{\rm BNV}^{-1}$, if (see Eq. (3.8))

$\langle\sigma v\rangle\gg 2\times
10^{-55}\left(\frac{10^{57}}{B_{i}}\right)\,\left(\frac{T_{\chi}}{m_{\chi}}\right)^{3/2}\left(\frac{10^{15}\,{\rm
g}/{\rm cm}^{3}}{\bar{\cal E}}\right)^{3/2}\left(\frac{\Gamma_{\rm
BNV}}{10^{-10}\,{\rm yr}^{-1}}\right)\,{\rm cm}^{3}\,{\rm s}^{-1},$ (3.13)

in which $T_{\chi}$ and $m_{\chi}$ have the same units. We can also find the
equilibrium value for $\chi$ number density at the core by combining the
definition of equilibrium number, $N_{\chi}^{\infty}=\sqrt{B_{i}\,\Gamma_{\rm
BNV}/C_{\rm ann}}$, with Eq. (3.11) to arrive at

$\frac{n_{\chi}^{\infty}(0)}{n_{\rm sat}}=8\times
10^{-16}\left(\frac{B_{57}\,\Gamma_{\rm BNV}}{10^{-10}\,{\rm
yr}^{-1}}\right)^{1/2}\left(\frac{10^{-26}\,{\rm cm}^{3}\,{\rm
s}^{-1}}{\langle\sigma v\rangle}\right)^{1/2}\left(\frac{\bar{\cal
E}}{10^{15}\,{\rm g}/{\rm
cm}^{3}}\right)^{3/4}\left(\frac{m_{\chi}}{T_{\chi}}\right)^{3/4},$ (3.14)

in which we defined $B_{57}\equiv B_{i}/10^{57}\sim\mathcal{O}(1)$. Assuming
the reference values in this equation, we can see that this ratio is about
$10^{-13}\,\left({\rm MeV}/T_{\chi}\right)^{3/4}$ for
$m_{\chi}\sim\mathcal{O}({\rm GeV})$. Therefore, we have shown that the self-
annihilation of $\chi$ can be very effective at keeping its concentration
negligible. This concludes the analysis of the necessary conditions on $\chi$
self-annihilation cross section imposed by Eq. (3.3). The explicit forms of
decay and annihilation rates for $\chi$ in terms of our model parameters are
presented in Sec. 2.1.

With the assumptions set forth in this section, we only need to specify the
EoS of hadronic matter to find the neutron star structure. Once the EoS is
specified, the Tolman-Oppenheimer-Volkoff (TOV) [62, 63] equations can be
integrated with the initial conditions $M(0)=0$ and ${\cal E}(0)={\cal E}_{c}$
up to the surface of the star, corresponding to $P(r)=0$. In other words, to
study BNV effects on neutron stars generated by a fixed EoS, we focus on the
unique family of stars, each parameterized by its central energy density
(${\cal E}_{c}$), known as the single-parameter sequence [20] of stars.

The baryon decay rate (per baryon) in a small volume ($V$) in the nuclear
matter (n.m.) rest frame ($\Gamma_{\rm nm}$) is defined by
$d(nV)/d\tau=-\Gamma_{\rm nm}\,n\,V$, in which $\tau$ is the fluid’s proper
time, and $n$ is the proper baryon number density. We can define a baryon
number-flux vector by $j^{\mu}=u^{\mu}n$ [64], in which $u^{\mu}$ is the four-
velocity of the fluid ($u^{\mu}u_{\mu}=1$) and use the definition of
$\Gamma_{\rm nm}$ to write $j^{\mu}_{;\mu}=-n\,\Gamma_{\rm nm}$, in which ‘;’
denotes the covariant derivative. We then use the relationship
$\sqrt{-g}\,j^{\mu}_{;\mu}=\left(\sqrt{-g}\,j^{\mu}\right)_{,\mu}$ [20], in
which ‘,’ denotes ordinary partial derivative and $g\equiv{\rm
det}|g_{\mu\nu}|$, to arrive at

$\frac{\partial}{\partial
t}\int\sqrt{-g}\,n\,u^{0}\,d^{3}x+\int\left(\sqrt{-g}\,n\,u^{i}\right)_{,i}\,d^{3}x=-\int
d^{3}x\sqrt{-g}\,n\Gamma_{\rm nm},$ (3.15)

in which $i=1,2,3$. For a static ($u^{i}=0$), spherically symmetric neutron
star with radius $R$, we can use the metric in Eq. (3.1) to simplify Eq.
(3.15) as

$\frac{\dot{B}}{4\pi}\equiv\frac{\partial}{\partial
t}\left[\int_{0}^{R}\left[1-\frac{2M(r)}{r}\right]^{-\frac{1}{2}}\,n(r)\,r^{2}\,dr\right]=-\int
e^{\nu(r)}\,\left[1-\frac{2M(r)}{r}\right]^{-\frac{1}{2}}\,\Gamma_{\rm
nm}(r)\,n(r)\,r^{2}\,dr,$ (3.16)

where $B$ is the total baryon number of the neutron star. We have used
$\sqrt{-g}=\exp(\nu(r)+\lambda(r))\,r^{2}\,\sin\theta$, with
$\exp(2\lambda(r))=(1-2M(r)/r)^{-1}$, and $M(r)$ is the total mass included
within radius $r$:

$M(r^{\prime})=4\pi\int_{0}^{r^{\prime}}{\cal E}(r)r^{2}dr.$ (3.17)

Given a particle physics model for BNV we can evaluate $\Gamma_{\rm nm}(r)$
and use Eq. (3.16) to find the resulting $\dot{B}$.

###  3.2 Framework

It was shown in Ref. [5] that the conditions in Sec. 3.1 are necessary for a
model-independent analysis of BNV effects on neutron star. These conditions
can be summarised as: (1) BNV is slower than chemical and dynamical responses
in the neutron star, and (2) the contributions to the EoS from any new
particles (e.g., $\chi$) are negligible (see Eq. (3.3)). The overall effect of
BNV within this framework is to relocate the neutron star along its single-
parameter sequence prescribed by the chosen baryon-number-conserving EoS. The
rate of change in any neutron star observable $\mathcal{O}$ as a result of
this quasi-equilibrium evolution can be written as

$\dot{\mathcal{O}}\equiv\left(\frac{d{\cal
E}_{c}}{dt}\right)\left(\frac{\partial\mathcal{O}}{\partial{\cal
E}_{c}}\right).$ (3.18)

Here, we ignore any possible dependence of $\mathcal{O}$ on the angular
velocity $\Omega$, i.e., we assume $\mathcal{O}$ evolves along a one-
dimensional trajectory with $\Omega=0$ on the general two-dimensional space
parameterized by ${\cal E}_{c}$ and $\Omega$. We can solve for $\dot{\cal E}$
in terms of the rate of baryon loss, $\dot{B}$, such that

$\frac{\dot{\mathcal{O}}}{\mathcal{O}}=\left(\frac{B}{\mathcal{O}}\times\frac{\partial_{{\cal
E}_{c}}\mathcal{O}}{\partial_{{\cal
E}_{c}}B}\right)\frac{\dot{B}}{B}\equiv-b(\mathcal{O})\times\Gamma_{\rm BNV},$
(3.19)

in which we defined the effective BNV rate $\Gamma_{\rm BNV}\equiv-\dot{B}/B$
and the dimensionless parameter $b(\mathcal{O})$ encodes the relative rate of
change in $\mathcal{O}$ with respect to $\Gamma_{\rm BNV}$. We pick hadronic
versions of the DS(CMF) EoS [21] that includes a crust [65] from the CompOSE
database [66]. The details of these EoS including their Lagrangians and
particle contents are given in Sec. 4.1. In order to evaluate $b(\mathcal{O})$
factors, we generate a sequence of neutron stars on a grid of ${\cal E}_{c}$
values, and then find the derivative of $\mathcal{O}$ using the central finite
difference method. The resulting $b(\mathcal{O})$ is plotted in Fig. 5 for
various observable quantities as a function of neutron star masses for
DS(CMF)-1 EoS.

Figure 5: (Color online) (a) The set of observable quantities ($O$): mass
$(M)$, radius $(R)$, baryon number $(B)$, and moment of inertia $(I)$ for a
sequence of neutron stars as a function of the central energy density (${\cal
E}_{c}$) relative to their canonical values ($O_{*}$):
$M_{\star}=1.4\,M_{\astrosun}$, $R_{\star}=12\,{\rm km}$, $B_{\star}=10^{57}$,
$I_{\star}=70\,({\rm M_{\astrosun}\,{km}^{2}})$ assuming DS(CMF)-1 EoS. (b)
The $b(\mathcal{O})$ factors for three parameters ($\mathcal{O}=M,R,I$) as a
function of neutron star masses for DS(CMF)-1 EoS. See Sec. 4.1 for more
details on our choice of EoS.

In Sec. 3.3, we use this formalism to show how changes in neutron star
parameters due to BNV would affect pulsar-binary orbital decay rates.

###  3.3 Observables

Baryon loss in pulsars may lead to observable effects on their individual
spin-down rate ($\dot{P}_{s}$), and their orbital period lengthening
($\dot{P}_{b}$) if they belong to a binary system [5]. The BNV modifications
to $\dot{P}_{s}$ are caused by the quasi-equilibrium changes in the moment of
inertia ($I$), and angular momentum loss due to light particles (e.g.,
$\phi_{B}$) escaping the pulsar. While the first contribution can be expressed
in a model-independent manner, the latter depends on the specific BNV model
and the masses of particles involved. Therefore, we focus our attention on BNV
modifications to $\dot{P}_{b}$, which can still be formulated in a model
independent way.

The energy loss due to BNV can modify the orbital period decay rate in a
binary system, assuming it is active in one or both of the components. This
energy loss can be written as [5]

$\dot{M}^{\rm
eff}\equiv\frac{d}{dt}\left(M+\frac{1}{2}I\Omega^{2}\right)=\underbrace{\,b(M)\left(\frac{\dot{B}}{B}\right)M+b(I)\left(\frac{\dot{B}}{B}\right)\left(\frac{2\pi^{2}I}{P_{s}^{2}}\right)\,}_{\textrm{BNV}}-\frac{4\pi^{2}I\dot{P_{s}}}{P_{s}^{3}},$
(3.20)

in which $b(M)$ and $b(I)$ are defined in Eq. (3.19), $P_{s}$ and
$\dot{P_{s}}$ are the observed pulsar spin period and its observed rate of
change respectively. Note that the rates of change in $I$ due to spin-down,
$(dI/d\Omega)\dot{\Omega}$, are negligible in the pulsars that we consider.
The relative rate of change in a binary period due to energy loss in its
components is given by [67, 68, 69]

$\left(\frac{\dot{P}_{b}}{P_{b}}\right)^{\dot{E}}=-2\left(\frac{\dot{M}_{1}^{\rm
eff}+\dot{M}_{2}^{\rm eff}}{M_{1}+M_{2}}\right),$ (3.21)

in which $1$ and $2$ refer to the components of the binary system. After
plugging Eq. (3.20) into (3.21), we get the following BNV and spin-down
contributions to the energy-loss term

$\displaystyle\left(\frac{\dot{P}_{b}}{P_{b}}\right)^{\rm BNV}=$
$\displaystyle\frac{-2}{M_{1}+M_{2}}\sum_{i=1,2}\left(\frac{\dot{B}_{i}}{B_{i}}\right)\left[b_{i}(M)M_{i}+b_{i}(I)\left(\frac{2\pi^{2}I_{i}}{P_{s,i}^{2}}\right)\right],$
(3.22) $\displaystyle\left(\frac{\dot{P}_{b}}{P_{b}}\right)^{\dot{\Omega}}=$
$\displaystyle\frac{8\pi^{2}}{M_{1}+M_{2}}\left(\frac{I_{1}\dot{P}_{s,1}}{P_{s,1}^{3}}+\frac{I_{2}\dot{P}_{s,2}}{P_{s,2}^{3}}\right).$
(3.23)

We should note that the second term in Eq. (3.22), which is due to changes in
the moment of inertia, is $\mathcal{O}(10^{-3})\times b(I)\,M_{\astrosun}$ for
J1614$-$2230 and even smaller for the other two systems considered in this
work. Given that $b(M)\approx 1$ and $|b(I)|\sim\mathcal{O}(1)$, and the spin-
down contributions from Eq. (3.23) are usually subdominant, we conclude that
our limits would be mainly controlled by the first term in Eq. (3.22). For
this reason, our inferred limits on $\dot{B}/B$ are not sensitive to the
specific choices of EoS. We can use the observed pulsar binary period decay
rate to limit the contributions from Eq. (3.22), but first we need to identify
other sources of binary orbital decay.

###  3.4 Interpretation

The dominant contributions to the observed relative rate of orbital period
decay can be written as [70]:

$\left(\frac{\dot{P}_{b}}{P_{b}}\right)^{\rm
obs}=\underbrace{\,\left(\frac{\dot{P}_{b}}{P_{b}}\right)^{\rm
GR}+\left(\frac{\dot{P}_{b}}{P_{b}}\right)^{\dot{E}}\,}_{\textrm{intrinsic}}+\left(\frac{\dot{P}_{b}}{P_{b}}\right)^{\rm
ext},$ (3.24)

in which the first term is due to gravitational radiation [71], and the third
term includes extrinsic effects, e.g., due to the relative motion of a binary
pulsar with respect to the solar system barycenter. The numerical values for
each of these contributions and the limits on $\dot{P}_{b}^{\dot{E}}$ (found
by subtracting the GR contribution, $\dot{P}_{b}^{\rm GR}$, from the intrinsic
orbital-period decay rate, $\dot{P}_{b}^{\rm int}\equiv\dot{P}_{b}^{\rm
obs}-\dot{P}_{b}^{\rm ext}$) are given in Table 1 for three binary systems.
Two of these systems J0348$+$0432 and J1614$-$2230) have heavy pulsars that
may contain hyperons [72], and the third one is a double pulsar system
(J0737$-$3039A/B) with an extremely high precision in its orbital parameters.

  1. 1.

PSR J0348$+$0432: A pulsar–white dwarf binary discovered in 2007 with the
Robert C. Byrd Green Bank Telescope [73] with an orbital period of about $2.4$
hr. We use the results from the analysis in Ref. [74], in which it was shown
that the kinematic, spin-down (Eq. (3.23)), and tidal
($\dot{P}_{b}^{T}\lessapprox 10^{-16}$) contributions to $\dot{P}_{b}$ are
negligible and the observed $\dot{P}_{b}$ should be mainly caused by the GW
emission. We use the value from Ref. [74] for the intrinsic period decay rate,
$\dot{P}_{b}^{\rm int}=-0.275(45)\times 10^{-12}$.

  2. 2.

PSR J1614$-$2230: A pulsar–white dwarf binary discovered in 2006 with the
Parkes radio telescope [75]. We use the Shapiro delay mass estimates from Ref.
[76], and the binary parameters from NANOGrav 12.5 yr data set [77] at $56323$
MJD. The observed value of $\dot{P}_{b}^{\rm obs}=1.57(13)\times 10^{-12}$ is
dominated by the Doppler shift due to the pulsar motion which is itself mainly
caused by the Shklovskii effect [78]:

$\dot{P}_{b}^{\rm Shk}=\frac{\mu^{2}\,d}{c}P_{b}=1.24(9)\times 10^{-12},$
(3.25)

in which we input the value for proper motion $\mu=32.4(5)\,{\rm
mas\,yr^{-1}}$, and used the parallax distance $d=0.65\pm 0.04\,{\rm kpc}$
[79]. We use Eq. (16) from Ref. [80] to estimate the contribution due to the
Galactic potential, namely,

$\left(\frac{\dot{P}_{b}}{P_{b}}\right)^{\rm
Gal}=-\frac{K_{z}|\sin(b)|}{c}-\frac{\Omega_{\astrosun}^{2}R_{\astrosun}}{c}\left(\cos(l)+\frac{\beta}{\beta^{2}+\sin^{2}(l)}\right)\cos(b),$
(3.26)

in which $\beta\equiv(d/R_{\astrosun})\cos(b)-\cos(l)$, $R_{\astrosun}=8.0(4)$
is the Sun’s Galactocentric distance, $\Omega_{\astrosun}=27.2(9)\,{\rm
km}\,{\rm s}^{-1}\,{\rm kpc}^{-1}$ is its Galactic angular velocity, $K_{z}$
is the vertical component of Galactic acceleration approximated by

$K_{z}\left(10^{-9}\,{\rm cm}\,{\rm s}^{-2}\right)\approx 2.27z_{\rm
kpc}+3.68(1-\exp(-4.31\,z_{\rm kpc})),$ (3.27)

for Galactic heights $z\equiv|d\,\sin(b)|\leq 1.5\,{\rm kpc}$. We use the
pulsar’s coordinates $(l,b)=(352.64\degree,20.19\degree)$ to find
$z=0.223(14)\,{\rm kpc}$, and $\dot{P}_{b}^{\rm Gal}=1(5)\times 10^{-14}$.
These extrinsic effects combine to yield $\dot{P}_{b}^{\rm
ext}=\dot{P}_{b}^{\rm Gal}+\dot{P}_{b}^{\rm Shk}=1.25(10)\times 10^{-12}$. Our
resulting estimate for the period derivative, $\dot{P}_{b}^{\rm
int}=0.32(16)\times 10^{-12}$, is positive at $2\sigma$ significance, pointing
to a possible underestimation of extrinsic effects and their errors. However,
we note that if, for example, we instead assume a negligible value for
$\dot{P}_{b}^{\rm int}\approx 0$ and double our error estimates, then we would
still obtain the same limits. We also evaluate the relatively small GW
contribution which for circular orbits is given by [71]

$\begin{split}\dot{P}_{b}^{\rm GW}&=-\frac{192\pi}{5}\left(\frac{2\pi
T_{\astrosun}}{P_{b}}\right)^{5/3}\frac{M_{p}M_{c}}{(M_{p}+M_{c})^{(1/3)}}=-4.17(4)\times
10^{-16},\end{split}$ (3.28)

in which we used the pulsar and white dwarf masses from Ref. [76],
$T_{\astrosun}=4.92549094\times 10^{-6}\,s$, and we neglected the small
eccentricity of the orbit $e=1.333(8)\times 10^{-6}$ [81]. In estimating
$\dot{P}_{b}^{\dot{\Omega}}$ using Eq. (3.23) we assumed the canonical value
$I=10^{45}\,{\rm g\,cm^{2}}$ for the pulsar’s moment of inertia.

  3. 3.

PSR J0737$-$3039A/B: A double pulsar discovered in 2003 [82], comprised of two
radio pulsars (A and B) with pulse periods of $22.7$ ms and $2.8$ ms,
respectively. We use the data from Ref. [83] and the inferred limits on BNV
contributions from Ref. [5].

Name | J0348$+$0432 | J1614$-$2230 | J0737$-$3039A/B
---|---|---|---
$M_{p}\,(M_{\astrosun})$ | $2.01(4)$ | $1.908(16)$ |  $1.338\,185(+12,-14)$ [A]
$M_{c}\,(M_{\astrosun})$ | $0.172(3)$ | $0.493(3)$ |  $1.248\,868(+13,-11)$ [B]
$P_{s}\,({\rm ms})$ | $39.122\,656\,901\,780\,6(5)$ | $3.150\,807\,655\,690\,7$ |  $22.699\,378\,986\,472\,78(9)$ [A]
$\dot{P}_{s}^{\rm obs}\,(10^{-18})$ | $0.240\,73(4)$ | $9.624\times 10^{-3}$ |  $1.760\,034\,9(6)$ [A]
$P_{b}\,({\rm days})$ | $0.102\,424\,062\,722(7)$ | $8.686\,619\,422\,56(5)$ | $0.102\,251\,559\,297\,3(10)$
$\dot{P}_{b}^{\rm obs}\,(10^{-12})$ | $-0.273(45)$ | $1.57(13)$ | $-1.247\,920(78)$
$\dot{P}_{b}^{\rm ext}\,(10^{-12})$ | $1.6(3)\times 10^{-3}$ | $1.25(10)$ | $-1.68(+11,-10)\times 10^{-4}$
$\dot{P}_{b}^{\rm int}\,(10^{-12})$ | $-0.275(45)$ | $0.32(16)$ | $-1.247\,752(79)$
$\dot{P}_{b}^{\rm GR}\,(10^{-12})$ | $-0.258(+8,-11)$ | $-4.17(4)\times 10^{-4}$ | $-1.247\,827(+6,-7)$
$(\frac{\dot{P}_{b}}{P_{b}})^{\dot{E}}_{2\sigma}\,({\rm yr^{-1}})$ | $2.7\times 10^{-10}$ | $2.7\times 10^{-11}$ | $8.3\times 10^{-13}$
$(\frac{\dot{P}_{b}}{P_{b}})^{\dot{\Omega}}\,({\rm yr^{-1}})$ | $<1.4\times 10^{-13}$ | $\approx 4.2\times 10^{-15}$ | $1.04(7)\times 10^{-13}$
$(\frac{\dot{P}_{b}}{P_{b}})^{\rm BNV}_{2\sigma}\,({\rm yr^{-1}})$ | $2.7\times 10^{-10}$ | $2.7\times 10^{-11}$ | $7.3\times 10^{-13}$
$(\frac{\dot{B}}{B})^{\rm BNV}_{2\sigma}\,({\rm yr^{-1}})$ | $1.8\times 10^{-10}$ | $2.0\times 10^{-11}$ | $4.0\times 10^{-13}$
Table 1: The relevant binary parameters for J0348$+$0432 [74], J1614$-$2230
[76, 77], and J0737$-$3039A/B [83]. See the discussion in Sec. 3.4 for more
details.

We can now translate the bounds on $(\dot{P}_{b}/P_{b})^{\rm BNV}$ from Table
1 to limits on $(\dot{B}/B)$ using Eq. (3.22), which are presented in the last
row of Table 1. In deriving these limits, we assumed that BNV is only active
in the pulsars. We also note that we can only infer a model-independent limit
on a linear combination of BNV in pulsars A and B of the double pulsar system
(J0737$-$3039A/B). However, we expect that the rates of BNV (per baryon) would
be about the same in both pulsars, i.e.,
$(\dot{B}_{1}/B_{1})\approx(\dot{B}_{2}/B_{2})$, since their masses are very
close and the composition of light neutron stars ought not change much over
$0.1\,M_{\astrosun}$. In Sec. 6, in which we adopt a specific BNV model ($\cal
B\to\chi\gamma$), our inferred limits on the mixing parameter
($\varepsilon_{{\cal B}\chi}$) are found by evaluating the individual BNV
rates in each of the two pulsars J0737$-$3039A and J0737$-$3039B, which we
then sum to compare to the observational limit on BNV in this system. We also
observe that changing between the DS(CMF) EoSs (see Table 2) induces variation
in, at most, the last significant digit in our limits (see the discussion
below Eq. (3.23)).

##  4 Dense Matter Considerations for Particle Processes

Different lines of evidence reveal that dense matter environments can be
discriminating probes of non-SM processes. For example, limits on
$\Lambda\to\chi\gamma$, as well as other decay channels with dark particles,
follow by noting that the duration of the observed neutrino pulse in SN 1987A
should not be significantly impacted by dark sector emission [16]. We, too,
have found severe limits on BNV from binary pulsar period lengthening, as
shown in Table 1. Here we sharpen such studies by computing particle processes
within a theoretical framework suitable to the description of the dense matter
in the interior of a neutron star.

To compute particle processes in dense matter we might first turn to chiral
effective theory to describe the low-energy interactions of such hadrons [84,
85]. At the simplest level, these studies exploit the symmetries of QCD to
systematize the interactions of mesons and baryons in a momentum expansion in
powers of $(Q/\Lambda_{\chi})$, in which $Q$ is the momentum or pion mass and
$\Lambda_{\chi}$ is the chiral-symmetry breaking scale ($\Lambda_{\chi}\approx
1\,{\rm GeV}$), with experiments fixing the value of the unknown low-energy
constants (LECs) that appear. This framework can also be extended to the
determination of the EoS of neutron stars [86, 87]. The empirical nature of
the LEC determinations limit the applicability of chiral effective theory to
densities no more than $2n_{\rm sat}$ [88]. Moreover, in neutron stars, the
central densities can easily exceed that of saturation density by a factor of
a few, making the nucleons relativistic. As a result, we turn to relativistic
mean-field (RMF) theory in hadronic degrees of freedom to describe the dense
matter at the core of a neutron star. In what follows we first describe how a
RMF treatment emerges from a simple, covariant quantum field theory
description of hadronic interactions before describing the specific chiral
mean-field (CMF) EoS that we employ for generating our numerical results,
showing how this specific choice maps onto the RMF treatment of the simpler
model. We then show how particle decays can be computed within that framework.

###  4.1 Modelling Dense Matter

A prototypical choice is the Walecka model [17, 18, 19], namely,

$\displaystyle{\cal L}_{\varphi/V}$ $\displaystyle=$
$\displaystyle{\bar{\psi}}[(i\not{\partial}-g_{V}\not{V})-(m_{N}-g_{s}\varphi)]\psi+\frac{1}{2}(\partial_{\mu}\varphi\partial^{\mu}\varphi-
m_{s}^{2}\varphi^{2})$ (4.1)
$\displaystyle-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}+\frac{1}{2}m_{v}V_{\mu}V^{\mu}+\delta{\cal
L}\,,$

where $F_{\mu\nu}=\partial_{\mu}V_{\nu}-\partial_{\nu}V_{\mu}$ and a
counterterm $\delta{\cal L}$, as the model is renormalizable. It is similar to
massive QED with a scalar extension and a conserved current (baryon number).
Both a neutral scalar meson ($\varphi$) and a neutral vector meson
($V^{\mu}$), describing the attractive and repulsive features, respectively,
of the nucleon-nucleon force appear. The equations of motions (EoMs) take the
form

$\displaystyle\left(\not{\partial}^{2}+m_{s}^{2}\right)\varphi(x)=g_{s}{\bar{\psi}}\psi$
(4.2)
$\displaystyle\partial_{\nu}F^{\nu\mu}+m_{v}^{2}V^{\mu}=g_{v}{\bar{\psi}}\gamma^{\mu}\psi$
(4.3)
$\displaystyle\left\\{\left[i\not{\partial}-g_{v}\not{V}(x)\right]-\left[m_{N}-g_{s}\varphi(x)\right]\right\\}\psi(x)=0.$
(4.4)

The EoMs are nonlinear and thus complicated. Working in the mean field limit
is grossly simplifying, however. That is, at high baryon number densities, the
sources for $\varphi(x)$ and $V^{\mu}(x)$ fields become large, and these field
operators can be replaced by their vacuum expectation values (VEV) in the n.m.
frame: $\varphi(x)\to\langle\varphi(x)\rangle\equiv\overline{\varphi}$, and
$V_{\mu}(x)\to\langle V_{\mu}(x)\rangle\equiv\delta_{\mu 0}\overline{V}_{0}$.
In doing this, we assume rotational invariance and note that in static uniform
matter, as in a neutron star, $\overline{\varphi}$ and $\overline{V}_{0}$
become constants that only depend on density. The solutions to Eq. (4.4) would
take the form of that of the free Dirac equation if the replacements
$k_{\mu}\to k^{*}_{\mu}\equiv k_{\mu}-g_{v}V_{\mu}$ and $m\to m^{*}\equiv
m-g_{s}\varphi_{0}$ are made. In other words, the medium effects in the RMF
limit are captured by a shift in the baryon momenta and masses. In
generalizing this result for broader use, we note that the Lagrangian of
interactions for a more realistic hadronic model would have more ingredients
(e.g., mesons). However, we would still be able to add up the scalar meson
VEVs that modify the baryon’s mass in a similar manner and denote the
effective baryon mass by $m^{*}$, independent of the specific scalar mesons in
our model. Similarly, we can combine all the contributions to the baryon’s
momentum from vector mesons and denote them by $\Sigma^{\mu}$, such that in
going from the vacuum to the in-medium formalism we would replace $k_{\mu}\to
k^{*}_{\mu}\equiv k_{\mu}-\Sigma_{\mu}$. Equipped with this result, we can
write the wave-function for a baryon with (canonical) four-momentum $k^{\mu}$
(in a uniform medium) as

$\psi(x)=e^{-ik\cdot x}u(k^{*},\lambda),$ (4.5)

in which $k^{*\mu}\equiv
k^{\mu}-\Sigma^{\mu}=\left\\{E^{*}(k^{*}),\vec{k}-\vec{\Sigma}\right\\}$ is
defined to be the kinetic four-momentum and the vector self-energy
($\Sigma^{\mu}$) is generated by the vector meson VEVs, with $\vec{\Sigma}=0$
in the n.m. frame. The time-component of $k^{*\mu}$ is defined by
$E^{*}(k^{*})\equiv\sqrt{{m^{*}}^{2}+|\vec{k}^{*}|^{2}}$, in which $m^{*}$ is
generated by the scalar meson VEVs. The baryon spinor $u(k^{*},\lambda)$
satisfies the Dirac equation

$\left(\not{k}^{*}-m^{*}\right)u(k^{*},\lambda)=0,$ (4.6)

which has the following solution in Dirac-Pauli representation

$u(k^{*},\lambda)=\sqrt{E^{*}(k^{*})+m^{*}}\begin{pmatrix}1\\\
\frac{\vec{\sigma}\cdot\vec{k}^{\,*}}{E^{*}(k^{*})+m^{*}}\end{pmatrix}\chi_{\lambda},$
(4.7)

in which $\vec{\sigma}$ contains the Pauli matrices, and $\chi_{\lambda}$ is
the Pauli spinor with $\chi_{\uparrow}=(1,0)^{T}$ and
$\chi_{\downarrow}=(0,1)^{T}$. Note that $u$ has a Lorentz-invariant
normalization given by $\overline{u}(k^{*},\lambda)u(k^{*},\lambda)=2m^{*}$.
The wave-function for antibaryons can be similarly constructed. The energy
spectrum of baryons ($k^{0}$) is given by

$E(k)=\sqrt{{m^{*}}^{2}+|\vec{k}-\vec{\Sigma}|^{2}}+\Sigma^{0};$ (4.8)

in the mean-field approximation, $\Sigma^{\mu}$ and $m^{*}$ do not depend on
$k^{\mu}$ but they do vary with density. The values for $m^{*}$ and
$\Sigma^{0}$ (in the n.m. frame) decrease and increase respectively (see Fig.
6) in such a way that the total energy of baryons in Eq. (4.8) increases at
higher densities. As we will see shortly, this brings about in-medium baryon
decays to particles that are heavier than the baryon’s vacuum mass since
$E(0)>m_{\cal B}$ at high densities. In general, the increase in the repulsion
between baryons in a RMF framework can be understood by comparing the time-
like component of vector (repulsive) interactions, which are proportional to
$u^{\dagger}u$, with scalar (attractive) interactions, which are parameterized
by $\overline{u}u=(m^{*}/E^{*})u^{\dagger}u$. As the density increases,
$m^{*}$ decreases and the strength of the attractive forces relative to the
repulsive ones diminishes [89]. However, we should note that having a highly
repulsive nuclear interaction at extremely high densities (compared to $n_{\rm
sat}$) is a reasonable expectation, regardless of the specific dense matter
formalism. Having explained the formalism utilized in this work, we now
describe the specific EoS that we use for generating our numerical results.

We choose an EoS based on a non-linear hadronic SU($3$) CMF model [90], in
which the baryonic degrees of freedom include nucleons ($n$, $p$), hyperons
($\Lambda$, $\Sigma$, $\Xi$) and the spin-3/2 resonances ($\Delta$,
$\Sigma^{*}$, $\Xi^{*}$, $\Omega$). These baryons interact via exchange of
scalar ($\sigma$, ${\bm{\delta}}$, $\zeta$, $\chi$) and vector mesons
(${\bm{\rho}}^{\mu}$, $\omega^{\mu}$, $\phi^{\mu}$), in which
${\bm{\rho}}^{\mu}$ and ${\bm{\delta}}$ are both isovectors. In the RMF limit,
the mesons become classical fields, and in the n.m. frame only the zeroth
components of vector mesons develop VEV. The Lagrangian density of the CMF
model is given by [21]

$\mathcal{L}=\mathcal{L}_{\rm Kin}+\mathcal{L}_{\rm Int}+\mathcal{L}_{\rm
Self}+\mathcal{L}_{\rm SB},$ (4.9)

in which $\mathcal{L}_{\rm Kin}$ contains the usual kinetic terms for baryons
and leptons, $\mathcal{L}_{\rm Int}$ is due to the baryon-meson interactions
which are given by

$\mathcal{L}_{\rm
Int}=-\sum_{i}\overline{\psi}_{i}\left(g_{i\omega}\gamma_{0}\langle\omega^{0}\rangle+g_{i\phi}\gamma_{0}\langle\phi^{0}\rangle+2g_{i\rho}\gamma_{0}I_{3i}\langle\rho_{3}^{0}\rangle+m_{i}^{*}\right)\psi_{i}\,.$
(4.10)

We note $\psi_{i}$ denotes a baryon of species $i$ with an effective mass
$m_{i}^{*}$ and an isospin 3-component $I_{3i}$, and the expectation value is
evaluated in the ground state. The last two terms in Eq. (4.9), i.e.,
$\mathcal{L}_{\rm Self}$ and $\mathcal{L}_{\rm SB}$, contain the self-
interactions of scalar and vector mesons and explicit chiral symmetry breaking
terms respectively. The explicit expressions are given in Eqs. (3), (4), and
(5) of Ref. [21]. The baryon effective masses are generated by the scalar
meson VEVs, except for a small explicit mass term $\delta m_{i}\sim 150$ MeV,
and are given by

$m_{i}^{*}=g_{i\sigma}\langle\sigma\rangle+2g_{i\delta}\langle\delta_{3}\rangle
I_{3i}+g_{i\zeta}\langle\zeta\rangle+\delta m_{i},$ (4.11)

in which $\delta_{3}$ is the isospin 3-component of ${\bm{\delta}}$. The time-
component of baryon self-energy is given by

$\Sigma^{0}_{i}=g_{i\omega}\langle\omega^{0}\rangle+g_{i\rho}\langle\rho^{0}_{3}\rangle
I_{3i}+g_{i\phi}\langle\phi^{0}\rangle.$ (4.12)

The numerical values for $m^{*}$ and $\Sigma^{0}$ are plotted in Fig. 6. We
note that the reduction of the effective baryon masses at high densities as
shown in Fig. 6 is due to chiral symmetry restoration at high densities. The
coupling constants are chosen [91, 92, 93] to reproduce the hadron vacuum
masses, the nuclear saturation properties (density $n_{\rm sat}=0.15\,{\rm
fm}^{-3}$, binding energy per nucleon $B/A=-16.00\,{\rm MeV}$, compressibility
$K=300\,{\rm MeV}$), the asymmetry energy ($E_{\rm sym}=30\,{\rm MeV}$), and
hyperon potentials. Furthermore, the pion and kaon decay constants constrain
the scalar meson VEVs.

This conventional approach in determining the coupling constants in RMF models
relies on an extrapolation from symmetric finite nuclei to infinite neutron
matter. We would like to contrast this with an alternative that we may wish to
employ in the future, which is based on fitting uniform pure neutron matter
properties determined through the use of chiral effective field theory [94].
The latter procedure involves fitting the RMF couplings with the synthetic
neutron matter data generated using Quantum Monte Carlo (QMC) many-body
methods [95], in addition to reproducing $n_{\rm sat}$, $B/A$, and $K$.

Our chosen class of EoS has a set of variations that depend on the degrees of
freedom that are included, and they are given in Table 2 for convenience. Thus
in order to explore how our results vary with EoS, we employ the choices given
there.

EoS | DoF | Add. Int. |  $L$ (MeV) | $U_{\Lambda}\,({\rm MeV})$ | $U_{\Sigma}\,({\rm MeV})$ | $U_{\Xi}\,({\rm MeV})$ | $U_{\Delta}\,({\rm MeV})$ | $M_{\rm max}\,(M_{\astrosun})$
---|---|---|---|---|---|---|---|---
$1$ | N$+$Y | – | $88$ | $-28$ | $5$ | $-18$ | – | $2.07$
$2$ | N | – | $88$ | $-28$ | $5$ | $-18$ | – | $2.13$
$3$ | N$+$Y | $\omega\rho$ | $75$ | $-28$ | $5$ | $-18$ | – | $2.00$
$4$ | N | $\omega\rho$ | $75$ | $-28$ | $5$ | $-18$ | – | $2.05$
$5$ | N$+$Y | $\omega\rho+\omega^{4}$ | $75$ | $-27$ | $6$ | $-17$ | – | $2.07$
$6$ | N | $\omega\rho+\omega^{4}$ | $75$ | $-27$ | $6$ | $-17$ | – | $2.11$
$7$ | N$+$Y$+\Delta$ | $\omega\rho+\omega^{4}$ | $75$ | $-27$ | $6$ | $-17$ | $-64$ | $2.07$
$8$ | N$+\Delta$ | $\omega\rho+\omega^{4}$ | $75$ | $-27$ | $6$ | $-17$ | $-64$ | $2.09$
Table 2: The set of CMF EoS variants taken from Ref. [92, 93]; we refer to
them as DS(CMF)-1 through DS(CMF)-8, respectively, in later use. The second
and third columns describe the degrees of freedom (DoF); nucleons (N),
hyperons (Y), and delta resonances ($\Delta$); and the additional vector
interactions (“Add. Int.”) beyond the standard terms ($\mathcal{L}_{\rm
Self}$) that are included for each EoS respectively. The fourth column
represents the assumed value for symmetry energy ($E_{\rm sym}$) slope ($L$).
The fifth to eighth columns are the single-particle hyperon potentials, and
the last column is the maximum neutron star mass ($M_{\rm max}$) that can be
generated.

The set of EoSs that we utilize has also been extended to include crusts based
on a zero-temperature unified EoS [65] at $\beta$–equilibrium with similar
values of the symmetry energy slope ($L$), in which the effective interactions
are Skyrme forces Rs [96] (EoS 1-2) and SkMP [97] (EoS 3-8) with cluster
energy functionals taken from Ref. [98].

Figure 6: (Color Online) The effective masses (a), vector self-energies (b),
and the energy of baryons at rest in the n.m. frame (c) as a function of
density in the DS (CMF)-1 EoS. The horizontal lines correspond to the vacuum
masses of baryons, and the vertical lines indicate the central number density
($n_{c}$) of the pulsars we consider in this work.

###  4.2 Medium Effects: Effective Masses and Beyond

In this section, we discuss some of the notable features that emerge in
studying processes in the medium, and make comparisons with the vacuum
formalism. We start with the quantization of baryon fields in the medium
followed by the rate and cross section calculation formalism. We then discuss
the electromagnetic form factors of the baryons that are needed for our
calculations in Sec. 5.2.

The presence of the baryon Fermi sea modifies the quantization procedure of
the baryon fields, $\psi(x)$, in medium [18] compared with the usual procedure
in vacuum [99]. Once the coefficients behind Fourier modes of $\psi(x)$ are
promoted to baryon creation ($a^{\dagger}(k)$) and annihilation ($a(k)$)
operators (likewise $b^{\dagger}(k)$ and $b(k)$ for antibaryons), we conclude
that the action of these operators on the medium ground state
$|\Omega\rangle$, which contains baryon levels filled to a Fermi momentum
($k_{F}$), should be given by

$\begin{split}b(k)|\Omega\rangle=&0\quad\forall\vec{k},\\\
a^{\dagger}(k)|\Omega\rangle=&0\quad|\vec{k}|<k_{F},\\\
a(k)|\Omega\rangle=&0\quad|\vec{k}|>k_{F}.\end{split}$ (4.13)

This leads to a different form (compared to vacuum) for the baryon propagator
which is given by [18]

$G(p)\equiv\left(\not{p}^{*}+m^{*}\right)\Bigg{\\{}\frac{1}{p^{*2}-m^{*2}+i\varepsilon}+2\pi
i\delta(p^{*2}-m^{*2})\theta(p^{*}_{0})\,\theta\Big{[}k_{F}^{2}+p^{*2}-\frac{\big{(}p^{*}_{\mu}B^{\mu}\big{)}^{2}}{B_{\mu}B^{\mu}}\Big{]}\Bigg{\\}},$
(4.14)

in which $\theta$ is the Heaviside step function, $B^{\mu}$ is the baryon
current density, which in the n.m. frame is given by $B^{\mu}_{nm}=\delta^{\mu
0}n_{B}$, and the second term in Eq. (4.14) allows for the propagation of
holes in the Fermi sea. Using this modified propagator and the spinors in Eq.
(4.7), one can derive Feynman rules [18] for calculating the amplitudes for
various processes (see Sec. 5.2). However, in calculating rates via phase
space integrals, we should first observe that an on-shell ($p^{*2}=m^{*2}$)
and positive energy ($p^{0}>0$) Lorentz-invariant integral over the four-
momentum is given by

$\begin{split}\int
d^{4}p\,\delta\left(\left(p^{\mu}-\Sigma^{\mu}\right)^{2}-m^{*2}\right)\theta(p^{0})f(p^{\mu})=&\int
d^{3}p\,\frac{f\left(\sqrt{|\vec{p}-\vec{\Sigma}|^{2}+m^{*2}}+\Sigma^{0},\vec{p}\right)}{2\sqrt{|\vec{p}-\vec{\Sigma}|^{2}+m^{*2}}}.\end{split}$
(4.15)

Therefore, we identify the Lorentz-invariant (on-shell) volume element in the
medium as $d^{3}p/2E^{*}(p)$. This means that the normalization factors in the
in-medium phase space integrals should contain $\left(2E^{*}\right)^{-1}$ in
place of the usual vacuum expression.

We also note that the velocity of a baryon is defined in terms of the kinetic
momentum as opposed to the canonical one, i.e., $v_{\mu}\equiv
k^{*}_{\mu}/E^{*}$. This velocity should be used for calculating the cross
section of two-body scattering involving a baryon (see App. C). We can
explicitly show this by performing an integration over the longitudinal
($\hat{z}$) components of the incident beams’ momenta
($\bar{k}^{z}_{\mathcal{A}}$ and $\bar{k}^{z}_{\mathcal{B}}$). Let us assume
for the moment that only one baryon ($\mathcal{B}$) is involved, in which case
we have (see Eq. (4.77) of Ref. [100])

$\begin{split}&\int
d\bar{k}^{z}_{\mathcal{A}}d\bar{k}^{z}_{\mathcal{B}}\,\delta\Big{(}\bar{k}^{z}_{\mathcal{A}}+\bar{k}^{z}_{\mathcal{B}}-\sum_{f}p_{f}^{z}\Big{)}\delta\Big{(}\bar{E}_{\mathcal{A}}+\bar{E}_{\mathcal{B}}-\sum_{f}E_{f}\Big{)}\\\
&=\int
d\bar{k}^{z}_{\mathcal{A}}\,\delta\Big{(}\sqrt{\bar{k}^{2}_{\mathcal{A}}+m_{\mathcal{A}}^{2}}+\sqrt{\big{(}\bar{k}_{\mathcal{B}}-\vec{\Sigma}_{\mathcal{B}}\big{)}^{2}+m_{\mathcal{B}}^{2}}+\Sigma^{0}_{\mathcal{B}}-\sum_{f}E_{f}\Big{)}\bigg{|}_{\bar{k}^{z}_{\mathcal{B}}=\sum
p_{f}^{z}-\bar{k}^{z}_{\mathcal{A}}}\\\
&={\left|\frac{\bar{k}^{z}_{\mathcal{A}}}{\bar{E}_{\mathcal{A}}}-\frac{\bar{k}^{z}_{\mathcal{B}}-\Sigma^{z}_{\mathcal{B}}}{\bar{E}^{*}_{\mathcal{B}}}\right|}^{-1}={\left|\frac{\bar{k}^{z}_{\mathcal{A}}}{\bar{E}_{\mathcal{A}}}-\frac{\bar{k}^{*z}_{\mathcal{B}}}{\bar{E}^{*}_{\mathcal{B}}}\right|}^{-1}\equiv{\left|v_{\mathcal{A}}-v_{\mathcal{B}}\right|}^{-1},\end{split}$
(4.16)

in which in the last line we are assuming $\bar{k}^{z}_{\mathcal{B}}=\sum
p_{f}^{z}-\bar{k}^{z}_{\mathcal{A}}$ and have identified the baryon velocity
using the kinetic momentum, such that
$\left|v_{\mathcal{A}}-v_{\mathcal{B}}\right|$ is the relative velocity of the
beams as viewed from the laboratory frame. The generalization to the case with
two baryons is straightforward. The fact that the velocity of a baryon is zero
when $\vec{k}^{\,*}=0$ could have also been deduced by inspecting the kinetic
energy component in Eq. (4.8). For this reason, the frame in which
$\vec{k}^{\,*}=0$ holds is called the center of velocity (c.v.) frame which is
distinct from the center of mass (c.m.) frame defined by $\vec{k}=0$.
Therefore, the decay rate of a baryon in an arbitrary frame ($\Gamma$) is
found by boosting ($\gamma$) the rate evaluated in the c.v. frame using
$\Gamma=\gamma^{-1}\,\Gamma_{\rm c.v.}$.

Since we study processes that involve electromagnetic (EM) interactions with
baryons, the generalization of EM form factors from the vacuum to within the
medium should be checked. The in-medium spinors in Eq. (4.7) are different
from their vacuum counterparts. Therefore, certain commonly used properties
(e.g., Gordon decomposition) in vacuum need to be reestablished. However, we
note that the general form of these interactions is determined by the
structure of Dirac algebra. While important for formulating our analyses, this
is slightly tangential to the broader narrative of this work; we thus relegate
the details to the Appendices, but we encourage the reader to study them
nonetheless. In App. A, we explicitly show that the vacuum EM vertex form can
be generalized to its in-medium form if one replaces $m\to m^{*}$, $p\to
p^{*}$, and identifies the electric charge and magnetic moment of a baryon
from the scattering amplitudes in the c.v. frame. Our numerical results in
Sec. 5 assume the vacuum values for the in-medium form factors $F_{1,2}^{*}$
of neutron and $\Lambda$. We also derive the non-relativistic limit of
baryon’s EM interactions and their elastic scattering formalism in App. B. We
present the calculations for in-medium Compton scattering in App. C, as a
demonstration of the RMF formalism utilized in this work.

##  5 Baryon Dark Decay Rates in Dense Matter

In this section, we develop the procedures for evaluating particle physics
processes, such as neutron decays and neutron-neutron scattering, in the
neutron-star medium. Our particular interest is in radiative decays such as
${\cal B}\to\chi\gamma$ in the core of the star. In the absence of a matter
environment, a common procedure, adopted in many contexts, is to assume the
mixing is weak and to redefine the fields, here $B^{i}$ and $\chi$ [101], so
that they no longer mix, and then to analyze $B^{i}\to\chi$ transitions in
that new basis. In Sec. 5.1, we show why and how this procedure can fail in
strongly interacting dense matter, and we argue for a Feynman diagram analysis
in its place. Subsequently, starting in Sec. 5.2, we show how the transition
rates can be evaluated explicitly and consider their implications.

###  5.1 General Considerations

To illuminate the essential points, we consider the possibility of $n$-$\chi$
mixing in a background field $\Sigma^{\mu}$, the vector self-energy of a
neutron in the neutron-star medium, which interacts with the neutron field
$\psi_{n}$, but not the $\chi$ field $\psi_{\chi}$. Thus we adopt the
following simple model:

$\mathcal{L}=\overline{\psi}_{n}\left(i\not{\partial}-\not{\Sigma}-m_{n}^{*}\right)\psi_{n}+\overline{\psi}_{\chi}\left(i\not{\partial}-m_{\chi}\right)\psi_{\chi}-\varepsilon\left(\overline{\psi}_{n}\psi_{\chi}+\overline{\psi}_{\chi}\psi_{n}\right)\,.$
(5.1)

Under a field redefinition, $\psi\to\psi^{\prime}$, prescribed by

$\begin{pmatrix}\psi_{n}^{\prime}\\\
\psi_{\chi}^{\prime}\end{pmatrix}=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}\psi_{n}\\\
\psi_{\chi}\end{pmatrix},$ (5.2)

Eq. (5.1) becomes

$\begin{split}\mathcal{L}^{\prime}=&\overline{\psi}_{n}^{\prime}\left(i\not{\partial}-\not{\Sigma}\cos^{2}\\!\theta-
m_{n}^{*}\cos^{2}\\!\theta-
m_{\chi}\sin^{2}\\!\theta\right)\psi_{n}^{\prime}\\\
&+\overline{\psi}_{\chi}^{\prime}\left(i\not{\partial}-\not{\Sigma}\sin^{2}\\!\theta-
m_{\chi}\cos^{2}\\!\theta-m_{n}\sin^{2}\\!\theta\right)\psi_{\chi}^{\prime}\\\
&+\overline{\psi}_{n}^{\prime}\left[\frac{\sin(2\theta)}{2}\left(m_{n}^{*}-m_{\chi}+\not{\Sigma}\right)-\varepsilon\cos(2\theta)\right]\psi_{\chi}^{\prime}\,.\end{split}$
(5.3)

If $\Sigma^{\mu}$ were absent, and with $\varepsilon$ real, then for
$\tan\left(2\theta\right)={2\varepsilon}/({m_{n}^{*}-m_{\chi}})$, $\cal
L^{\prime}$ describes two decoupled fields with a modified energy spectrum.
These fields can then map to the asymptotic (“in” and “out”) states needed to
define the $S$-matrix [102]. To do this, any interactions with these fields
should vanish as $t\to\pm\infty$. For the neutron (and other SM baryons), we
note that the effect of the vector self-energy can be absorbed into the
definition of a modified single-particle spinor, as discussed in Sec. 4.1, and
thus suitable “in” and “out” states can still be constructed. In the current
case, $\Sigma^{\mu}$ mediates an interaction between the rotated $n$ and
$\chi$ fields, putting the utility of our field redefinition procedure into
question. After all, even in the mean-field limit, $\Sigma^{0}$ can greatly
exceed the $n$ and $\chi$ masses at the high densities reached within a
neutron star, and it cannot vanish as $t\to\pm\infty$, since we work within a
medium of infinite extent. Since $\Sigma^{\mu}$ is not a Lorentz scalar, we
cannot extend our field redefinition approach to include it. Therefore, there
would seem to be no advantage to following a field redefinition approach in
neutron matter. Moreover, in the small mixing limit
($\varepsilon\ll|m_{n}^{*}-m_{\chi}|$), the mass ($n^{\prime}$,
$\chi^{\prime}$) and interaction ($n$, $\chi$) eigenstates are nearly the
same. Working with Eq. (5.1), we can treat
$\varepsilon\overline{\psi}_{n}\psi_{\chi}$ as a tiny interaction that
mediates $n\leftrightarrow\chi$ transitions within perturbation theory. This
Feynman diagram analysis, through the in-medium baryon propagator, Eq. (4.14),
naturally includes the impact of momentum dependence and of the neutron self-
energy on $n$-$\chi$ mixing. We emphasize that both effects are absent in the
field redefinition procedure. As a result, too, we do not have large
enhancements in our predictions should the in-medium neutron and $\chi$ states
become degenerate in energy — the imaginary part of the neutron self-energy
effectively eliminates that possibility. Nevertheless, $n$-$\chi$ mixing
within the neutron-star medium could potentially lead to effects not possible
in terrestrial experiments, and we consider those possibilities more carefully
in Sec. 5.3.

###  5.2 Dark Decay Rate Estimates

We now turn to the explicit evaluation of rates of particle processes within
the neutron-star medium, with a particular focus on dark decay rates. As long
known, the background field associated with matter leads to a spontaneously
breaking of Lorentz symmetry, but as a consequence of our Lorentz covariant
description, discussed in Sec. 4.1, our expressions always have definite
Lorentz transformation properties. In what follows, we exploit our freedom to
choose a frame to simplify our analysis.

Generally, processes of the form ${\cal B}+\\{X\\}\to\chi+\\{Y\\}$ lead to the
following rate of change of the local baryon density $n_{B}$ (with respect to
the proper time, $\tau$, referenced to that spacetime point):

$\displaystyle\dfrac{dn_{B}}{d\tau}$ $\displaystyle=-\int d\Pi_{\cal
B}\biggl{(}\prod_{\\{X\\}}d\Pi_{X}\biggr{)}d\Pi_{\chi}\biggl{(}\prod_{\\{Y\\}}d\Pi_{Y}\biggr{)}$
(5.4) $\displaystyle\quad\qquad\times f_{\cal B}(\vec{p}_{\cal
B})\biggl{(}\prod_{\\{X\\}}f_{X}(\vec{p}_{X})\biggr{)}\Bigl{(}1-f_{\mathcal{\chi}}(\vec{k}_{\chi})\Bigr{)}\biggl{(}\prod_{\\{Y\\}}\Bigr{[}1\pm
f_{Y}(\vec{k}_{Y})\Bigr{]}\biggr{)}$
$\displaystyle\quad\qquad\times|\mathcal{M}|^{2}\times(2\pi)^{4}\delta^{(4)}\biggl{(}p_{\cal
B}+\sum_{\\{X\\}}p_{X}-k_{\chi}-\sum_{\\{y\\}}k_{Y}\biggr{)},$

where $\\{X\\}$ ($\\{Y\\}$) is some set of other states in the initial (final)
state — which may be empty. Moreover,
$d\Pi_{i}=d^{3}\vec{p}_{i}/[(2\pi)^{3}(2E^{*}_{i})]$ is the Lorentz-invariant
phase space measure, $f(\vec{p})$ are the species-dependent occupation
numbers444The occupation factor for $Y$ depends on whether or not $Y$ is a
boson ($+$) or fermion ($-$)., and $|\mathcal{M}|^{2}$ is the _spin-summed_
(as opposed to _spin-averaged_) squared matrix element. We denote final-state
momenta with $k_{i}$ instead of $p_{i}$. Consistent with our assumption that
there is no appreciable background of $\chi$, we set its occupation factor
$f_{\mathcal{\chi}}(\vec{k}_{\chi})$ to zero. All baryonic species abide by
zero-temperature Fermi distributions characterized by distinct Fermi momenta
$p_{F,{\cal B}}$.

We briefly discuss important qualitative features of the evaluation of Eq.
(5.4) for the decay process ${\cal B}\to\chi\gamma$ and present the
corresponding results. We relegate details of the calculation to App. E. We
work in the interaction basis, so that the decay proceeds via the Feynman
diagram containing the $n-\chi$ interaction and the baryon magnetic dipole
moment operator, which we write as

$\mathcal{O}_{{\cal B}\gamma}=\frac{g_{\cal B}e}{8m_{\cal
B}^{*}}\overline{{\cal B}}\sigma^{\mu\nu}{\cal B}F_{\mu\nu},$ (5.5)

noting $g_{n}=3.826$ and $g_{\Lambda}=-1.226$ [50]. This computation is made
in a background mean-field of neutron matter, and the associated decay
amplitude, as developed in Sec. 4.1, is determined by replacing the canonical
momenta of the in-vacuum computation with kinetic momenta as per Eq. (4.5).
Labeling canonical momenta as ${\cal B}(p_{\cal
B})\to\chi(k_{\chi})+\gamma(k_{\gamma})$, the corresponding spin-summed
squared matrix element is

$|\mathcal{M}|^{2}=\frac{\varepsilon_{{\cal B}\chi}^{2}g_{\cal
B}^{2}e^{2}}{2(m_{\cal B}^{*})^{2}}\left[(p_{\cal B}^{*}\cdot
k_{\chi})+m_{\cal B}^{*}m_{\chi}\right]\,,$ (5.6)

noting that both the kinetic momentum $p^{*}_{\cal B}$ and the in-medium mass
$m_{\cal B}^{*}$ appear. It then remains to insert this into Eq. (5.4) and
integrate. The integral takes the form

$\frac{dn_{\cal B}}{d\tau}=-\int\frac{d^{3}\vec{p}_{\cal
B}}{(2\pi)^{3}(2E^{*}_{\cal
B})}\frac{d^{3}\vec{k}_{\chi}}{(2\pi)^{3}(2E_{\chi})}\frac{d^{3}\vec{k}_{\gamma}}{(2\pi)^{3}(2E_{\gamma})}f_{\cal
B}(\vec{p}_{\cal
B})\times|\mathcal{M}|^{2}\times(2\pi)^{4}\delta^{(4)}\left(p_{\cal
B}-k_{\chi}-k_{\gamma}\right)\,,$ (5.7)

and it can be computed in different ways. The Lorentz-invariance of each
measure $d\Pi_{i}$ affords the opportunity of performing different parts of
the integration in different frames. We note that the integration over
$d\Pi_{\cal B}$ is simplest in the n.m. frame: there, we have $f_{\cal
B}(\vec{p}_{\cal B}^{\,\rm(n.m.)})=\Theta(p_{F,{\cal B}}-|\vec{p}_{\cal
B}^{\,\rm(n.m.)}|)$, and the integrand is isotropic. Contrariwise, it is
simplest to evaluate the integration over $d\Pi_{\chi}~{}d\Pi_{\gamma}$ in the
c.v. frame. Moreover, as discussed in Sec. 4.2, the width of an individual
baryon is most simply interpreted in its respective c.v. frame, since the
baryon is not moving. The baryon width in the c.v. frame takes the form

$\Gamma_{\rm c.v.}(|\vec{p}_{\cal B}^{\,\rm(n.m.)}|)=\frac{1}{2m_{\cal
B}^{*}}\int\frac{d^{3}\vec{k}_{\chi}}{(2\pi)^{3}(2E_{\chi})}\frac{d^{3}\vec{k}_{\gamma}}{(2\pi)^{3}(2E_{\gamma})}\times\frac{1}{2}|\mathcal{M}|^{2}\times(2\pi)^{4}\delta^{(4)}\left(p_{\cal
B}-k_{\chi}-k_{\gamma}\right)\,,\,$ (5.8)

where the argument of $\Gamma_{\rm c.v.}$ follows from our earlier frame
choice.555The additional factor of $\frac{1}{2}$ in Eq. (5.8) arises because
$|\mathcal{M}|^{2}$ has been spin-summed and not spin-averaged. Henceforth we
abbreviate $p_{B}\equiv|\vec{p}_{\cal B}^{\,\rm(n.m.)}|$. We have $p_{\cal
B}^{*,\rm{(c.v.)}}=(m_{\cal B}^{*},0)$ and $p_{\cal
B}^{\rm(c.v.)}=(m_{n}^{*}+\Sigma_{\cal B}^{\rm{(c.v.)},0},\vec{\Sigma}_{\cal
B}^{\rm(c.v.)})$, with $\Sigma_{\cal B}^{\rm(c.v.)}$ the baryon vector self-
energy in the c.v. frame. The results in different frames are connected by
Lorentz boosts, yielding

$\Gamma_{\rm c.m.}(p_{B})=\left(\frac{m_{\cal B}^{*}}{E_{\cal
B}^{*,\rm{(c.m.)}}}\right)\Gamma_{\rm c.v.}(p_{B})=\left(\frac{E_{\cal
B}^{*,\rm{(n.m.)}}}{E_{\cal B}^{*,\rm{(c.m.)}}}\right)\Gamma_{\rm
n.m.}(p_{B})\,.$ (5.9)

The total rate of baryon loss in the n.m. frame is then given by integrating
over all baryons in the local fluid, accounting for the contraction of their
individual widths by a factor $\gamma^{-1}=m_{\cal B}^{*}/E_{\cal B}^{*}$,
with $E_{\cal B}^{*}=\sqrt{|\vec{p}_{\cal B}^{\,\rm(n.m.)}|^{2}+(m_{\cal
B}^{*})^{2}}$ the n.m.-frame energy associated with the kinetic momentum:

$\dfrac{dn_{B}}{d\tau}=-2\times\int_{0}^{p_{F,B}}\frac{p_{B}^{2}dp_{B}}{2\pi^{2}}\gamma^{-1}\Gamma_{\rm
c.v.}(p_{B})\,.$ (5.10)

The prefactor of 2 comes from the baryon’s two spin degrees of freedom, and
the factor of $2\pi^{2}$ in the denominator ensures that if the function
$\gamma^{-1}\Gamma_{\rm c.v.}(p_{B})$ were a constant, then the result would
be $(\gamma^{-1}\Gamma_{\rm c.v.})\times n_{B}$. Combining Eqs. (5.8) and
(5.10), we find that this procedure leads to the same result as the direct
evaluation of Eq. (5.7) shown in App. E. In either case, we arrive at the
following result:

$\displaystyle\dfrac{dn_{\cal B}}{d\tau}=-\frac{\varepsilon_{B\chi}^{2}g_{\cal
B}^{2}e^{2}}{128\pi^{3}}(m_{\cal B}^{*})^{2}\int_{1}^{x_{F}}$ $\displaystyle
dx\,\sqrt{x^{2}-1}\times\frac{1+2x\sigma+\sigma^{2}-\mu^{2}}{(1+2x\sigma+\sigma^{2})^{2}}$
$\displaystyle\times\left[(1+2x\sigma+\sigma^{2})(1+x\sigma+2\mu)+\mu^{2}(1+x\sigma)\right]\,,$
(5.11)

in which

$\displaystyle x\equiv\frac{E_{\cal B}^{*,\rm{(n.m.)}}}{m_{\cal
B}^{*}},\quad\quad x_{F}\equiv\frac{E_{F,{\cal B}}^{*,\rm{(n.m.)}}}{m_{\cal
B}^{*}},\quad\quad\sigma\equiv\frac{\Sigma_{\cal B}^{\rm{(n.m.)},0}}{m_{\cal
B}^{*}},\quad\quad\mu\equiv\frac{m_{\chi}}{m_{\cal B}^{*}}\,.$ (5.12)

The corresponding c.m.-frame single baryon decay rate $\Gamma_{\rm
c.m.}(p_{B})$ is given in Eq. (E.18).

Figure 7: (Color Online) The proper rates for ${\cal B}\to\chi\gamma$ decays,
$(-dn_{\cal B}/d\tau)$ (in fm-3 s-1), for neutrons (a) and $\Lambda$s (b)
assuming DS(CMF)-1 EoS. In either panel, we fix the corresponding
$\varepsilon_{B\chi}$ to be $10^{-16}$ MeV. Note that the color scales are
different between the two panels. See text for additional details. Figure 8:
(Color Online) The baryon proper decay (${\cal B}\to\chi\gamma$) rate
$(-dn_{\cal B}/d\tau)$ (per unit volume) assuming $\varepsilon_{{\cal
B}\chi}=10^{-16}$ MeV for (a) neutrons and (b) $\Lambda$s as a function of
radius for four pulsars using DS(CMF)-1 EoS. The numbers next to curves
(colors) indicate different values of $m_{\chi}=\\{0.5,1.0,1.2,1.3\\}\,{\rm
GeV}$.

We illustrate the rates of ${\cal B}\to\chi\gamma$ as a function of $m_{\chi}$
and $n$ in Fig. 7, for both neutrons (a) and $\Lambda$s (b). This calculation
is for the DS(CMF)-1 EoS, but the results are qualitatively similar for the
other EoS in this family.666Of course, the EoS that do not contain hyperons
will not lead to $\Lambda\to\chi\gamma$ decays within neutron stars. The
vertical axes have been normalized to the value of nuclear saturation density
in this EoS, $n_{\rm sat}=0.15$ fm-3. The respective color scales are shown at
right, assuming $\varepsilon_{{\cal B}\chi}=10^{-16}$ MeV; the units are fm-3
s-1 and we emphasize that these rates scale as $\varepsilon_{{\cal
B}\chi}^{2}$. To guide the eye, we have also added black contours every
quarter order of magnitude. Solid contours correspond to integer numbers; dot-
dashed contours correspond to half-integer numbers; and dotted lines
correspond to quarter-integer numbers. In either panel, the dashed vertical
line indicates the vacuum mass of the corresponding baryon. The dotted
horizontal line corresponds to the central density of the heaviest stable
neutron star within this EoS, corresponding to $M_{\rm TOV}\approx
2.07\,M_{\odot}$; the region above this line has been grayed out because these
densities do not occur in a stable neutron star. Similarly, the dot-dashed
horizontal lines correspond to the central densities of neutron stars with the
masses of J0348+0432, J1614–2230 and J0737–3039A/B. We observe that when
$n\to\chi\gamma$ is operative, it is almost always numerically larger than the
rate of $\Lambda\to\chi\gamma$ (for
$\varepsilon_{n\chi}=\varepsilon_{\Lambda\chi}$). This is a simple consequence
of larger neutron number fractions at these densities, and the two rates often
differ by several orders of magnitude. However, $\Lambda$s have a further
reach in $m_{\chi}$ when they are present than neutrons do, owing to the
larger total energy of $\Lambda$s in neutron matter.

###  5.3 Medium-Enabled Dark Decay Processes

It was shown in Sec. 4.1 that baryons in neutron stars have a lower effective
mass ($m_{\cal B}^{*}$) and a higher self-energy ($\Sigma^{0}_{\cal B}$) at
higher densities (see Fig. 6), but their overall energy can be much higher
than their vacuum rest mass ($m_{\cal B}$). In order to illustrate this for a
heavy neutron star, we plot the baryon rest-energies ($E_{\cal B}^{0}\equiv
E_{\cal B}(p=0)$ in the n.m. frame) for PSR J0348$+$0432 as a function of
radius in Fig. 9. We can see that baryon decays containing a final state
$\chi$ with $m_{\chi}>m_{{\cal B}}$, which would be forbidden in vacuum, can
occur at the core of heavy neutron stars. This enables a novel way of
analyzing models with $m_{\chi}$ values for which nuclear and vacuum decays
are kinematically forbidden. Furthermore, constraints derived from heavy
neutron stars can still be applicable in the vicinity of $m_{\chi}\approx
m_{\cal B}$ and beyond that. This should be contrasted with limits derived
from processes in vacuum and within nuclei, which diminish at $m_{\chi}\approx
m_{\cal B}$ or even at much lower values of $m_{\chi}$ due to the binding
energy and possible energy cuts on the final states. For example, when
inferring limits from $n\to\chi\gamma$ via detection of $\gamma$ there is an
energy cut $E_{\gamma}^{\rm min}$ [103], which means $m_{\chi}$ values larger
than $m_{n}-E_{\gamma}^{\rm min}$ cannot be constrained.

Figure 9: (Color Online) The particle composition (a), and the energy of
baryons at rest in the n.m. frame (b) as a function of radius in PSR
J0348$+$0432 assuming the DS (CMF)-1 EoS. The horizontal lines correspond to
the vacuum masses of baryons. We choose models for which we expect the steady-
state admixture of dark states to be completely negligible.

#### 1\. Spontaneous ${\cal B}\to\chi$ Conversion

The existence of $\chi$ raises the possibility that the baryons to which they
couple might undergo spontaneous conversion to $\chi$ in the neutron-star
medium as they propagate. Such an effect could prove loosely analogous to
empirically observed matter-enhanced neutrino oscillations [104] or to the
possibility of neutron-antineutron oscillations [105, 106, 107], breaking
baryon number by two units. In the latter case the presence of external
interactions from matter or magnetic fields modify the energy of the $n$ and
${\bar{n}}$ differently, severely reducing the spontaneous oscillation
probability for a fixed source of new physics [108], and the cross-section for
scattering-mediated $n$-${\bar{n}}$ conversion is also very small [109]. In
this section, we note the distinct features of ${\cal B}$-$\chi$ conversion.

The essential physics is thus: ${\cal B}$ and $\chi$ constitute a two-level
quantum system. As we have noted in Sec. 5.1, if the coupling
$\varepsilon_{{\cal B}\chi}$ is nonzero, then ${\cal B}$ and $\chi$ constitute
the interaction basis, whereas the eigenstates of this Hamiltonian, which we
term $f_{1}$ and $f_{2}$ for this discussion, constitute the mass basis.
Formally, the strong interactions that operate in neutron matter only ever
produce $n$ — this is what it means for ${\cal B}$ to be an interaction
eigenstate. This ${\cal B}$ is, however, a coherent superposition of $f_{1}$
and $f_{2}$ at the moment of its creation. The subsequent evolution of this
coherent wavepacket depends on the details of the ${\cal B}-\chi$ system.
These details are discussed in depth in App. D; we pick out the most relevant
results as they pertain to this discussion.

The Hamiltonian that describes our two-state system depends on the local
environment: the total energy of the baryon depends on the density through
$m_{\cal B}^{*}$ and $\Sigma_{\cal B}$, and baryons with different n.m.-frame
momenta will mix differently with $\chi$ because Lorentz invariance is
spontaneously broken by the background. There exists a _resonance_ in this
system wherever the condition, which follows from energy-momentum conservation
of the canonical momenta,

$\sqrt{m_{\chi}^{2}+|\vec{p}^{\,\rm(n.m.)}|^{2}}\approx\sqrt{(m_{\cal
B}^{*})^{2}+|\vec{p}^{\,\rm(n.m.)}|^{2}}+\Sigma_{\cal B}^{0,{\rm(n.m.)}}$
(5.13)

is satisfied. We expect that this condition will occur for at most one value
of the (magnitude of the) baryon momentum for a given density. Moreover, Eq.
(5.13) cannot be satisfied if $\Sigma_{\cal B}^{0,{\rm(c.v.)}}$ is _complex_.
In what follows, we set this latter possibility aside, because, as we will
see, other effects act to suppress the likelihood of ${\cal B}$-$\chi$
conversion. At resonance, the offset between the interaction and mass bases is
maximal, corresponding to a mixing angle of $\theta=45^{\circ}$; however, if
Eq. (5.13) is violated by more than a few times $\varepsilon_{{\cal B}\chi}$,
then the mixing angle is parametrically small: $\theta\sim\varepsilon_{{\cal
B}\chi}^{2}/(\omega^{(+)}_{\cal
B}-\omega^{(+)}_{\chi})^{2}\equiv(\varepsilon_{{\cal
B}\chi}/\delta\omega^{(+)})^{2}$, where $\omega^{(+)}_{\chi}$ and
$\omega^{(+)}_{\cal B}$ are the left- and right-hand sides of Eq. (5.13),
respectively (cf. Eqs. (D.14) and (D.15)).

We first consider what happens when the system is not close to resonance. In
this case, the eigenvalues of the system are given by Eqs. (D.25) and (D.27),
which are very nearly given by $\omega^{(+)}_{\chi}$ and $\omega^{(+)}_{\cal
B}$ up to $\mathcal{O}(\varepsilon_{{\cal B}\chi}^{2})$ corrections. If the
system is far from resonance, then these eigenvalues are well separated. As a
result, the ${\cal B}$ states produced in scattering processes will
essentially immediately decohere into their component $f_{1}$ and $f_{2}$
with, respectively, probabilities of $\cos^{2}\theta$ and $\sin^{2}\theta$. As
such, the state that emerges from the scattering process manifests as either
$f_{1}$ with probability $\cos^{2}\theta\sim 1$ or $f_{2}$ with probability
$\sin^{2}\theta\sim(\varepsilon_{{\cal B}\chi}/\delta\omega)^{2}$, and the
latter may be vanishingly small — and thus so would be any yield in $\chi$.

The situation is richer if the state is close to resonance. In this case, if
the canonical momentum of the baryon is fixed by Eq. (D.17), then the
eigenvalues of the Hamiltonian are better given by Eq. (D.32). They are nearly
identical to each other, but they are split by a small factor,

$\Delta\omega^{*}=2\varepsilon_{{\cal B}\chi}\Sigma_{\cal B}^{0,({\rm
n.m.})}\sqrt{\frac{\left(m_{\cal B}^{*}+m_{\chi}\right)^{2}-(\Sigma_{\cal
B}^{0,({\rm n.m.})})^{2}}{\left[(m_{\cal
B}^{*})^{2}-m_{\chi}^{2}\right]^{2}-(\Sigma_{\cal B}^{0,({\rm
n.m.})})^{4}}}+\mathcal{O}(\varepsilon_{{\cal B}\chi}^{3})\,.$ (5.14)

This means that when ${\cal B}$ is produced in some strong interaction, the
wavepacket containing $f_{1}$ and $f_{2}$ may remain coherent over relatively
long timescales. This is analogous to how neutrino mass eigenstates remain
coherent as they propagate in terrestrial oscillation experiments, despite
being formed in an interaction eigenstate.777We point the interested reader to
Ref. [110] for a discussion about the role of (de)coherence in understanding
neutrino oscillations specifically, as well as Ref. [111] for a comparative
analysis of neutrino oscillations with adiabatic conversion. As in the case of
neutrino oscillations, the $f_{1}$ and $f_{2}$ components of the ${\cal B}$
state generically evolve with different phases; over time, this leads to
nonzero overlap between the evolved state and either $\cal{B}$ or $\chi$. The
state is then measured, in a sense, at its next interaction some time $t$
later, either by its environment or by some experimental apparatus. It is
appropriate, in this case, to invoke the concept of an oscillation
probability; this is estimated by

$P_{{\cal
B}\to\chi}(t)=\sin^{2}2\theta\times\sin^{2}\left[\frac{(\Delta\omega^{*})t}{2}\right]\,.$
(5.15)

When the state is observed, however, it collapses to the combination of
$f_{1}$ and $f_{2}$ appropriate to either ${\cal B}$ or $\chi$ with
probabilities given by Eq. (5.15), and the process repeats for further
interactions. While the oscillations have a large amplitude
($\sin^{2}2\theta\sim\mathcal{O}(1)$) in this regime, the probability to
convert will remain small if the time between successive measurements $\delta
t_{\rm meas}$ is small, in the sense $(\Delta\omega^{*})(\delta t_{\rm
meas})\ll 1$. This is precisely the quantum Zeno effect [112, 113].

It remains to determine the timescale of the interactions in the nuclear
medium in order to estimate the rate of ${\cal B}\to\chi$ conversions. We
estimate this to be the light time of the mean interparticle separation around
nuclear saturation density: $\delta t_{\rm strong}\sim n^{-1/3}_{\rm
sat}c^{-1}\sim\mathcal{O}(10^{-23})$ s. For a benchmark value
$\varepsilon_{{\cal B}\chi}=10^{-16}$ MeV, the argument of the latter sine
function in Eq. (5.15) is $\sim\mathcal{O}(10^{-39})$ MeV s
$\sim\mathcal{O}(10^{-18})$; this is safely approximated as small, and we see
that the quantum Zeno effect is indeed operative under these conditions.
Therefore, even if the mixing angle is large, we estimate the probability to
be

$P_{B\to\chi}(\delta t_{\rm
strong})\sim\mathcal{O}(10^{-36})\times\left(\frac{\varepsilon_{B\chi}}{10^{-16}{\rm\,MeV}}\right)^{2}\times\left(\frac{\delta
t_{\rm strong}}{10^{-23}{\rm\,s}}\right)^{2}.$ (5.16)

This implies an approximate per-baryon conversion rate of

$R_{B\to\chi}(\delta t_{\rm
strong})\sim\mathcal{O}(10^{-13}){\rm\,s^{-1}}\times\left(\frac{\varepsilon_{B\chi}}{10^{-16}{\rm\,MeV}}\right)^{2}\times\left(\frac{\delta
t_{\rm strong}}{10^{-23}{\rm\,s}}\right).$ (5.17)

One might expect that this would multiply the large density of baryons to
yield a macroscopically relevant rate. However, the near-resonance region
occupies a thin shell (parametrically of width $\sim\varepsilon_{{\cal
B}\chi}$) within the baryon Fermi sphere; the fraction of baryons relevant for
this phenomenon is fantastically small, even in the best case scenario. Thus
we summarize by emphasizing that we do not expect ${\cal B}-\chi$ conversion
to be a phenomenologically relevant mechanism for the production of $\chi$.

Figure 10: (Color Online) The volume-integrated rates for ${\cal
B}\to\chi\gamma$ decays, $(-d{\cal B}/dt)$ (in s-1) for neutrons (a) and
$\Lambda$s (b) assuming the DS(CMF)-1 EoS. As in Fig. 7, we have fixed
$\varepsilon_{{\cal B}\chi}=10^{-16}$ MeV. The right panels (c & d) show the
relationship between the neutron star mass and its central density.  Figure
11: (Color Online) The total baryon loss rate per baryon $(-\dot{{\cal
B}}/{\cal B})$ due to ${\cal B}\to\chi\gamma$ decays assuming
$\varepsilon_{{\cal B}\chi}=10^{-16}$ MeV for neutron and $\Lambda$ as a
function of $m_{\chi}$ for four pulsars using the DS(CMF)-1 EoS.

###  5.4 Total Rates

In this section, we report the total baryon decay rates that emerge after
integrating our earlier results over the structure of a neutron star with a
given central density, $n_{\rm c}$. For example, in Fig. 10, we show the rates
that result from integrating the local BNV rates in Fig. 7 over the neutron
star volume using Eq. (3.16) and report these results as a function of
$m_{\chi}$ and $n_{c}$. Panel (a) is for neutron decays, while panel (b) is
for $\Lambda$ decays; similar to Fig. 7, the contours correspond to constant
(base-ten log of the) integrated rate of $\cal{B}\to\chi\gamma$. We have again
fixed $\varepsilon_{{\cal B}\chi}=10^{-16}$ MeV, and note that the results are
in s-1. We have coded the black contours in the same way as in Fig. 7, and we
have again indicated the central densities of J0348+0432, J1614–2230, and
J0737–3039A/B within this EoS. The right panels, (c) and (d), contextualize
these results by showing the neutron star mass, $M_{\rm NS}$, on the
horizontal axis as a function of the central density on the vertical axis.
Note that Figs. 7 and 10 together imply that J0737–3039A/B are both too light
to contain hyperons.

##  6 Inferred Limits on Baryon Dark Decays

Figure 12: (Color Online) Exclusion constraints at 2$\sigma$ on the $n-\chi$
(upper panel) and $\Lambda-\chi$ (lower panel) mixing parameters as functions
of $m_{\chi}$ for the DS(CMF)-1 EoS. The blue, red and green curves correspond
to the $2\sigma$ limits derived from J0348+0432, J1614–2230, and
J0737–3039A/B, respectively. The dashed black curves correspond to the
combined limits, realized as per the discussion in the text. The vertical
dashed lines indicate the in-vacuum baryon mass in each case.

We now turn to the task of assessing the limits on the ${\cal B}-\chi$ mixing
parameters that emerge from our numerical assessment of the stellar-volume-
integrated baryon dark decay rates, as shown in Fig. 11, and the macroscopic
baryon number loss limits we have determined from astrophysical observations
and their analysis. The latter, namely, are limits on anomalous binary-pulsar
period lengthening, to which we refer as “binary spin-down,” and they are
given in Table 1. We show the limits we find for each astrophysical system as
well as that associated with a final combined limit. To make our presentation
more compact, we first discuss how the individual limits on
$\varepsilon_{{\cal B}\chi}$ can be combined before showing all of these
results. Note, too, that since our constraint depends on the square of
$\varepsilon_{{\cal B}\chi}$ that its sign is left unconstrained — we choose
$\varepsilon_{{\cal B}\gamma}>0$ in reporting our limits.

### Combining Individual Limits

Here we briefly describe our statistical procedure for combining limits on
$\varepsilon_{{\cal B}\chi}$ derived from different pulsar binary systems. The
limits we show have implicitly been determined as contours of constant
$\chi^{2}(m_{\chi},\,\varepsilon_{{\cal B}\chi})$. Our assumed-true hypothesis
is that rate of BNV-induced binary spin-down vanishes in these systems, so we
have $\chi^{2}=0$ for $\varepsilon_{{\cal B}\chi}=0$. As such, each $\chi^{2}$
function is generically of the form

$\chi^{2}(m_{\chi},\varepsilon_{{\cal
B}\chi})=\left(\frac{\dot{P}_{b}}{P_{b}}(m_{\chi},\,\varepsilon_{0})\right)^{2}\frac{(\varepsilon_{{\cal
B}\chi}/\varepsilon_{0})^{4}}{\sigma^{2}}\equiv
F(m_{\chi})\times\varepsilon_{{\cal B}\chi}^{4}\,.$ (6.1)

The first equality follows from the fact that
$\dot{P}_{b}/P_{b}\propto\varepsilon_{{\cal B}\chi}^{2}$, noting Eq. (3.22),
and we emphasize that $F$ is a function of $m_{\chi}$ only. The limits we have
shown correspond to $\chi^{2}=c$;888For two degrees of freedom, the $2\sigma$
exclusion curves we have shown correspond to $c=6.18$. we call the resulting
curve $\widetilde{\varepsilon}(m_{\chi})$. From this, we determine

$F(m_{\chi})=\frac{c}{\widetilde{\varepsilon}(m_{\chi})^{4}};$ (6.2)

this allows to determine the $\chi^{2}$ function over the entire parameter
space.

The combined limit, then, corresponds to the contour along which the sums of
the individual $\chi^{2}$ functions also equals $c$. Using the definitions
above, we determine the combined limit $\widetilde{\varepsilon}_{\rm
comb}(m_{\chi})$ as follows:

$\displaystyle\chi^{2}_{\rm comb}(m_{\chi},\varepsilon_{{\cal
B}\chi})=\sum_{i}F_{i}(m_{\chi})\varepsilon_{{\cal
B}\chi}^{4}=\sum_{i}c\left(\dfrac{\varepsilon_{{\cal
B}\chi}}{\widetilde{\varepsilon}_{i}(m_{\chi})}\right)^{4}$ (6.3)
$\displaystyle\chi^{2}_{\rm comb}\bigl{(}m_{\chi},\widetilde{\varepsilon}_{\rm
comb}(m_{\chi})\bigr{)}=c\implies\widetilde{\varepsilon}_{\rm
comb}(m_{\chi})=\Bigl{(}\sum_{i}\widetilde{\varepsilon}_{i}(m_{\chi})^{-4}\Bigr{)}^{-1/4}$
(6.4)

This discussion has assumed that all $\widetilde{\varepsilon}_{i}$ are defined
at the same level $c$, and that the desired combined limit is also at $c$.
This result can be generalized for distinct individual significances $c_{i}$
and combined significance $C$:

$\widetilde{\varepsilon}_{\rm
comb}=\Bigl{(}\sum_{i}c_{i}/C\times\widetilde{\varepsilon}_{i}(m_{\chi})^{-4}\Bigr{)}^{-1/4}.$
(6.5)

We show our individual pulsar limits as well as our combined limits, realized
via our described procedure, for the DS(CMF)-1 EoS in Fig. 12.

Figure 13: (Color Online) Combined exclusion constraints at 2$\sigma$ on the
$n-\chi$ (upper panel) and $\Lambda-\chi$ (lower panel) mixing parameters as
functions of $m_{\chi}$ for the eight equations of state in the DS(CMF)
family.

Figure 13 depicts our results for the constraints on $\varepsilon_{n\chi}$
(upper panel) and $\varepsilon_{\Lambda\chi}$ (lower panel) as functions of
$m_{\chi}$, calculated for each of the eight EoS in the DS(CMF) family.
Equations of state that do not include hyperons are indicated with dashed
curves in the upper panel. We also note that the DS(CMF)-3 EoS formally cannot
support a neutron star with a mass of 2.01$M_{\odot}$ — its maximum TOV mass
is 2.00$M_{\odot}$. However, this is within 1$\sigma$ of the observed mass of
J0348+0432; we therefore elect to include it in this figure, but instead show
the constraint derived for this maximal neutron star. This constraint has been
shown in dot-dashing to indicate that it is qualitatively different from the
others.

We underscore that we have fixed the masses of these neutron stars to their
best-fit values to construct these limits. A more statistically complete
analysis would propagate the uncertainty in the inferred masses of the
observed pulsars into the determinations of their central densities (within
the context of a given EoS), and thus into the predicted baryon loss and
binary spin-down rates. The mass uncertainties on J0348+0432 and J1614–2230
are $\mathcal{O}(10^{-2})M_{\odot}$; we anticipate that there would be
$\mathcal{O}(1)$ corrections to the limits whenever these are the only
operative constraints, though the orders of magnitude are expected to be
correct. That said, the uncertainties on the masses of J0737–3039A/B are
$\mathcal{O}(10^{-5})M_{\odot}$ [83], so that we expect the limits on
$\varepsilon_{n\chi}$ for $m_{\chi}\lesssim 1000$ MeV, noting Fig. 12, to be
quite robust.

Figure 14: (Color Online) Exclusion limits at 2$\sigma$ on the vacuum
branching fraction for ${\cal B}\to\chi\gamma$ for neutrons (upper panel) and
$\Lambda$s (lower panel). The results for each EoS are color-coded as in Fig.
13. We also show constraints from KamLAND [114], SuperKamiokande [115], and
BESIII [116].

In Fig. 14, we reinterpret our constraints on $\varepsilon_{{\cal B}\chi}$ as
constraints on the branching fractions for ${\cal B}\to\chi\gamma$ in vacuum
and contrast them against laboratory constraints, with neutrons ($\Lambda$s)
in the upper (lower) panel. For neutrons, we also show the KamLAND constraint
on invisible neutron decay [114] in red and the SuperKamiokande constraint on
$n\to\nu\gamma$ [115] in light blue. These are as much as twenty orders of
magnitude stronger than the constraints we have derived, but we note that
these are only operative up to $m_{\chi}$ = 920 and 827 MeV, respectively.
This is a result of experimental cuts — heavier $\chi$s result in less
energetic photons in the decay, and eventually these become too soft to be
meaningfully detected. We emphasize, in particular, that these experiments
cannot probe the region $m_{\chi}>m_{n}$; while they are more powerful when
they are operative, they are fundamentally constrained in ways that
astrophysical probes of new physics are not. For $\Lambda$s, we show the
constraint on invisible decays from BESIII [116] in dark cyan. In this case,
we find the opposite result: pulsar binaries are able to probe this branching
ratio as much as twenty orders of magnitude more severely than laboratory
constraints! The caveat is that this requires hyperons to appear in neutron
stars, which is still a matter of debate, simply because EoSs without hyperons
exist that confront current observational data successfully. However, if
hyperons appear in an appreciable amount in these objects, then one can expect
vast improvements on laboratory searches.

Figure 15: (Color Online) Exclusion limits at 2$\sigma$ on the vacuum
branching fraction for $n\to\chi\gamma$ as per Fig. 14, in the particular
$\chi$ mass region pertinent to an explanation of the neutron lifetime
anomaly. Additional constraints and expected limits have been included as
detailed in the text, after Ref. [34].

The upper panel of Fig. 14 is incomplete in that there are additional
constraints around $m_{\chi}\approx m_{n}$, a region that has become of
interest in recent years as a result of tests of new-physics explanations [10]
of the neutron lifetime anomaly [12]. We examine this region more closely in
Fig. 15; panel (a) casts these searches in terms of constraints on
$\varepsilon_{n\chi}$, while panel (b) casts them in terms of constraints on
Br$(n\to\chi\gamma)$. We show in blue the estimated constraint from a direct
search for $n\to\chi\gamma$ using ultra-cold neutrons (UCN) [103], and in
green we show a constraint from Borexino from searches for hydrogen decay,
both from Ref. [34]. We also show the curve along which the free hydrogen
lifetime is supposed to be $\tau_{H}=10^{32}$ s in dashed gold, also from Ref.
[34]. (The constraints from Ref. [34] are reported at 90% CL, though the
differences between those and limits at $2\sigma$ should be very small given
the ranges shown in the figure.) Clearly, neutron stars are more sensitive to
these decays than these (would-be) laboratory constraints by many orders of
magnitude.

It was noted in Ref. [10] that the existence of $\chi$ can destabilize nuclear
matter, including 9Be. This constraint was calculated more precisely in Ref.
[49], assuming that the lifetime of 9Be is longer than $3\times 10^{9}$ years
to account for the presence of 9Be in old, metal-poor stars [117]. This
constraint is shown in red in Fig. 15 and is competitive with (if not dominant
to) our neutron star constraints in the region of its operation,
$m_{\chi}<937.993$ MeV. We note that other probes of dark decays of nuclei
with low neutron separation energies have been discussed in, e.g., Ref. [47].
Particular attention has been paid to decays of 11Be, with experimental
efforts underway at CERN-ISOLDE [118] and ISAC-TRIUMF [119], though we are
unaware of any efforts to interpret these experimental results as constraints
on new physics. As a side note, it is curious that there are no laboratory
constraints, as far as we can tell, on the lifetime of 9Be. We find the
arguments about the presence of 9Be in old stars compelling and agree that
this is a valid constraint, but we are surprised, frankly, that the lifetime
is only constrained at the billion-year scale. While experimentalists of yore
would have had little reason to interrogate the stability of 9Be – or indeed,
any species thought to be stable in the SM – we regard the observation that
the stability of these systems has not been tested in a detailed way in the
laboratory as a potentially promising avenue for constraining new physics.

We conclude by noting that Ref. [49] has also presented constraints on
$n\to\chi\gamma$ from cosmology and from neutron star cooling. The former is a
combination of constraints coming from modifications to Big Bang
Nucleosynthesis (BBN) and the Cosmic Microwave Background (CMB); this
treatment includes the reverse decay $\chi\to n\gamma$ when $m_{\chi}>m_{n}$,
and so constrains the region shown. However, in their calculations, $\chi$ is
assumed to constitute (at least some of) the dark matter. This is unlike our
framework, in which we introduced more new states ($\xi$ and $\phi_{B}$) to
prevent overaccumulation of $\chi$. Therefore, the limits they derive from BBN
and CMB do not apply here, though we agree that this would be an interesting
and important avenue to explore. The neutron star cooling constraint derived
there, however, makes very rough assumptions about how heat from decays is
deposited into the neutron star, with the implicit assumption that increases
in the temperature of the core of the neutron star lead to commensurate
increases in the observed effective temperature. Yet thermal transport and
cooling in neutron stars demands careful investigation; for instance, BNV
decays lead to $\beta$-disequilibrium, which leads to neutrino cooling via
(direct and modified) Urca processes, which impact how the energy released in
the decays is deposited back into the SM fluid. While we agree that while old,
cold neutron stars should constrain this model, the details are intricate and
expected to be sufficiently impactful that we decline to include such
constraints here.

##  7 Implications for Models of Baryogenesis and Dark Matter

The prospect of explaining the origins of both the dark matter abundance and
the cosmic baryon asymmetry within a single dynamical framework is a beguiling
one. Different possibilities have existed for some time, and many share a
common feature: there is a dark-sector baryon that carries baryon number and
into which SM baryons can decay. A particularly intriguing variant is that of
$B$ mesogenesis [14, 29, 120]. It proceeds in the early universe from late
time, out of equilibrium production of $B$ mesons (with equal fractions of $b$
and ${\bar{b}}$ quarks) that evolve under SM CP-violating processes before
decaying to a SM baryon and a dark fermion carrying the opposite sign of
baryon number. Thus: no new sources of CP violation of the SM are required;
the baryon number of the universe is conserved — it is just sequestered into
visible and dark sectors with opposite baryon number; and it occurs late in
the history of the universe in that occurs after the QCD phase transition,
making it possible to realize hadronic states, and before the epoch of big-
bang nucleosynthesis. Finally, it is an example of a testable mechanism of
baryogenesis [25], in that its essential features are subject to direct
experimental investigation. Particularly, its reliance on the SM mechanism of
CP violation (albeit new CPV sources could enter) implies that the branching
ratios of $B$ mesons in SM baryons and the dark fermion (antibaryon) cannot be
too small, with the expectation that the branching fractions can roughly be no
less than ${\rm Br}(B_{s,d}^{0}\to\chi{\cal B})\gtrsim 10^{-5}$ or ${\rm
Br}(B^{+}\to\chi{\cal B}^{(+)})\gtrsim 10^{-6}$ [120]. The expected
theoretical window in $\chi$ mass is $0.94\,\rm GeV<m_{\chi}<4.34\,\rm GeV$
[29]. Studies from Belle [121] and BaBar [122] limit the available parameter
space in the mass region of $1-4.4\,\rm GeV$, and it is anticipated that the
remaining parameter space can be probed at Belle-II [122]. This model is
particularly close to the model we study, in both its visible and hidden-
sector components. In this paper we have established severe limits on the
$\varepsilon_{n\chi}$ and $\varepsilon_{\Lambda\chi}$ mixing parameters for
$\chi$ masses satisfying $m_{\chi}\lesssim 1400\,\rm MeV$, as shown in Fig.
13. In this mass region and for the regions of hidden-sector parameter space
we have chosen, our limits constrain the flavor structure of models of
$B$-mesogenesis, and we now turn to those and their implications.

Different UV completions of $B$-mesogenesis models fare differently in light
of our constraints. Here we consider versions in which only one extra particle
is needed. For example, in Ref. [14], a color-triplet, SU(2)L singlet scalar
with the SM quantum numbers $(3,1,-1/3)$ is used, though a scalar of form
$(3,1,+2/3)$ [29] or a vector of form $(3,2,-1/6)$ [16] are noted
alternatives. We do not consider this list exhaustive. The two scalars are
just the leptoquarks we have noted in Sec. 2: $S_{1}^{\ast}$ and
$\bar{S}_{1}^{*}$ [10, 15]. The phenomenology of these specific models has
been studied, and in order to explain the baryon asymmetry, the dark matter
abundance, and all empirical constraints, including those on $|\Delta F|=2$
meson mixing, a rich flavor pattern of couplings to quarks is needed [29].

To determine the implications of our constraints, we first note the structure
of the Lagrangian for each UV completion, following Ref. [16], though we write
our 2-spinors as in Ref. [123] and employ the conventions given there.
Denoting the new scalars as $Y_{\rm Y}$ and the new vector as $X^{\mu}$, we
have

$\displaystyle{\cal L}_{Y_{\frac{2}{3}}}\supset-
y_{d_{a}d_{b}}\epsilon_{\alpha\beta\gamma}Y^{\alpha}_{\frac{2}{3}}d^{\beta}_{a}d^{\gamma}_{b}-y_{\chi
u_{c}}Y^{\alpha\,*}_{\frac{2}{3}}\chi u^{\alpha}_{c}+\rm h.c.\,,$ (7.1)
$\displaystyle{\cal L}_{Y_{\\!{-\frac{1}{3}}}}\supset-
y_{u_{a}d_{b}}\epsilon_{\alpha\beta\gamma}Y^{\alpha}_{-\frac{1}{3}}u^{\beta}_{a}d^{\gamma}_{b}-y_{\chi
d_{c}}Y^{\alpha\,*}_{-\frac{1}{3}}\chi
d^{\alpha}_{c}-y_{Q_{a}Q_{b}}\epsilon_{\alpha\beta\gamma}Y^{\alpha}_{-\frac{1}{3}}\left(Q^{\beta}_{a}\varepsilon
Q^{\gamma}_{b}\right)+\rm h.c.\,,$ (7.2) $\displaystyle{\cal L}_{X}\supset-
y_{Q_{a}d_{b}}\epsilon_{\alpha\beta\gamma}\left(X^{\alpha}_{\mu}\varepsilon
Q^{\beta}_{a}\right)\sigma^{\mu}d^{\gamma}_{b}-y_{\chi
Q_{c}}\left(X^{\dagger\,\alpha}_{\mu}Q^{\alpha}_{c}\right)\sigma^{\mu}\chi+\rm
h.c.\,,$ (7.3)

where $\varepsilon$ is an antisymmetric tensor in the two-spinor indices and
$\chi$ is a right-handed field. With the ${\cal B}$ assignments of $-2/3$ for
the scalars $Y_{\frac{2}{3}}$ and $Y_{-\frac{1}{3}}$ and ${\cal B}=1$ for
$\chi$, the noted interactions conserve baryon number. In Refs. [29, 15]
$y_{Q_{a}Q_{b}}$ (for each $a,b$) is taken to be zero. The color structure of
the first term of Eq. (7.1) requires that the product of $d$-like quarks be
antisymmetric in the generation indices $a,b$, which follows because we have
assumed the scalar is a color triplet. As for the last case, the vector
$X^{\mu}$ can be written in two-spinor form as [16]

$X^{\mu}=\left(\stackrel{{\scriptstyle{Y^{\mu}_{\frac{2}{3}}}}}{{{}_{Y^{\mu}_{\\!{-\frac{1}{3}}}}}}\right)$
(7.4)

and thus through Eq. (7.3) we see that both scalars couple to left-handed
quarks. We have defined our scalar-fermion couplings in the flavor basis,
rather than the mass basis, but in the case of couplings to right-handed
quarks no distinction needs be made. However, in the case of couplings to
left-handed quarks we need to rotate the fields to the mass basis, to parallel
the treatment of the charged weak current in the SM. As a result, a flavor
diagonal coupling to a left-handed quark of a single flavor can engender a
contribution to a flavor-changing neutral current (FCNC). In the example of
$Z^{\prime}$ models, satisfying FCNC constraints with a large $Z^{\prime}$
coupling requires nearly flavor-universal couplings [124], where we note that
in the flavor universal limit the unitary structure of the CKM matrix makes
the FCNC couplings vanish. We will see that this effect does not appear here
because our scalars do not ever couple to two left-handed quarks of the same
flavor. Replacing a left-handed flavor state $d^{i}$ with a combination of
mass states via $V_{ij}d_{j}$, with $V$ the CKM matrix, we see that the
$X^{\mu}$ completion does lead to a FCNC of form [16]

${\cal L}_{X\,;\,\rm FCNC}\supset-
y_{Q_{a}d_{b}}Y_{\frac{2}{3}\,\mu}V_{aa^{\prime}}\bar{d}^{c}_{a^{\prime}}\gamma^{\mu}P_{R}d_{b}\,,$
(7.5)

where we have employed 4-component notation. This interaction engenders not
only $|\Delta F|=2$ meson-mixing but also structures such as
$B_{(s)}\to{\bar{K}}$ or $B_{(s)}\to\pi^{0}$ at tree level, which can be
probed through $B$ decay studies. We also see explicitly that the structure of
the vertex does not require a flavor universal coupling to control the size of
the effect. Thus there are no particular flavor conspiracies in satisfying the
$|\Delta F|=2$ constraints, and to determine the impact of the constraints we
have found on the mixing parameters $\varepsilon_{n\chi}$ and
$\varepsilon_{\Lambda\gamma}$ on these models, it suffices to consider the
contributions to these quantities from the scalar-fermion couplings with a
particular UV complete model.

Figure 16:  (Color Online) Feynman diagrams for $n\to\chi\gamma$ and
$\Lambda\to\chi\gamma$ decays as mediated by the baryon-number-carrying
scalars $Y_{\frac{2}{3}}$ and $Y_{-\frac{1}{3}}$ as discussed in text, after
Ref. [15].

Considering, then, the flavor structure of Eq. (7.1) we see that
$n\to\chi\gamma$ cannot occur at tree level, and a loop graph with $W$ and
$Y_{\rm Y}$ exchange is needed to generate the process [15]. The opposite
situation is true for $\Lambda\to\chi\gamma$, with Eq. (7.1) and Eq. (7.2)
yielding that process at tree level and one loop level, respectively. The
pertinent Feynman diagrams are illustrated in Fig. 16, replacing the
illustration of Fig.(1). Noting Eqs. (2.5) and (2.7), it is apparent that the
mixing parameters $\varepsilon_{n\chi}$ and $\varepsilon_{\Lambda\chi}$ depend
very differently on the underlying scalar-fermion couplings in the two cases —
we refer to Ref. [15] for explicit expressions. In particular, the one-loop
diagrams bring in a coupling to the $b$ quark as well, with the following
combinations of couplings:

$\displaystyle y_{db}y_{\chi u}\,;\quad y_{sb}y_{\chi u}$ (7.6) $\displaystyle
y_{db}y_{\chi c}\,;\quad y_{sb}y_{\chi c}$ (7.7) $\displaystyle y_{db}y_{\chi
t}\,;\quad y_{sb}y_{\chi t}$ (7.8)

each of which could saturate the bound we have found for
$\varepsilon_{n\chi}$. In regards to the mechanism of $B$-mesogenesis,
operators with the flavor combinations $\chi bud$, $\chi bud$, $\chi bcd$, and
$\chi bcs$ are pertinent, and they take one of three forms [29]

${\cal\theta}^{(1)}_{ij}=(\chi
b)(u_{i}d_{j})\,,\quad{\cal\theta}^{(2)}_{ij}=(\chi
d_{j})(u_{i}b)\,,\quad{\cal\theta}^{(3)}_{ij}=(\chi u_{i})(d_{j}b)\,,$ (7.9)

where $i\in d,s$, $j\in u,c$, and the colors have been contracted to form a
color singlet in each case. Taking the couplings in Eq. (7.8) one at a time,
we find that saturating our $\varepsilon_{n\chi}$ constraint we have found
limits the coefficient of each of the ${\cal\theta}^{(3)}_{ij}$ operators to
be powers of ten smaller than that needed for $B$-mesogenesis to be successful
[29]. We emphasize, however, that this is particular to the mass window in
$\chi$ and region of hidden-sector parameter space we have noted. For the
$Y_{\frac{2}{3}}$ scalar, those are the operators that would act — and thus we
have ruled out this specific model for $B$-mesogenesis under the conditions we
have noted. The other UV completions we have considered are not similarly
constrained, because the $\varepsilon_{\chi\Lambda}$ constraints limit just
the flavor combinations $y_{\chi b}y_{du}$ and $y_{\chi b}y_{dc}$ pertinent to
$B$-mesogenesis — the other flavor combinations associated with
${\cal\theta}^{(1)}_{ij}$ and ${\cal\theta}^{(2)}_{ij}$ remain unconstrained
despite the severity of our limits.

##  8 Summary

BNV has not yet been observed in terrestrial experiments, and its deep ties to
explanations of the observationally well-established cosmic baryon asymmetry
[2] argue persuasively for its investigation on broader fronts. Previously, we
have considered how it might eventually be discovered through precision
measurements of neutron star observables, particularly those of changes in the
binary-pulsar period, familiar from tests of general relativity [5]. Thus far
we have found limits, and they are macroscopic ones, in that they emerge from
the consideration of a neutron star as a whole. Such constraints miss a
concrete connection to particle physics, and it is badly needed: regardless of
whether we continue to constrain or, finally, discern the existence of BNV (in
contradistinction to a failure of general relativity) from these studies,
further theoretical progress on the problem of BNV requires constraints on the
particle physics models of BNV themselves. In this paper, we have developed
just such a connection, using a concrete description of the neutron star
interior based on a relativistic mean-field theory in hadronic degrees of
freedom [17, 18, 19] that successfully confronts existing macroscopic
properties of neutron stars [21]. Within this context, we have developed how
to assess the rates for BNV particle processes in dense matter, and we present
explicit rates for benchmark processes, particularly ${\cal B}\to\chi\gamma$,
considering its rate both at local points within a neutron star as well as its
volume rate after integration over the structure of the entire star, up to its
crust. Although our in-medium formalism is germane to the evaluation of any
particle process in the dense medium of a neutron star, the focus of this
paper — noting current sensitivities — is that of apparent BNV through baryon
decays to hidden-sector particles. Finally, with this in place, we match the
computed rate to our inferred limits on anomalous binary-period lengthening,
i.e., how the binary itself spins down, to set one-sided limits at 2$\sigma$
on the mixing parameters $\varepsilon_{{\cal B}\chi}$, for individual binary-
pulsar systems, as well as a combined limit for all of the studied systems.

As a result of these studies, we discover that neutron stars open new windows