Bose–Einstein Condensate Dark Matter That Involves Composites

Abstract

Improving the Bose–Einstein condensate model of dark matter through the repulsive three-particle interaction to better reproduce observables such as rotation curves reveals both different thermodynamic phases and few-particle correlations. 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 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 model 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 obtained results can be extended to the models of multicomponent dark matter consisting of composites formed by particles with a large scattering length.

1. Introduction

The current observable restrictions on the properties of dark matter still admit the construction of theoretical models with different constituents [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]. Adopting the concept of dark matter consisting of light bosons in a condensate state [15,16,17,18,19,20,21,22,23,24,25,26,27], we might also borrow the known properties of atomic Bose–Einstein condensates in gaseous or liquid form [28,29,30] reproduced in the laboratory.

In the case of weakly interacting structureless bosons with or without specifications 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,21,25]. 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,24,25] 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 [31,32,33] 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 structural transformations. Much attention is paid to its study in the literature [34,35,36] in connection with experimental exploration of the Bose condensate. Theoretically, processes similar to the ones 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 [37]. Thus, the assumptions proposed give rise to many implications, and we would like to dwell on a few. It is worth noting that dark matter made of several constituents is already being actively studied [38,39].

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 [40]. Various scenarios for the formation of composites from three particles are considered in Section 4 within the framework of chemical kinetics [41] 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 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 [25], such a model reveals a 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 at describing their special features made in the literature [5,21,24], 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.

We first refer to our exploration of the model elaborated in [25], where the Bose–Einstein condensate (BEC) is described by a real function $\psi \left(r\right)$ of the radial variable $r=\left|\mathbf{r}\right|$ in a ball $B=\{\mathbf{r}\in {\mathbb{R}}^3|\left|\mathbf{r}\right|\le R\}$. The main equation of our model is generated by the functional $\mathsf{\Gamma }$, dependent on a chemical potential $\mu$, and is supplemented by the Poisson equation for gravitational field $4\pi GmV\left(r\right)$:

$$\begin{array}{cc}\hfill \mathsf{\Gamma }=& 4\pi {\int }_0^R\left[\frac{{\hslash }^2}{2m}{\left({\partial }_r\psi \left(r\right)\right)}^2-\mu {\psi }^2\left(r\right)+4\pi G{m}^2{\psi }^2\left(r\right)V\left(r\right)\right.\hfill \\ & \left.+\frac{U_2}{2}{\psi }^4\left(r\right)+\frac{U_3}{3}{\psi }^6\left(r\right)\right]r^2\mathrm{d}r,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill {\Delta }_{r}V\left(r\right)& =& {\psi }^2\left(r\right).\hfill \end{array}$$

(2)

Let two- and three-particle interactions be characterized by the quantities:

$$U_2=\frac{4\pi {\hslash }^{2}a}{m},U_3=g{U}_2{\lambda }^3,\lambda =\frac{\hslash }{mc},$$

(3)

where $a>0$ is the s-scattering length of pairwise interaction; $\lambda$ is the Compton wave length of the particle with mass m; g is a dimensionless parameter introduced for convenience; ${\Delta }_r$ and ${\Delta }_r^{-1}$ are the radial part of the Laplace operator and its inverse.

Note that the relativistic nature of three-particle interaction because of dependence of $\lambda$ and $U_3$ on the speed of light c is pointed out in [5].

Applicability of the local form of interaction requires that:

$$n_0{a}^3<1,n_0{\lambda }^3<1,$$

(4)

where $n_0$ is the maximal particle density in the system.

The vanishing of variation $\delta \mathsf{\Gamma }$ with respect to $\delta \psi$ provides the stationary Gross–Pitaevskii equation with the gravitational and three-particle interactions:

$$\left[-\frac{{\hslash }^2}{2m}{\Delta }_r-\mu +4\pi G{m}^{2}V\left(r\right)\right]\psi +U_2{\psi }^3+U_3{\psi }^5=0,\psi \left(R\right)=0,$$

(5)

where R is defined as the (first) zero of $\psi \left(r\right)$ and determines a core boundary of BEC.

The further description of the system is conveniently carried out in terms of dimensionless wave function $\chi$ and coordinate $\xi$:

$$\psi \left(r\right)=\sqrt{n_0}\chi \left(\xi \right),r=r_0\xi ,R=r_0{\xi }_B,$$

(6)

where $n_0$ and $r_0$ set the scale of the particle density and radius.

A normalization of the dimensionless wave function $\chi \left(\xi \right)$ defines the mean particle density n in the total volume $4\pi R^3/3$:

$$\sigma \equiv \frac{3}{{\xi }_B^3}{\int }_0^{{\xi }_B}{\chi }^2\left(\xi \right){\xi }^2\mathrm{d}\xi =\frac{n}{n_0}.$$

(7)

Empirically, the value of $\sigma$ is less than unity [24,25].

The wave function $\chi \left(\xi \right)$ determines the local particle density $\eta \left(\xi \right)={\chi }^2\left(\xi \right)$ and is found from the field equation

$$\frac{1}{2}({\Delta }_{\xi }+\nu )\chi -A\chi \mathsf{\Phi }-Q{\chi }^3-B{\chi }^5=0,$$

(8)

where dimensionless parameters of the model are:

$$Q=4\pi r_0^{2}a{n}_0,A=Q\frac{G{m}^3{r}_0^2}{{\hslash }^{2}a},B=gQ{\lambda }^3{n}_0,{\epsilon }_0=\frac{{\hslash }^2}{m{r}_0^2}{n}_0.$$

(9)

The effective chemical potential $\nu$ plays the role of free parameter absorbing the gravitational potential in the center (at $\xi =0$):

$$\upsilon =-A{\int }_0^{{\xi }_B}\eta \left(\xi \right)\xi \mathrm{d}\xi ,$$

(10)

and it leads to the following expression for $\mathsf{\Phi }\left(\xi \right)$:

$$\mathsf{\Phi }\left(\xi \right)=-\frac{1}{\xi }{\int }_0^{\xi }\eta \left(s\right)s^2\mathrm{d}s+{\int }_0^{\xi }\eta \left(s\right)s\mathrm{d}s,{\Delta }_{\xi }\mathsf{\Phi }\left(\xi \right)=\eta \left(\xi \right).$$

(11)

As usual [28], we associate $r_0$ with the healing length:

$$r_0=\frac{1}{\sqrt{4\pi a{n}_0}},$$

(12)

