Neutron Star Constraints on Neutron Dark Decays

Abstract

Motivated by the neutron lifetime puzzle, it is proposed that neutrons may decay into new states yet to be observed. We review the neutron star constraints on dark fermions carrying unit baryon number with masses around 939 MeV, and discuss the interaction strengths required for the new particle. The possibility of neutrons decaying into three dark fermions is investigated. While up to six flavors of dark quarks with masses around 313 MeV can be compatible with massive pulsars, any such exotic states lighter than about 270 MeV are excluded by the existence of low-mass neutron stars around ∼$1.2M_{\odot }$. Light dark quarks in the allowed mass range may form a halo surrounding normal neutron stars. We discuss the potential observable signatures of the halo during binary neutron star mergers.

1. Introduction

Neutrons and protons are fundamental building blocks of the visible universe. Although not considered as elementary particles, they are windows to the world below femtometer scales that is challenging and costly to probe. Understanding their properties is essential in testing the Standard Model (SM) and could also be instrumental in searching for signatures of physics beyond the SM. For instance, the observed incredibly long lifetime of protons has put severe constraints on supersymmetric and grand unified theories that predict appreciable baryon number violating processes at low energies [1,2,3].

Unlike their electrically charged counterpart, neutrons are more challenging to study in terrestrial laboratories. This difficulty is especially highlighted in recent years by the neutron decay anomaly. It refers to a small yet distinctive discrepancy among the values of neutron lifetime ${\tau }_n$ measured with two classes of techniques [3,4]. In the so-called “beam” method, fluxes of beta-decay products, protons, and electrons are measured; in the “bottle” method, one continually monitors an ensemble of ultra-cold neutrons in a trap, and the lifetime is obtained by fitting the survival probability as a function of time.

Attempts by various groups yielded averaged values of ${\tau }_n^{\mathrm{bottle}}=878.4\pm 0.5\mathrm{s}$ [5,6,7,8,9,10] for the bottle method, and ${\tau }_n^{\mathrm{beam}}=888.0\pm 2.0\mathrm{s}$ [11,12,13] for beam measurements. This tension, on the order of ∼$10\mathrm{s}$, has already grown above the $4\sigma$ level [3,4]. Since the bottle method measures the total decay width, whereas, in the beam method, only the beta-decay width is measured, the shorter lifetime reported using the bottle method could indicate a hitherto unknown decay channel of neutrons [14].

On the theory side, current Standard Model calculations of the neutron lifetime do not reveal particular preferences for either the beam or the bottle results [14,15]. For instance, if one fix the CKM matrix element $|V_{ud}|$ to the Particle Data Group average [3], although a few measurements of the axial vector coupling inferred larger values of $g_A\gtrsim 1.275$ [16,17,18], favoring the bottle average, recent lattice calculations [19,20] and several early experiments [21,22] point to the lower side $g_A\lesssim 1.270$ predicting ${\tau }_n$ closer to the beam results.

As the neutron lifetime puzzle continues haunting the community, one shall not overlook where neutrons abound. Containing over $\sim {10}^{57}$ baryons, neutron stars (NSs) are gigantic nuclei rich in neutrons. With masses on the order of solar mass $M_{\odot }\approx 1.99\times {10}^{30}$ kg and radii around 10 km, neutron stars are compact objects of extreme densities. Their compactness $C=GM/R{c}^2\sim 0.1\text{–}0.3$ (a parameter related to the escape velocity) is large, only second to black holes ($C_{\mathrm{BH}}=0.5$), suggesting an environment of extreme gravity. Due to the strong gravitational field, baryons are compressed to densities $n_B$ a few times higher than the nuclear saturation density $n_0=0.16{\mathrm{fm}}^{-3}$, the typical density encountered near the center of nuclei. These conditions make neutron stars ideal laboratories to study new degrees of freedom that mix with neutrons.

The remainder of this paper is organized as follows. In Section 2, we briefly review proposals for neutron dark decays. Section 3 is a short introduction to the physics of neutron stars, with a focus on the equation of state and the neutron star mass–radius relations. In Section 4, we review existing neutron star constraints on dark neutrons [23,24,25,26,27], and discuss neutron decays into dark quarks in Section 5. Below, we have outlined the natural unit system in which $\hslash =G=c=1$ is adopted.

2. Dark Neutrons and the Neutron Decay Anomaly

We begin by examining the general features of neutron dark decay scenarios required to explain the puzzle surrounding ${\tau }_n$. New decay channels are possible at tree levels if the SM neutron couples to another spin one-half operator X, where X could be an intermediate state or is composed of the decay products. The final state may consist of several fermionic and possibly bosonic fields

$$n\to X_1{X}_2\dots ,$$

(1)

where at least one of them is not the content of the Standard Model. Based on kinematic considerations, the dark decay may proceed if the sum of the masses of final state particles does not exceed that of neutrons ${\sum }_i{m}_{X_i}<m_n$. Additionally, the incredibly long lifetime of protons suggests the proton decay process $p\to e^{+}{\nu }_e{X}_1{X}_2\dots$ must be forbidden. This can be achieved via kinematics when ${\sum }_i{m}_{X_i}>m_p-m_e=937.76\mathrm{MeV}$. A slightly stronger bound is given by the stability of nuclei, and the strongest of such requirements comes from the weakly bounded ${}^9$Be nucleus, which translates into ${\sum }_i{m}_{X_i}>937.90\mathrm{MeV}$ [14]. All in all, neutron dark decay scenarios can be viable solutions to the lifetime puzzle provided that the sum of final state masses lies in the narrow range

$$937.90\mathrm{MeV}<\sum _i{m}_{X_i}<939.57\mathrm{MeV}.$$

(2)

A number of models that can generate the operator $nX$ and predict Equation (1) have been proposed [14,23,26]. In this section, for concreteness, we shall focus on the possibility in which the final state consists of the SM photon $\gamma$ and a Dirac fermion $\chi$ carrying unit baryon number [14,23], i.e.,

$$n\to \chi \gamma .$$

(3)

According to Equation (2), the mass of this dark fermion $\chi$ is constrained to the interval

$$937.90\mathrm{MeV}<m_{\chi }<939.57\mathrm{MeV}.$$

(4)

We further note that, if $m_{\chi }<m_p+m_e=938.78\mathrm{MeV}$, $\chi$ is stable, as guaranteed by the conservation of baryon number and electric charge, and is a well-motivated dark matter candidate.

The dark decay Equation (3) may arise from the effective Lagrangian

$$\begin{array}{cc}{\mathcal{L}}_{\mathrm{eff}}=\hfill & \overline{n}\left(i\cancel{\partial }-m_n+{\mu }_n{\sigma }^{\mu \nu }{F}_{\mu \nu }\right)n\hfill \\ & +\overline{\chi }\left(i\cancel{\partial }-m_{\chi }\right)\chi -\delta \left(\overline{\chi }n+\overline{n}\chi \right),\hfill \end{array}$$

