Neutrino bound states and bound systems

Yukawa interactions of neutrinos with a new light scalar boson $\phi$ can lead to formation of stable bound states and bound systems of many neutrinos ($\nu$-clusters). For allowed values of the coupling $y$ and the scalar mass $m_\phi$, the bound state of two neutrinos would have the size larger than $10^{12}$ cm. Bound states with sub-cm sizes are possible for keV scale sterile neutrinos with coupling $y>10^{-4}$. For $\nu$-clusters we study in detail the properties of final stable configurations. If there is an efficient cooling mechanism, these configurations are in the state of degenerate Fermi gas. We formulate and solve equations of the density distributions in $\nu$-clusters. In the non-relativistic case, they are reduced to the Lane-Emden equation. We find that (i) stable configurations exist for any number of neutrinos, $N$; (ii) there is a maximal central density $\sim 10^9$ cm$^{-3}$ determined by the neutrino mass; (iii) for a given $m_\phi$ there is a minimal value of $Ny^3$ for which stable configurations can be formed; (iv) for a given strength of interaction, $S_\phi = (ym_\nu/m_\phi)^2$, the minimal radius of $\nu$-clusters exists. We discuss the formation of $\nu$-clusters from relic neutrino background in the process of expansion and cooling of the Universe. One possibility realized for $S_\phi>700$ is the development of instabilities in the $\nu$-background at $T<m_\nu$ which leads to its fragmentation. For allowed values of $y$, cooling of $\nu$-clusters due to $\phi$-bremsstrahlung and neutrino annihilation is negligible. The sizes of $\nu$-clusters may range from $\sim$ km to $\sim 5$ Mpc. Formation of clusters affects perspectives of detection of relic neutrinos.

In recent years neutrinos and new physics at low energy scales became one of the most active areas of research. In particular, neutrino interactions with new and very light bosons have been actively explored, with considerations covering e.g. oscillation phenomena [1][2][3][4][5][6] and various astrophysical and cosmological consequences [7][8][9][10][11][12]. Yet another possible consequence of such interactions (the subject of this paper) is the formation of bound states and bound systems of neutrinos.
The exploration of possible existence of N -body bound states of neutrinos (neutrino clouds, balls, stars, or halos) has a long story. In 1964 Markov proposed the existence of "neutrino superstars" formed by massive neutrinos due to the usual gravity [13]. The system is similar to neutron stars and can be described as the degenerate Fermi gas in thermodynamical and mechanical equilibrium. Substituting the mass of neutron, m n , by the mass of neutrino, m n → m ν in results for neutron stars, Markov obtained the radius of a neutrino star: which corresponds to the Oppenheimer-Volkoff limit. In comparison to the radius of neutron star, R increases by a factor of β 2 where β ≡ m n m ν = 1.9 · 10 10 , and for numerical estimation we take m ν = 0.05 eV. Correspondingly, the mass of neutrino star increases and the central density decreases as For m ν ≈ m e used by Markov, the size, the mass and the number density of the superstars would be R = 10 12 cm, M = 10 6 M and n ν 10 29 cm −3 correspondingly. For the present values of neutrino mass (0.1 eV), the relations (1) and (2) give R 5 · 10 26 cm, M ν = 4 · 10 20 M , n ν (0) 10 8 cm −3 .
The radius is only 20 times smaller than the size of observable Universe (the present Hubble scale), and the total number of neutrinos is 10 times larger than the number of relic neutrinos in whole the Universe with the standard density. The bound states of much smaller sizes and higher densities with substantial implications for cosmology and structure formation as well as for laboratory experiments require new physics beyond the Standard Model.
Another way to reduce the size and mass of a neutrino star is to consider new long-range interactions between neutrinos which are much stronger than gravity. In the case of Yukawa forces due to a new light scalar boson φ with the coupling constant to neutrinos, y, the Newton coupling G N should be substituted in (1) by As a result, we obtain Notice that this expression does not depend on the mass of mediator m φ as long as m φ 1/R. For y = 10 −7 and m ν = 0.05 eV, one finds from (3) √ G ν = 0.56 MeV −1 (that is, 22 orders of magnitude larger than √ G N ), and correspondingly Eq. (4) gives R = 0.7 km.
Properties and formation of the ν-clusters due to new scalar interactions were studied in [17,18] using equations of Quantum Hadrodynamics. The central notion is the effective neutrino massm ν in the background of degenerate neutrino gas and scalar field. The equation form ν was obtained in the limit of static and uniform medium. This non-linear equation depends on the Fermi momentum k F and the strength of interactions which is defined as Here ω is the number of degrees of freedom. For large densities before neutrino clustering, the effective mass is close to zero [17]:m At small densities:m ν → m ν . Usingm ν , the energy density of neutrinos, ρ ν , and the energy density of the scalar field, ρ φ , were computed in [17] which gives the total energy per neutrino here n F = k 3 F /6π 2 . For large enough strength G φ the energy tot as a function of k F acquires a minimum with tot < m ν around k min F ≈ 0.8m ν . Clustering of neutrinos starts when k F decreases due to the expansion of the Universe down to It is also assumed that with the expansion (decrease of n), the kinetic energy could be converted into the increasing effective massm ν and the gas became strongly degenerate.
To describe finite size objects, the static equations of motion were used with non-zero spatial derivatives of the fields, and consequently,m ν . It was assumed that the distribution of clusters on the number of neutrinos N follows the density fluctuations and has the Harrison-Zeldovich spectrum. In [17], the neutrino mass m ν = 13 eV motivated by the Tritium experiment anomaly was used. Very small mediator masses m φ ∼ 10 −17 eV and couplings g ∼ 10 −14 were taken. Consequently, the sizes of clusters equal ∼ 10 13 cm, which is roughly the size of the Earth orbit, and the central densities could reach 10 15 cm −3 .
Later the bound systems of heavy Dirac fermions ("nuggets") due to Yukawa couplings with light or massless scalar were explored [19][20][21][22]. With fermion masses m D ∼ 100 GeV and coupling α φ ∼ (0.01−0.1) such systems were applied to the Asymmetric Dark Matter (ADM). The description of nuggets in [20], [21] is similar to that of neutrino stars [17] 1 . The effective mass of fermion was introduced as in [17] and equations for its dependence on r was derived. The system of equations for the scalar field and Fermi momentum of ADM particles p F (r) (and therefore their density) was presented. In [20] the equations were solved by introducing an ansatz for p F (r), thus reducing the system to a single equation. In [21] the complete system of equations was solved numerically. Dependence of properties of nuggets on N , in particular R(N ), were studied. It was established in [20] that with the increase of N the radius first decreases, reaches minimum and then increases. The increase is associated to transition from non-relativistic to relativistic regime. While in [20] the mediator was assumed to be massless, in [21] dependence of characteristics of nuggets on m φ was studied and, in particular, the case of "saturation", which corresponds to R > 1/m φ , was noticed.
Due to strongly different values of masses and coupling constant the formation and cosmological consequences of nuggets are completely different from those of neutrino stars. Also observational methods and perspectives of detection are different. Our paper is continuation of the work in [17] for neutrinos. Where possible, we compare our results with those for ADM.
We perform a systematic study of various aspects of physics of neutrino clusters. The paper contains a number of derivations and auxiliary material which will be useful for further studies and applications.
In this paper, we revisit the possibility of formation of the neutrino bound states and bound systems due to neutrino interactions with very light scalar bosons. First the bound states of two neutrinos are considered. We formulate conditions for their existence and find their parameters such as eigenstates and eigenvalues. For usual neutrinos and allowed values of y and m φ , the sizes of these bound states are of the astronomical distances which has no sense. For sterile neutrinos with mass m ν ≥ 10 keV and y > 10 −4 the size can be in sub-cm range and formation of such states may have implications for dark matter.
Then we study properties and formation of the N -neutrino bound systems -the ν-clusters with a sufficiently large N , so that the system can be quantitatively described using Fermi-Dirac statistics. We compute characteristics of possible stable final configurations of the ν-clusters as function of N -mostly in the state of degenerate Fermi gas (ground state) but also with thermal distributions. Energy loss of the ν-clusters is estimated. Depending on N (we explore the entire possible range) the clusters can be in non-relativistic or relativistic regimes and we study the transition between these two regimes. In the non-relativistic case, using equations of motions for ν and φ, we rederive the Lane-Emden equation for the neutrino density. We obtain qualitative analytic results for parameters of the neutrino cluster. The equation for the relativistic case is derived which is similar to the one used in [17]. We use presently known values of neutrino masses which are two orders of magnitude smaller than in [17], and consequently, parameters of the clouds are different.
Understanding the formation of neutrino clusters from homogeneous relativistic neutrino gas requires numerical simulations. In this connection we present some relevant analytic results. One possible mechanism of formation is the development of instabilities in the uniform neutrino background, which leads to fragmentation and further redistribution of neutrino density, thus approaching final degenerate configurations. In some parameter space, one can use an analogy between Yukawa forces and gravity and therefore compare it with the formation of dark matter halos due to gravity.
The paper is organized as follows. We start with re-derivation of the Yukawa potential in Sec. II, where we adopt the approach developed in Ref. [5] and generalize it to the relativistic case. In Sec. III, we study the two-body bound states in the framework of quantum mechanics. In Sec. IV, the N -body bound systems are considered. The results of Sec. III and Sec. IV apply to the non-relativistic neutrinos, while a study of the relativistic case is presented in Sec. V. In Sec. VI, we speculate on possible mechanisms of formation of the ν-clusters which include fragmentation of the cosmic neutrino background induced by the expansion of the Universe. Further evolution and formation of final configurations are outlined. Energy loss of the clusters is estimated. Analogy with formation of DM halos is outlined. Finally, we conclude in Sec. VII. Various details and explanations are presented in the appendices.