that leads immediately to $Q=1$ in (9). Under the diluteness conditions (4) that guarantee applicability of the given formalism, the following inequality holds:

$$r_0^3{n}_0\gg 1.$$

(13)

Let us relate the characteristics of dark particles with the parameters of our model, using the typical mass density value ${\rho }_0={10}^{-20}\mathrm{kg}/{\mathrm{m}}^3$, that gives:

$$n_0=\frac{{\rho }_0}{m}.$$

(14)

Then, substituting it and (12) into (9), we obtain the relations:

$$\begin{array}{ccc}& & m={\left(\frac{g}{B}\right)}^{1/4}{m}_0,m_0\equiv {\left(\frac{{\hslash }^3}{c^3}{\rho }_0\right)}^{1/4}=0.081\mathrm{eV}/c^2,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & a={\left(\frac{g}{4\pi AB}\right)}^{1/2}{\ell }_{\mathrm{P}},{\ell }_{\mathrm{P}}\equiv \sqrt{\frac{G\hslash }{c^3}}=1.62\times {10}^{-35}\mathrm{m},\hfill \end{array}$$

(16)

$$\begin{array}{ccc}& & r_0={\left(\frac{B}{g}\right)}^{1/8}{\left(\frac{A}{4\pi }\right)}^{1/4}L,L\equiv \sqrt{\frac{m_0}{{\rho }_0{\ell }_{\mathrm{P}}}}=9.44\times {10}^8\mathrm{m}.\hfill \end{array}$$

(17)

It is required that $m{c}^2\le 1\mathrm{eV}$ according to [20]. For definiteness, let us fix $g=4.73\times {10}^{-83}$ and:

$$Q=1,A={\pi }^2\simeq 9.87,B=20,$$

(18)

to obtain that $a=2.23\times {10}^{-78}\mathrm{m}$, $m={10}^{-22}\mathrm{eV}/c^2$, $r_0=0.818\mathrm{kpc}$.

Of course, we admit a more suitable parametrization for concrete situations. In addition, the core size $R=r_0{\xi }_B$ and the mass density $\rho ={\rho }_0\sigma$ 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) and (11) are numerically integrated under the following initial conditions: ${\chi }^{\prime }\left(0\right)=0$ and ${\chi }^{″}\left(0\right)<0$. As argued in [25], the finite initial value $\chi \left(0\right)=z$ lies in the range $z_1<z<z_2$, where:

$$z_1={\left[\sqrt{{\left(\frac{3Q}{10B}\right)}^2+\frac{\nu }{10B}}-\frac{3Q}{10B}\right]}^{1/2},z_2={\left[\sqrt{{\left(\frac{Q}{2B}\right)}^2+\frac{\nu }{2B}}-\frac{Q}{2B}\right]}^{1/2},$$

(19)

and it should be the minimal positive solution of the equation:

$$\frac{A}{2}{z}^2-\left(5B{z}^4+3Q{z}^2-\frac{\nu }{2}\right)\left(\frac{\nu }{2}-Q{z}^2-B{z}^4\right)=0.$$

(20)

Absence of such a solution z at given $(A,Q,B,\nu )$ means that $\chi \left(\xi \right)=0$ everywhere. It happens at $\nu <{\nu }_{\min }$ where ${\nu }_{\min }(A,Q,B)$ is found numerically.

Studying the statistical properties of matter, one uses the part of oscillating wave function $\chi \left(\xi \right)$ in the range $\xi \in [0;{\xi }_B]$, limited by its first zero $\chi \left({\xi }_B\right)=0$. The value of ${\xi }_B$ determines the total radius $R=r_0{\xi }_B$ of the system and depends on the set of parameters $(A,Q,B,\nu )$. Let us write down the internal energy density $\mathcal{E}$ and the local pressure $\mathcal{P}$ in dimensionless units:

$$\mathcal{E}\left(\xi \right)=\frac{Q}{2}{\eta }^2\left(\xi \right)+\frac{B}{3}{\eta }^3\left(\xi \right),\mathcal{P}\left(\xi \right)=\frac{Q}{2}{\eta }^2\left(\xi \right)+\frac{2}{3}B{\eta }^3\left(\xi \right).$$

(21)

The statistical mean of any local characteristic is defined by the integral:

$$\langle (\dots )\rangle =\frac{3}{{\xi }_B^3}{\int }_0^{{\xi }_B}(\dots ){\xi }^2\mathrm{d}\xi .$$

(22)

In general, there are two limiting regimes of the model defined by (8) and (11): (i) the Thomas–Fermi (TF) approximation when the terms of quantum fluctuations such as ${\left({\partial }_{\xi }\chi \right)}^2$ and ${\Delta }_{\xi }\chi$ are neglected; (ii) the regime of the standing spherical wave when the chemical potential $\nu$ 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 [21], it is useful to indicate its basic equations and the values of the parameters. The set of this model equation contains Equation (11) and:

$$\frac{\nu }{2}-A\mathsf{\Phi }-Q\eta =0.$$

(23)

Solution to these equations is given (see [21,23]) by:

$$\eta \left(\xi \right)={\eta }_c{j}_0\left(\kappa \xi \right),\kappa =\sqrt{\frac{A}{Q}},{\eta }_c=\frac{\nu }{2Q},$$

(24)

where $j_0\left(z\right)=sinz/z$ is the spherical Bessel function. In fact, the central density ${\eta }_c$ is regarded as a free parameter of the model instead of $\nu$.

Limiting the system radius by the value $r=r_0$ (or ${\xi }_B=1$), one derives that $\kappa =\pi$. Hence, for arbitrary $\nu >0$, the parameters of such a model can be fixed as $Q=1$ and $A={\pi }^2$. This confirms our choice of the value of A in (18).

Note that close parameter values, $Q\sim 1.36$ and $A=10$, were used in [25]. However, taking into account the quantum fluctuations ${\Delta }_{\xi }\chi$ and the three-particle interaction regulated by parameter B, the system radius ${\xi }_B>1$ and the central density ${\eta }_c={\chi }^2\left(0\right)$ reveal nonlinear dependence on $\nu$ (see Figure 1 below and Figure 3 in [24]).

Let us note the relation $r_0(a,m)$:

$$r_0=\pi \sqrt{\frac{{\hslash }^{2}a}{G{m}^3}},$$

(25)

which is discussed in detail in [21] 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:

$$({\Delta }_{\xi }+k^2)\chi \left(\xi \right)=0,k^2=\nu -w>0,$$

