Neutron Stars with Baryon Number Violation, Probing Dark Sectors

The neutron lifetime anomaly has been used to motivate the introduction of new physics with hidden-sector particles coupled to baryon number, and on which neutron stars provide powerful constraints. Although the neutron lifetime anomaly may eventually prove to be of mundane origin, we use it as motivation for a broader review of the ways that baryon number violation, be it real or apparent, and dark sectors can intertwine and how neutron star observables, both present and future, can constrain them.


Introduction
stars are remarkable for their very existence: they are the densest objects known in the observable Universe. With an upper mass limit of about 2 M and a typical radius of about 12 km, their central density can be in excess of a few times nuclear matter saturation density, where n sat ≈ 0.16 fm −3 [1][2][3][4][5][6]. A typical neutron star can contain in excess of 10 57 baryons [2]. Thus a neutron star is an exquisitely sensitive environment in which to study the possibility of new sources of baryon number violation (BNV).
That a neutron star (NS) can serve as a "graveyard" of different theoretical extensions of the Standard Model (SM), particularly in regards to suggested solutions to the dark matter problem, is long known [7,8]. In these scenarios, the capture of the suggested dark matter candidate by a NS alters the latter so severely that the existence of the model tested is precluded by that of the NS. The environment of a proto-NS is also long-known to be a sensitive discriminant of light new physics, such as axions [9][10][11][12][13][14][15] or dark photons [16][17][18], through its impact on the observed cooling of the star. Yet the advent of gravitational wave (GW) observations of compact object mergers [19], and other observational facilities for the realization of multi-messenger probes of these objects, offer increasingly sensitive probes of X X X X X X X Figure 1: In this schematic, we update earlier notions regarding the interior structure of a neutron star [2] to include not only the possibility of hyperonic matter in the inner regions of the star, but also the possibility of either a quark or a mixed quark/hadron phase, as well as the possibility of new matter or force mediators, "X", that impact its structure, terrestrial constraints on neutron decay to dark final states (Sec. 2). Baryonnumber violation impacts the thermodynamics of the star (Sec. 3) and, be it apparent (Sec. 4), explicit (Sec. 5), or spontaneous (Sec. 6), acts as a source of X. For X of still lighter mass, the cooling of the neutron star, as well as its merger dynamics, can be modified. The nature of the innermost region of the neutron star is unknown, but a dark-matter core or a dark-mediator condensate figure among the possibilities (Sec. 7). Neutron stars typically have masses of O(1 M ) and radii of O(10 km).
latter are also subject to observational studies. Finally, we turn to a brief assessment of the broader ways in which dark sectors can impact neutron star observables before offering our final summary.

Setting the Stage -The Neutron Lifetime Anomaly
Measuring the neutron lifetime with ever-increasing precision has been the ongoing work of decades [20]. This work has been motivated, in part, by anticipating the needs of precision cosmology: from the recognition that the neutron lifetime not only fixes the effective weakinteraction rate in standard Big-Bang nucleosynthesis [36,37], but it also represents the largest experimental uncertainty in the prediction of the 4 He yield [38]. Over the last decade, or so, a significant disagreement has appeared in the determination of the neutron lifetime via two distinct methods [20]. This is the neutron lifetime anomaly -namely, that the neutron lifetime determined from the detection of its decay products [39][40][41], as studied in neutron beam experiments, is longer than that inferred from counting the surviving neutrons [42][43][44][45][46][47][48][49][50][51], as studied in neutron bottle or trap experiments. That the lifetime inferred from "counting the living" is smaller than that from "counting the dead" is evocative [20] -perhaps the neutron decays to exotic, weakly-coupled final states and that could explain the difference. To our knowledge, the first work along these lines is due to Berezhiani [52]. 3 Fornal and Grinstein have developed new-physics models particular to the anomaly, yielding exotic final states in which no proton appears but containing particles that carry baryon number; and they have noted that these models can be probed through ancillary empirical tests [21], which have been made [53][54][55], with null results thus far. It is quite possible that the anomaly could yet be explained through a subtle combination of experimental systematic effects in either or both types of experiments. In addition, different lines of evidence suggest that the entirety of the anomaly would not reasonably arise from new-physics effects. Powerful constraints come from the existence of neutron stars [22][23][24], as well as from the connection to precision measurements of β-decay correlations in the Standard Model [35]. Yet the possibility of new-physics effects remain, and we use this as a springboard to consider the interconnections between neutron-star physics with new-physics models that contain baryon number, and its violation, and dark sectors in a broad way.
The perspective on the interconnections between the neutron lifetime, β-decay observables, and SM radiative corrections (RCs) has shifted, due to new theoretical and experimental results, and we pause to consider these developments before proceeding to the main body of our article.

Constraints from empirical studies of neutron β decay within the SM
Precision measurements of β-decay observables, along with accurate calculations of electroweak RCs, yield precision tests of the SM [56][57][58][59][60][61][62]. For example, the unitarity test stemming from the first row of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, namely, is the most precise known [30]. A significantly nonzero value of Σ CKM would establish the existence of physics beyond the SM, and there has been much discussion of the accuracy of the ingredients needed to determine V ud and V us , since |V ub | 2 ∼ O(10 −5 ) and is thus relatively negligible. An important ingredient is the electroweak RCs, which includes the evaluation of the γW box diagram and in which non-perturbative effects appear to be significantly larger [63] than earlier estimated [62]. An updated analysis [64] of the latter agrees with the sense of the shift but finds a result intermediate to that of Refs. [62,63]. Moreover, both electromagnetic and isospin corrections are also key to determining consistent values of |V us | from K 2 and K 3 decays [65,66], with the current scale factor of S = 2.7 [30] seemingly indicating the need for further theoretical and experimental work [67][68][69][70]. Here we focus on the V − A structure of neutron decay in the SM, and the concomitant ties between its 3 See references therein.
observables, as this constrains the possibility that the neutron lifetime anomaly comes from new physics [35] -and we refer to Refs. [71][72][73][74] for reviews. To that end, we consider [35] 1 where τ n is the neutron lifetime, G F is the Fermi constant determined from µ decay after QED radiative corrections are subtracted, G F = 1.1663787(6) × 10 −5 GeV −2 [75], m e is the electron mass, m e = 0.5109989461 (31) MeV [30], g A is the axial vector coupling constant of the nucleon, 4 δ RC is the electroweak radiative correction, and f is the statistical rate function [76]. The last follows from the allowed phase space, the recoil corrections assessed in the isospin symmetric limit, and the Coulomb correction in the e−p final state as encoded in the Fermi function -and it has been reevaluated to yield f = 1.6887(1) [61]. In the SM, the combination is tightly constrained, because η is 4908.6(1.9) s [62], or 4903.6(1.0) s [63], 4905.7(1.5) s [64], depending on the calculation of δ RC used. The uncertainties are dominated by that in δ RC ; thus we use that to determine the reported errors in η. Using Eq. (3), different measurements of the neutron lifetime at 68% confidence level (C.L.) give the diagonal bands in the |V ud | versus g A plot shown in Fig. 2. We combine statistical and systematic errors in quadrature, assuming uncorrelated errors, to realize the bands shown in Fig. 2.
We now turn to the inputs used to generate Fig. 2. We consider the most precise beam lifetime result, 887.7 ± 1.2 (stat) ± 1.9 (sys) s [41], the most recent PDG average of the neutron lifetime determined from bottle and trap experiments, 879.4 ± 0.6 s [30], and the most precisely measured neutron lifetime, determined in a magnetic trap experiment, 877.75 ± 0.28 (stat) +0. 22 −0. 16 (sys) s [51]. The numerical difference between the beam and bottle/trap measurements constitutes the neutron lifetime anomaly. The latest magnetic trap result [51] is compatible with an earlier measurement with the same method [50] but is not included in the PDG average. We pull it out for explicit study because of the possibility of significantly large and/or underestimated systematic errors in the older experiments; e.g., the scale factor in the average reported by the PDG is 1.6 [30].
We also include the value of g A from both measurements of neutron decay correlation coefficients and computations within lattice QCD, as well as the value of |V ud | from superallowed 0 + → 0 + nuclear decays, in which the effect of the axial vector current enters in RCs, which can be modified by nuclear structure as well [81]. The empirical determination of g A comes from that of λ ≡ |g A /g V |, where the interpretation of the measured A and a correlation coefficients in terms of λ requires the application of radiative and recoil corrections, and in the latter additional hadronic matrix elements appear [82,83]. The weak magnetism  Figure 2: The SM relationship between the CKM element |V ud | and the axial coupling constant g A , with different values of the neutron lifetime and estimated RCs, with values taken at 68% C.L. throughout. We consider the beam lifetime [41] (cyan), as well as the 2020 PDG average of bottle/trap lifetimes [30] (purple) and the latest magnetic trap result [77] (green), with the RCs of Refs. [62][63][64], respectively, applied in each case, realized from the top to the bottom throughout. Assessments of g A from neutron β-decay from both decay-correlation measurements and lattice QCD calculations are also shown, as is the value of |V ud | from 0 + → 0 + nuclear β decays, which is also sensitive to the precise value of the RCs. The lattice values of g A are the 2021 FLAG average for N f = 2 + 1 + 1 flavors [78] and the 2018 CalLat result [79]. The decay correlation determinations of |λ| are from the 2020 PDG compilation [30] and from Ref. [80] (ochre). We refer to the text for all details.
contribution therein, in the isospin limit, is fixed by the determination of the isovector nucleon magnetic moment [30], and the matrix elements that are nonzero only if the u and d quarks differ in mass [84] are set to zero. We note that the PDG average is λ = 1.2754 (13) with a scale factor of 2.7, whereas the most precise determination, from that of the A decay correlation in n decay, is λ = 1.27641(56) [80]. This last result is consistent with the other two most precise determinations of λ, which are also determined from A [85,86]. In what follows, following common practice, we employ λ as g A , which is strongly supported by an analysis of its RCs [87]. As for the lattice QCD results, we note the recent FLAG average from simulations with N f = 2 + 1 + 1 flavors, g A = 1.246 (28), as well as the most precise lattice result, g A = 1.271 (13) [79]. Both calculations have significant errors, but the latter result is compatible with the most precise empirical determination of λ. We note that the g 2 form factor, which vanishes in the isospin limit, and/or the possibility of scalar and/or tensor currents can make the two assessments differ [73,88,89]. We refer to Refs. [73,88] for complete expressions for the hadronic vector and axial-vector currents. These contri-leptons and photons and the hyperons do not. The conservation of baryon number and electric charge are used as constraints in determining the state of matter inside neutron stars. This is achieved by finding the ground state of electrically neutral (Q = 0), cold (T = 0) matter at a given fixed baryon number density 5 (n B ), i.e., by defining a function and minimizing it with respect to the individual number densities (n i ), i.e., ∂Φ/∂n i = 0. The sum in Eq. (4) is over all of the particles present in the matter, α and β are the Lagrange multipliers enforcing the electric charge and baryon number conservation respectively, ε i is the energy density that depends on the density of each species j, Q i is the electric charge, and B i is the baryon number of a particle of type i. The two constraints (α and β) relate the chemical potentials of particles present in a neutron star at chemical equilibrium. The chemical potential for a particle of type i (if present and in equilibrium with the matter inside the star) is given by µ i = B i µ n − Q i µ e , in which µ e and µ n are the chemical potentials for the electron and neutron, respectively. In general a violation of baryon number conservation would change the chemical equilibrium and the composition of the star in a model-dependent manner. However, as we show later (3.1), a model-independent probe is feasible for a class of sufficiently slow BNV processes (τ BNV τ weak , τ hyd ) which act as small baryon perturbations over time. In response to each of these out-of-equilibrium perturbations, the star regains its chemical equilibrium using the standard (faster) baryon-number conserving (BNC) reactions (e.g., weak interactions) and ends up with a lower total baryon number. Here τ weak is the timescale for weak interactions in the neutron star medium (e.g., Urca reactions) which may be different from the free neutron lifetime (τ n ), and τ hyd is the time needed to adjust to and maintain hydrostatic equilibrium. We will explain these timescales further in Sec. 3.1. We will study the generic effects of BNV in this section, and defer a discussion of specific models to the following sections.
We assume that the neutron star matter (without BNV) has spherical symmetry. 6 This spherical symmetry would remain intact after the inclusion of BNV processes, because BNV processes would be sourced by the matter already present in the star. Furthermore, we work in a quasi-static regime in which the changes to the metric (g µν ) are very slow in time. This warrants the use of the line element for a static spherically symmetric system [97] in which ν(r), λ(r) are solutions to the Einstein field equations [98], G µν = −8πGT µν , in which G µν is Einstein's tensor, G is the gravitational constant, and T µν is the stress-energy tensor. We use a geometric unit system in which the speed of light (c) and G are both set to unity, i.e., G = c = 1. For a perfect fluid T µν has the form in which p and ε are the local pressure and energy density of the fluid respectively, and u µ is the 4-velocity of the fluid, which has zero 3-velocity (u i = 0, u 0 = 0) in a static star. The time component of u µ is calculated (from the normalization condition: g µν u ν u µ = 1) to be Therefore, for a static perfect fluid the only non-zero components of the stress-energy tensor are given by Moreover, λ(r) is found to be g 11 (r) = − exp(2λ(r)) = −(1 − 2M (r)/r) −1 , with M (r) being the total mass included within radius r: Further simplification of the Einstein field equations yields a differential equation for the pressure inside the star: Eq. (9) together with Eq. (10) are known as the Tolman-Oppenheimer-Volkoff (TOV) [99,100] equations. The pressure (p) and energy density (ε) are in general functions of the number density of baryons in the rest-frame of the fluid (n) and temperature (T ). Neutron stars cool down to T E Fermi ≈ 30 MeV within a few seconds after formation. Therefore, the thermal contribution to the pressure and energy density can be neglected, i.e., p(n, T ) = p(n, T = 0) = p(n) and ε(n, T ) = ε(n, T = 0) = ε(n). We can then deduce both ε and p from the knowledge of n.
The TOV equations (9) and (10) can be integrated with the initial conditions M (0) = 0 and ε(0) = ε c up to p(r) = 0 (surface of the star). Therefore, for any given equation of state (EoS), there is a unique family of stars parameterized by the central energy density (ε c ) also known as the single parameter sequence [101] of stars. We note in passing that in the case of a rotating neutron star [102], or a neutron star with a dark matter core [103,104], extra parameters in addition to ε c are needed to describe the star uniquely. We discuss the possibility of generalizing our analysis to the rotating case in Sec. 3.2.
The baryon number current is given in terms of the fluid velocity (u µ ) and the baryon number density (n) by [105] j µ (r) = n(r)u µ . Bearing in mind that the invariant 4-volume is given by √ −g d 4 x (g ≡ det|g µν |), the total baryon number in a static, spherically symmetric neutron star is given by [101] in which we used √ −g = exp(ν(r) + λ(r)) r 2 sin θ.