II. YUKAWA POTENTIAL FOR NEUTRINOS
In general, to form a bound state, the range of the force, R force , needs to be longer than the de Broglie wavelength of the neutrino, λ ν : For a non-relativistic neutrino λ ν = (m ν v) −1 , where v is the neutrino velocity, and R force ∼ 1/n φ is given by the inverse of the mediator mass. Therefore Eq. (5) can be written as m φ < m ν v, which implies that mediator should be lighter than the neutrino to form a bound state 2 . Let us consider a light real scalar boson φ interacting with neutrinos We assume m φ m ν , so that ν can be treated as point-like particles while φ as a classical field. This is valid as long as we are not concerned with hard scattering processes. In this regime, the effect of the φ field on ν particles can be described using Yukawa potential. From Eq. (6), we obtain the Euler-Lagrange equations of motion (EOM) for ν and φ: According to Eq. (7) the effect of the φ field on the motion of ν can be accounted by shifting its mass: Here V is the potential, andm ν is treated as the effective mass of neutrino in background. At quantum level νν is replaced by the expectation value computed as (see Appendix A): where E p = m 2 ν + |p| 2 and f (t, x, p) is the neutrino distribution function normalized in such a way that the number density, n, and the total number of ν, N , are given by In the non-relativistic case,m ν /E p ≈ 1, Eq. (10) gives We assume a static distribution, so that ∂ t φ = 0. Replacing ψψ by ψψ ≈ n in Eq. (8) and using relation φ = V /y, (see Eq. (9)) we obtain equation for the potential For spherically symmetric distribution of neutrinos, n = n(r), where r is the distance from the center, Eq. (12) can be solved analytically [5], giving For n = N δ 3 (x) (all neutrinos are in origin) the solution reproduces the Yukawa potential Note that V is the potential experienced by a given neutrino, and for two neutrinos separated by a distance r, the total binding energy is given by Eq. (14) with N = 1 (see Appendix B for details).
For relativistic neutrinos, if ∂ t φ can be neglected 3 , one should use general expression in (10) and the equation for the potential becomes In the non-relativistic caseñ → n.
In this paper we focus on the real scalar field. Pseudo-scalar mediators (e.g. the majoron) cause spindependent forces, and the corresponding potential can be found in [23]. For a two-body system, the potential might be able to bind the pair if the mediator mass is sufficiently light. In an N -body system with large N , the effect of spin-dependent forces is suppressed by spin averaging.
In the case of a complex scalar, our analysis can be applied to its real part, while the effect of its imaginary part, which is a pseudo-scalar, can be neglected if the two parts have the same mass. In models with spontaneous symmetry breaking such as the one in [23], only the pseudo-scalar component is massless and therefore leads to long-range effects.

III. TWO-NEUTRINO BOUND STATES
In this section, we consider bound states of two non-relativistic neutrinos due to the attractive φ potential. For DM particles this was considered in a number of papers before (see [19] and references therein).
The strength of the interaction required to form bound states can be estimated in the following way. If the wave function of the bound state is localized in the region of radius R, then according to the Heisenberg uncertainty principle the momentum of neutrino is p ∼ 1/R, and consequently, the kinetic energy, E kin ≈ p 2 /2m ν , equals The condition of neutrino trapping in the potential V reads Plugging E kin from Eq. (16) and V from Eq. (14) with N = 1 into this equation we obtain Maximal value of the r.h.s. of this equation equals y 2 /(4πem φ ), which is achieved at R = m −1 φ . Then, according to (18), the sufficient condition for existence of bound state can be written as where λ can be interpreted as strength of interaction. Quantitative study below gives λ > 0.84. The condition can be rewritten as the upper bound on m φ :

eV.
A. Quantitative study Quantitative treatment of neutrino bound states requires solution of the two-body Schrödinger equation. Let r 1 and r 2 be the coordinates of two neutrinos. The potential depends on the relative distance between neutrinos r ≡ r 1 − r 2 only. Consequently, the wave function of bound state depends on r, ν = ν(r), and it obeys the Schrödinger equation with V (r) = V (r, N = 1) given in Eq. (14), r = |r| and E being the energy eigenvalue. The numeric coefficient of the kinetic term takes into account the reduced mass of neutrino (see Appendix C). Since the potential V (r) is independent of the direction, the wave function can be factorized as   (19). In all the panels we take the orbital momentum l = 0. For λ = 0.5, since no bound states can be formed, we plot the wave function with the lowest energy.
Here u(r) is a radial component and Y m l (θ, ϕ) is a spherical harmonic function. Plugging Eq. (21) into Eq. (20), we obtain the equation for u(r): For zero angular momentum, l = 0, the equation simplifies further: We solve this equation numerically (see details in Appendix C). u(r) for different values of λ are shown in Fig. 1. According to this figure there is no bound state for λ = 0.5: the solution with the lowest energy level does not converge to zero when r → ∞. For other specific values, λ = 0.9, 3.6, and 9.0, the figure shows one, two, and three bound states. With further increase of λ, the number of allowed energy levels of increases. Varying λ we find that for • λ < 0.84, no bound states can be formed; • 0.84 < λ < 3.2, there is one bound state (the 1s state); • 3.2 < λ < 7.2, two bound states with l = 0 (the 2s state) can be formed; • λ > 7.2, three or more bound states exist, including states with l = 0.

B. Coulomb limit
The limit of m φ → 0 (λ → ∞) corresponds to a Coulomb-like potential, and in this case, the Schrödinger equation (22) can be solved analytically. There is an infinite number of levels and the wave function of ground state equals where is the radius of the bound state, and the corresponding eigenvalue equals This E can be interpreted as the binding energy of the state. For non-zero m φ such that R 0 m −1 φ (equivalent to λ 1), the range of the force is much larger than the radius of the bound state in the Coulomb limit. In this case, Eq. (24) and (25) can be used as a good approximation. The correction to the Coulomb limit result due to non-zero m φ can be computed using the variation principle (see below).
Assuming solution in the form Eq. (23): u(r) = 2a 3/2 re −r/a with a being a free parameter, we find Minimization of H with respect to a, that is d H /da = 0, gives (16πλ + am ν y 2 ) 3 − 32πλ 2 am ν y 2 (16πλ + 3am ν y 2 ) = 0, which is a cubic equation with respect to a. The cubic equation has one negative solution which can be ignored, and two others which can be either both complex (in this case the two solutions are conjugate to each other) or both real and positive (in this case the two solutions are identical). Only the later case leads to a physical solution and requires that the discriminant is positive: From this inequality we obtain which is about 13% larger than the exact value 0.84 obtained in Sec. III A. In [19] the bound λ ≥ 0.84 was also obtained. Using the variation principle in this way may not give an accurate wave function, but usually leads to a very accurate result for eigenvalues [24]. Note that here we do not require that λ −1 ∝ m φ is perturbatively small. For very small λ −1 the result is expected to be more accurate. The series expansion of the positive physical solution of Eq. (27) in λ −1 gives Therefore, the radius of the system R = a = x/y 2 m ν equals The eigenvalue is corrected according to (26) as The lowest order expressions for R and binding energy (28) coincide with those in [19]. Numerically, the radius equals This is much larger than the distance between relic neutrinos (without clustering) ∼ 0.16 cm, and therefore existence of bound states of such size has no practical sense. For two-body bound states, it is impossible to enter the relativistic regime. Indeed, according to Eq. (25) or Eq. (28) the kinetic energy being of the order of binding energy equals y 4 m ν /64π 2 . Therefore the relativistic case, E kin /m ν ≥ 1 implies y (64π 2 ) 1/4 which would not only violate the perturbativity bound but also be excluded by various experimental limits -see Appendix D.
The results of this section may have applications for sterile neutrinos with mass m ν ≥ 10 keV. These neutrinos could compose dark matter of the Universe. For m ν ≥ 10 keV and y = 10 −3 we find from (29) R = 0.05 cm.

IV. N -NEUTRINO BOUND SYSTEMS
Let us consider N neutrinos trapped in the Yukawa potential with N being large enough, so that neutrinos form a statistical system which is described by the Fermi-Dirac distribution function Here E is the total neutrino energy, µ is the chemical potential and T is the temperature. Recall that for this distribution the number density, n, the energy density, ρ, and pressure P of the system equal respectively The ground state of this system (T → 0) is the degenerate Fermi gas with distribution f (p) = 1, if E < µ, and f (p) = 0 if E > µ according to Eq. (30). (The case of finite T will be discussed in Sec. VI). Consequently, in the degenerate state and the total number of neutrinos N for uniform distribution in the sphere of radius R equals Eq. (31) applies to both relativistic and non-relativistic cases. For ρ and P, the non-relativistic results are while the relativistic results can be obtained by computing the integrals numerically.