(26)

where w is a correction due to small interaction, which is not derived here.

This equation has a simple analytic solution:

$$\chi \left(\xi \right)={\chi }_0{j}_0\left(\pi \frac{\xi }{{\xi }_B}\right),{\xi }_B=\frac{\pi }{\sqrt{\nu -w}}.$$

(27)

Merging of the approximate function ${\xi }_B\left(\nu \right)$ and the numerically found one at large $\nu$ is shown in Figure 1 for the fixed values of Q and A.

In addition, the form of wave function (27) guarantees that $\langle {\eta }^n\left(\xi \right)\rangle ={\eta }_c^2{c}_n$, where ${\eta }_c\equiv {\chi }_0^2$, and the numerical coefficients $c_n$ are independent of $\nu \gg {\nu }_{\min }$. Although the exact situation is only qualitatively the same as we shall see below, this effect excludes a possibility of the composites formation because of the destruction of many-particle states at both large $\nu$ and particle momentum $k=\sqrt{\langle {({\partial }_{\xi }\chi )}^2\rangle /\sigma }$. As it is easily seen, the averages of ${\eta }^n$ with similar properties are obtained for the solution (24) in the Thomas–Fermi approximation.

Referring to [25], we especially note that the successful attempts at reproducing rotational curves of dwarf galaxies require the use of $\nu \simeq {\nu }_{\min }$ that can be associated with a ground state of the model (8) and (11).

Moreover, the two regimes of “small” and ”large” $\nu$ are naturally separated by condition $\mu =0$ for the physical chemical potential or its dimensionless value $u\equiv \nu /2-A\upsilon =0$. In other words, the regime of large $\nu$ corresponds to $u>0$, while the bound states are described at $u<0$. As we have seen in [24,25], the point $u=0$ determines the 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 $\sigma$ as a function of the central density ${\eta }_c$, at large values of which the transformation of particles can occur.

As already shown in [21] for models (23) and (24) with ${\xi }_B=1$, the relation between $\sigma$ and ${\eta }_c$ is given by

$$\sigma =\frac{3}{{\pi }^2}{\eta }_c.$$

(28)

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) and (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 ${\eta }_c^{*}\simeq 0.1294$.

At ${\eta }_c\ge {\eta }_c^{*}$, there are two branches of the function $\sigma \left({\eta }_c\right)$ 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 ${\xi }_B$ on ${\eta }_c$.

We see in Figure 2b that an increase in the system size ${\xi }_B$ along the red branch is accompanied by a decrease in the corresponding mean density in Figure 2a. That is, the red sectors in Figure 2a,b describe quite dilute matter.

Since the blue curve in Figure 2b describes a system with a larger radius ${\xi }_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 ${\eta }_c$ grows, the blue curve ${\xi }_B\left({\eta }_c\right)$ tends to ${\xi }_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 test our assumptions, we introduce the ratios of the dimensionless means:

$$C_1=\frac{\langle {\eta }^2\rangle }{{\langle \eta \rangle }^2},C_2=\frac{\langle {\eta }^3\rangle }{{\langle \eta \rangle }^3},C_3=\frac{\langle {\eta }^3\rangle }{{\langle \eta \rangle }^2\langle \eta \rangle },$$

(29)

where the particle density $\eta \left(\xi \right)={\chi }^2\left(\xi \right)$ is formed from the solution $\chi \left(\xi \right)$ of (8).

To examine the behavior of these quantities, found with the same parameters $(A,Q,B)$ as before, but varying $\nu$, we use the mean internal pressure P (see [25]):

$$P=\frac{Q}{2}\langle {\eta }^2\rangle +\frac{2}{3}B\langle {\eta }^3\rangle .$$

(30)

The quotients $C_{1,2,3}\left(P\right)$, presented in Figure 3, reveal the splitting and equivalence:

$$\langle {\eta }^2\rangle \sim {\langle \eta \rangle }^2,\langle {\eta }^3\rangle \sim {\langle \eta \rangle }^3,\langle {\eta }^3\rangle \sim \langle {\eta }^2\rangle \langle \eta \rangle ,$$

(31)

described by an almost horizontal plateau at relatively high pressure $P>1.8\times {10}^{-3}$. (To obtain the value of pressure and energy density in physical units, the dimensionless quantity should be multiplied by ${\epsilon }_0={\hslash }^2{n}_0/\left(m{r}_0^2\right)$, see (9)). This means that two- and three-particle configurations behave as two and three isolated particles, respectively. The reasons for 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 [24,25].

It is natural to associate the instability of many-particle configurations (composites) with quantum fluctuations in the system. Indeed, the quantity $\tau =\langle {({\partial }_{\xi }\chi )}^2\rangle$ undergoes a jump in Figure 4a when passing to the pressure interval of interest. Moreover, at $P>1.8\times {10}^{-3}$, the almost linear proportionality between P and $\tau$ 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. In addition, the inequality $n{r}_0^3\gg 1$ is fulfilled in BEC due to large healing length $r_0$, which determines the BEC core radius [28]. Despite the strong correlations characterized by $r_0$, the diluteness condition $n{a}^3<1$ persists.

3. 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:

$$\begin{array}{ccc}& {\psi }_1={|3\rangle }_{\mathrm{M}}{\otimes |0\rangle }_{\mathrm{D}}{\otimes |0\rangle }_{\mathrm{T}},{\psi }_2={|1\rangle }_{\mathrm{M}}{\otimes |1\rangle }_{\mathrm{D}}\otimes {|0\rangle }_{\mathrm{T}},& \\ & {\psi }_3={|0\rangle }_{\mathrm{M}}{\otimes |0\rangle }_{\mathrm{D}}\otimes {|1\rangle }_{\mathrm{T}},& \hfill \end{array}$$

(32)

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 $\mathsf{\Psi }\left(t\right)$:

$$\mathsf{\Psi }\left(t\right)=\sum _{i=1}^3{a}_i\left(t\right){\psi }_i,\sum _{i=1}^3{\left|a_i\left(t\right)\right|}^2=1,$$

(33)

that evolves in time t according to the Schrödinger equation with initial state (at $t=0$) assumed to be ${\mathsf{\Psi }}_0={\psi }_1$.

Then we define the fractions ${\mathsf{\Omega }}_i\left(t\right)={\left|a_i\left(t\right)\right|}^2$ and analyze their means ${\overline{\mathsf{\Omega }}}_i$ over a period of time.

The evolution of the initial state:

$$|{\mathsf{\Psi }}_0\rangle =\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)\stackrel{t}{\to }|\mathsf{\Psi }\left(t\right)\rangle =\left(\begin{array}{c}{a}_1\left(t\right)\\ a_2\left(t\right)\\ a_3\left(t\right)\end{array}\right),$$

