Condensed dark matter with a Yukawa interaction

We explore the possible phases of a condensed dark matter (DM) candidate taken to be in the form of a fermion with a Yukawa coupling to a scalar particle, at zero temperature but at finite density. This theory essentially depends on only four parameters, the Yukawa coupling, the fermion mass, the scalar mediator mass, and the DM density. At low fermion densities we delimit the Bardeen-Cooper-Schrieffer (BCS), Bose-Einstein Condensate (BEC) and crossover phases as a function of model parameters using the notion of scattering length. We further study the BCS phase by consistently including emergent effects such as the scalar density condensate and superfluid gaps. Within the mean field approximation, we derive the consistent set of gap equations, retaining their momentum dependence, and valid in both the non-relativistic and relativistic regimes. We present numerical solutions to the set of gap equations, in particular when the mediator mass is smaller and larger than the DM mass. Finally, we discuss the equation of state (EoS) and possible astrophysical implications for asymmetric DM.


I. INTRODUCTION
The nature of Dark Matter (DM) remains to a large extent a mystery. Like ordinary matter, it is possible that it consists of an asymmetric population of particles [1]. Additionally, DM may be in a cold and dense state in several regions of the Universe, for instance at the core of galaxies and dwarf galaxies [2][3][4][5] or, if captured, inside of dense and cold stars [6][7][8] or in the form of nuggets of DM [9,10]. It is even considered that DM particles can condense and make dark stars [11,12]. Much like ordinary condensed matter, it is then conceivable that DM particles manifest complex, emergent behavior at low densities, like superfluidity. This may have various implications, from the formation of proto-DM halos to the dissipationless transport of heat, to vortex formation during DM halo mergers, or simply for the equation of state of DM, which is essential for the study of compact DM object formation and evolution. Incidentally, superfluid DM is a possibility that is much studied, usually in the form of fundamental scalar DM, see e.g. [13][14][15]. Here we focus on the less studied possibility that condensed fermionic DM particles can be superfluid, see e.g. [16,17].
More precisely, we explore the possible superfluid phases of a very simple and yet quite rich asymmetric DM model. It consists of a DM Dirac fermion ψ with a Yukawa coupling to a spin 0 (φ, a real scalar) mediator. This theory rests essentially on 3 parameters on top of the DM density: the DM mass, the mass of the mediator and the Yukawa coupling. For simplicity, we neglect a possible self-coupling of the mediator. Since spin 0 boson exchange is attractive, at low densities and de-pending on the parameters, this fermionic DM may either form di-molecules that could be in a Bose-Einstein condensate (BEC phase), or Bardeen-Cooper-Schrieffer pairs (BCS phase). The thermodynamic description of the BEC and BCS phases are qualitatively different, since in the BEC phase there are true di-fermionic bound states, see e.g. [18]. But, in either cases, such DM should be superfluid at low enough temperatures. The problem turns out to be quite complex, so our aim here will be to set the ground for more phenomenological studies. In particular, we will focus on the thermodynamic properties of this simple DM model at zero temperature and will delineate which phases can be formed as a function of model parameters and DM density. This will allow us to address a question regarding the gravitational collapse of such DM, if it is accumulated at the core of neutron stars [19,20]. This requires the determination of the overall normalization of the gaps and hence necessitates a more careful treatment that goes beyond the standard, textbook BCS approximations.
Our plan is as follow. In section II we setup the model and establish which phases (BCS or BEC) DM may form at low temperatures, depending on the model parameters. To do so, we characterize the strength of Yukawa interaction at finite density using the concept of scattering length (see also appendix B). Next, in section III, we focus on the case where DM particles are in a BCS phase. We derive the consistent set of gap equations, taking into account the change of DM mass at finite density due to the formation of a scalar condensate. The gaps equations are solved numerically in section III D and, where possible, comparison is made with analytical solutions. We show that in the large DM density limit the system is generically in the BCS phase. We next discuss some possible implications for DM phenomenology, ranging from DM in neutron stars to dwarf galaxies in section IV. Finally, we draw conclusions in section V. Several technical results which may be of interest to a broad audience can be found in the appendices A -F.

II. LOW DENSITY PHASES OF A YUKAWA THEORY
Our starting point is a degenerate gas of Dirac fermions ψ (the asymmetric DM) with a Yukawa coupling g to a real scalar φ (the mediator) (1) The fermion ψ and the mediator φ are singlets of the Standard Model (SM) but ψ is charged under a global dark U (1) symmetry. Therefore µ is the chemical potential conjugate to DM fermion number N in a volume V , corresponding to a DM density n = N/V ≡ ψ γ 0 ψ . The expectation value . . . is on the ground state of the system, here taken to be at finite fermion density but at zero temperature. For a degenerate gas of free fermions g → 0, the chemical potential is equal to the Fermi energy E F , µ = E F ≡ m 2 + k 2 F and n = k 3 F /3π 2 . In eq. (1), m and m φ denote the bare fermion and boson masses at zero density. Both are modified in a medium and, in particular, at finite density. The most dramatic effect is the change of the fermion mass. In this case, physically, the scalar operatorψψ has a non-zero mean, n s = ψ ψ > 0 [21]. In turn, n s sources the scalar field due to its Yukawa interactions with the fermions δL δφ = 0 → m 2 φ φ + g ψ ψ = 0 , which consequently changes the mass of the DM see fig. 1 for a diagrammatic representation. We have expressed the fact that n s is itself a function of the effective fermion mass. 1 This effect implies that the effective fermion mass decreases with increasing fermion density. For a given value of Yukawa coupling and the mediator mass, at low enough densities, n s is close to the DM fermion number and m * ≈ m. However, as the density increases , the in-medium mass m * decreases and eventually tends to zero (m * → 0), see appendix A. 2 1 See also ref. [22,23] for an application of this effect to SM neutrino clustering. 2 Inclusion of self coupling term of the mediator in the scalar po-ψψ ψ ψ Figure 1: The tadpole is non-vanishing at finite density. This corresponds to a scalar density condensate ψ ψ which modifies the effective mass of the fermions in the medium.
The change of the fermion mass is of course not the only effect of Yukawa interactions. If the attraction due to spin 0 particle exchange is strong enough, fermions can be bound into di-molecules at low densities. These bound states behaves as bosons and can condense into a Bose-Einstein condensate (BEC) at low enough T . Nevertheless, even if attraction is weak such that no true bound state can be formed, Cooper instability can lead to the formation of Bardeen-Cooper-Schrieffer (BCS) pairs. The distinction between the two possibilities is not sharp and which of these situation is realized at low densities depends on the parameters of the theory. The transition between the BEC and BCS states as these parameters vary is continuous or a crossover [25].
In the non-relativistic limit, the nature of the low density phase can be qualitatively understood by examining the scattering length (a) of the fermions. While this is part of the standard toolbox of condensed matter physics, it is less common in the high energy physics literature. A notable and interesting exception is ref. [26], in which the scattering length is put forward as a convenient tool to discuss the properties of interacting DM. The most relevant feature of the scattering length is that, while in the Born approximation a Born 0 for an attractive potential, a > 0 signals the possibility of forming a bound state [27] and so a BEC. For convenience, we summarize the basic and relevant properties of the scattering length in appendix B.
For a Yukawa interaction, the s-wave scattering length for the singlet, spin 0 channel reads where δ 0 (k) is the s-wave phase shift. Obtaining this quantity requires solving the Schrödinger equation for the tential (λφ 4 ) qualitatively changes the picture in two ways. The equilibrium values of in-medium fermion mass gets bounded by below, i.e. m * → m √ λ in the limit of large λ. The second effect is related to the effective mass of the mediator, which would scale as m eff φ = m φ 1 + 2V ( φ )/m 2 φ φ 2 [24].
scattering problem. A much used approach, valid in the limit of large particle separation compared to the effective range of the interaction, is to approximate the attractive potential by a contact interaction, a delta function in other words, and to re-express the scattering problem in terms of the scattering length [28]. Applied to a degenerate system of fermions at finite density with Fermi momentum k F , the phases can then be characterized in terms of the dimensionless parameter k F a, with a dilute system corresponding to |k F a| → 0, while the sign indicates the phase of the system [29][30][31]. Consequently, the BEC phase corresponds to small positive values of k F a, physically corresponding to a scattering length smaller than the particle separation, while the BCS regime is characterized by small negative values of k F a. Finally, the crossover regime is realized for large absolute values of k F a. This corresponds formally to a diverging scattering length |a|, and is called a unitary Fermi gas in the literature [32].
In this work, so as to map the model parameters of the theory (1) to the possible phases of condensed DM, we go beyond the above contact interaction approximation and compute directly the scattering length by solving the Schrödinger equation. We solve the scattering problem using the numerical method proposed in [26] which is also summarized in appendix B for reference. In doing so, we can determine parameters of the Yukawa theory for which DM particles are clearly in the BCS (large negative k F a) or in the BEC (large positive k F a) phases. This result is depicted in fig. 2 as function of the dimensionless parameters β = αm/m φ with α = g 2 /4π and the ratio of length scales k F /m φ . The Born approximation corresponds to β 1 (also denoted b in the DM literature [33]), so we can expect the onset of bound state formation (and thus BEC phases) to be around β 1. However, the sign of a changes each time a new bound state channel opens, so the relation between the possible phases and the parameters is complex. The other parameter (k F /m φ ) is simply a measure of the mean particle separation over the range of the Yukawa potential, large k F /m φ corresponding to large densities. In fig. 2, the red shaded regions indicate the BCS phase, which we define to correspond to (k F a) −1 < −1, a value motivated by the results obtained based on the contact interaction approximation, see e.g. [30]. The cyan shaded regions are characterized by (k F a) −1 > 1 and are delimiting the BEC phase. The gray shaded area show the intermediate crossover phase, −1 < (k F a) −1 < 1. The unitarity limit is reached when k F |a| → ∞, indicating crossover regime at all densities, which is seen as a feature in fig. 2. Further, for finite densities, the cases of anti-resonance k F a → 0 are captured by the peaks delimiting BEC to BCS transitions.
We conclude that, for fixed k F /m φ , the BCS, BEC, and crossover regimes alternate as we change the value of β. However, for large values of β, i.e. large coupling or small mediator mass limit, the system is more likely to be in a BEC or, at large densities, in crossover regimes. Whereas, in the opposite limit, the interaction is not strong enough to form bound states and the DM fermions are expected to be in a BCS phase at low temperatures. Figure 2: Contours of (k F a) −1 . Red shaded regions are characterized by (k F a) −1 < −1 indicating BCS regime. In the cyan shaded regions (k F a) −1 > 1 indicating BEC regime, and the gray regions correspond to possible BEC-BCS crossover with −1 < (k F a) −1 < 1.
The above considerations are based on a very simple criteria, essentially the sign of the scattering length, a low momentum/short range parameter, indicative of the possible outcomes of a system of fermions at relatively low DM densities k F m. Notice that there only two possibilities for fixed parameters, either BCS or BEC at low densities with a crossover at increasing densities. This behaviour is manifest in the approximation of a contact interaction, see [30]. These considerations however neglect the fact that the effective mass of the fermion m * changes, an effect which is particularly dramatic at larger densities. To take that effect in a self-consistent way, in the next section we derive the gap equations for the BCS phase in complete generality by retaining the mediator mass and momentum dependence of the anomalous propagators. On that basis, we shall argue that, regardless of the phase at low densities, the fermions are in a pure BCS phase at large densities, corresponding to k F m * . That only a BCS phase is possible in the relativistic regime is probably due to the fact that it is not possible to form bound state of relativistic particles with a Yukawa interaction.