A. Qualitative estimation
Let us consider N neutrinos uniformly distributed in a sphere of radius R. The Yukawa attraction produces a pressure, P Yuk , which for m φ = 0 can be estimated as where V = 4 3 πR 3 is the volume and the quantity in the brackets is the total potential energy. Static configuration corresponds to the equilibrium between this Yukawa pressure and the pressure of degenerate fermion gas (Pauli repulsion): Using Eq. (33) and Eq. (34) this equality can be rewritten as Inserting Eq. (32) in Eq. (35) we find Alternatively, we can express p F in terms of N and obtain Plugging R from Eq. (37) back into Eq. (32), we obtain a relation between N and p F : Notice that for N → 2, Eq. (37) recovers approximately the result in Eq. (24) for the two-body system. A noteworthy feature of Eq. (37) is that the radius decreases with N : R ∝ N −1/3 . This implies that increasing N leads to a more compact cluster. This dependence is valid only in the non-relativistic regime. When N is too large, the Fermi momentum p F which increases with N according to Eq. (38) will exceed m ν and the system enters the relativistic regime. As we will show in Sec. V, in the relativistic regime, the attractive force becomes suppressed, causing R to increases with N .

B. Quantitative study
In an object of finite size, the density and pressure depend on distance from the center of object. The dependence can be obtained using the differential (local) form of the equality of forces, which is given by the equation of hydrostatic equilibrium. For degenerate neutrino gas it is reduced to the Lane-Emden (L-E) equation [25]. The L-E equation can be also derived from equations of motion of the fields (see Appendix F).
The equation of hydrostatic equilibrium expresses the equality of the repulsive force due to the degenerate pressure, F deg , and the Yukawa attractive force, F Yuk , acting on a unit of volume at distance r from the center: For non-relativistic neutrinos, P deg is given in Eq. (33), and using Eq. (31), P deg can be rewritten as The Yukawa force equals where is the number of neutrinos inside a sphere of radius r. Inserting Eqs. (40) and (41) into Eq. (39) gives the equation for n(r): Dividing this equation by n/r 2 and then differentiating by r, we obtain where Eq. (42) is nothing but the Lane-Emden equation [25] for specific values of parameters. In particular, γ, defined as the ratio of powers of n in the r.h.s. and the l.h.s. of Eq. (42), equals γ = 1/(2/3) = 3/2. For this value the N -body system has a finite size R defined by n(R) = 0.
Integrating n(r) with r ranging from 0 to R gives the total number of neutrinos: Combining Eqs. (43) and (44), we obtain expression for the radius which shows that R decreases with N , in agreement with the qualitative estimations in Sec. IV A. The results in Eqs. (43), (44) and (45) have the same parametric dependence as Eqs. (36), (38) and (37) and differ by numerical factors 2.5, 2.3 and 3 correspondingly. Notice that for the Fermi momentum p F 0 ∼ m ν the Eq. (43) gives R = 19.9/ym ν , which practically coincides with (4) obtained from rescaling of parameters of neutron star. The coincidence is not accidental since parameters of neutron star correspond to p F ∼ m n . Other characteristics of a ν-cluster for p F 0 ∼ m ν are M ν = 3 · 10 −11 g (for y = 10 −7 ) and n 0 = 2.6 · 10 8 cm −3 . The central density is determined according to Eq. (31) as In Eq. (42) both y 2 and κ can be absorbed by the redefinition of r: In terms of r , Eq. (42) can be written as the universal form and its solution depends on the boundary condition n(0) which is free parameter restricted from above by p F 0 < m ν . Thus, solution for different values of y and m ν can be obtained from rescaling of r: where r 7 ≡ r(y = 10 −7 ). The radius of a cluster is rescaled as in (47). For y = 10 −14 this equation gives R ∼ 10 12 cm, i.e. orders of magnitude smaller than the Earth orbit. For the ν-cluster made of ν 2 (in the case of mass hierarchy) the radius is 6 times larger. The mass of such a cluster equals 3 · 10 10 g. The central density, and consequently, the mass of ν-cluster is an independent parameter restricted by the neutrino mass. This consideration matches results in Eqs. (43) and (44).

C. Non-degenerate ν-clusters
As we will show in Sec. I the energy loss (cooling) of ν-clusters due to φ bremsstrahlung is negligible. Therefore the final configuration of degenerate neutrino gas may not be achieved. In this connection we consider the non-degenerate neutrino system characterized by a non-zero temperature T . The number density equals and I 3 is related to the Riemann zeta function: For simplicity we assume that neutrinos are distributed uniformly in the sphere of radius R T . Then If neutrinos are non-relativistic, the total kinetic energy can be computed as where I 5 = 45ζ(5)/2 = 23.33 [see (48)] and we used Eq. (49) for R T . The potential energy of a cluster equals or eliminating R T : The ratio of the kinetic and potential energies can be written as The ratio ξ E in Eq. (53) depends on the temperature T . Therefore, given a certain value of ξ E , it allows to determine T : For T → p F /3 it matches Eq. (44), numerically: Then, according to (49) the radius equals and numerically These equations match Eq. (45). The stable configuration of a cluster corresponds to the hydrostatic equilibrium P T = −P Yuk . In the nonrelativistic case, P T = (2/3)ρ T and therefore from Eq. (50) we have E K = 3/2V P T . So, P T = (2/3)E K /V . The Yukawa pressure equals P Yuk = E V /V . From this we find that the hydrostatic equilibrium corresponds The ratio of the degenerate ν-cluster radius in Eq. (45), denoted by R deg , to the radius R T , equals where the last equality is for the hydrostatic equilibrium. The ratio does not depend on other parameters. So, if the evolution proceeds via the formation of a cluster in hydrostatic equilibrium, achieving final configuration require that the radius decreases by a factor of two. Our consideration of the thermal case is oversimplified: applying the Hydrostatic equilibrium locally leads to a distribution of neutrinos in a cluster with a larger radius. Therefore, the density of a non-degenerate cluster equals and numerically n = 6.2 · 10 6 cm −3 y 10 −7

A. Equations for degenerate neutrino cluster
For a sufficiently large N , neutrinos become relativistic (p F m ν ) in the center of cluster, while near the surface they are almost at rest (p F = 0). For relativistic neutrinos the Yukawa potential is suppressed by the factorm ν /E p in Eq. (15). Due to this suppression, the pressure of degenerate gas is able to equilibrate the Yukawa attraction for arbitrarily large N , in contrast to the gravity case.
The system is described by Eq. (15) and relativistic generalization of equilibrium condition between the Yukawa attraction and the degenerate neutrino gas repulsion, Eq. (39). To avoid potential ambiguities in the definition of forces in the relativistic regime, instead of F Yuk and F deg we will use the derivatives dφ/dr and dp F /dr which represent these two forces and are well-defined quantities in the relativistic case. The equilibrium condition can be derived from the equations of motion for the fields (see Appendix E) giving (r < R) which is the relativistic analogy of equality F Yuk = F deg . In turn, Eq. (15) for the scalar field can be rewritten as is the effective neutrino density andm ν is introduced in Eq. (9). Let us underline that since the effective neutrino mass is given bym ν , the relativistic regime is determined by the condition p F /m ν 1 which can differ substantially from p F /m ν 1.
Using Eq. (59), we can rewrite the system of Eqs. (57) and (58) in terms ofm ν and p F : Eqs. (60) and (61) is a system of differential equations for p F andm ν with the boundary conditions in the center:m Eq. (61) can be integrated from r = 0 to a given r: Insertingm ν from this equation into Eq. (60) we obtain a single equation for p F , which can be solved numerically 4 and using n(r) = p 3 F (r)/6π 2 (which still holds for relativistic neutrinos) to compute the number density.
In the differential equations, the evolution ofm ν represents the evolution of the scalar field while p F determines the neutrino density. Therefore the boundary condition form ν is essentially the condition for the field φ. Note that p F 0 /m 0 is the measure of relativistic character of the system. Hence p F 0 /m 0 is used as the main input condition which determines the rest, whilem 0 is determined by a self-consistent condition below.
As in the case of the non-relativistic Lane-Emden equation, at certain r = R, the Fermi momentum p F (and hence the density) drops down to zero: and this defines the radius of a cluster. Since the field φ extends beyond a star radius, r > R, the effective massm ν differs from m ν at r > R. But it should coincide with the vacuum mass m ν at r → ∞. To obtain a self-consistent solution, we need to tunem 0 to matchm ν (r → ∞) with m ν , i.e.m 0 is not a free parameter but a parameter determined by the neutrino mass at r → ∞.
The system of Eq. (60) and (61) is similar to that for nuggets in [21]. Indeed, differentiating Eq. (5) of [21] with respect to r gives exactly Eq. (61). Eq. (60) corresponds (up to a factor of 1/π) to Eq. (6) in [21]. However, the boundary conditions (7) there differs from our boundary conditions determined by Eq. (62) and by the aforementioned matching ofm ν (r → ∞). The boundary condition in the relativistic regime together with the matching is non-trivial as already noticed in [17]. This condition can affect the existence of solutions.
Eqs. (60) and (61) can be reduced to the Lane-Emden equation in the non-relativistic limit. Indeed, in this limitñ ≈ n and for m φ = 0 we obtain from (60) which is the spherically symmetric case of (60). From (57), neglecting V = yφ in comparison to m ν we find Inserting this expression for the derivative into Eq. (63) gives Eq. (42).

