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

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


INTRODUCTION
The cosmic excess of baryons over antibaryons is well established [1], but the theoretical mechanism by which it is produced is not. The essential theoretical ingredients are thought to be known: baryon number violation (BNV), along with C and CP violation, must all be present in a non-equilibrium environment [2]. Thus BNV would seem to play an essential role, though in the Standard Model (SM) BNV is thought to occur appreciably only at extremely high temperature [3,4] -and the existence of BNV at low energies has as yet to be established. In this paper we continue our scrutiny of such effects through observations of neutron stars, which contain enormous reservoirs of baryons. In earlier work we identified sensitive limits on BNV through the interpretation of precise observations of energy loss in isolated neutron stars and in neutron-star binary systems [5]. These studies limit the baryonnumber-violating effects that occur across the entirety of a neutron star. In this sense they are macroscopic limits. In this paper we interpret these limits in a microscopic way, in that we develop a framework in which they can be translated to limits on the parameters of particular particle physics models that generate baryon-number-violating effects.
The particular models to which our studies are most sensitive are those in which baryons decay or otherwise transform to dark-sector fermions, of O(1 GeV) in mass, that carry baryon number. In such cases BNV becomes an apparent, rather than explicit, effect, because the dark-sector particles are unobserved, even if baryon number is not broken. Although the existence of dark matter is certainly established through astrometric observations, both its nature and origin continue to be open questions. It is possible that the origins of dark matter and of the cosmic baryon asymmetry are related, so that the loosely similar value of the cosmic baryon and dark matter energy densities today may follow from a single underlying model [6]. The possibility of baryons that connect to hidden-sector baryons of comparable mass figure in many such explanations. In this paper we constrain this possibility through the study of neutron and hyperon transitions to final states with dark baryons in the neutron star. To our knowledge, an in-depth, quantitative study of non-SM processes within dense nuclear matter has not previously been realized 1 , and its execution necessitates much care.
The existence of neutron stars of about 2 M ⊙ in mass speaks to central densities in excess of three times nuclear matter saturation density, so that in this paper we employ relativistic mean field theory in baryonic degrees of freedom for our dense matter description, as its accuracy should improve with increasing density -and thus it should work best at the core of the star. We note that a neutron star may become a hybrid star, i.e., one with a quarkbased core predicated by a finite-density quark-hadron phase transition, if it is sufficiently heavy, and this possibility can also be constructed within this framework [9]. Transitions to dark baryons could also occur within the quark-based core, though we will set aside this possibility in this paper -and revisit it only in offering an assessment of our uncertainties in our concluding summary.
The broader possibility of dark decays of the neutron has been noted in explanation [10,11] of the long-standing neutron lifetime anomaly [12], in which the lifetime inferred from counting surviving neutrons is significantly different from that inferred from counting the protons subsequent to ordinary neutron decay. Although the discrepancy may arise from experimental effects, the possibility that dark decays contribute to it in some measure is a continuing possibility [5]. In this paper we provide severe limits on the flavor structure of possible new-physics models with dark baryonic sectors, such as Refs. [10,[13][14][15][16], that arise from the interpretation of neutron-star energy loss constraints we developed in Ref. [5]. In this paper we flesh out the general assumptions of that earlier analysis and note how the specific models we consider can satisfy them.
Let us conclude our introduction with a brief outline of the body of our paper. In Sec. 2 we detail the models of baryon dark decays that we are able to constrain through our neutron star studies, and we note how they are distinct from models that we cannot. We also compute baryon dark decay rates in vacuum, for later reference, as well as dark baryon 1 Albeit studies of exotic light particle emission in dense matter, which possesses simplifying aspects, are of long standing [7] and continue to be investigated [8].
removal rates, because our analysis assumes that SM dynamics determine the response of the star in the presence of BNV. In Sec. 3 we consider macroscopic baryon number violation in neutron stars, revisiting our earlier work [5] and fleshing out constraints following from its assumptions in greater detail. In Sec. 4 we develop how to evaluate particle processes within dense matter, employing RMF, as developed in Refs. [17][18][19], to describe the neutron star medium in β-equilibrium [20,21]. In this context, uncertainties in our description of the dense medium are captured through variations in the equation of state (EoS). With these developments in hand, we evaluate particle processes within our framework for the dense nuclear medium of a neutron star in Sec. 5 and use our macroscopic limits on BNV from Sec. 3 to report limits on the parameters of the microscopic models we consider in Sec. 6.
In Sec. 7 we consider the implications of our results for models of dark-sector baryogenesis and dark matter, and we offer a summary and outlook in Sec. 8.

PARTICLE PHYSICS MODELS OF BARYON DARK DECAYS
The possibility of hadronic processes with dark-sector particles naturally emerges in models that explain both the origin of dark matter and the cosmic baryon asymmetry, particularly if the dark sector candidate carries a baryonic charge [13,14,22]. Although it has long been thought that dark matter could also be described as a relic asymmetry [23,24], in these models, rather, the two problems are solved simultaneously [6]. More recently, highly testable scenarios [25] have been developed [16,[26][27][28][29][30], and we probe their flavor structure through the studies of this paper -and in Sec. 7 we consider the implications of the constraints that we find. Since the dark-sector particles are presumably SM gauge singlets, they could be light in mass, potentially with masses comparable to that of the known hadrons, and yet have escaped experimental detection thus far.
Our current discussion is loosely inspired by models connected to explanations of the neutron lifetime anomaly [10,16,31], with neutrons decaying to a dark baryon with a photon or an e + e − pair. Models with similar content have been considered for broader purposes [15,[32][33][34], and alternative solutions have also been noted [35][36][37]. The dark channels in the various models would impact the determined bottle lifetime, with a mirror neutron model [35] serving as a rare exception. There, neutron-to-mirror-neutron conversion occurs in a strong magnetic field, impacting the ability to detect protons in the beam-lifetime experiment. This last possibility has been excluded as a complete explanation of the anomaly by a direct experimental search [38]. We note that models that would explain the anomaly through neutron disappearance or decay to dark-sector final states can also be constrained by the close empirical agreement of the neutron lifetime with its measured A decay correlation as interpreted in the SM [5,39,40]. This agreement limits the branching ratio on such exotic processes to [41]: Br(n → exotics) < 0.16% (95% one-sided C.L.) , (2.1) where we note that the neutron lifetime anomaly is roughly a 1% effect [10] 2 . Direct experimental limits on n → χγ [45] and n → χe + e − [46] decays also exist, removing ranges of parameter space as an explanation of the anomaly. We will be able to set much more severe limits through our studies, where we note the limit on Λ → χγ from SN1987 for reference [16]. We regard the neutron lifetime anomaly as a motivation for further investigation of baryon dark decays, with new limits constraining the manner in which the co-genesis of dark matter and the cosmic baryon asymmetry could possibly occur. We now turn to the development of models of dark baryon decays.
Following Ref. [33], we introduce a Dirac fermion χ with baryon number B = +1 which interacts with SM quarks via the generic form where i, j, k are generational indices, Q and q denote a left-handed quark doublet and a right-handed quark, respectively -and color and Lorentz indices are left implicit. Such interactions can generate both decay and scattering processes involving dark final states, which we consider closely in this paper. First, though, we address their flavor structure.
We could neglect this possibility altogether, dropping all subscript dependence, but simple, renormalizable models that produce Eq. (2.2), at energies below the mass scale of their new physics, show that strong flavor sensitivity can nevertheless exist. Turning to models with leptoquarks [10,15], we consider colored scalars S 1 andS 1 transforming as (3, 1, 1/3) 3 and (3, 1, −2/3), respectively, under the SM gauge groups and SM invariant scalar-fermion interactions. Non-trivial flavor structure follows from the choice of leptoquark in that S 1 2 The most precise measurement of the A correlation coefficient yields the ratio of the axial-vector to vector coupling constants |λ| = 1.27641(56) [42], but recent measurements of the a correlation do not completely fit this picture, yielding |λ| = 1.2677 (28) [43] and |λ| = 1.2796 (62) [44]. 3 This variant was first considered in Ref. [10].
can mediate both n → χγ and Λ → χγ decay at tree level, whereasS 1 can mediate Λ → χγ at tree level but to mediate n → χγ requires a one-loop process with W ± exchange as well [15]. Thus in this paper we strive to probe both n → χγ and Λ → χγ decay processes.
These models also readily generate proton decay [10,15,33], noting p → χπ + or p → χK + decay as examples, so that the possible range of χ masses is rather restricted as a result.
We note that the stability of the 9 Be nucleus [10], particularly stability against 9 Be → χαα decay [47], requires m χ ≥ 0.937993 GeV , (2.3) slightly in excess of the proton stability constraint m χ > m p −m e , and that atomic hydrogen is stable if m χ > m p +m e = 0.93878 GeV [34]. If either constraint were not satisfied, then the empirical limit on the pertinent lifetime would bound the parameters of the model. Within the SM both systems are absolutely stable, yet empirical tests of that, with a determined lifetime as an outcome, should be possible. We note H lifetime estimates, made finite through a model with a suitably light χ, are made in Ref. [34]. Moreover, the radiative decay H → νχγ, which is subdominant relative to H → νχ, can be probed through measurements at Borexino [34,48]. Similar expectations follow from violating Eq. (2.3) -and a concrete estimate of the 9 Be lifetime can be found in Ref. [49].
In what follows we ignore the possible chiral structure of the quark-χ couplings and simply consider [16] L ⊃ u i d j d k χ Λ 2 + h.c. (2.4) Since the quarks carry electric charge, we have, at the energy scales for which baryonic degrees of freedom are pertinent, L n =n i/ ∂ − m n + g n e 8m n σ αβ F αβ n +χ(i/ ∂ − m χ )χ + ε nχ (nχ +χn) , (2.5) noting g n = −3.826 is the g-factor of the neutron [50]. This form also holds for the Λ upon the replacement n → Λ, taking g Λ = −1.22. After redefining the fields to remove the mixing term in Eq. (2.5), then if ε ≪ m n − m χ , with m χ < m n , we have [10,15] L n→χγ = g n e 8m n ε nχ m n − m χχ σ αβ F αβ n , (2.6) though potentially this operator could also stem from a distinct higher-energy source. Generally, the interaction of Eq. (2.4) can also generate transitions to dark baryon states with mesons, such as the decays n → χ + meson or Λ → χ + meson. Ref. [16] uses chiral effective theory [51] to relate the possibilities. We eschew this path because chiral effective theory ceases to be valid if the density of the neutron star medium much exceeds that of nuclear matter saturation density. Since our particular purpose is to set limits on microscopic models given BNV limits determined from observables associated with the entire neutron star, we set aside the study of final states containing both dark and hadronic degrees of freedom.
They are distinct from the final states we do study, and cancellations cannot occur. We thus expect that including these additional decays with hadrons can only make our limits more severe, though the inclusion of hadronic channels would make our estimates less sure. At low energies, the magnetic interaction of Eq. (2.6), employed in Refs. [10,15], can be used to compute n → χγ or Λ → χγ. A pertinent Feynman diagram is illustrated in Fig. 1.
Denoting B as either n or Λ, the total decay rate for B → χγ is given by in agreement with Ref. [15]. The mixing parameter ε Bχ follows once the UV model is given, and it is what we constrain through the analysis of this paper.
To determine the impact of these microscopic processes on the neutron star requires further model building. Thus far, at low energies we have a dark baryon χ, which we take to be a massive Dirac fermion. If it is a stable particle, then it can also be a dark matter candidate. If so, then it may already exist within the material that collapsed to form the protoneutron star, though likely only in small amounts, and through dark decays or adsorption on the star it may accumulate within the star. If it is able to give up its kinetic energy, then it may settle in the core of the star, ultimately impacting its properties and evolution. There are many processes in which χ could participate, though the interactions with baryons are severely limited by the cold, degenerate nature of the interior of the neutron star. In principle, given the nχγ and nχπ 0 effective interactions in the models we have noted, and using N to denote either a neutron or a proton, χ could (i) be produced via nN → χN scattering, (ii) interact elastically with another nucleon via a nN intermediate state, (iii) be formed via the annihilation nn → χχ or (iv) it can decay via χ → p + e − +ν e if it is heavy enough. The reverses of the reactions in (i) and (iii) could also occur. Pauli-blocking effects associated with the cold, dense neutron medium strongly suppress all of the reactions in which nucleons appear in the final state. Moreover, χ − N elastic scattering is further suppressed in that it occurs at O(ε 2 nχ ) at amplitude level. We note Fig. 2 for an illustration. Given this, and our interest in limiting BNV within the star in a model-independent way, implying that the response of the star to BNV ought be controlled by SM dynamics, we think that ensuring χ disappearance is important. Thus we consider two different pathways to do just that. In the first, we add χ-lepton interactions [5], which intrinsically break baryon number and are intrinsically very poorly constrained. We would also want the rate for χ decay to be no less of that for χ production. This path, however, is potentially subject to severe constraints from proton decay experiments. For example, we could have χ → e + e − ν or χ → 3ν, and these channels could give rise to proton decay via an off-shell χ * state as in (2.8) (Exotic proton decays of just this ilk also emerge in models with quark and lepton compositeness [52].) Admittedly, this process, as well as the p → π + 3ν channel, may evade severe constraints due to the particular nature of existing |∆B| = 1 searches, both because of the final-states studied and the cuts on the final-state particle momenta needed to control backgrounds. Although this path could prove to be viable, we favor an alternate choice: we will allow χ to decay to other dark particles. A simple realization of this is given by [14] L dark ⊃ y dχ ϕ B ξ + h.c. , (2.9) where ϕ B is a complex scalar with B = +1 and ξ is a Majorana fermion -and both are dark-matter candidates. Introducing a Z 2 symmetry, so that χ, n, and p are all Z 2 even, but ϕ B and ξ are Z 2 odd, we see Eq. (2.9) is the only surviving interaction that traces to the visible sector, with n − ξ oscillations, say, forbidden by the Z 2 symmetry. One interesting consequence of this new path is that dark decays can be induced in the scattering with either ϕ B or ξ in the initial state, as developed in Ref. [53] and illustrated in Fig. 3. A similar mechanism, considered in the context of the neutron lifetime anomaly, has been studied in Ref. [36]. The same process can destabilize the proton, with |∆B| = 1 experimental studies constraining the model parameters [53]. We note that theχϕ B ξ interaction can also induce χχ annihilation, as noted and illustrated in decay from the star, yet ϕ B could potentially accumulate in its core -and impact the survival of the neutron star [55]. If we suppose, rather, that ϕ B is light enough to escape the star, then that outcome can be avoided. We now turn to the explicit evaluation of proceses that can remove χ from the neutron star.  Fig. 3. Alternatively, χχ annihilation via ϕ B exchange in tchannel would yield a ξξ final state, which could ultimately rematerialize as aχχ pair.

A. Dark Baryon Removal Rates
If the masses of ξ and ϕ B sum to less than the mass of χ, then the decay χ → ξϕ B is operative. Using Eq. (2.9) and Refs. [54,56], we calculate the width of this decay to be However, if this decay is operative and if m ϕ B + m ξ < m p − m π , then this allows for proton decay via p + → π + ξϕ B . We avoid potentially running afoul of these constraints by insisting that this decay not be operative and thus require m ξ > m χ . Instead, we focus on possible annihilation processes of χ, where we have assumed that only ϕ B is lighter than χ. Adopting the same tools to compute χχ → ϕ B ϕ B we have: We note that this cross section goes to zero as m ξ → 0. This must occur, so that this outcome serves as a non-trivial check of our procedure. Our cross section result does not depends on whether the scalar is real or complex, but its interpretation does. If the scalar is real, it cannot carry baryon-number, and χχ annihilation to scalars would then break B by two units. This can only occur if m ξ has a nonzero baryon-number-violating mass. Thus its rate vanishes if m ξ does.
We would like to understand how these annihilation processes operate within a neutron star. As we will see in Sec. 3 A, these cross sections would need to be averaged over the true distribution of χs produced in baryon decays within the star. Generically, χs need not be distributed thermally, and the process of thermalization would require self-interactions, which do not appear at tree level in our simple model. The problem of χ transport in the neutron star is beyond the scope of this paper, so that we assume that the thermally averaged cross section is a reasonable estimate of what the true averaged cross section would be.
We proceed by employing pertinent results from the seminal Ref. [57]. The thermally averaged cross section ⟨σv⟩ is given formally by where T χ is the χ temperature (which is generically nonzero and may be different from the temperature of the rest of the neutron star) and K 1,2 are modified Bessel functions of the second kind. This expression assumes that it is appropriate to describe the χ fluid as abiding by a Maxwell-Boltzmann distribution; it would be inappropriate to apply this expression to a cold, degenerate population of χ, but such a population does not occur in our framework.
To perform the thermal averaging, we expand σ(s) × v in powers of ϵ ≡ s/(4m 2 χ ) − 1: 2.13) this requires that v = 2 ϵ(1 + ϵ)/(1+2ϵ). In the limit in which the χ fluid is nonrelativistic, the thermally averaged cross section can we written in terms of the coefficients a (n) as follows: This prescription is expected to be valid as long as T χ ≲ 3m χ [57]. For χχ → ϕ B ϕ B , we find the leading-order contribution to the thermally averaged cross section in T χ to be Since the a (0) term vanishes, we conclude that the s-wave annihilation contribution vanishes, resulting in a suppression at low temperatures. We expect our χs to have a nonzero average kinetic energy from decays, so we do not expect to encounter a scenario in which these annihilations are completely quenched by the low energies of their parents, but it is an interesting feature to note.
We conclude by noting some relevant qualitative features of this model. Since χ selfinteractions do arise at the one-loop level as a result of interactions with ϕ B and ξ, we can expect the χ population would thermalize, but that timescale is likely slow relative to that of their annihilation to scalars. There are many more interesting phenomenological consequences of this model that one could explore, but for our purposes, it is enough to assume that the masses and coupling conspire such that χ can be removed from neutron stars quickly enough that our formalism is valid.

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

A. Assumptions
The structure of a neutron star can be approximated by a static and spherically symmetric metric (g µν ) with a line element given by [58] in which ν(r), λ(r) are solutions to the Einstein field equations [59], G µν = −8πGT µν , in which G µν is Einstein's tensor, G is the gravitational constant, and T µν is the stress-energy tensor. The rotation effects on the neutron star structure, which are O (Ω 2 / (G M/R 3 )) [60], amount to less than 3% for the fastest spinning pulsar (J1614−2230) that we consider in this work. Furthermore, the inclusion of quasi-static BNV processes, which are sourced by the matter in the star, would keep the spherical symmetry intact and changes to the metric (g µν ) very slow in time, such that the use of Eq. (3.1) is warranted.
We also assume that the medium in neutron star can be described by a perfect fluid with as the only nonzero components of the stress-energy tensor in which P and E are the local pressure and energy density of the fluid respectively which in general depend on the local baryon number density (n) and temperature (T ) via the EoS. In the standard picture, neutron stars cool down to internal temperatures T ≲ 10 11 K ≪ E F ≲ GeV within a minute after formation [61], such that the thermal contribution to the pressure and energy density can be neglected. The neutron star fluid can then be described as a cold degenerate Fermi gas at β-equilibrium. The existing terrestrial constraints on neutron dark decay, Eq. (2.1), along with the BNV limits we find in Table I, show that BNV rates should be slow with respect to other dynamical processes in the neutron star. We have also devised a model in which χ, the dark baryon-like particle, can be removed efficiently from the star. Thus we expect the deviations from a degenerate state at β-equilibrium due to BNV should be negligibly small, and we leave a more detailed study of possible thermal effects on neutron stars from BNV to future work.
In order to be able to apply our model-independent formalism [5], we are going to focus on a subset of models in which the dark contributions to the EoS are negligible relative to the energy density and pressure of the visible sector. In other words, we demand that the following (local) conditions hold throughout the neutron star at all times, which can be equivalently written as a condition on the local number density of χ: n χ (r) ≪ n(r). This means that χ has to decay or annihilate either back to the visible sector or to some other dark particles that can escape the neutron star. We assume that χ participates in self-annihilation to lighter dark particles that can escape the neutron star (see Sec. 2 A for more details).
We can express the condition n χ (r) ≪ n(r) in terms of the BNV rate, Γ BNV , and the annihilation cross section that is averaged over χ distribution, which we denote by ⟨σv⟩. We note that the exact distribution of χ in the neutron star can in principle be found by solving the Boltzmann transport equation in the star, but this is not practical for our estimation purposes. We instead consider two scenarios for χ: one in which the annihilation rate is much faster than the self-interactions which help establish a thermal equilibrium, and another in which self-interactions of χ are much faster than its annihilation rate.
We first consider the scenario in which dark particles have a non-thermal distribution at the time of their annihilation. If we ignore the effects due to radial redistribution of χ after their production and prior to their annihilation, their number density (n χ ) would approximately satisfyṅ (3.4) in which n i (t) is the decaying baryon number density which we take to be constant on short timescales, and ⟨σv⟩ is the annihilation cross section averaged over the distribution of χ.
The asymptotic value for χ number density (at times t ≫ 1/ n i Γ BNV ⟨σv⟩) is then equal to n ∞ χ = n i Γ BNV /⟨σv⟩, which relative to the local baryon number density n(r) is given by in which f i (r) ≡ n i (r)/n(r) < 1 is the fraction of baryon i relative to the total baryon number density, n sat = 0.15 fm −3 is the nuclear saturation density, and we used the scale of the canonical weak-scale cross section (10 −26 cm 3 s −1 ) for comparison. We can see that this ratio is negligible for the reference values in this equation if ⟨σv⟩ ≫ 10 −56 cm 3 s −1 .
We can generalize Eq. (3.4) to scenarios in which the redistribution of χ's , after their production and prior to their annihilation, is not negligible, by noting that the total χpopulation satisfiesṄ (3.6) in which B i (t) is the number of decaying baryons of type i, and C ann is the annihilation rate per particle, such that the total annihilation rate is identified as Γ ann ≡ C ann N 2 χ /2. We are interested in short timescales during which B i (t) can be taken as a constant (t ≪ Γ −1 BNV ). In this case, the solutions to Eq. (3.6), assuming N χ (0) = 0, are given by in which the timescale for achieving an equilibrium between the production and annihilation of χ (Ṅ χ (τ ∞ ) ≈ 0) can be identified as τ ∞ = 1/ √ B i Γ BNV C ann , which can be achieved for The total number of χs can then be approximated by its equilibrium value given by N ∞ χ = B i Γ BNV /C ann . We can see that if the condition in Eq. (3.8) holds, then N ∞ χ ≪ B i . We now calculate C ann in the scenario in which the annihilation rate of χ is slower than its self-interaction rate, and the χ's are distributed spherically with an average radius of R χ , according to Boltzmann distribution. Using the virial theorem and assuming a radially uniform distribution of background neutron star matter (over R χ ) with an average energy densityĒ, we can write 9) in which k B is the Boltzmann constant, and T χ is the dark sector temperature. The total annihilation rate (Γ ann ) and N χ can then be evaluated as (3.11) in which ⟨σv⟩ is the thermally averaged annihilation cross section. Using the definition of C ann we have 12) and an equilibrium between the production and annihilation can be achieved on timescales in which T χ and m χ have the same units. We can also find the equilibrium value for χ number density at the core by combining the definition of equilibrium number, N ∞ χ = B i Γ BNV /C ann , with Eq. The baryon decay rate (per baryon) in a small volume (V ) in the nuclear matter (n.m.) rest frame (Γ nm ) is defined by d(nV )/dτ = −Γ nm n V , in which τ is the fluid's proper time, and n is the proper baryon number density. We can define a baryon number-flux vector by j µ = u µ n [64], in which u µ is the four-velocity of the fluid (u µ u µ = 1) and use the definition of Γ nm to write j µ ;µ = −n Γ nm , in which ';' denotes the covariant derivative. We then use the relationship √ −g j µ ;µ = ( √ −g j µ ) ,µ [20], in which ',' denotes ordinary partial derivative and g ≡ det|g µν |, to arrive at where B is the total baryon number of the neutron star. We have used √ −g = exp(ν(r) + λ(r)) r 2 sin θ, with exp(2λ(r)) = (1 − 2M (r)/r) −1 , and M (r) is the total mass included within radius r: Given a particle physics model for BNV we can evaluate Γ nm (r) and use Eq. (3.16) to find the resultingḂ.