III. BCS PHASES
We have so far described the low density phases of the Yukawa theory with no reference to the details of superfluidity, such as the bound-state wave-functions [25]. In this section, we study in detail the BCS gap matrix ∆ ≡ ψ cψ , with ψ c being the charge conjugate of ψ, see fig. 3. We analyze its Dirac structure so as to decompose them into simple gap functions, derive the corresponding gap equations and finally solve them, taking into account their interplay with other in-medium density effects, in particular the scalar condensate and its impact on the DM mass. We aim at showing the evolution of the gaps with varying density, from the non-relativistic to the ultra-relativistic regimes, an aspect that has not yet been studied as far as we know. This implies that we will focus our study on the red regions shown in fig. 2, i.e. for small values of β, in the non-relativistic regime. However, we will argue that, for fixed β, a BEC phase becomes a true BCS phase at large densities.  The techniques to derive the gap equations for a potentially relativistic fermionic system at finite density are well-established. In this section, we first analyze the Dirac structure of the gap matrix, ∆, which is a 4 × 4 matrix (for one flavor of DM) which can be decomposed using the Dirac matrices. We follow the approach set up in [34] and in particular [35], in which the Yukawa theory is studied in detail in the massless limit, m = 0, m φ = 0. Then to derive the gap equations, we follow the energy functional approach of [36,37], (see also [38] for the Yukawa theory), which is based on the so-called Hubbard-Stratonovich transform. In particular, we were inspired by [39] to take into account the variation of the DM mass through the scalar condensate.

A. Gap Dirac structure
We will consider a general ansatz for the 4 × 4 gap matrix ∆ = ψ cψ . We work in the rest frame of the fermion gas, which is assumed to be infinite and homogeneous. In this case, one can be easily convinced that the ∆ matrix can be written as a sum of up to 8 translation invariant terms that can be expressed using the Clifford basis of matrices built upon the Dirac γ µ and γ 5 [34,35]. As the Yukawa theory preserves parity, we expect that the gap matrix is also parity symmetric. Also, the ground state is expected to be rotationally invariant. This implies that pairing of fermions should be in J P = 0 + channel. This allows us to express the gap matrix ∆ in terms of only 3 gap functions (the gaps in the sequel): The task is then to determine self-consistently these gaps, together with the scalar condensate While this is not a new problem, to our knowledge it has not been worked out in the framework of the Yukawa theory. We will show that the ∆ i 's strongly depend on n s though the effective fermion mass m * , while the dependence of the scalar condensate on the gaps is mild, see below.

B. Quasi-particle dispersion relations
The next step is to examine the spectrum of fermionic excitations (quasi-particles) near the Fermi surface, see appendix D 1 and [30]. For a free degenerate gas, the dispersion relations are simply with ω k = √ k 2 + m 2 and the subscript −(+) corresponding to fermion (anti-fermion) excitations. These dispersion relations can be read off directly from Lagrangian (1). Simply put, this means that it costs little energy to create a fermionic excitation near the Fermi surface, with v F = d /dk| k F = k F /µ and k > k F for a particlelike excitation while k < k F for a hole-like excitation. The linearity of the dispersion relation for particle/hole excitations near the Fermi surface is valid both in the non-relativistic and relativistic regimes. For antiparticles, however, the cost is at least + ≈ 2µ. If the gaps are non-vanishing, it is tedious to derive the dispersion relations but the final result can be approximated by the following fairly simple expression where again −(+) corresponds to particle (anti-particle) excitations, see appendix D 1 for details. We have assumed that the gaps are smaller than the chemical potential, which we expect to be the case in the BCS phase. As far as we could judge, our ansatz is consistent with results derived in [34], albeit with a distinct approach. It generalizes the results presented in [35], where the focus was only on the ultra-relativistic regime, m = 0 in which case only ∆ 1,2 ∆ 3 are relevant. As the gap functions are momentum dependent, in principle the dispersion relation eq. (8) interpolates between the non-relativistic and relativistic regimes. This also goes beyond the simple ansatz put forward in [30], where only ∆ 1 pairing channel is retained in the non-relativistic regime. This is consistent with setting ∆ 2 = 0, whose contribution to the dispersion relation is suppressed if k m. This leaves open the possible relevance of ∆ 3 in this regime. Notice however that its contribution is also related to a possible distinction between particle and antiparticle excitations.
This motivates the introduction of the new gap func-tions∆ We expect and will verify thatκ, which is orthogonal tõ ∆ ± − ∆ 1 , is relatively small compared to ∆ ± in the nonrelativistic and ultra-relativistic regimes. Indeed, in the relativistic limit, eq. (9) suggests that∆ ± ≈ ∆ 1 ± ∆ 2 and |κ| ≈ −∆ 3 . Whereas in the non-relativistic limit, ∆ ± ≈ ∆ 1 ± ∆ 3 and |κ| ≈ ∆ 2 . In either case,κ should be relatively small. This intuition is supported by the numerical results, see appendix E 3. So in the sequel, we simply setκ to zero.