B. Solution for massless mediator
Following the above setup, we first solve the differential equations (60) and (61) with m φ = 0. In Fig. 3, we show the obtained density n(r) and the effective densityñ(r) profiles for several values of p F 0 /m 0 . Recall that the effective density is the number density that includes the suppression factorm ν /E p . This factor depends on p F , and decreases toward the center. As a result, suppression in the center is stronger and the maximum ofñ(r) occurs at certain point away from r = 0. This effect becomes stronger with the increase of Using the obtained density profiles we computed various characteristics of the neutrino clusters for y = 10 −7 and m φ = 0 (see Tab Using the numbers in the second and third lines we find p F 0 /m ν (the 4th line), which in turn, gives the central density of a ν-cluster: As N increases, the central density n 0 (also p F 0 ) first increases, reaches maximum at N y 3 = 1.5 · 10 2 or N = 1.5 · 10 23 and then it decreases. The maximum corresponds to p F 0 0.6 m ν , that is, to the point of transition between the non-relativistic and relativistic regimes. The maximal density is determined by the mass of neutrino only: The corresponding values of radius of a cluster are shown in the last line of Tab. I. On the contrary, as N increases, the radius R first (in non-relativistic range) decreases, reaches the minimum at N ≈ 1.5 · 10 23 and then increases, as shown in Fig. 4. The minimal radius (64) is about 20% smaller than the non-relativistic result in Eq. (43) for p F = 0.3m ν , which can still be accurately computed using the non-relativistic approximation according to the difference between the dashed and solid curves in Fig. 3. Notice that dependence n(r) substantially differs from n(r) = p 3 F (r)/6π 2 with p F (r) taken from the ansatzes of [20]. For some distances r, the difference is given by a factor of 5 to 6. At the same time there is good agreement of our n(r) shown in Fig. 3 with that obtained in [21].
The effective neutrino mass at the border of the cluster can be obtained asm Hereñ(r) should be computed according to Eq. (59). For m φ = 0 the expression simplifies tõ Dependence of N and R on y can be obtained using rescaling: r → r = ry, m φ → m φ = m φ /y -see Eq. (G1) in Appendix G. The rescaling means that relative effects of the density gradient and mass of φ in the equation form ν do not change with y, if m φ /y is given at a fixed value. We can write where R 7 and N 7 are the radius and density at y = 10 −7 given in the Table I. From Eq. (65), we find e.g. that at y = 10 −14 and parameters in the last column of the Table I, R = 2.4 · 10 12 cm and N = 2.3 · 10 45 , which corresponds to M = 3.3 · 10 12 g. For m φ = 0, the distribution in r simply dilates with y, and N changes correspondingly. According to Tab. II, for small N the field φ in a cluster is weak, thereforem ν ≈ m ν and p F 0 /m 0 ≈ p F 0 /m ν . With increase of N the field increases,m 0 decreases, i.e., the medium suppresses the effective neutrino mass. N 7 = 1.5 · 10 23 is critical value between the non-relativistic and relativistic cases. Further increasing N , the dependence of R and n 0 on N becomes opposite: R increases, while n 0 decreases.

C. Neutrino clusters for nonzero m φ
The dependence of properties of a ν-cluster on m φ is rather complicated. With increase of m φ the radius of scalar interaction becomes smaller and hence the interactions become weaker. Therefore, in general, one expects that the binding effect weakens and the system becomes less compact for a fixed N . As a result, the radius R increases with increase of m φ , which indeed is realized in the non-relativistic case.
In Fig. 4 we show the dependence of the radius R on N for different values of m φ . For a given m φ the dependence has a minimum R min . As in the massless case, with increase of N the radius first decreases, reaches minimum and then increases in the relativistic regime. This behavior is explained by the fact that in the relativistic case the attractive force is suppressed bym ν /E, and with increase of N the effective mass m ν decreases.  The main features of dependence of R on N are discussed below.
(i) Existence of a lower bound on N . At certain N the lines become vertical. The bound is absent for m φ = 0. For each nonzero m φ , there is a lower bound, N min , below which it is impossible to make the neutrinos confined. Since the potential energy (attraction) increases with N , to support bound system one should have large enough number of neutrinos. N min increases with increase of m φ , that is, with decrease of radius of interactions of individual neutrinos 1/m φ . The smaller the radius, the larger N should be to keep the system bounded.
The dependence of R on N in Fig. 4 differs from that in Fig. 2 of [21]. In the region of N around minimal radius and above that for m φ = 0 radius in Fig. 4 is larger than the radius in [20] and [21] and the difference increases with decrease of N . The difference is due to that [21] used an analytical solution for φ which is only valid whenñ is a constant.
Furthermore, the curves R = R(N, m φ ) of [21] do not extend far below N that correspods to the minimum of R, since the number N becomes small. Indeed, in [21] the constant α φ ∼ (0.01 − 0.1) corresponds to y 3 = (0.045 − 1.4). For N ≥ 10 2 this gives N y 3 ≥ (4.5 − 140). Our results extend to much smaller values of N y 3 , where important feature of existence of minimal value of N for m φ = 0 is realized.
The dependence of R on m φ for fixed N is also similar here to those in [21]. In the relativistic branch the radius decreases with the increase of m φ .
For m φ = 0 (blue line) the dependence of R on N is well reproduced by Eq. (45) in the non-relativistic case. Notice that in this case the scale of R is related to the neutrino mass.
(ii) Existence of a minimal radius for a given m φ . The minimum corresponds to p F 0 ∼m ν (i.e. to the effective mass of neutrinos which can be much smaller than the vacuum mass). This minimal radius, R min , in the m φ = 0 case is determined by Eq. (64). For nonzero m φ , it becomes larger because nonzero m φ suppresses the Yukawa potential exponentially.
With increase of m φ the minimum R min shifts to larger N . The value R min itself changes weakly. For large m φ the minimum of R(N ) does not correspond tom ν ∼ p F 0 , but still occurs close to transition from NR to UR regimes. Because of the shift of the curves for fixed N in the NR regime the radius increases with increase of m φ , while in the UR case it decreases. The shift can be interpreted in the following way: For fixed R the volume per neutrino is determined by m φ : . Therefore with increase of m φ , the volume per neutrino decreases and the total number of N in a cluster increases.
There is maximal value of m φ for which solution with finite radius exists. With approaching of m φ to this value the curve further shifts to larger N and minimum moves up. For m φ > m max φ the radius of interaction becomes too short to bound the neutrino system and the system disintegrates. The upper bound on m φ corresponds to the lower bound on strength of interaction: which can be compared with the bound G φ > 3.27 found in [17] from the condition of existence of minimum of total energy per neutrino for infinite medium. Notice that our bound is obtained from condition of stability of finite configuration. Let us compare the radius R with the radius of interactions R φ = 1/m φ . For small m φ there is the range of values of N in which R < 1/m φ . It is possible to obtain R = R φ at two different values of N . For m φ = 0.03ym ν the equality R = R φ is fulfilled at definite value N = 10 2 /y 3 . For m φ > 0.03ym ν , we have R > R φ for all N and the difference between them increases. R can be orders of magnitude larger than R φ and therefore in general 1/m φ does not determine the size of the cloud and fragmentation scale.
We denote by X and Y the values at horizontal and vertical axes of Fig. 4 correspondingly: The coordinate Y can be written as Y = R/R φ S −1/2 φ . Therefore the radius of a cluster over the radius of interactions equals  110) and R/R φ = 5 − 11. Thus, with decrease of strength the radius increases and becomes substantially larger than R φ . This regime was observed in [21] and called "saturation". Still 1/m φ determines the scale of sizes of ν-clusters within an order of magnitude.
In Fig. 5 we show the dependence of p F 0 (and consequently, the central density) on N for different values of m φ . For fixed m φ , p F 0 as a function of N has a maximum. With increase of m φ , this maximum shifts to larger N and increases. However, this increase slows down and is restricted by p F 0 /m ν 0.9. Explicitly, for this value With increase of m φ the central density decreases at small N (the NR regime) and it increases at larger N (UR regime). This anti-correlates with change of radius (for fixed N ) in Fig. 4. The maximal density (corresponding to the maximal p F 0 ) increases with m φ and shifts to larger N . The anticorrelation is due to the relation where b ≈ 0.2 − 0.5 is smaller than 1, since the density distribution is not uniform and the density is substantially smaller than n F 0 in peripheral regions. Thus, the overall change of characteristics of cluster (for fixed N ) with m φ is the following: In the NR case central density decreases and radius increases (cluster becomes less compact). In the UR case, the radius varies non-monotonically with m φ , as shown by the relativistic branches of the curves in Fig. 4. These characteristics could be important for cluster formation processes.
For fixed m φ the bound system exists for N above certain minimal value. With the increase of N the central density (and p F 0 ) increases quickly and therefore the radius of cluster decreases. Then the density reaches its maximum and with further increase of N the density decreases. Correspondingly the radius of cluster increases.
Notice that strengths of interactions in the case of two neutrino bound states λ [see Eq. (19)] and for the bound neutrino system [Eq. (66)] have different functional dependence on the parameters. From Eqs. (19) and (66) we obtain the following relation between the two strengths: For y ≤ 10 −7 and λ 2 ≥ 0.7 required to have 2ν-bound state we find S φ > 6.3 · 10 16 . That is, if 2ν-bound states exist also ν-clusters should exist. The opposite does not work: For S φ ∼ 10 2 we have λ ∼ 1.6 · 10 −15 , and so the 2ν bound states are not possible. Both the bound states and bound systems exist for very small masses of scalars. If m φ < m cr φ (for fixed y), the strength increases and according to Fig. 4, the ranges of X and Y expand. Correspondingly, N = Xy −3 = 10 25 X/X cr , and therefore clusters with much smaller number of neutrinos are possible. E.g. for S φ = 10 4 , X can be as small as 0.2, and therefore N = 2 · 10 20 . X becomes a free parameter in the interval 0.2 − 10 4 . The allowed range of Y expands in both directions: Y can be as small as 0.3Y cr , and therefore the radius reduces down to 0.7 km. The maximal value of Y can be much bigger than Y cr and it increases with increase of strength (decrease of m φ ).
For S φ = 10 4 (orange line in Fig. 4) Y = (30 − 380) and correspondingly, R = (0.6 − 7.6) km. If also y decreases and m φ decreases, so that S φ does not change, the line Y = Y (X) is the same, but parameters of cluster increase as N ∝ y −3 , R ∝ y −1 .
For fixed y and m φ stable bound systems are situated along the corresponding lines S φ = const in Figs. 4 and 5. This, in turn, fixes the interval of X and Y (hence N and R): As an example, for y = 10 −20 and m φ = 10 −23 eV we obtain S These relations are valid for m φ 1/R. With decrease of y (and fixed N ), R increases as 1/y and the density decreases as y 3/2 .