(5)

where ${\mu }_n=-1.91e/m_n=-0.3{\mathrm{GeV}}^{-1}$ is the neutron magnetic dipole moment. The $n-\chi$ coupling $\delta$ can be generated by the dimension-6 operator $\frac{1}{{\Lambda }^2}\overline{\chi }udd$. UV completions of Equation (5) suggest that the value of this coupling is $\delta$∼$\frac{0.01{\mathrm{GeV}}^3}{{\Lambda }^2}$ [28,29], where $\Lambda$ is the characteristic energy of the UV theory. Taking $\Lambda$ to be around the electroweak scale $\Lambda$∼$100\mathrm{GeV}$, we see that the $n-\chi$ coupling is expected to be small $\delta$∼${10}^{-3}\mathrm{MeV}$. In particular, it is much smaller than the mass difference $\delta \ll \Delta m\equiv m_n-m_{\chi }$.

Diagonalizing the mass terms in Equation (5) is straightforward, and it leads to a mixing between n and $\chi$. To leading order in the mixing angle $\theta =\delta /\Delta m$ and the mass eigenstates are obtained by taking $n\to n+\theta \chi$, and $\chi \to \chi -\theta n$. In the mass eigenbasis, the operator responsible for the dark decay can be identified as

$${\mu }_n\theta \overline{\chi }{\sigma }^{\mu \nu }n{F}_{\mu \nu }.$$

(6)

It predicts $n\to \chi \gamma$ with a partial width

$${\Gamma }_{n\to \chi \gamma }=\frac{{\mu }_n^2{\theta }^2{m}_n^3}{16\pi }{\left(1-\frac{m_{\chi }^2}{m_n^2}\right)}^3\simeq \frac{{\mu }_n^2{\theta }^2\Delta m^3}{2\pi }.$$

(7)

Given a total width of ${\Gamma }_n=1/{\tau }_n^{\mathrm{bottle}}={\left(879.6\mathrm{s}\right)}^{-1}$, the branching ratio for the dark decay into $\chi \gamma$ is

$${\mathrm{Br}}_{n\to \chi \gamma }=0.01{\left(\frac{\Delta m}{1\mathrm{MeV}}\right)}^3{\left(\frac{\theta }{7\times {10}^{-10}}\right)}^2.$$

(8)

A search for photons in the expected energy range was conducted shortly after the publication of reference [14], where the branching ratio is constrained to be less than $\sim {10}^{-3}$ for $E_{\gamma }>m_n-m_p-m_e\approx 0.78\mathrm{MeV}$ [30], thereby ruling out this $\chi$ as a dark matter candidate, or that $n\to \chi \gamma$ as a complete solution to the neutron lifetime puzzle if $\chi$ is stable. We shall demonstrate that neutron star constraints on neutron dark decays are complementary to and stronger than this terrestrial bound, as they can probe branching ratios much less than required to explain the lifetime puzzle and do not require visible states in the decay product.

Additional models are proposed in which neutrons decay into two dark states [14,27],

$$n\to \chi \varphi ,$$

(9)

where $\chi$ is a fermion and $\varphi$ is a boson.

In reference [27], baryon number $B=1$ is assigned to the fermion $\chi$, whereas in reference [14], the scalar $\varphi$ carries unit baryon number. We will briefly review the neutron star constraints for this scenario in Section 4.4, and comment on future directions. In Section 5, the possibility of neutron decaying into three dark fermions $n\to {\chi }_i{\chi }_j{\chi }_k$ is investigated [31], and its implications for neutron stars examined.

3. Neutron Stars

Here, we briefly review the physics of neutron stars relevant for the discussion on neutron dark decay.

3.1. The Mass–Radius Relation

The hydrostatic equilibrium of spherically symmetric configurations in general relativity was first solved by Tolman and by Oppenheimer and Volkov almost a century ago [32,33]. The resulting equations that govern compact star structures, dubbed the TOV equations, are

$$\begin{array}{cc}\frac{\mathrm{d}P}{\mathrm{d}r}\hfill & =-\frac{M\mathcal{E}}{r}\left[1+\frac{P}{\mathcal{E}}\right]\left[1+\frac{4\pi r^{2}p}{M}\right]{\left[1-\frac{2M}{r}\right]}^{-1},\hfill \\ \frac{\mathrm{d}M}{\mathrm{d}r}\hfill & =4\pi r^2\mathcal{E}.\hfill \end{array}$$

(10)

Above, M is the enclosed mass, P is the pressure, and $\mathcal{E}$ the energy density of the stellar matter.

Equation (10) is written in a suggestive form to reveal its resemblance with its Newtonian counterpart. The weak field, non-relativistic limit is obtained upon dropping all three terms in the bracket in the first line. In this limit, the equations can be understood through analyzing the forces acting on an infinitesimally thin slab of the star. Equilibrium is achieved when the difference in pressure across the slab balances the gravitational force per unit area sourced by this thin slab. Even though the relativistic corrections in Equation (10) can be as large as factors of a few, the heuristic “pressure gradients balance gravity” argument remains useful in construing global structures of relativistic stars.

Once an equation of states (EOS), a relation between the pressure and energy density $\mathcal{E}\left(P\right)$, is specified, the TOV equations can be solved numerically as follows. For given boundary condition $P(r=0)=P_c$, specified at the center of the star $r=0$ where $M\left(r\right)=0$, Equation (10) can be integrated all the way up to the surface $r=R$, defined by $P(r=R)=0$. The total stellar mass is then found as $M=M(r=R)$. A family of parametric solutions $R\left(P_c\right)$ and $M\left(P_c\right)$ can thus be generated by varying the central pressure $P_c$, from which the mass–radius relation is obtained. Figure 1 shows a few such curves. The EOSs that underlie these curves are discussed below.

Ubiquitous to all mass–radius relations is the existence of maximum masses, beyond which points increasing the central energy density will only make lighter stars. That gravitational binding energy grows faster than the baryonic contribution signifies gravitational instability. In other words, the stars on the branch beyond the limiting mass are unstable, and are not shown in Figure 1.

Known as the TOV limit, the value of the maximally attainable neutron star mass depends on the details of the underlying EOS. Based on the “pressure gradients balance gravity” argument discussed earlier, one could expect that EOSs predicting higher pressure at given energy densities are more likely to support heavier neutron stars. An EOS that has higher pressure at given energy densities is usually referred to as “stiff”, and those with smaller pressure are typically called “soft”. The correlation between neutron star maximum mass and pressure of the matter content suggests that observations of massive neutron stars could provide valuable insights into the microscopic physics governing the dense stellar medium. This avenue of research has been explored extensively in the nuclear physics community in studying neutron-rich matter [34]. Here, instead, we aim to use this connection to assess the possibilities of dark contents in neutron stars.