C. Gap and scalar condensate equations
With the dispersion relation eq. (8) at hand, we now derive the gap equations in the mean-field approximation. To do so, we follow a variational approach based on the so-called Hubbard-Stratonovich transformation [40,41]. In this framework, a potential gap function, say ∆ i , is introduced as an auxiliary field in the expression for the free energy, Ω = −T ln Z/V = −p(µ, T ) at finite fermion density (and temperature), i.e. Ω → Ω[∆ i ]. Minimizing the free energy with respect to the gap function correspond to the true ground state of the system at finite density and, by the same token, a gap equation for the corresponding gap function. An extra complication is the presence of the scalar condensate n s . To take this into account, we follow the approach put forward in [39]. Alternatively, one can derive the gap equations in the mean-field approximation using the the Wick theorem, see e.g. [42], p.441. Finally, recall that n = −(∂Ω/∂µ) T (here assumed to be at T = 0).
Altogether, the gap functions and scalar condensate are thus determined imposing Thus we obtain the following equations, which we have expressed using the∆ ± combination of gap functions, .
Here, the gap equation forκ is omitted as it is expected to be subdominant in all regimes, see appendix D 2 for the complete set of gap equations. While having a relatively simple structure, these expressions call for some explanations. Consider first the equation for n s . It is easy to verify that, in the limit of zero gap, it reduces to eq. (A1) with an integral over the volume of the Fermi sphere of radius k F . Note that we have renormalized the expression for n s by subtracting the vacuum part; hence the −1 in the first term of eq. (12). We further find that n s is momentum independent, and the presence of non-zero gaps manifests as the second term in the equation. This term is proportional toκ and represents a parametrically minor correction to the gap-less expression. Now consider the equation for the gap∆ ± , eq. (13). The limit m φ k F m, corresponds to the case of contact interactions and relativistic fermions. In this case the above system reduces to the well known simple BCS gap equation [30] see also appendix D 2. Note however that in general the gap functions are momentum dependent. The standard BCS approximation assumes that the gap is momentum independent and hence constant. This requires the introduction of a cut-off on the integral over momentum. The  cut-off may be dictated by physical considerations, for instance it may be set by the Debye screening length in the case of ordinary superconductivity. Such cut-off also sets the overall normalization of the gap. In the case of Yukawa interactions, it is not clear a priori which cut-off scale can be introduced. Furthermore, we want to have a handle on the overall scale of the gap functions. In such a situation, the best way forward is to keep the full momentum dependence of the gap functions, which by construction should vanish in the limit of large momentum. So, the advantage of our consistent set of gap equations is that we can in principle integrate the above equations all the way to infinity without having the need to introduce spurious cut-off dependence. Finally note that we have not considered possible in-medium corrections to the mediator mass m φ . These corrections are O(gµ), and contributes positively to its mass. We neglect them for simplification, but one should keep in mind that for large coupling or massless mediator, they may play an important role [35].

D. Solutions to the gap equations
In this section, we present the numerical solutions to the set of gap equations, together with the scalar condensate. We use a method called matrix inversion and developed in ref. [43] (see appendix E 2). We present our findings in fig. 4, for various representative values of parameters of the theory, i.e. the Yukawa coupling and mediator mass. In the left panel we show the results for both heavy and moderately heavy mediator masses of m φ = 5 m (in blue) and 0.5 m (in red), for g = 3 (2) in solid (dashed), respectively, including the effects of scalar density condensate. For the case of heavy mediator we recover the familiar BCS solution [30], see appendix E 3. In the case of the moderately heavy mediator, we obtain solutions which are parametrically different from the BCS case. We interpret this to be due to the fact that interactions are not point-like (or contact interaction) contrary to the standard BCS approximation. In particular, we show that the gaps can be substantially larger at relatively small densities. Nevertheless, the exponential fall off, typical of BCS gap solutions, is recovered at moderate densities. At very large densities or, in other words, in the relativistic limit, we find that the solutions to the gap equations do not depend on the mediator mass. For the whole range of DM densities, we find that the gap is always substantially smaller than the kinetic energy at the Fermi surface, see the dotted line labelled F − m.
Here, we also clearly see the impact of the scalar density condensate at intermediate densities, where the solution behaves similarly to a relativistic system. For comparison, the gray dotted lines show the evolution of the gaps neglecting the change of the effective mass m * . In other words, the effect of the shift of the mass of the DM is to drive more rapidly the system into the relativistic regime.
In the right panel we present the results for light mediator with m φ = 0.13 m for the same values of couplings as in the left panel. As expected, the solution in the high density regime behaves parametrically as in the left panel. However, at low densities the situation could be drastically different from the left panel depending on the value of g. This is best understood through the low density phase diagram of the Yukawa theory put forward in fig. 2. For the case g = 3 (2), corresponding to β = 5.5 (2.4), we see that the system is in BCS (BEC) phase at low densities. As we increase the density, both the cases are in a crossover regime at k F ≈ |a −1 | = 0.07 m. Interestingly, for this choice parameters the scattering length turns out to be the same in magnitude but opposite in sign. In light of these observations, we can now understand the low density regions shown in the right panel of fig. 4. For g = 2, at low densities, we argue that the system is in the BEC phase. The solution yields the gap to be constant and much larger than k 2 F /2m. Although we get a solution for the gap, it does not represent a small perturbation to the Fermi surface, i.e. the chemical potential is no longer given by k 2 F + m 2 but should evolve towards the binding energy of the would-be DM molecules [25]. To obtain the correct solutions in this regime, we would need to simultaneously solve for both the gaps (and the scalar condensate) and for the chemical potential. This would require a more sophisticated analytical and numerical approach than the one we have considered, so we leave this particular situation for possible future works. Whereas, for g = 3, at low densities the gap is exponentially suppressed, indicative of the non-relativistic BCS phase. As we approach densities close to k F ≈ |a −1 | = 0.07 m, the system goes to the crossover regime; for which we do not present any solution and it is shown as the shaded region. Regardless, at very large densities, the system becomes relativistic (this is further enhanced by the decrease of the effective mass) and, as the formation of true bound state becomes impossible, the system makes a transition to the relativistic BCS phase, a feature which seems to be novel.

