Asymmetric Dark Stars and Neutron Star Stability

We consider gravitationally bound states of asymmetric dark matter (ADM stars), and the impact of ADM capture on the stability of neutron stars. We derive and interpret the equation of state for ADM with both attractive and repulsive interactions, and solve the Tolman-Oppenheimer-Volkoff equations to find equilibrium sequences and maximum masses of ADM stars. Gravitational wave searches can utilize our solutions to model exotic compact objects (ECOs). Our results for attractive interactions differ substantially from those in the literature, where fermionic ADM with attractive self-interactions was employed to destabilize neutron stars more effectively than non-interacting fermionic ADM. By contrast, we argue that fermionic ADM with an attractive force is no more effective in destabilizing neutron stars than fermionic ADM with no self-interactions.


INTRODUCTION
Work on hidden sector dark matter has exploded over the last decade [1]. One conclusion of this work is that even modest extensions of the standard paradigm of dark matter-as a single, stable, weakly interacting particle coupling only via Standard Model forces-to include the dynamics of dark forces can easily change the cosmology, astrophysics and terrestrial signatures of the dark sector [2]. A leading example of this is theories of asymmetric dark matter (ADM), where coupling to dark forces arises naturally as a means to annihilate the symmetric abundance of dark matter [3][4][5][6], similar to the annihilation of electron-positron pairs to photons in the early Universe. When dark forces are present, the cosmology of DM is generically modified due to self-interactions [4,[7][8][9][10][11][12].
When dark ADM forces are sufficiently strong and attractive, bound states can form, similar to the formation of nuclei in the Standard Model [13][14][15]. If the dark sector simultaneously lacks a repulsive long range force (the analogue of the photon), very large states-nuggets-are generically synthesized [13,[16][17][18][19].
The same dynamics that leads to nugget formation in the early Universe can also lead to the formation of ADM stars in the late Universe, via condensation arising from radiation of dark force mediators or small nugget fragments [19]. Self-interacting or not, ADM may also collect in neutron stars. If the ADM is a scalar particle, a black hole can form, destroying the parent neutron star [20][21][22][23][24][25][26][27]. If the ADM is a fermion, Fermi degeneracy pressure tends to stabilize the ADM, though in principle attractive selfinteractions can help to overcome degeneracy pressure.
The primary results in this paper are a self-consistent set of mass-radius relationships and stability bounds of exotic compact objects (ECOs) comprised of fermionic ADM with attractive and/or repulsive self-interactions. We also consider the impact of such fermionic ADM on neutron star stability. We focus on a model with a single stable Dirac spin-1/2 fermion, X, as the dark matter candidate, with attractive self-interactions mediated by a real scalar, φ, and repulsive self-interactions by a vector, We solve the Tolman-Oppenheimer-Volkoff (TOV) equations for this theory to determine the stability against gravitational collapse. We will argue that for ECOs composed of fermionic constituents with arbitrary selfinteractions, the smallest maximum stable ECO size for a given mass scale per constituent, m X , is approximately realized in this model when only a scalar mediator with negligible potential is present. Then we will see that this minimum is of the same order of magnitude as Landau's estimate for Fermi degeneracy supported matter: N max ∼ M 3 Pl /m 3 X , M max ∼ M 3 Pl /m 2 X [28]. We find in particular that spin-1/2 ADM with an attractive force never collapses neutron stars (NSs) over their lifetime in our Universe unless (perhaps) the fermionic constituents are heavier than order 10 6 GeV. This means that fermionic ADM with an attractive force does not in general solve the missing pulsar problem [29][30][31], have limits from imploding NSs [32], or lead to non-primordial solar mass black holes [33]. We also find a different equilibrium sequence for stars made of self-attractive ADM than derived elsewhere [34]. The most important difference between our work and previous treatments is use of a fully consistent equation of state (EoS), instead of utilizing a Yukawa potential valid only in the non-relativistic and low-density limit in the case of scalar-mediated interactions. 2 This crucially changes both the impact on NS stability and the ADM star equilibrium sequence. Our fully relativistic treatment extends to a fully general relativistic treatment of ECOs composed of two possibly interacting but separately conserved constituents-in our case baryonic matter and fermionic ADM.
The potential structure and stability of fermionic dark matter ECOs [34,35] and of ADM-admixed NSs [36][37][38][39][40][41][42] has been examined before, and there has been renewed interest in such objects in the context of gravitational wave observations [43][44][45][46][47]. This work gives a comprehensive account of the effect of interactions, including the first correct treatment of attractive self-interactions that cause binding, and baryon-ADM interactions in the case of admixed stars.
The outline of this paper is as follows. In Sec. II we specify and interpret the EoS for spin-1/2 dark matter with attractive and/or repulsive scalar-and/or vectormediated self-interactions. In Sec. III we find and interpret the sequence of gravitationally stable stars composed of such matter. Then in Sec. IV, employing results from the previous section, we argue that ADM with spin-1/2 constituents smaller than about 10 6 GeV cannot collapse NSs, regardless of whether the constituents are selfinteracting. Finally, in Sec. V we conclude. Appendix A explains non-generic features of ADM at the cusp of being self-bound. Appendix B lays out the general relativistic equations appropriate for determining structure and gravitational stability of static stars composed of multiple separately conserved, possibly interacting, components. It details methods we used to obtain numerical solutions that confirm the less technical arguments presented in Sec. IV.

II. EQUATION OF STATE FOR SELF-INTERACTING SPIN-1/2 ADM
The EoS for fermionic matter described by Eq. 1 is given by [48][49][50][51][52] where The effective Dirac mass (c.f. Eq. 1) is m * = m X (1 − ϕ) and the number density is given by n = 2 The equations above assume zero temperature, though generalization to nonzero temperature is straightforward and has been worked out in the context of the σ-ω model (see e.g. [51,52]). The equations are derived in the mean field limit, where scalar and vector fields are approximated by their mean values. Additionally the mean fields are assumed to be static and spatially uniform. This last assumption is inconsistent with solutions to the general relativistic equilibrium equations when order one variations in star density occur over length scales comparable to or smaller than the force range, 1/m φ , 1/m V . For example with C 2 φ fixed, the approximation breaks down in the decoupling limit, α φ → 0. We explored the transition where spatial uniformity breaks down for self-bound matter in [15].
Eqs. 2 and 3 are related through the thermodynamic relation where µ = ∂ ∂n is chemical potential. Note that rest energy per constituent, /n, is necessarily minimized when p = 0. For large enough attractive interactions, there are solutions where p = 0, ∂p ∂n > 0, and the binding energy per particle is positive (m X − /n > 0), meaning that large stable self-bound states exist, elsewhere called nuggets [14,15,18,19]. Fig. 1 shows the rest energy per particle as a function of number density for C 2 φ = 0.5, 10, C 2 V = 0, and V (φ) = 0. The solid/dotted lines show Eq. 2 while the "semi-relativistic" dashed lines show the energy computed assuming = kin + Y with potential energy given by the nonrelativistic expression , with indices i and j running over particles, and kin = 2 3 , as was done explicitly or implicitly in Refs. [29][30][31][32][33][34], for example. As number density grows and the constituents grow more relativistic, the fully relativistic and semi-relativistic expressions diverge. Were it correct, the semi-relativistic expression would imply that self-attractive ADM is microscopically unstable such that the energy per constituent becomes negative at high density and is unbounded below. By contrast, the fully relativistic expression for energy per constituent remains positive but can develop a local or global minimum at nonzero density. This happens because Fermi pressure overcomes the attractive force at high density and the pressure grows again. More specifically, as density grows, φ grows, initially decreasing pressure, but = 0.5 (thin) and 10 (thick) using a semi-relativistic treatment (dashed) alongside the expressions obtained using relativistic mean field theory (solid). When C 2 φ = 10, ADM is self-bound; the dotted line shows the analytic EoS in the density range below the density of self-bound nuggets, where the matter is unstable to coalescence.
simultaneously decreasing the effective Dirac mass, accelerating the growth of Fermi degeneracy pressure.
For C 2 φ > 1.09 there is a global minimum in /n, with energy per constituent at this minimum less than m X (the value at n = 0), indicating the existence of large stable bound states, or nuggets, 3 that form in principle without the aid of gravity. The analog of surface tension causes large states to form spheres [15]. The large point on the C 2 φ = 10 curve in Fig. 1 represents the saturation density-the density of large self-bound nuggets. The dotted curve lies in a density region with mostly negative pressure, representing the instability of ADM to condensation into large, dense nuggets.
Numerical solutions to Eqs. 2 and 5 for 0.840 < C 2 φ < 1.09 reveal a local minimum with /n > m X , indicating a phase change at positive pressure. At this pressure, zero-temperature matter jumps from a gas-like state at low density to a liquid-like state at higher density. The liquid state can be realized only with the help of another force-for example with the aid of gravity in the core of a star. See Appendix A 1 for further discussion. With the EoS for ADM self-interacting through scalar and vector exchange in hand, we now explore the structure of self-gravitating objects composed of such matter.

III. ADM STAR STABILITY AND EQUILIBRIUM SEQUENCE
In the previous section, we showed that self-attractive fermionic ADM is microscopically stable. Our main ob- 3 For large enough C 2 φ , local and global minima also exist when C 2 V = 0, V (φ) = 0. See e.g. [15], [19].
jective here is to pinpoint when cold ADM stars become gravitationally unstable. In particular we will show that the maximum stable mass cannot deviate much below Landau's estimate, M fermion max ∼ M 3 Pl /m 2 X [28] (see also [53]). An implication is that self-attractive fermionic ADM cannot seed collapse of neutron stars more efficiently than non-self-interacting fermionic ADM.
We restrict our attention to spherically symmetric compact objects such that the Tolman-Oppenheimer-Volkoff (TOV) equation governing such objects reads, dp dr = − (p + )(GM (r)/r + 4πGr 2 p) r(1 − 2GM (r)/r) (7) or equivalently where r is the radial coordinate, M (r) = 4π r 0 r 2 dr , c 2 s = dp d is squared sound speed that characterizes the stiffness of matter, and we have used Eq. 6. Given an EoS relating and p, this single integro-differential equation can be solved for any given choice of central energy density by integrating out from r = 0 to the edge of the star r = R where p(R) = 0. Eqs. 2 and 3 represent a parametric EoS, {p(n), (n)}, which leads to an integrodifferential equation for number density, n, as a function of r. The number of constituent particles in the star is given by N X = n(r)/ 1 − 2GM (r)/r d 3 r.
The gravitational stability of ADM stars is calculated from the equilibrium configurations, which are solutions to the TOV equations for a given central density. Maxima in mass as a function of central density indicate a transition from stable to unstable [53]. Figs. 2, 3, and 4 show the mass and radius of solutions to the TOV equations for spin-1/2 ADM with various self-interaction strengths. Figs. 2 and 3 show similar results to those in [34,35] for repulsive interactions. However, for m X ∼ {10 GeV, 100 GeV, 1 TeV} the attractive interactions corresponding to m φ = 10 MeV, α φ ∈ {10 −2 , 10 −3 , 10 −4 } increase the maximum mass (and maximum number of constituents, N X ) of a stable gravitationally bound ECO relative to the case of non-interacting fermions of the same mass, in contrast to the results in [34]. For comparison, the magenta region in Fig. 2 includes equilibrium configurations for NS matter that can support a 2M star, satisfy the 90% confidence level constraint on tidal deformability from the NS binary inspiral gravitational wave observation, GW170817 [56,57], and are consistent with known limits on the baryonic matter EoS at (low) nuclear densities and at very high densities [55].
The attractive self-interaction parameters shown in Figs. 2 and 3 all correspond to and thus to the case of strong self-binding; the growth of radius with mass until near the gravitational stability endpoint (here marked with asterisks) is characteristic of compact objects made of self-bound matter such    Table I. Asterisks mark the stable equilibrium sequence endpoints corresponding to the global maximum in M as a function of central density. The gray region corresponds to R < Rs = 2GM and the cyan contour represents maximum compactness, GM/R = 0.354. Compare Fig. 3 in [34]; the repulsive case agrees but the attractive case dramatically differs due to differences in the EoS. Right-hand plots show energy density for a benchmark star (marked with in the left-hand plots) as a function of distance from its center. The cutoffs at finite density in the bottom right-hand plot indicate the discontinuity in energy density at r = R due to self-boundedness. For comparison, the middle plot also shows the equilibrium sequence for a sample NS matter EoS (magenta) consistent with NS observations to date-the HB EoS as defined in [54]. The NS benchmark is a 1.5M star, and the shaded magenta region is digitized from [55].  as those of hypothetical self-bound strange quark matter stars [58]. The increased gravitational stability of these objects relative to their non-interacting counterparts stems from their effectively stiffer EoS due to enhanced Fermi degeneracy pressure all the way out to the edge of the star. For self-interactions satisfying C 2 φ 1 and negligible scalar potential, we have identified universal formulas for N max , M max and GM R max . We report these in Table I alongside analogous relations for the cases of no interactions, purely repulsive interactions, and attractive-repulsive interactions with The formulae work well when C 2 i 10. In the attractive case, the asymptotic values may alternatively be written N attractive wherem X is the zero-pressure chemical potential, equivalent to the average mass per constituent of large nuggets. [15,19]. An alternative explanation of the enhanced gravitational stability 4 Our formula for Mmax in the purely repulsive case matches that reported in [35]. of matter with large attractive self-interactions is that these interactions decrease the effective constituent mass scale, so the Landau limit M max ∼ M 3 Pl /m 2 still applies but with m to be interpreted asm X wherem X < m X .
A scalar potential term of the form λφ 4 with λ > 0 tends to limitm X from below [15,19], reducing the possibility of substantially raising M max for fixed m X . In general, such a term tends to push c 2 s closer to its form with no interactions, indicating that scalar potentials tend to push the maximum mass and other equilibrium sequence characteristics closer to that for no interactions. Fig. 4 shows equilibrium sequences for matter with more moderate (C 2 φ , C 2 V 10) attractive, repulsive, or both attractive and repulsive self-interactions side-byside with squared sound speed, c 2 s = dp d , characterizing the stiffness of the matter. The figure demonstrates that softer equations of state lead to smaller maximum masses and vice versa: the larger the density range where c 2 s lies below the no-interactions curve, the smaller the maximum mass relative to the no-interactions case and vice versa. In the top left plot, we see that any softening of the EoS relative to the no-interactions EoS accelerates the approach to the high density limit, c 2 s ∼ 1/3. no interactions attractive both repulsive 21   TABLE I. Mass, M , number of constituents, N , and compactness, GM R , for the static spherically symmetric maximum-mass stars made of spin-1/2 matter with large attractive, repulsive, both attractive and repulsive, and no interactions. The dimensionless In the attractive case, the matter is strongly self-bound with the chemical potential at zero pressure equal to µ|p=0 =mX = mX (2/C 2 φ ) 1/4 . In all other cases shown the matter is not self-bound so µ|p=0 = mX . For quick reference, note: And in the top right figure, comparing purple to blue, we again see that attractive interactions soften the EoS at low densities but stiffen it a larger densities, accelerating the approach to the high density limit when a vector is present, c 2 s ∼ 1.
In all examples, the more extreme the softening at lower densities, the more extreme the stiffening at higher densities. This can also be seen analytically as follows.
The chemical potential for the model Eq. 1 is given by with m * (k F ) determined through Eq. 5. With a scalar potential guaranteeing positive energy density and therefore microscopic stability, one can show that m * approaches zero in the large-density limit; and larger C 2 φ drives m * to zero faster while a quartic potential term moderates the decrease. For spin-1/2 matter in general, c 2 s = dp d = 1 3 d ln µ d ln k F by Eq. 6 and n = k 3 F /3π 2 . In the large density limit, vector repulsion dominates pressure and µ ∼ k 3 F /m 2 X if a vector is present, or fermi pressure dominates and µ ∼ k F if the vector is absent; thus c 2 s → 1 or 1/3 with or without vector repulsion, respectively. Furthermore since attractive interactions cause m * to decrease, though this initially drives c 2 s below its value absent the attractive interactions, it also accelerates the approach to the asymptotic limit. This explains the tendency of attractive interactions to soften the matter at lower densities and stiffen it at higher densities.
The left plots in Fig. 4 represent matter with attractive self-interactions, including examples of non-self-bound matter (C 2 φ = 0.5) and self-bound matter (C 2 φ = 10). The equilibrium sequence for C 2 φ = 0.5 begins at low central density and large radius following the no-interactions sequence, but ends at lower maximum mass and smaller radius. By contrast, since self-bound matter fuses before reaching a cool state, the C 2 φ = 10 equilibrium sequence begins at relatively large densities and small radii, and never meets up with the no-interactions sequence. For C 2 φ < 0.516 and C 2 φ > 1.09, the first maximum in M as a function of central density is the global maximum.
In the purely attractive, V (φ) = 0 case, the equi-librium solutions for 0.516 < C 2 φ < 1.09 hug the nointeractions sequence at low central density and large radius, and then develop a local maximum at low compactness before reaching the global maximum at larger compactness, indicating the existence of two separate sequences analogous to the white dwarf and NS sequences (see e.g. [52,59]). We further discuss this range in Appendix A 2, but note that the first local maximum occurs in the range 0.
Pl /m 2 Xless than a factor of ten lower than the no-interactions global maximum of 0.38M 3 Pl /m 2 X -except in the narrow range 1.05 < C 2 φ < 1.09.
For purely attractive interactions and fixed m X , the global maxima in the entire C 2 φ range satisfy M max > 0.23M 3 Pl /m 2 X , N max > 0.24M 3 Pl /m 3 X , and (GM/R) max > 0.092 with the bounds saturated when C 2 φ = 0.45, C 2 φ = 0.42, and C 2 φ = 0.26, respectively. In each case, the parameter decreases from the no interactions value to the value at the minimum, and then increases monotonically toward the asymptotic value in Table I. In the case of equal strength attractive and repulsive interactions, as the interactions become more extreme with C 2 φ = C 2 V 1, the speed of sound curve approaches a step function, c 2 s ∼ θ(n − n crit ). Matter with such behavior is thought to produce the theoretically most compact stars, and indeed we find that GM R max ∼ 0.354, the posited maximum in the literature assuming subluminal sound speed [60,61], for also  Table I.) We now consider self-interacting fermionic ADM more generally. If the cost of softening matter at a given density through attractive interactions is accelerating the approach to the asymptotic limit, generally, then the softness of microscopically stable fermionic matter is limited, and therefore the amount that self-interactions can reduce the maximum stable mass below that for free fermionic matter is limited. On this basis, we conjecture that M max 0.1M 3 Pl /m 2 X , N max 0.1M 3 Pl /m 3 X holds true for spin-1/2 ADM, generally. s = dp d , as a function of (kF /mX ) 3 for spin-1/2 dark matter with attractive (red), repulsive (blue), both attractive and repulsive (purple), and no (thick black) self-interactions. The strength of attractive and repulsive self-interactions is characterized by C 2 φ and C 2 V , respectively, as defined in Eq. 4. Here kF is Fermi momentum and X number density is n = k 3 F /3π 2 . The dotted section of the red c 2 s curve corresponds to the dotted region in Fig. 1, where matter is unstable to condensation into large self-bound states with density marked by the dot. Bottom: Mass and radius of static, spherically symmetric stars composed of such matter, representing the equilibrium sequence. The curves cut off at the maximum-mass gravitationally stable stars, denoted with asterisks. Gray regions correspond to (R < Rs = 2GM ). The cyan boundary marks GM R = 0.354, corresponding to the theoretically most compact non-black-hole objects [60,61].

IV. IMPLICATIONS FOR NEUTRON STAR COLLAPSE
So far we have focused on ADM-only stars, including pinpointing the maximum mass of gravitationally stable self-attractive ADM stars. Based on a semi-relativistic calculation, Refs. [29][30][31][32][33] have claimed that ADM with an attractive force, captured inside of NSs, can destabilize and destroy them. Here we argue that relativistic effects stabilize the NS over most of the parameter space, and destabilization can occur only for very large ADM mass of order PeV.

A. ADM Capture
The amount of ADM captured in a NS in a time t is at most the amount that impacts the NS [20], where b max = 2GM/R 1−2GM/R R vDM is the impact parameter corresponding to DM that just scrapes the surface of the NS at closest approach. 5 Here the energy density ρ DM and velocity scale v DM are to be taken asymptotically far away from the NS. For typical NSs, GM/R ∼ 0.2 and R ∼ 10km leading to The upper bound is realized only when on average 100% of the DM passing through the NS deposits enough energy to be captured. For DM with mass of order GeV to PeV, the minimum required baryon-DM cross section for this to occur is order 2×10 −45 cm 2 [20]. For m X GeV, Fermi blocking suppresses the scattering [27], and for m X PeV, multiple scatters are required for gravitational capture [62] and therefore greater cross sections are required to realize the upper bound in Eq. 11. If the ADM is bound in large composite states at the time of capture, additional considerations apply [63]. In general since the density and velocity scales entering Eq. 11 are not vastly different from the fiducial values, the amount of ADM captured is a tiny fraction of the total mass of a NS (order 1.5M ). One can speculate about other ways to realize ADM-admixed NSs stars with a much larger fraction of ADM than can be collected gradually through capture over the NS lifetime. One interesting possibility is copious production and capture of dark matter in the core-collapse supernova of the NS's progenitor [44].

B. ADM-Admixed Neutron Stars
For ADM captured by a NS over its lifetime to induce its collapse, a self-gravitating ADM star of the same mass as the captured ADM must itself be unstable to collapse. Based on this observation, we argue that capture of spin-1/2 ADM-self-interacting or not-with constituent masses smaller than order 10 6 GeV by NSs cannot in any circumstances induce collapse.
As detailed in Appendix B, we solved the general relativistic equilibrium equations for NS matter admixed with cold ADM as modeled in Sec. II. In general, the maximum number of ADM constituents, (N X | N b =0 ) max , possible in a stable ADM-admixed NS with fixed baryon number can be smaller than the maximum number of ADM constituents in an ADM-only star (N X | N b =0 ) max made of the same kind of ADM. However, we find that any appreciable differences between ( ) max is comparable to a solar mass. Including baryon-ADM interactions does not affect this conclusion. And including thermal effects only increases stability. Now assume the amount of captured ADM is a small fraction, of the NS mass, as predicted by Eq. 11, and that the NS is not already teetering at its stability bound with mass greater than 2M . Then as long as an ADM-only star of size N X = N X cap does not collapse, neither will an ADM-admixed NS with the same amount and type of ADM. From our treatment of fermionic ADM-only stars, since taking M NS ∼ 1.5M we find Assuming maximally efficient ADM capture over 10 billion years and galactic ADM densities ρ DM ∼ 1-100 GeV/cm 3 and the velocity scale v DM ∼ 200 km/s, we find m X, collapse 2 × 10 5 -2 × 10 6 GeV.
We have not yet argued that ADM can collapse NSs, but rather only that ADM with spin-1/2 constituent mass m X < (10 f ) −1/2 GeV cannot collapse NSs if f 1. Now we consider whether ADM can collapse NSs in certain cases when m X > (10 f ) −1/2 GeV. More specifically, if an ADM-only star with constituents N X = N X cap is unstable to gravitational collapse, then is an ADM-admixed NS with N X = N X cap also unstable? If there is a process to cool the ADM sufficiently that it self-gravitates, the answer is yes. Let us briefly consider the ADM cooling process after capture. Ref. [64] estimates the cross section needed to maximize ADM capture, saturating Eq. 11, and the time required for the ADM to deposit most of its kinetic energy in the NS. When m X 10 6 GeV and Eq. 11 is saturated, for example, a captured ADM particle deposits most of its kinetic energy in less than a day. The ADM heats the NS, which will be detectable by the next generation of infrared telescopes [64]. Again if the ADM-baryon cross section is anywhere near the level required to maximize ADM capture, then according to the estimates in [23], thermalization of the ADM with baryonic matter happens quickly. Accounting for NS heating through ADM capture, the maximum blackbody temperature of NSs near Earth is expected to be 1750 K (0.15 eV) [64], which is a minuscule fraction of the Fermi momentum of nearcollapse ADM with m X 10 6 GeV. The ADM indeed cools to effectively zero temperature and for ADM with M ADM-only max M , once the amount of captured ADM approaches this maximum, the ADM core density within the NS far exceeds the baryon density. The ADM core self-gravitates and its structure is unaffected by the baryonic matter. When N X cap > N ADM-only max , the ADM core is unstable to collapse.

V. CONCLUSIONS AND OUTLOOK
We have argued that the capture of spin-1/2 ADM by neutron stars cannot lead to their implosion unless the mass of the spin-1/2 ADM constituents exceeds approximately one PeV. This includes ADM with attractive or other varieties of self-interactions. Thus the existence of old NSs can only set limits on ADM-baryon cross sections for spin-1/2 ADM constituent masses larger than about one PeV. Once the next generation of infrared telescopes comes online, limits from dark kinetic heating of NSs in this mass range may compete with any such limits [64,65]. If there is a positive detection of dark kinetic heating of NSs, the existence of old neutron stars will provide complementary information on possible models of ADM with m X PeV.
After deriving and interpreting the equation of state for cold spin-1/2 ADM self-interacting through scalar and/or vector mediators, we found solutions to the general relativistic gravitational equilibrium equations for cold static spherically symmetric ADM stars and identified the maximum size of gravitationally stable ADM stars. We found formulas for this maximum stable size in the case of strong self-interactions (see Table I). We also found that the maximum size at fixed m X generically does not drop below M star max = 0.1M 2 Pl /m 2 X , and we conjectured that a similar limit holds for spin-1/2 ADM with arbitrary self-interactions.
If fermionic ADM stars are realized in our Universe, they might be detectable by gravitational wave observatories [66,67] in the event of mergers with other compact objects [43]. The masses and radii of these objects can be drastically different from those of NSs, though not necessarily so. As compared to a NS binary merger, we expect the electromagnetic signature of a merger involving at least one ADM star to be a smoking gun signal of the difference even if the gravitational wave form does not reveal the presence of the ECO. As compared to a black hole binary merger, the waveform will be modified due to the ADM star's spatial structure and tidal deformability. We leave the prospects of gravitational wave detection of ADM stars for future work. Purely self-attractive ADM with C 2 φ > 1.09 is selfbound, not relying on gravity for its boundedness. Here we focus on ADM that is not quite self-bound.
Zero temperature matter with purely attractive interactions of strength 0.840 < C 2 φ < 1.09 and V (φ) = 0, as modeled in Sec. II, undergoes a phase change at positive pressure and number density. Further, when 0.516 < C 2 φ < 1.09, the equilibrium sequence obtained by solving the TOV equations has an additional unstable region as compared to the sequences described in the main text. Rest energy per constituent, sound speed, and gravitational equilibrium mass and density configurations are shown for this range of couplings in Fig. 5. The dotted regions show unphysical solutions to Eqs. 2-5, signaled by a local minimum in energy per constituent as a function of number density and a negative squared sound speed, c 2 s = dp d . The dashed regions in the bottom plots represent unphysical solutions to the TOV equations, separating two separate gravitationally stable equilibrium sequences that emerge due to a temporary stall in the growth of sound speed with density. In the next two subsections we discuss how the physical solutions are constructed, first examining the equation of state, and then considering the equilibrium sequence obtained from the TOV equations.

Phase Change and Maxwell's Construction for the Equation of State
We first consider the rest energy per constituent and speed of sound, shown in the upper two panels of Fig. 5. The local minimum in the analytic /n, C 2 φ = 1 curve at nonzero density indicates a phase change. Matter lying near the local maximum in /n can lower its energy per constituent (and pressure) by condensing, and matter between the local maximum and minimum has negative pressure. As seen in the upper right panel of Fig. 5, sound speed also becomes imaginary near the local maximum. These features all signal the unphysical nature of the analytic EoS in this density domain. The physical EoS is obtained by choosing endpoints n A and n B that lie at equal pressure on either side of the density region with ∂p/∂n < 0 such that ∂ ∂n A = ∂ ∂n B = B − A n B −n A . This construction is equivalent to Maxwell's construction in standard thermodynamics [53], and is shown by the solid red line labeled "physical" in Fig. 5. Matter at densities less than n A exists in a stable gas-like state and matter at densities greater than n B exists in a stable liquid-like state. When there is a phase change, we use "liquid" and "gas" to refer to the high-density and low-density phases, respectively. Matter in the density region between n A and n B coexists in two different density states (liquid and gas). The phase change occurs at nonzero pressure; to realize this nonzero pressure, some other force must be applied-for example, gravity. The liquid state can exist in the cores of gravitationally bound stars; in static spherically symmetric stars an abrupt transition from liquid core to gas crust occurs at a given radius.

Two Equilibrium Sequences when
0.516 < C 2 φ < 1.09 Next we consider solutions to the TOV equations. For purely attractive interactions and 0.516 < C 2 φ < 1.09, the sequence of solutions to the TOV equations contains a local maximum in M as a function of central density before hitting a global maximum, as shown in the bottom right panel of Fig. 5. The local maximum at lower density occurs at larger radius and lower compactness, as shown in the bottom left panel. This same behavior occurs for cold catalyzed matter due to a relative softening of the EoS near the neutron drip density; near this density the speed of sound temporarily decreases (softening) before increasing again (stiffening) as a function of density (see e.g. [52,53,59]). Correspondingly, a local maximum and minimum in star mass as a function of central density develops near the neutron drip density, matching onto the stability endpoint of the white dwarf sequence and the beginning of the NS sequence, respectively. In Fig. 5, analogs of the white dwarf sequence are represented by solid lines stretching from low central density and larger radius up to the local maxima marked by diamonds. Analogs of the NS sequence are represented by solid lines that stretch from the local minima in mass to the global maxima at higher central densities and smaller radii, marked by asterisks.
Another way of characterizing the feature of cold catalyzed matter that leads to the separate white dwarf and NS sequences is that its sound speed growth (or stiffening) temporarily stalls near the neutron drip density. This is similar to what we see in the top right panel of the figure. For C 2 φ = 0.52, the stall is moderate, with c 2 s continuing to increase but at a lower rate near densities n ∼ 0.2 m 3 X 3π 2 . Correspondingly, as seen in the bottom right panel, the growth of gravitational mass with central density near n(0) ∼ 0.2 m 3 X 3π 2 stalls so much that a local maximum develops, signaling an instability to contraction. The dashed gray region lying between local maximum and minimum represents gravitationally unstable solutions to the TOV equations. In this case, the instability is relatively mild, with stability taking hold again at only slightly higher central densities. For C 2 φ = 0.75, the stall in sound speed growth is more severe. Correspondingly, the the minimum in star mass as a function of central density is deeper, and the density gap between equilibrium sequences is larger. The stall grows more severe with increasing C 2 φ . As demonstrated by the C 2 φ = 1 curve of the bottom left panel, for matter with a phase change, a cusp develops near the local maximum in the TOV solutions for M (R) and is associated with the discontinuity in density because of the phase change. Solutions to the TOV equations with central densities corresponding to the coexistence density range do not exist; there is a gap in the mass versus central density curve where n A < n(0) < n B , marked by the dotted red line in the lower right panel. This entire gap is mapped onto the cusp in the lower left panel. Such behavior also occurs in models of compact stars that include QCD phase changes, see e.g. [68,69].
Understanding the final states of stars or ADM cores within NSs that reach the white dwarf-like stability endpoint-be they NS-like or black holes-is beyond the scope of this work. But we remark that the white dwarf-like mass endpoints lie within a factor of ten below the no-interactions global stability endpoint except in the range very near the transition to self-bound matter, 1.05 < C 2 φ < 1.09.

Appendix B: ADM-Admixed Neutron Stars and Baryon-ADM Interactions
In the bulk of the paper we focused on ADM-only stars. Here we lay out the equations for baryon-ADM admixed stars. Then based on our solutions to these equations, we argue that, unless the mass of the ADM and baryons in the star is of the same order, the amount of ADM that destabilizes a NS is little modified from the amount that destabilizes an ADM-only star.

Gravitational Equilibrium Equations
To investigate ADM-admixed NSs, we need the analog of the TOV equations for two interacting but separately conserved matter species. The TOV equations are equivalent to extremizing mass, with baryon number, held fixed [53]. To generalize to multiple conserved species, we extremize M with each conserved species number separately held fixed through the method of Lagrange multipliers. That is, for ADM-baryon stars, extremize the functional F = M − µ b N b − µ X N X , treating baryon density, n b , and ADM density, n X , as independent functions. The Lagrange multipliers µ i are gravityinclusive chemical potentials. The derivation is worked out for N such species in detail in Ref. [70]. One finds, where µ i is the gravity-inclusive chemical potential of species i, p = i ∂ ∂ni n i − , and ν is a metric function defined by |g tt | = e 2ν . The TOV equation, Eq. 8, is equivalent to i n i d dr µ i = 0, with ν eliminated through Eq. B4.
The total mass M , constituent numbers N X , N b , and radius R, are determined by the equilibrium equations, Blue dots on the C 2 φ = 1 curve mark the matching points for Maxwell's construction; the matter is in its liquid phase at densities n > nB, in its gas phase at densities n < nA, and coexisting in both phases in between. The dotted red lines in the top panels represent the unphysical analytic EoS in the coexistence density range. In a star, the phases do not coexist but rather density abruptly jumps from nA to nB at a given radius, and equilibrium solutions with central density nA < n(0) < nB do not exist; this region is marked with a dotted red line in the bottom right panel, which maps onto the cusp in the bottom left panel. The asterisks mark the maximum mass stable star in the higher-density NS-like sequence while diamonds mark that of the lower-density white dwarf-like sequence. Dashed gray lines indicate gravitationally unstable equilibrium solutions between the two sequences. See the text for further detail.
Eqs. B3-B4, along with the equations defining M and N i , Eqs. B1-B2, for given central baryon and ADM number densities, {n b (0), n X (0)}. 6 We interpret the largest stable star at fixed baryon number, N b0 , to correspond to the first local maximum in M as a function of {n b (0), n X (0)}, subject to the constraint N b = N b0 . Since equilibrium configurations satisfy dM = µ b dN b + µ X dN X = 0, and since µ i are finite for any equilibrium 6 We will discuss a caveat when significant baryon-ADM interactions are present. configuration, with fixed N b = N b0 , M and N X attain their maxima at the same {n b (0), n X (0)}. 7 When the ADM constitutes a tiny fraction of the star by mass, it is numerically easier to identify the maximum in N X .
Absent DM-baryon interactions, the number density of a given species affects the other only through the total gravitational mass and pressure. In this case, Eq. B3 and Eq. B4 are equivalent to where p i = n i ∂ i ∂ni − i . Given a large hierarchy between the densities of the two species, as we will detail below, even a weak DM-baryon interaction can dramatically affect the density profile of the subdominant species where the two species overlap; in this case it is important to use Eqs. B3-B4 rather than Eq. B5. Furthermore, in such cases, we find there can be multiple distinct equilibrium configurations corresponding to zero central density of the component with lighter constituents. In the next section we describe how to include baryon-ADM interactions before discussing our numerical solutions to the gravitational equilibrium equations with and without ADM-baryon interactions in Sec. B 3.

Modeling Baryon-ADM Interactions
Using relativistic mean field theory and the same techniques used in the context of the σ-ω model of nuclear physics (see e.g. [51,52]), we find that a vector-mediated baryon-ADM interaction gives rise to an interaction energy density and pressure where g b and g X are the vector-nucleon and vector-ADM coupling constants, respectively, and m A is the vector mediator mass. The low energy elastic nucleon-X scattering section is given by σ bX = with µ bX the reduced mass, so Eq. B6 is alternatively written, where m b ≈ GeV is the nucleon mass scale. The total energy density in an admixed star is = b + X + I with b independent of n X , and X independent of n b , and similarly for pressure. 8 Consider m X GeV. Current direct detection constraints on σ bX in this range are σ bX 10 −46 m X 100 GeV cm 2 ≈ 10 −19 m X 100 GeV GeV −2 [71][72][73]. Baryon number densities (energy densities) toward the centers of NSs are order 10 −2 GeV 3 ( GeV 4 ). Thus given interaction strengths near current direct detection limits, interaction energy density is comparable to baryon energy density when m X GeV n X 10 10 GeV 3 . The number density of free fermionic ADM reaches order m 3 X in the cores of near-collapse stars, implying the bound can be satisfied for m X TeV dark matter. The hierarchy I , b X with I b naturally occurs for m X TeV when ADM-baryon interactions are near current direct detection limits. In this case, both the gravitational and non-gravitational ADM-baryon interactions affect the baryon density profile in the small admixed core, while the dark matter density profile is unaffected by the baryonic matter.
Conversely, direct detection bounds on sub-GeV dark matter become weak, and furthermore for ADM masses much smaller than a GeV, we expect X number (energy) densities to be much smaller than order GeV 3 ( GeV 4 ). In this case ADM-baryon interactions can be important in determining density profiles of sub-GeV dark matter within ADM-admixed NSs.

Results: Numerical Solutions to the Equilibrium Equations
The left-hand plot in Fig. 6 shows contours of constant N X , N b , and M as functions of central (r = 0) ADM and baryon number densities for solutions to the admixed star equations, Eqs. B3 and B4, with noninteracting zero-temperature, m X = 100 m b ADM, where m b ≡ 939.5 MeV is the nucleon mass scale. The baryonic matter was modeled with the HB EoS as described in [74] and shown in Fig. 2, though this detail is unimportant. The important points are: • The constant N X contours are independent of central baryon density, and N X max at fixed baryon number is the same as for an ADM-only star, demonstrating that baryonic matter affects neither the structure of the ADM core nor its stability endpoint.
• The maximum M and N b at any fixed value of N X are the same as a baryons-only star up to negligible fractions, demonstrating that the NS matter stability endpoint is unaffected by the ADM. This is because the ADM's contribution to the total mass of the NS is small: m X N X ≪ M . The curvature of the M , N b contours indicates that the ADM affects the baryon density profile at the center of the star. C.f. Fig. 7.
• The X number density and Fermi momentum for solutions near the stability endpoints, where N X = (N X | N b ) max , are order 10 −2 m 3 X and m X , respectively-much greater than the baryon density, and also relativistic so that even if the NS is relatively warm, the zero-temperature approximation used for the ADM EoS is still valid.
We checked that these features also hold for selfinteracting ADM when M ADM-only star max M . Thus for f 1, Eq. 14 is accurate.
By contrast, the right-hand plot of Fig. 6 shows similar contours when m X = m b . The shape and overlap of M , N b , and N X contours along the diagonal from near the bottom left to top right are highly interdependent because here the ADM and baryonic matter densities are similar and both components contribute similarly to the total star mass. For N b fixed at its value corresponding to a 1.5M non-admixed NS (N b = 0.66N bmax ), the maximum-mass stable admixed star corresponds to M ∼ 1.7M and N X ∼ 0.4N X max . The amount of m X = m b ADM that destabilizes the NS is smaller than the amount that destabilizes an ADM-only star. But the amount by mass is a sizable fraction of the remainder of baryonic matter that would destabilize a baryons-only NS, and by Eq. 11 far exceeds the amount that can be captured.
The main message conveyed by Fig. 6 is: small amounts of fermionic ADM by mass (m X N X M ) cannot induce collapse of a NS unless the ADM is concentrated in a very dense ( X GeV 4 ) core. Due to Fermi degeneracy pressure, this can occur only when m X GeV. And in this case, the amount of ADM that destabilizes the NS is negligibly modified from the amount that destabilizes an ADM-only star.
We find that dense ADM cores can dramatically affect the density profile of baryonic matter overlapping the core, even though as noted above the maximum mass, radius, and baryon number at fixed ADM number is affected by negligible fractions. This is shown in Fig. 7, with baryon-ADM interactions as modeled in Sec. B 2 present or not. With even moderate σ bX well below direct detection constraints on large m X ADM, there are classes of solutions with vanishing n b at r = 0-because of the high ADM densities, the repulsive interaction between ADM and baryonic matter wins over gravitational attraction and expels the baryonic matter from most of the core. In these cases we scan solutions by setting the ADM density at r = 0 and the baryon density at the edge of the ADM core. We checked that including the baryon-ADM interaction at levels allowed by direct detection does not affect the conclusions described in the previous paragraphs. Finally we note that it could be interesting to examine the structure of NSs admixed with sub-GeV constituent mass ADM given repulsive ADM-baryon interactionsthe ADM could be concentrated at the outer edge of the NS and beyond, which could lead to larger-than-naivelyexpected effects on the tidal deformability and/or moment of inertia of the NS for a given total amount of collected ADM. . These plots demonstrate that, unless the amount of ADM is of the same order or larger than the baryonic matter by mass (mX NX m b N b ), the amount of ADM that destabilizes a NS is little modified from the amount that destabilizes an ADM-only star. Configurations lying below the maximum in N b (at fixed NX ) and to the left of the maximum in NX (at fixed N b ) are stable. Other configurations are gravitationally unstable. The ADM is modeled as a free fermi gas and the baryonic matter through a spliced polytrope, HB as in Fig. 2. Note that the maximum mass of a baryons-only star for the same model Since the number densities are normalized by m 3 i and m b ∼ GeV, the ADM number density in most of the core at small radius is more than 10 9 times greater than typical baryon number densities. The total star mass is 1.5M and NX = 0.396(M Pl /mX ) 3 . The left-hand plots assume no ADM-baryon interactions while the middle and right-hand plots assume a repulsive ADM-baryon interaction as described in Sec. B 2 such that σ bX = 4 × 10 −49 cm 2 and 4 × 10 −47 cm 2 , respectively. The EoS used for baryonic matter is as in Fig. 6.