Bose-Einstein Condensate Dark Matter That Involves Composites

By improving the Bose-Einstein condensate model of dark matter through the repulsive three-particle interaction to better reproduce observables such as rotation curves, both different thermodynamic phases and few-particle correlations are revealed. Using the numerically found solutions of the Gross-Pitaevskii equation for averaging the products of local densities and for calculating thermodynamic functions at zero temperature, it is shown that the few-particle correlations imply a first-order phase transition and are reduced to the product of single-particle averages with a simultaneous increase in pressure, density, and quantum fluctuations. Under given conditions, dark matter exhibits rather the properties of an ideal gas with an effective temperature determined by quantum fluctuations. Characteristics of oscillations between bound and unbound states of three particles are estimated within a simple random walk approach to qualitatively models the instability of particle complexes. On the other hand, the density-dependent conditions for the formation of composites are analyzed using chemical kinetics without specifying the bonds formed. The obtain results can be extended to the models of multicomponent dark matter consisting of composites formed by particles with a large scattering length.

In the case of weakly interacting structureless bosons with or without specification of their sort, the description of the macroscopic properties of the dark matter halo of galaxies is usually based on the Gross-Pitaevskii equation under the imposed conditions of diluteness and a large healing length. Indeed, such models for the halo of dwarf galaxies give a good qualitative description with a very short (pairwise) scattering length of ultralight bosons [5,20,24]. It turns out that attempts to achieve better agreement between theoretical estimates and astronomical observables by taking into account additional terms of few-particle interactions in the Gross-Pitaevskii equation result in the existence of various thermodynamic phases of dark matter [5,23,24], as well as nontrivial correlations associated with few-particle complexes (composites).
On the other hand, there are experimental data on the formation of Bose condensate droplets [30,31,32] from a gas of alkali atoms with a long scattering length, which is accompanied by the processes of particle recombination, three-particle losses, and the formation of molecules that are technically suppressed in the laboratory. This regime, associated with the unitary limit of infinite scattering length, implies a strong interaction between particles leading to the structural transformations. Much attention is paid to its study in the literature [33,34,35] in connection with experimental exploration of the Bose condensate. Theoretically, similar processes like noted above cannot in principal be ruled out in dark matter, which could eventually be a (multicomponent) mixture of different composites with a resulting short scattering length, as required by current observations. In addition, they could take place under the extreme conditions in the early Universe [36]. Thus, the assumptions proposed give rise to many implications, and we would like to dwell on a few. It is worth to note that dark matter, made of several constituents, is already being actively studied [37,38].
Summing up, we outline the tasks of this paper, which is devoted to the analysis of density-dependent conditions for the formation of composites using rather simple models.
In Section 2, considering a model based on the Gross-Pitaevskii equation that takes into account two-and three-particle interactions, we try to obtain the admissibility conditions for the occurrence of particle complexes by estimating few-particle correlators over a relatively wide range of density, pressure, and quantum fluctuations at zero temperature, i.e., quantities that affect the transfer momentum in an inhomogeneous one-component system of bosons. In Section 3, we describe the time-oscillatory instability of two-and three-particle states that arises in a dense medium using a simple random walk model [39]. Various scenarios for the formation of composites from three particles are considered in Section 4 within the framework of chemical kinetics [40] without clarifying the nature of the bonds. We operate only with the scattering length and the density of the components as the most universal characteristics, in order to extract the necessary conditions for the existence of those. We end with the Discussion where we present our conclusions on the obtained results.
2 Self-Gravitating Bose-Einstein Condensate: The Role of Few-Particle Interactions We start with an analysis of the equilibrium properties of halo dark matter, described by the Gross-Pitaevskii equation, where both two-and three-particle repulsive interactions are taken into account. As established in [24], such a model reveals two-phase structure, represented by dilute and denser (liquidlike) states, and the first-order phase transition. It is intended to describe the resulting state of dark matter after a long-term evolution under certain conditions, although only some of the physical options provided by this model can be justified. After many attempts of description, with their special features, made in the literature [20,5,23], here we can expect and generally conclude that the resulting constituents of dark matter possess a vanishingly small (pair) s-scattering length a. This does not contradict the development of models describing particles with a longer scattering length, because their composites (with a short scattering length a) are allowed to exist in this "final" stage.
So, we first refer to our exploration of the model elaborated in [24], where the Bose-Einstein condensate (BEC) is described by a real function ψ(r) of the radial variable r = |r| in a ball B = {r ∈ R 3 | |r| ≤ R}. The main equation of our model is generated by the functional Γ, depended on a chemical potential µ, and is supplemented by the Poisson equation for gravitational field 4πGmV (r): Let two-and three-particle interactions be characterized by the quantities: where a > 0 is the s-scattering length of pairwise interaction; λ is the Compton wave length of the particle with mass m; g is a some dimensionless parameter introduced for convenience; ∆ r and ∆ −1 r are the radial part of the Laplace operator and its inverse.
Note that the relativistic nature of three-particle interaction, because of dependence of λ and U 3 on speed of light c, is pointed out in [5].
Applicability of the local form of interaction requires that where n 0 is the maximal particle density in the system. Vanishing of variation δΓ with respect to δψ gives the stationary Gross-Pitaevskii equation with the gravitational and three-particle interactions: where R is defined as the (first) zero of ψ(r) and determines a core boundary of BEC. The further description of the system is conveniently carried out in terms of dimensionless wave function χ and coordinate ξ: where n 0 and r 0 set the scale of the particle density and radius. A normalization of the dimensionless wave function χ(ξ) defines the mean particle density n in the total volume 4πR 3 /3: Empirically, the value of σ is less than unity [23,24]. The wave function χ(ξ) determines the local particle density η(ξ) = χ 2 (ξ) and is found from the field equation where dimensionless parameters of the model are The effective chemical potential ν plays the role of free parameter absorbing the gravitational potential in the center (at ξ = 0): and it leads to the following expression for Φ(ξ): As usual [27], we associate r 0 with the healing length: what leads immediately to Q = 1 in (9). Under the diluteness conditions (4) which guarantee applicability of the given formalism, the following inequality holds: Let us relate the characteristics of dark particles with the parameters of our model, using the typical mass density value ρ 0 = 10 −20 kg/m 3 , that gives Then, substituting it and (12) into (9), we obtain the relations It is required that mc 2 ≤ 1 eV accordingly to [19]. For definiteness, let us fix g = 4.73 × 10 −83 and Q = 1, A = π 2 9.87, B = 20, to obtain that a = 2.23 × 10 −78 m, m = 10 −22 eV/c 2 , r 0 = 0.818 kpc. Of course, we admit a more suitable parametrization for concrete situations. Besides, the core size R = r 0 ξ B and the mass density ρ = ρ 0 σ should be found by solving the field equations for given restrictions. However, here we focus on analyzing the physical effects without applying the model to specific objects.
Equations (8), (11) are numerically integrated under the following initial conditions: χ (0) = 0 and χ (0) < 0. As it was argued in [24], the finite initial value χ(0) = z lies in the range z 1 < z < z 2 , where and it should be minimal positive solution of equation: The absence of such a solution z at given (A, Q, B, ν) means that χ(ξ) = 0 everywhere. It happens at ν < ν min where ν min (A, Q, B) is found numerically. Studying the statistical properties of matter, one uses the part of oscillating wave function χ(ξ) in the range ξ ∈ [0; ξ B ], limited by its first zero χ(ξ B ) = 0. The value of ξ B determines the total radius R = r 0 ξ B of the system and depends on the set of parameters (A, Q, B, ν). Let us write down the internal energy density E and the local pressure P in dimensionless units: The statistical mean of any local characteristic is defined by the integral: In general, there are two limiting regimes of the model (8), (11): i) the Thomas-Fermi (TF) approximation, when the terms of quantum fluctuations like (∂ ξ χ) 2 and ∆ ξ χ are neglected; ii) the regime of the standing spherical wave, when the chemical potential ν is large enough and/or the interaction terms are negligibly small in (8).
Since our model generalizes the well-studied model with pair interaction (B = 0) in the Thomas-Fermi approximation [20], it is useful to indicate its basic equations and the values of the parameters. The set of this model equations contains Equation (11) and Solution to these equations is given (see [20,22]) by where j 0 (z) = sin z/z is the spherical Bessel function. In fact, the central density η c is regarded as a free parameter of the model instead of ν.
Limiting the system radius by the value r = r 0 (or ξ B = 1), one derives that κ = π. Hence, for arbitrary ν > 0, the parameters of such a model can be fixed as Q = 1 and A = π 2 . This confirms our choice of the value of A in (18).
Note that close parameter values, Q ∼ 1.36 and A = 10, were used in [24]. However, taking into account the quantum fluctuations ∆ ξ χ and the three-particle interaction regulated by parameter B, the system radius ξ B > 1 and the central density η c = χ 2 (0) reveal nonlinear dependence on ν (see Figure 1 below and Figure 3 in [23]).
Let us note the relation r 0 (a, m): which is discussed in detail in [20] to outline physically reasonable values of a and m by describing the dark matter of dwarf galaxies.
In contrast with the Thomas-Fermi approximation, equation (8) admits the limiting regime of the standing wave, implied by the equation: where w is a correction due to small interaction, which is not derived here. This equation has a simple analytic solution: Obvious merging of the approximate function ξ B (ν) and the numerically found one at large ν is shown in Figure 1 for the fixed values of Q and A.
Besides, the form of wave function (27) guarantees that η n (ξ) = η 2 c c n , where η c ≡ χ 2 0 , and the numeric coefficients c n are independent of ν ν min . Although the exact situation is only qualitatively the same, as we shall see in Figure 3 below, this effect excludes a possibility of the composites formation because of destruction of many-particle states at large both ν and particle momentum k = (∂ ξ χ) 2 /σ. As it is easily seen, the averages of η n with similar properties are obtained for the solution (24) in the Thomas-Fermi approximation.
Referring to [24], we especially note that the successful attempts of reproducing rotational curves of dwarf galaxies require the use of ν ν min , what can be associated with a ground state of the model (8), (11).
Moreover, the two regimes of "small" and "large" ν are naturally separated by condition µ = 0 for the physical chemical potential or its dimensionless value u ≡ ν/2 − Aυ = 0. In other words, the regime of large ν corresponds to u > 0, while the bound states are described at u < 0. As we have seen in [24,23], the point u = 0 determines parameters of the first-order phase transition in the system with three-particle interaction.
After matching the parameters of our model and the reference one, we would also like to compare some dependencies. We start by examining the mean particle density σ as a function of the central density η c , at large values of which the transformation of particles can occur.
As already shown in [20] for model (23)-(24) with ξ B = 1, the relation between σ and η c is given by Figure 2: Mean particle density σ (a) and system radius ξ B (b) as functions of central density η c in two models. In both panels, the orange line corresponds to the model with pair interaction in the Thomas-Fermi approximation. The red and blue branches of the single curve correspond to the model that takes into account both quantum fluctuations and three-particle interaction.
A completely different relation is expected for our model (8). The dependencies of interest are shown in Figure 2a.
The first difference between the model with pair interaction in the Thomas-Fermi approximation (23)-(24) and our model which takes into account both quantum fluctuations and three-particle interaction is the presence of the minimal admissible central density in the latter, which for the given model parameters is equal to η * c 0.1294. At η c ≥ η * c , there are two branches of the function σ(η c ) for the model (8), which are colored in blue and red in Figure 2a. The behavior of the blue branch is similar to the orange line described by (28), while the red curve looks abnormal due to three-particle repulsion. To explain the discrepancy between the mean density curves colored in blue and orange as well as the red branch, let us turn to the dependence of the size ξ B on η c .
We see in Figure 2b that an increase in the system size ξ B along red branch is accompanied by a decrease in the corresponding mean density in Figure 2a. That is, the red sectors in Figures 2a and 2b describe quite dilute matter.
Since the blue curve in Figure 2b describes a system with a larger radius ξ B than the orange line for the reference model, this is consistent with the lower density in blue compared to the mean density in orange in Figure 2a. As η c grows, the blue curve ξ B (η c ) tends to ξ B = 1. Therefore, the blue curves do show a denser state of matter.
On the other hand, the existence of composites in this dense matter becomes problematic due to the mutual influence leading to their dissociation into particles. Assuming that the interactions still guarantee the formation of composites for a short time, which is similar to oscillations, we will consider such a process in the next Section using the model of quantum random walk.
To examine the behavior of these quantities, found with the same parameters (A, Q, B) as before, but varying ν, we use the mean internal pressure P (see [24]): The quotients C 1,2,3 (P ), presented in Figure 3, reveal the splitting and equivalence: described by almost horizontal plateau at relatively high pressure P > 1.8 × 10 −3 . (To obtain the value of pressure and energy density in physical units, the dimensionless quantity should be multiplied by ε 0 = 2 n 0 /(mr 2 0 ), see (9).) This means that two and three-particle configurations behave as two and three isolated particles, respectively. The reasons of such a behavior are explained above by considering the standing wave solution (27). The transition to this regime is, generally speaking, a first-order phase transition described in [23,24]. It is natural to associate the instability of many-particle configurations (composites) with quantum fluctuations in the system. Indeed, the quantity τ = (∂ ξ χ) 2 undergoes a jump in Figure 4a when passing to the pressure interval of interest. Moreover, at P > 1.8×10 −3 , the almost linear proportionality between P and τ resembles the relation between pressure and temperature of an ideal gas.
In summary, we can conclude that the high pressure regime does not provide the formation of composites in BEC, although the particle density increases with the pressure, as shown in Figure 4b. Besides, the inequality nr 3 0 1 is fulfilled in BEC due to large healing length r 0 , which determines the BEC core radius [27]. Despite the strong correlations, characterized by r 0 , the diluteness condition na 3 < 1 persists.