To date, three measurements of massive pulsars have been made with high confidence and good accuracy. These are $1.97\pm 0.04M_{\odot }$ for J1614-2230 [35], $2.01\pm 0.04M_{\odot }$ for J0348+0432 [36], and $2.08\pm 0.7M_{\odot }$ for J0740+6620 [37]. All of the pulsars being measured reside in binary systems. The rather accurate mass measurements are made possible thanks to relativistic orbital effects which are used to break degeneracies among the Keplerian elements. Below, we shall see these observations are powerful in setting limits on dark particles that mix with neutrons.

3.2. Equation of State (EOS)

It is often stated in popular science books and introductory-level texts that neutron stars are supported by the degeneracy pressure of neutrons. This is somewhat misleading. In fact, it was known in the very early days by Oppenheimer and Volkov [32] that, if the degeneracy pressure is the only source of stabilization, neutron stars cannot exist at all.

Although the ideal fermi gas is not appropriate to describe dense nuclear matter, it will play a central role in the discussion below surrounding neutron dark decays. Here, we present a brief review. The EOS for arbitrary relativistic ideal Fermi gas (with degeneracy 2 due to spin) is given by

$$\begin{array}{cc}\mathcal{E}\hfill & =\frac{1}{{\pi }^2}{\int }_0^{k_F}\sqrt{k^2+m^2}{k}^2\mathrm{d}k,\hfill \\ P\hfill & =-\frac{\mathrm{d}(\mathcal{E}/n)}{\mathrm{d}(1/n)}=\frac{1}{3{\pi }^2}{\int }_0^{k_F}\frac{3k^4\mathrm{d}k}{\sqrt{k^2+m^2}},\hfill \end{array}$$

(11)

where $n=k_F^3/\left(3{\pi }^2\right)$ is the fermion number density, $k_F$ is the Fermi momentum, and m denotes the mass of the fermion. These integrals can be evaluated analytically, leading to closed-form expressions

$$\begin{array}{cc}\mathcal{E}\hfill & =\frac{m^4}{8\pi }\left[x\sqrt{x^2+1}\left(2x^2+1\right)-{sinh}^{-1}\left(x\right)\right],\hfill \\ P\hfill & =\frac{m^4}{24\pi }\left[x\left(2x^2-3\right)\sqrt{x^2+1}+3{sinh}^{-1}\left(x\right)\right],\hfill \end{array}$$

(12)

where we have defined the dimensionless quantity $x=k_F/m$.

The resulting mass–radius relation is shown in green in Figure 1. Clearly, the predicted maximum neutron star mass sits around $0.7M_{\odot }$, far below even the Chandrasekhar limit. Such low masses are not only incompatible with observations of heavy pulsars but also in stark contradiction with the neutron star formation mechanism. The main takeaway here is that the degeneracy pressure alone is not adequate to support any realistic neutron stars, and that the short-range nuclear repulsive forces are crucial in determining even the global structure of neutron stars.

As discussed in the previous section, the Fermi gas EOS cannot support heavy neutron stars because it is too soft, i.e., it predicts (much) lower pressure at given energy densities compared to more realistic models (to be discussed below). This is evident in Figure 2. A useful measure of the “stiffness” of an EOS is the isentropic speed of sound. It is related to the slope of the EOS as

$$c_s\equiv \sqrt{C_s}=\sqrt{\frac{\mathrm{d}P}{\mathrm{d}\mathcal{E}}}.$$

(13)

For the ideal Fermi gas EOS, the sound speed squared takes the simple form

$$C_s=\frac{1}{3}{\left[1+{\displaystyle \frac{m^2}{{\left(3{\pi }^{2}n\right)}^{2/3}}}\right]}^{-1}.$$

(14)

It gradually rises with density and asymptotes to the ultra-relativistic limit $1/3$. Note that, depending on the mass of the fermion, the asymptotic value $1/3$ may be achieved at densities much higher than or comparable to values (around a few times $n_0$) relevant for the interior of neutron stars (see Figure 3). This is the key motivation for Section 5.

We shall now turn to more realistic nuclear EOSs. Although it remains a daunting task to describe neutron star matter, especially as one gets closer to the central regions of the star, much progress has been made to improve our understanding of the EOS below about twice the nuclear saturation density $n_B\le 2n_0$. In particular, the chiral effective field theory ($\chi \mathrm{EFT}$) [38,39,40,41,42] developed over the past three decades has produced fruitful predictions for neutron-rich matter in this density range [43,44]. It is an effective low-energy description of QCD in terms of explicit nucleonic degrees of freedom. All the operators consistent with the chiral symmetry are taken into account, and power counting schemes are adopted to organize them into a series in nucleon momenta. This not only enables systematically improvable calculations but also provides uncertainty estimates through order-by-order comparisons. For instance, the recent many-body calculations based on chiral potentials up to the next-to-next-to-next-to-leading order (N3LO) suggest that the pressure of pure neutron matter is about 19 MeV${\mathrm{fm}}^{-3}$ at $n_B=2n_0$ with an uncertainty of about $30\%$ at the $1\sigma$ level [43,44,45].

In this work, we adopt a family of EOSs for nuclear matter in beta equilibrium [45,46,47] based on the aforementioned state-of-art N3LO $\chi \mathrm{EFT}$ calculations. This model accounts for both the central values and the correlated uncertainties through a simple and accurate parameterization [45]. For baryon number densities above around $1.5n_0$, we shift to an EOS with a constant speed of sound, i.e.,

$$P=C_s\mathcal{E}+P_0,$$

(15)

where the constant $P_0$ is chosen so that the pressure is continuous. Since causality requires the speed of sound to be no greater than one, the limiting case $C_s=C_{\max }=1$ is the stiffest possible nuclear EOS (assuming $\chi \mathrm{EFT}$ calculations are valid up to $n_B=1.5n_0$). This is shown in Figure 2 as the blue curve labelled “stiff”. Its prediction for the TOV limit, around $3.5M_{\odot }$, places an upper bound on how heavy neutron stars can be. Below, when neutron dark decay channels are turned on, we shall show that new fermionic states carrying baryon number drastically soften the EOS, which leads to lowered neutron star maximum mass. This extreme EOS will provide the most conservative yet robust bounds on such dark states, as the results are insensitive to the uncertainties associated with the poorly understood high-density QCD phase.

Figure 1 and Figure 2 show the EOSs discussed so far and their predictions for the mass-radius relations. The EOS labelled APR, calculated by Akmal, Pandharipande, and Ravenhal [48], has been widely used in the literature to describe neutron stars and it is included here as a reference. The curve labeled DD2 is based on the relativistic mean field calculation reported in reference [49]. All of these SM EOSs support two-solar-mass neutron stars, and will serve as the baseline when dark decay channels are turned on.

4. Neutron Star Constraints on Dark Neutrons

Given the long lifetime of observed neutron stars $t_{\mathrm{NS}}\approx {10}^6\text{–}{10}^8$ years, one could expect that any dark decay scenario aiming to explain the neutron lifetime puzzle, which in vacuum has a characteristic time scale

