Dark Matter interpretation of the neutron decay anomaly

We add to the Standard Model a new fermion $\chi$ with minimal baryon number 1/3. Neutron decay $n \to \chi\chi\chi$ into non-relativistic $\chi$ can account for the neutron decay anomaly, compatibly with bounds from neutron stars. $\chi$ can be Dark Matter, and its cosmological abundance can be generated by freeze-in dominated at $T \sim m_n$. The associated processes $n \to \chi\chi\chi \gamma$, hydrogen decay ${\rm H}\to \chi\chi \chi\nu(\gamma)$ and DM-induced neutron disappearance $\bar\chi n \to \chi \chi (\gamma)$ have rates below experimental bounds and can be of interest for future experiments.

These measurements contradict at about 4.6σ level the Standard Model (SM) that predicts Γ tot n = Γ β n up to sub-leading channels that can be neglected. If confirmed, they indicate new physics. It has been proposed that could be due to an extra decay channel of the neutron that does not produce protons. The new decay mode needs to be nearly invisible with a branching ratio around 1%. Various authors proposed a n → χγ decay into a new neutral fermion χ with mass M slightly below the neutron mass so that E γ ≈ m n − M is small [12,13] (see also [14] and [15,16]). However, this interpretation has the following problems: 1. First, if M < m p + m e the new particle χ does not promptly decay back to SM charged particles and can be Dark Matter; but in this range E γ > m n − m p − m e is large enough that the extra decay is visible enough that it has been tested and disfavoured [17,18,13]. An alternative that can explain both the neutron anomaly and Dark Matter (DM) would be welcome.
2. Second, n ↔ χγ would lead to χ thermalization inside neutron stars, as they are older than a Myr, while the neutron decay anomaly needs the much short time-scale of eq. (1). Thermalized χ soften the equation of state ℘(ρ) of neutron stars. Free χ add more energy density ρ than pressure ℘ reducing the maximal mass of neutron stars below 0.7M sun , like in the 1939 Oppenheimer and Volkoff's calculation that assumed free neutrons, neglecting their nuclear self-repulsion [19].
3. Third, the precise SM prediction of the neutron decay rate, possibly Γ SM n = 1/(878.7 ± 0.6) s [20,21], agrees with the bottle experiments disfavouring new physics contributions that increase the neutron decay rate. This might indicate that the neutron decay anomaly could just be due to some under-estimated systematic uncertainty of 'beam' experiments, or that it needs alternative interpretations where the neutron oscillates (with a resonant enhancement thanks to magnetic fields in the apparata) into a composite twin of the neutron [21][22][23]. 1 However, the quoted precise SM prediction for Γ SM n employs determinations of V ud and measurements of the neutron axial coupling constant g A that are potentially problematic at the claimed level of precision. In particular, the pre-2002 measurements of g A gave instead a Γ SM n value close to the beam experiments [20]. We thereby proceed ignoring this possible third problem.
A proposed solution to the second problem is adding some new light mediator such that χ undergoes repulsive interactions (with itself or with neutrons) stronger than the QCD repulsion among neutrons: this would make energetically favourable to avoid producing χ in thermal equilibrium inside a neutron star, stiffening the equation of state [34][35][36][37][38]. However this requires QCD-size scattering cross sections that, if χ is DM, risk conflicting with the bullet-cluster bound on DM interactions. Again, an alternative that can explain both the neutron anomaly and Dark Matter would be welcome.
In section 2 we discuss how more general models where new physics opens different neutron decay channels affect neutron stars depending on the conserved baryon number B χ assigned to χ. In section 3 we show that the minimal choice B χ = 1/3 leads to n → χχχ decays into free DM χ that give a milder modification of neutron stars, compatible with observations. In section 4 we elaborate on the theory behind and compute related processes. In section 5 we show that χ can be DM, with abundance possibly generated via freeze-in. In section 6 we show that χ DM has unusual signals compatible with current bounds. Conclusions are presented in section 7.