General conditions
The exact consequences of BNV processes for neutron stars depend on the modeling of neutron stars' structure (in the absence of BNV) and the particle physics model producing those specific BNV reactions. Although the details of neutron star models may change the numerical results in this section, we expect that their order of magnitude and qualitative behaviour remain intact. On the other hand, the particle physics modeling of BNV could have drastic effects, and a generic study would require imposing some simplifying assumptions on the BNV models. We attempt to find a minimal set of conditions that makes such a broad investigation viable. There are two major effects that the inclusion of a BNV process in the star can generate. The first one is caused by the relaxation of the baryon number conservation constraint (i.e., β = 0 in Eq. (4)). The system is now allowed to transition into more energetically favorable states subject to electric charge conservation (only). The EoS would be different from the standard BNC EoS, and this could cause drastic changes to the composition of the matter inside stars. As we will explain below, we are interested in slow BNV processes for which this effect is eliminated and the BNC EoS is revived. The second effect is because of the production of new particles in the BNV process that are not otherwise present in the star. The modification in the pressure and energy density of the matter (i.e., the EoS) would depend on the specific final states produced in the BNV process. For example, a fermionic final state would exert a Fermi pressure and its production would be Pauli suppressed, whereas a scalar final state (with negligible self-interaction) would significantly reduce the pressure of the system. This is an obstacle to our model-independent analysis objective. Accordingly, we set forth the following essential condition for the BNV processes that we consider in this section: The final states are either already present (via BNC processes) or if they are not already present, then they maintain a negligible contribution to the EoS.
In cases with new particles in the final states, the above condition can be realized if: 1. The new final state particles participate in annihilation or decay channels to yield particles already present plus neutrinos and photons. 2. Their production rate (Γ BNV ) is much less than their elimination rate via annihilation Γ ann or decay Γ dec .
We extend the above constraints by also demanding that: The BNV rate(s) are slower than the weak-interaction processes that they activate in the neutron star in response to their presence, i.e., Γ BNV < Γ weak .
We should elaborate on the nature of these responses and their timescales which vary greatly depending on the EoS, mass, and temperature of the neutron star. Reactions are generally suppressed by the small phase space available to fermions in chemical equilibrium in a cold, degenerate state. This is because the fraction of fermions on the edge of their Fermi surface that can undergo inelastic scattering is about ∼ k B T /E F 1. The slow BNV processes perturb the chemical equilibrium in the star. This imbalance in chemical equilibrium temporarily activates or enhances BNC reactions until a new chemical equilibrium is achieved. The exact timescale of the response to BNV reactions by the weak processes (τ weak ) would depend on the specific reaction, and the temperature of the star. For a general estimate of the timescales involved, let us consider the Urca processes [106] which are extensively studied in the context of neutron star cooling theories [107][108][109]. Direct Urca processes involve baryon -decays and lepton ( = e, µ) capture: in which B 1,2 denotes nucleons or hyperons. The nucleonic direct Urca reactions (B 1,2 = n, p) would be active if the proton fraction is above a minimum threshold (n p n n /8) [110,111], which is possible at the inner-core of a heavy neutron star with supranuclear densities (ε 2ε nucl ). The hyperonic direct Urca processes (B 1,2 = Λ, Σ − , . . .) can occur [112] when the neutron chemical potential (µ n ) surpasses the energy of the lowest state of a Λ, and if µ n + µ e is greater than the energy of the lowest state of a Σ − in the neutron star. At lower densities (e.g., in the outer core), the direct Urca processes may be suppressed, but the following modified Urca reactions would still occur [113,114]: Since these modified Urca reactions have two extra degenerate fermions, their rate is suppressed by a factor of (k B T /E F ) 2 compared to the direct processes. In the presence of so-called β-disequilibrium and processes involving the neutron, the sign of the disequilbrium parameter δµ ≡ µ n − µ p − µ e determines which of the reactions in Eq. (12) dominates, so that if δµ > 0, β decay occurs and if δµ < 0, electron capture occurs. Thus for δµ > 0, which yields netν e production, the rates for direct and modified Urca processes in a simple npe model (i.e., a degenerate Fermi gas consisting of neutrons, protons and electrons) are given by [115] Γ Urca = 8.86 × 10 31 n e n sat 1/3 Γ mod Urca = 5.91 × 10 23 n e n sat 1/3 in which n sat = 0.16 fm −3 , T 9 ≡ T /(10 9 K), and the dimensionless functions G d and G m are defined as For δµ < 0, yielding ν e production, the rates evaluate to the same numerical value as in the δµ > 0 case, implying the replacement δµ → −δµ. The various Urca timescales are plotted in Fig. 3 as a function of |δµ|/T 9 . We take the baryon number density to be n ≈ 0.5 fm −3 , and (n e /n) = (n p /n) ≈ 11.2%, n n /n ≈ 88.8%. We can see that the Urca rates are highly sensitive to δµ. We expect BNV reactions to generate β-disequilibrium of order |δµ| ≈ 10 − 100 MeV from kinematics. The other important timescale is the hydrodynamical relaxation time of the star (τ hyd ) for regaining hydrostatic equilibrium. We can approximate this timescale by finding the period of small oscillations of a uniform Newtonian fluid in hydrostatic equilibrium: in which ε(r) = ε, m(r) = (4/3)πr 3 ε. The period of small adiabatic oscillations is given by , in which Γ 1 ≡ (n/p) (∂p/∂n) s is the mean adiabatic index. The hydrodynamical timescale (τ hyd ≡ 2π/ω) is then given by in which G is the Newton's constant, Γ 1 ∈ [2, 3] [116] and we used ε ∈ [10 14 , 10 15 ] (g/cm 3 ). The second condition on BNV (i.e., Γ BNV < Γ weak ), permits the use of the standard BNC EoS. In essence, because baryon number is conserved on shorter timescales (compared to τ BNV ), the instantaneous states of neutron stars are governed by the same BNC EoS before the inclusion of the BNV reactions. Nevertheless, the star will be slowly leaking baryon number and its structure will change over time.
To better clarify these statements, let us consider a specific example: a process with |∆B| = 2 such as nn → e − e + . In this case the positrons would start annihilating with the electrons that are already present at a rate Γ ann ≈ n e − σ e + e − c, in which n e − and σ e + e − are the electron number density, and electron-positron annihilation cross section to two photons and c = 1 is the speed of light. In general, if the particles produced in the final states have a cross section σ ≈ 10 −43 cm 2 , and their counterpart in the annihilation has a density of at least 10 −3 fm −3 , then we get Γ ann ≈ 10 3 s −1 . As long as the rate for the BNV processes under consideration is smaller than this annihilation (or decay) rate Γ BNV Γ ann , our assumptions in this section are valid. In other words, we are considering the cases in which, through a chain of reactions, the final products end up being those particles that are already present in the star plus photons and neutrinos. Neutron stars become transparent to low-energy (E ν MeV) neutrinos as they cool down to temperatures T MeV. The mean-free path of electron-neutrinos (ν e ) with E ν E F (e) is given by [114] λ νe ≈ 5 × 10 6 km ε nucl ε in which ε nucl = 3.7 × 10 14 g/cm 3 . This mean-free path is much larger than the typical neutron star radius (∼ 10 km) for E ν 5 MeV. The mean-free path of other neutrino flavors would be even larger in the SM. In comparison, photons have a much shorter meanfree path, and would deposit most of their energies before they can escape, which would result in heating of the neutron star.

Effects of Slow BNV perturbations
If the BNV processes are much slower than the weak interaction rates, then their sole effect is to change the baryon number density inside the star (n) and perturb the system out of equilibrium (ñ). 7 The system will respond by adjusting the densities of each species via reactions that conserve baryon number and electric charge. Therefore, the final equilibrium state (n ) of the star would be the same as a star with a lower baryon number (B < B) or equivalently a lower central energy density (ε c < ε c ) from the same single parameter sequence (Table 1). The total baryon number after the perturbation generated by the BNV process (B ) is (from Eq. (11)) 7 The perturbations are Eulerian, i.e., they are changes measured by an observer at a fixed point (t, r, θ, ϕ).

13
in which R and M are the equilibrium radius and mass of the star after the perturbation.
Since the structure of a spherically symmetric neutron star is fully determined (via Eq. (9) and Eq. (10)) by an EoS and the value of the energy density at the origin (ε c ), changes in the total baryon number δB = B − B can be uniquely mapped onto δε c . Therefore, we can quantify the changes in a neutron star observable O along the single parameter sequence as If we include the effects of rotations, then two parameters are needed: (ε c , ω c ), in which ω c ≡ Ω − ω(r), Ω is the angular frequency of the star, and ω(r) is the frequency of the local inertial frame [102]. Since BNV does not change the angular momentum of the star, L(ε c , ω c ) = L (ε c , ω c ) and B (ε c , ω c ) = B(ε c , ω c ) + δB can still be solved for a unique set (δε c , δω c ). Therefore, the generalization to the rotating case is possible but we will consider the static non-rotating scenario for simplicity. Given the relative rate of change in the baryon numberḂ/B ≡ (dB/dt)/B, we find the relative rate of change in the observable O to bė The most stringent bound on BNV rates can be set with a careful choice of O. The relative change in O is given by δO/O ∝ d ln O/dε c , such that quantities that are more sensitive to ε c will have a greater variation. At the same time, an observable quantity (O) that can be measured to a higher precision would yield a better constraint. Thus far we have presented a global examination of a neutron star's macroscopic observable quantities in terms of a given rate of change in the total baryon number (Ḃ). Alternatively, we can inspect the local properties of Eq. (22) in terms of the baryon number density (n(r)). Since neutron stars are compact objects with n(R) n(r < R), the baryon conservation condition (Eq. (22)) can be expanded to leading order in the perturbation as in which δM (r) is the variation of the mass function (Eq. (9)) and it is given by with δε(r) = (dε/dn)δn(r). In principle the integral equation in (25) should be solved for δn(r) givenδn(r) and in a consistent manner with the TOV equations (Eq. (10)). The functionδn(r) can be calculated from the rates of the BNV processes, which implicitly depend on time (t) and radial position (r). Specifically we havẽ in which f (n) is the relative fraction of the decaying species (i) in the BNV process, i.e., n i (r, t) = f (n) × n(r, t), and ∆B is the change in the baryon number per each decay. We expect a BNV process to slow down and eventually halt once the density falls below the threshold needed for that specific BNV reaction. For example, if only hyperons are involved in the BNV process, then as the baryon number density drops below a certain threshold the star will essentially be depleted of hyperons and the BNV reaction stops. If we assume that the conditions from Sec. 3.1 are satisfied by the BNV processes, then n(r, t) is simply given by the solutions to the one-parameter sequence. Specifically, let us assume that BNV is only active in regions of the star that have a baryon number density greater than a certain threshold, i.e., n BNV > N n sat = 0.16N fm −3 , in which N is a ratio of order unity. In other words, we assume the form Γ BNV (n) = 0 n < N n sat Γ BNV n ≥ N n sat Starting from any point (ε i c ) along the sequence in ε c with a known density (n(r)) and mass profile (M (r)), we can find the change in the total baryon number of the star δB using Eq. (27) and We can then find the δε c corresponding to this δB, and repeat this process for the new point (ε f c = ε i c + δε c ) until the central density n(r = 0) falls below the threshold (n BNV ), and the BNV process is deactivated. We have adopted the hyperonic EoS H3 (K = 300 MeV, m /m = 0.70, x σ = 0.60) from Ref. [117] and plotted the results for f = ∆B = 1 and Γ BNV = 10 −10 yr −1 in Fig. 4. We considered a set of four BNV threshold densities n BNV > {1, 2, 3, 5}n sat for illustration. As expected, BNV processes with a lower density threshold have a much more significant effect because the BNV is active in a larger region inside the star. We expect the BNV effects to be most notable for a timescale T ∼ 1/Γ BNV , and we see that the epoch between 10 8 (yr) T 10 10 (yr) shows the fastest evolution for the neutron star. Similar results can be found for other choices of parameters f , ∆B, and Γ BNV via proper scaling of the timescale T .