VI. FORMATION OF NEUTRINO CLUSTERS
Let us consider the formation of neutrino clusters from the uniform background of relic neutrinos in the course of cooling and expansion of the Universe.
We find that for the allowed values of y and m ν , the typical time of cooling via φ bremsstrahlung and via annihilation νν → φφ are much larger than the age of the Universe (see Appendices I and J). Therefore the formation of neutrino clusters differs from the formation of usual stars. 5 As we will show in subsections VI C and VI D, for the interaction strength above a certain critical value it has a character of phase transition at which the kinetic energy is transformed into the increase of the effective neutrino mass. The mechanism was first suggested in [17] and here we explore it further.
The phase transition leads to fragmentation of the cosmic neutrino background. Here we present some simplified arguments and analytic estimations. A complete solution of the formation problem would require numerical simulation of evolution of the ν background. Such a simulation is beyond the scope of this paper.

A. Maximal density and fragmentation of the ν background
The key parameter that determines formation of ν-clusters is the central density given by Eq. (68). According to Fig. 5, there is a maximal value of p F 0 , and consequently, a maximal value of the density n max 0 which is determined by the mass of neutrino in Eq. (68). After the fragmentation of a uniform relic ν-background (see below), the central density of a cluster stopped decreasing. Consequently, at the epoch of fragmentation the density was n frag ≤ n max 0 . Thus, we can take n frag ≤ 10 9 cm −3 .
In turn, this density determines the epoch of fragmentation: where n rel = 113 is the number density of one neutrino mass state that would be in the present epoch in the absence of clustering. That is, formation of ν-clusters starts after the epoch of recombination (z = 1100, T = 0.1 eV). At this epoch (z f +1 = 215) the neutrino average momentum was p ∼ 0.02 eV. So, fragmentation could start when p ∼ m ν , which is not accidental since n max 0 is determined by m ν only.

B. Evolution of the effective neutrino mass
Dependence of the effective massm ν on T plays crucial role in the formation of neutrino bound systems. If neutrinos are the only source of the field φ, for the allowed coupling constant φ cannot be thermalized in the early Universe and the amount of φ particles is negligible. Evolution of the scalar field in expanding universe is described by Eq.
and νν is given in (10). Thus, the expansion term regularizes the solution for m φ → 0. The energy density in the field: In what follows we assume that H m φ , so that expansion can be treated as slow change of density and temperature in the solution of equation without expansion.
Neglecting the Hubble expansion term, Eq. (72) gives for uniform medium From Eq. (74), we obtain the equation for the effective mass of neutrino: Here we use the thermal distribution with T = 0, which is relevant for T > m ν , while in [17] the distribution of degenerate Fermi gas was used. Qualitatively, the results are similar. In Fig. 6, we show dependence ofm ν on T for different values of strength of interaction. In the ultrarelativistic limit, T m, we obtainm where in general, the integrals I n have been introduced in Eq. (48). From Eq. (75) we find That is, the effective massm ν ∝ T −2 decreases as T increases and becomes negligible in the early Universe.
In the non-relativistic case, T m ν , the solutions is That is, as T → 0,m ν → m ν . These analytic results agree well with results of numerical computations shown in Fig. 6. 6 Notice that results for the degenerate case can be obtained by substituting Note that the relativistic regime is determined bym ν and not m ν : T >m ν . The effective massm ν increases as 1/T 2 till T ≈m ν . According to Eq. (76) this equality gives Below this temperature the growth ofm is first faster and then slower than 1/T 2 . The change of derivative occurs at maximum of the relative energy in the scalar field (see below).m(T ) converges to the non-relativistic regime atm rel ν (T ) = m ν . This gives, according to Eq. (76),

C. Energy of the ν-φ system and the dip
Let us compute the energy density in the scalar field, ρ φ , in neutrinos, ρ ν , and the total energy density defined as as a function of T . Expansion of the Universe can be accounted for by considering the average energy per neutrino: where the number density of neutrinos equals The energy density in the scalar field is given by Here φ can be expressed via the effective mass, φ = (m ν − m ν )/y, and therefore So, ρ φ is directly determined by the deviation of the effective mass of neutrino from the vacuum mass. The energy per neutrino equals Using the known form ofm ν (T ), Eq. (75), we compute the energy density φ . In the ultra-relativistic limit, according to Eq. (82), the energy density in φ converges to the constant: where In the non-relativistic limit According to Eq. (84), as the Universe cools down, φ first increases and then decreases. The maximum of φ is achieved when rel φ ≈ nr φ , which happens at χ = 2, or Hence the maximum is max = m ν /2. When S φ increases (correspondingly m φ decreases), the maximum max does not change significantly, but shifts to smaller T . In Fig. 7, we show the dependence of φ /m ν on T (dashed lines). As can be easily checked, the analytic expressions above agrees well with the numerical results.
The energy density of neutrinos is given by It can be rewritten as Here the second integral coincides with integral in νν and φ. The dependence of ρ ν (m ν ) on T can be found explicitly.
In the relativistic case, we have and rel ν = as expected. The dependence of ν on T is shown in Fig. 7 (dotted lines). Solid lines in Fig. 7 correspond to the total energy (80). Scalar contributions shift the minimum of tot to slightly larger T with respect to the minimum of ν .
The key feature of the tot dependence on T is development of the dip at some temperature below m ν , denoted by T dip .
According to Fig. 7, tot becomes smaller than m ν when T m ν /3. As T decreases below T dip , the energy tot converges to m ν from below. The dip essentially coincides with the transition region form ν (Fig. 6). Without the Yukawa interaction (or if the Yukawa interaction is too weak, as illustrated by the red line for S φ = 0.04 −2 = 625 in Fig. 7), tot would decrease as tot ∝ T in the relativistic case; at T m ν /3, the curve flattens and tot converges to m ν . For the entire range it would always be tot > m ν which corresponds to unbounded neutrinos. For larger S φ (other curves in Fig. 7), the dip appears. The appearance of the dip with tot < m ν means that neutrinos can form bound systems with the average binding energy equal to m ν − tot . As the interaction strength increases, the dip shifts to lower temperatures and becomes deeper. Correspondingly the binding energy becomes stronger.
The critical value of the strength needed for existence of the dip, S c φ , is determined by the condition that temperature of beginning of deviation of rel ν from linear decrease equals T rel ≈ m ν /3. This temperature coincides with T rel of deviation ofm ν from its 1/T 2 behavior. The condition (78) can be rewritten as Then for T rel /m ν ≈ 1/3 it gives S c φ = 580, which is close to S φ = 625 for the red curve in Fig. 7. In actual numerical calculations, we find that the dip appears roughly at S c φ ≈ 700.