Overview of possible models
To motivate the model we will propose, we first classify models where SM particles couple to one new particle, by describing how baryon number B and lepton number L can be assigned to the new particle such that they are conserved, in order to satisfy the strong phenomenological bounds on B, L violations.
The L and B charges of the new particle determine its spin needed to couple to SM particles via operators with lowest dimension. Since in the SM only fermions carry lepton and baryon number, the new particle must have charge (−1) 3B+L under the Z 2 symmetry that flips signs of all fermions and leaves bosons unchanged. Table 1 lists the main possibilities, focusing on the minimal spin: either a fermion χ with spin 1/2, or a boson φ with spin 0 (or, equivalently, two fermions χχ). See [39][40][41][42][43] for lists of SM operators, partially relevant for our study.
The apparently minimal choice (upper row of table 1) arises when a fermion χ carries baryon number B χ = 1 and vanishing lepton number L χ = 0. This is the model of [12,13] where the neutron decays as n → χγ. Models that involve a boson φ not protected by Fermi  repulsion worsen the neutron star problem. Models where a fermion χ carries both lepton and baryon number do not avoid the neutron star problem, as neutrinos freely escape from neutron stars. Independently of the specific particle-physics interactions, the χ chemical potential µ χ in thermal equilibrium in a neutron star is fixed in terms of the chemical potential of conserved charges (baryon number and vanishing electric charge) as Eq. (2) shows that all models affect neutron stars in a qualitatively similar way. Eq. (2) also shows how important quantitative differences can arise: reducing µ χ allows to substantially reduce the impact on neutron stars. This is achieved noticing that a more minimal choice of quantum numbers exists, given that baryon number can be fractional, like the electric charge.
In the more minimal model χ carries baryon number B χ = 1/3 and L χ = 0 (second row of table 1). So neutrons decay as n → χχχ if the new particle is light enough, M < m n /3. This neutron decay can have a small enough impact on neutron stars for the same reason why ordinary neutron decay n → peν e has a small impact: the neutron decays to particles light enough that their Fermi repulsion is big enough, without the need of introducing extra repulsive interactions. We thereby focus on this model. The DM mass M is strongly restricted as illustrated in fig. 1, where we plot the following key kinematical thresholds: • One needs M < m n /3 ≈ 313. 19 MeV so that n → χχχ is kinematically open.
• One needs M > (m n − E Be )/3 ≈ 312.63 MeV so that nuclear decays into χχχ are kinematically closed, where E Be = 1.664 MeV since the strongest bound comes from 8 Be [12]. The neutron decay anomaly can be explained in the mass range highlighted in green in fig. 1, while it cannot be explained in the red shaded range.