IV. ASTROPHYSICAL CONSTRAINTS AND IMPLICATIONS
DM self interaction constraints. -In the context of DM self interactions, in fig. 5, we show the phase diagram of the Yukawa theory in the k F /m − m φ /m plane for a dark matter candidate of mass m = 1 GeV for g = 3 (left) g = 2 (right). The gray dotted line corresponds to k F = m φ . We overlay the constraints on the dark matter self-interaction cross section at the scale of the Bullet cluster, requiring σ/m 1 cm 2 /g at a velocity of v = 2000 km/s. For such DM candidate, we find that a m > 20 is excluded. The corresponding excluded mediator mass range m φ is shaded in gray. As could be expected, the unitarity regime are excluded, as is most of the very light mediator regime. The very fine viable intervals of m φ correspond to vanishing self-interactions, i.e. a → 0. This nicely illustrates the possibility for a dark sector to manifest emergent phenomenon like superfluidity, while being not entirely excluded by self interaction constraints on DM.
DM in neutron stars. -We now turn to discuss the impact of superfluid gaps on the fate of captured asymmetric DM in the core of neutron stars (NS). More precisely, DM particles in the galactic halo can be efficiently captured by NS [6,44]. For any given neutron star of mass M and radius R, the maximal rate of DM accretion is given by the so-called geometric rate corresponding to scattering cross section πR 2 /N ∼ 10 −45 cm 2 , where N is the total number of target neutrons in the NS. The maximum number of accumulated DM in the NS over its life time of 10 Gyrs is given by After DM thermalization with the medium [7,[45][46][47], the DM thermal sphere can collapse to a black hole. For non-interacting fermion DM particles, this happens as soon as the number of accumulated particles exceeds the Chandrasekhar limit given by N ch ≈ M 3 pl /m 3 . In ref. [19], it was suggested that in the presence of attractive self-interactions, such as in the Yukawa theory, the above Chandrasekhar limit could be parametrically reduced to ∼ (m φ /m √ α) 3 N ch . Soon after, ref. [20] pointed out that as the DM thermal sphere shrinks the DM density increases therefore the emergence of scalar density condensate should be accounted for consistently in the calculation. Therefore, DM particles become relativistic at smaller densities than in the free case and the Chandrasekhar limit is not parametrically modified.
In this section we comment on the impact of superfluid gaps on the above scenario. To this end, the order parameter ∆ should be evaluated. As an example, we choose the model parameters m = 200 GeV, m φ = 1 MeV, and α = 10 −3 , as exemplified in ref. [19]. We find that at low densities, for these parameters, the system finds itself in the BEC phase, and consequently the gaps are much larger than the chemical potential. However, at higher densities, fermions become relativistic and the system transitions to the relativistic BCS phase. At this stage we can estimate the values of the gap using our relativistic BCS analytical approximation ). This yields the numerical value of the gap to be exponentially suppressed and thus negligible. While the BEC phase may change the equation of state at low DM densities [25], we conclude that at large densities, in the BCS phase this has little impact on the equation of state of the DM with a Yukawa interaction and thus the conclusions of [20] are essentially unchanged.
Condensed DM Halos. -The possibility that DM may be in a degenerate fermionic gas phase and form a core at the center of galaxies and dwarf galaxies has been put forward in the literature, see e.g. [3][4][5], albeit for noninteracting particles. More recently, possible superfluid phase for DM in the context of dwarf galaxies was studied in ref. [17], in a model of ultra-light dark-QCD matter parameterized by rescaling the known QCD CFL phase [48]. It is well known that the Yukawa theory shares qualitative features at finite densities with quark matter [35]. We now consider the possibility that DM particles that make up the halos of galaxies (e.g. dwarf spheroidal galaxies) could manifest emergent phenomena described in our work, in a simpler model such as the Yukawa theory. We first present the equation of state (EoS) for the Yukawa theory 3 and then showcase a possible realization of condensed DM halos at galactic scales.
In the left panel of fig. 6 we present the EoS in dimensionless units, for the case of free degenerate fermions (dashed), and for attractive self-interactions (solid), see also appendix F for further details. The strength of interaction in this case is most conveniently characterized by the parameter C 2 φ = 4/3 αm 2 /(3π 2 m 2 φ ). The curves shown in red, cyan, and magenta, correspond to values of C 2 φ = 3, 10 and 100, respectively. At relativistic densities (ρ m 4 ), each of the EoS obtained in the presence of attractive interactions reproduce that of the relativistic Fermi gas, P = ρ/3. However, in the non-relativistic regime (ρ m 4 ) significant departures from the free case appear. Due to the emergence of the scalar density condensate, at large values of C 2 φ , m * m, and hence the relation P ∼ ρ is valid down to densities ρ ∼ m 4 * . Additionally, the presence of attractive self-interactions lead to the appearance of a critical saturation density analogous to the nuclear saturation density, at which P vanishes for a finite ρ. Close to the saturation density, the EoS is very stiff and the system becomes incompressible. Interestingly, for such densities, the pressure of the interacting gas is much smaller than that of the free case. This would lead to more compact self-gravitating configurations of the DM cloud, when compared to the free case. It is also seen that, as the system contracts and ρ becomes larger (for example during its evolution towards gravitational collapse) the pressure quickly rises above the values of the free case, hence hindering collapse, as described above.
In the right panel of fig. 6, we present the mass-radius relationship for a few representative cases depicted in the left panel. We have obtained these configurations by solving the Tolman-Oppenheimer-Volkoff equations with the EoS presented in the left panel. Particularly, we have considered the free degenerate EoS (C 2 φ = 0) and the mildly interacting EoS (C 2 φ = 3). As is well known, for free degenerate fermions M ∝ R −3 [3,4,49], and such DM candidates cannot simultaneously reproduce typical DM halos of both dwarf spheroidal and large galaxies. However, if DM has sufficient attractive self-interactions, as exemplified in the Yukawa theory, stable solutions exist for incompressible configurations which behave as M ∝ R 3 . An estimate of the M − R relationship is given by As shown in the figure, such configurations could lead to halo masses and sizes that are both typical of dwarfs (blue shaded regions) and large galaxies (green shaded regions). Due to their incompressible nature, both halos would exhibit cored profiles. It is interesting to note that perhaps only the inner cores of halos are in such configurations, while the outer regions could be at temperatures high enough for the scalar condensate to not play any role. A comprehensive analysis of such phenomena could, for example, be captured by a piece-wise EoS correctly interpolating between the cold, compact core and a hot, less dense outer halo. Light fermionic dark matter candidates are clearly challenging to reconcile with the formation of cosmic structures. For m < 2 keV, the Fermi pressure is itself sufficient to modify substantially the matter power spectrum [5,50]. Nonetheless, if the DM particles are bound in di-fermionic molecule, for example in the BEC phase, at early time, such bounds could be alleviated.

V. CONCLUSION AND PERSPECTIVES
We have developed a formalism to study possible phases in the Yukawa theory. To this end, we have brought together concepts originating from the theory of non-relativistic scattering, the physics of high density nuclear matter, and phenomena that are realized in condensed matter systems, and have presented our findings in the context of particle dark matter at finite density. By computing the scattering length in the Yukawa theory we have found that DM could be in BCS, BEC or crossover phase, depending on the parameters of the theory at low densities. We have further explored the BCS phase. To do so we have computed superfluid energy gaps, including the effects of the scalar density condensate; by retaining the mediator mass dependence and deriving the consistent set of UV finite gap equations for the three most dominant pairing channels. We have found that at large densities (or ultra-relativistic limit) DM finds itself in the BCS phase, irrespective of the phase realized at low densities.
As an application of our formalism we have examined whether captured asymmetric fermion DM in the core of neutron stars can collapse to a black hole in models with attractive DM self-interactions, for DM mass ∼ GeV. Even in the presence of large attractive interactions, it was pointed in ref. [20] that collapse to a black hole is impeded due to the emergence of scalar density condensate, and that the Chandrasekhar limit remains the same parametrically as that for free degenerate fermions. In this work, we have studied the interplay between the scalar density condensate and superfluid gaps consistently. We have found that further accounting for the superfluid gaps does not qualitatively change the conclusions drawn in ref. [20], since the gap contribution to the pressure is small and positive, with scaling ∝ µ 2 ∆ 2 .
By consistently computing the EoS in the Yukawa theory, we have also explored the possibility that halo DM could realize cored density profiles due to emergent density effects, at both small and large scales, corresponding to dwarf galaxies and other large galaxies. We have found that the impact of the scalar density condensate is important and dominates over that of superfluid gaps.
There are various directions for further phenomenological enquiry. In this work we have focused on identifying the possible phases of condensed DM with a Yukawa interaction, on self-consistently determining the gaps and the scalar condensate and the impact of the latter on the energetics of the system. The next natural step would be to examine the transport properties of such systems.
On astrophysical scales high density regions of DM are found, e.g. in dwarf spheroidal galaxies, or at the center of more massive galaxies. Examining velocity dispersion data of several Milky-Way dwarf galaxies would in principle provide a window to determine the phase of DM. Moreover, since superfluidity is realized through breaking of the global U (1) we expect formation of strings and vortices. Observations of galaxy-galaxy mergers and/or dwarf galaxy-host galaxy mergers might provide a way to probe DM phases. Finally, depending on how dark sector couples to the visible sector dark stars or hybrid stars could exist in the Universe. Gravitational wave observations of mergers of such objects with black holes, neutron stars, or white dwarfs, would pave the way in testing the scenario through waveform template fitting. The scalar condensate at finite density n s = ψ ψ has peculiar properties [21]. In particular, it leads to the vanishing of the effective fermion mass at large densities.
A simple application of the decomposition of ψ in terms of creation/destruction operators gives for k F m, n s ≈ n = ψ γ 0 ψ = k 3 F /3π 2 while for k F m, n s ≈ mk 2 F /2π 2 . This scalar condensate leads to a modification of the fermion mass through m * = m − g 2 /m 2 φ n s (m). Replacing m by m * in the argument of n s gives a self-consistent, mean-field approximation for m * . Alternatively, one can start from the grand partition function of the Yukawa theory eq. (1) where Ω is the (Landau) thermodynamical potential Ω = −p(T, µ). At non-zero temperature β = 1/T , integrating over the fermions and replacing the field φ by a constant mean-field value φ 0 gives One recognizes contributions from fermions and antifermions. Note also, ω = ω k (m * ) with m * = m − gφ 0 . In the limit T → 0, the integral reduces to the first term which can be written as There is an infinite contribution from vacuum which can be renormalized away, see below. Then, minimizing Ω with respect to φ gives Within the same mean-field approximation, the energy density of the Yukawa gas is written as with the free gas contribution * and the Fermi energy E * F = k 2 F + m 2 * . Now we examine the energy density above. In the non-relativistic limit, k F m, this gives Notice that a factor 1/2 comes from the 1/2m 2 φ φ 2 0 term. This reproduces the leading term derived in [19] from the potential energy of a sphere of N particles with a short range attractive interaction, here R is the radius of the sphere, F G is the free Fermi gas energy density, and α = g 2 /4π. At very small k F , the energy density ∝ n. If gm/m φ is large enough, it becomes negative. However (at large densities) when m * approaches zero. The energy density is then as a relativistic gas, We conclude this appendix with a brief note about the vacuum contribution to n s . The above discussion is for a renormalized n s so that n s = 0 for T and µ = 0. The vacuum contribution is potentially infinite being given by for some cut-off scale Λ 2 . This leads to a negative n s and so a positive contribution to the fermion mass. Following the logic of above leads to the NJL mechanism of chiral symmetry breaking, in which case the effective dressed mass is where G is some effective coupling (dimension 1/M 2 ). The chiral limit is m = 0 and a scalar condensate contributes dynamically to yield a non-zero fermion mass, m * ∝ ψ ψ 0 . At finite density, the contribution to n s is positive and non-zero if and only if m = 0, i.e., the opposite to the vacuum case so to speak.