(34)

is given by the Schrödinger equation with Hamiltonian $\widehat{H}={\widehat{H}}_0+\widehat{U}$, where ${\widehat{H}}_0=\mathrm{diag}({\epsilon }_1,{\epsilon }_2,{\epsilon }_3)$:

$$\mathrm{i}\hslash \frac{\mathrm{d}}{\mathrm{d}t}|\mathsf{\Psi }\rangle =\widehat{H}|\mathsf{\Psi }\rangle .$$

(35)

Actually, the matrix $\widehat{H}$ with constant elements determines quantum random walk (a Markov process). Analytically, the solution of (35) for the wave function can be immediately represented (and found) in the form:

$$|\mathsf{\Psi }\left(t\right)\rangle =\widehat{S}\left(t\right)|{\mathsf{\Psi }}_0\rangle ,\widehat{S}\left(t\right)=exp\left(-\frac{\mathrm{i}}{\hslash }\widehat{H}t\right).$$

(36)

Neglecting spatial evolution of localized particles and time delay in each state, we omit ${\widehat{H}}_0$ and remain with the matrix of admissible transitions:

$$\widehat{U}=\left(\begin{array}{ccc}0& U_{-1}& U_{-3}\\ U_1& 0& U_{-2}\\ U_3& U_2& 0\end{array}\right),$$

(37)

where the quantities $U_{\pm i}$ play a role of transition amplitudes, and we put $U_{-i}=U_i\in \mathbb{R}$ to make $\widehat{U}$ self-adjoint with real elements.

The conserved energy of the system is:

$$E\equiv \langle \mathsf{\Psi }\left(t\right)|\widehat{H}|\mathsf{\Psi }\left(t\right)\rangle =U_1{a}_2^{*}\left(t\right)a_1\left(t\right)+U_2{a}_3^{*}\left(t\right)a_2\left(t\right)+U_3{a}_1^{*}\left(t\right)a_3\left(t\right)+\mathrm{c}.\mathrm{c}.$$

(38)

Introducing auxiliary quantities:

$$U^{\left(3\right)}=U_1{U}_2{U}_3,U^{\left(2\right)}=\frac{1}{3}(U_1^2+U_2^2+U_3^2),$$

(39)

three eigenvalues of $\widehat{U}$, which are roots of equation:

$$0=det\left(\widehat{U}-\lambda \widehat{I}\right)=-{\lambda }^3+3\lambda U^{\left(2\right)}+2U^{\left(3\right)},$$

(40)

are found to be:

$$\begin{array}{ccc}& {\lambda }_1=A+B,{\lambda }_2={\mathrm{e}}^{-2\pi \mathrm{i}/3}A+{\mathrm{e}}^{2\pi \mathrm{i}/3}B,& \hfill \\ & {\lambda }_3={\mathrm{e}}^{2\pi \mathrm{i}/3}A+{\mathrm{e}}^{-2\pi \mathrm{i}/3}B,\end{array}$$

(41)

where:

$$A={\left[U^{\left(3\right)}+\sqrt{{\left(U^{\left(3\right)}\right)}^2-{\left(U^{\left(2\right)}\right)}^3}\right]}^{1/3},B=\frac{U^{\left(2\right)}}{A}.$$

(42)

The general solution to Equation (35) is presented in the form:

$$|\mathsf{\Psi }\left(t\right)\rangle =\sum _{i=1}^3{C}_i{\mathbf{u}}_{i}exp\left(-\mathrm{i}\frac{{\lambda }_{i}t}{\hslash }\right),{\mathbf{u}}_i=\left(\begin{array}{c}{\lambda }_i^2-U_2^2\\ {\lambda }_i{U}_1+U_2{U}_3\\ {\lambda }_i{U}_3+U_1{U}_2\end{array}\right),$$

(43)

where ${\mathbf{u}}_i$ are the non-normalized eigenvectors of matrix $\widehat{U}$; $C_i$ are the constants determined by initial condition $|\mathsf{\Psi }\left(0\right)\rangle =|{\mathsf{\Psi }}_0\rangle$ from (34):

$$C_1=\mathrm{i}\frac{\sqrt{3}}{9}\frac{{\lambda }_2-{\lambda }_3}{A^3-B^3},C_2=\mathrm{i}\frac{\sqrt{3}}{18}\frac{{\lambda }_3-{\lambda }_2-3{\lambda }_1}{A^3-B^3},C_3=\mathrm{i}\frac{\sqrt{3}}{18}\frac{{\lambda }_3-{\lambda }_2+3{\lambda }_1}{A^3-B^3}.$$

(44)

We omit the study of time dependence in a three-parameter case (43) and restrict ourselves to a particular case, setting $U_3=0$ in (42), when:

$$A=3^{-1/2}{\mathrm{e}}^{\mathrm{i}\pi /6}\sqrt{U_1^2+U_2^2}.$$

We easily obtain the spectrum:

$${\lambda }_1=\sqrt{U_1^2+U_2^2},{\lambda }_2=0,{\lambda }_3=-\sqrt{U_1^2+U_2^2}.$$

(45)

Physically, these quantities determine the energy levels in the system and indicate the presence of the unbound state of 3M (three monomers) with ${\lambda }_1$, the bound one of T with ${\lambda }_3$, and a transient one with ${\lambda }_2$, associated with M+D.

It is useful to parametrize $U_1$ and $U_2$ by s and T as:

$$U_1=Jcoss,U_2=Jsins,J=\hslash \omega =\frac{2\pi \hslash }{T}.$$

(46)

Parameter T defines the time scale of the problem, while the quantities $coss/T$ and $sins/T$ represent the rates of the processes.

Under the imposed conditions, the components of the matrix $\widehat{S}\left(t\right)$ are easily calculated and yield the instant fractions:

$$\begin{array}{ccc}& {\mathsf{\Omega }}_{\mathrm{M}}(t,s)={\left[{cos}^{2}scos\left(\omega t\right)+{sin}^{2}s\right]}^2,{\mathsf{\Omega }}_{\mathrm{D}}(t,s)={\left[cosssin\left(\omega t\right)\right]}^2,& \\ & {\mathsf{\Omega }}_{\mathrm{T}}(t,s)={\left[cosssins\left(cos\left(\omega t\right)-1\right)\right]}^2.& \hfill \end{array}$$