$${\tau }_{n\to \chi \dots }=1/\left(\frac{1}{{\tau }_n^{\mathrm{bottle}}}-\frac{1}{{\tau }_n^{\mathrm{beam}}}\right)\approx 0.9\mathrm{day}$$

(16)

will be efficient in bringing the decay products into equilibrium with the star. This general expectation is indeed true, despite complications due to finite density effects. In a dense baryon environment, each neutron feels strong interactions with a large number of nucleons in the background, and thus the excitation spectrum of neutrons is modified. This leads to a suppressed in-medium mixing angle. For the model described in Equation (5), the n-$\chi$ mixing at finite density is given by

$$\tilde{\theta }=\frac{\delta }{\sqrt{{\widetilde{\Delta m}}^2+{\Sigma }_i^2}}$$

(17)

where $\widetilde{\Delta m}=\Delta m+{\Sigma }_r$, and ${\Sigma }_r$ and ${\Sigma }_i$ are the real and imaginary parts of the neutron dispersion relation, respectively. Since ${\Sigma }_r$ and ${\Sigma }_i$ are on the order of 10–100 MeV in neutron stars [50], the suppression on the mixing angle $\tilde{\theta }/\theta$ is expected to be around 0.01–0.1. The enhanced energy level splitting $\widetilde{\Delta m}$, on the other hand, enlarges the available phase space for $n\to \chi \dots$ in the star. Since the production rate in Equation (7) is proportional to the mixing angle squared and the energy level difference cubed, it is safe to assume that $\chi$ will come into equilibrium on a timescale much shorter than the typical age of neutron stars as long as $\delta \gtrsim {10}^{-17}$ GeV, much less than the value required to explain the neutron lifetime puzzle. In other words, even if the decay anomaly is resolved within the nuclear experiment community, neutron star considerations remain relevant as they are sensitive to $n-\chi$ mixings several orders of magnitude less than being probed at terrestrial laboratories.

Once in equilibrium, the large baryon chemical potential will source a considerable population of dark particles carrying baryon number. Depending on the dark decay products, the implications of those exotic states on neutron star structures will differ. Below we will examine the impacts of three types of dark decay modes on neutron star mass–radius relations.

4.1. Non-Interacting Dark Neutrons

In this section, we consider a Dirac fermion $\chi$ that carries unit baryon number. $\chi$ could be the fermion in Equation (5), or it can be the final fermionic state in Equation (9). We note that the search for the decay $n\to \chi \gamma$ already set a limit on this particular channel. For $937.90\mathrm{MeV}<m_{\chi }<938.78\mathrm{MeV}$, the branching ratio is constrained to be less than ${10}^{-3}$ [30], a level insufficient to explain the neutron decay anomaly on its own. Neutron stars can place stronger constraints on such a dark fermion with a mass close to that of nucleons, independent of other accompanying decay products [23,24,25].

As discussed earlier, in equilibrium, the dark neutron $\chi$ will be sourced by the baryon chemical potential ${\mu }_B$. In the absence of strong repulsive interactions, as is the case for models Equation (5) and in reference [14], the $\chi$’s are described by the Fermi gas EOS Equation (12). The chemical potential for this dark neutron gas is given by

$${\mu }_{\chi }=\sqrt{k_{\mathrm{F}}^2+m_{\chi }^2},\mathrm{where}{k}_F={\left(3{\pi }^2{n}_{\chi }\right)}^{1/3}.$$

(18)

The equilibrium dark neutron number density $n_{\chi }$ can be obtained by solving the condition

$${\mu }_{\chi }={\mu }_B=(P_{\mathrm{nucl}}+{\mathcal{E}}_{\mathrm{nucl}})/n_{\mathrm{nucl}}$$

(19)

for given nuclear EOSs in beta-equilibrium where $P_{\mathrm{nucl}}$, ${\mathcal{E}}_{\mathrm{nucl}}$, and $n_{\mathrm{nucl}}$ denote the pressure, energy density, and number density, respectively. The hybrid EOS appropriate for neutron stars admixed with dark neutron $\chi$’s then follows as

$$\begin{array}{cc}{P}_{\mathrm{tot}}\hfill & =P_{\mathrm{nucl}}+P_{\chi },\hfill \end{array}$$

(20)

$$\begin{array}{cc}{\mathcal{E}}_{\mathrm{tot}}\hfill & ={\mathcal{E}}_{\mathrm{nucl}}+{\mathcal{E}}_{\chi }.\hfill \end{array}$$

(21)

The reduction in pressure and, consequently, the maximum mass due to $\chi$’s is striking. Even for the maximally stiff EOS, the presence of non-interacting dark neutrons reduces the maximum mass to values well below observed neutron star masses. Thus, a dark neutron with $m_{\chi }\simeq m_n$ and weak interactions is robustly excluded [23,24,25]. The drastic softening of the EOS and the reduction to the TOV limit is universal across all nuclear models. This is because, in the absence of repulsive interaction among $\chi$’s it is energetically favorable to convert neutrons into $\chi$’s. Indeed, for the hybrid star based on the stiffest possible nuclear EOS, there are almost eight times more $\chi$’s than neutrons at the center of the heaviest star. Consequently, the hybrid EOS is dominated by that of the weakly interacting $\chi$’s, leading to TOV limits decreased to around $0.7M_{\odot }$, the value predicted by the Fermi gas EOS. The takeaway here, again, is that microscopic interactions are critical in supporting neutron stars.

4.2. Mirror Neutrons

A simple scenario in which dark neutrons may interact substantially is the possibility that they are mirror copies of the SM neutrons [51,52,53]. Here, we will assume that the mirror neutron is almost identical to its SM counterpart, except its mass $m_{\chi }$ is slightly less than $m_n$, and falls in the narrow range in Equation (4). In this case, we have $P_{\chi }=P_n$ and ${\mathcal{E}}_{\chi }={\mathcal{E}}_n$. We may further ignore the small fraction of protons in the star and arrive at $P_{\chi }=P_{\mathrm{nucl}}$ and ${\mathcal{E}}_{\chi }={\mathcal{E}}_{\mathrm{nucl}}$. The resulting mass–radius relations are shown in Figure 4. Here, the reduction to the neutron star maximum mass persists, but the effects are less prominent since $n_{\chi }=n_n$ and since short-range repulsions experienced by $\chi$’s provide additional stabilization. While hybrid stars constructed based on the stiffest possible nuclear EOS can still reach around $2.4M_{\odot }$, those assuming APR EOS as the underlying nuclear model can no longer satisfy the massive pulsar bound. We find that the maximum speed of sound squared for neutron star matter needs to be greater than about 0.6 ($C_s\gtrsim 0.6$ in Equation (15)) for mirror neutrons to be a viable neutron decay product.

4.3. Self-Interacting Dark Neutrons