Appendix B: Scattering in Yukawa theory
The scattering length is a basic quantity used to study the onset of BEC and, more generally, low energy scattering processes. It is less commonly encountered in DM literature, see however [26] for applications to DM self-scattering in galaxies and dwarf galaxies. It is particularly useful to characterize in simple terms complicated, possibly strongly interacting scattering systems in terms of a physical observable such as the scattering length.
The wave-function in the CM frame of two non-relativistic particles is of form with scattering amplitude f (θ) and differential cross section dσ/dΩ = |f (θ)| 2 . In the case of s-wave scattering, which is dominant for many systems (e.g. a pair of identical fermions in a singlet spin configuration), the scattering amplitude is constant f (θ) = f 0 . In terms of the s-wave phase shift δ 0 , Then, for low energy scattering k → 0, the wave function can be written as where the scattering length (a) is given by, In this limit σ = 4πa 2 . For a > 0, the vanishing of the wave-function signals the possibility of bound state formation [27]. For instance, for a neutron-proton (n-p) pair in a spin S = 1 triplet state, the scattering length is a ≈ +7 fm corresponding to the formation of a deuteron bound state. There is no n-p bound state with in the spin singlet state S = 0, but the s-wave scattering is very large (∼ resonant) at low energies, corresponding to a ≈ −20 fm. A resonance corresponds to a phase shift δ 0 (k res ) = π/2 and through eq. (B4), to a change of sign of the scattering length. Formally, the scattering length is infinite if δ 0 → π/2 as k → 0. In the condensed matter literature such systems are called unitary Fermi liquids. Similarly, for a neutron-neutron (n-n) interactions in the s-wave, a ≈ −17 fm; the nuclear interaction between two neutron is attractive, but not large enough to form a real 2-body bound state.
The relevance of all of this is that an attractive interaction, for instance through a Yukawa coupling may, depending on the parameters of the theory, lead or not to formation of bound state at low energies. Thus corresponding to a BEC or a BCS state, respectively, see e.g. [30]. In the following we will reproduce the results of scattering length in Yukawa theory by solving the Schrödinger equation (following ref. [26]) and use these findings to study phases of Yukawa theory.
The interaction Yukawa potential in co-ordinate space takes the form where α = g 2 /(4π) and m φ is the mass of the force mediator. The phase shifts due to DM self scattering are obtained by solving the Schrödinger equation for the radial wavefunction (R l,k ) through the equation with boundary condition rR l,k = 0 at r = 0. As discussed in detail in ref. [26], with change of variables, the above equation could be re-expressed in terms of spherical Hankel function of the first kind h l . This results in the following equation for the phase shift which is to be solved with boundary conditions δ l,k (0) = 0 and δ l,k → δ l at r → ∞. Having determined the phase shift, the s-wave scattering length is obtained by For comparison we also compute the scattering length for Huelthen potential. The interaction potential in coordinate space in this case reads with δ = 2ζ(3) m φ [51]. The s-wave scattering length for this potential can be derived by solving the Schrödinger equation as above. However, for the Huelthen potential analytical solution exists for the scattering length and the effective range of the interaction [26,51,52], which are given by Here, η = α m/δ and ψ (n) (z) are polygamma functions of order n and γ is Euler-Mascheroni constant. In fig. 7 (left) we present the scattering length in units of the mediator mass (m φ ) as a function of dimensionless variable β = αm/m φ , for Yukawa (Huelthen) potential in solid (dashed) lines. Blue (red) color indicates instances where the scattering length is positive (negative). As expected, we find that the Huelthen potential could be a good approximation for the Yukawa potential. In fig. 7 (right) we show 1/(k F a) as a function of α, for the case m φ /m = 0.13 with k F /m = 10 −2 , for Yukawa (Huelthen) potential in black (gray) colors. We conclude that the Huelthen potential is qualitatively accurate enough to describe the unitary limit.

BCS-BEC crossover
We end this appendix by recalling the well known phenomenon of BCS-BEC crossover in the limit of contact interactions. This transition is best parameterized in terms of the s-wave scattering length of the theory in the nonrelativistic limit at finite density. For such a system, within the constant gap approximation, the equation for energy gap and the chemical potential for given a scattering length should be solved simultaneously [30]. The gap equations take the form with There two classes of asymptotic solution to eq. (B12) depending on the sign and magnitude of the dimensionless parameter k F a. For large and negative values of (k F a) −1 , the chemical potential µ is close to the Fermi energy and the gap is exponentially small as could be expected in the BCS theory. and their condensation, often referred to as the BEC phase. While the limit (k F a) −1 → 0 corresponds to the unitary limit, where the phase of matter is neither in BCS nor in BEC phases, called the crossover regime. Quantitatively, (k F a) −1 < −1 delimits well the BCS phase from the crossover. The BEC and crossover regimes are similarly delimited by the inequality (k F a) −1 > 1 [30,53]. The BEC-BCS phase transition is illustrated in fig. 8. In theories with light mediators the interactions can nevertheless be considered to be short-range at low enough densities. Analogous to contact interactions, the scattering length in these theories can be used to understand the phase diagram of nonrelativistic systems.

Appendix C: Condensation in the Yukawa theory
Cooper pairing according to BCS theory in metals results from an extremely small effective attractive interaction between electrons sourced by background lattice of positively charged ions. In weakly coupled theories Cooper pairing is a collective phenomena in which the Cooper pairs are spatially separated at distance scales larger than the average separation between the constituent fermions of the system. In contrast to di-fermionic molecules that are bosonic, physically Cooper pairs are thought to be quasi-particles that could be called a boson. Formation of Cooper pairs could also be sourced by fundamental attractive interaction between fermions. To this end, we consider a simple model where fermion (ψ) interacts by the exchange of scalar boson (φ) through an attractive Yukawa potential, given by where Γ = I (iγ 5 ) for scalar (pseudo scalar) couplings. The chemical potential is denoted by µ and g (> 0) is the coupling constant. The bare fermion (boson) mass is denoted by m (m φ ). The fermion ψ and the mediator φ are singlets under the SM gauge group. However ψ is charged under a dark global U (1) symmetry conserving dark fermion number. The spontaneous breaking of this global symmetry results in modes that propagate with energies smaller than Fermi energy, the so-called Cooper pairs. Consequently the system exhibits superfluidity (and not superconductivity) at finite density and temperature [35]. The Dirac structure is analogous to single-flavour, single-colour limit of QCD. This simple structure already captures essential features and dynamics of superfluidity that could arise in models that have particles charged under more non-trivial groups [30].
The study of fermion-fermion condensate (gaps) and fermion-anti-fermion condensate (scalar condensate) is well established through the computation of partition function Z(µ, T, V ), within the mean-field approximation. The partition function and the action reads with the inverse propagators defined as D −1 (x, y) = δ 4 (x − y) ∂ µ ∂ µ + m 2 φ and G −1 0 (x, y) = δ 4 (x − y) i / ∂ + γ 0 µ − m , for scalar mediator and fermionic particles, respectively. We follow the standard finite temperature field theory notation for space-time integral, abbreviated as x ≡ 1/T 0 dτ d 3 x [54]. The model presented here is reminiscent of the so-called σ − ω model of dense nuclear matter interactions [21]. Broadly speaking, protons and neutrons feel a long-range attractive force, arising from the exchange of a scalar meson (σ), and a short-range repulsive force, due to exchanges of a vector boson (ω). Thus, eq. (C1) corresponds to the case of purely attractive interactions, i.e when the ω meson plays little role. The features of our "σ model" often shares qualitative similarities that appears in the description of formation of large bound states of nuclei which appear beyond a particular value of saturation density, and with molecular interactions via van der Waals forces.
We proceed by first integrating out the mediator (φ). Since superfluidity is the phenomenon at the surface of the Fermi-sea the most relevant scale in the problem is identified to be the Fermi-momentum (k F ). Before we begin, however, it is instructive to ponder whether it is justified to integrate out the mediator also when m φ is the smallest scale in the problem (i.e. m φ k F ). The answer is provided in [55]. For completeness we reiterate the arguments here and specify by what we mean by integrating out the mediator in the context of superfluidity, and the resulting computational advantage. The partition function in eq. (C2) sums over all possible configurations for the fields. Shifting the scalar field arbitrarily does not change any physical quantities. For example, we can shift the scalar field the following way: After integrating over all the field configurations of Dφ we get the new effective action for the fermions The effective approximation is reasonable to describe the formation of cooper pairs since the mediator φ is never produced on-shell around the Fermi surface. In other words the momenta exchanged between fermions that would form a cooper pair is always in the t-channel. Note that s-channel diagram is absent due to the presence of a finite chemical potential. The scattering is causal and space-like (q 2 < 0). Thus, for extremely large m φ k F , m the bosonic propagator in momentum space scales simply as m −2 φ . Whereas, in the limit m φ k F , m the propagator scales as t −1 . The exchange momenta can never be much larger than the Fermi momentum due to degeneracy, as the fermions very close to the Fermi surface effectively participate in scattering. In this work we keep the momentum dependence of the bosonic propagator explicitly and examine in detail how the superfluid gap depends on m φ and momentum of the gaps. The above procedure clearly brings the fermion bilinears in quadratic form, thus easing the computation of partition function. Before we get into details in the next section, we first examine the possible combination of bilinear that can get expectation values.

