Yukawa Bound States of a Large Number of Fermions

We consider the bound state problem for a field theory that contains a Dirac fermion $\chi$ that Yukawa couples to a (light) scalar field $\phi$. We are interested in bound states with a large number $N$ of $\chi$ particles. A Fermi gas model is used to numerically determine the dependence of the radius $R$ of these bound states on $N$ and also the dependence of the binding energy on $N$. Since scalar interactions with relativistic $\chi$'s are suppressed two regimes emerge. For modest values of $N$ the state is composed of non-relativistic $\chi$ particles. In this regime as $N$ increases $R$ decreases. Eventually the core region becomes relativistic and the size of the state starts to increase as $N$ increases. As a result, for fixed Yukawa coupling and $\chi$ mass, there is a minimum sized state that occurs roughly at the value of $N$ where the core region first becomes relativistic. We also compute an elastic scattering form factor that can be relevant for direct detection if the dark matter is composed of such $\chi$ particles.


Introduction
One of the mysteries of modern physics is the composition of the dark matter (DM). Various extensions of the standard model (SM) with dark matter candidates have been proposed and studied. A popular scenario has the Higgs boson mediating the connection between DM and the SM. This setup may be testable at LHC experiments which explore Higgs physics. A very simple candidate for the DM is a Dirac fermion (which we denote by χ) that is a singlet under the SM gauge group. If χ interacts with the SM sector through the Higgs boson, the lowest operator has dimension 5 which is not renormalizable. § If the cutoff scale is higher than that of the relevant physical processes, one can use this effective operator description. On the other hand, if the cutoff scale is lower than the DM mass, working with a ultraviolet (UV) completion is necessary. The simplest UV completion has an additional gauge singlet scalar (which we denote by φ) that interacts with the DM through a Yukawa coupling g χ , and interacts with the Higgs boson via renormalizable couplings. For m φ < m χ , the Yukawa coupling allows the DM to efficiently annihilate in the early universe. This very simple dark matter sector and has been previously studied in the literature [1][2][3]. The scenario where DM self interacts with a light mediator has also drawn astrophysical interest [4][5][6][7][8] because of its possible implications for structure formation on small scales.
One interesting aspect of the DM in this model is that for small enough scalar φ mass, the DM particles have a rich spectrum of χ particle bound states. Two-body bound states were studied in [9] and in the non-relativistic regime bound states with N 1 χ particles were studied in [10]. For small values of the dark Yukawa coupling fine structure constant α χ = g 2 χ /(4π) and moderately large values of N α −3/2 χ these bound states are non-relativistic. It was found that the size of these states decreases as N is increased.
As N increases the χ particles become more relativistic and eventually the methods used in [10] are no longer applicable. The nature of the bound states for these larger N 's is not known. The purpose of this paper is to fill in this gap by providing an understanding of the bound states with a large number of fermions where relativistic physics is important. While the main motivation for such states existing in nature is dark matter, the results of this paper are not dependent on that physical interpretation.
One important difference between scalar interactions and vector interactions with fermions is that scalar Yukawa couplings give rise to interactions that are suppressed at large fermion energies. Because of this, we find that the character of the bound states changes as one enters the relativistic regime for the core region. They no longer get smaller as N increases but rather start grow in size. Eventually these bound states become so large that the screening of the Yukawa potential by the scalar mass cannot be ignored. Larger stable bound states with more particles do not exist.
In this paper we use both analytical and numerical methods to study the binding energy and the size of the Yukawa bound states as a function of N . We present our results in a number of § We neglect the renormalizable operatorLHχ which explicitly breaks the symmetry that stabilizes the DM.
Including such interaction with a tiny coupling could lead to the decaying DM scenario and indirect detection signals.
plots and provide a detailed discussion of the methods used. Finally we provide a calculation of an elastic form factor that may be relevant for direct detection experiments if the χ particles are the dark matter.