Quantum Random Walk Based Treatment
Based on the predictions of the previous Section, we would like to consider the instability of two-and three-particle states in dense matter within the framework of a simple probabilistic model. In the case of an elementary system of three particles (monomers, denoted M), which are confined within the radius of interactions leading to the formation of a dimer (D) and a trimer (T), there is a possibility of oscillations, namely, a continuous transition between three states: representing the unbound state and two bound ones. In the absence of external factors, the reversibility of microscopic processes is supported by the equality of the probabilities of forming and decaying of the composites. Of course, the known lifetimes of each (metastable) state could point out the dominant configuration of the particles. However, having permanently interacting particles, we will neglect the probability of their adhesion. Nevertheless, it seems possible to estimate the probabilities of states over the period of oscillations. Thus, we consider a toy model with three identical monomers, which is described by the wave function Ψ(t): that evolves in time t according to the Schrödinger equation with initial state (at t = 0) assumed to be Ψ 0 = ψ 1 . Then we define the fractions Ω i (t) = |a i (t)| 2 and analyze their means Ω i over a period of time.
Evolution of the initial state, is given by the Schrödinger equation with HamiltonianĤ =Ĥ 0 +Û , whereĤ 0 = diag(ε 1 , ε 2 , ε 3 ): Actually, the matrixĤ with constant elements determines quantum random walk (a Markov process). Analytically, solution of (35) for the wave function can be immediately represented (and found) in the form: Neglecting spatial evolution of localized particles and time delay in each state, we omitĤ 0 and remain with the matrix of admissible transitions: where the quantities U ±i play a role of transition amplitudes, and we put U −i = U i ∈ R to makeÛ self-adjoint with real elements. The conserved energy of the system is Introducing auxiliary quantities: three eigenvalues ofÛ , which are roots of equation are found to be where General solution to equation (35) is presented in the form: where u i are the non-normalized eigenvectors of matrixÛ ; C i are the constants, determined by initial condition |Ψ(0) = |Ψ 0 from (34): We omit the study of time dependence in three-parameter case (43) and restrict ourselves to a particular case, setting U 3 = 0 in (42), when We easily obtain the spectrum: Physically, these quantities determine the energy levels in the system and indicate the presence of the unbound state of 3M (three monomers) with λ 1 , the bound one of T with λ 3 , and a transient one with λ 2 , associated with M+D.
It is useful to parametrize U 1 and U 2 by s and T as To obtain this, we can also use in (43) that Note that configuration with Ω M = 0 and Ω T > Ω D is realized at the time moment (during the first period) when s tends to π/4 from below.
Averaging Ω α (t) over the period T , we arrive at Comparison of these fractions for different s is made in Figure 5. We see a slight dominance of Ω T in the range s c < s < π/4, while Ω D prevails for s < s c , where s c = arcsin √ 3/3 3.13π/16, and Ω M (s c ) = Ω D (s c ) = Ω T (s c ) = 1/3.
In summary, we note that in a three-particle system with permanent oscillations between quantum states there exists the possibility of the composites dominance in a long-time picture. This effect is due to the tuning of interaction parameters. In general, the probability of detecting monomers (of dark matter) increases at higher density and pressure in the model with few-particle interactions.