Scalar density condensate
The number density of ψ is conserved. Therefore the mediator φ (scalar) will couple to the number density. Diagramatically this is shown in fig. 1. In terms of Wick contraction, we have where the contractions are taken one at a time. This in effect modifies the free fermion propagator through the correction to the bare fermion mass. Thus the in-medium mass of the fermion could be written as m * = m − Σ, with where the minus sign appears due to the fermion loop (see fig. 1). The self energy Σ is related to the scalar density through Σ = g 2 n s /m 2 φ . Finally, note that if the mediator were to interact with pseudo scalar couplings to the fermions (i.e. Γ = iγ 5 ) Σ is identically zero. Thus, we do not expect any in-medium correction to the bare mass for pseudo scalar interactions.

Fermion-Fermion condensate: Cooper pairs
As mentioned above, in the presence of attractive interactions (however small), the Fermi surface is unstable and generically leads to mixing between particles and holes at the edge of the surface that manifests as Cooper pairs. Diagramatically, this is shown in fig. 3, where the blob represents the gap function. In contrast to fermion-anti-fermion condensate (scalar density condensate), Cooper pairs in a superfluid are condensates of fermion-fermion pairs. The BCS contribution to this condensate are obtained by the following Wick contraction of the 4-point function in eq. (C4), Thus the object that will get an expectation value is of the form ψ(x)ψ(y), i.e. fermions at two different co-ordinates behave collectively in the same way. Clearly, the product ψψ (orψψ) is not well defined. With the introduction of charge conjugate spinor (ψ c ≡ Cψ T ), the Cooper pairs can now be written as ψψ c , where C = iγ 2 γ 0 . In the following appendix we describe the Dirac structure of the gap and derive the gap equations.

Appendix D: The consistent set of gap equations
We outline the derivation of the set of consistent gap equations from the free energy (Ω) starting from the effective action eq. (C4). The four-fermion interaction part in eq. (C4) describes all possible 2→2 scattering processes involving the fermions. Due to the biquadratic nature of this term we cannot perform integrations over ψ andψ to obtain the partition function. One possible way to correctly describe the superfluid ground state of the system is to approximate the the above biquadratic term as a product of a fermion bilinear and a fermion condensate, the so-called mean field approximation. To this end, we rewrite the four-fermion interaction in eq. (C4) in terms of the usual spinor ψ and its charge conjugated counterpart ψ c (x) = Cψ T (x), as described in [30,35]. In the mean field limit, the fermion bilinears ψ c (x)ψ(y) (ψ(y)ψ c (x)) andψ(x)ψ(x) (ψ c (x)ψ c (x)) are approximated by their expectation values and small fluctuations around their expectation values, respectively. Additionally the mean fields are assumed to be static. Assuming D(x, y) = D(y, x), we havē The following fermion bilinears can get expectation values The parametric and function forms of Σ and ∆ + have been computed individually in several works, most notably in refs. [20,35]. In this work we want to consistently capture condensation in multiple channels (see below). In the mean field approximation, we capture one-particle irreducible (1PI) diagrams. The most general ansatz for the structure of ∆ in the rest frame of the medium must be consistent with Grassmanian nature of the fermions, i.e. they should not vanish upon anti-symmetrization. Depending on the model at hand, parity and helicity impose restrictions on the allowed gap structure. In general, the ansatz can contain up to eight terms (expressed in the basis of linear combination of γ matrices) that are translationally invariant [34,35]. Assuming parity is not violated and that the condensation is in J = 0 + channel, one possible ansatz for the fermion pair Dirac structure (that would dominate) and the scalar condensate would have the form Of the possible eight structures only four structures apply to the case of scalar mediators as they are parity conserving. As we have only one type of fermionic particle ψ in the fundamental theory, only three gap structures are most significant, and contribution from ∆ 7 is identically zero.
As mentioned before, we include all channels of condensation within a single framework. To this end, we suitably modify and closely follow the imaginary time coherent state path integral formalism as presented in [39,56,57] and simultaneously include all possible channels of condensation. Having identified the relevant gap structures we introduce fermionic (bosonic) auxiliary fields Φ ± αβ (Σ αβ ) that captures the degrees of freedom of the low energy theory via Hubbard-Stratonovich transformation (HST).
We introduce Hubbard-Stratonovich auxiliary fields ∆,∆ and ρ(≡ n s ). The transformation is carried out by inserting two 'fat-identites' into the partition function. The first for the scalar condensate being For the fermionic HS fields the following Gaussian integral is introduced Rewriting eq. (C1) in terms of the auxiliary fields results in the following interaction Lagrangian in co-ordinate space Notice that now the Lagrangian (D6) is quadratic in ψ, φ thus enabling integration over the fundamental fields ψ, φ. In other words we can trade integration over fundamental fields for an integration over auxiliary fields. In other words, setting the auxiliary fields to their vevs (i.e. the mean-field approximation) allows for computing the path integral of a theory of free particles, albeit with anomalous (i.e. off-diagonal) propagators. To this end, we go to momentum space by Fourier transforming the fields that are in co-ordinate space. We now treat the action to be a functional of charge conjugated spinor ψ c as well (cf. [30]), resulting in the following action in Nambu-Gorkov space (i.e. 8 × 8 matrix) with the new spinors defined asΨ The inverse propagator reads We now have all the ingredients at hand to write down the free energy of the theory through the partition function that reads Inserting our ansatz for the gap structure eq.(D3) in the above and exploiting the fact that the determinant of a matrix is the product of its eigenvalues ( i (k)), the Helmholtz free energy reduces to the form We note that at this level the difference between the contribution of the scalar condensate and the gap functions is factor four. By setting the HST fields to their expectation values a stationary phase analysis can be performed [36,37] and the dynamics of the low energy theory that describe condensed DM can be examined. The advantage of HST is now manifest. The expression for the free energy of the system above includes the contribution of 1PI diagrams only. The contribution of Ω is UV finite upon subtraction of the free energy of fermion in the non-interacting theory. The scalar density condensate contribution Ω Σ is also finite owing to their momentum dependence. However, the contribution from the gaps Ω Φ are UV divergent which implies divergent pressure. To remedy this problem one should consistently include 2PI contributions [48], which is the CJT formalism. Upon inclusion of these diagrams the free energy contribution of the gaps reads

Dispersion Relation
The roots of the determinant of the inverse propagator are four fold degenerate. The dispersion relations take the form (D17) The effective mass is denoted by m * = m − Σ(0) and the energy is ω 2 = m 2 * + |k| 2 . The dispersion relation for particles (anti-particles) corresponds to − ( + ). In the above expression all the gap functions are momentum dependent. Admittedly not much physics can be extracted from the above expression. However, we know that the gap functions are not larger than the chemical potential upon including the momentum dependence in the energy gaps [35]. We make use of this fact and expand the dispersion relation in ∆ i / √ µω and a-posteriori check that they are consistent. The approximated dispersion is of the form We further note that the last two terms approach zero as m * → 0 or |k| → 0. We believe neglecting the last two terms is nevertheless a good approximation as all the relevant dynamics are encoded in the first few terms. This results in the following much simpler dispersion relation For convenience we introduce the parameters As stated in the main text, when m * → 0, the gap∆ ± ≈ ∆ 1 ± ∆ 2 and |κ| ≈ −∆ 3 . In the non-relativistic limit as k → 0, the gap∆ ± ≈ ∆ 1 ± ∆ 3 and |κ| ≈ ∆ 2 . The gap parameter∆ − is the same as the order parameter d obtained in Bailin and Love, see eq. (3.54a), evaluated at the Fermi momentum k F .