(47)

To obtain this, we can also use in (43) that:

$$C_1{\mathbf{u}}_1=\frac{1}{2}\left(\begin{array}{c}{cos}^{2}s\\ coss\\ cosssins\end{array}\right),C_2{\mathbf{u}}_2=\left(\begin{array}{c}{sin}^{2}s\\ 0\\ -cosssins\end{array}\right),C_3{\mathbf{u}}_3=\frac{1}{2}\left(\begin{array}{c}{cos}^{2}s\\ -coss\\ cosssins\end{array}\right).$$

(48)

Note that the configuration with ${\mathsf{\Omega }}_{\mathrm{M}}=0$ and ${\mathsf{\Omega }}_{\mathrm{T}}>{\mathsf{\Omega }}_{\mathrm{D}}$ is realized at the time moment (during the first period):

$$t_1=\frac{T}{2}-\frac{T}{2\pi }arccos\left({tan}^{2}s\right),$$

(49)

when s tends to $\pi /4$ from below.

Averaging ${\mathsf{\Omega }}_{\alpha }\left(t\right)$ over the period T, we arrive at:

$$\begin{array}{ccc}& & {\overline{\mathsf{\Omega }}}_{\mathrm{M}}\left(s\right)=\frac{3}{2}{cos}^{4}s-2{cos}^{2}s+1,\hfill \\ & & {\overline{\mathsf{\Omega }}}_{\mathrm{D}}\left(s\right)=\frac{1}{2}{cos}^{2}s,{\overline{\mathsf{\Omega }}}_{\mathrm{T}}\left(s\right)=\frac{3}{2}{cos}^{2}s{sin}^{2}s.\hfill \end{array}$$

(50)

Comparison of these fractions for different s is made in Figure 5. We see a slight dominance of ${\overline{\mathsf{\Omega }}}_{\mathrm{T}}$ in the range $s_c<s<\pi /4$, while ${\overline{\mathsf{\Omega }}}_{\mathrm{D}}$ prevails for $s<s_c$, where $s_c=\arcsin \left(\sqrt{3}/3\right)\simeq 3.13\pi /16$, and ${\overline{\mathsf{\Omega }}}_{\mathrm{M}}\left(s_c\right)={\overline{\mathsf{\Omega }}}_{\mathrm{D}}\left(s_c\right)={\overline{\mathsf{\Omega }}}_{\mathrm{T}}\left(s_c\right)=1/3$.

Additionally, let us indicate the changes in the energy spectrum when ${\widehat{H}}_0\ne 0$. Given ${\widehat{H}}_0=\mathrm{diag}(\mu ,\mu ,\mu )$, we get a shift in (45): ${\lambda }_i\to {\lambda }_i+\mu$. Including ${\widehat{H}}_0=\mathrm{diag}(0,-2Jsinh\left(a\right),0)$ leads to rescaling: ${\lambda }_1=J{\mathrm{e}}^{-a}$, ${\lambda }_2=0$, ${\lambda }_3=-J{\mathrm{e}}^a$ (see (46)). Dynamics in these cases will be somewhat more complicated, and we omit it.

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.

4. 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.

An 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 constant. This allows us to explore the properties of quasi-homogeneous subsystems. However, by 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:

$$n_0\simeq 5.61·{10}^{15}{\mathrm{m}}^{-3}{\left[\frac{m{c}^2}{1\mathrm{eV}}\right]}^{-1}\left[\frac{{\rho }_0}{{10}^{-20}\mathrm{kg}{\mathrm{m}}^{-3}}\right],$$

(51)

which is determined by the individual bosonic monomer mass $m\lesssim 1$ eV and the typical mass density ${\rho }_0$.

Let us consider possible reactions occurring simultaneously in each subsystem:

$$2\mathrm{M}\underset{k_{-1}}{\overset{k_1}{\rightleftharpoons }}\mathrm{D},\mathrm{M}+\mathrm{D}\underset{k_{-2}}{\overset{k_2}{\rightleftharpoons }}\mathrm{T},3\mathrm{M}\underset{k_{-3}}{\overset{k_3}{\rightleftharpoons }}\mathrm{T},$$

(52)

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_{\alpha }$ (the numbers of particles of sort $\alpha =\mathrm{M},\mathrm{D},\mathrm{T}$ in a fixed volume) can be described according to the rules of chemical kinetics [41]:

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}{C}_{\mathrm{M}}}{\mathrm{d}t}& =& -2k_1{V}_{2\mathrm{M}}{C}_{\mathrm{M}}^2+2k_{-1}{C}_{\mathrm{D}}-k_2{V}_{\mathrm{MD}}{C}_{\mathrm{M}}{C}_{\mathrm{D}}\hfill \\ & & +k_{-2}{C}_{\mathrm{T}}-3k_3{V}_{3\mathrm{M}}^2{C}_{\mathrm{M}}^3+3k_{-3}{C}_{\mathrm{T}},\hfill \\ \hfill \frac{\mathrm{d}{C}_{\mathrm{D}}}{\mathrm{d}t}& =& k_1{V}_{2\mathrm{M}}{C}_{\mathrm{M}}^2-k_{-1}{C}_{\mathrm{D}}-k_2{V}_{\mathrm{MD}}{C}_{\mathrm{M}}{C}_{\mathrm{D}}+k_{-2}{C}_{\mathrm{T}},\hfill \\ \hfill \frac{\mathrm{d}{C}_{\mathrm{T}}}{\mathrm{d}t}& =& k_2{V}_{\mathrm{MD}}{C}_{\mathrm{M}}{C}_{\mathrm{D}}-k_{-2}{C}_{\mathrm{T}}+k_3{V}_{3\mathrm{M}}^2{C}_{\mathrm{M}}^3-k_{-3}{C}_{\mathrm{T}},\hfill \end{array}$$

(53)

with the initial conditions at $t=0$: $C_{\mathrm{M}}^{\left(0\right)}>0$, $C_{\mathrm{D},\mathrm{T}}^{\left(0\right)}=0$; for $V_{2\mathrm{M}}$, $V_{3\mathrm{M}}$, and $V_{\mathrm{MD}}$, see next proposition.

The above set of equations for $C_{\alpha }$ conserves in time the total number of monomers in the subsystem of a fixed volume:

$$C_{\mathrm{M}}\left(t\right)+2C_{\mathrm{D}}\left(t\right)+3C_{\mathrm{T}}\left(t\right)=\mathrm{const},$$

(54)

and implies that the elementary act depends on constants $k_{\pm i}$ and interaction volumes $V_s\simeq a_s^3$, where $a_{2\mathrm{M}}$, $a_{3\mathrm{M}}$ and $a_{\mathrm{MD}}$ are the scattering lengths (for s-wave channel, for instance). It is natural to require conditions for successfully performing three forward reactions:

$$\begin{array}{ccc}& & V_{2\mathrm{M}}{C}_{\mathrm{M}}\ge 1,V_{3\mathrm{M}}{C}_{\mathrm{M}}\ge 1,\hfill \\ & & V_{\mathrm{MD}}{C}_{\mathrm{M}}\ge 1\mathrm{and}{V}_{\mathrm{MD}}{C}_{\mathrm{D}}\ge 1.\hfill \end{array}$$

(55)

These conditions 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_{\alpha }$ with the local particle density far from the core (of galactic halo) center to satisfy the imposed conditions (55), we require that $n_0{a}_s^3\gg 1$ for the scattering lengths $a_s$ in the model.

Strictly speaking, the validity of condition $n{a}_s^3>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:

$$K_1=k_1{V}_{2\mathrm{M}},K_2=k_2{V}_{\mathrm{MD}},K_3=k_3{V}_{3\mathrm{M}}^2,K_{-i}=k_{-i}.$$

It is convenient to describe the subsystem evolution in terms of the fractions:

$${\mathsf{\Omega }}_{\alpha }=\frac{c_{\alpha }{C}_{\alpha }}{C_{\mathrm{M}}+c_{\mathrm{D}}{C}_{\mathrm{D}}+c_{\mathrm{T}}{C}_{\mathrm{T}}},{\mathsf{\Omega }}_{\mathrm{M}}+{\mathsf{\Omega }}_{\mathrm{D}}+{\mathsf{\Omega }}_{\mathrm{T}}=1,$$

(56)

where $c_{\mathrm{M}}=1$, $c_{\mathrm{D}}=2$, $c_{\mathrm{T}}=3$.

Let us analyze the equilibrium fractions ${\mathsf{\Omega }}_{\alpha }^{\mathrm{eq}}$ when $\mathrm{d}{C}_{\alpha }/\mathrm{d}t=0$. For equilibrium concentrations (dropping the label ‘eq’), we find $C_{\mathrm{D}}$ and $C_{\mathrm{T}}$ as functions of $C_{\mathrm{M}}$:

$$\begin{array}{ccc}& & C_{\mathrm{D}}=C_{\mathrm{M}}^2{f}_{\mathrm{D}}\left(C_{\mathrm{M}}\right),C_{\mathrm{T}}=C_{\mathrm{M}}^3{f}_{\mathrm{T}}\left(C_{\mathrm{M}}\right),\hfill \end{array}$$

(57)

$$\begin{array}{ccc}& & f_{\alpha }\left(x\right)=\frac{a_{\alpha }^{\left(11\right)}x+a_{\alpha }^{\left(12\right)}}{a_{\alpha }^{\left(21\right)}x+a_{\alpha }^{\left(22\right)}},A_{\alpha }=\left(\begin{array}{cc}{a}_{\alpha }^{\left(11\right)}& a_{\alpha }^{\left(12\right)}\\ a_{\alpha }^{\left(21\right)}& a_{\alpha }^{\left(22\right)}\end{array}\right),\hfill \end{array}$$

(58)

where the elements of matrix $A_{\alpha }$ are formed from $K_{\pm i}$:

$$\begin{array}{ccc}& & A_{\mathrm{D}}=\left(\begin{array}{cc}{K}_3{K}_{-2}& K_1(K_{-2}+K_{-3})\\ K_2{K}_{-3}& K_{-1}(K_{-2}+K_{-3})\end{array}\right),\hfill \end{array}$$

(59)

$$\begin{array}{ccc}& & A_{\mathrm{T}}=\left(\begin{array}{cc}{K}_2{K}_3& K_1{K}_2+K_3{K}_{-1}\\ K_2{K}_{-3}& K_{-1}(K_{-2}+K_{-3})\end{array}\right).\hfill \end{array}$$

(60)

For a subsystem in a fixed (unit) volume, the initial concentration $C_{\mathrm{M}}^{\left(0\right)}$ is related to the equilibrium one $C_{\mathrm{M}}$:

$$C_{\mathrm{M}}+2C_{\mathrm{M}}^2{f}_{\mathrm{D}}\left(C_{\mathrm{M}}\right)+3C_{\mathrm{M}}^3{f}_{\mathrm{T}}\left(C_{\mathrm{M}}\right)=C_{\mathrm{M}}^{\left(0\right)}.$$

(61)

It means that the initial state $C_{\mathrm{M}}^{\left(0\right)}$ can be determined by using the $C_{\mathrm{M}}$.

Note the possibility of generalizing our considerations by specifying the concentration $C_{\alpha }^{\left(q\right)}$ of composites of sort $\alpha$ in certain quantum state q at a given temperature T, that is, the occupation of energy level $E_{\alpha }^{\left(q\right)}$ with the degeneracy $g_q$. For the processes that do not violate the equilibrium distribution of particles over states, one writes [41]:

$$C_{\alpha }^{\left(q\right)}=C_{\alpha }{p}_{\alpha }\left(q\right),p_{\alpha }\left(q\right)=\frac{g_q}{Z_{\alpha }}exp\left(-\frac{E_{\alpha }^{\left(q\right)}}{T}\right),\sum _q{p}_{\alpha }\left(q\right)=1,$$

(62)

where $Z_{\alpha }$ is the partition function over states, and $C_{\alpha }$ is the total concentration. From Equation (53), with account of Boltzmann factors $p_{\alpha }\left(q\right)$, one can deduce the dependence of the reaction rates $k_{\pm i}$ on nonzero temperature.

At finite $K_{\pm i}$, there are two limiting regimes: (1) $C_{\mathrm{M}}\to 0$ (low concentration of monomers) admits approximation: $C_{\mathrm{M}}^{\left(0\right)}=C_{\mathrm{M}}+O\left(C_{\mathrm{M}}^2\right)$, and