Constraining BNV from neutron star observations:
We calculated the expression in Eq. (24) for a sequence of neutron stars with the hyperonic EoS H3 from Ref. [117].  Fig. 5. Note that the sequence becomes unstable (dotted curves) where dM/dε = 0 [119]. The rates of change for mass, radius, and moment of inertia from the expression in Eq. (24) are plotted on the right side of Fig. 5. We can see that for almost all the points along the sequence the ratio in Eq. (24) is of O(1). This indicates that relative changes in a neutron star's observable would be at the same order as the relative changes in the baryon number. The exceptions to this statement . We have used the hyperonic EoS H3 from [117] together with the BPS EoS [118] for lower densities.
are in the beginning and at the end of the sequence where B has extrema, i.e., dB/dε c = 0. Although the exact locations of these extrema depend on the EoS, their existence is independent of it. Close to these points an infinitesimal change in the baryon number would result in a substantial variation in other observable quantities. As a consequence, the existence of BNV processes as we have outlined would shift the heavier neutron stars away from the maximum mass, and make light neutron stars close to the minimum mass unstable. The main question regarding the destabilizing effects of BNV on low-mass neutron stars is how light neutron stars in nature can be. In order to answer this question, it is helpful to consider different aspects of neutron stars including their theoretical mass limit, observation, birth, and evolution. The absolute (theoretical) minimum mass for a stable cold neutron star is about M min = 0.1 M [118]. If the mass of a neutron star in hydrostatic equilibrium decreases below M min , then it will become unstable and explode [120]. This will result in a burst of hard X-rays and soft gamma rays with a total energy of 10 43 −10 47 erg [120]. On the observation side, the current minimum neutron star mass is observed to be in the range 1 -1.1 M [4,121], albeit with significant errors, with one of the lightest neutron stars observed, PSR J0453+1559, having a more precisely determined mass of M = 1.174(4) M [121]. As for their birth, neutron stars can either be directly born in explosive death of a massive star (> 8M ) [122,123], or as a result of white dwarf accretion-induced collapses [124]. Typical core-collapse supernova creation scenarios for neutron stars predict that their masses must be at least 1 M [3,125]. Furthermore, the minimum mass of hot proto-neutron stars We have chosen the hyperonic EoS H3 from [117] together with the BPS EoS [118] for lower densities.
predicted by the models considered in Ref. [126] and formed from supernovae is in the range 0.89 -1.13 M . Therefore, the core-collapse supernova paradigm would appear to impose a lower limit on neutron star masses at birth [4,127]. As an example of a non-standard scenario, the effects of dark matter accretion by a massive white dwarf and the subsequent core collapse into a light proto-neutron star is studied in Ref. [128]. They show that a dark matter admixed core with 0.01 M of dark matter can lower the minimum mass to 1 M . On a distinct note, the discovery of ZTF J190132.9+145808.7 [129], which is a nearby (41 pc) white dwarf with a mass ∼ 1.3M , opens up the possibility of observing a new formation channel [130] for neutron stars, which may yield relatively lighter masses. Therefore, in the absence of a mass-loss mechanism, it is unlikely to have isolated neutron stars with masses below 1 M i.e., they evolve from the right side to the left. The arrows indicate if an observable is increasing (green) or decreasing (blue). Typical core-collapse supernova scenarios generate neutron stars with masses above ≈ M . Those with masses below M min ≈ 0.1 M and above M max ≈ 1.78 M are unstable. We note that the value of the latter depends on the adopted EoS.
its radius (making it even more susceptible to mass loss in the binary), this accretion can be self-accelerating and lead to the explosion of the low-mass component [131,132]. The existence of BNV processes would only increase this acceleration, as in Fig. 6. If BNV is active in light isolated neutron stars, then it can lower their masses below 1 M , and even possibly result in a mass as small as 0.1 M , in which case it would lead to an explosion which we may be able to detect. However, we should note that if a baryon number density of at least n sat is required for the BNV reactions to occur, it seems implausible for the mass to decrease below ∼ 0.4 -0.5 M (see Fig. 4). In any case, the observation of light isolated neutron stars with M < 1 M would point to a mass-loss mechanism which could be caused by BNV, or else a modification of neutron star genesis theories would be required. Moreover, the continued lack of such observations can be used to constrain BNV processes in neutron stars.
As pulsars emit radiation they lose rotational energy and their spin periods (P s = 2π/Ω) slow down over time [133]. An accurate knowledge of this mechanism along with pulsar timing data can then be used to constrain additional non-standard contributions to this slow-down. A spherically symmetric BNV sink within the neutron star would not affect its angular momentum (L = IΩ), but as we showed in Fig. 5, the moment of inertia (I) changes and as a result the spin period would also change. For heavier pulsars (M 1.7 M ) BNV will have a spin-down effect, whereas for lighter pulsars it will cause a spin-up effect, as indicated in Fig. 6. Assuming that the BNV processes (only) are responsible for a change in a pulsar's spin period, we haveİ in which forḂ < 0, the plus (minus) sign corresponds to spin-up (down), and the last equality comes from |(İ/I)/(Ḃ/B)| ∼ O(1). The spin period and its first derivative have been measured to a remarkable precision for many pulsars. Taking advantage of these precise measurements would require a robust understanding of the pulsar spin-down due to magnetic dipole radiation. An independent measurement of the magnetic field would make it possible to separate the exotic (e.g., BNV) from the standard electromagnetic contributions to the spin-down rate. This is not currently possible and so the extreme precision inṖ s /P s can not be utilized. However, the total value can still be used for comparison, i.e., Alternatively, stronger limits (compared to the pulsar spin-down bounds) can be inferred from the decay rate of a binary pulsar's orbital period (Ṗ b ), which can be used to constrain changes in its components' parameters (e.g., in their masses [134]). The observed relative rate of orbital period decay comprises of various intrinsic and extrinsic terms with the following dominant contributions [135]: The first term is due to gravitational radiation [136], the second term is due to mass-energy loss, and the third term includes the extrinsic effects such as Doppler effects caused by the relative acceleration (due to the Galactic gravitational potential) of a binary pulsar with respect to the solar system. The rate of change in P b due to gravitational radiation (Ṗ GR b ) can be expanded as a series in powers of (v/c) 2 in the post-Newtonian (PN) approximation. To leading order (2.5PN), 8 P GR b is given by the quadrupole formula in Ref. [136], and the next higher correction (3.5PN) is calculated in Ref. [137] which would be needed for more accurate values ofṖ obs b (e.g., in the case of J0737−3039A/B [138]). We set limits onṖĖ b by subtracting this GR contribution (Ṗ GR b ) from the intrinsic orbital-period decay rate,Ṗ int b ≡Ṗ obs b −Ṗ ext b , for three binary pulsar examples ( Table 2): 1. PSR B1913+16: This binary system (Hulse-Taylor binary) is the first binary pulsar ever discovered [139], and consists of a neutron star (M c = 1.39 M ) and a pulsar (M p = 1.44 M ) with a pulse period of 59 ms. We use the results from the analysis in Ref. [140] which is based on timing measurements performed over the last 35 years. 2. PSR J0737−3039A/B: The only known double pulsar was discovered in 2003 [141], and is comprised of two radio pulsars (A and B) with masses M A = 1.34 M , M B = 1.25 M , and with pulse periods of 22.7 ms and 2.8 ms, respectively. We use Ref. [138] which is based on data acquired over 16 years of observation. As a result of the increased accuracy in measurements, the higher-order GR corrections (3.5PN) toṖ GR b , and contribution toṖĖ b from the spin-down of pulsar A are added [138]. 3. PSR J1713+0747: This binary system was discovered in 1993 [142], and it contains a 4.6 ms radio pulsar with M = 1.3 M and a companion white dwarf with M = Name J0737−3039A/B B1913+16 J1713+0747 Table 2: The values of the orbital period (P b ), its intrinsic decay rate (Ṗ int b ), and the gravitational wave radiation contributions (Ṗ GR b ) to it for J0737-3039A/B [138], B1913+16 (Hulse-Taylor) [140], and J1713+0747 [143] binary pulsars. The bound on ( 0.29 M [143]. A sub-microsecond precision is achieved at measuring its pulse time of arrivals [143] owing to the short spin period and its narrow profile. It has a much longer orbital period (P b = 67.8 day) compared to the other two binaries considered here. This is why the inferred limit on (Ṗ b /P b )Ė is one order of magnitude better than the limit from the Hulse-Taylor binary despite the higher precision in the latter.
The 2σ (98% C.L.) bound on (Ṗ b /P b )Ė for each of these binaries is listed in Table 2. We now elaborate on the BNV contributions, via changes in M and I, to (Ṗ b /P b )Ė.
We begin by noting that since the mass loss due to BNV is spherically symmetric, and it appears in the form of photons and neutrinos, implying that we have the very high velocity ejecta needed, we can thus apply the Jean's mode of mass ejection [144]. In this mode the relative rate of change in the binary period is given by [145,146] in which M 1,2 are the masses for each of the components in the binary system, andṀ eff 1,2 is their respective mass loss which, by virtue of Einstein's mass-energy equivalence, can be written asṀ in which we suppressed indices 1, 2. The first term is due to a direct (rest) mass loss, which could be caused by BNV. The second term is due to a change in the moment of inertia (I) which has a direct contribution from BNV (İ BNV ) and an indirect contribution as a result of changes in the angular velocity (Ω). We assume that the latter effect, i.e., I Ω = (dI/dΩ)Ω, is negligible. The third term is the energy loss due to the pulsar spindown which arises from both BNV and electromagnetic radiation. Therefore, after defining for an observable (O), we can rewrite Eq. (32) in terms of the observed pulsar spin periods (P s ), and its observed rate of change (Ṗ s ) aṡ in which the values for η (M,I) can be read from Fig. 5. Equation (31) can then be written in terms of the two separate contributions from BNV and spin-down effects (Ω) as with each of the contributions given by We now turn to the assessment of the period decay rate due to pulsar spin-downṖΩ b for the different binary systems of interest. As we will see, this is only relevant for PSR J0737−3039A/B. In that case, the rates of orbital decay for pulsars A and B are given by 2.3 × 10 −17 I 45 A and 6.3 × 10 −21 I 45 B respectively [147], in which I 45 ≡ I/(10 45 g cm 2 ), [138], and the negligible contribution from pulsar B is ignored, giving the result reported in Table 2. In the case of the Hulse-Taylor binary we use the estimated value for the pulsar (1), and the 68% C.L. bound on the companion neutron star (2) from Ref. [135]: We see that the combined effect is two orders of magnitude smaller than the limit on (Ṗ b /P b )Ė shown in Table 2 and can be safely ignored. In the case of PSR J1713+0747, we estimate the spin-down effects using Eq. (36) with I 45 ≈ 1 [148],Ṗ s = 8.96(3)×10 −21 [143] and ignore the companion white dwarf's contribution as it is much lighter in mass. We see that the spin-down effect is about two orders of magnitude smaller than the limit on (Ṗ b /P b )Ė for this binary and thus it can be neglected.
With these estimates, we subtract the spin-down contributions (Ω) from the energyloss term in Eq. (34) to find 98% C.L. limits on the BNV contributions and record them in Table 2. In principle, given the values for all of the parameters in Eq. (35), we would be able to infer limits on a linear combination of (Ḃ/B) 1,2 , but not on each of them individually. One may attempt to resolve this degeneracy by measuring parameters other than binary period decay rate, such as the individual pulsar spin-down rate. Unfortunately, as we have already mentioned, our understanding of the electromagnetic contributions to the pulsar spin-down rate are not as precise as the GR contributions to the binary period decay rate, making the overall separation of the contributions challenging.
For a general estimate, however, we can make theoretical assumptions about the nature of BNV processes to resolve this issue. In the case of PSR J1713+0747, we assume that BNV would be only active in the pulsar and not the white dwarf companion. For the other two binaries, given that the masses of binary components are close to each other, we assume thatḂ 1 /B 1 ≈Ḃ 2 /B 2 ≡Ḃ/B. The exact value of η coefficients depends on the adopted EoS, but since they are of order unity (Fig. 5), we can approximate them as η (M ) ≈ η (I) ≈ 1. Equation (35) is then simplified as The second term in both cases is more than two orders of magnitude smaller than 1 and can be safely neglected. We translate the bounds on (Ṗ b /P b ) BNV into limits onḂ/B and report our results in the last row of Table 2. We can write an expression for the derivative of baryon number,Ḃ = f × B × Γ BNV , with f being the proportion of the baryons involved in the BNV process. If we assume that most of the matter inside the neutron star has densities above the required BNV threshold we would have f ≈ 1 for neutrons, and f ≈ 10 −3 for hyperons. The limit on Γ BNV is then given by in which α = 0.4, 7, 1 for PSR J0737−3039A/B, the Hulse-Taylor binary, and PSR J1713+0747 respectively. This indicates that if BNV is active throughout any of these binary pulsars, its rate must be less than one per 10 10 yr, i.e., the characteristic lifetime for a typical pulsar. 9 In this section, we have studied the generic consequences of BNV processes in neutron stars. Our only assumptions have been that the rates for such processes are slower than the weak interactions, and their final products will ultimately, either directly or via a cascade of interactions, turn into the particles already present in the neutron star and with structure dictated by the standard BNC EoS. We have demonstrated that these processes can relocate the neutron stars along their one-parameter sequence away from the maximum mass configuration and have analyzed the rate at which this can occur. We have also shown that observations of neutron star properties such as the orbital periods of pulsar binaries can lead to stringent constraints on this generic class of BNV processes.
In the following sections we consider specific BNV mechanisms, as well as model realizations thereof, and their effects on neutron star physics. Three types of processes are of interest. The first concerns models in which baryon number is not violated but transferred to a hidden sector. Since these baryon-number-carrying particles are not observed, these processes appear to be violating baryon number. We consider this apparent BNV and its implications in the next section. In subsequent sections we consider explicit and spontaneous BNV in turn.