Fermions with a Yukawa Interaction
In this paper we study bound states of N Dirac fermions χ interacting through the exchange of a light scalar mediator φ. We are interested in understanding the properties of such states at large N , § in particular their binding energy and the dependence of the size of the states R with the number of particles. Exchange of the light mediator results in an attractive potential between the χ particles. The range of the potential is set by the mass of the scalar m φ . Since the Yukawa potential is screened at distances much larger than 1/m φ , stable bound states with radius R so large that Rm φ 1 do not exist. Similarly, for small enough R such that Rm φ 1, m φ can be neglected. Going forward in this paper we will work in the latter regime neglecting m φ . At the quantum field theory level we have a Dirac fermion χ with mass m χ interacting with a scalar φ through the Lagrange density Without loss of generality we take the Yukawa coupling g χ to be positive. The Lagrange density in Eq.
(1) has a shift invariance where g χ φ → g χ φ+c and at the same time the parameter m χ → m χ +c.
Here c is a constant. In this paper we determine the static properties of bound states with a large number of χ particles. So we match this field theory onto a classical theory of N χ particles interacting with a scalar field φ by substituting in the usual particle action the shift invariant mass where x i (t) is the coordinate of the i'th χ particle. Then the shift invariant particle Lagrangian becomes The canonical momenta are, The equation of motion for the scalar field is, (5) § We call such states nuggets.
For solutions that go to zero at spatial infinity Eq. (5) is equivalent to the integral equation, Note that φ appears on the right hand side of the above differential and integral equations implicitly through the dependence of m(x) on it.
The Hamiltonian for this system is where p i is the momentum of the i-th χ particle and the sum of i is over all the N particles. The second term comes from the scalar part of the field theory Lagrangian. We integrated by parts to put this part of the Lagrangian into a form where the equations of motion can be used. Integrating by parts is not consistent with the φ shift symmetry since it assumes the field vanishes at spatial infinity. That is why the shift symmetry is not manifest in the second term in Eq. (7). Schematically for a single massless χ the Hamiltonian in Eq. (7) has the form, At the field theory level each relativistic χ particle interaction with the scalar field φ is suppressed by 1/p. This property of scalar interactions with very energetic fermions causes each factor of the scalar field φ, in the classical particle Hamiltonian above to also come suppressed by 1/p. § Since we are interested in a large number of fermions, N >> 1, throughout this paper we replace sums over particles by an integral over a phase space density f (r, p) of a Fermi gas, i.e., We also assume spherical symmetry. This corresponds to having the χ particles confined to a coordinate sphere of radius R, and at each spatial point having a Fermi sea filled to a Fermi momentum p F (r) that can depend on the radial coordinate r, where the factor of 2 is from the spin degrees of freedom. In this case, the total number of particles Converting the sums over i in Eq. (7) to phase space integrations we find that the total energy of such a Fermi gas with N χ particles and radius R is § This is very different from N particles of charge q and mass m integrating via an electric field, E = −∇φ (em) .

Then the static Hamiltonian is
where the integrations over momenta gave rise to the functions For small z, i(z) ∼ z 3 /3 while for large z, i(z) ∼ z 2 /2. Inside the nugget the scalar field φ(r) satisfies the differential equation, At the origin the scalar field satisfies the boundary condition φ (0) = 0. Outside the nugget Differentiating this gives the boundary condition, With m φ neglected the integral over the angle between the Fermi momentum and the position of the particles in the nugget can be done in Eq. (6) leaving just the radial integral, Differentiating the above equation with respect to r gives for r < R, As was noted earlier the effective mass, m(r) ≡ m χ + g χ φ(r). The field φ(r) is negative and gets larger in magnitude as one goes towards the center of the nugget at r = 0. Suppose at the center m(0) is positive then the above differential equation indicates that m(r) increases as one increases r towards the surface of the nugget. Conversely if m(0) is negative then it decreases as r increases towards the surface of the nugget. Hence the effective mass never changes sign inside the nugget.
In practice for the situations we consider this means that m(r) is always positive. The properties of a nugget can be calculated once the r-dependence of the Fermi momentum p F (r) is known. In general, it can be determined from the hydrostatic equilibrium. We describe this method in Sec. 3 but find it difficult to implement numerically. So in this paper we take a more heuristic approach making an ansatz for the form of p F (r) and then use Eq. (10) to write p F (r) as a function of r, N and R. For each N and R the non-linear differential equation for φ(r) (i.e., Eq. (14)) is solved using the boundary conditions on the radial derivative of φ at the origin and the surface of the nugget. This is used to compute E(N, R) which is then minimized with respect to R at fixed N to find the nugget size and binding energy. In the non-relativistic regime where the equations of hydrostatic equilibrium can be solved this heuristic approach gives reasonable results for the two ansatzs for the Fermi momentum that we will make.
A very simple ansatz the Fermi momentum inside the nugget is a constant independent of r p F (r) = 9πN 4 Constant Fermi momentum in the nugget is not compatible with the physical condition that the Fermi momentum be continuous at the surface of the nugget. A one-parameter family of physically more reasonable ansatzs with this property is 3 Hydrostatic Equilibrium In this section, we briefly discuss the approach to the Yukawa bound state problem based on hydrostatic equilibrium. As discussed in Sec. 2, the key is the knowledge of the Fermi momentum as a function of the position, p F (r). It can be determined by the the energy-momentum conservation law ∂ µ T µν = 0, where T µν is the stress-energy tensor.
In the static situation we are dealing, the above equation becomes ∇ k T kl = 0, and the spatial components are Conservation of stress-energy, ∇ k T kl = 0, then implies the first order integral-differential equation, where p F (r) = dp F (r)/dr. Since ∇ 2 φ(r) is positive we see that p F (r) is negative. Hence the Fermi momentum increases as r decreases.
Converting this integral-differential equation into a second order differential equation one obtains, 1 r 2 d dr as the equation of hydrostatic equilibrium. This is a second order differential equation and the boundary conditions are p F (0) = 0 and p F (R) = 0.
In general, the equation of hydrostatic equilibrium Eq. (22) cannot be solved independentlyit is coupled to the Laplacian equation Eq. (14), and the two differential equations must be solved together in order to determine p F (r) and φ(r). For each radius R, p F (r) fixes the χ number N through Eq. (10), and determines the total energy E through Eq. (11).
In the non-relativistic limit p F (r) m(r), the above equation can be simplified. The number density of χ particles is, n(r) = p F (r) 3 /(3π 2 ) and in this limit the Laplacian of φ is, ∇ 2 φ(r) = g χ n(r). Introducing the effective pressure for non-relativistic Fermi gas, p(r) = p F (r) 5 /(15π 2 m(r)), the equation for hydrostatic equilibrium takes the more recognizable form, Under the further assumption of weak field, gφ m χ , the solution to the equation of hydrostatic equilibrium is known and this solution was discussed in [11]. The application to dark matter bound states was discussed in [10]. We find that in the the non-relativistic weak field regime, the choice a = 1/2 in Eq. (19) is a very good approximation and as was remarked on in [10] even the constant Fermi-momentum ansatz in Eq. (18) gives results that have the correct qualitative behavior.