Another possibility to obtain sizable repulsion among dark neutrons is by charging $\chi$ under a global $U{\left(1\right)}_{A^{\prime }}$ [23,26]. This leads to a class of well-motivated models that have rich phenomenology if the new gauge boson mixes with the Standard Model photon $A^{\mu }$ [54]. The effective Lagrangian involving the new gauge boson $A^{\prime }$ of mass $m_V$ that complements Equation (5) is given by

$${\mathcal{L}}_{A^{\prime }}=-\frac{1}{4}{F}^{\prime \mu \nu }{F}_{\mu \nu }^{\prime }-\frac{1}{2}{m}_V^2{A}^{\prime \mu }{A}_{\mu }^{\prime }-\frac{ϵ}{2}{F}^{\mu \nu }{F}_{\mu \nu }^{\prime }+g\overline{\chi }{\cancel{A}}^{\prime }\chi .$$

(22)

Above, $F^{\prime \mu \nu }={\partial }^{\mu }{A}^{\prime \nu }-{\partial }^{\nu }{A}^{\prime \mu }$, and $ϵ$ is the kinetic mixing between the SM photon and the dark photon $A^{\prime }$. This mixing could open up another decay channel $n\to \chi A^{\prime }$ for sufficiently light $A^{\prime }$. The coupling g is responsible for the repulsive interactions among $\chi$’s mediated by the vector boson $A^{\prime }$. At the relativistic mean field level, it leads to the following contributions to the EOS in addition to those of the Fermi gas (e.g., references [55,56])

$$\begin{array}{cc}\Delta P_{\chi }\hfill & =\Delta \mathcal{E}\chi =\frac{1}{2}\frac{g^2}{m_V^2}{n}_{\chi }^2,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\Delta {\mu }_{\chi }\hfill & =\frac{g^2}{m_V^2}{n}_{\chi }.\hfill \end{array}$$

(24)

Once again, besides making the EOS of $\chi$’s stiffer, the interaction also raises the cost of converting n’s to $\chi$’s, as can be seen from Equation (23). This helps to keep the population of $\chi$’s in check. A back-of-the-envelope estimate suggests that mediator mass on the order of $m_V$ ∼ 100 MeV can achieve this by matching the self-energy of neutrons:

$$\Delta {\mu }_{\chi }\approx 10\mathrm{MeV}\left(\frac{g^2}{0.1}\right){\left(\frac{100\mathrm{MeV}}{m_V}\right)}^2\left(\frac{n_{\chi }}{n_0}\right).$$

Figure 5 shows the required dark neutron interaction strengths as functions of TOV limits for the stiff and the APR nuclear EOSs. Assuming the coupling $g=1$, the observed massive pulsar would require a mediator lighter than about 100 MeV. If a three-solar-mass pulsar is found in the future, the mediator mass will not be allowed to exceed about 30 MeV. Such light gauge bosons may be subject to constraints from CMB, BBN, stellar cooling, etc., [26] and could have interesting consequences in other contexts, for example, explaining the DM small-scale structure puzzles [57].

4.4. Bosonic Dark Decay Products

In reference [27], it is proposed that appreciable repulsion may arise between neutrons and the dark fermion $\chi$ in Equation (9), and can therefore evade the neutron star maximum mass bound discussed above. However, we note that the boson $\varphi$ in the decay final state will thermalize with the $\chi$’s and neutrons and, therefore, this scenario is subject to constraints imposed by black hole formation. The general idea is that, in the absence of the Pauli exclusion principle, the bosons will settle into a small sphere near the center of neutron stars once thermalized. The characteristic sizes of this boson cloud

$$R_{\varphi }\sim 10\mathrm{m}{\left(\frac{T}{{10}^{8}K}\right)}^{1/2}{\left(\frac{m_{\varphi }}{\mathrm{GeV}}\right)}^{1/2}$$

are expected to be small. Without strong repulsions, the gravitational collapse of the boson gas to a black hole is inevitable since the critical mass is tiny

$$M_{\mathrm{crit}}\simeq m_{\mathrm{pl}}^2/m_{\varphi }\approx \left(\frac{\mathrm{GeV}}{m_{\varphi }}\right){10}^{-19}{M}_{\odot }.$$

Such consideration has already been used to place constraints elsewhere, such as WIMP capture by neutron stars [59] and dark bosons produced in supernovae [60]. In the context of neutron dark decay, it is expected to impose even stronger bounds given the large number of bosons produced in neutron dark decay. Furthermore, for sufficiently strong boson self-interactions, the $n-\chi$ repulsion may not be required to stabilize heavy neutron stars.

It is interesting to note that, if the scalar carries baryon number, as is the case in reference [14], and that their self-repulsion is adequate to evade the black hole formation bound, a new baryon number superfluid phase could emerge inside NSs. For a similar scenario that gives rise to a lepton number superfluid, see [60]. If, on the other hand, the scalars are weakly interacting but are light (sub-MeV), they could leave the star, carrying away considerable energies. This extra cooling channel will impact the early evolutions of neutron stars, a scenario potentially constrained by the cooling trajectories of young pulsars [61]. A detailed study on bosonic dark decay byproducts is underway and will be reported in future work.

5. Neutron Star Constraints on Dark Quarks

Yet another possibility is that neutrons may decay into three dark fermions $n\to {\chi }_1{\chi }_2{\chi }_3$, each carrying baryon number 1/3. The fermions in the final state can be identical ${\chi }_1={\chi }_2={\chi }_3$ or distinct. The only requirement is that the sum of their masses lies in the range Equation (2) in order to explain the neutron decay anomaly while, at the same time, to not destabilize ${}^9$Be. Due to their resemblance to quarks in the Standard Model, we refer to these fermions as dark quarks.

5.1. Equal Masses

We first study the simplest cases where all the dark quarks share the same mass $m_{{\chi }_1}=m_{{\chi }_2}=m_{{\chi }_3}\approx m_n/3$, while still allowing ${\chi }_{1,2,3}$ to be distinct states. Implementing the dark quark admixed EOS is straightforward by noting that the dark quark chemical potential is only 1/3 of that of baryons. Taking ${\mu }_{\chi }={\mu }_B/3$ in Equation (19), we obtain

$${\mu }_{\chi }={\mu }_B/3=(P_{\mathrm{nucl}}+{\mathcal{E}}_{\mathrm{nucl}})/\left(3n_{\mathrm{nucl}}\right).$$

(25)

Equation (18) can then be used to calculate the equilibrium dark quark densities, from which Equation (12) yields the EOSs for each flavor.

$$\begin{array}{cc}{P}_{\mathrm{tot}}\hfill & =P_{\mathrm{nucl}}+{\displaystyle \sum _{i=1}^{N_f}}{P}_{{\chi }_i},\hfill \end{array}$$

(26)

$$\begin{array}{cc}{\mathcal{E}}_{\mathrm{tot}}\hfill & ={\mathcal{E}}_{\mathrm{nucl}}+{\displaystyle \sum _{i=1}^{N_f}}{\mathcal{E}}_{{\chi }_i},\hfill \end{array}$$

(27)

where the subscript i is the dark quark flavor index.