Implications of apparent BNV
The possibility of apparent baryon-number violation, with the concomitant notion that baryon number can also be carried by particles of a hidden sector, emerges naturally from the idea that ordinary baryonic matter and dark matter could share a common origin [149][150][151], since the dark matter relic density is within a factor of a few of that of ordinary matter. Thus the dark matter relic density would be set by its cosmic relic asymmetry much as in the case of the relic density of ordinary baryonic matter.
The earliest examples of such asymmetric dark matter models were in the context of technicolor models [150,151], and a number of mechanisms to explain the cosmic genesis of dark and visible matter have since been proposed -and we direct the reader to the review of Ref. [152] for a succinct summary.
We have noted that the neutron lifetime anomaly could, in principle, be resolved through the existence of a dark decay channel of the neutron, in which one or more or all particles in the final state are uncharged under the SM gauge groups [21,153]. There are many such possibilities. The neutron, e.g., could oscillate to a mirror, or dark, neutron [154] or a neutron could be destroyed through its interaction with asymmetric DM [155] -or it could decay to an exotic final state [21]. Alternatively, if a neutron star were to capture a Q-ball [156,157], noting that B-or L-carrying Q-balls can appear in supersymmetric extensions of the SM [158], then neutrons can be consumed by the Q-ball, increasing its baryonic charge, shortening the lifetime of the star [159]. Although these scenarios are all realizations of neutron disappearance, we say that the baryon number violation is apparent, rather than explicit, because baryon number remains unbroken. That is, in these scenarios the concept of baryon number has been generalized, as appropriate, to include particles in a hidden sector [152]. This construct is convenient in that (1) it permits the appearance of exotic decays of the neutron without incurring proton decay [21], (2) it enables models of cosmic baryogenesis through visible-hidden sector interactions [152,160], and (3) it can also provide a natural way of stabilizing a DM candidate. Thus the possibility of apparent baryon violation can be used more broadly, and we note Ref. [161] for a study of the possible decay channels in the presence of light particles with B or L, such as p → π + χ decay, which searches for p → π +ν limit severely, as we show in the next section. Earlier work concerning the possibility of dark neutron decay [162] via a mass-dimension-six "neutron portal" such as uddχ c L /Λ 2 n [160] gives rise to decays such as n → χγ and n → χZ d with Z d → e + e − , though the predicted rates are too slow to explain the neutron lifetime anomaly [21]. It strikes us that apparent baryon violating processes stand out among the BNV processes we consider in that they are not immediately required to occur very much more slowly than ordinary weak processes within the SM. Thus in this section, we focus on models of this ilk and consider not only how the existence of neutron stars constrains them but also how such constraints can potentially be evaded.
We open our discussion with a recap of made-to-measure models of the neutron lifetime anomaly [21], and we refer to Ref. [153] for a review. A model with operator O that mediates an exotic decay of the neutron can potentially admit proton decay vis-à-vis the same operator, through p → n + e − +ν e , with the virtual neutron n * decaying via the exotic decay channel. This is at odds with proton lifetime constraints, as we detail in Sec. 5, but it can be eliminated altogether if the mass of the exotic final state M f exceeds that of m p − m e , with nuclear stability imposing a still stricter constraint [21], namely that 937.993 MeV < M f < 939.565 MeV to prevent 9 Be → f αα decay [163]. The exotic final states contain at least one hidden sector particle χ, and χ can be a dark matter candidate if its mass is less than m p + m e . Nuclear decays open the opportunity of studying neutron decays to invisible final states, and we return to this point later in this section.
Alternatively, e.g., a Φ with (3, 1) 2/3 would also work [21], though it cannot couple to two first-generation d-like quarks; n → χγ can, however, appear via the strange quark content of the neutron [21] or through a one-loop weak process [165]. In contrast, the realization of the invisible decay n → χφ, with the particle content thus far noted, requires the introduction of another Dirac fermionχ as well [21]. Finally the decay n → χχχ with B χ = 1/3 can be realized at the nucleon level via [164] where B χ = 1/3 and Γ is a combination of V ± A interactions. For Dirac fermions χ, a choice which makes "short-cut" |∆B| = 2 transitions via dark neutron decay impossible, empirical constraints from direct searches for Φ at colliders, allow a loci of parameters that would permit the resolution of the neutron lifetime anomaly with new physics, that is, a dark branching ratio of ∼ 1% [21]. However, direct searches for n → χγ [53] and n → χe + e − [54,55] decay processes significantly constrain the allowed possibilities, particularly if χ is light enough to be stable. In particular, the study of Ref. [53] does not constrain E γ < 0.782 MeV. However, studies at Borexino [166] can also be employed to the same end, limiting Br(n → χγ) < 10 −4 , as noted by Ref. [167], thus precluding this particular solution to the anomaly, though dark H decay is still possible [168,169]. In addition, dark neutron decays, at the strength to explain the anomaly, particularly that of n → χγ, can render a massive neutron star unstable [22][23][24], limiting its maximum mass to 0.8 M , which is inconsistent with observations [139,170]. If χ has repulsive self-interactions, then such effects can make dark neutron decays to final states with χ energetically less favorable and ultimately permitting neutron star masses that are consistent with observations. We note Refs. [22,171] for models in which such self-interaction effects has been studied. Interestingly, Ref. [164] has shown that it is possible to evade these constraints if only the dark decay mode n → χχχ is permitted -indicating that a new physics solution to the anomaly would be possible, though the β-decay constraints studied in Sec. 2 would seem to limit its role.
Thinking broadly, we emphasize that dark decay models [21,161,162,165,172], even if they ultimately make a negligible contribution to the neutron lifetime anomaly, nevertheless allow for much larger apparent BNV effects than that permitted from direct searches for explicit BNV. For example, the minimal dark sector model with Φ in the (3, 1) 2/3 representation can mediate n → χγ up to the ∼ 10 −6 level, with a Φ at the TeV scale, opening the possibility for its discovery at the LHC [165]. Moreover, this model permits Λ → χγ, for which there are no direct constraints, and it is possible to trade the size of n → χγ for Λ → χγ, or vice versa, given the existing constraint from D 0 −D 0 oscillations [165]. In the current context, the notion of differing rates for n → χγ and Λ → χγ, predict very different evolutionary effects within the neutron star, as we have studied in Sec. 3. In particular, the possibility of Λ → χγ decay is only appreciable at central densities for which a Λ population appears, whereas n → χγ decay has no such requirement. Neutron decays into dark final states can also potentially be discovered or constrained through the study of nuclear decays [21,163,173]. One possibility concerns the study of β-delayed proton emission in 11 Be decay, 11 Be(βp) [163], though the estimated n * → χ rate in the nucleus (the γ is not needed in the nuclear process) appears to be significantly larger than the empirical width of 11 Be, significantly constraining this solution to the neutron lifetime anomaly [173]. It is nevertheless the case that a surprisingly large branching ratio for quasi-neutron-like decay has been inferred from the detection of 10 Be in 11 Be → 10 Bepe −ν e decay [174], perhaps the neutron also decays invisibly in this process [163]. This process has been investigated further, with direct observation of the final-state protons [175] confirming the size of the branching ratio from the earlier indirect result [174], even if a subsequent experiment [176] in the manner of Ref. [174] fails to do so. It has been noted that a new resonance state could explain the large branching ratio, and this appears to be possible theoretically, both from direct study of the resonance properties [177] and from a study within halo-nucleus effective field theory [178]. Without a resonance, the decay rate would be very challenging to explain [179,180]. These studies constrain the size of a possible dark decay of the neutron in reference to the neutron lifetime anomaly, but have little impact on the broader possibilities we consider here.
Alternative explanations for the neutron lifetime anomaly come from the possibility of dark-matter-neutron interactions in the bottle experiments [181], the existence of a dark Z d [89], or from the possibility of neutron-mirror-neutron mixing [52,154,182], where we note Ref. [183] for a theoretical review of the latter set of models. The last possibility is limited by neutron star heating constraints [184,185], direct experimental searches [186][187][188][189][190][191][192], as well as pulsar timing studies [193], though the possibility of "hidden" magnetic fields [52,194] makes for a large phenomenological parameter space to explore [195][196][197].
We now turn to the study of explicit BNV.

Implications of explicit BNV
Baryon number is only an accidental symmetry of the Standard Model: given its particle content, the complete set of possible renormalizable interactions conserves this quantity without requiring it a priori. This concept has been invoked to explain the apparent stability of matter, analogous to how the electron is stabilized by electric charge conservation [198][199][200][201]. Neutron decays into new states, considered in the previous section, apparently violate B, but this accidental symmetry can be readily generalized to incorporate the new particles so that they do not. However, a generic new physics scenario need not conserve baryon number. Proton decay into SM states, in particular, is a generic consequence of grand unified theories (GUTs) [202,203] and of models of supersymmetry [204,205]; see, e.g., Refs. [206][207][208][209][210]. Given that baryon number must ostensibly be violated on some level in order to generate the observed baryon asymmetry of the universe (BAU) [25], it is pertinent to consider how this might appear at low energies. 10 We adopt a model-independent perspective on low-energy baryon number violation, considering the operators, comprised solely of SM fields, that generate some subset of BNV processes with no assumed relationships between them. This approach dates to Refs. [216][217][218][219][220] and can be recast in the language of Standard Model Effective Field Theory (SMEFT) [221][222][223][224][225]. If heavy new particles exist, then these can only manifest at low energies through the tower of operators that they engender, including nonrenormalizable ones. These operators must be invariant under the Lorentz group and the SM gauge group, but need not preserve accidental symmetries. The nonrenormalizable part of the SMEFT Lagrangian, L NR , can be generically expressed as follows: where O (d) is an operator with mass-dimension d. The c i are defined to be dimensionless; the factor 1/Λ d−4 is required by dimensional consistency and reflects the intuition that higherscale physics contributes with lesser strength to low-energy processes, i.e., that ultraviolet physics decouples from physics in the infrared [216]. On one hand, operators at high massdimension will be suppressed relative to those with lower mass-dimension if they share a common scale Λ. On the other, processes with small rates may be connected to relatively low-scale physics if they are dominantly generated by higher mass-dimension operators. We categorize these operators by the extent to which they break B, i.e., the ∆B of the operator, starting with operators with |∆B| = 1 before moving on to those with |∆B| = 2. Processes with larger violations have received significantly less theoretical and experimental interest, though some constraints do exist [226][227][228][229]. 11 We note, however, that these operators would be suppressed by relatively high powers of the scale Λ, allowing for small rates of baryon-number violation to be connected to relatively low-scale new physics.

Processes with |∆B| = 1
These first arise at mass-dimension six 12 within SMEFT and are schematically of the form L where q represents a quark operator and represents a lepton operator; the sum over i represents the sum over possible chirality and gauge structures. It is this class of operators that can give rise to proton decay and nonstandard neutron decay. Restricting to the first quark generation, there are four possible operators for each generation of lepton (= e, µ): On the left, we write the operators before electroweak symmetry is broken in terms of two-component (Weyl) spinors: the left-handed quark doublet, Q; the left-handed lepton doublet, L ; the right-handed singlet up quark, u; the right-handed singlet down quark, d; 11 We note electroweak sphalerons, nonperturbative gauge configurations that can convert baryon number into lepton number (and vice versa) in units of three. These processes are only operative at high temperatures and naturally allow for BNV in the early universe; these are inoperative at the energy scales we consider, and we refer to Ref. [230] for a detailed study in SM lattice gauge theory, though we also note the exploration of possible exceptions [231][232][233]. 12 One often reads about contributions that arise at dimension five [204,205] within supersymmetric models, not six. This is in a framework in which the superpartners are dynamical -if one were to integrate these out, as would be the case in SMEFT, then the resulting effective interactions would be no lower than dimension six. Figure 7: Constraints on the partial lifetimes of representative two-body BNV nucleon decays: p → + π 0 [234], p → ν π + [235], p → + η [236], p → + γ [237], n → + π − [236], n → ν π 0 [235], n → ν η [238] and n → ν γ [239], with = e, µ. All values represent 90% C.L. upper limits. and the right-handed singlet charged lepton e ). Fermion fields enclosed in parentheses have been grouped into Lorentz scalar bilinears and ε is the antisymmetric tensor, which is used to form weak singlets. On the right, we write them in the broken phase in terms of the more familiar four-component spinors; ψ c = ψ T C is the charge conjugate of the field ψ and P R,L are the standard projection operators. Color contractions have been left implicit; only the totally antisymmetric combination may be used. The task becomes constraining the coefficients c i /Λ 2 |∆B|=1 for these operators. Constraints on the proton lifetime date back to the 1950s [207] and the current limit on the total lifetime is 3.6 × 10 29 yr [240], though some partial lifetimes have been measured at the level O(10 32 − 10 34 ) yr; we refer the reader to the comprehensive tabulation of partial widths in Ref. [241] and references therein. A subset of two-body decay lifetimes are shown in Fig. 7. One can interpret these limits in terms of the operators in Eq. (45); we briefly outline this analysis, following the formalism established in, e.g., Refs. [242][243][244]. The decay width for N → L Π, with N a nucleon, L an antilepton of flavor = e, µ and Π a pseudoscalar meson, is given by 13 where m L , p L , E L are the antilepton mass, three-momentum and energy, respectively, c ( ) χχ 13 We write Λ instead of the more cumbersome Λ |∆B|=1 here.

28
is the coefficient of the operator O ( ) χχ and W χχ 0,1 are form factors for the matrix element LΠ|O ( ) χχ |N , which are calculated on the lattice. If L = ν, then contributions proportional to c χR must be taken to zero. We use the nucleon-pion form factors from lattice calculations with a physical pion mass of Ref. [244]; nucleon-eta form factors are taken from calculations with unphysical pion masses in Ref. [243]. Moreover, we calculate the decay widths for N → Lγ to find where Q N , µ N are the nucleon charge and magnetic moment (the latter in units of the Bohr magneton) and α, β are low-energy constants that characterize the proton-to-vacuum transitions, which we obtain from Ref. [244].
The results of this analysis are shown in Fig. 8. Allowing for one operator to be active at a time for simplicity and fixing the corresponding c Consequently, |∆B| = 1 processes are severely constrained; this is sufficient to rule out the simplest GUT models. We note, however, that there may be cancellations among the c ( ) χχ , suppressing rates of some subset of decay modes; this requires fine-tuning, but is logically possible. However, the nonobservation of at most one channel can be explained by such a cancellation 14 -we are inexorably led to the conclusion that proton decay at mass-dimension six requires the scale of new physics to be not far below the Planck scale. That said, while the tower of operators starts at d = 6, all higher mass dimensions also include |∆B| = 1 operators. It is logically possible that dimension-six contributions to proton decay could be suppressed in a given model and that, say, dimension-eight contributions could dominate. These latter contributions will depend on higher powers of 1/Λ |∆B|=1 , thereby allowing this scale to be significantly less than the ∼ 10 16 GeV figure we have found.
Nucleons may also decay into final states with kaons; there exist models [245][246][247] in which symmetries are invoked to suppress decays to first-generation fermions, thereby leading to decays involving the heavier mesons. These channels are only slightly less constrained than the strongest first-generation final states [30]. There are six operators that can engender these decays for each generation of lepton. The constraints on the associated scale turn Figure 8: Representative constraints on operators with |∆B| = 1. We employ the lattice matrix elements, form factors, and low-energy constants of Refs. [243,244]. Note that the operators O ( ) χR do not generate decays to antineutrinos; these are thus not constrained by these final states. out to be one or two orders of magnitude weaker than those shown above. We omit these due to space constraints, but, irrespective of the sensitivity of the experimental limits, we emphasize that the inferred scale of these operators need not be the same as the scale of the first-generation operators considered previously. In particular, if the coefficients of SMEFT, determined by the true high-scale theory of Nature, depend on quark flavor, then it might be reasonable to expect that processes involving the light quarks would be suppressed, and that new-physics effects are more likely to appear in heavy quark systems.
We note an important theorem regarding SMEFT [248,249]: operators with even (odd) mass dimension induce changes to B −L in even (odd) multiples of two. The mass-dimension six operators in Eq. (45) thus conserve B − L. However, it need not be the case that proton decay conserves B − L: there exist four BNV operators at mass-dimension seven that all violate B − L by two units. The importance of the (non)conservation of B − L in proton decay was recognized as early as the late 1970s [250], where it was noted that B −L breaking would imply the existence of a new scale between the weak and GUT scales. Here, we note that if proton decay were both B-and B − L-violating -say, via p → e − π + π + -then the dependence on 1/Λ 3 allows for the scale of proton decay to be much smaller than the values obtained in our analysis. We coarsely estimate this as This scale is far larger than directly accessible experimental energy scales, but is many orders of magnitude suppressed relative to typical GUT and Planck scales. Lastly, we note that oneunit violations of B can manifest in processes besides nucleon decays: scattering processes such as pe − → e + e − or pn → e + n may be generated at even mass dimension, whereas processes such as nn → pe − , which violate B − L, may occur at odd mass dimension. We emphasize that there is no innate reason why baryon number should necessarily be violated by only one unit at a time. The study of |∆B| = 1 operators and the processes they engender is motivated by the observation that in the absence of non-SM states that there are no other means by which protons can decay. It would be a remarkable signature, to be sure, but it is not guaranteed that proton decay should occur if B is explicitly violated -the true theory of Nature may be such that it is dramatically suppressed relative to the intrinsic BNV scale or that it is for some reason completely forbidden. While the large scales discussed in this subsection may be discouraging from the point of view of phenomenological relevance, it is critical to keep in mind that proton decay is not the only way in which BNV can manifest at low energies. In the remainder of this section, we aim to demonstrate that B violation may still be accessible in both the laboratory and in astrophysical settings and we turn to the consideration of processes with |∆B| = 2.