B. Framework
It was shown in Ref. [5] that the conditions in Sec Here, we ignore any possible dependence of O on the angular velocity Ω, i.e., we assume O evolves along a one-dimensional trajectory with Ω = 0 on the general two-dimensional space parameterized by E c and Ω. We can solve forĖ in terms of the rate of baryon loss,Ḃ, such in which we defined the effective BNV rate Γ BNV ≡ −Ḃ/B and the dimensionless parameter b(O) encodes the relative rate of change in O with respect to Γ BNV . We pick hadronic versions of the DS(CMF) EoS [21] that includes a crust [65] from the CompOSE database [66]. The details of these EoS including their Lagrangians and particle contents are given in

C. Observables
Baryon loss in pulsars may lead to observable effects on their individual spin-down rate (Ṗ s ), and their orbital period lengthening (Ṗ b ) if they belong to a binary system [5]. The BNV modifications toṖ s are caused by the quasi-equilibrium changes in the moment of inertia (I), and angular momentum loss due to light particles (e.g., ϕ B ) escaping the pulsar.
While the first contribution can be expressed in a model-independent manner, the latter depends on the specific BNV model and the masses of particles involved. Therefore, we focus our attention on BNV modifications toṖ b , which can still be formulated in a model independent way.
The energy loss due to BNV can modify the orbital period decay rate in a binary system, assuming it is active in one or both of the components. This energy loss can be written in which b(M ) and b(I) are defined in Eq. (3.19), P s andṖ s are the observed pulsar spin period and its observed rate of change respectively. Note that the rates of change in I due to spin-down, (dI/dΩ)Ω, are negligible in the pulsars that we consider. The relative rate of change in a binary period due to energy loss in its components is given by [67][68][69] in which 1 and 2 refer to the components of the binary system. After plugging Eq. (3.20) into (3.21), we get the following BNV and spin-down contributions to the energy-loss term We should note that the second term in Eq.