The resulting mass–radius relations are shown in Figure 6. While the APR EOS does not appear to support the dark quark hypothesis, for the maximally stiff nuclear EOS, dark quarks are compatible with current observations unless there are more than six species with masses equal to $m_n/3\approx 313$ MeV. Note that, although the final state in this scenario can only accommodate no more than three species, if additional flavors of light dark quarks exist (see below), they will be sourced by the baryon chemical potential and become relevant. Since the baryon chemical potential in the most massive neutron stars can reach ${\mu }_B\simeq 2.3$ GeV, any dark quark carrying baryon number $1/3$ can populate neutron stars if it is lighter than about 800 MeV.

The key insight into why weakly interacting dark quarks work is by noting the $m_{\chi }^2$ dependence in the speed of sound squared for the Fermi gas EOS Equation (14). Reducing the fermion mass enables the asymptotic value $C_s=1/3$ to be approached faster at lower densities. Figure 3 shows the speed of sound for selected values of $m_{\chi }$. At $m_{\chi }\approx 300$ MeV, the dark quark matter is endowed with $C_s$ a few times higher than that of the weakly interacting dark baryon gas at densities above $\sim n_0$. We therefore expect stiffer hybrid EOSs for dark quarks compared to dark neutrons. Furthermore, the zero-temperature thermodynamic identity $C_s^{-1}=\mathrm{d}logn/\mathrm{d}log\mu$ suggests that, at fixed chemical potentials, higher sound speeds imply lower densities. Hence, in the dense background of SM baryons that source ${\mu }_B$, the population of dark particles would be more restrained if they are considerably lighter than dark neutrons. Both of these factors allow dark quarks to support massive pulsars, even in the absence of strong repulsion.

5.2. Unequal Masses

One may further relax the assumption of equal dark quark masses. This would imply the existence of new states lighter than about $m_n/3$ carrying baryon number 1/3. In the absence of additional symmetries, the lightest dark quark(s) would be the ground state in vacuum into which heavier flavors decay on a time scale comparable to Equation (16). The equilibrium configuration in a dense baryon background is specified by the baryon chemical potential ${\mu }_{{\chi }_i}={\mu }_{{\chi }_j}={\mu }_{{\chi }_k}={\mu }_B/3$, as in the previous section. Unlike the equal mass case, however, additional considerations on neutron star stability arise here because the dark quarks are no longer contained within the SM baryon surface.

The outer layer of a cold neutron star, known as the outer crust, is a solid primarily composed of nuclei and electrons [62,63]. With increasing depth, the nuclei become increasingly neutron-rich. The top layer of the outer crust is mostly ${}^{56}\mathrm{Fe}$ and is characterized by a vanishing density and pressure (ignoring the atmosphere). Since the binding energy of ${}^{56}\mathrm{Fe}$ is about $8\mathrm{MeV}$ per baryon, the dark quark chemical potential there is ${\mu }_{\mathrm{threshold}}\approx (m_n-8\mathrm{MeV})/3=310\mathrm{MeV}$. According to Equations (18) and (25), this chemical potential would source a considerable population of dark quarks if their mass is less than this threshold value:

$$n_{\chi }\left({\mu }_{\mathrm{threshold}}\right)\approx \frac{{[{(m_n/3)}^2-m_{\chi }^2]}^{\frac{3}{2}}}{3{\pi }^2}\approx n_0{\left[1-{\left(\frac{3m_{\chi }}{m_n}\right)}^2\right]}^{\frac{3}{2}}$$

(28)

where the approximation ${\mu }_{\mathrm{threshold}}\approx m_n/3$ is implied. This nonzero density of dark quarks at the iron surface implies a finite pressure that would push the dark quark gas outward until a gradient of pressure is established outside, forming a halo surrounding the SM neutron star.

We note that dark quarks do not “drip” out of isolated stable nuclei. This is guaranteed by Equation (2). Their leakage outside neutron stars is a gravitational effect, as baryon numbers are redistributed across the star to minimize the total energy. The time scale associated with (re-)establishing hydrostatic equilibrium can be estimated as $\Delta t$∼$\sqrt{R^3/M}$ which, for typical neutron stars, is on the order of milliseconds. Since it is much shorter than the $n-\chi$ conversion time scale in Equation (16), the hydrostatic equilibrium is always maintained as the star readjusts itself. Once the chemical equilibrium is achieved, the true ground state is characterized by a constant ${\mu }_B$. In the presence of gravitational fields, this constancy condition reads (e.g., reference [64])

$${\mu }_B\left(r\right)\sqrt{g_{tt}\left(r\right)}=\mathrm{constant}$$

(29)

where ${\mu }_B\left(r\right)$ is the chemical potential in the local frame at radial coordinate r, and $g_{tt}\left(r\right)=e^{\nu \left(r\right)}$ is the temporal component of the metric tensor. This metric function is given by the equation

$$\frac{\mathrm{d}\nu }{\mathrm{d}r}=-\frac{\mathrm{d}P/\mathrm{d}r}{P+\mathcal{E}}$$

(30)

which can be solved alongside Equation (10) with the boundary condition $\nu (r=R)=log(1-\frac{2M}{R})$, i.e., matching the Schwarzschild metric at the surface of the star. In practice, this boundary condition may be imposed after integrating Equation (30) using an arbitrary value at the center of the star $\nu (r=0)$ and noting that Equation (30) is shift-invariant. The metric function for a $1.4M_{\odot }$ hybrid star is shown in Figure 7 as a solid black line.

To determine the structure of hybrid stars in the presence of dark quark halos, one may adopt the iterative approach based on Equation (29) described in reference [65]. Alternatively, simplifications are possible by noting that the dark quarks are only charged under the baryon number, so the equilibrium state of matter is uniquely specified by ${\mu }_B$. The simplified procedure is as follows. We first obtain the number density $n_{{\chi }_i}\left({\mu }_B\right)$ for each flavor down to baryon chemical potential ${\mu }_B=m_{{\chi }_i}/3$ via Equations (18) and (26), and compute the EOSs for every species according to Equation (12). The full EOS is then obtained by summing all contributions due to SM baryons and $\chi$’s via Equation (26), which, once supplemented to the TOV equation, would yield the ground state of dark quark admixed neutron stars.