D. Fragmentation
Here we consider S φ > S c φ so that the dip occurs. Above T dip , cosmic neutrinos evolve as a uniform background with decreasing density and effective temperature. For T < T dip , further cosmological expansion assuming the uniform neutrino background would imply that the energy of the system increases. Since there is no injection of energy into the ν-φ system from outside, it is energetically profitable that the uniform neutrino background fragments into parts (clusters) with the temperature T ≈ T dip , which determines also the corresponding number density. Further evolution of fragments will proceed, subject to the dynamics of ν clusters. The expansion of the Universe then increases the distance between fragments.
The highest temperature of fragmentation corresponds to T f ∼ 0.3m ν ∼ 0.03 eV. This temperature is achieved at redshift z f = 200 when the number density of neutrinos was n f = 9 · 10 8 cm −3 . (86) In the process of formation of final configuration this density could further increase in the central parts of objects. This value agrees well with the one in Sec. VI A. The biggest possible object after fragmentation, which satisfies minimal (free) energy condition, would have the size D 0 ≤ 1/2D U (z f ) in one dimension, where D U (z f ) is the size of horizon in the epoch of fragmentation, D U (0) = 4.2 Gpc. Therefore in three dimensions one would expect that 2 3 = 8 such objects could be formed, and for estimations we will use ∼ 10 objects. At z f = 200 we find D 0 = 10 Mpc. In the present epoch distance between the objects is z f times larger, that is ∼ 2000 Mpc.
Using the number density density (86) and the radius R f = 5 Mpc we obtain the total number of neutrinos in a cluster Correspondingly, the mass is M = 4 · 10 17 M . These biggest structures can be realized for m φ = 4 · 10 −30 eV and y = 4 · 10 −27 .
If m φ and y are much larger, these structures do not satisfy conditions for stable bound systems. The number of neutrinos and radii should be much smaller. One can explore two possible ways of formation of small scale structures: 1. The neutrino background first fragments into the biggest structures, thus satisfying the energy requirement. These structures are unstable and further fragment into smaller parts. The secondary fragmentation can be triggered by density fluctuation of the neutrino background itself or via gravitational interactions with already existing structures like the Dark Matter halos.
If the secondary fragmentation is slow, the overall final structure could consist of superclusters (originated from primary fragmentation) and clusters within superclusters.
2. The uniform neutrino background may immediately fragment into small structures (close to final stable bound states) with R = O(1/m φ ). These structures may interact among themselves.
As we established, for given m φ and y (i.e. fixed strength) there is a range of stable configurations with different N and R. E.g. for S −1/2 φ = 0.01, N can be within 4 orders of magnitude. What is the distribution of clusters with respect to N ? The answer depends on details of formation, in particular, on the surface effects. One can consider several options: (i) a flat distribution within the allowed interval; (ii) a peak at N = X(n F 0 )/y 3 , where n F 0 is the density at the beginning of fragmentation which, in turn, is determined by the strength, etc.
Considering T min as the temperature of the beginning of fragmentation we can find using Figs. 6 and 7 the physical conditions of medium which determine the initial states of pre-clusters. The kinetic energy per neutrino is (in the unit of m ν ) kin = 3.14T min /m ν . The kinetic energy turns out to be much bigger thañ m(T min ): kin m(T min ).
So neutrinos are ultra-relativistic and therefore our usage of the relativistic distribution is justified. Correspondingly, the total energy in neutrinos nearly coincides with the kinetic energy: ν ≈ kin .
The number density of neutrinos is given by n min 0 = 2 · 10 6 cm −3 which would correspond to the Fermi momentum p F /m ν = 0.1. This momentum is smaller than the kinetic energy, so the cooling is needed to reach degenerate configuration.
Further results and statements should be considered as a possible trend since to construct the figures we used approximation of uniform infinite medium and relativistic distributions of neutrinos over momenta. Below T dip , the dynamics of finite size objects should be included.
As T deceases below T min , the kinetic energy of neutrinos decreases while the effective massm increases. This means that the kinetic energy transforms into the increase of the mass as well as into the energy of scalar field. Thus, for T /m ν = 0.03 we obtain kin = 0.09,m/m ν = 0.03, i.e. the ratio kin m ν /m = 3; for T /m ν = 0.02 the ratio becomes 0.3. According to Figs. 6 and 7, the minimum of total neutrino energy, ν , is achieved at a temperature lower than T dip . In the minimum, for S −1/2 φ = 10 −3 , ν = 0.1 andm/m ν = 0.07, i.e. neutrinos become nonrelativistic.
As T further decreases, both ν andm quickly increase withm → ν . This indicates that the kinetic energy is transformed intom as well as into φ . One can guess that the phase transition ends when φ reaches maximum which corresponds to the maximal potential energy. This happens at a slightly lower temperature than the temperature of ν reaching its minimum.
The period of (primary) fragmentation can be estimated as the time interval which corresponds to the width of the dip: where t ≈ 2 3H(T ) (for matter dominated era) and H is roughly proportional to T 2 . The specific form can be determined by H 2 = 8πρ total /(3m 2 pl ) where ρ total is the total energy density of the Universe and m pl is the Planck mass. The two temperatures T 1 and T 2 correspond to the maximum of the dashed curve in Fig. 7 and the minimum of the solid curve respectively. This gives e.g. ∆t = 0.014τ U for S −1/2 φ = 10 −3 and ∆t = 0.24τ U for S −1/2 φ = 10 −4 , where τ U is the present age of the Universe. Note that the estimations based on Fig. 7 may not be reliable since in the region of temperatures below T dip neutrinos become non-relativistic and the dynamics of clusters starts to dominate.
According fo Fig. 7, T dip decreases as S φ increases. This means that for larger S φ , fragmentation starts at later epochs and lower densities. For S

E. Virialization
In the range of S φ = 70 − 700, the dip disappears and the phase transition is absent. Such interaction strengths still allow for stable bound systems of neutrinos. In this case, ν clusters might be formed from the growth of local over-density to the non-linear regime and then go also through virialization, similar to the formation of dark matter halos.
Quantitatively, results for Yukawa forces (if m φ is sufficiently small) can be obtained from the gravitational results by replacing where G is the gravitational constant and M = m ν N is the total mass. The key difference from the gravitational case is that the strength of Yukawa interaction is many orders of magnitude larger than the strength of gravity. With this analogy we can consider the following picture of structure formation. 1. The initial state is the neutrino background with some primordial or initial density perturbations (primordial clouds) with ∆n p /n p 1 and sizes R p . These perturbations may be related to the density perturbation of DM in the Universe, since neutrino structures are formed latter.
2. As the Universe expands, the size of clouds and distance between them increase simultaneously as the scale factor R ∝ d ∝ a. Correspondingly, the density decreases as a −3 and ∆n/n is constant. This is the linear regime.
3. At certain epoch, when T < 0.3m ν the evolution becomes non-linear: the attraction becomes more efficient than expansion and the increase of the size of clouds becomes slower than the expansion of the Universe. This happens provided that the potential energy of the cloud is substantially larger than the total kinetic energy. Then the distance between two clouds increases faster that size of clouds. Also the decrease of the density slows down. 4. At some point (for which we denote the scale factor by a max ) the density stops decreasing, and after that (a > a max ), the system starts collapsing. In this way a bound system is formed and its evolution is largely determined by internal factors. The distance between two clouds continues increasing.
Since the system is collisionless and there is no energy loss (emission of particles or and waves) the contraction has a character of virialization 7 , which starts from the stage that kinetic energy of the neutrino system becomes significantly smaller than the potential energy: During virialization, the kinetic energy increases (the effective T increases, although the distribution may not be thermal). The process stops at a = a vir when virial equilibrium is achieved, which also corresponds to hydrostatic equilibrium. The radius decreases by a factor of 2 so that the density increases by a factor of 8. Approaching the equilibrium the cluster may oscillate. At this point, the formation of the cluster is accomplished.
In what follows we will make some estimations and check how this picture matches the final configurations of clusters we obtained before. We will discuss bounds on initial parameters of the perturbations , formation procedure, epochs of different phases, etc, which follow from this matching.
The total energy is conserved during virialization. When the system reaches virial equilibrium, the virial theorem implies that Here E vir K and E vir V are the total kinetic and potential energies after virialization. Thus, the ratio of energies increases.
Although quantitative studies of virialization involve N -body simulation, the orders of magnitude of various characteristics can be obtained from simplified considerations. For gravitational virialization, it is known that 8 the virialization time (roughly the time of collapse) driven by gravity equals where E K and E V are the total kinetic and potential energies (E K + E V < 0) of the system. The time of m φ -driven virialization can be obtained from Eq. (92) by making the replacement (89): The contraction stops when virial equilibrium is achieved and we assume that the equilibrium has the same form as in the gravity case (91). Since the total energy is conserved during virialization, The total energy in the beginning of virialization can be written as where in the last equality we have taken ξ E = 1/3. Inserting this total energy with |E V | given in (52) into Eq. (93), we obtain the time of virialization A system of collisionless (i.e. no additional interactions except for gravity) particles can collapse from a large cold distribution (or more strictly, a distribution with low mean kinetic energy since it is usually not a thermal distribution) to a smaller and hotter one. This processes is known as virialization. 8 See, e.g., http://www.astro.yale.edu/vdbosch/astro610_lecture8.pdf or Ref. [26].
which does not depend on N . Numerically, we get t vir 3.4 · 10 −6 sec 10 −7 y which implies that for T /m ν = 1, the virialization takes t vir 3.4 · 10 −6 sec. As aforementioned, after the virialization process, the radius decreases by a factor of 2 and correspondingly, the density increases by a factor of 8. Hence in the example we discussed at the end of the previous subsection, this gives R vir = 1.3 km, and n vir = 7 · 10 10 cm −3 .
As a result of virialization, the distribution reaches virial equilibrium with the virialized energies given by A halo of neutrinos is formed. In the absence of energy loss, its number density remains constant without being diluted by a −3 (a is the scale factor) due to the cosmological expansion. That is, a halo of high-density cosmic neutrinos would be frozen till today. Let us estimate temperatures of the ν-clusters at different phases of formation. We can assume that deviation from the linear regime occurs at T ∼ m ν . Also we assume that before collapse the state of the system can be approximately described by thermal equilibrium and the border effects can be neglected. Therefore, we can use results obtained in sect. IV C. The expression for temperature (54) can be rewritten as Let us take y 3 N = 1 (see implications in the Fig. 4, 5) then T /m ν = 1 (the scale factor a N L ) corresponds to ξ E = 189. For ξ E = 1/3 we obtain T /m ν = 1.8 · 10 −3 . For y 3 N = 32, T /m ν = 1 corresponds to ξ E = 19 and ξ E = 1/3 is achieved at T /m ν = 1.8 · 10 −2 . If y 3 N = 10 3 , T /m ν = 1 is realized at ξ E = 1.9, and ξ E = 1/3 at T /m ν = 0.18.