Perturbation Theory
The main purpose of this section is to show how perturbation theory in the coupling breaks down for large N . The problem with perturbation theory does not depend on the explicit ansatz for the r dependence of the Fermi momentum and so in this section we use the very simple constant Fermi momentum ansatz, p F (r) = p F θ(R − r). Expanding in powers of g χ where the term with subscript n is of order g 2n+1 χ .
Neglecting m φ we have from eq.(6) that, Using our ansatz for p F (r) this becomes, where the function i(z) was defined in the previous section. At the next order, which with our ansatz for p F (r) becomes, where the function  With the perturbative solution for φ the total energy E(N, R) can be obtained using Eq. (11). Minimizing it with respect to R (at fixed N ) yields the radius of the ground state of N χ particles and the binding energy of that state. We find as one increases N the radius first shrinks and then expands as shown in Fig. 1 where we used m χ = 100GeV and α χ equal to 0.1 and 0.01.
The hierarchy φ 0 φ 1 is satisfied throughout the nugget. Despite this perturbation theory breaks down for large N because the field becomes very large towards the core of the nugget driving m(r) from positive to negative values (the dashed curves in Fig. 2) which is in conflict with the general results we derived earlier.

Numerical Approach
Since solving the hydrostatic equilibrium equations is too difficult, and perturbation theory is not valid for large N because the field φ gets too large, we adopt the method described in Sec. 2. For fixed N , the Laplacian equation (14) is solved numerically, for different choices of R. We present results for both the constant and power law (with a = 1/2) ansatzs for p F (r) given in Eqs. (18) and (19) respectively. The energy in Eq. (11) is minimized as a function of R to determine the radius and binding energy of a nugget containing N χ particles.
The model parameters that determine the physics of nuggets are, m χ , α χ = g 2 χ /(4π) and the mediator mass m φ . Our ansatzs for the dependence of the Fermi momentum on the radius do not introduce any new dimensionful parameters once they are normalized to the number of particles. We are neglecting m φ here and so the only dimensionful parameter is m χ . Hence we work at the fixed value, m χ = 100GeV. Using dimensional analysis we can determine the dependence of physical quantities on m χ , for example R ∝ 1/m χ . To capture the trends with α χ we display our results for two values 0.1 and 0.01.
In Fig. 3 and 4, we show m(r) and p F (r) throughout the nugget, for different values of N . We find that for small N the χ particles are non-relativistic throughout the bound state. For larger N , the χ particles are not ultra relativistic near the surface, but are ultra relativistic in the core region. The effective mass m(r) becomes small near the core but does not change sign.
In Fig. 5, we plot the nugget radius R and binding energy per particle, as a function of the number of particles N , for the same sets of parameters as Fig. 3. Note that as N increases, the radius R first shrinks and then grows. At larger N , more of the χ particles are relativistic, the Yukawa interactions among these particles are m/p suppressed, and the Fermi pressure pushes the minimal energy configuration to larger R. The binding energy per particle ε increases monotonically with N , indicating that the nuggets are stable against breaking up into smaller pieces. The behavior of R and ε as functions of N are the main results of this paper. The qualitative behavior of our results for the N dependence of the nugget radius and binding energy do not depend on the the particular ansatz for the Fermi momentum chosen. The most striking feature, that the radius first decreases with N and then increases, occurs for both ansatzs for p F (r) that we used and is even present if perturbation in the coupling α χ is used to determine φ(r).
Throughout this paper, we have neglected self-interactions of the φ field. We can estimate the range of self-couplings λ for which this is a reasonable approximation. Suppose the Hamiltonian in Eq. (7) contains an additional term δH[φ] = d 3 xλφ(x) 4 /4!. Using our solution for φ, we show in Fig. 6 (for α χ = 0.1 and 0.01) the largest values of λ for which this new contribution does not to exceed the total energy we calculated before, Here the a = 1/2 power law ansatz for p F (r) in Eq. (19) was used. Clearly for larger N this requires smaller λ for the self interaction contribution to be negligible. A non-zero λ term draws the φ field closer to 0, thus it tends to reduce the nugget radius R.   To close this section, we discuss the possible existence of nuggets with radius R much larger than the screening length 1/m φ . It seems intuitively clear that the screening forbids such states but it is not difficult to address this issue somewhat quantitatively using the methods developed in this paper. Assuming the constant Fermi momentum ansatz p F ∼ N 1/3 /R, that the χ particles are non-relativistic and that the weak field limit is appropriate, the sum of the total kinetic and potential energy for a nugget with R 1/m φ is Here a and b are positive constants. The first term in Eq. (32) arises from the kinetic energy and the second from the potential energy. For fixed N , the above expression is minimized at R = 0 or ∞ and hence stable nuggets with R 1/m φ do not exist.