Figure 7 shows the profile of an equilibrium configuration based on the APR EOS, assuming $N_f=2$ and that the mass of the lightest state is $m_{{\chi }_1}=280\mathrm{MeV}$. The mass of the other flavor can either be (1) $m_2=m_n-2m_1$ if $n\to {\chi }_1{\chi }_1{\chi }_2$ or (2) $m_2=(m_n-m_1)/2$ if $n\to {\chi }_1{\chi }_2{\chi }_2$, where, for brevity, we have defined $m_i\equiv m_{{\chi }_i}$. The resulting stellar structure is quite similar for these two cases and we picked the second possibility in Figure 7, since ${\chi }_2$ would populate the star in greater abundance thanks to the lower threshold chemical potential ${\mu }_{\mathrm{threshold}}=3m_2$. Although the ${\chi }_1$ halo extends to a large radius of about 18 km, the surface of SM baryons which determines thermal X-ray emissions sits at a radius of about 12 km, consistent with astrophysical observations (Spins are ignored. The hybrid stars most likely experience differential rotations due to the weak coupling between the halo and SM baryons, and the halo would remain intact in the presence of rapidly spinning SM cores. We defer discussions pertaining non-zero angular momentum to future work). And despite its large size, the halo only accounts for a small fraction of the total baryon numbers in the star. The baryonic mass of SM degrees of freedom is about $1.5M_{\odot }$, whereas for ${\chi }_2$, the share is close to $0.01M_{\odot }$ and for ${\chi }_1$ is around $0.07M_{\odot }$. Of the latter, just about $20\%$ resides in the halo.

Ignoring the thermal effects, the initial configuration can be found by searching for the cold neutron star composed solely of SM baryons that contain the same total baryon number as the hybrid one. Profiles of this normal neutron star are depicted as dotted lines in Figure 7. Upon comparing the solid and the dotted curves, we see that the effects of light dark quarks are mainly twofold. On one hand, the presence of dark quarks inside the iron surface softens the EOS, which leads to a small contraction of the SM sphere and a slight increase in the central density, thereby deepening the gravitational potential (shown in black) near the center and raising the gravitational binding energy by a few percent. This is similar to scenarios in previous sections and the upshot is mild reductions in the maximum attainable masses of hybrid stars. On the other hand, the light state ${\chi }_1$ provides a means to transport baryon numbers outside the iron surface. As noted earlier, the baryonic mass of the halo is about $0.01M_{\odot }$ and has negligible impacts on the spacetime geometry. Indeed, the $tt$ component of the metric tensor (solid black line) is almost identical to that of the Schwarzschild spacetime geometry $g_{tt}\left(r\right)=1-2M/r$ outside the SM portion of the star (dotted black line, very close to and hard to distinguish from the solid black). However, since number densities at the iron surface in Equation (28) can be quite large for very light dark quarks, the leakage of baryon numbers would be so severe that SM components do not survive. We now examine this possibility and demonstrate that low-mass neutron stars put stringent lower limits on dark quark masses.

Figure 8 shows the mass–radius relations for a few hybrid EOS models accommodating dark quark halos. As expected, maximum masses are moderately reduced, an effect most sensitive to $N_f$ and that closely resembles those seen in equal mass scenarios (Figure 6). With decreasing stellar mass (and decreasing total baryonic mass), the total number of dark quarks decreases at first as ${\mu }_B$ in the core drops. But, as the halo grows in size, dark quarks account for a larger share of the total baryon number as they are transferred to the halo. Eventually, the whole star becomes unbounded, leaving the minimum attainable mass above $\sim M_{\odot }$ where dark quarks make up ∼20–30% of the baryonic contents. For instance, if the mass of the lightest state $m_1\equiv min\left\{m_{\chi }\right\}=260$ MeV, equilibrium configurations below $\sim 1.5\text{–}1.8M_{\odot }$ are forbidden.

It is helpful to estimate the size of the dark quark halo to better understand its impact on minimum masses. Since the halo only accounts for a small fraction of the matter contents, the gravitational potential is dominated by SM baryons. To a good approximation, the spacetime outside the iron surface $R_{★}$ is the Schwarzschild geometry defined by the enclosed mass $M_{★}\approx M$. In this gravitational potential, the constancy of the chemical potential in Equation (29) suggests that ${\mu }_B$ at $R_{★}$ and at the edge of the halo $R_{\chi }\equiv R$ are related as $m_n\sqrt{g_{tt}\left(R_{★}\right)}=3m_{\chi }\sqrt{g_{tt}\left(R_{\chi }\right)}$ where $g_{tt}\left(r\right)=1-2M/r$. It follows that

$$R_{\chi }=\frac{2M}{1-\left(1-2M/R_{★}\right){\left(3m_{\chi }/m_n\right)}^{-2}}.$$

(31)

For $m_{\chi }=m_n/3$, we recover $R_{\chi }=R_{★}$ as expected. With decreasing $m_{\chi }$, the radii of the halo $R_{\chi }$ increase rapidly until they blow up when the denominator vanishes. Requiring that $R_{\chi }$ remains finite leads to

$$m_{\chi }\ge \frac{m_n}{3}\sqrt{1-\frac{2M}{R_{★}}}.$$

(32)

Note that the right-hand side is the asymptotic dark quark chemical potential ${\mu }_{\chi }(r\to \infty )={\mu }_B(r\to \infty )/3$. The origin of the lower limit on hybrid star masses is now clear. To have bound states, the SM core needs to be sufficiently massive to provide the adequate redshift. When it fails to do so, it becomes energetically favorable to get rid of all the SM baryons via $n\to {\chi }_i{\chi }_j{\chi }_k$ as the dense core expands and dissolves into the halo. The lowest energy configuration would then consist solely of $\chi$’s and would have much lower densities at the center where the baryon chemical potential is below the threshold ∼931 MeV. Indeed, such dark quark stars are stable solutions of the TOV equation. For $100\mathrm{MeV}\lesssim {\mathrm{m}}_{\chi }\lesssim {\mathrm{m}}_{\mathrm{n}}/3$, they predict large radii R ∼ 100–5000 km in the mass range of interest.

The bound Equation (32) is depicted in Figure 9 as the solid black line, where we took $R_{★}$∼14 km. The minimum masses of hybrid stars predicted by EOSs discussed previously are also presented. It is quite remarkable that the minimum masses are raised above ∼$1.2M_{\odot }$ for $min\left\{m_{\chi }\right\}\lesssim 270\mathrm{MeV}$, regardless of the underlying SM baryonic EOSs and the number or the mass of heavy dark quarks. This is in strong tension with most pulsars observed to date as they are near the Chandrashekhar limit ∼1.3–1.4$M_{\odot }$ [66,67]. Among the handful of reliably measured pulsars, the secondary component of the neutron star–neutron star binary system J1453+1559 is the lightest with $M=1.174\pm 0.004M_{\odot }$ [68]. The theoretical lower limits on neutron star masses stem from considerations regarding their formation via core-collapse supernovae. While most cold EOSs predict very low masses around ∼0.1–0.3$M_{\odot }$, finite temperature lepton-rich models generally do not support bound states below ∼$M_{\odot }$. Numerical simulations suggest that the minimum is around 1.15–1.20$M_{\odot }$ [69,70], consistent with J1453+1559. Everything considered, we conclude that dark quarks need to be heavier than about $m_{\chi }\gtrsim 270\mathrm{MeV}$ to accommodate low-mass neutron stars.