Bounds from neutron stars
Neutron stars are described by the Tolman-Oppenheimer-Volkoff (TOV) equations [19,44] d℘ for the pressure ℘(r) at radius r and for the mass M(r) inside radius r. The TOV equations describe spherical hydrostatic equilibrium in general relativity and can be solved for any equation of state that gives the density ρ in terms of ℘, by starting with an arbitrary pressure at the center at r = 0 where M(0) = 0, and by evolving outwards until finding the neutron star radius r = R at which ρ(R) = 0. Repeating this procedure predicts the relation between the radius R and the total mass M, in particular giving a maximal mass M above which neutron stars become unstable and collapse into black holes. In the SM this maximal mass is around two solar masses, compatibly with observations of neutron stars around this mass.
Adding new stable particles χ to neutrons and to the other sub-dominant SM particles, the total energy density and pressure become ρ = ρ n + ρ χ and ℘ = ℘ n + ℘ χ . The equation of state with the two components in equilibrium is computed as follows. We assume that the χ particles are free fermions with mass M and g = 2 degrees of freedom, so that their levels are occupied up to their Fermi momentum p χ . Their number density and pressure are Thermal equilibrium of n ↔ χχχ relates the χ chemical potential µ χ = M 2 + p 2 χ to the chemical potential of neutrons as in eq. (2), i.e. µ χ = µ n /3.  [45], the dotted curve is ℘ n = κ 0 ρ 2 n ); adding n ↔ χχχ (blue curves); adding n ↔ χγ (red curves). Right: the corresponding relation between the radius and mass of neutron stars. This shows that the observed neutron stars with mass around two solar masses are compatible with n ↔ χχχ, but not with n ↔ χγ. The solutions below the peaks at smaller radii are unstable.
Concerning the SM particles (mostly neutrons), their equation of state is significantly affected by nuclear repulsion, and different approximations give somehow different results [46-48, 45, 49, 50]. We adopt the following computations: 1. A first, rough but simple, equation of state that includes nuclear effects is ℘ n = κ 0 ρ 2 n with κ 0 ≈ 52/ GeV 4 [46]. The thermodynamic relation dn n /n n = dρ n /(ρ n + ℘ n ) of eq. (5) determines n n up to a constant, fixed imposing µ n = dρ n /dn n → m n as n n → 0. The result is ρ n = m n n n 1 − m n n n κ 0 , This equation of state is shown by the black dashed curve in fig. 2a.
equilibrium implies the relation µ i + µ j = µ i + µ j among the chemical potentials, independently of particlephysics details. We recall the thermodynamic relations at T = 0 used later. Particles are added with energy E = µ at the Fermi sphere, so dρ = µ dn. The pressure ℘ is then given by having used U = N E and N = nV . 2. The possibly precise BSk24 equation of state from [45], shown by the black curve in fig. 2a, agrees well with neutron star observations. We use its numerical form, and present a rough analytic Taylor approximation: ρ n = a 1 n n + a 2 n 2 n + a 3 n 3 n with a 1 ≈ m n , a 2 ≈ 16/ GeV 2 and a 3 ≈ 4800/ GeV 5 . Then thermodynamic relations imply µ n = m n + 2a 2 n n + 3a 3 n 2 n and ℘ n = a 2 n 2 n + 2a 3 n 3 n .
The red curves in fig. 2a show how the equations of state computed in the SM get significantly softened if n ↔ χγ is in thermal equilibrium (i.e. µ n = µ χ ) with M = m n . The blue curves in fig. 2a show how thermal equilibrium of n ↔ χχχ (i.e. µ n = 3µ χ ) with M = m n /3 leads to a milder modification of the equation of state, that remains almost as hard as in the SM. As a result, fig. 2b shows that the relation between the neutron star mass M and radius R in the presence of n ↔ χχχ remains close enough to the SM limit (in particular allowing for observed neutron stars with mass M ≈ 2M sun ), unlike what happens if n ↔ χγ is in thermal equilibrium. In particular, while n ↔ χγ reduces the maximal mass of neutron stars in contradiction with data, n ↔ χχχ leads to a mild reduction comparable to current SM uncertainties. The neutron star radius is reduced in a similarly mild way compatible with data. Different SM computations lead to maximal neutron star masses between 1.8M sun and 2.6M sun and minimal radii between 10 km and 14 km [47] in apparent agreement with data. Thereby, if the SM can account for observed neutron stars, its n → χχχ extension can too. More precise future computations and observations might be able to test the mild difference. So far we considered N = 1 generation of χ particles. Fig. 3 shows the maximal neutron star mass as function of the number N of χ generations, assuming for simplicity they all have the same mass. We see that n ↔ χχχ leads to such a small reduction that also N > 1 is allowed. The SM corresponds to N = 0.