Chemical Kinetics Consideration
For the sake of simplicity, let us assume that the elementary constituents of dark matter are identical spinless particles (monomers, or M). Due to an interaction between them, the nature of which is not specified herein, we assume the ability of formation of dimers (D) and trimers (T) in the ground state. Therefore, we restrict ourselves to processes involving no more than three monomers, taking into account also those contained in dimers and trimers. Attempting to describe the formation/decay of dimers and trimers, as well as their fractions, we clarify some physical conditions. Account of the gravitational interaction, resulting in an inhomogeneous distribution of particles, implies that the statistical characteristics of the matter should, in principle, depend on space. We assume that the particle density is highest in a core and tends to zero at halo periphery. However, in regions of space that are small in astronomical scale, but large enough to apply the statistical approach, the gravitational potential (together with the chemical one) can be considered as constant one. This allows us to explore the properties of quasi-homogeneous subsystems. However, simplifying the problem by introducing regions of homogeneity with different particle densities, the dominance of certain processes does depend on the ratio of the interparticle distance to the radius of interaction (or scattering lengths).
The particle density in the core can be estimated as which is determined by the individual bosonic monomer mass m 1 eV and the typical mass density ρ 0 . Let us consider possible reactions occurring simultaneously in each subsystem: where the positive constants k i and k −i (i = 1, 2, 3) are related to the rates of the corresponding reactions and have the dimension of inverse time.
The time evolution of concentrations C α (the numbers of particles of sort α = M, D, T in a fixed volume) can be described accordingly to the rules of chemical kinetics [40]: with the initial conditions at t = 0: D,T = 0; for V 2M , V 3M , and V MD , see next proposition. The above set of equations for C α conserves in time the total number of monomers in the subsystem of a fixed volume, and implies that the elementary act depends on constants k ±i and interaction volumes V s a 3 s , where a 2M , a 3M and a MD are the scattering lengths (for s-wave channel, in instance). It is natural to require conditions for the successful performing of three forward reactions: They indicate the presence of at least one particle in each interaction volume V s . Failure to meet one of these conditions may lead to the exclusion of the corresponding process from consideration. Associating C α with the local particle density far from the core (of galactic halo) center, in order to satisfy the imposed conditions (55), we require that n 0 a 3 s 1 for the scattering lengths a s in the model. Strictly speaking, the validity of condition na 3 s > 1 depends on both n and a s . However, as shown above, an increase in monomer density n (due to squeezing, for example) correlates with the destruction of composites. Therefore, there remains the only possibility of considering interactions with longer scattering lengths a s . Next, we evaluate the contribution of dimers and trimers to dark matter.
To simplify further notations, we will also use It is convenient to describe the subsystem evolution in terms of the fractions: , where c M = 1, c D = 2, c T = 3.
Let us analyze the equilibrium fractions Ω eq α when dC α /dt = 0. For equilibrium concentrations (dropping the label 'eq'), we find C D and C T as functions of C M : where the elements of matrix A α are formed from K ±i : .
For a subsystem in a fixed (unit) volume, the initial concentration C M is related to the equilibrium one C M : That means, the initial state C M can be determined by using the C M . Note the possibility to generalize our considerations by specifying the concentration C (q) α of composites of sort α in certain quantum state q at a given temperature T , that is, the occupation of energy level E (q) α with the degeneracy g q . For the processes that do not violate the equilibrium distribution of particles over states, one writes [40]: where Z α is the partition function over states, and C α is the total concentration. From equation (53), with account of Boltzmann factors p α (q), one can deduce the dependence of the reaction rates k ±i on non-zero temperature. At finite K ±i , there are two limiting regimes: i) C M → 0 (low concentration of monomers) admits approximation: C (0) and , Ω eq T = 3 ii) asymptotic expressions at large C M yield: , where C 0 denotes a real solution among three roots C m (m = 0, ±1) of a cubic equation: To provide Ω eq T ≥ 0.5, we demand 3C 2 M ≥ κ 1 κ 2 + 2κ 2 C M [see (65)], or, equivalently, Therefore, C M a 3 MD > k −2 /k 2 , where k −2 /k 2 ≤ 1 to permit the synthesis of T. Thus, condition (74) is not stronger than (55).
Summarizing, the applied analysis enables the leading role of composites under the derived conditions (55) and the definite values of monomer density C  (69) and (74)]. However, the mechanism of the formation of composites needs detailed study.