D. Interpretation
The dominant contributions to the observed relative rate of orbital period decay can be written as [70]: in which the first term is due to gravitational radiation [71], and the third term includes extrinsic effects, e.g., due to the relative motion of a binary pulsar with respect to the solar system barycenter. The numerical values for each of these contributions and the limits oṅ PĖ b (found by subtracting the GR contribution,Ṗ GR b , from the intrinsic orbital-period decay rate,Ṗ int b ≡Ṗ obs b −Ṗ ext b ) are given in Table I for three binary systems. Two of these systems J0348+0432 and J1614−2230) have heavy pulsars that may contain hyperons [72], and the third one is a double pulsar system (J0737−3039A/B) with an extremely high precision in its orbital parameters.
1. PSR J0348+0432: A pulsar-white dwarf binary discovered in 2007 with the Robert C. Byrd Green Bank Telescope [73] with an orbital period of about 2.4 hr. We use the results from the analysis in Ref. [74], in which it was shown that the kinematic, spin-down (Eq. (3.23)), and tidal (Ṗ T b ⪅ 10 −16 ) contributions toṖ b are negligible and the observedṖ b should be mainly caused by the GW emission. We use the value from Ref. [74] for the intrinsic period decay rate,Ṗ int b = −0.275(45) × 10 −12 .
2. PSR J1614−2230: A pulsar-white dwarf binary discovered in 2006 with the Parkes radio telescope [75]. We use the Shapiro delay mass estimates from Ref. [76], and the binary parameters from NANOGrav 12.5 yr data set [77] at 56323 MJD. The observed value ofṖ obs b = 1.57(13) × 10 −12 is dominated by the Doppler shift due to the pulsar motion which is itself mainly caused by the Shklovskii effect [78]: in which we input the value for proper motion µ = 32.4(5) mas yr −1 , and used the parallax distance d = 0.65 ± 0.04 kpc [79]. We use Eq. (16) from Ref. [80] to estimate the contribution due to the Galactic potential, namely, for the period derivative,Ṗ int b = 0.32(16)×10 −12 , is positive at 2σ significance, pointing to a possible underestimation of extrinsic effects and their errors. However, we note that if, for example, we instead assume a negligible value forṖ int b ≈ 0 and double our error estimates, then we would still obtain the same limits. We also evaluate the relatively small GW contribution which for circular orbits is given by [71] in which we used the pulsar and white dwarf masses from Ref. [76], T ⊙ = 4.92549094× 10 −6 s, and we neglected the small eccentricity of the orbit e = 1.333(8) × 10 −6 [81].
In estimatingṖΩ b using Eq. (3.23) we assumed the canonical value I = 10 45 g cm 2 for the pulsar's moment of inertia.
3. PSR J0737−3039A/B: A double pulsar discovered in 2003 [82], comprised of two radio pulsars (A and B) with pulse periods of 22.7 ms and 2.8 ms, respectively. We use the data from Ref. [83] and the inferred limits on BNV contributions from Ref. [5].  We can now translate the bounds on (Ṗ b /P b ) BNV from Table I to limits on (Ḃ/B) using Eq. (3.22), which are presented in the last row of Table I. In deriving these limits, we assumed that BNV is only active in the pulsars. We also note that we can only infer a model-independent limit on a linear combination of BNV in pulsars A and B of the double pulsar system (J0737−3039A/B). However, we expect that the rates of BNV (per baryon) would be about the same in both pulsars, i.e., ( In Sec. 6, in which we adopt a specific BNV model (B → χγ), our inferred limits on the mixing parameter (ε Bχ ) are found by evaluating the individual BNV rates in each of the two pulsars J0737−3039A and J0737−3039B, which we then sum to compare to the observational limit on BNV in this system. We also observe that changing between the DS(CMF) EoSs (see Table II) induces variation in, at most, the last significant digit in our limits (see the discussion below Eq. (3.23)).

DENSE MATTER CONSIDERATIONS FOR PARTICLE PROCESSES
Different lines of evidence reveal that dense matter environments can be discriminating probes of non-SM processes. For example, limits on Λ → χγ, as well as other decay channels with dark particles, follow by noting that the duration of the observed neutrino pulse in SN 1987A should not be significantly impacted by dark sector emission [16]. We, too, have found severe limits on BNV from binary pulsar period lengthening, as shown in Table I. Here we sharpen such studies by computing particle processes within a theoretical framework suitable to the description of the dense matter in the interior of a neutron star.
To compute particle processes in dense matter we might first turn to chiral effective theory to describe the low-energy interactions of such hadrons [84,85]. At the simplest level, these studies exploit the symmetries of QCD to systematize the interactions of mesons and baryons in a momentum expansion in powers of (Q/Λ χ ), in which Q is the momentum or pion mass and Λ χ is the chiral-symmetry breaking scale (Λ χ ≈ 1 GeV), with experiments fixing the value of the unknown low-energy constants (LECs) that appear. This framework can also be extended to the determination of the EoS of neutron stars [86,87]. The empirical nature of the LEC determinations limit the applicability of chiral effective theory to densities no more than 2n sat [88]. Moreover, in neutron stars, the central densities can easily exceed that of saturation density by a factor of a few, making the nucleons relativistic. As a result, we turn to relativistic mean-field (RMF) theory in hadronic degrees of freedom to describe the dense matter at the core of a neutron star. In what follows we first describe how a RMF treatment emerges from a simple, covariant quantum field theory description of hadronic interactions before describing the specific chiral mean-field (CMF) EoS that we employ for generating our numerical results, showing how this specific choice maps onto the RMF treatment of the simpler model. We then show how particle decays can be computed within that framework.

A. Modelling Dense Matter
A prototypical choice is the Walecka model [17][18][19], namely, where F µν = ∂ µ V ν − ∂ ν V µ and a counterterm δL, as the model is renormalizable. It is similar to massive QED with a scalar extension and a conserved current (baryon number).
Both a neutral scalar meson (φ) and a neutral vector meson (V µ ), describing the attractive and repulsive features, respectively, of the nucleon-nucleon force appear. The equations of motions (EoMs) take the form The EoMs are nonlinear and thus complicated. Working in the mean field limit is grossly simplifying, however. That is, at high baryon number densities, the sources for φ(x) and V µ (x) fields become large, and these field operators can be replaced by their vacuum expec- In doing this, we assume rotational invariance and note that in static uniform matter, as in a neutron star, φ and V 0 become constants that only depend on density. The solutions to Eq. (4.4) would take the form of that of the free Dirac equation if the replacements In other words, the medium effects in the RMF limit are captured by a shift in the baryon momenta and masses. In generalizing this result for broader use, we note that the Lagrangian of interactions for a more realistic hadronic model would have more ingredients (e.g., mesons). However, we would still be able to add up the scalar meson VEVs that modify the baryon's mass in a similar manner and denote the effective baryon mass by m * , independent of the specific scalar mesons in our model. Similarly, we can combine all the contributions to the baryon's momentum from vector mesons and denote them by Σ µ , such that in going from the vacuum to the in-medium formalism we would replace k µ → k * µ ≡ k µ − Σ µ . Equipped with this result, we can write the wave-function for a baryon with (canonical) four-momentum k µ (in a uniform medium) as Σ is defined to be the kinetic four-momentum and the vector self-energy (Σ µ ) is generated by the vector meson VEVs, with ⃗ Σ = 0 in the n.m. frame. The time-component of k * µ is defined by E * (k * ) ≡ m * 2 + | ⃗ k * | 2 , in which m * is generated by the scalar meson VEVs. The baryon spinor u(k * , λ) satisfies the Dirac which has the following solution in Dirac-Pauli representation in which ⃗ σ contains the Pauli matrices, and χ λ is the Pauli spinor with χ ↑ = (1, 0) T and Note that u has a Lorentz-invariant normalization given by u(k * , λ)u(k * , λ) = 2m * . The wave-function for antibaryons can be similarly constructed. The energy spectrum of baryons (k 0 ) is given by in the mean-field approximation, Σ µ and m * do not depend on k µ but they do vary with density. The values for m * and Σ 0 (in the n.m. frame) decrease and increase respectively (see Fig. 6) in such a way that the total energy of baryons in Eq. (4.8) increases at higher densities. As we will see shortly, this brings about in-medium baryon decays to particles that are heavier than the baryon's vacuum mass since E(0) > m B at high densities. In general, the increase in the repulsion between baryons in a RMF framework can be understood by comparing the time-like component of vector (repulsive) interactions, which are proportional to u † u, with scalar (attractive) interactions, which are parameterized by uu = (m * /E * )u † u.
As the density increases, m * decreases and the strength of the attractive forces relative to the repulsive ones diminishes [89]. However, we should note that having a highly repulsive nuclear interaction at extremely high densities (compared to n sat ) is a reasonable expectation, regardless of the specific dense matter formalism. Having explained the formalism utilized in this work, we now describe the specific EoS that we use for generating our numerical results.
We choose an EoS based on a non-linear hadronic SU(3) CMF model [90], in which the baryonic degrees of freedom include nucleons (n, p), hyperons (Λ, Σ, Ξ) and the spin-3/2 resonances (∆, Σ * , Ξ * , Ω). These baryons interact via exchange of scalar (σ, δ, ζ, χ) and vector mesons (ρ µ , ω µ , ϕ µ ), in which ρ µ and δ are both isovectors. In the RMF limit, the mesons become classical fields, and in the n.m. frame only the zeroth components of vector mesons develop VEV. The Lagrangian density of the CMF model is given by [21] L = L Kin + L Int + L Self + L SB , (4.9) in which L Kin contains the usual kinetic terms for baryons and leptons, L Int is due to the baryon-meson interactions which are given by We note ψ i denotes a baryon of species i with an effective mass m * i and an isospin 3component I 3i , and the expectation value is evaluated in the ground state. The last two terms in Eq. (4.9), i.e., L Self and L SB , contain the self-interactions of scalar and vector mesons and explicit chiral symmetry breaking terms respectively. The explicit expressions are given in Eqs. (3), (4), and (5) of Ref. [21]. The baryon effective masses are generated by the scalar meson VEVs, except for a small explicit mass term δm i ∼ 150 MeV, and are given by The time-component of baryon self-energy is given by The numerical values for m * and Σ 0 are plotted in Fig. 6. We note that the reduction of the effective baryon masses at high densities as shown in Fig. 6 is due to chiral symmetry restoration at high densities. The coupling constants are chosen [91][92][93]  This conventional approach in determining the coupling constants in RMF models relies on an extrapolation from symmetric finite nuclei to infinite neutron matter. We would like to contrast this with an alternative that we may wish to employ in the future, which is based on fitting uniform pure neutron matter properties determined through the use of chiral effective field theory [94]. The latter procedure involves fitting the RMF couplings with the synthetic neutron matter data generated using Quantum Monte Carlo (QMC) many-body methods [95], in addition to reproducing n sat , B/A, and K. , and delta resonances (∆); and the additional vector interactions ("Add. Int.") beyond the standard terms (L Self ) that are included for each EoS respectively. The fourth column represents the assumed value for symmetry energy (E sym ) slope (L). The fifth to eighth columns are the single-particle hyperon potentials, and the last column is the maximum neutron star mass (M max ) that can be generated.
Our chosen class of EoS has a set of variations that depend on the degrees of freedom that are included, and they are given in Table II  In this section, we discuss some of the notable features that emerge in studying processes in the medium, and make comparisons with the vacuum formalism. We start with the quantization of baryon fields in the medium followed by the rate and cross section calculation formalism. We then discuss the electromagnetic form factors of the baryons that are needed The presence of the baryon Fermi sea modifies the quantization procedure of the baryon fields, ψ(x), in medium [18] compared with the usual procedure in vacuum [99]. Once the coefficients behind Fourier modes of ψ(x) are promoted to baryon creation (a † (k)) and annihilation (a(k)) operators (likewise b † (k) and b(k) for antibaryons), we conclude that the action of these operators on the medium ground state |Ω⟩, which contains baryon levels filled to a Fermi momentum (k F ), should be given by This leads to a different form (compared to vacuum) for the baryon propagator which is given by [18] in which θ is the Heaviside step function, B µ is the baryon current density, which in the n.m. frame is given by B µ nm = δ µ0 n B , and the second term in Eq. (4.14) allows for the propagation of holes in the Fermi sea. Using this modified propagator and the spinors in Eq. (4.7), one can derive Feynman rules [18] for calculating the amplitudes for various processes (see Sec. 5 B). However, in calculating rates via phase space integrals, we should first observe that an on-shell (p * 2 = m * 2 ) and positive energy (p 0 > 0) Lorentz-invariant integral over the four-momentum is given by Therefore, we identify the Lorentz-invariant (on-shell) volume element in the medium as . This means that the normalization factors in the in-medium phase space integrals should contain (2E * ) −1 in place of the usual vacuum expression.
We also note that the velocity of a baryon is defined in terms of the kinetic momentum as opposed to the canonical one, i.e., v µ ≡ k * µ /E * . This velocity should be used for calculating the cross section of two-body scattering involving a baryon (see App. C). We can explicitly show this by performing an integration over the longitudinal (ẑ) components of the incident beams' momenta (k z A andk z B ). Let us assume for the moment that only one baryon (B) is involved, in which case we have (see Eq. (4.77) of Ref. [100]) 16) in which in the last line we are assumingk z B = p z f −k z A and have identified the baryon velocity using the kinetic momentum, such that |v A − v B | is the relative velocity of the beams as viewed from the laboratory frame. The generalization to the case with two baryons is straightforward. The fact that the velocity of a baryon is zero when ⃗ k * = 0 could have also been deduced by inspecting the kinetic energy component in Eq. (4.8). For this reason, the frame in which ⃗ k * = 0 holds is called the center of velocity (c.v.) frame which is distinct from the center of mass (c.m.) frame defined by ⃗ k = 0. Therefore, the decay rate of a baryon in an arbitrary frame (Γ) is found by boosting (γ) the rate evaluated in the c.v. frame using Since we study processes that involve electromagnetic (EM) interactions with baryons, the generalization of EM form factors from the vacuum to within the medium should be checked. The in-medium spinors in Eq. (4.7) are different from their vacuum counterparts. Therefore, certain commonly used properties (e.g., Gordon decomposition) in vacuum need to be reestablished. However, we note that the general form of these interactions is determined by the structure of Dirac algebra. While important for formulating our analyses, this is slightly tangential to the broader narrative of this work; we thus relegate the details to the Appendices, but we encourage the reader to study them nonetheless. In App. A, we explicitly show that the vacuum EM vertex form can be generalized to its in-medium form if one replaces m → m * , p → p * , and identifies the electric charge and magnetic moment of a baryon from the scattering amplitudes in the c.v. frame. Our numerical results in Sec. 5 assume the vacuum values for the in-medium form factors F * 1,2 of neutron and Λ. We also derive the non-relativistic limit of baryon's EM interactions and their elastic scattering formalism in App. B. We present the calculations for in-medium Compton scattering in App. C, as a demonstration of the RMF formalism utilized in this work.

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

A. General Considerations
To illuminate the essential points, we consider the possibility of n-χ mixing in a background field Σ µ , the vector self-energy of a neutron in the neutron-star medium, which interacts with the neutron field ψ n , but not the χ field ψ χ . Thus we adopt the following simple model: Under a field redefinition, ψ → ψ ′ , prescribed by Eq. (5.1) becomes If Σ µ were absent, and with ε real, then for tan (2θ) = 2ε/(m * n − m χ ), L ′ describes two decoupled fields with a modified energy spectrum. These fields can then map to the asymptotic ("in" and "out") states needed to define the S-matrix [102]. To do this, any interactions with these fields should vanish as t → ±∞. For the neutron (and other SM baryons), we note that the effect of the vector self-energy can be absorbed into the definition of a modified single-particle spinor, as discussed in Sec. 4 A, and thus suitable "in" and "out" states can still be constructed. In the current case, Σ µ mediates an interaction between the rotated n and χ fields, putting the utility of our field redefinition procedure into question. After all, even in the mean-field limit, Σ 0 can greatly exceed the n and χ masses at the high densities reached within a neutron star, and it cannot vanish as t → ±∞, since we work within a medium of infinite extent. Since Σ µ is not a Lorentz scalar, we cannot extend our field redefinition approach to include it. Therefore, there would seem to be no advantage to following a field redefinition approach in neutron matter. Moreover, in the small mixing limit (ε ≪ |m * n − m χ |), the mass (n ′ , χ ′ ) and interaction (n, χ) eigenstates are nearly the same. Working with Eq. (5.1), we can treat εψ n ψ χ as a tiny interaction that mediates n ↔ χ transitions within perturbation theory. This Feynman diagram analysis, through the in-medium baryon propagator, Eq. (4.14), naturally includes the impact of momentum dependence and of the neutron self-energy on n-χ mixing. We emphasize that both effects are absent in the field redefinition procedure. As a result, too, we do not have large enhancements in our predictions should the in-medium neutron and χ states become degenerate in energy -the imaginary part of the neutron self-energy effectively eliminates that possibility.
Nevertheless, n-χ mixing within the neutron-star medium could potentially lead to effects not possible in terrestrial experiments, and we consider those possibilities more carefully in Sec. 5 C.

B. Dark Decay Rate Estimates
We now turn to the explicit evaluation of rates of particle processes within the neutronstar medium, with a particular focus on dark decay rates. As long known, the background field associated with matter leads to a spontaneously breaking of Lorentz symmetry, but as a consequence of our Lorentz covariant description, discussed in Sec. 4 A, our expressions always have definite Lorentz transformation properties. In what follows, we exploit our freedom to choose a frame to simplify our analysis.
Generally, processes of the form B+{X} → χ+{Y } lead to the following rate of change of the local baryon density n B (with respect to the proper time, τ , referenced to that spacetime point): is some set of other states in the initial (final) state -which may be empty.
are the species-dependent occupation numbers 4 , and |M| 2 is the spin-summed (as opposed to spin-averaged ) squared matrix element. We denote final-state momenta with k i instead of p i . Consistent with our assumption that there is no appreciable background of χ, we set its occupation factor f χ ( ⃗ k χ ) to zero. All baryonic species abide by zero-temperature Fermi distributions characterized by distinct Fermi momenta p F,B .
We briefly discuss important qualitative features of the evaluation of Eq. (5.4) for the decay process B → χγ and present the corresponding results. We relegate details of the calculation to App. E. We work in the interaction basis, so that the decay proceeds via the Feynman diagram containing the n − χ interaction and the baryon magnetic dipole moment operator, which we write as noting g n = 3.826 and g Λ = −1.226 [50]. This computation is made in a background meanfield of neutron matter, and the associated decay amplitude, as developed in Sec. 4 A, is determined by replacing the canonical momenta of the in-vacuum computation with kinetic momenta as per Eq. (4.5). Labeling canonical momenta as B(p B ) → χ(k χ ) + γ(k γ ), the corresponding spin-summed squared matrix element is where the argument of Γ c.v. follows from our earlier frame choice. 5 Henceforth we abbreviate The prefactor of 2 comes from the baryon's two spin degrees of freedom, and the factor of we arrive at the following result: in which This is a simple consequence of larger neutron number fractions at these densities, and the two rates often differ by several orders of magnitude. However, Λs have a further reach in m χ when they are present than neutrons do, owing to the larger total energy of Λs in neutron matter. 6 Of course, the EoS that do not contain hyperons will not lead to Λ → χγ decays within neutron stars. MeV. Note that the color scales are different between the two panels. See text for additional details.

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

Spontaneous B → χ Conversion
The existence of χ raises the possibility that the baryons to which they couple might undergo spontaneous conversion to χ in the neutron-star medium as they propagate. Such an effect could prove loosely analogous to empirically observed matter-enhanced neutrino oscillations [104] or to the possibility of neutron-antineutron oscillations [105][106][107], breaking baryon number by two units. In the latter case the presence of external interactions from matter or magnetic fields modify the energy of the n andn differently, severely reducing the spontaneous oscillation probability for a fixed source of new physics [108], and the crosssection for scattering-mediated n-n conversion is also very small [109]. In this section, we note the distinct features of B-χ conversion.
The essential physics is thus: B and χ constitute a two-level quantum system. As we have noted in Sec. 5 A, if the coupling ε Bχ is nonzero, then B and χ constitute the interaction basis, whereas the eigenstates of this Hamiltonian, which we term f 1 and f 2 for this discussion, constitute the mass basis. Formally, the strong interactions that operate in neutron matter only ever produce n -this is what it means for B to be an interaction eigenstate. This B is, however, a coherent superposition of f 1 and f 2 at the moment of its creation. The subsequent evolution of this coherent wavepacket depends on the details of the B −χ system.
These details are discussed in depth in App. D; we pick out the most relevant results as they pertain to this discussion.
The Hamiltonian that describes our two-state system depends on the local environment: the total energy of the baryon depends on the density through m * B and Σ B , and baryons with different n.m.-frame momenta will mix differently with χ because Lorentz invariance is spontaneously broken by the background. There exists a resonance in this system wherever the condition, which follows from energy-momentum conservation of the canonical momenta, is satisfied. We expect that this condition will occur for at most one value of the (magnitude of the) baryon momentum for a given density. Moreover, Eq. up to O(ε 2 Bχ ) corrections. If the system is far from resonance, then these eigenvalues are well separated. As a result, the B states produced in scattering processes will essentially immediately decohere into their component f 1 and f 2 with, respectively, probabilities of cos 2 θ and sin 2 θ. As such, the state that emerges from the scattering process manifests as either f 1 with probability cos 2 θ ∼ 1 or f 2 with probability sin 2 θ ∼ (ε Bχ /δω) 2 , and the latter may be vanishingly small -and thus so would be any yield in χ. This means that when B is produced in some strong interaction, the wavepacket containing f 1 and f 2 may remain coherent over relatively long timescales. This is analogous to how neutrino mass eigenstates remain coherent as they propagate in terrestrial oscillation experiments, despite being formed in an interaction eigenstate. 7 As in the case of neutrino oscillations, the f 1 and f 2 components of the B state generically evolve with different phases; over time, this leads to nonzero overlap between the evolved state and either B or χ. The state is then measured, in a sense, at its next interaction some time t later, either by its environment or by some experimental apparatus. It is appropriate, in this case, to invoke the concept of an oscillation probability; this is estimated by When the state is observed, however, it collapses to the combination of f 1 and f 2 appropriate to either B or χ with probabilities given by Eq. (5.15), and the process repeats for further interactions. While the oscillations have a large amplitude (sin 2 2θ ∼ O(1)) in this regime, the probability to convert will remain small if the time between successive measurements δt meas is small, in the sense (∆ω * )(δt meas ) ≪ 1. This is precisely the quantum Zeno effect [112,113].
It remains to determine the timescale of the interactions in the nuclear medium in order to estimate the rate of B → χ conversions. We estimate this to be the light time of the mean interparticle separation around nuclear saturation density: δt strong ∼ n One might expect that this would multiply the large density of baryons to yield a macroscopically relevant rate. However, the near-resonance region occupies a thin shell (parametrically of width ∼ ε Bχ ) within the baryon Fermi sphere; the fraction of baryons relevant for this phenomenon is fantastically small, even in the best case scenario. Thus we summarize by emphasizing that we do not expect B − χ conversion to be a phenomenologically relevant mechanism for the production of χ.

D. Total Rates
In this section, we report the total baryon decay rates that emerge after integrating our earlier results over the structure of a neutron star with a given central density, n c . For example, in Fig. 10, we show the rates that result from integrating the local BNV rates in

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

Combining Individual Limits
Here we briefly describe our statistical procedure for combining limits on ε Bχ derived from different pulsar binary systems. The limits we show have implicitly been determined as contours of constant χ 2 (m χ , ε Bχ ). Our assumed-true hypothesis is that rate of BNVinduced binary spin-down vanishes in these systems, so we have χ 2 = 0 for ε Bχ = 0. As such, each χ 2 function is generically of the form The first equality follows from the fact thatṖ b /P b ∝ ε 2 Bχ , noting Eq. (3.22), and we emphasize that F is a function of m χ only. The limits we have shown correspond to χ 2 = c; 8 we call the resulting curve ε(m χ ). From this, we determine this allows to determine the χ 2 function over the entire parameter space.
The combined limit, then, corresponds to the contour along which the sums of the individual χ 2 functions also equals c. Using the definitions above, we determine the combined limit ε comb (m χ ) as follows: This discussion has assumed that all ε i are defined at the same level c, and that the desired combined limit is also at c. This result can be generalized for distinct individual significances c i and combined significance C: We show our individual pulsar limits as well as our combined limits, realized via our described procedure, for the DS(CMF)-1 EoS in Fig. 12. show the constraint derived for this maximal neutron star. This constraint has been shown in dot-dashing to indicate that it is qualitatively different from the others.
We underscore that we have fixed the masses of these neutron stars to their best-fit values to construct these limits. A more statistically complete analysis would propagate the We also show constraints from KamLAND [114], SuperKamiokande [115], and BESIII [116].
light blue. These are as much as twenty orders of magnitude stronger than the constraints we have derived, but we note that these are only operative up to m χ = 920 and 827 MeV, respectively. This is a result of experimental cuts -heavier χs result in less energetic photons in the decay, and eventually these become too soft to be meaningfully detected.
We emphasize, in particular, that these experiments cannot probe the region m χ > m n ; while they are more powerful when they are operative, they are fundamentally constrained in ways that astrophysical probes of new physics are not. For Λs, we show the constraint on invisible decays from BESIII [116] in dark cyan. In this case, we find the opposite result: pulsar binaries are able to probe this branching ratio as much as twenty orders of magnitude more severely than laboratory constraints! The caveat is that this requires hyperons to appear in neutron stars, which is still a matter of debate, simply because EoSs without hyperons exist that confront current observational data successfully. However, if hyperons appear in an appreciable amount in these objects, then one can expect vast improvements on laboratory searches.
The upper panel of Fig. 14 is incomplete in that there are additional constraints around m χ ≈ m n , a region that has become of interest in recent years as a result of tests of newphysics explanations [10] of the neutron lifetime anomaly [12]. We examine this region more closely in Fig. 15; panel (a) casts these searches in terms of constraints on ε nχ , while panel (b) casts them in terms of constraints on Br(n → χγ). We show in blue the estimated constraint from a direct search for n → χγ using ultra-cold neutrons (UCN) [103], and in green we show a constraint from Borexino from searches for hydrogen decay, both from Ref. [34]. We also show the curve along which the free hydrogen lifetime is supposed to be τ H = 10 32 s in dashed gold, also from Ref. [34]. (The constraints from Ref. [34] are reported at 90% CL, though the differences between those and limits at 2σ should be very small given the ranges shown in the figure.) Clearly, neutron stars are more sensitive to these decays than these (would-be) laboratory constraints by many orders of magnitude.
It was noted in Ref. [10] that the existence of χ can destabilize nuclear matter, including 9 Be. This constraint was calculated more precisely in Ref. [49], assuming that the lifetime of 9 Be is longer than 3 × 10 9 years to account for the presence of 9 Be in old, metal-poor stars [117]. This constraint is shown in red in Fig. 15 and is competitive with (if not dominant to) our neutron star constraints in the region of its operation, m χ < 937.993 MeV. We note that other probes of dark decays of nuclei with low neutron separation energies have been discussed in, e.g., Ref. [47]. Particular attention has been paid to decays of 11 Be, with experimental efforts underway at CERN-ISOLDE [118] and ISAC-TRIUMF [119], though we are unaware of any efforts to interpret these experimental results as constraints on new physics. As a side note, it is curious that there are no laboratory constraints, as far as we can tell, on the lifetime of 9 Be. We find the arguments about the presence of 9 Be in old stars compelling and agree that this is a valid constraint, but we are surprised, frankly, that the lifetime is only constrained at the billion-year scale. While experimentalists of yore would have had little reason to interrogate the stability of 9 Be -or indeed, any species thought to be stable in the SM -we regard the observation that the stability of these systems has not been tested in a detailed way in the laboratory as a potentially promising avenue for constraining new physics.
We conclude by noting that Ref. [49] has also presented constraints on n → χγ from cosmology and from neutron star cooling. The former is a combination of constraints coming from modifications to Big Bang Nucleosynthesis (BBN) and the Cosmic Microwave Background (CMB); this treatment includes the reverse decay χ → nγ when m χ > m n , and so constrains the region shown. However, in their calculations, χ is assumed to constitute (at least some of) the dark matter. This is unlike our framework, in which we introduced more new states (ξ and ϕ B ) to prevent overaccumulation of χ. Therefore, the limits they derive from BBN and CMB do not apply here, though we agree that this would be an interesting and important avenue to explore. The neutron star cooling constraint derived there, however, makes very rough assumptions about how heat from decays is deposited into the neutron star, with the implicit assumption that increases in the temperature of the core of the neutron star lead to commensurate increases in the observed effective temperature. Yet thermal transport and cooling in neutron stars demands careful investigation; for instance, BNV decays lead to β-disequilibrium, which leads to neutrino cooling via (direct and modified) Urca processes, which impact how the energy released in the decays is deposited back into the SM fluid. While we agree that while old, cold neutron stars should constrain this model, the details are intricate and expected to be sufficiently impactful that we decline to include such constraints here.

IMPLICATIONS FOR MODELS OF BARYOGENESIS AND DARK MATTER
The prospect of explaining the origins of both the dark matter abundance and the cosmic baryon asymmetry within a single dynamical framework is a beguiling one. Different possibilities have existed for some time, and many share a common feature: there is a dark-sector baryon that carries baryon number and into which SM baryons can decay. A particularly intriguing variant is that of B mesogenesis [14,29,120]. It proceeds in the early universe thesis. Finally, it is an example of a testable mechanism of baryogenesis [25], in that its essential features are subject to direct experimental investigation. Particularly, its reliance on the SM mechanism of CP violation (albeit new CPV sources could enter) implies that the branching ratios of B mesons in SM baryons and the dark fermion (antibaryon) cannot be too small, with the expectation that the branching fractions can roughly be no less than Br(B 0 s,d → χB) ≳ 10 −5 or Br(B + → χB (+) ) ≳ 10 −6 [120]. The expected theoretical window in χ mass is 0.94 GeV < m χ < 4.34 GeV [29]. Studies from Belle [121] and BaBar [122] limit the available parameter space in the mass region of 1 − 4.4 GeV, and it is anticipated that the remaining parameter space can be probed at Belle-II [122]. This model is particularly close to the model we study, in both its visible and hidden-sector components. In this paper we have established severe limits on the ε nχ and ε Λχ mixing parameters for χ masses satisfying m χ ≲ 1400 MeV, as shown in Fig. 13. In this mass region and for the regions of hidden-sector parameter space we have chosen, our limits constrain the flavor structure of models of B-mesogenesis, and we now turn to those and their implications.
Different UV completions of B-mesogenesis models fare differently in light of our constraints. Here we consider versions in which only one extra particle is needed. For example, in Ref. [14], a color-triplet, SU(2) L singlet scalar with the SM quantum numbers (3, 1, −1/3) is used, though a scalar of form (3, 1, +2/3) [29] or a vector of form (3, 2, −1/6) [16] are noted alternatives. We do not consider this list exhaustive. The two scalars are just the leptoquarks we have noted in Sec. 2: S * 1 andS * 1 [10,15]. The phenomenology of these specific models has been studied, and in order to explain the baryon asymmetry, the dark matter abundance, and all empirical constraints, including those on |∆F | = 2 meson mixing, a rich flavor pattern of couplings to quarks is needed [29].
To determine the implications of our constraints, we first note the structure of the Lagrangian for each UV completion, following Ref. [16], though we write our 2-spinors as in Ref. [123] and employ the conventions given there. Denoting the new scalars as Y Y and the new vector as X µ , we have where ε is an antisymmetric tensor in the two-spinor indices and χ is a right-handed field.
With the B assignments of −2/3 for the scalars Y 2 3 and Y − 1 3 and B = 1 for χ, the noted interactions conserve baryon number. In Refs. [15,29] y QaQ b (for each a, b) is taken to be zero. The color structure of the first term of Eq. (7.1) requires that the product of dlike quarks be antisymmetric in the generation indices a, b, which follows because we have assumed the scalar is a color triplet. As for the last case, the vector X µ can be written in two-spinor form as [16] and thus through Eq. (7.3) we see that both scalars couple to left-handed quarks. We have defined our scalar-fermion couplings in the flavor basis, rather than the mass basis, but in the case of couplings to right-handed quarks no distinction needs be made. However, in the case of couplings to left-handed quarks we need to rotate the fields to the mass basis, to parallel the treatment of the charged weak current in the SM. As a result, a flavor diagonal coupling to a left-handed quark of a single flavor can engender a contribution to a flavor-changing neutral current (FCNC). In the example of Z ′ models, satisfying FCNC constraints with a large Z ′ coupling requires nearly flavor-universal couplings [124], where we note that in the flavor universal limit the unitary structure of the CKM matrix makes the FCNC couplings vanish. We will see that this effect does not appear here because our scalars do not ever couple to two left-handed quarks of the same flavor. Replacing a left-handed flavor state d i with a combination of mass states via V ij d j , with V the CKM matrix, we see that the X µ completion does lead to a FCNC of form [16] where we have employed 4-component notation. This interaction engenders not only |∆F | = 2 meson-mixing but also structures such as B (s) →K or B (s) → π 0 at tree level, which can be probed through B decay studies. We also see explicitly that the structure of the vertex does not require a flavor universal coupling to control the size of the effect. Thus there are no particular flavor conspiracies in satisfying the |∆F | = 2 constraints, and to determine the impact of the constraints we have found on the mixing parameters ε nχ and ε Λγ on these models, it suffices to consider the contributions to these quantities from the scalar-fermion couplings with a particular UV complete model.
Considering, then, the flavor structure of Eq. (7.1) we see that n → χγ cannot occur at tree level, and a loop graph with W and Y Y exchange is needed to generate the process [15].
The opposite situation is true for Λ → χγ, with Eq. (7.1) and Eq. (7.2) yielding that process at tree level and one loop level, respectively. The pertinent Feynman diagrams are illustrated in Fig. 16, replacing the illustration of Fig.(1). Noting Eqs. (2.5) and (2.7), it is apparent that the mixing parameters ε nχ and ε Λχ depend very differently on the underlying scalar-fermion couplings in the two cases -we refer to Ref. [15] for explicit expressions.
In particular, the one-loop diagrams bring in a coupling to the b quark as well, with the following combinations of couplings: y db y χu ; y sb y χu (7.6) y db y χc ; y sb y χc (7.7) y db y χt ; y sb y χt (7.8) each of which could saturate the bound we have found for ε nχ . In regards to the mechanism of B-mesogenesis, operators with the flavor combinations χbud, χbud, χbcd, and χbcs are as discussed in text, after Ref. [15].
pertinent, and they take one of three forms [29] θ (1) where i ∈ d, s, j ∈ u, c, and the colors have been contracted to form a color singlet in each case. Taking the couplings in Eq. (7.8) one at a time, we find that saturating our ε nχ constraint we have found limits the coefficient of each of the θ ij operators to be powers of ten smaller than that needed for B-mesogenesis to be successful [29]. We emphasize, however, that this is particular to the mass window in χ and region of hidden-sector parameter space we have noted. For the Y2 3 scalar, those are the operators that would act -and thus we have ruled out this specific model for B-mesogenesis under the conditions we have noted.
The other UV completions we have considered are not similarly constrained, because the ε χΛ constraints limit just the flavor combinations y χb y du and y χb y dc pertinent to B-mesogenesis -the other flavor combinations associated with θ (1) ij and θ (2) ij remain unconstrained despite the severity of our limits.

SUMMARY
BNV has not yet been observed in terrestrial experiments, and its deep ties to explanations of the observationally well-established cosmic baryon asymmetry [2] argue persuasively for its investigation on broader fronts. Previously, we have considered how it might eventually be discovered through precision measurements of neutron star observables, particularly those of changes in the binary-pulsar period, familiar from tests of general relativity [5].
Thus far we have found limits, and they are macroscopic ones, in that they emerge from the consideration of a neutron star as a whole. Such constraints miss a concrete connection to particle physics, and it is badly needed: regardless of whether we continue to constrain or, finally, discern the existence of BNV (in contradistinction to a failure of general relativity) from these studies, further theoretical progress on the problem of BNV requires constraints on the particle physics models of BNV themselves. In this paper, we have developed just such a connection, using a concrete description of the neutron star interior based on a relativistic mean-field theory in hadronic degrees of freedom [17][18][19] that successfully confronts existing macroscopic properties of neutron stars [21]. Within this context, we have developed how to assess the rates for BNV particle processes in dense matter, and we present explicit rates for benchmark processes, particularly B → χγ, considering its rate both at local points within a neutron star as well as its volume rate after integration over the structure of the entire star, up to its crust. Although our in-medium formalism is germane to the evaluation of any particle process in the dense medium of a neutron star, the focus of this paper -noting current sensitivities -is that of apparent BNV through baryon decays to hidden-sector particles. Finally, with this in place, we match the computed rate to our inferred limits on anomalous binary-period lengthening, i.e., how the binary itself spins down, to set one-sided limits at 2σ on the mixing parameters ε Bχ , for individual binary-pulsar systems, as well as a combined limit for all of the studied systems.
As a result of these studies, we discover that neutron stars open new windows on the study of BNV, probing m χ parameter space not accessible to terrestrial nucleon decay experiments, due to experimental limitations in the detection of a final-state photon. More than this, the dense nuclear medium admits the study of regions for which m χ exceeds the vacuum mass of the nucleon, as well as the possibility of probing strange baryon decays. Our final limits are reported in Figs. 14 and 15. We observe that in the regions of parameter space to which proton decay (nuclear stability) experiments are sensitive [114,115], they exceed the limits we set by nearly twenty orders of magnitude. In contrast, however, our neutron star limits exceed the sensitivity of those from terrestrial Λ and neutron β-decay experiments by a comparably large amount. Let us emphasize that our limits are likely upper bounds, and hence are conservative, in that they are determined by the electromagnetic decay B → χγ alone, although the particle physics models we study do admit the possibility of B → χ + meson(s) decays as well. This latter set of decays has no reason to be negligible compared to the electromagnetic decays in rate -and we note Ref. [16] for specific examples computed within (in-vacuum) chiral EFT [51]. As a result, we would expect larger B decay rates for fixed ε Bχ , but the challenges in realizing a suitable theoretical assessment of the hadronic channels prompt the conservative approach we have espoused in this paper.
We now turn to an assessment of the limitations in our approach. One key question concerns the largest value of ε Bχ , ε max Bχ , we can possibly limit with our formalism, in which the SM drives the dynamical response of the neutron star to BNV. (In our work, dark-sector interactions drive the removal of χ, so that the neutron star survival constraints on the mass of m χ noted in Refs. [101,125,126] do not operate.) We believe a realistic assessment of ε max Bχ requires a study of neutron star heating from relatively fast rates of BNV, the complexities of which lie beyond the scope of this paper. We note, however, the outcomes of terrestrial neutron β-decay searches [45], shown in Fig. 15, as well as limits arising from constraints due to the charged-current structure of the SM [41], noted in Eq. (2.1). Since n → χγ does not derive from a SM weak process in any way, a Br(n → χγ) limit of O(10 −3 ) implies a limit on ε nχ of O(10 −9 )! Thus we think these limits are severe enough that determining ε max Bχ precisely is not an immediate concern, but, rather, an important topic for future investigation.
Another potential limitation may be our use of a relativistic mean-field theory framework [17][18][19] in which to describe the nuclear medium within a neutron star. This approach is computationally tractable and readily allows for the treatment of more sophisticated models of the nucleon-nucleon interaction than those in which it was first devised. We have employed the chiral SU(3) hadronic model of Refs. [9,21,91] in this paper. This is admittedly a model that is not QCD, and our ability to assess the errors predicated by this choice is rather limited. We have, however, studied how our results change within a family of EoSs, namely DS(CMF) 1-8 EoSs [92,93], to which it can be connected. Moreover, frankly, there is no other alternative for the treatment of dense nuclear matter, though this may ultimately change [127]. We note that the use of chiral effective theory has been championed in this regard [87], but its applicability does not stretch much beyond that of nuclear saturation density. In the future, it may be advantageous to consider EoSs that blend the chiral effective field theory and relativistic mean-field theory approaches [94]. Nevertheless, given our interest in order-of-magnitude estimates, we believe that our choice is also reasonably realistic.
Different paths beckon as opportunities for future work. We believe that studies of neutron star heating from BNV is important not only to discerning the limits of our existing formalism, but also, crucially, to interpreting what a significant observation of anomalous binary spin down might mean. It strikes us that theoretical heating studies and concomitant observational studies of neutron star cooling may be the only tangible way to tell a failure of general relativity, in some undetermined way, from BNV.
As for other possibilities, we could consider how our results could change if the neutron star were a hybrid star, containing a quark core [9], or how viable models with a significant χ admixture in the neutron star (albeit constrained by Eq. (2.1) [41]), such as that of Ref. [37], could be addressed through modifications of our formalism. As for future terrestrial experiments that could complement the studies of this paper, it strikes us that empirical studies of the lifetime of SM-stable composites, such as atomic hydrogen, or of the 9 Be nucleus, could yield fruitful results.
in which σ µν ≡ (i/2)[γ µ , γ ν ] and γ µ are the usual Dirac matrices. The in-medium Gordon decomposition is then given by in which we defined q ν ≡ p ′ ν − p ν . The general form of a vector interaction vertex, Γ µ , can be written as in which A, B, C, D are functions of scalar quantities (e.g., q 2 ). Applying the Ward identity, q µ Γ µ = 0, plus p ′ * 2 = p * 2 = m * 2 and p ′ 2 − p 2 = 2q · Σ, yields C = 0 and 2B = D. The electromagnetic vertex factor can then be written as in which F * 1,2 are in principle distinct from their vacuum counterparts F 1,2 . We now show how the electric charge can be identified in the scattering amplitude of a baryon from a Coulomb potential A µ = (Φ(x), ⃗ 0). Employing equations u(k * , λ)u(k * , λ) = 2m * and u(k * , λ)γ 0 u(k * , λ) = 2E * (k * ), this amplitude can be written as in which E * = m * 2 + (⃗ p * ) 2 , and χ is the Pauli spinor. The electric charge (q) can then be identified, by considering this scattering in the c.v. frame (⃗ p * = 0), as q = Similarly, we can identify the magnetic moment from the scattering amplitude of a baryon from a static magnetic field potential A µ = (0, ⃗ A) at small momentum transfers (q 2 ≈ 0), which is given by The first term inside the bracket can be written as in which σ i are the Pauli matrices, and χ, η represent the spin states. This expression can be further simplified using σ i σ j = δ ij + iϵ ijk σ k , such that The F 2 term in the scattering amplitude (A.6) already contains a factor of q, and so we can evaluate it using the leading order expansion of the spinors in the non-relativistic limit (⃗ p * ≪ m * ), which is given by u(⃗ p * = 0) = √ 2m * (χ, 0) T . We also note that such that the spin-dependent contribution from Eq. (A.10), i.e., u(p ′ * )(σ i0 q 0 )u(p * ) is proportional to q 0 q j , which is subdominant to other terms. The term from Eq. (A.9), i.e., u(p ′ * )(σ ij q j )u(p * ) is given by The amplitude in Eq. (A.6) can then be written as (note q j = −q j ) in which we defined the magnetic field byB k ≡ −iϵ ijk q iÃj , spin by ⃗ S ≡ (1/2)χ † ⃗ ση, and the baryon g-factor can be identified as g * = 2 [F * 1 (0) + F * 2 (0)].

Appendix B: Nonrelativistic Limit of In-Medium Scattering
In this appendix we study the non-relativistic (NR) limit of the RMF model, and derive the elastic scattering formalism in the Born approximation. Since the medium effects in RMF formalism resemble an electromagnetic interaction with a constant EM background field given by eA µ → Σ µ , it is instructive to consider the NR limit of baryon EM interactions in medium. We explicitly show that the NR limit of the modified Dirac (Eq. (4.6)) solutions under the influence of EM interactions, reduces to the two-component Pauli spin theory, with replacements m → m * , eΦ → eΦ + Σ 0 , e ⃗ A → e ⃗ A + ⃗ Σ, in which Σ 0 and ⃗ Σ are the self-energies due to the medium effects, e is the baryon electric charge, with Φ and ⃗ A as the scalar and vector EM potentials respectively. We start from the Schrodinger equation, which can be written by denoting the Dirac wave-function (ψ) in two-component notation [128], Using the definition (φ,χ) = exp(−im * t) (φ, χ), we can rewrite Eq. (B.1) as We note that in the NR limit, in which kinetic and interaction energies are much smaller than m * , the second component χ is subdominant to the first component φ and is approximately given by We also arrive at the Pauli equation governing the first component (φ): This expression can be further simplified for a weak uniform magnetic in which ⃗ p * = ⃗ p − ⃗ Σ is the kinetic three momentum, and ⃗ L * = ⃗ r × ⃗ p * and ⃗ S = ⃗ σ/2 are the baryon's kinetic orbital angular momentum and spin respectively. Note that in the n.m.
We now construct the elastic scattering formalism off of an arbitrary potential (V ) in with solutions of the form which can also be written in a more symmetric way in terms of ⃗ p * . If we orient our coordinates such that ⃗ Σ.⃗ x > 0, then for a positive ⃗ p (⃗ p.⃗ x > 0) the first term is a plane wave moving to the right and the second term is a wave moving to the left. Therefore, we pick the first term for incident waves in the elastic scattering problem. Let H 0 be the Hamiltonian used in Eq. (B.4) (with Φ = ⃗ A = 0), and |k (+) ⟩ be the state that satisfies the following Schrodinger equation in the presence of a potential V then, |k (+) ⟩ can be found from the Lippmann-Schwinger equation: The momentum representation of operator G ≡ (E − H 0 + iε) −1 is given by and the position space representation is given by in which we performed the angular integration in the second line, and the complex contour integration in the third line. To characterize the scattering problem at r → ∞ we approximate the above expression for (r ′ /r) → 0 using R = |⃗ r − ⃗ r ′ | ≈ r −r · ⃗ r ′ , such that We now write the asymptotic form of the Lippmann-Schwinger equation in position space in which ψ k (⃗ r ) ≡ ⟨⃗ r | ⃗ k (+) ⟩ and φ k (⃗ r ) ≡ (2π) −3/2 exp i ⃗ k · ⃗ r . The exponential outside of the integral in the second term is an ellipsoidal wave (stretched along ⃗ Σ) which becomes spherical in the n.m. frame ( ⃗ Σ = 0). The exponent inside the integral is a vector pointing in the direction of | ⃗ k − ⃗ Σ |r + ⃗ Σ, which reduces to the familiar kr term in the n.m. frame. We can see that the gradient of ellipsoidal surface is equal to the vector in the exponent inside the integral since which suggests that the exponent ⃗ k ′ ≡ | ⃗ k − ⃗ Σ |r + ⃗ Σ is the momentum of scattered particle in the direction of an observer at r. Note that the kinetic energy of the scattered particle is given by 16) and the scattering is indeed elastic. We can therefore deduce the scattering amplitude by which is the Fourier transform of the potential in the Born approximation.