Possible theories
Having established that n ↔ χχχ is compatible with neutron star bounds, we study its possible theory. Since conserved χ number is needed, χ must be a complex fermion(s) described by Dirac spinor(s) Ψ = (χ L ,χ R ) containing two Weyl spinors χ L and χ R . The n → χχχ decay arises from 4-fermion effective operators of the form at nucleon level. This unusual operator involves the charge-conjugated field Ψ c = CΨ T such that the neutron decays into three χ particles and noχ anti-particles. In the relevant nonrelativistic limit M < ∼ m n /3 the decay rate is given by where A is the n ↔ χχχ decay amplitude. A minimal non-relativistic suppression |A | 2 nr ≡ (m n /Λ n ) 4 with ∼ 1 − 3M/m n necessarily arises if the decay is into N = 1 χ generation. This can be seen from non-relativistic quantum mechanics: if the χ momenta vanish, the Pauli principle demands the fully anti-symmetric product of three χ spinors (heavy-quark effective theory can be adapted to obtain a systematic expansion [51]). On the other hand, decays involving N > 1 generations of χ particles allow for ∼ 1. We thereby estimate A lower Λ χn ≈ 30 TeV is needed if ∼ 1 − 3M/m n .
The rate of the related visible decay channel with an extra γ produced from the neutron magnetic moment interaction is estimated as and is safely below the experimental bounds [17].
In the present model χ can be DM while preserving hydrogen stability: hydrogen decay is kinematically open only in the part of the parameter space with lower M , as shown in fig. 1. In such a case, the hydrogen decay rate can be estimated taking into account weak interactions of neutrons, where |ψ(0)| 2 ≈ α 3 m 3 e /π is the inverse atomic volume, and E ν < ∼ m p +m e −3M . Since this dominant hydrogen decay mode is fully invisible, the experimental bound on its rate is comparable to the inverse universe age and largely satisfied. A mildly stronger bound Γ H→χχχν e < ∼ 1/(10 14 yr) arises considering electron ep capture in the sun [52]. Adding an extra photon gives a subdominant but visible hydrogen decay mode. Its rate Γ H→χχχν e γ ≈ αE 2 γ 4πm 2 n Γ H→χχχν e ≈ 1 10 39 yr (12) is compatible with the bound from Borexino, Γ H→χχχν e γ < ∼ 1/(10 28 yr) [18,53].
Coming back to theory, the nucleon-level operator of eq. (7) can arise from 6-fermion quarklevel operators invariant under SM gauge interactions, such as where Q = (u, d) L , d R , u R are the SM quarks in Weyl notation and we omitted χ chiralities, χ ∼ χ L ∼χ R and Lorentz indices that can be contracted in multiple ways. Lattice computations of the nuclear matrix element 0|(ud) L,R d L,R |n = β L,R n find |β L,R | ≈ 0.014 GeV 3 ∼ Λ 3 QCD [54], and chiral perturbation theory allows to compute interactions with extra pions or photons [55].
In view of their large dimension 9, the operators in eq. (13) are strongly suppressed by Λ qχ . The neutron-level operator of eq. (7) is obtained with coefficient 1/Λ 2 χn = β L,R /Λ 5 χq , so that the neutron decay anomaly is reproduced for low Λ χq ∼ 30 GeV. Since Λ χq ∼ 30 GeV is below the weak scale, a mediator below the weak scale is needed. This can be achieved, compatibly with collider bounds, adding for example a neutral fermion n , singlet under all SM gauge interactions, coupled via dimension-6 operators as n χχχ/Λ 2 χn +n ddu/Λ 2 nn .
Baryon number is conserved assigning B n = 1 to n . In turn, the 4-fermion operators in eq. (14) can be mediated by renormalizable couplings of bosons, such as a W R or a lepto-quark, and by their dark-sector analogous. One gets Λ 5 n = M n Λ 2 χn Λ 2 nn , so that the 4-fermion operator n ddu involving SM particles can be suppressed by a Λ nn above the weak scale, up to a few TeV, while keeping M n > 1.5m n in order to avoid conflicting with neutron star bounds [36], and keeping the dark interaction Λ χn above the QCD scale in order to avoid conflicting with bounds on DM self-interactions, σ/M < ∼ 10 4 / GeV 3 [56,57].

Dark Matter cosmological abundance
To match the observed cosmological DM abundance, the χ particles must have abundance Y ≡ n/s = 0.44 eV/M = 1.4 10 −9 , where n is the χ +χ number density and s = 2π 3 d SM T 3 /45 is the entropy density of the thermal bath with d SM (T ) degrees of freedom at temperature T .
The neutron decay anomaly is reproduced for Λ χn v, so the interaction rates of χ particles at T ∼ m n are below electroweak rates and thereby below the Hubble rate: thermal freeze-out is not possible. 3 For the same reason, no equilibration of asymmetries happens at T < ∼ m n .
We thereby consider freeze-in, assuming vanishing initial χ abundance and computing the multiple contributions it receives from particle-physics processes. We avoid entering into modeldependent mediator issues, and explore the effect of the interaction with neutrons in eq. (7), for example by assuming that the reheating temperature is below the mediator masses. At leading order, neutrons and anti-neutrons that decay at T ∼ m n close to thermal equilibrium 4 contribute to the DM abundance as where R ∼ e −Λ QCD /m n corrects the freeze-in formula taking into account that neutrons only form after the QCD phase transition at T < ∼ Λ QCD . This Y is about 10 4 smaller than the cosmological DM abundance. Three different effects provide the needed enhancement.
1. First, thermal scatterings such as π 0 n ↔ χχχ and γn ↔ χχχ arise at higher orders in the QCD or QED coupling g, and avoid the non-relativistic suppression of n ↔ χχχ. This suppression was computed in eq. (9): two powers of 1 − 3M/m n ≈ few × 10 −3 arise from the phase space, and possibly one extra power arises from the squared amplitude. Thereby, at the relevant temperature T ∼ m n , the scattering rates are enhanced by a factor of order g 2 (T /E Be ) 2−3 ∼ 10 4−6 compared to the Γ n→χχχ rate. Indeed these finitetemperature scatterings can be partially accounted by a n ↔ χχχ rate enhanced by taking into account the thermal contribution to the neutron squared mass, m 2 n + O(g 2 T 2 ).
2. Second, n * → χχχ decays of the excited neutron n * with mass m n * ≈ 1.44 GeV and of other QCD resonances similarly avoid the non-relativistic suppression of n → χχχ. Dedicated lattice computations are needed to confirm that these resonances have matrix elements similar to those of the neutron. We expect that their effect is numerically similar to thermal scatterings at point 1.
3. Third, the processesχn ↔ χχ and χn ↔χχ have similar phase-space unsuppressed rates, and increase the number of DM particles that have already been produced, addinḡ χ anti-particles. Their effect is numerically comparable to those at previous points.
A precise computation is not possible since non-perturbative QCD effects determine the rates at point 2. Furthermore, larger but model-dependent contributions to freeze-in can arise at higher temperature, if the reheating temperature is higher than Λ QCD .