$$\begin{array}{ccc}& & {\mathsf{\Omega }}_{\mathrm{M}}^{\mathrm{eq}}=1-2\frac{K_1}{K_{-1}}{C}_{\mathrm{M}}+O\left(C_{\mathrm{M}}^2\right),\hfill \\ & & {\mathsf{\Omega }}_{\mathrm{D}}^{\mathrm{eq}}=2\frac{K_1}{K_{-1}}{C}_{\mathrm{M}}+O\left(C_{\mathrm{M}}^2\right),\hfill \\ & & {\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}=3\frac{K_1{K}_2+K_3{K}_{-1}}{K_{-1}(K_{-2}+K_{-3})}{C}_{\mathrm{M}}^2+O\left(C_{\mathrm{M}}^3\right);\hfill \end{array}$$

(63)

(2) asymptotic expressions at large $C_{\mathrm{M}}$ yield:

$$\begin{array}{ccc}& & C_{\mathrm{M}}^{\left(0\right)}=\frac{3K_3}{K_{-3}}{C}_{\mathrm{M}}^3+o\left(C_{\mathrm{M}}^2\right),\hfill \\ & & {\mathsf{\Omega }}_{\mathrm{M}}^{\mathrm{eq}}=\frac{K_{-3}}{3K_3}{C}_{\mathrm{M}}^{-2}+O\left(C_{\mathrm{M}}^{-3}\right),{\mathsf{\Omega }}_{\mathrm{D}}^{\mathrm{eq}}=\frac{2K_{-2}}{3K_2}{C}_{\mathrm{M}}^{-1}+O\left(C_{\mathrm{M}}^{-2}\right),\hfill \\ & & {\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}=1-\frac{2K_{-2}}{3K_2}{C}_{\mathrm{M}}^{-1}+O\left(C_{\mathrm{M}}^{-2}\right).\hfill \end{array}$$

(64)

The formulas (63) and (64) are obtained by expanding up to the first nontrivial terms. Thus, the normalization ${\mathsf{\Omega }}_{\mathrm{M}}^{\mathrm{eq}}+{\mathsf{\Omega }}_{\mathrm{D}}^{\mathrm{eq}}+{\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}=1$ is valid for ${\mathsf{\Omega }}_{\alpha }^{\mathrm{eq}}$ in the same approximation.

As seen from (63) and (64), the case (1) indicates a rather low probability of obtaining composites D and T, while the mode (2) shows the dominant role of trimers (due to reaction $\mathrm{M}+\mathrm{D}\to \mathrm{T}$) and the suppression of the monomer number because of converting $3\mathrm{M}\to \mathrm{T}$. In this case we also see that $C_{\mathrm{M}}\sim {\left(C_{\mathrm{M}}^{\left(0\right)}\right)}^{1/3}$.

Since the admissible value of ${\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}$ varies in the interval $(0;1)$, we can obtain ${\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}$ higher than the expected threshold ${\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{thr}}\in (0;1)$ by requiring $C_{\mathrm{M}}>C_{\mathrm{M}}^{\mathrm{thr}}$. For a given ${\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{thr}}$, the value of $C_{\mathrm{M}}^{\mathrm{thr}}$ is determined by a real positive root z of the quartic equation:

$$3z^3{f}_{\mathrm{T}}\left(z\right)\left(1-{\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{thr}}\right)-{\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{thr}}\left(z+2z^2{f}_{\mathrm{D}}\left(z\right)\right)=0.$$

(65)

A similar problem arises when expressing $C_{\mathrm{M}}$ through $C_{\mathrm{M}}^{\left(0\right)}$. In addition, one can also find the threshold value of $C_{\mathrm{M}}^{\left(0\right)}$ when ${\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}={\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{thr}}$. To demonstrate some analytical solutions, we restrict ourselves to simpler special cases.

Let $K_{-1}=K_{-2}=K_3=0$. It means that the allowed processes are:

$$3\mathrm{M}\stackrel{k_1}{\to }\mathrm{M}+\mathrm{D}\stackrel{k_2}{\to }\mathrm{T}\stackrel{k_{-3}}{\to }3\mathrm{M}.$$

(66)

The composites fractions in such a model are given by:

$$\begin{array}{ccc}& & {\mathsf{\Omega }}_{\mathrm{M}}^{\mathrm{eq}}=\frac{K_2{K}_{-3}}{K_2{K}_{-3}+2K_1{K}_{-3}+3K_1{K}_2{C}_{\mathrm{M}}},\hfill \\ & & {\mathsf{\Omega }}_{\mathrm{D}}^{\mathrm{eq}}=\frac{2K_1{K}_{-3}}{K_2{K}_{-3}+2K_1{K}_{-3}+3K_1{K}_2{C}_{\mathrm{M}}},\hfill \\ & & {\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}=\frac{3C_{\mathrm{M}}{K}_1{K}_2}{K_2{K}_{-3}+2K_1{K}_{-3}+3K_1{K}_2{C}_{\mathrm{M}}},C_{\mathrm{M}}={\mathcal{C}}_{+}\left(C_{\mathrm{M}}^{\left(0\right)}\right),\hfill \end{array}$$

(67)

where ${\mathcal{C}}_{+}$ is a positive solution among the two solutions ${\mathcal{C}}_{\pm }$ to the quadratic equation obtained from (61) under imposed restrictions:

$${\mathcal{C}}_{\pm }=\sqrt{\frac{K_{-3}}{3K_1}}\left[\pm \sqrt{C_{\mathrm{M}}^{\left(0\right)}+\mathcal{C}}-\sqrt{\mathcal{C}}\right],\mathcal{C}=\frac{K_{-3}}{12K_1}{\left(1+2\frac{K_1}{K_2}\right)}^2.$$

(68)

Thus, one has that $C_{\mathrm{M}}\sim \sqrt{C_{\mathrm{M}}^{\left(0\right)}}$ within this model.

To guarantee ${\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}\ge {\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{thr}}$ for ${\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{thr}}=0.5$, we arrive at the condition:

$$C_{\mathrm{M}}\ge \frac{K_{-3}}{3K_1}+\frac{2K_{-3}}{3K_2}\mathrm{or}{C}_{\mathrm{M}}^{\left(0\right)}=4\mathcal{C}\gamma (\gamma +1),\gamma \ge 1,$$

(69)

with definite $\gamma$.

Rewriting this as:

$$\frac{1}{3}\frac{k_{-3}}{k_1}\frac{1}{C_{\mathrm{M}}{a}_{2\mathrm{M}}^3}+\frac{2}{3}\frac{k_{-3}}{k_2}\frac{1}{C_{\mathrm{M}}{a}_{\mathrm{MD}}^3}\le 1,$$

(70)

one has two auxiliary inequalities: $C_{\mathrm{M}}{a}_{2\mathrm{M}}^3\ge k_{-3}/k_1$ and $C_{\mathrm{M}}{a}_{\mathrm{MD}}^3\ge k_{-3}/k_2$, which are obviously satisfied at $k_{-3}/k_1\le 1$, $k_{-3}/k_2\le 1$, and $C_{\mathrm{M}}{a}_s^3\ge 1$ from (55).

Now let $K_3=K_{-3}=0$ in (52). Denoting ${\varkappa }_i=K_{-i}/K_i$, we obtain the fractions:

$$\begin{array}{ccc}& & {\mathsf{\Omega }}_{\mathrm{M}}^{\mathrm{eq}}=\frac{{\varkappa }_1{\varkappa }_2}{{\varkappa }_1{\varkappa }_2+2{\varkappa }_2{C}_{\mathrm{M}}+3C_{\mathrm{M}}^2},{\mathsf{\Omega }}_{\mathrm{D}}^{\mathrm{eq}}=\frac{2{\varkappa }_2{C}_{\mathrm{M}}}{{\varkappa }_1{\varkappa }_2+2{\varkappa }_2{C}_{\mathrm{M}}+3C_{\mathrm{M}}^2},\hfill \\ & & {\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}=\frac{3C_{\mathrm{M}}^2}{{\varkappa }_1{\varkappa }_2+2{\varkappa }_2{C}_{\mathrm{M}}+3C_{\mathrm{M}}^2},C_{\mathrm{M}}={\mathcal{C}}_0\left(C_{\mathrm{M}}^{\left(0\right)}\right),\hfill \end{array}$$

(71)

where ${\mathcal{C}}_0$ denotes a real solution among three roots ${\mathcal{C}}_m$ ($m=0,\pm 1$) of a cubic equation:

$$\begin{array}{ccc}& & {\mathcal{C}}_m=\mathcal{A}{\mathrm{e}}^{\mathrm{i}\frac{2\pi }{3}m}-\frac{{\mathcal{A}}^{\left(2\right)}}{\mathcal{A}}{\mathrm{e}}^{-\mathrm{i}\frac{2\pi }{3}m}-\frac{2}{9}{\varkappa }_2,\mathcal{A}={\left[{\mathcal{A}}^{\left(3\right)}+\sqrt{{\left({\mathcal{A}}^{\left(3\right)}\right)}^2+{\left({\mathcal{A}}^{\left(2\right)}\right)}^3}\right]}^{1/3},\hfill \end{array}$$

(72)

$$\begin{array}{ccc}& & {\mathcal{A}}^{\left(2\right)}=\frac{{\varkappa }_2}{9}\left({\varkappa }_1-\frac{4}{9}{\varkappa }_2\right),{\mathcal{A}}^{\left(3\right)}=\frac{{\varkappa }_1{\varkappa }_2}{6}\left[C_{\mathrm{M}}^{\left(0\right)}+\frac{2}{9}{\varkappa }_2\left(1-\frac{8}{27}\frac{{\varkappa }_2}{{\varkappa }_1}\right)\right].\hfill \end{array}$$

(73)

To ensure that ${\mathsf{\Omega }}_{\mathrm{T}}^{\mathrm{eq}}\ge 0.5$, we demand $3C_{\mathrm{M}}^2\ge {\varkappa }_1{\varkappa }_2+2{\varkappa }_2{C}_{\mathrm{M}}$ (see (65)), or, equivalently,

$$C_{\mathrm{M}}\ge \frac{{\varkappa }_2}{3}\left(1+\sqrt{1+3\frac{{\varkappa }_1}{{\varkappa }_2}}\right).$$

(74)

It follows that $C_{\mathrm{M}}{a}_{\mathrm{MD}}^3>k_{-2}/k_2$, where $k_{-2}/k_2\le 1$ to permit the synthesis of T. Thus, condition (74) is not stronger than (55).

To summarize, the applied analysis enables the leading role of composites under the derived conditions (55) and the definite values of monomer density $C_{\mathrm{M}}^{\left(0\right)}$ and $C_{\mathrm{M}}$ (see (69) and (74)). However, the mechanism of the formation of composites needs detailed study.

5. 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, which is based on the Gross–Pitaevskii equation extended by an additional three-particle interaction, we find that few-particle complexes can be formed at relatively low density and pressure when quantum fluctuations and the transfer momentum is small enough and cannot destroy the condensate and the desired structures. Such a conclusion is deduced on the basis of calculating the averages of the powers of the local density $\eta$ and thermodynamic functions at zero temperature. However, increasing the pressure (and density), the multiparticle averages are reduced to the product of single-particle averages that determine the dependence on thermodynamic parameters. This allows us 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 [25]. In addition, as found in ref. [25], 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 three-particle interaction arises due to quantum corrections in a system with pairwise interaction and the use of renormalization group calculations [34,42,43,44]:

$$\begin{array}{ccc}& & \mathcal{E}=\frac{\mathcal{Q}\left(n\right)}{2}{n}^2+\frac{\mathcal{B}\left(n\right)}{3}{n}^3,\hfill \end{array}$$

(75)

$$\begin{array}{ccc}& & \mathcal{Q}\left(n\right)=\frac{4\pi {\hslash }^{2}a}{m}\left(1+\frac{128}{15\sqrt{\pi }}\sqrt{n{a}^3}\right),\hfill \end{array}$$

(76)

$$\begin{array}{ccc}& & \mathcal{B}\left(n\right)=\frac{16\pi {\hslash }^2{a}^4}{m}\left(4\pi -3\sqrt{3}\right)\left[ln\left(n{a}^3\right)+2d\right].\hfill \end{array}$$

(77)

a Here n is the particle density of homogeneous system of bosons, a is their pairwise s-scattering length, and d is a numeric parameter.

Note that the internal energy $\mathcal{E}$ is proportional to the inverse of the mass of (dark matter) 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 [45,46]. This is motivated, in particular, by the study of Bose condensates of alkali atoms, the strong interaction between which can cause molecule 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 an appropriate mechanism for forming the trimers from dark bosons. The existence of such an effect is successfully confirmed by experiments with cold-atom systems [47] and is investigated [36]. The 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.