Gap equations and scalar condensate
We use a variational approach to derive the gap equations. Upon setting the auxiliary fields to their expectation values we minimize the free energy eq. (D15) with respect to the gap parameters, i.e. non-trivial energy gaps would be the solution to The gap equations take the form Few comments are in order: We need to perform the Matsubara sum in k 0 = i(2n + 1)πT . As we are interested in stationary solutions, we suppose that Σ and ∆ depend only on three momentum k. Furthermore, as in Pisarski and Rischke [35], we assume that the exchanged boson has zero energy, i.e. p 0 = k 0 = 0. This removes the k 0 dependence of the gaps. This approximation is justified, as k ∼ µ. For large energy exchange of bosons, the fermions lie far away from the Fermi surface (which implies that it is too expensive to have an energy gap). Note that, these considerations are reminiscent of Eliashberg theory of superconductivity where the energy dependence of the gap is retained while the momentum dependence is neglected [58]. The final set of gap equations in the T → 0 limit reads .
Note we have renormalized the expression for the scalar density condensate by subtracting the vacuum part (hence the −1). We find that the condensate contribution is momentum independent. The advantage of having a consistent set of gap equations is that we can integrate the above equations all the way to infinity without having the need to introduce spurious cut-off dependence.
The BCS limit of the gap equations We now consider the limit m φ k F , m, i.e. when the mediator mass is the largest scale. We aim to recover the well-known BCS equation from the above set of equations. In the relativistic limit (m * → 0), we find that the gap channels ∆ 2 and ∆ 3 being identically zero in the heavy mediator, relativistic limit. Within the BCS approximation, the gaps are assumed to be momentum independent. Such an assumption signifies possible UV divergence for the gap equation. In the non-relativistic limit (ω → m), we find that with∆ + being identically zero. In this limit, two pairing channels ∆ 1 and ∆ 3 dominantly contributes to∆ − . The gap equation for∆ − is, up to a factor 2, the BCS gap equation. Note thatκ retains an overall momentum dependencẽ κ ∝ p/m, hence it is suppresed compared to ∆ − .
The above equations were evaluated in the limit T → 0 corresponding to temperatures well below the critical temperature of phase transition. More generally, the gap equations including the dependence of finite non-zero temperature takes the form the form In the above we have introduced a cut-off momentum k 2 c = 2m(E c + µ). This way the rhs of the renormalized gap equation converges as k c → ∞ and E c µ. In obtaining the above we have also assumed m∆ k 2 F and m∆ k F k 2 F , with k F . Thus, the solution to eq. (E1) is Adding the condensate contribution does not bring about further complications. The Fermi momentum is defined as k 2 F /(2m) = µ + g 2 /m 2 φ n, with fermion number density n. For this case, one finds a parametrically similar solution as above when we proceed following eq. (B12).
The solution in the relativistic limit is similar and reads (E4) These ansatz are exactly applicable only in the case of strictly contact interactions. For all other cases the full set of gap equations should be solved numerically.

Structure of momentum dependent gap equations and their numerics
Momentum dependent gap equations are of the form with the kernel function denoted by K(p, k ) and the dispersion relation written as − (k ). The task is to solve for ∆(p).
In mathematics, the above equation goes by the name non-linear Fredholm equations of the first kind (NLFEFK). It is well known that a direct brute-force solution to the above equation is very sensitive to the initial guess value and consequently unreliable.
In the context of superfluidity in high density neutron matter where such energy gap equations are often encountered, new techniques have been proposed to overcome the numerical difficulty imposed by eq. (E5). One such technique is discussed in refs. [43,60,61], where eq. (E5) is rewritten in the form of non-linear Fredholm equations of the second kind (NLFESK), which is more tractable and offers numerical stability in all regimes without making severe physical assumptions.
NLFESK take the form where G(x, t) is the kernel function and the non-linearity is hidden in the function H(y(t)). In the context of gap equation the first step is to recast eq. (E5) in the form eq. (E6), as suggested by ref. [43]. The idea is that the kernel function K p,k ≡ K(p, k ) can be written as a sum of variable separable parts close to the Fermi surface (i.e. at p = k = k F ), and a remainder part that is non-zero only away from the Fermi surface (i.e. the excitations/deformations of the sphere). More strictly K p,k ≥ 0 or K p,k ≤ 0 ∀ p, k ∈ [0, ∞] (as they are our integration limits). Furthermore, F (x), y(x) : [0, ∞] → R. These conditions seem to be trivially satisfied by the gap equations we have. The ansatz can be written as where the separable function is φ, the remainder function is W p,k that vanishes on the surface, and K (k F ,k F ) = K(p = k F , k = k F ). Using the above two conditions in eq. (E7), we obtain the following Thus the gap equation eq. (E5) becomes then we have To obtain the last equation (constraint equation) we have set p = k F in the reformulated gap equation above. As suggested by refs [43,60,61], it is numerically convenient to divide both sides by the normalization of the gap at Fermi momentum and recast the equations in terms of dimensionless variables, by introducing the so-called shape function y(p) = ∆(p)/∆(k F ),