Processes with |∆B| = 2
We turn now to |∆B| = 2 operators. These first appear in SMEFT at mass-dimension nine and are of the form The Lagrangian in Eq. (53) differs from Eq. (45) in that the former does not conserve B − L. This is important because if |∆B| = 2 processes were subject to the same high scale as |∆B| = 1 processes, then the suppression by higher powers in Λ would imply that the former would be impossible to observe, for all practical purposes. However, as we have noted, if B − L were broken, then this would imply that new dynamics must be operative between the weak and unification scales [250]. One may construct models [250][251][252][253][254][255][256][257][258][259][260][261] that invoke a new scale at which B − L is broken, such that |∆B| = 2 processes are generated while generating small (or vanishing) contributions to proton decay. If Nature realized this scheme, then it would be quite a blow to the GUT paradigm -the introduction of new states far below the GUT scale may spoil the promise of unification, or at least render such a framework hopelessly unpredictive. Still, the promise of lowering the scale of new physics to (potentially) accessible energies makes this a tantalizing possibility. This invites us to take seriously the possibility that Λ |∆B|=2 might be much smaller than the assumed GUT scale deduced from experimental limits on p decay. A crucial observation is that |∆B| = 2 operators do not contribute to free nucleon decay, 15 immediately eliminating a powerful set of potential constraints. The most well-studied effect of these operators at low energies is to mix neutrons and antineutrons, allowing for spontaneous oscillations between the two. Ref. [263] provides a comprehensive review of this subject (see also Refs. [264][265][266]); here, we only briefly touch upon the most phenomenologically relevant details. The operators of Eq. (53) generate the following contribution to the nucleon-level Lagrangian: with n the neutron field and n c = Cn T its charge conjugate. This induces nondiagonal Hamiltonian matrix elements for the n − n system. Consequently, the probability P nn (t) for a neutron to manifest as an antineutron after a propagation time t is where τ n is the free neutron lifetime. We identify (δm) −1 as the timescale of nn oscillations; we employ the standard notation τ nn = (δm) −1 [263][264][265][266]. Searches for free nn oscillations have been performed at nuclear reactors dating back nearly forty years [267][268][269]; the leading limit on the oscillation timescale using free neutrons comes from the ILL [270]: To interpret this result in terms of the scale Λ |∆B|=2 in Eq. (53), we note that δm can be parametrically written as for C ∼ O(1) and Λ QCD = 1 GeV, this implies Λ |∆B|=2 O(100) TeV. One can improve this estimate using more precise calculations of the n − n matrix elements, including those of lattice QCD [271,272] or phenomenological models such as the MIT bag model [220]. These are useful for interpreting specific UV model predictions, but do not lead to qualitatively different conclusions, for our purposes. Therefore, we eschew an effective-operator analysis of the sort we have discussed for |∆B| = 1 processes. It has been nearly 30 years since a new search for free nn oscillations has been performed. That said, the coming decade-plus promises an improvement of more than an order of magnitude in τ nn (∼ 3 × 10 9 s) from NNBAR at the European Spallation Source [265]; this implies a sensitivity to Λ |∆B|=2 at the PeV scale. This scheme can be made arbitrarily more complicated: • Explicit and apparent BNV are not mutually exclusive. If mirror neutrons exist and |∆B| = 2 dynamics are operative, then one can consider simultaneous oscillations among all four of n, n, n , n . It has been proposed that nn oscillations may operate through a "shortcut" [273] through the mirror neutron, via n → n /n → n: these transitions may occur even if the nn matrix element is small. This framework can be difficult to test due to the ill-constrained Hamiltonian involving the interactions of the mirror neutron. Proposals exist to probe this scenario at the European Spallation Source [265], at Oak Ridge National Laboratory [274], and at the Paul Scherrer Institut [275].
• Eq. (55) assumes that the n and n are exactly degenerate in the experimental environment. There are several reasons why this would not be the case; among these are the presence of external matter and magnetic fields [276,277], which lift the degeneracy and suppress oscillations. While these suppress oscillations, they may also stimulate |∆B| = 2 phenomena via other mechanisms, such as neutron-antineutron conversion via scattering, such as ne − → ne − . This may be generated either by long-distance contributions (as in, e.g., Ref. [278]) or through short-distance contributions in SMEFT -these processes can select different subsets of operators from those that generate nn oscillations directly. The latter come with higher inverse powers of the new scale; this may allow for Λ |∆B|=2 to be lowered to the O(1 − 10) TeV scale without running aground of existing constraints.
While |∆B| = 2 operators do not engender new nucleon decays, these can render some nuclei unstable through dinucleon decays, N N → X, where X is some state comprised of leptons, mesons and photons [279]. The nuclear decay lifetimes can be related to the free oscillation timescale via [280] τ dinuc. = R dinuc.,nn × τ 2 nn , where R dinuc.,nn accounts for the effects of nuclear structure: ambient nuclear matter breaks the degeneracy between n and n, thereby suppressing the probability for a given neutron to convert. Searches for these decays have been performed using deuterium [281], 16 note that this is stronger than the exclusion from searches with free neutrons. DUNE proposes to measure the dinucleon decay lifetime of 40 Ar at the level of 6.45 × 10 32 s −1 for a 400 kt·year exposure [290]; using R Ar dinuc.,nn = 0.56 × 10 23 s −1 [291], this implies sensitivity to τ nn 6.0 × 10 8 s.  Figure 9: A sketch of a dimension-13 operator contribution to nn oscillations. On the left, we show a cartoon of an operator of the form (udd)(udd)(hee). The Higgs field h is required by gauge invariancethis operator must flip the electron chirality in order to be Lorentz invariant, so that the Higgs is required to preserve all SM charges. On the right, we have closed the electron and Higgs legs to form a two-loop nn diagram. This diagram is suppressed by the electron Yukawa coupling, y e ≈ 2 × 10 −6 (denoted by the white circle), as well as two loop factors, (16π 2 ) −2 . Therefore, for some scale of new physics Λ |∆B|=2 , nn oscillations would be suppressed relative to naïve expectations.
An important caveat is that the dinucleon decay lifetime depends on how |∆B| = 2 phenomena manifest at low energies. The values presented above assume that the dominant mechanism is nn oscillations, resulting in annihilation with a spectator nucleon. In truth, the connection between τ dinuc. and τ nn depends on the (assumed) relative strengths of all operators that can generate |∆B| = 2 processes -including those that do not generate nn oscillations at tree level. As an example of this, suppose the leading contribution to the SMEFT Lagrangian from high-scale BNV physics were of the form L ∼ 1 Λ 9 (13) (udd)(udd)(hee); (60) in other words, assume that some UV model gives rise to this dimension-13 contribution at tree level, but that the typical dimension-nine operators were, for some reason, absent. If so, then the process nn → e + e − occurs at tree level, but nn oscillations only arise at two-loop level; this is demonstrated in Fig. 9. We estimate the contribution of this operator to δm to be where y e ≈ 2 × 10 −6 is the electron Yukawa coupling. The prefactor is significantly less than unity; if one tried to associate this with a dimension-nine operator as in Eq. (57), then one would misinterpret the scale of the associated physics by about an order of magnitude -to wit, current limits on τ nn imply Λ (13) ∼ O(10TeV) as opposed O(100TeV. Within the nuclear medium, τ dinucleon would receive contributions both from nn oscillations with subsequent annihilation, and from the direct reaction nn → e + e − . There is no reason to expect, a priori, that the latter should be suppressed relative to the former. We can generalize Eq. (58) to the precise relationship between R dinuc.,nn and R dinuc.,nn→e + e − depends on the details of nuclear structure and on the precise connection between nn → e + e − and nn oscillations. Moreover, |∆B| = 2 operators that conserve B −L do appear starting at mass-dimension 12. 16 This class of processes is interesting because they contain contributions that are second order in the |∆B| = 1 operators discussed in the previous subsection. If, for instance, p → e + π 0 were present, then processes such as pp → e + e + and e − p → e + p [292], as well as hydrogen-antihydrogen oscillations [293][294][295], must be, too. However, the rates or cross sections for these processes scale as Λ −8 |∆B|=1 ; given the strong constraints on the associated scale, these must be suppressed to a fantastic degree. However, while the existence of p → e + π 0 is sufficient to generate these processes, it is not necessary -the model building and phenomenology become much richer if one considers more general mechanisms of BNV. In particular, if these were connected to physics that innately provides for violations of baryon number by two units while conserving B − L, then these may be operative in the absence of proton decay. As a concrete realization of these ideas, we highlight the minimal scalar models studied in Ref. [296]. The catch is that because the associated operators are dimension-12, the rates must scale as Λ −16 |∆B|=2 . However, the upside is that nonobservation of any of these processes results in a constraint on Λ only around the TeV scale [292]; put more optimistically, there is still room for new, BNV physics to exist in a way that can be discovered at future collider experiments.

Connections to Lepton Number Violation
We have seen that |∆B| = 2 operators may also violate B − L. This correspondence is not necessary, because the operator need not have an odd mass dimension [248]. Both nν → nν scattering and hydrogen-antihydrogen oscillations, e.g., arise from dimension-12 operators that violate both B and L by two units but preserve B − L; this suggests that there may be connections between BNV processes and particular lepton-number violating (LNV) processes. In particular, we are free to wonder if nn oscillations could possibly be connected to Majorana neutrino masses, similar to how the Schechter-Valle theorem connects neutrinoless double-β decay (0νββ) to Majorana neutrino masses [298], because the Lagrangian in Eq. with, e.g., n → e − π + → e − e + ν and F = e + e − [297]; b) with n →νF , a |∆(B − L)| = 0 process, with, e.g., n → e + π − → e + e −ν and F = e + e − (Ref. [297] considers the process p → e + π 0 , which is related by isospin.); c) with FF oscillations noting that F is restricted to have B − L = 0 only. If F = e − p, then, e.g., e − p → e +p must occur for nn oscillations to imply the existence of 0νββ [296], or vice versa. In the context of a NS, these connections give rise to nn → νν (νν), or nn → F F . See the text for further discussion.
operator analysis of Ref. [297] determines that the observation of any two of p → e + π 0 , n → e − π + , and nn oscillations would also show that 0νββ decay can occur. In contrast, in Ref. [296], minimal scalar models from Ref. [256] are used to show that the observation of a |∆B| = |∆L| = 2 scattering process, along with that of nn oscillations, would also show that 0νββ can occur -and if 0νββ decay can occur, then a Majorana neutrino mass also exists [298]. These ideas can be generalized to realize three distinct mechanisms for connecting neutron and neutrino Majorana masses, illustrated in Fig. 10. That is, nn oscillations can combine with either n → νF decay, or n →νF decay, where F is a state of SM particles with total B = L = 0 and zero electric charge in each case, or with F −F oscillations, where F has B − L = 0, to give rise to a Majorana neutrino mass. In the first instance, we could have n → e − π + → e − e + ν, so that F = e + e − , but F = µ + µ − or F = π, η, or simply γγ are also possible, as is F = K. Analogously, in the second case, we could have n → e + π − → e − e +ν , to yield the same set of F we have just enumerated. In the last case, with SM neutron β decay we have F = e − p, so that BSM physics is required to yield e − p → e +p [296]. Ultimately, we see that the various connection processes would be signalled by the appearance of dineutron decay to F F orνν final states.
In conclusion, we emphasize that BNV processes can only be connected to LNV processes through operators that violate both. Although the connections we have noted could be discoverable, the resulting rates for LNV from BNV and vice versa are expected to be extremely small [297].