F. Neutrino structure of the Universe
The formation of the neutrino clusters leads to neutrino structure in the Universe with overdense areas and voids. For simplicity we consider that all the clouds have the same number of neutrinos N .
In the present epoch the distance between neutrinos without clustering equals If neutrinos cluster into clouds with N neutrino in each, the distance between (centers of) these clusters is Taking, e.g., N = 6 · 10 22 , we obtain d = 78 km which is much larger than the radius of each cloud, R = 0.62 km given by Eq. (64).
The existence of such a structure could have an impact on the direct detection of relic neutrinos in future experiments such as the PTOLEMY [27,28]. The result also depends on the motion of the clusters with respect to the Earth. The motion of the Earth and the solar system will average the effect of these small structures. According to Eq. (45), in the non-relativistic regime the radius of cluster decreases with increase of N as N −1/3 . Therefore the ratio of spacing (distance) and radius increases as ∝ N 2/3 : Taking the above example (d = 78 km, R = 0.62 km), we have d/R = 126. So, voids occupy most of the space. With decrease of y, the radius of cluster increases and N also increases (for fixed maximal density).
As a result, d/R = const which is determined by the neutrino mass.
Let us estimate how small scale clusters may affect the direct detection of relic neutrinos assuming that clouds reached their final configurations. The distance needed for a detector to travel through to meet a ν-cluster, in average, can be estimated as L = (σ c n c ) −1 , where σ c ≈ πR 2 is the cross section of the cluster and n c = d −3 is the number density of clusters. Consequently, Numerically, in our example, we obtain from (96) L = 4 · 10 5 km. Since the circumference of the Earth is 4 · 10 4 km, the Earth detector will cross the cluster every 10 days in average. Here we have not taken into account motion of clusters.
If the time of observation is much longer than 10 days, the total rate of events should be the same with and without clustering of neutrinos. This can be seen by noticing that the interaction rate when the detector is in the cloud is enhanced by a factor of n c /n 0 while the rate of the detector being in a sphere of neutrinos is reduced by a factor of n c /n 0 . Hence the time-averaged result will be the same as in the non-clustered case.
If the clusters are very big (e.g. comparable with the solar system or more) which is realized for very small coupling constant y, the rate will be modified if the detector is in a void or in the overdense region of cluster.