Results
As mentioned previously, we solve the full gap equation keeping the momentum dependence using techniques described in ref. [43,60] (see above). We systematically present results for the full set of gaps below. We begin with the discussion of the solution of gap equations in the absence of scalar density condensate. These results are shown in fig. 9. All but the right panels are the same as in fig. 4 (see gray dotted curves) but where the effect of the scalar density condensate is neglected, i.e. Σ = 0. In the top left panel we show in this scenario the numerical results for heavy mediator masses m φ = 5 m (blue curves) and moderately heavy mediator masses m φ = 0.5 m (red curve), for g = 3 (2) in solid (dashed) respectively, neglecting the effect of the scalar density condensate. In blue dash-dotted (dotted) is given, for g = 3 (2), the following analytical ansatz for which good agreement is found with the numerical results in the heavy mediator mass regime. The first exponential factor was derived in ref. [35] under the assumption of massless fermions, the second exponential factor is the well-known BCS ansatz. In the top right panel, we show the solutions for the gaps for each individual pairing channel, as well as along the two orthogonal directions,∆ and κ. As discussed in the main text (section III), we find that not all the pairing channels contribute at all densities. At small densities∆ − is the largest gap and it is dominated by the combination of pairing channels ∆ 1 − ∆ 3 . At large densities however∆ + dominates through the combination of pairing channels ∆ 1 + ∆ 2 . These results are consistent with expectations, even though the transition from one regime to the other had not been worked out so far. Also, as anticipated, at all densities κ is subdominant, ∆ 3 (∆ 2 ) being small in the relativistic (resp. non-relativistic) limit. Finally, we note that the dominant gap combination,∆, smoothly evolves as the system changes from the non-relativistic regime to the relativistic regime. The bottom panels focus on the light mediator case with m φ /m = 0.13. In the left panel we present the gaps evaluated at k = k F for g = 3 (2) in solid (dashed). In the right panel we show the momentum dependence of the gap∆ − (k), the so-called shape function, at representative values of densities, k F = 0.01, 1, 10, represented by markers. At large momentum (k k F ), the momentum dependence of the gap is found to be ∼ k −1.85 , ensuring self-consistent UV convergence of the gap equations. At low densities, the gap can exhibit one (or multiple) change of signs, and this is exemplified by the dip of the shape function for the lowest shown density, k F = 0.01. For such choice of parameters, the gap is negative at large momentum, but nevertheless goes smoothly to 0 as k → ∞.  The effective mass m * as a function of density is shown in fig. 10, for the same model parameters considered in the above figures and in the main text. We follow the same color scheme. The non-relativistic regime k F /m < 1 and relativistic regime k F /m > 1 are delimited by the vertical gray dashed line. For the lowest effective coupling C 2 φ , i.e. comparatively large m φ /m and low g, the effective mass m * deviates from the bare mass m in the relativistic regime. The effect of the scalar density condensate on the gaps is thus expected to be small, see the blue curves along their gray dotted companion in the right panel of fig. 4. At low densities, as m φ /m becomes smaller than unity, for the studied choice of the coupling, large impact on the gaps due to the scalar density condensate are expected. This is in particular shown by the differences between the red curves along their gray dotted companion in the left panel of fig. 4. Similar effect is also seen in the right panel of fig. 4.
In fig. 11, the focus is put on the dependence of the shape functions on the scalar density condensate for light mediator masses. The bottom left panel is the same as the right panel of fig. 4, the markers represent the densities at which the shape function is given on the bottom right panel. As before, the gray dotted lines correspond to the case of no scalar density condensate, Σ = 0. For g = 3 and m φ /m = 0.13, due to the change in the effective mass, at moderate density k F /m = 1, the system is however relativistic as k F /m * 1. In turn, the shape function for k F /m = 1 closely resembles that for k F /m = 10, albeit shifted to lower k F . Note that this is dissimilar to the case of Σ = 0. The top left panel differs from the top right panel of fig. 9 by the fact that in the limit k F → ∞, ∆ 3 ,κ → 0, as both gaps are proportional to the effective mass m * instead of the bare mass m. This vindicates the simplification we proposed in the expression for the dispersion relation, see also section III B. Building upon the solutions obtained for the gaps and condensate, we can now determine the free energy of the system of fermions and study its phase diagram. We evaluate them in the limit of zero temperature. As the interaction between the fermions is attractive, a gas to liquid transition is expected at large fermion densities [21]. The free energy of the system when superfluid gaps are only a small correction is given by [54] where ω * = m 2 * + p 2 . The pressure P is simply minus the free energy and takes the following form at zero temperature [20] where we have denoted ϕ = Σ/m and C 2 φ = 4 αm 2 χ /(3πm 2 φ ). The adimensional condensate renormalizes the mass as m * = m (1 − ϕ) and C 2 φ is the effective coupling strength. The energy density (not to be confused with the dispersion relation ± ) at zero temperature is then found through where the number density of particles is given by n = k 3 F /3π 2 . Alternatively, the energy density of the fermions is simply given by the integral over all momentum of the Fermi-Dirac distribution weighted by the energy of each mode. The chemical potential is then related to k F by µ 2 = m 2 * + k 2 F where m * is the effective mass of the fermions in the medium. Subsequently, ϕ is determined through the following equation This equation is the same as that of eq. (D27) when the gaps are zero. Furthermore, it is easy to accommodate for repulsive interactions by introducing a vector boson. Such simplified system has been treated extensively in the nuclear physics literature (the so-called σ − ω model) in order to understand the properties of nuclear matter at very high densities. Neutrons, protons interact attractively by exchange of the σ meson and repulsively by exchange of the vector meson ω. Condensation of σ renormalizes the nucleon mass, condensation of ω (more precisely, the 0 th component of ω µ by virtue of rotational symmetry) renormalizes the chemical potential. By fitting the couplings g σ,ω to measured quantities, some characteristics of nuclear matter can be reproduced. Interestingly, this system is similar in spirit to a Lennard-Jones-type potential and some features are expected to be shared between the nuclear matter and more mundane matter. For the scenario that is treated here, i.e. DM, the same parallel can be made, with the attractive interaction rising from the exchange of the scalar mediator and the degenerate Fermi pressure at relativistic Fermi momentum balancing it. Thermodynamic quantities of the system are obtained for a given m and C 2 φ by solving eq. (F4) for ϕ at every k F and then computing P (k F ) and (k F ). The pressure and energy density of a free gas (no interactions) is recovered by taking the limit C 2 φ → 0 which, in turn, sends ϕ → 0 and the in-medium mass becomes the bare mass, m * → m. It is interesting to first consider the binding energy per fermion as a function of density. This broadly paints the behavior of the system. In fig. 12 we show /mn as a function of the adimensional density (k F /m) 3 for a few illustrative values of the effective coupling. For C 2 φ 1 the non-interacting case is recovered. For C 2 φ > 1.09, the binding energy becomes bigger than the rest mass of the fermions and therefore a global minima develops at high density. This signals that the matter in the system can clump into stable bound states of high densities, the so-called nuggets [9,24], which are identified as the liquid phase. To understand better the behaviour of the system at different DM densities and coupling, the isothermal curves of the pressure P as a function of volume v = 1/n are shown in fig. 13; such figures are reminiscent of textbook treatment of the van der Waals equation of state, see for example [62]. The inverse coupling τ = 1/C 2 φ plays a role analogous to the temperature. Indeed, at high τ (small C 2 φ ) the system can only be in the gas phase whereas at small τ (high C 2 φ ) the system can be put in a high density phase. The conjugate variable to τ is identified to be the square of the scalar density ϕ 2 , as seen from the expression of free energy above. For C 2 φ > 0.841, the pressure P (v) start developing a local minima for small volumes, and for C 2 φ > 0.885 the minima becomes global with the pressure dropping below 0. Therefore for strong enough interactions, some parts of the P (v) curves have dP/dv > 0, meaning that the compressibility of the system is negative and that these regions are unstable under any small perturbations. If prepared in such states, the system is expected to collapse to a mixture of low and high density matter at the same pressure.
Constructing the physical EoS is a matter of finding the equilibrium states and of identifying the metastable and unstable configurations. Let us summarize how the physical EoS is built from the standard Maxwell constructions [62]. The Gibbs-Duhem relation dµ = −sdT + vdP provides, at constant temperature (in our case, it would be constant τ ), a relation between the chemical potentials between two states A and B taken to be, say, at low and high density respectively and the area under the curve v (P ) between those two states which, of course, depend on the coupling. For appropriately chosen states A and B, it possible for the chemical potential of, say, B, to be bigger, lower or equal to the chemical potential of A. As can be seen on the right of fig. 13, many different state share the same equal pressure P but different volumes v. The essence of the Maxwell construction is, for an EoS P (v), to identify the two states A and B such that they are at equal pressure and at the same chemical potential. Using eq. (F5), this amounts to looking at two points on the P (v) curve such that the areas between the curve and a constant P line (which is the pressure of both the state A and B) is the same below and above. In fig. 13, the green curve (C 2 φ = 0.9) showcases this construction. In the right panel of fig. 13, we show µ (P ) for the free case, and for C 2 φ = 0.9, 1.09, 1.12, respectively. Following the Maxwell construction of the physical EoS, the physical iso-coupling for C 2 φ = 0.9 (1.09) is explicitly shown and are presented as the solid green (red) line respectively, while the unphysical iso-coupling in shown as dashed lines. As the coupling grows, the gas to liquid phase transition happens at smaller and smaller densities. Interestingly, for C 2 φ > 1.09 the Maxwell construction cannot be made as there are no pairs of states with different densities but with the same µ and P . Technically, for large enough coupling, the area between the P < 0 region of the curve and P = 0 horizontal axis is bigger than the area under the curve of the region where P > 0, even if integrated up to infinite volume. Amusingly, a similar feature is encountered in a fluid analogy to black hole thermodynamics [63].
As mentioned before, the EoS in this model closely resembles that of the Lennard-Jones potential and so a similarity between both phase diagrams is expected. This is mostly conveniently seen in P − τ plane which is shown in fig. 14 (left panel) using the effective coupling C 2 φ = 4 αm 2 χ /(3πm 2 φ ) and the number density n/m 3 . We now interpret the phase diagram presented on the right panel of fig. 14. Given a coupling C 2 φ , a system of DM particles at density n/m 3 and in an equilibrium configuration could behave as a 'gas' (blue regions) or as a 'liquid' (orange regions) or coexist as a mixture of the two (green regions) or in the so-called fluid phase (purple regions), respectively. It is also possible for the system to have no equilibrium configuration strictly speaking but to be instead in metastable (very long lived) or unstable states; the region where such behavior occurs is show in hatched red. In the context of our DM model, say captured at the core of neutron stars, we can assume that such a situation could not be realized.
At this point, it may be instructive to distinguish what we mean by fluid, gas and liquid phases. As is evident from the phase diagram, the gas phase corresponds to low density configuration of the system at moderate C 2 φ ≈ 1. Therefore, the expected behavior is very similar to a non-interacting degenerate Fermi gas. For similar and larger values of the coupling, the liquid phase is characterized by very large densities. Near the saturation density, indicated by a solid orange line, the EoS of the liquid phase closely resembles that of an incompressible fluid [20]. The intermediate coexistence region is reached as we increase the density for similar values of coupling. This transition from the gas phase happens to be of first order [54]. Furthermore, in the so-called fluid phase there are no phase transitions, and consequently the system smoothly goes from the limit of non-interacting Fermi gas to a relativistic one.
Finally, we briefly address the possible impact of the superfluid phase on this phase diagram. We find that the presence of superfluid phase only gives as small correction to the pressure (as the gaps are parametrically smaller than the kinetic energy at the Fermi surface). More precisely, the gaps contributes positively to the total pressure ∝ µ 2 ∆ 2 (k F ). Qualitatively, the correction to picture in fig. 14 is modest.