Effects in Neutron Stars
For a neutron star containing ∼ O(10 57 ) neutrons with a BNV decay lifetime of ∼ O(10 30 ) years, corresponding to the most optimistic limit on the proton lifetime [240], one would expect ∼ O(10 27 ) such decays to occur in a year. This relatively small rate is unlikely to lead to observable changes in its macroscopic observables (mass, radius, moment of inertia, etc.), and any possible signature is expected to be too weak to observe. In addition to decays, though, there may be |∆B| = 1 scattering processes operative within the star, e.g., ne − → νe − . While these processes would only expedite the conversion of neutrons into energy, it is unlikely that they would lead to observable signatures, either. Fundamentally, this is a consequence of the stringent constraints from proton decay experiments: the limits on Λ |∆B|=1 are just too strong to allow for observable consequences at such low energies. Recall from our discussion in Sec. 3, however, that binary spin-down considerations tolerate much larger rates of BNV: taking α ≈ 1, f ≈ 1 (for nucleons) in Eq. (41) yields an upper limit of Γ BNV 10 −12 yr −1 . We therefore allow ourselves to contemplate rates this large by allowing for violations of B by more than one unit, though apparent BNV could also act.
How |∆B| = 2 physics manifests in neutron stars is expected to be qualitatively different from how it would manifest in the laboratory. Matter effects break the degeneracy between neutrons and antineutrons in neutron star matter, quenching free nn oscillations. Therefore, we do not expect to be able to associate any phenomena directly with the timescale τ nn for free neutrons. This is something of a pity -the fact that pulsars with characteristic ages of O(10 8 − 10 10 ) years have been observed (see, e.g., the pulsar catalogs of Refs. [299][300][301]) would surely be able to probe phenomena that occur on the timescale τ nn , which is only constrained to be O(3) years [270]. However, it seems this is just not the universe that we occupy. This is also the case in nuclei, as discussed above, but matter effects are stronger at the supranuclear densities present in the cores of neutron stars. As such, the rate of nn oscillations will be dramatically suppressed in neutron stars, even compared to nuclei.
These |∆B| = 2 operators can manifest as two broad classes of processes: • Processes that destroy two nucleons. These include processes such as nn → 2γ, 3ν, etc. This is completely analogous to dinucleon decay, discussed above.
• Processes that convert nucleons to antinucleons. These include scattering processes such as e − n → e − n.
These are not mutually exclusive -the process nn → nn belongs to both. In the remainder of this discussion, we focus on processes involving neutrons. Protons are present in the star at the ∼ O(10%) level; these may also participate in BNV interactions with the ambient matter, e.g., e − p → e + p, but we expect neutrons to dominate the dynamics. On one hand, |∆B| = 2 processes should scale with δm (see Eq. (54)) at the amplitude level, assuming that they depend on the same underlying mechanism of B violation. The partial lifetime for any such process thus generically scales as R NS depends both on properties of the nuclear medium and, as in the case of dinucleon decays in nuclei, on which BNV processes are operative within the NS. On the other, R NS need not be of the same order of magnitude as its counterpart in nuclear matter, ∼ O(10 23 ) s −1 , given the larger densities, the lower proton fraction and the requirement of charge neutrality, allowing charged species, such as electrons, muons, pions, and kaons to appear in appreciable numbers as well.
We can place weak limits on the size of τ NS , and thus on R NS . Based on our upper limit from Eq. One would need a more complete analysis -including the effects of the neutron star EoS and a concrete set of BNV interactions -to derive more robust limits than this, but this estimate speaks to the power that the mere existence of cold, old neutron stars has on constraining this sort of new physics.
Let us emphasize a key finding: the present limit on the temperature of PSR J2144-3933 gives a stronger constraint than that derived from the nonobservation of anomalous binary spin-down. Indeed, if B violation were operative and saturated the upper limit on the rate in Eq. (41), then this would imply an energy production rate of ∼ O(10 34 − 10 35 ) erg s −1 from the destruction of baryons; this corresponds to ∼ O(10 − 100) L . If most (or all) of this energy is trapped by the NS, then this would lead to significant heating: the corresponding asymptotic temperature would be O(10 6 ) K, in clear violation of the limit on J2144-3933. This constraint is only operative, however, if the products of the BNV reaction cannot escape the star. For species that interact electromagnetically, this is clearly the case. For neutrinos, this is less clear, a priori. However, recall that in the discussion surrounding Eq. (21), we concluded that a ν e with energy E ν 5 MeV is likely to scatter on its way out of the core, where it is most likely to have been produced -neutrinos produced in BNV processes have much higher energies than this, so that even these would need to deposit most of their energy before they can escape. As such, it would only be possible for BNV rates to be this large if the decay products were invisible to the rest of the star, as in our discussion of apparent B violation in Sec. 4. Put more sharply, decay processes such as n → χχχ [164] would evade constraints from neutron star heating.
Explicit BNV processes invariably lead to production of (anti)neutrinos from a combination of three sources: 1. Those produced directly in the BNV reaction. 2. Those produced by the weak reactions that restore chemical equilibrium after some BNV process has disrupted it (i.e., Urca reactions). 3. Those emitted as a result of the heating of the star, via processes such as N N → N N νν.
Those of the first category are directly sensitive to the relationship between B and L violation. To wit, if L is conserved in B-violating processes, then these should produce neutrinos and antineutrinos in equal numbers. If, however, B is violated in such a way that B − L is conserved, then these reactions must produce antineutrinos in excess of neutrinos. The dependence of the second category on L (non)conservation is more difficult to pin down -it depends on the precise connection between the BNV processes and the weak disequilibrium they engender. 17 The third type are, by assumption, SM processes, which conserve L by default; neutrinos and antineutrinos must be emitted in equal numbers. We will discuss the observability of (anti)neutrino signals associated with B violation in more detail in Sec. 8. Lastly, we note that BNV processes need not restrict themselves to protons and neutrons in such an environment. At high densities, hyperonic degrees of freedom may emerge within the cores of neutron stars; for recent reviews on the subject, see Refs. [303][304][305]. These may also participate in BNV interactions, either among themselves or with the nucleons [279], that are, at best, poorly constrained. These cannot be directly probed in the laboratory, and only indirect comparisons with the operators controlling the nonobservation of nn oscillations or N N → kaons can be formed. Even less well understood are the contributions of the spin-3/2 ∆ resonances, which may also appear in the cores of neutron stars [303] (see also Ref. [306], which aggregates predictions of neutron star properties for several EoS involving ∆s), and might decay to antinucleons (e.g., ∆ → N π). Appreciable amounts of non-nucleonic hadrons within neutron stars may yet allow for fantastic signatures of B violation at the extremes of matter.

Implications of spontaneous BNV
If baryon number were a gauge symmetry of nature instead of an accidental one, then null results from fifth force searches (e.g., [307][308][309]) imply that it is unlikely to be a symmetry of the vacuum that we occupy -but it could be spontaneously broken. This notion has existed nearly as long as the concept of baryon number itself . The mechanism of this breaking may leave an observable imprint on the low-energy spectrum of the theory. If an analogue of the Higgs mechanism within the SM or pion condensation in QCD were operative within the baryon-number sector, then there might be (fundamental or composite) scalars present. Even if the mechanism of breaking is not directly observable at low energies, the low-energy spectrum of the theory may contain more than simply the gauge boson. If one were to embed baryon number into a larger gauge group [332][333][334], or incorporate quark flavor into the gauge structure [335], then there may exist additional gauge states or fermions. In this article, we focus on the existence of a single, new gauge boson. New scalars and fermions are interesting in their own right, but these can exist irrespective of whether or not baryon number is gauged -a massive new gauge state would be a smoking-gun signature of the gauging (and breaking) of B. We will call this state X: its gauge coupling is g X , and its mass is m X .