Appendix C: In-Medium Compton Scattering
In this section we evaluate the Compton scattering cross section of baryons, B(p 1 ) + γ(k 1 ) → B(p 2 ) + γ(k 2 ) (see Fig. 17) in neutron star medium, denoting the photon and baryon energies by ω 1,2 and E 1,2 respectively. We first note that the second term in the k 1 k 2 baryon propagator defined in Eq. (4.14) vanishes since and similarly, it can be shown that (p * 1 − k 2 ) 2 − (m * B ) 2 is strictly negative. The amplitude for the diagrams shown in Fig. 17 can then be written as in which in which F * 1,2 are the in-medium form factors. The interaction term in the amplitude can be simplified using which follows from ϵ µ (k 1 )k µ 1 = ϵ µ (k 2 )k µ 2 = 0. The spin-averaged squared amplitudes simplify to and with the cross-term given by , (C. 8) in which We now define the following Mandelstam variables such that s * + t * + u * = 2(m * B ) 2 . We suppress the superscripts ("*") of m * B , F * 1,2 in some of the following equations for convenience. The averaged amplitude-squared can be written as 16) in which we note that I and IV are related via (s ↔ u) replacement. Equation (C.13) can then be written as .

(C.17)
We now consider the Compton scattering in the rest (c.v.) frame of B(p 1 ) (see Fig. 18), in which ⃗ p * 1 = 0. We first note that the relationship k 1 · k 2 = p * 1 · (k 1 − k 2 ), written in the c.v. frame, yields ω 1 ω 2 (1 − cos θ) = m * B (ω 1 − ω 2 ). We then arrive at the following kinematics in the c.v. frame which resembles the familiar Compton's formula. We use these kinematical relationships to write Eq. (C.13) in terms of ω 1 and the scattering angle (θ) in the c.v. frame as The phase space integrals over the final states (see Eq. (4.15)) can be written as in which dΩ 2 = d cos(θ) dϕ is the differential solid angle of ⃗ k 2 in the c.v. frame, and f B (⃗ p 2 ) is the Pauli blocking factor for the outgoing baryon. The shape of the Fermi surface in a general frame (such as c.v.) changes from being spherical to an ellipsoid. The general form of f B (⃗ p 2 ) in an arbitrary frame is given by and noting v * B = 0, v * A = 1 in our chosen frame (c.v.), the in-medium Compton scattering differential cross section can be written as in which we recover the Klein-Nishina [100,129] formula if we set f B (⃗ p 2 ) = 0, F 1 = e, F 2 = 0 and replace m * B by m e .