Application to Dark Matter
If the χ particles are the DM, then for a range of parameters nuggets can form in the early universe [10], shortly after the (free) χ's freeze out. Such dark nuggets can be cosmologically abundant if the DM relic density is dominated by χ particles, i.e., the DM is asymmetric. In the example shown in Fig. 5, it is possible to have a bound state containing thousands of weak scale χ particlesa nugget with mass of order ∼ 10 3 TeV. Usually such supermassive DM cannot be fully thermalized in the early universe due to its small annihilation cross section. However, our analysis of Yukawa bound states provides a possible way to have a supermassive and thermal DM candidate.    The rate for their direct detection is determined partly by the elastic form factor F (q 2 ), r 2 dr sin(qr) qr (m(r)) 3 i(p F (r)/m(r)) .
For q = 0, the naive expectation for the form factor is F (0) = N . However relativistic particle couplings to the scalar mediator are suppressed and this causes F (0) to be less than N . This feature becomes more prominent for large N since then more of the χ particles are relativistic. The form factor falls as the momentum transfer q 2 increases because the scattering becomes less coherent. These features are shown in Fig. 7. This form factor will also be relevant for elastic DM-DM scattering through virtual φ particle exchange.
In this paper we have focussed on the lowest energy bound state of N χ particles. Of course there are excitations and so inelastic processes are possible. We have not explored the form factors relevant in that case but expect that the inelastic channels are suppressed at small q because of Pauli blocking.
Indirect detection signals of DM bound states are also very important. Since the binding energy per χ particle grows with the number N , the χ's are more deeply bounded in larger nuggets. There will be a release of energy by emission of a mediator φ (either real or virtual) when a free χ particle is captured by a nugget or when a small nugget captured by a large one. In the model discussed in Ref. [10], a real mediator φ materializes as a SM µ + µ − or ππ final state. This could offer interesting signals for indirect DM searches using cosmic rays.

Outlook
We have studied bound states of a large number of Dirac fermions χ interacting through exchange of a light scalar field that they are Yukawa coupled to. For very large N the cores of these objects contain very relativistic χ's and the size of the state R increases with N . That is in contrast to the smaller N regime where the χ particles are non relativistic and the size of the state shrinks as N increases. There are several extensions of this work that are worth pursuing.
We made a number of approximations in order to draw these conclusions that are worthy of further exploration. For example we used a simple Fermi gas model for the fermions. For strong enough coupling pairing of fermions and Bose-Einstein condensation may occur. Also we neglected the scalar self coupling. A preliminary estimate we made shows that in some cases the scalar self coupling must be very small for this to be a good approximation and so it would be nice to extend our analysis to include the scalar self coupling. We have shown that nuggets with Rm φ 1 are not stable. It would be interesting to explore in more detail the properties of nuggets in the regime Rm φ ∼ 1.
Finally for the interpretation of χ particles as dark matter it is important to estimate the fraction of dark matter that ends up in bound states with N > 2, § i.e., the analog of big-bang nucleosynthesis § The cosmological production of two body bound states was considered in [10].