Laboratory Constraints on a New Gauge Boson
We are particularly interested in the effect of X on two-nucleon interactions. The new vector state necessarily generates a repulsive contribution between nucleons; how this contribution compares to the strong nuclear force depends on m X . At one extreme, the new force could be light enough that its range exceeds that of the nuclear force, the latter being set by m −1 π ∼ O(1 fm). Constraints on a new gauge boson in the range m X ∈ [10 −3 , 10 9 ] eV from nucleon dynamics are shown in Fig. 11. In particular, we show the following: 18 • A selection of fifth force searches have been aggregated in Ref. [309]. We have converted Figs. 8 and 9 of this reference from the α − λ parameter space to g 2 X /4π − m X and plotted the result in dark blue.
• Neutron scattering and neutron optics constraints are taken from Fig. 2 of Ref. [336] (see also Ref. [343]), and are shown in black and gray, respectively.
• This new interaction would also change energy levels, relative to QED predictions, of the antiproton-helium (p-He) bound state. The constraint from Ref. [337] is shown in blue.
• New forces also change the charge radii and binding energies of nuclei. Ref. [338] studies the effects of new nuclear-range interactions on 48 Ca, 120 Sn, and 208 Pb; their constraint is shown in light blue.
• This interaction would modify the long-distance potential between two protons in the Sun, thereby altering the rate of solar fusion. This could then (1) change the inferred age of the Sun to be inconsistent with the age of the solar system, and (2) modify solar neutrino production to be inconsistent with observations. Constraints have been derived in Ref. [339]. These are shown in shades of pink, corresponding to different proton energies: 10 (light), 50 (medium) and 100 (dark) keV. Figure 11: Constraints on the coupling and mass of a new boson associated with gauged baryon number, U (1) B . Aggregate constraints from fifth-force searches (dark blue) are taken from the review of Ref. [309]. Neutron scattering (black) and neutron optics (gray) constraints are from Ref. [336]. The p-He constraint (blue) is from Ref. [337]. The nuclear charge constraint (light blue) is from Ref. [338]. Constraints from solar fusion [339] are shown in various shades of pink; see the text for details. Constraints from the anomalous cooling and trapping from SN1987A are shown in orange [340]; we emphasize that only the shaded region is excluded. Anomalous cooling constraints from NS1987A [341] and Cas A [342] are given in red and dark red, respectively. The gray-shaded region corresponds to a new force whose range is comparable to the strong nuclear force; constraints in this region are presented in more detail in Fig. 12.
At the other extreme, the new state might be so heavy -that is, short ranged -that it cannot contribute in any meaningful way to low-energy nucleon processes. If this is the case, then these interactions contribute to contact terms in the chiral Lagrangian, i.e., to the so-called low-energy constants (LECs) of chiral effective field theory (χEFT) [344][345][346][347].
One would expect this if m X Λ QCD ∼ O(1) GeV -it may not make sense to talk about nucleons for momentum exchanges much larger than this scale, but these heavy states can be probed at colliders [324,327,348].
In between these domains -that is, for m π m X Λ QCD -is the regime in which the new interaction is not so short-range that it can be integrated out of the N N force, but not so long-ranged that its contributions can be clearly distinguished from nuclear forces. This latter aspect is particularly confounding, because we lack a precise, first-principles description of the forces between nucleons, though lattice QCD may yet provide one [349].Historically, interactions among nucleons have been described using phenomenological models [350][351][352][353]: one introduces a set of interactions with unknown coefficients that are fixed by low-energy nuclear data, including N N phase shifts and deuterium data. Even with more modern approaches such as χEFT, the unknown coefficients of the theory must be fit to data in order to put the framework to use. If the new degrees of freedom are not explicitly included in the N N potential, then their contributions could be inadvertently subsumed into some other part of the interaction -these may contribute to the short-range pieces of, e.g., the Argonne v 18 potential [352]. In the context of χEFT, these would contribute to LECs; however, if these new states are not too much heavier than the pion, then one would expect these to become dynamical for not-too-large momentum exchanges, even if they do not appear for low-momentum exchanges. This would manifest as an apparent inconsistency in the effective theory.
One further complication is that if two-body nuclear forces are poorly understood from first principles, then three-body forces (and beyond) are even more so -and these are not negligible. As with two-body forces, various phenomenological prescriptions exist for their inclusion [354][355][356][357][358][359][360][361] and they are naturally included within χEFT [344,345,347,[362][363][364]. However, these contribute to the uncertainty in these potentials, further complicating extractions of new-physics contributions to nuclear processes. When discussing the impact of new states on neutron star structure, we assume that the new interaction is abelian, so that it does not innately generate new three-body interactions, making the SM contributions to N N N forces overwhelmingly dominant.
Of course, the effects of such a new interaction would not be localized to nucleons. In particular, this new state vector contributes to radiative decays of light mesons [330,335,370]. These decays would essentially be two-step processes: 1. The meson decays radiatively to γ and X, e.g., η → γX. The computed rate of this process is scheme dependent: rates calculated at the quark level are different than rates calculated in, e.g., the vector meson dominance (VMD) scheme; see App. A.1 of Ref. [330]. 2. The X then decays into some observable final state. At tree level, X can decay to, e.g., π 0 γ or 3π. If there exists some nontrivial kinetic mixing with γ, then X may also decay into dilepton pairs, e + e − or µ + µ − , even though these are uncharged under B. Additionally, because X only couples to isoscalar currents, tree-level decays to π + π − are absent, barring either (1) nonzero kinetic mixing, or (2) more complicated gauge structures.
This width can then be compared against the width for the decay η ( ) → γγ. Constraints of this sort are shown in Fig. 12. The black curves are constraints derived from η → π 0 γγ (solid), η → π 0 γγ (dashed), η → ηγγ (dotted) and η → π + π − π 0 γ (dot-dashed); we collect the branching fractions in Table 3. Constraints may also be derived from decays of vector mesons [330]; these turn out to be weaker than pseudoscalar meson constraints, so we will not consider these further. On one hand, these new contributions increase the decay rates relative to pure-SM production. Conservative constraints can be derived by insisting that the new contribution not exceed the total observed value at some significance; this is what has been presented in Refs. [330,335,370]. One could try to strengthen these constraints by comparing with theoretical predictions of the rates [371][372][373][374][375] -in fact, there is ostensibly a discrepancy between measurement and the state-of-the-art prediction for η → π 0 γγ in Ref. [375]. However, Figure 12: Constraints on a sub-GeV vector boson, adapted from Refs. [330,331,335]. Black curves are from rare decays of pseudoscalar bosons: η → π 0 γγ (solid) [30,365], η → π 0 γγ (dashed) [366], η → ηγγ (dotted) [367] and η → π + π − π 0 γ (dot-dashed) [368]. The gray, double-dot-dashed curve represents the projected sensitivity to X photoproduction (γp → pX) for the full design luminosity of GlueX Phase IV [369]. The dark orange line is derived from the decay width for Υ(1S) → hadrons [317][318][319]; the light orange line is from the hadronic decay width for ψ(1S) [319]. The red-pink-white band represents the loci of points for which the new interaction increases the maximum neutron star mass by 0.1 − 0.5 M ; see Sec. 6.2 for details.
one should be circumspect in so doing: Ref. [375] employs a combination of the VMD and linear σ models to calculate decay widths under the assumption of isospin conservation with empirically-derived meson coupling constants. It is difficult to assess the systematic uncertainty incurred by these model choices; a dead-reckoning between theory and measurement may not be entirely reliable.
On the other hand, one can sidestep the issue of theoretical predictions and use a more robust experimental observable to constrain the existence of such a state. If X is a narrow state, then its decays induce narrow features in the invariant-mass distributions of the meson decays -in other words, bumps. It has been previously proposed to use bump hunts to probe new gauge states in Ref. [376] in the context of electron and proton beam-dump experiments; the same principles apply here, but with slightly different experimental config-urations. One might be able to conduct such a search with the upcoming JLab Eta Factory experiment [377,378] or with the REDTOP proposal [77,379,379]. Moreover, X can also be produced via γp → pX at GlueX [369], where such a search is also possible [331,335].
Also shown in Fig. 12 are constraints from hadronic decays of Υ(1S) [317][318][319] (dark orange) and of ψ(1S) [319] (light orange). These are sufficient to rule out new physics at the nuclear scale with strength comparable to electromagnetism (g 2 X /4π ≈ 1/137). However, we note that one can circumvent these processes [335] by insisting that the new physics only couples to first-generation baryon number, B 1 . In so doing, one must rederive the radiative pseudoscalar meson constraints assuming no mixing to the strange quark; the results are shown in Fig. 13. Apart from its tree-level couplings to quarks, the new vector state can kinetically mix with the SM photon. Limits on the kinetic mixing parameter ε have been compiled in, e.g., Refs. [380,381]. Usually, however, limits on ε are derived from searches for minimal dark photons, in which the new vector only couples to the SM through this kinetic mixing. One must reinterpret these constraints with the tree-level couplings to quarks from the outset; Fig. 6 of Ref. [380] has recast these searches in terms of limits on g X and m X . However, these limits assume that the kinetic mixing is given by ε = e 2 /(4π) 2 -otherwise, none of the constraints would be operative. We will not discuss these constraints in depth, but we note that, in the region 100 MeV m X 1 GeV, the kinetic mixing is most strongly probed by searches at LHCb for dimuon final states [382,383].
We also note that baryon number is anomalous within the SM -it is a symmetry of the Lagrangian, but not of the corresponding action. This is an acceptable state of affairs for global symmetries, but must be remedied for gauge symmetries by introducing additional fermions. From a model-building perspective, there is significant freedom in choosing how to resolve the anomalies, but in general, the existence of new fermions charged under baryon number (or some generalization thereof) can be probed at colliders [329]. If these new fermions are heavier than the electroweak scale, then integrating them out of the theory at low energies leads to three-gauge-boson interaction terms, XBB, where B is the gauge Decay Channel Observed Branching Fraction Limit [330,335,370] η → π 0 γγ (2.56 ± 0.22) × 10 −4 [30,365] Table 3: The observed branching fractions for radiative η ( ) decays and the assumed upper limit on contributions from a gauge boson of U (1) B , given at the 2σ level (98% C.L.). The limit for η → π 0 γγ is taken from the PDG average η decay measurements; we note that the measurement from CrystalBall@AGS [365] was the basis of the analyses in Refs. [330,370]. The observed branching fraction in the η → π + π − π 0 γ row corresponds to η → ωγ. boson of hypercharge in this context [384][385][386][387][388][389][390]. These interactions enhance the emission of longitudinal X in decays such as Z → Xγ [391][392][393]. Aside from this, there are also terms involving the charged W s, XW W ; these give rise to nonstandard flavor-changing neutral currents (FCNCs) such as b → sX at the quark level or B → KX at the hadron level [391,392]. Limits from these anomalous decays, however, depend on ε, so we have not shown them in Figs. 11-13. We further note that FCNCs can also appear in models with generationdependent couplings, such as the U (1) B 1 model discussed above: if the three left-handed quark doublets are charged differently under the interaction, then K − K or B − B mixing contributions at odds with experimental constraints are induced. Ref. [329] estimates that in the absence of new fermions, the couplings should satisfy g X |z is the charge of the first-generation (third-generation), left-handed quark doublet under the new interaction (in units of g X ). The only ways in which large couplings can exist without including additional fermions would be (1) to charge the third-generation quarks, thereby invoking the Υ(1S) constraint once again, or (2) to charge only the righthanded quarks, which are not constrained by FCNCs. However, in this latter scenario, the anomaly-cancellation conditions are sufficiently difficult to satisfy that it is all but required to invoke new fermions. The new states could, in principle, constitute some part of the dark matter; if so, then dark matter direct detection experiments predicated on nuclear recoils might encounter an observable scattering rate. We conclude this subsection by contrasting the scenario of gauged U (1) B symmetry with that of the related U (1) B−L . Both symmetries are anomalous in the SM, but U (1) B−L can be easily rendered nonanomalous by introducing three right-handed neutrinos. This property has made U (1) B−L an attractive candidate for a gauge symmetry of nature, and it has been studied extensively as a result; see, e.g., Refs. [380,394,395] and references therein. Because U (1) B−L models couple to charged leptons at tree level, this incurs strong constraints from processes involving electrons and positrons, to which U (1) B models are not subject (in the absence of kinetic mixing). However, constraints derived for U (1) B from, e.g., fifth force searches will also apply to U (1) B−L , albeit with some ∼ O(1 − 10) differences; see, e.g., Fig. 5 of Ref. [396]. To wit, the charge neutrality of matter implies that the numbers of protons and electrons must be equal; this implies that the B charge of some laboratory probe cannot be less than its B − L charge, though these are of the same order of magnitude.

Effects in Neutron Stars -Heavy X
If new, repulsive contributions are incorporated into the N N potential, then this stiffens the nuclear EoS -for a given (number) density of baryons, there is more energy density and pressure in a given fluid element than if these were not present. Conversely, attractive contributions soften the EoS. Stiffening the EoS has two primary effects, for our purposes here: 1. Neutron stars can be more massive. The increase in the energy density of a given fluid element, relative to some nominal prediction, partially offsets some of the gravitational binding energy of the neutron star, resulting in a heavier star. The heaviest confirmed neutron star, PSR J0740+6620, was initially determined to have a mass 2.14 +0.10 −0.09 M [397], though this has since been refined to 2.08 ± 0.07 M [398]; any candidate EoS must be stiff enough to support neutron stars at least this heavy. 2. Neutron stars have larger radii. As the EoS is stiffened, the increase in pressure makes nuclear matter harder to compress. Therefore, a neutron star of a fixed gravitational mass will be physically larger with a stiffer EoS.
Measurements of this sort have been performed for PSR J0030+0451 [407,408] and PSR J0740+6620 [409,410]; these indicate that the EoS is relatively stiff, supporting radii in the range 12-13 km over a wide range of possible masses (see, for instance, Fig. 11 of Ref. [409]). In Fig. 14, we show the neutron star mass-radius relationship for the Akmal-Pandharipande-Ravenhall (APR) EoS [411] -a representative EoS with nucleonic degrees of freedom based on the Argonne v 18 two-body potential [352] and the Urbana IX three-body potential [358], including relativistic effects -in blue. In orange, we have added a new vector interaction with coupling strength g 2 X /4π = 1 and mass m X = 600 MeV. For each curve, the thin, Figure 14: The neutron star mass-radius diagram. The blue line shows the predicted mass-radius relationship for the APR EoS [411]; the orange line adds to this a new, repulsive interaction with m X = 600 MeV and g 2 X = 4π. The thin, dashed portion of either curve represents the points for which the sound speed in the core of the star exceeds c; these points are unphysical. The pink and cyan regions, respectively, represent inferences for J0030+0451 [408] and J0740+6620 [410] from NICER and XMM-Newton; the dark (light) shading corresponds to 68% (95%) C.R. The red violin plot represents the posterior on the radius of a 1.4M neutron star, calculated in Ref. [6]. The black curve represents the mass-radius relation for black holes (i.e., the Schwarzschild radius), while the gray curve represents a constraint from causality [412].
dashed region denotes that the sound speed in the core of the star is greater than c; this is unphysical, and reflects that this EoS should not be applied for large central densities (n ∼ 5 − 10 × n sat , where n sat ≈ 0.16 fm −3 is the empirical nuclear saturation density). As described above, the mass-radius curve is shifted to larger radii, as reflected by the horizontal arrows, and the maximum neutron star mass is increased, as reflected by the vertical arrow. This figure also shows the 68% (dark) and 95% (light) credible regions (C.R.) for J0030+0451 [408] and J0470+6620 [410] in cyan and pink shading, respectively. Moreover, the dark red violin plot depicts the posterior probability distribution [6] for the radius of a 1.4M neutron star, R 1.4 , conditioned on a combination of data ("d") from (1) heavy pulsar masses, (2) the binary neutron star gravitational wave events GW170817 [413] and GW190425 [414], and (c) the same NICER observations of J0030+0451 and J0740+6620. The red shading represents the 90% C.R. 19 Taken together, these observations indicate that the neutron star EoS is required to be relatively stiff compared to the landscape of possible equations of state. The APR EoS that we have shown as an example is certainly compatible with existing data, but as more data become available in the coming decades from gravitational wave observatories, the nuclear EoS may end up being remarkably stiff. New repulsive interactions are one possible route to this outcome.
One effect not represented in Fig. 14 is rotation: for a fixed baryonic mass, a rotating neutron star will have a larger gravitational mass than a nonrotating one [415][416][417]. For maximal uniform rotation, this effect may reach the level of ∼ 20% [417]. Nonuniform rotation can lead to even larger enhancements [418][419][420], but these configurations are not expected to be stable on long timescales. The fastest known NS spin is that of J1748-2446ad [421,422], determined to be 716 Hz; the mass of this object is unknown, so it is unclear how close this object is to maximal rotation, and thus how large the contribution of its rotation to its total mass are. However, Ref. [423] has calculated how the mass-radius relationship for neutron stars with this rotational frequency differs from that of nonrotating stars; see Fig. 10 of that reference. The effect of rotation is to enhance the (equatorial) radius of the star, with the largest changes occurring for lower masses (around ≈ 1M ). If the spin of the NS were unknown, then rotations can mimic the effect of an intrinsic stiffening of the EoS; this is a potentially important systematic in making extractions of the EoS.
It has long been suggested that new physics may be operative within neutron starsindeed, the appearance of a new, long-range force would be a realization of modified gravity [424,425]. Some of the earliest work on explicitly determining the effects of a new interaction on neutron star structure comes from Ref. [426]. There, the authors determined that the new boson would modify the structure of neutron stars if the new coupling g and boson mass M satisfied The arguments of this paper are, however, approximate. In particular, this estimated limit arises from a comparison with the strength of the omega-exchange potential, g 2 ω /M 2 ω ≈ 400 − 500 GeV −2 [427]. Moreover, this treatment does not account for in-medium effects, which break the naïve dependence on g 2 /M 2 . Still, this work reflects the intuition that adding in a repulsive new interaction allows for larger, heavier neutron stars. Later work would show that, as expected, the precise modifications depend on the treatment of the baseline EoS considered; see, e.g., Refs. [428][429][430][431][432].
Recently, Ref. [335] presented an analysis in which a new two-nucleon force was studied in pure neutron matter accounting for in-medium effects using Brueckner-Bethe-Goldstone theory (see, for instance, Refs. [433][434][435] for details). In particular, this new interaction was studied in conjunction with the Argonne v 18 potential [352]; three-body and relativistic effects were studied through comparisons with the APR EoS [411]. Within this framework, the shift in the maximum neutron star mass, ∆M TOV , was calculated as a function of the new boson mass m X and coupling g X ; the results are shown in Figs. 12 and 13 for the cases of gauged B and B 1 , respectively. Specifically, the red-pink-white band sweeps over contours of constant ∆M TOV , with the red side corresponding to ∆M TOV = 0.1M and the pink side to ∆M TOV = 0.5M . As is evident from the figures, if the new boson has a mass in the range O(10 2 − 10 3 ) MeV, then the new coupling must be fairly strong, g 2 X /4π ∼ O(0.1 − 10) to have a marked impact, though it can also be rather weaker than the strong force in this region. Nevertheless these findings are crudely compatible with the expectation shown in Eq. (64). Most interestingly, if these effects were operative at this level, then they would contribute significantly to the rare meson decays we have discussed, allowing for a decisive test of this scenario.
We reiterate a point made previously. As of this writing, all usable nucleonic potentials are necessarily phenomenological -our knowledge of the nuclear force is fundamentally reliant on data. This state of affairs is necessary to make progress, but it is technically insufficient to add in a new interaction without considering how this changes the underlying parameters of the description. A useful analogy to this situation is muon decay. If one introduces a new interaction that contributes to µ − → e − ν µ ν e , then one cannot constrain this interaction through modifications to the rate. 20 This is because muon decays are used to define the Fermi constant, G F -this is one of a handful of parameters in the Standard Model that one must simply measure. Instead, one should probe consistency between a set of observables that all depend on G F , even if these do not directly couple to the new interactions. Therefore, we should aspire to a similar scheme in the context of nucleon interactions: one must formally adjust the parameters of one's phenomenological prescription to accommodate the new interaction(s) and test for consistency among some subset of independent data. As such, it might be difficult for two-or many-nucleon interactions to definitively rule out the existence of new interactions -but they can tell us where in parameter space to look.
Another issue on which we remark is the so-called masquerade problem [436]. In its original framing, this reflects the observation that hybrid stars, comprised of both nucleonic and quark degrees of freedom, may have a mass-radius relationship that can be emulated by a pure-nucleonic EoS. More broadly, a number of new phenomena may combine with the strong nuclear interactions in such a way that it is not possible to distinguish among them on the basis of the mass-radius relationship (though it may be possible to break this degeneracy using the g-modes of the NS [437][438][439].) One of the conceptual advantages of introducing new interactions among nucleons, however, is that these can be directly probed in the laboratory. Indeed, we have seen that if new interactions are strong enough to meaningfully alter neutron star structure, then these can be emphatically probed with, e.g., rare meson decays. This is qualitatively different from, say, the presence of critical phenomena, such as a transition to hyperon matter or quark matter: these latter scenarios are extremely difficult to probe in the laboratory -if not altogether impossible. Thus the 20 We emphasize that we are assuming this new interaction is subdominant to the SM. falsifiability of the former scheme is very much an asset, irrespective of whether the new interactions are connected to gauged B or not.
We also remark that the presence of new interactions is not mutually exclusive with the appearance of other critical phenomena. Of particular interest is the so-called hyperon puzzle (see, e.g., Refs. [303][304][305]). In Sec. 5, we noted that for heavy neutron stars, the nucleon chemical potentials become large enough that it is energetically favorable to produce strange baryons (Λ, Σ, etc.). This significantly softens of the EoS to such an extent that, in the absence of new ingredients, EoSs with strange baryons can fail to support compact stars heavier than ∼ 2M . It has been suggested that the appearance of hyperons can be pushed to higher densities and that their contribution to the EoS can be stiffened through (1) hyperon-hyperon repulsion, or (2) three-body forces involving hyperons (or both), among other possibilities; see Sec. 3.1 of Ref. [305] for a review of the literature on these subjects. 21 A new vector interaction between quarks may provide this additional stiffness. Additionally, the presence of the new interactions may modify the onset of certain critical phenomena by modifying, e.g., the relative energetics of the nucleonic and quark phases. This issue requires an in-depth treatment that is, to our knowledge, currently lacking in the literature.

Effects in Neutron Stars -(Ultra-)Light X
If the new boson is light -roughly below the MeV scale -then it may be produced onshell in the neutron star. In this case, it may exist as a real state with finite energy within the neutron star, beyond simply contributing to the potential of nuclear matter. The effects of this new state on the long-term evolution of the neutron star depend on its coupling to the nucleons. If the couplings are too weak, then the new boson will be produced too infrequently to have a meaningful affect on the star, as, e.g., on its cooling. As the coupling strength is increased, the state may begin to overcool the star: enough bosons are produced to transport significant amounts of energy out of the core. This energy would then manifest as either • a flux of the new state, which one could hope to observe directly; or • a flux of SM particles, produced via, e.g., decays of X or conversion to photons in the star's magnetosphere [441][442][443].
If, however, the new state is too strongly coupled, then the new bosons are produced copiously but are trapped by the medium -their mean-free path is too short to transport energy out of the star. In this case, there would be some accumulation of the new state within the nuclear medium. While these may modify the EoS, couplings in this regime do not lead to anomalous cooling. As an illustration of these ideas in a similar context, we show a constraint on the existence of a new gauge boson of U (1) B derived from observations of the neutrino signal of SN1987A from Ref. [340] in orange in Fig. 11. Where the constraint is operative, the lower bound represents the constraint from overcooling via X emission: if too much energy had been radiated in the form of X, then the neutrino signal would have been substantially reduced. The upper bound comes from trapping of X: if X is sufficiently strongly coupled, then it is produced and reabsorbed in the explosion, so that it does not prevent neutrinos from transporting energy out of the supernova. Constraints on anomalous neutron star cooling from new bosons have been derived for NS1987A (the remnant of SN1987A) [341,342] and Cassiopeia A (Cas A) [342]; Ref. [444] presents constraints derived from SGR 0418+5729, Swift J1822.3-1606, and 1E 2259+586. The precise form of the constraints depends on the scenario considered; typically, these studies are framed in terms of either (1) a dark photon that mixes kinetically with the SM photon, or (2) a new gauge state for U (1) B−L . The latter is closer to the scenario of gauged B than the former, but we note that gauged U (1) B−L contains a tree-level coupling to the electron, which changes the constraint. Still, neutron stars are largely (though not entirely) insensitive to the differences between B and B − L, so we expect constraints on the latter to be representative of constraints on the former. The constraints from NS1987A [341] and Cas A [342] are shown in red and dark red, respectively, in Fig. 11; we have omitted the constraints from Ref. [444] because these depend on the kinetic mixing between the new gauge state and the SM photon, or are otherwise subdominant. Cas A is the stronger of the two, setting a limit of g 2 X /4π O(10 −26 ) for m X 100 keV. 22 These energy-loss arguments are quite general and are not restricted to any particular source of BNV, nor to any particular astrophysical environment (though the form of the constraint is, of course, model dependent). To compare how these constraints fit in to the broader landscape of searches for new bosons, we refer the reader to, e.g., Refs. [380,381,394,445].
If instead the new boson is ultra-light -far below the eV scale -then its interaction range may be of macroscopic size. In this mass regime, the constraints discussed above prevent the new interaction for making observable changes to the structure of the neutron star or to its cooling, but there is a mass regime in which the new state can mediate the interactions between astrophysical objects. Constraints on an ultra-light gauge boson are shown in Fig. 15. The dark blue line represents the same set of constrains on a fifth force shown in Fig. 11, adapted from Ref. [309]. To these, we have added the exclusion from the Eöt-Wash torsion balance experiment from Ref. [308] in purple. The vertical gray bands represent constraints from black hole superradiance [446,447] -ultra-light bosonic fields can condense in the region just outside of a black hole, thereby sapping it of its energy and angular momentum. Constraints from observations of astrophysical-scale black holes have been derived in Refs. [448][449][450][451]; see also Ref. [381]. We also note in passing that the conjecture that gravity should be the weakest force [452], here loosely interpreted as m X g X M Planck , is satisfied throughout the parameter space that we have shown.
We have also included constraints from probes of nonstandard compact object interactions, on which we now elaborate. Firstly, if such a state exists, then it can induce modifications to the gravitational waveform of two merging compact objects. As we have Figure 15: Constraints on an ultra-light new gauge boson. The dark blue line is the same set of fifth force constraints from Ref. [309] as in Fig. 11, continued down to lower masses. The purple curve is the constraint from the Eöt-Wash collaboration [308]. The gray regions are excluded from nonobservations of black hole spin-down via superradiance [448][449][450][451] (see also Ref. [381]). The dotted brown and pink curves represent the sensitivities of the third-generation gravitational wave observatories Einstein Telescope and Cosmic Explorer, respectively, to modifications to inspiral waveforms, adapted from Ref. [453]. The dark and light orange contours represent constraints from anomalous energy loss via X radiation from B1913+16 and J0737−3039A/B, respectively. The dot-dashed contours are projected sensitivities from LIGO (light green) and LISA (dark green) under the assumption that the gauge boson of U (1) B is the dark matter. discussed, a new vector state induces a repulsive interaction between nucleons. For ultralight mediators, the contributions of the component baryons add coherently; neutron stars, having B ∼ O(10 57 ), will experience new interactions on macroscopic scales. These changes to the forces between neutron stars will manifest as apparent deviations from GR, inducing shifts to the gravitational waveform. To wit, a new vector interaction induces a long-range repulsion between neutron stars, thereby drawing out the merger over longer times. These shifts can be observable with gravitational wave interferometry [453][454][455][456]. The dashed lines in Fig. 15 represent the projected sensitivities of third-generation gravitational wave observatories to non-Newtonian contributions to binary NS mergers, adapted from Ref. [453]; brown is for the Einstein Telescope and pink is for the Cosmic Explorer. This reference explicitly concerns scenarios in which the non-Newtonian contribution arises from forces between dark matter cores, but we have reinterpreted their results in the context of gauged U (1) B . While this mechanism would be an exquisite probe of dark forces, in the case of gauged baryon number, the parameter space of interest is already excluded at high significance, speaking to the exquisite sensitivity of the Eöt-Wash measurements.
Secondly, if the two compact objects have different ratios of B/M , then the orbiting stars constitute a time-varying dipole of baryon number. This results in radiation of X as long as its mass is below the orbital frequency of the dipole [457], and, as mentioned above, expediting the merger in addition to modifying the inspiral waveform. Interestingly, this effect is present even if the compact objects are not close to merging. In that case, this manifests as a nonstandard contribution toṖ b /P b . We note Refs. [458,459] as particularly interesting realizations of this idea: the authors have studied several such systems in order to constrain a boson associated with gauged U (1) Lµ−Lτ , as opposed to U (1) B . We have adapted the analyses performed in these works for the case of U (1) B for the Hulse-Taylor binary (B1913+16) and for J0737−3039A/B, previously discussed in Sec. 3; we briefly summarize our calculation below. The quantity of interest is the ratio of the time-averaged energy loss from X emission, Ė X , and that from GW emission, Ė GW [457][458][459]: where M 1,2 are the masses of the compact objects, B 1,2 are their respective total baryon numbers, is the eccentricity of the binary, and J n (x) is the n th -order Bessel function. These observables have been summarized in  the APR EoS [411]. We then compare against the ratio of the intrinsic and GR-predicted spin-down rates of the binaries,Ṗ GR ḃ using the values presented in Table 2. Results for B1913+16 and J0737−3039A/B are shown in dark orange and light orange, respectively, in Fig. 15; 23 we have not included J1713+0747 as a part of this study, since its poorly constrainedṖ GR b /Ṗ int b ratio does not lend itself to a competitive constraint. Similar to the projected sensitivities to modifications to the inspiral waveform, these regions are entirely excluded by fifth force searches -precision studies of binary spin-down are well suited for some classes of new mediators, but are seemingly overwhelmed by terrestrial experiments whenever tree-level couplings to npe matter are present.
We note an interesting, complementary capability of the global GW observation program. If the new vector is ultra-light, then it may constitute some or all of the dark matter [460][461][462][463], though its mass is required to be no smaller than O(10 −22 − 10 −21 ) eV [464][465][466][467] (see also Ref. [468] and references therein). Such a light dark matter candidate would coherently oscillate over long length scales, producing a nearly monochromatic, stochastic gravitational wave signature. Ref. [469] has studied the sensitivity of LIGO and LISA to such a signature; Fig. 15 reproduces their findings for two years of observation, specific to the case of gauged U (1) B , in light green and dark green, respectively. These curves are presented in dot-dashing to remind the reader that these sensitivities assume that the new state constitutes the entirety of the dark matter. We note analyses of LIGO O1 data performed in Ref. [470,471]; the resulting exclusions are between one and two orders of magnitude weaker than the projection shown here. Similar projections have been made for other future GW observatories [472,473] as well as pulsar timing arrays [474,475]; see also Sec. V of Ref. [476].
We conclude this section by synthesizing some of the possibilities we have discussed in Fig. 16, where we show a generic NS-NS and NS-BH merger. We have not yet commented on how X might participate in such a merger; we briefly sketch this here. The green sinusoids represent the emission of X, which may occur (1) as a result of the time-varying B dipole, either before or during the collision, or (2) through thermal processes operative in the hot remnant. If X is not too strongly coupled, then either of these effects might lead to increased cooling rate of the remnant, or may lead to novel electromagnetic signatures. For instance, Ref. [478] has recently presented a calculation of the rate of dark photon production and

NS-NS NS-BH
Metastable NS Accretion Disk NS BH Figure 16: Evolutionary steps in a NS-NS or a NS-BH merger, after Ref. [477], with the possibility of X emission, as indicated by the green squiggles. Ultra-light X emission can appear throughout the merger event and can modify the inspiral wave form. In some portions of parameter space, X can be trapped in the interior, leading to anomalous cooling or transport of SM particles. If X stiffens the EoS of nuclear matter, then intermediate, metastable super-massive or hyper-massive neutron stars could potentially be heavier and longer-lived. These various modifications can modify the likelihood of the merger event following a particular evolutionary path and can also modify multi-messenger signals of the merger event at later times.
decay after merger; this produces a flash of gamma radiation potentially as energetic as ∼ O(10 46 ) ergs within the first second. More broadly, the new state may alter the cooling of the remnant in an observable way, or modify the kilonova signature in the hours, days and weeks after the event [479]. Moreover, if X is in the regime in which it can modify the N N potential, then the additional stiffness can modify the dynamics of the merger itself, or modify the intermediate state(s) that can appear in the merger. If the remnant is a supermassive or hyper-massive neutron star, then additional stiffness in the EoS may allow for heavier remnants to persist for longer times before collapsing to either a neutron star or black hole. There remains a significant amount of uncovered territory in understanding how these new states can manifest in extreme astrophysical environments; the age of multi-messenger astronomy promises a new avenue by which to study light, weakly coupled physics.

Other imprints of DM physics on NSs and their mergers
In this article we have delved into the ways different manifestations of BNV can impact the structure and evolution of a NS. Yet there are still broader ways in which dark sector and dense matter physics can intersect, and we offer a brief survey of the sweep of the possibilities here. These should largely be available to dark sector particles that carry baryon number, but could be possible even if they do not. Here we distinguish between DM and hidden-sector force mediators, where we assume for simplicity that such mediators are not sufficiently long-lived to be a component of DM.
Generally, DM can be produced within the NS, or it can accumulate within the NS through capture onto the star [480]. We consider each scenario in turn. The DM production rate can be extremely slow, in which case chemical equilibrium is never attained, as studied in Ref. [32], or it can be fast enough so that DM reaches chemical equilibrium with baryonic matter. The models proposed for the neutron lifetime anomaly are of this latter class, and they are severely constrained by the existence of NSs.
In the case that DM is captured on the NS, the exothermic nature of the reaction leads to DM thermalization and with the possibility of DM annihilation as well, can lead to significant heating of the star [481]. In such stars, if DM annihilation does not occur, or if DM carries an internal quantum number such as baryon number, then a DM core can form, and the possibility that this can induce collapse to a black hole serves as a constraint on the DM model, as studied in the case of bosonic Antisymmetric Dark Matter (ADM) models [8,482,483]. If the DM is also self-interacting, then the collapse into a black hole can be avoided, and asymmetric dark stars can form [484]. In this later case, the capture and accretion of ordinary baryonic matter on these dark objects could presumably occur, forming an outer skin of baryonic matter.
If significant amounts of DM can be either produced or accumulated in the NS, then the structure of the NS can be modified. In the presence of a NS with a dark core, the star can be become more compact with a smaller M max , modifying the M − R relationship [32,104]. This scenario has also been discussed in the context of the possibility of mirror neutron dark cores [103,485,486]. Moreover, new observables can appear in NS mergers, such as additional peaks in the postmerger frequency spectrum [487], and the tidal deformability (Λ) can decrease (increase) in the case of dark cores [32] (halos [488]); see also Ref. [489]. Moreover, the DM content of a NS can impact its EoS, so that, e.g., an EoS that was ruled out from the GW observation of the upper bound on Λ could be revived if a small admixture, say 5%, of DM were present. On the other hand, an assessment of the minimal value of Λ, as in Λ min ≈ 400 from Ref. [490], can also act as a constraint on the EoS with a dark core admixture.
The existence of dark, or hidden, sector mediators can also impact the structure and evolution of a NS. The mediator can couple to quarks and modify the EoS, making it either stiffer, as we have described in Sec 6.2, or softer, with concomitant implications for the tidal deformability. The existence of the mediator can also impact the nature of the critical phenomena that may occur with increasing density. In addition to this, Ref. [491] has considered the possibility of a dark lepton condensate at the core of a NS, noting that this can modify the neutrino transport properties in and evolution of the star.

Summary
The advent of the gravitational era has opened the nearby cosmos to us in new ways. We have considered the mechanisms by which BNV can exist and how such effects can combine with the physics of hidden sectors, whose inner workings are presumably key to the resolution of the dark-matter problem, or not to realize new ways of probing this physics through the study of neutron stars. Motivated by the neutron lifetime anomaly, we have observed that it is possible for apparent BNV to appear at rates not very much slower than ∼ 1% of the neutron lifetime and that these possibilities can be constrained through astrometric measurements. We emphasize that there is vast difference between the "global" BNV limits we have set in Eq. (41) from measurements of the decaying orbital period of binaries with at least one neutron star, presuming our assumptions of Sec. 3 hold, and the far more stringent constraints that appear on single-nucleon or dinucleon decays through either direct or indirect searches. These disparate limits are compatible in that only a fraction of the "star" may be active in regards to BNV processes, as particular local densities may be required for them to occur. Nevertheless, we have concluded that BNV, both real and apparent, is somewhat slower than the effective weak scale within the environment of a neutron star, thus making our studies of the thermodynamics of neutron stars with BNV both viable and concrete. Consequently the observables whence we can realize new insights into BNV and dark sectors, which emerge from our study, are: • In the presence of BNV, the distribution of neutron star masses to be found through gravitational wave studies can be expected to change with lookback time. We note that over the local volume available to us with present and next-generation gravitational wave detectors, the population of stars available to form neutron stars should differ little, making shifts in the mass distribution of the ensemble sensitive to the possibility of BNV effects, albeit likely apparent ones.
• It is possible that BNV can produce unbearably light neutron stars, leading to explosions with detectable signatures in X-rays or soft gamma rays [120]. This may be difficult to realize, however, as common mechanisms of neutron star formation favor roughly O(1 M ) stars and as BNV may become inefficient, due to its possible density dependence, in the lightest mass neutron stars.
• We can hope to detect neutron stars of sub-solar mass. This may speak to new mechanisms for neutron star formation and possibly, too, to BNV.
• The possibility of compact objects that are bright in X-ray or neutrino emission may allow the detection of these objects individually, or more probably, through their additional contributions to the diffuse supernova neutrino background, which may soon be detected at Super-K [492]. This effect has also been suggested from considerations of binary-star interactions [493].
• BNV processes with final states involving photons, mesons, and charged leptons (e.g., n → e ∓ π ± ) can be expected to dump all of their energy back into the NS, raising the temperature of the NS to a potentially detectable level, given upcoming observational possibilities, both in X-ray and the optical [494]. Ground-based follow-up optical studies of targets of opportunity from gravitational wave observations, given their expected sensitivity [495], may also yield new surprises. Put more pithily, old neutron stars should be cold; if they are not, then this would really be quite a coup.
We look forward to future discoveries!