Appendix D: Fermion Mixing in Dense Matter
In this appendix we evaluate the eigenvalues of a system consisting of a neutral baryon (B) and a dark fermion (χ) with a mixing term between them, in the context of RMF framework.
We suppress the superscript ("*") in baryon's effective mass (m * B ) for convenience. The Note that the baryon current J µ B ≡ ψ B γ µ ψ B satisfies and as expected is not conserved anymore, instead, the combined current J µ ≡ J µ B + J µ χ is conserved. The conserved energy-momentum prescribed by the Noether's theorem [100,130] is given by We expand each of the fields in terms of four modes ω ± 1,2 ( ⃗ k) as in which ψ(x, t) stands for ψ B,χ , a 1,2 and b 1,2 are the annihilation operators for particles and anti-particles, ω stands for ω( ⃗ k), we note the inequality ω( ⃗ k) ̸ = ω(− ⃗ k) if ⃗ Σ B ̸ = 0 (see Eq. (4.8)), and the fact that in the presence of medium (Σ 0 B ̸ = 0) the particle and antiparticle energies are not equal anymore (e.g., see Eq. (2.40) of [89]). The coefficients α and in which we note the definition P µ = i∂ µ = (H, − ⃗ P ). We now plug in the field expansion from Eq. (D.7) into Eq. (D.11) to arrive at the equation governing the spectrum of ω ± 1,2 modes: O(δ 2 ) : We can see that the first order term breaks the degeneracy by splitting the energies.
where we note and we define the quantity k * χ ≡ k χ − Σ B = p * B − k γ . Note in the neutron star medium that energy-momentum conservation of the total canonical momentum still holds: p µ B = k µ γ + k µ χ . Consideration of the kinematics show that we need only consider the first term in the full baryon propagator given in Eq. (4.14).

Integrated Rates
We now address full integral over phase space, (E. 6) In the main text, we presented the rate as an integral over the baryon Fermi sphere of the dilated widths of individual baryons. Here, we will contrast this approach with a more straightforward evaluation of this integral and demonstrate that these yield consistent results, as expected.
Our first step in the evaluation of the rate is to separate the integrals over the χ and γ phase spaces and evaluate these first: We have simplified the first integral by noting that it only depends on the magnitude of the three-momentum |⃗ p B | ≡ p, and that we only integrate within the neutron Fermi sphere. We tackle this second integral by computing it in the c.m. frame of the decaying neutron. We note, however, that the matrix element depends on p * B , which has a nonvanishing spatial component, though we will find that this is not relevant for the ultimate evaluation of the integral.
(E. 18) We note that if the self-energy were to vanish (σ = 0) we would recover the vacuum decay rate reported in Eq. (2.7).