Dark Matter signals
Theχn ↔ χχ process at point 3 leads to unusual DM direct detection signals: it is kinematically open today and the non-relativisticχ and n in its initial state produce relativistic χ with energy E = 2m n /3. Theχn ↔ χχ cross section avoids the non-relativistic suppression of n → χχχ and is thereby estimated as This is below current bounds from direct detection experiments for the value of Λ χn motivated by the neutron decay anomaly. No nuclear enhancement arises, since the energies involved in χn ↔ χχ are higher than the tens of keV produced by usual DM-induced nuclear recoils. For the same reason, bigger experiments dedicated to more energetic signals produced by neutrinos and proton decay have better sensitivity. Following their practice, we convert eq. (16) into the event rate per neutron as having assumed an equal number of χ andχ. We thereby have the following signals: 1. DM scatterings that lead to the disappearance of a neutron into invisible DM. These scatterings effectively make ordinary matter radioactive with life-time given by eq. (17), since a neutron that disappears within a nucleus leaves a hole, triggering nuclear deexcitations and decays. For example DM scatterings convert water 16 that is thereby satisfied.
Concerning indirect detection signals, extra interactions not needed by the model might lead to χχ annihilations into pairs of SM particles.

Conclusions
Past literature showed that the interpretation of the neutron decay anomaly in terms of n → χγ decays into new free invisible particles χ with baryon number B χ = 1 implies their thermalization inside neutron stars, and that the modified equation of state is in contradiction with observed neutron stars.
We found that the neutron star issue is more general, affecting all possible models listed in table 1 where we considered generic χ charges under conserved baryon and lepton numbers, B χ and L χ . Independently of particle-physics details, the χ chemical potential in thermal equilibrium is µ χ = B χ µ n . Neutron star bounds are found to be compatible with the minimal possibility that χ carries the smallest baryon number B χ = 1/3, so that the neutron decay anomaly can be accounted by n → χχχ decays. The resulting modification of neutron stars, shown in fig. 2, is comparable with SM uncertainties. Fig. 3 shows that neutron stars get modified in a specific and mild enough way even in the presence of multiple generations of χ particles, thereby providing a signal for improved future observations and computations.
In section 4 we presented possible theories for n → χχχ, showing that associated processes are compatible with bounds from colliders, on n → χχχγ decays, and on hydrogen decays into χχχν e or χχχν e γ (kinematically open only in a part of parameter space shown in fig. 1).
We next found that χ can be Dark Matter. In section 5 we showed that the cosmological DM abundance can be possibly matched by freeze-in production of χ after the QCD phase transition at T ∼ m n . Three processes contribute: π 0 n ↔ χχχ and γn ↔ χχχ thermal scatterings; decays of excited neutrons; χn ↔χχ andχn ↔ χχ scatterings. In section 6 we found thatχn ↔ χχ provide unusual DM signals with interesting rates below present bounds. Such DM scatterings lead to neutron disappearance, that effectively makes matter (such as water, carbon, etc) radioactive, giving signatures similar to those of neutron invisible decays. Predictions are compatible with bounds from DM detection experiments and from neutrino experiments such as SNO and KamLand. Furthermore, the higher order related processes with extra visible photons,χn → χχγ andχn → χχπ 0 → χχγγ, are compatible with Super-Kamiokande bounds.