Additional impacts of dark quark halos may arise in binary neutron star mergers [65]. During the late inspiral phase, component stars develop quadruple deformations due to tidal interactions within the binary. This deviation from point-particle dynamics accelerates the coalescence and leaves unique imprints on the gravitational wave emissions [71]. Within the Post-Newtonian framework, the tidal effect is captured by a dimensionless parameter commonly referred to as the tidal deformability $\Lambda$ [72]. The detection of gravitational waves from the binary neutron merger event GW170817 [73] already provided valuable insights into the $\Lambda$’s of the stars involved. Ref. [65] demonstrated that halos of mass $M_{\chi }$ $\sim {10}^{-3}$–${10}^{-4}{M}_{\odot }$ and radius $R_{\chi }$ ∼ 30–200 km formed by generic sub-GeV dark matter could significantly enhance tidal interactions. Here, the dark quark halos are more massive $M_{\chi }$ ∼ $\mathcal{O}(0.01$–$0.1M_{\odot })$ and can achieve similar enhancement at smaller radii $R_{1.4M_{\odot }}\lesssim 50$ km.

The tidal deformability of $1.4M_{\odot }$ dark quark admixed neutron stars is shown Figure 10. The calculation of $\Lambda$ is described in detail in reference [65]. ${\Lambda }_{1.4M_{\odot }}$ grows rapidly with decreasing dark quark masses. It is interesting to note that, in the limit $m_{\chi }\to m_n/3$, the tidal deformability of hybrid stars does not converge to and is lower than predictions in the absence of $\chi$’s (marked by colored arrows). This is expected since hybrid stars are slightly denser and smaller compared to their SM-only counterparts. Analyses of GW170817 generally favor low values of ${\Lambda }_{1.4M_{\odot }}$ and have led to the putative bound ${\Lambda }_{1.4M_{\odot }}\lesssim 600\text{–}800$ (e.g., references [73,74,75,76,77]). This translates to a lower bound on dark quark masses $m_1\gtrsim 270\text{–}290\mathrm{MeV}$. However, we note that, during the last few cycles of the inspiral phase, surfaces of dark halos would begin to touch and that the description of tidal effects in terms of the single tidal deformability parameter is expected to break down. The halos could be ejected as they merge, leaving the bare SM components behind. Since this occurs at most ∼ seconds before SM baryons merge over a short dynamical timescale, $n-\chi$ conversions would be effectively frozen and would not have the time (Equation (16)) to reestablish new halos. The bottom line is that higher values of ${\Lambda }_{1.4M_{\odot }}$ due to dark quark halos may still be compatible with gravitational waves from GW170817. Detailed studies of the binary evolution are only possible with the help of self-consistent hydrodynamic simulations accounting for dark halos, which will be reported in future work.

A final possibility is the existence of a dark analog of the SM QCD (e.g., references [78,79]). In this scenario, the sum of the masses of $\chi$’s need not be greater than about $m_p-m_e$. Protons and nuclei can remain stable provided that the $\chi$’s are confined to dark baryons in vacuum. Thus, they behave similarly to mirror neutrons in terrestrial experiments. In the dense neutron medium in the interior of neutron stars, however, the dark deconfinement may trigger one (additional) phase transition. Since the zero-temperature phase diagram of SM QCD remains largely a mystery, little can be said about this potential dark QCD phase transition, other than that, if the nature of the transition is first-order, the jump in energy density cannot be too large ($\Delta \mathcal{E}\le$ GeV/fm${}^3$). Transport properties might be useful in searching signatures of a dark QCD sector in neutron stars. We acknowledge that considerable fine-tuning may be required for such scenarios, and a careful investigation of the possibilities outlined above will be reported in future work.

6. Conclusions

Neutron stars are ideal laboratories for testing exotic neutron decays. Their long lifetime and large neutron contents make them sensitive to partial widths several orders of magnitude less than required to explain the current neutron lifetime puzzle [23,80,81]. Sourced by the large baryon chemical potentials, which can reach ∼$2.3$ GeV inside neutron stars, dark decay products carrying baryon number would populate dense neutron-rich matter with sizable abundance and affect the hydrostatic equilibrium configurations. Measurements of neutron star masses are fruitful in constraining the properties of potential dark decay products. These are summarized in

Table 1. Summary of neutron star constraints on neutron dark decays. The bound based on the tidal deformability (last row) is speculative as numerical simulations are required to clarify the evolution of halos during binary inspirals.

Exotic Decay	Sub-Model	Constraint	Observable	Section
Dark Neutron	non-interacting	excluded	$M_{\max }\gtrsim 2.0M_{\odot }$	Section 4.1
$n\to \chi \dots$	mirror neutron	$C_{s,\max }\gtrsim 0.6$		Section 4.2
Fermion $B=1$	self-interacting	$g/m_V\gtrsim 0.01 {\mathrm{MeV}}^{-1}$		Section 4.3
Dark Quark	equal masses	$N_f\le 6$	$M_{\max }\gtrsim 2.0M_{\odot }$	Section 5.1
$n\to {\chi }_i{\chi }_j{\chi }_k\dots$	unequal masses	$min\left\{m_{\chi }\right\}\gtrsim 270$ MeV	$M_{\min }\lesssim 1.2M_{\odot }$	Section 5.2
Fermions $B=1/3$		$min\left\{m_{\chi }\right\}\gtrsim 270$ MeV	${\Lambda }_{1.4M_{\odot }}\lesssim 800$

.

New fermions carrying unit baryon number with a mass close to that of neutrons was first proposed in reference [14] as an explanation for the neutron lifetime puzzle. The fermionic states drastically soften the equation of states and significantly reduce the maximum mass of neutron stars. The existence of two-solar-mass neutron stars thus places strong and robust requirements on the interaction strengths that dark neutrons must possess. If the decay product contains exotic bosonic degrees of freedom, even stronger interactions are required to ensure the bosons do not collapse and form black holes.

The neutron lifetime puzzle may also be explained by novel decays into three dark fermions each carrying baryon number $B=1/3$. We refer to these as dark quarks. Contrary to the dark neutron scenario where $m_{\chi }\simeq m_n$, massive neutron stars can remain stable even if the dark quarks do not experience sizable repulsion, as their underlying Fermi gas EOS supports a higher speed of sound at lower densities, akin to the weakly interacting quark EOSs postulated for dense matter within the Standard Model. Assuming dark quarks have equal masses $m_{{\chi }_i}\approx m_n/3$, up to six flavors of them can be compatible with $2M_{\odot }$ neutron stars. The equal mass case appears to be first discussed in ref. [31], and later in ref. [82]. In scenarios where the dark quark masses differ, light states with $m_{\chi }\lesssim m_n/3$ would form dark halos surrounding SM baryons, extending to large radii in the range ∼20–100 km. The stability of low-mass hybrid stars in the presence of dark quark halos strongly depends on $min\left\{m_{\chi }\right\}$, and the existence of pulsars around 1.1–1.2$M_{\odot }$ rules out dark quarks lighter than about 270 MeV. A dark quark halo could also affect the evolution of binary neutron star mergers, and we will report in future work detailed numerical simulations including halos.