Discussion
In this paper, we have considered three models that deal with the problem of formation and stability of the two-and three-particle composites in dark matter.
Starting with the Bose-condensate dark matter model with respective self-interaction on the base of the Gross-Pitaevskii equation extended by an additional three-particle interaction, we find that the fewparticle complexes can be formed at relatively low density and pressure, when quantum fluctuations and the transfer momentum are small enough and cannot destroy the condensate and the desired structures. Such a conclusion is deduced on the base of calculating the averages of the powers of the local density η and thermodynamic functions at zero temperature. However, increasing the pressure (and density), the multiparticle averages are reduced to the product of single-particle averages which determine the dependence on thermodynamic parameters. This enables to relate the thermodynamic functions by relations for an ideal gas, where the mean value of quantum fluctuations is regarded as the effective temperature of the inhomogeneous system.
Note that a similar splitting of many-particle correlators into single-particle functions is observed in the model with pairwise interaction in the Thomas-Fermi approximation. This emphasizes the importance of the model with three-particle interaction along with explicit accounting of quantum fluctuations.
We associate the revealed changes in properties with the first-order phase transition described in our previous work [24]. Besides, as found in Ref. [24], a description of rotational curves using this model is realized in the vicinity of the smallest value of the chemical potential (the ground state) and at a short pairwise scattering length a, the magnitude of which we reproduced above. Then, we can compare our expression (21) for the internal energy with the well-known one, where the term of threeparticle interaction arises due to quantum corrections in a system with pairwise interaction and the use of renormalization group calculations [33,41,42,43]: Here n is the particle density of homogeneous system of bosons, a is their pairwise s-scattering length, d is a numeric parameter. Note that the internal energy E is proportional to inverse of the mass of (dark matte) particle. However, we also admit that the origin of the three-body interaction can be explained by relativistic effects [5,6] or by the very nature of dark matter particles.
It is interesting to note the attempts to generalize the expression (75) in the case of the unitary regime with an infinitely long scattering length a [44,45]. This is motivated, in particular, by the study of Bose condensates of alkali atoms, the strong interaction between which can cause the molecules synthesis, recombination, and three-particle losses. Hypothetically, similar processes can result in multicomponent structure of dark matter , if the resulting composites are stable enough and have a short scattering length. Although, due to the greater mass, the Compton wavelength of the composite would already become shorter than that of individual constituents.
Within the framework of an efficient model of chemical kinetics in Section 4, we show that a long scattering length, but not a high density, is necessary for the formation of few-particle complexes. Even without knowledge of the nature of the bonds, it is clear that the appearance of dimers and trimers therein cannot be compared with the synthesis of molecules from atoms that are complicated by themselves. Moreover, there is no evidence that the processes under consideration are similar to the production of baryons: quark-antiquark mesons and three-quark nucleons or hyperons.
We assume that the Efimov effect could be appropriate mechanism of forming the trimers from dark bosons. The existence of such an effect is successfully confirmed by experiments with cold-atom systems [46] and is investigated [35]. Comparison of manifestations of this effect in the two different situations shows a number of similar features.
To simulate unstable composites consisting of no more than three particles, we consider a random walk model in Section 3. By adjusting the coupling constants, it is able to obtain configurations with the dominance of composites over long time intervals (long oscillation period). This model is quite useful to visualize the decay of correlators in the Bose condensate with a three-body interaction.
We expect that detailed quantum analysis may provide results applicable to the description of dark matter, which will be published elsewhere.