VII. CONCLUSIONS
Yukawa interactions of neutrinos with a new light scalar boson φ are widely considered recently. One of the generic consequences of such interactions is the potential existence of stable bound states and bound systems of many neutrinos (ν-clusters). Our study addresses these questions: whether and when do bound systems exist? what are the characteristic parameters of these systems and how are they determined? what are the observational (cosmological, laboratory, etc.) consequences of ν-clusters?
We find that for two-neutrino bound states the Yukawa coupling y and the mass of the scalar m φ should satisfy the bound λ ≡ y 2 m ν /(8πm φ ) > 0.84. Eigenstates and eigenvalues of a 2ν bound system are obtained by solving the Schrödinger equation numerically. Analytic results (including the critical value of λ as well as the binding energy and the radius of the system) can be obtained approximately using the variation principle. The binding energy in a two-neutrino state is always small compared to the neutrino mass due to smallness of y. For allowed values of couplings y and masses of scalar m φ , the bound state of two neutrinos would have the size greater than 10 12 cm. Bound states of sub-cm size are possible for 10 keV sterile neutrinos with coupling y > 10 −4 . As long as the Yukawa coupling is below the perturbativity bound, neutrinos in a two-neutrino bound state are non-relativistic.
For N -neutrino bound systems, the ground state corresponds to a system of degenerate Fermi gas in which the Yukawa attraction is balanced by the Fermi gas pressure. For non-relativistic neutrinos the density distribution in the ν-cluster is described by the Lane-Emden equation. We re-derive this equation and solve it numerically. As a feature of the Lane-Emden equation, the solution always has a finite radius R. The radius R decreases as N increases in the non-relativistic regime.
Unlike 2ν bound states, the N -neutrino bound system can enter the relativistic regime with sufficiently large N . We have derived a set of equations to describe the ν-clusters in relativistic regime. In the nonrelativistic limit they reduce to the L-E equation. In the relativistic regime, the Yukawa attraction acquires the m ν /E suppression which leads to a number of interesting consequences: • The effective number density of neutrinos,ñ, is suppressed. It deviates from the conventional number density n significantly when the neutrino momentum is large; • The Fermi degenerate pressure is able to balance the attraction for any N , thus preventing the collapse of the system in contrast to gravity; • The radius of system increases with N , which implies that relativistic N -body bound states cannot be much more compact than the non-relativistic ones.
• As N increases the radius of a ν-cluster first decreases (in the non-relativistic regime) and then increases (in the relativistic regime), attaining a minimum when p F 0 is comparable to the effective neutrino mass.
• There is a maximal central density ∼ 10 9 cm −3 , which is determined by the neutrino mass.
• For a given m φ there is a minimal value of N y 3 for which stable configurations can be formed. This minimal value increases as the interaction strength S φ ≡ • For a given interaction strength, there is a minimal radius of the cluster.
• The minimal interaction strength (maximal m φ for fixed y) required to bind neutrinos in a ν-cluster is S φ 70.
We have investigated possible formation of the ν-clusters from the uniform relic neutrino background as the Universe cools down.
• For the interaction strength S φ 700, a dip (see Fig. 7) appears in the evolution of the total energy of the ν-φ system at some temperature below m ν . In this case, the formation of neutrino structures has a character of phase transition. The dip leads to the instability of the uniform ν background, causing its fragmentation when T further decreases. The typical size of final fragments is R ∼ (0.1 − 10)/m φ .
A detailed picture of fragmentation should be obtained by numerical simulations, which is beyond the scope of this work.
• For 70 S φ 700, the dip is absent, and hence fragmentation does not occur, though stable neutrino clusters can exist. In this case neutrino clusters might be formed via the growth of local density perturbations followed by virialization, in a way similar to the formation of DM halos, though the energy loss of neutrino clusters due to φ emission and neutrino annihilation is negligible. In this case, a new cooling mechanism would be required.
• For S φ 70, stable neutrino clusters cannot exist.
If clusters are formed at redshift z f (which can be as large as 200) and retain the same size until today, the neutrino structure of the universe would show up as ν-clusters which are z 3 f times denser than the homogeneous neutrino background and voids which are z f times larger than clusters. This leads to important implications for the detection of relic neutrinos.
While conclusion of formation of the neutrino clusters via fragmentation at S φ 700 is robust, details of fragmentation, the distribution of clusters with respect to N , and the evolution of clusters require further studies. the probability of this particle appearing at position x is P(x) ≡ |Ψ(x)| 2 and we have chosen w(p) in a way that which is always possible as long as ∆x∆p 1. In this way, we can say that the particle is located approximately at x 0 with an approximate momentum p 0 .
Next, let us consider the normalization: Substituting Eq. (A3) in P(x) ≡ |Ψ(x)| 2 and then in Eq. (A5), we have where it is useful to note that With the above setup, we can evaluate the mean value of ψψ and ψ † ψ in the presence of the single particle state |w . First, note that Hence for w|ψ † ψ|w , using Eq. (A1) and Eq. (A8), we have where the spacetime-generalized wavefunction, Ψ(x), is defined similar to Eq. (A3) except that e ip·x is replaced by e −ip·x . Similarly, for w|ψψ|w , we obtain As has been formulated in Eq. (A4), the particle at t = 0 is localized in the region x ∈ [x 0 − ∆x, x 0 + ∆x]. So it is useful to inspect the average value of w|ψψ|w and w|ψ † ψ|w . At t = 0, we have which is exactly the particle number density of a single particle distributed in the small volume ∆x 3 . If ψ † is replaced with ψ in Eq. (A11), at the last step we would have u s (p)u s (p) → 2m ψ [see Eq. (A7)] and hence , which in the non-relativistic regime is equivalent to the number density but in the relativistic regime is suppressed.
The above calculation can be straightforwardly generalized to N particles with different momentum and difference spins. Each of them has an independent wave-packet function w i (p). In this case, we sum over the contributions together: which can be regarded as the momentum distribution function of the N particles. Then the summation gives Appendix B: Potential issues on the binding energy and the field energy There is a potential issue regarding the binding energy of two particles in Sec. II: for two particles, the binding energy does not contain a factor of two. To address this issue, here we discuss the two-body problem from both the classical and modern (field theory) perspectives.
In the classical picture, we assume that each particle can be infinitely split to smaller particles. Denote the two particles as A and B. If we could remove an infinitely small percentage ( , 1) of particle B to r = ∞, the energy cost would be y 2 4π 1 r e −rm φ (due to the attraction of A) plus the energy cost to overcome the attraction of B. The latter part will be canceled out when all the small fractions of B are re-assembled at r = ∞. So the total energy cost (i.e. binding energy) to separate A and B is Hence we conclude that the total binding energy for a two-body system should take N = 1 in Eq. (14). In addition to the classical picture, from the perspective of field theory, an equivalent interpretation is that the binding energy is stored in the field φ. According to Eq. (6), the Hamiltonian for static φ is For the case of two particles, we have where x A and x B are coordinates of A and B, and r A,B = |x| 2 − |x A,B | 2 . Although H computed in this way is divergent, by comparing two cases with |x A − x B | = r and |x A − x B | = ∞, the difference of H is finite. Plugging Eq. (B4) into Eq. (B2) and using (∇φ) 2 d 3 x = − φ∇ 2 φd 3 x, it is straightforward to obtain the energy difference: This is consistent with Eq. (B1) in the classical picture.
Note that here we assume that V depends on r only. Eqs. (C6) and (C7) allows us to separate the motion of R with that of r: In this work we are only concerned with the part of ψ(r), which can be expanded using spherical harmonics. Plugging Eq. (21) to we obtain − 1 2µ DefineṼ ≡ 2µV (r) + l(l + 1) Then the equation is simplified to − u (r) +Ṽ (r)u(r) =Ẽu(r).
Given any form ofṼ (r), in principle Eq. (C12) can be numerically solved. Since r varies from zero to infinity, it is actually more convenient to make the following variable transformation: so that v(ω) is a function on a finite interval ω ∈ [0, π/2). Here R in principle can be an arbitrary length scale but in order to improve the efficiency of the numerical method it should be set at the typical length scale of the wave function that is being inspected. After the transformation, Eq. (C12) becomes This equation can be solved by identifying it as the problem of find eigenvalues and eigenvectors of a large matrix. If we divide the interval [0, π/2) into N segments and ω takes one of these values: where dω = π 2N , then d 2 dω 2 can be regarded as an operator in this discrete space, represented by the (N + 1) × (N + 1) matrix below: AndṼ (ω) and cos 4 ω are represented by diagonal matrices: cos 4 ω → diag cos 4 (0), cos 4 (dω), cos 4 (2dω), · · · , cos 4 ( π 2 ) .
In the large-N limit, by solving eigenvalues and eigenvectors ofH defined below we can obtain energy levels and wave functions. In laboratory experiments, a light scalar that exclusively couples to neutrinos can be constrained by particle physics processes such as meson decay (π ± , K ± ) [29,[35][36][37], Z invisible decay [31], and double beta (ββ) decay [38][39][40][41][42][43]. Since the typical momentum transfer scales of these processes varying from 10 2 GeV to MeV are well above the mass of φ considered in this work, the constraints from these processes on the coupling y are almost independent of m φ when m φ is lighter than m ν . Most of these particle decay bounds only apply to the electron and/or muon flavor, except for Z invisible decay which is relevant to all flavors (For couplings to ν τ , τ decay is less restrictive than Z decay [31]). In addition, the ββ decay bound is only valid when φ is coupled to two electron neutrinos instead of one neutrino and one anti-neutrino. The specific values of aforementioned bounds are summarized in Tab. II. High precision measurements of neutrino scattering in principle could be sensitive to such interactions via bremsstrahlung or loop processes, which compared to the standard processes are typically suppressed by a factor of y 2 /16π 2 ∼ 10 −2 y 2 . Current neutrino scattering experiments can only measure the cross sections at percent-level precision [44,45], corresponding to an upper limit of y 1, which is negligibly weak. Long-range force searches (also referred to as the fifth force searches in the literature) based on precision test of gravitational laws (e.g. the inverse-square law, weak equivalent principle) are sensitive to light mediators coupled to normal matter (electrons or baryons) [46][47][48]. Although in the presence of φ-neutrino couplings, φ-electron couplings can be induced at the one-loop level [49,50], such loop-induced couplings are typically highly suppressed by the large hierarchy between neutrino masses and the electroweak scale. Hence bounds from long-range force searches can be neglected in this work.
Astrophysical and cosmological bounds, by contrast, are much more restrictive. Observations of neutrinos from the supernova event SN1987A have been employed in the literature to set various constraints from different considerations. The most commonly considered effect is energy loss, which leads to the constraint y < 3 × 10 −7 or 2 × 10 −5 < y < 3 × 10 −4 [33]. If the new interaction is able to convert ν e to ν e or to neutrinos of other flavors, then another effect known as deleptonization can also put a strong bound, y < 2 × 10 −6 [34]. Other supernova bounds [51][52][53][54] assuming the presence of quark couplings, heavy masses, or electromagnetic decays do not apply here. In addition, for y > 4.6×10 −6 the scalar could be thermalized in the early universe, causing observable effects in the Big Bang Nucleosynthesis (BBN) [9]. For smaller y, the scalar cannot be thermalized but could affect the evolution of cosmological perturbations and hence be constrained by CMB data [32]. Theses bounds are also summarized in Tab. II.

Appendix E: Equivalence between chemical equilibrium and force balance
In this appendix, we show that the balance of the two forces, F Yuk = F deg in Eq. (39), is equivalent to chemical equilibrium of the system (equilibrium of particle numbers in the presence of external forces). Indeed, the chemical equilibrium of system means that where µ is the chemical potential in the Fermi-Dirac distribution f = [e (E−µ)/T + 1] −1 with the neutrino energy E modified by φ: The Fermi energy, E F , is defined as Eq. (E2) with p = p F . Since it is the maximal energy of a particle in the distribution, E F = µ, Eq. (E1) is equivalent to dE F /dr = 0. So we have or y(m ν + yφ) dφ dr + p F dp F dr = 0.
Let us show that in nonrelativistic limit dφ dr and dp F dr can be related to the Yukawa force and the degenerate pressure respectively.
To make use of the classic concept of forces, we need to take the non-relativistic limit in which yφ m ν so that yφ can be neglected in the first term of (E4). The gradient dφ dr is related to the Yukawa force as For non-relativistic neutrinos without interactions, we have E = m 2 ν + p 2 ≈ m ν + p 2 2mν and hence E − µ ≤ 0 corresponds to m ν + p 2 2mν − µ 0. When p = p F , E − µ ≤ 0 is saturated. In the presence of the φ-induced potential V m ν , according to Eq. (9), E − µ ≤ 0 corresponds to m ν + p 2 2mν + V − µ 0. Hence it is convenient to define an effective chemical potential to absorb m ν and V : When µ eff is positive, we have p F determined by Note that µ eff sometimes can be negative, for which we should take p F = 0: Consider Eq. (12) where the r.h.s depends on the local number density n(x), which according to Eq. (31) is determined by the local Fermi momentum p F (x), which according to Eq. (F2) or (F3) depends on V . So eventually the r.h.s of Eq. (12), which as a source term generates the potential V , is determined by V itself: It is more convenient to express Eq. (F4) in terms of µ eff using Eq. (F1): For spherically symmetric distributions and m φ = 0, it becomes whereκ ≡ y 2 (2m ν ) 3/2 6π 2 , γ = 3 2 . (F7) When µ eff is expressed in terms of n, one obtains Eq. (42).
Appendix G: Numerical details of solving Eqs. (60) and (61) It is useful to note that we can remove the dependence of the differential equation on y by defining ∇ ≡ ∇/y, r ≡ yr, m φ ≡ m φ /y, and rewrite Eq. (60) as where Eq. (61) is equivalent tom wherem R ism ν at r = R (or r = R ≡ Ry). This allows us to substitutem ν (r) → m 2 R − p 2 F (r) in Eqs. (G2) and express it as an equation of p F (r).
The initial values in Eq. (62) are determined as follows. We note thatm 0 cannot be set arbitrarily because at r → ∞ it should match m ν . Technically, for an arbitrarilym 0 , quite often one gets a solutions with p F (r) monotonically increasing or oscillating as r increases, which is certainly not a physical solution. We find that in practice it is more convenient to usem R instead ofm 0 . Hence, in our code implementation, we first set values of (m R , p F 0 ) and then usem 0 = m 2 R − p 2 F 0 to obtain (m 0 , p F 0 ). After solving the differential equations for 0 < r < R, one can obtain N ≡ N y 3 andm ν (r → ∞) using to the following integrals: With the above setup, we implement the code in the following way: With the above result of |M| 2 , the total cross section of annihilation can be computed by Here q 2 = (p 3 − p 1 ) 2 ≈ m 2 ν (−1 + vc θ varies in the range m 2 ν (−1 − v) q 2 m 2 ν (−1 + v). Integrating q 2 in the kinetically allowed domain, we obtain the result in Eqs. (J1) and (J2).
Appendix I: Energy loss due to φ bremsstrahlung Due to φ bremsstrahlung and neutrino annihilation, the energy loss can cause further contraction of the virialized halo. Let us first consider the bremsstrahlung process in two-neutrino collision: ν + ν → ν + ν + φ with φ exchange in t channel. The differential cross-section can be estimated as where q is the momentum carried by the t-channel φ mediator and v is the relative velocity of two colliding neutrinos. The energy taken away by φ equals roughly E φ |q 2 |/2m ν , so that according to (I1) the energy loss is proportional to Here q 2 max and q 2 min are the maximal and minimal values of |q 2 |. The former is determined by kinematics, q 2 max ∼ (m ν v) 2 , and for the latter we took q 2 min ≈ 0. Consequently, the rate of energy loss of a given neutrino can be estimated as If q 2 r −2 mean , a given neutrino interacts coherently with neutrinos within the radius r coherent ∼ 1/p , and the differential cross section can be enhanced by a factor ofÑ 2 , whereÑ ∼ 4 3 πnr 3 coherent is the number of neutrinos within the coherence volume. On the other hand, takingÑ neutrinos as a target, the effective number density of scatterers is reduced to n/Ñ . As a result, the rate in Eq. (I3) still increases by a factor of NÑ = O(10).
The energy loss rate per neutrino reads The time it takes for the system to cool down and lose half of its kinetic energy can be estimated by where E vir K is determined in Eq. (95). We take the virialized radius R vir ∼ R in /2 where R in is given Eq. (55). Then using n = N/( 4 3 πR 3 vir ) and v 2 = E , which is much larger than the age of the Universe τ universe = 4.35 × 10 17 sec. dp νν The solution is where p 0 F is the Fermi momentum in the beginning of the cooling phase. The annihilation lifetime τ anni can be defined as the time during which the Fermi momentum reduces by a factor of 2: p F (τ anni ) = p 0 F /2, which corresponds to the density being reduced by a factor of 8. Then from (J7) and (J8) we obtain