Phases of the Bose-Einstein condensate dark matter model with both two- and three-particle interactions

In this paper we further elaborate on the Bose-Einstein condensate (BEC) dark matter model extended in our preceding work [Phys. Rev. D 102, 083510 (2020)] by the inclusion of 6th order (or three-particle) repulsive self-interaction term. Herein, our goal is to complete the picture through adding to the model the 4th order repulsive self-interaction. The results of our analysis confirm the following: while in the preceding work the two-phase structure and the possibility of first-order phase transition was established, here we demonstrate that with the two self-interactions involved, the nontrivial phase structure of the enriched model remains intact. For this to hold, we study the conditions which the parameters of the model, including the interaction parameters, should satisfy. As a by-product and in order to provide some illustration, we obtain the rotation curves and the (bipartite) entanglement entropy for the case of particular dwarf galaxy.


Introduction
Although the concept of dark matter (DM) is a widely accepted one, its precise nature is still escaping. There exist vast multitude of different approaches and models, among which the modeling of DM as Bose-Einstein condensate (BEC), see e.g. [1,2,3,4], the overviews [5,6] and many others, finds each time more and more support. Positions of BEC model of DM were especially enforced after the works [7,8,9] which demonstrated ability of BEC DM to avoid the core-cusp and the gravitational collapse [10] problems. Also, the model gives quite successful description [7,11,12,13] of the rotation curves of a number of galaxies, at least dwarf and the low surface brightness ones.
Nevertheless, even within this well-elaborated model there also exist some tensions and issues which can be improved. To this end, there is in particular a possibility to apply appropriate tools from the powerful and efficient theory of deformations. Namely, the µ-deformed analog of Bose-gas model developed in [14], with the so-called µ-calculus as a base, has demonstrated clearly the following preferable features: (i) the evaluated mass of DM halo appears more realistic; (ii) the obtained critical temperature of condensation of µ-Bose gas T (µ) C depends on the deformation parameter µ, µ > 0, and is higher [15] than the usual T C ; (iii) µ-deformation based description of the rotation curves [16] fits better than the curves inferred within the ordinary BEC model.
It is also worth to mention the recent work [17] which uses the concept of deformed spatial commutation relations for scalar field in order to develop a class of generalization of the Bosecondensate DM model. Such an extension has good potential to achieve improvements.
Another line of extension of the BEC model of DM, developed recently in [18], involves sixth order (or 3-particle) self-interaction term ψ 6 . Due to presence of the latter, the modified model manifests nontrivial phase structure: there exist two distinct phases, certain region of instability, and the possibility of first-order phase transition.
In order to make the extended model even more complete, it is natural to include, besides the ψ 6 , also the two-particle self-interaction encoded in the term ψ 4 . Analysis of such "doubly-nonlinear" extension of [18] and BEC model of DM is the goal of the present paper.
The structure of the paper is the following. Necessary details of the model are given in Section 2, and main part that involves obtaining thermodynamic functions and their key properties is presented in Section 3. In Sections 4 and 5 we consider briefly, again for the situation of presence of both two-and three-particle interactions, the rotation curves of selected galaxy and the respective bipartite entanglement entropy of two centrally symmetric regions of the halo of this same galaxy. In our concluding section, we present a discussion of the results.

The Model
We describe the Bose-Einstein condensate by real function ψ(r) of radial variable r = |r|, using a constant chemical potentialμ. Our study is based on the energy functional Γ in a ball B = {r ∈ R 3 | |r| ≤ R} and the Poisson equation: where ∆ r is the radial part of Laplace operator that acts as R being the radius of the ball where the matter is located. Focusing here on the effects of interparticle interactions, we leave aside the slow rotation of the condensate [13], which can be taken into account through the chemical potential [17].
For convenience, let us introduce dimensionless variables: Here χ(ξ) is a real dimensionless scalar field; ̺ 0 and r 0 characterize typical measures of the central particle density and the system size, respectively. Thus, we arrive at where ξ B = R/r 0 , whereas ∆ ξ and ∆ −1 ξ are given by (3), (4) in terms of ξ replacing r. To estimate the range of parameter values, we turn to astrophysical situations. Since we suggest to take into account the three-particle interaction in relatively dense DM of light bosons with masses of the order of 10 −22 eV c −2 , ranges of the parameters can be found by considering the DM of galactic cores with a central mass density ρ 0 = m̺ 0 of the order of 10 −20 kg m −3 , in the region of radius r 0 smaller than 1 kpc. Then, extracting r 0 from the definition (5) of the measure of gravitational interaction A, we can adopt that A ∼ 10 [18]. It is clear that the gravity results from integral effect of a whole system. Unlike, (thermo)dynamics of internal processes is determined by repulsive interactions among bosons, represented by the parameters Q and B under the condition B > A. The role of pairwise interaction, controlled by Q, is assumed to be comparable with the effect of gravity.
In order to find a decreasing solution χ(ξ) for admissible ξ with a finite initial value χ 0 = χ(0) < ∞, we first impose χ ′ (0) = 0 and then formulate the conditions which allow to fix χ 0 , through expanding χ(ξ) = χ 0 + C 2 ξ 2 + . . . at ξ → 0. On substituting that in (12), (13), the following algebraic equations result: Combining the two, we find that the value χ 0 should satisfy the equation S(A * , B * , Q * , ν, χ 0 ) = 0, where plus the condition 2C 2 = χ ′′ (0) ≤ 0. Taken altogether, these constrain χ 0 as z 1 < χ 0 < z 2 , where Technically, the search for the initial value χ 0 of the model which takes into account the pair interaction is similar to the problem with three-particle interaction only [18]. Likewise, we notice three regimes (for given A * , Q * , B * and ν): 1) no solution for χ 0 that, in Eq. (12), leads to χ(ξ) = 0; 2) single solution χ 0 that corresponds to a minimal admissible value ν min from which the system starts to evolve; 3) pair of (positive) solutions for χ 0 , when we should choose a minimal one, because the other leads to divergent χ(ξ). Usually, for fixed (A * , Q * , B * ), but increasing ν, the indicated sequence of all three options occurs. It is useful to analyze the system from the quantum-mechanical point of view. Equation (12) for ξ ≤ ξ B can be conveniently rewritten in the Schrödinger form to describe the scattering of a particle whose wave function 2 is taken as f (ξ) = cχ(ξ) (c is a normalization depended on a total number of particles) in the potential V eff : where potentials V 2 and V 3 come from two and three-particle interactions, respectively. Substituting the found solution χ(ξ), we can see that the form of potential V eff depends also on a chemical potential u (or parameter ν).
In Fig. 1 we have used the values A = 10, B = 20 and Q = 1.36. The latter one plays the role of "critical" value (its sense is seen in Fig. 2 below, with explanations at the end of the next Section). We relate the particular forms of V eff given in Fig. 1a with the physical situations depicted in Fig. 2 below. Namely, green curve in Fig. 1a is chosen for liquid-like (dense) state in Fig. 2, when the three-particle interaction contributes to a hard-core part of potential at small ξ. Red curve is constructed in the vicinity of the critical point of the first-order phase transition, when the 2 See Fig. 1b for the behavior of χ(ξ). potential V eff is similar to the harmonic trap. Black curve corresponds to gaseous (dilute) state in the effective (truncated) gravitational potential, when the kinetic energy term dominates (u > 0).
In the range ξ ≥ ξ B we come to the problem of a particle in the gravitational field (see the dashed curves in Fig. 1a) created by the system of N particles: where is the total number of particles within the ball ξ ≤ ξ B , which determines the total mass. At this stage a wave number k is ambiguous. Oscillating and decaying solution to this equation (for any real k and pure imaginary κ and s) is given as: Here M µ,ν (z) and W µ,ν (z) are the Whittaker functions. This solution might be interpreted as the gravitational capture of a dark matter particle (a radial wave with energy k 2 /2) outside the dark matter ball, when ξ → ∞. In the case of initially resting particle with k = 0, the solution is described in terms of Bessel functions J 1 (z) and Y 1 (z).
Thus, the account of interaction confirms the possibility of bound states, which can manifest themselves in the form of different thermodynamic phases. Further, all global quantities of the model, computed at fixed A, Q and B, are supposed to be functions of free parameter ν. Therefore, dependence, say, of a on b should be treated in parametric form: a(b) = {(b(ν), a(ν))|ν ≥ ν min }.

Thermodynamic Quantities and Two Phases
To study the macroscopic properties, let us define an effective chemical potential µ(ξ) [19], which includes the gravitational potential jointly with the term of quantum fluctuations, and replaces further the constant chemical potential u, that is Accordingly to the equation of motion (9), µ determines χ as For finding macroscopic characteristics, we appeal to the thermodynamic relations at T = 0, using a local particle density η(ξ) = χ 2 (ξ): where functions p(ξ) and ε(ξ) determine the (dimensionless) mean pressure P and the internal energy E: Hereafter, ξ 3 B /3 represents the volume of the system. Therefore, we need to integrate the Gibbs-Duhem relation (29) and then to substitute p(ξ) into the Euler relation (30) in order to find ε(ξ). This way, the explicit expressions are obtained: which give us the equation of state by inserting the solution η(ξ) of (12)- (14). The (local) internal pressure p(ξ) behaves spatially with radius accoding to We do not solve it because all necessary fields are found explicitly from Eqs. (12)- (14). Mean particle density is Dependence of mean internal pressure P , see (31), on this density σ is presented in Fig. 2 for various parameters of interaction. The curves obtained numerically reveal the presence of two stable phases of dark matter with ∂P/∂σ > 0: the dilute one (∂ 2 P/∂σ 2 < 0) and the denser liquidlike one (∂ 2 P/∂σ 2 > 0). The characteristic points of the phase diagram are determined from the conditions: ∂P/∂σ = 0 and ∂ 2 P/∂σ 2 = 0. The simultaneous fulfillment of these conditions at a single point determines the critical point, which belongs to the black curve in Fig. 2, and inferring of which is one of main tasks here. The existence of a critical point of a first-order phase transition is essentially related to the competition between gravitational and pair interactions, controlled by the parameter Q.
While the presence of metastable states, linked with the two-extrema behavior at Q > 1.36 and shown by the orange curve, clearly indicates mixing of the gaseous and liquid phases, a detailed description of the transition between the two phases at Q < 1.36 requires the use of additional thermodynamic characteristics, as is already done in [18] by means of the perturbation pressure Π ν for Q = 0.
It is important to emphasize that herein we reveal a discontinuous behavior (jump between two phases) of the density σ with the change in the internal pressure P . The existence of such a regime, as we show, is allowed at Q > 1.36 for A = 10 and B = 20. Although the existence of a denser, liquidlike phase of DM (in a relatively small region of halo) is possible at Q = 0, the condition Q > 0 is required in describing the galactic DM halos [7], even without taking into account the three-particle interaction.
Note that a continuous change in the parameter Q can lead to the intersection of different curves P (σ), that indicates the possibility of realizing one thermodynamic state by fixation of different sets of parameters. Such ambiguity in the set of parameters can complicate the interpretation and reproduction of observables.

The Rotation Curves
Important information about dark matter is extracted and verified from the rotation curves of galaxies. As our model is aimed to study the processes in dark matter in the galaxy cores, that is, in relatively small regions of space with a noticeable density, its direct application to the description of rotation curves should be limited to dwarf galaxies (or other dark matter dominated compact galaxies). One of the possibilities for describing larger (halo-type) objects is to extend the model, similarly to what was proposed in [18]. Besides, the lack of accounting for the rigid rotation of matter in the model does not allow us to describe the curves of rotating galaxies. Nevertheless, we consider it important and interesting to demonstrate the capabilities of our model, in which we accurately took into account quantum fluctuations in the condensate and (both two-and) threeparticle interactions. Although the extensions of the model are indicated in [18,17], the included effects make it possible to further validate the Bose-condensate approach, to compare with and supplement the results of [7].
Thus, the tangential velocity v of a test particle moving in the spherically symmetric DM halo can be represented as where ρ(r) = m|ψ(r)| 2 is a mass density such that ρ(R) = 0 at R = r 0 ξ B . Let us illustrate the model predictions for the rotation curves of the M81 galaxy, shown in Fig. 3. Note especially the free parameters of the model A, B, Q, ν, ρ 0 that we use for fitting, as well as the restricting characteristics: the total mass M of dark matter and the halo radius R.
The galaxy M81 with no rigid rotation gives us the most striking proof of the Bose-condensate approach, as noted in [7]. Let us use this example to emphasize main features of the description of rotation curves.
First of all, note it looks rather difficult to indicate unambiguously the parameters for the rotation curve: the same dependence v(r) can be realized for different sets of the parameters. This is a consequence of the symmetry of the model equations. For this reason, we give only graphs that are interesting for physics, and omit the mathematics of identifying the symmetries.
In our example, it is worth to note a possibility to construct the curve colored in cyan in Fig. 3 with Q = 0, when the pair interaction is absent. At the same time, the best fitted dependencies,  represented by other curves, always require Q > 0. This implies that the repulsive pair interaction should be taken into account in realistic models of dark matter in (dwarf) galaxies [11,20].

Entanglement Entropy
The possibility of estimating the entanglement entropy in a general Bose condensate was developed, in particular, in [21]. The idea of using this entropy as a criterion for differentiating the BEC dark matter from cold dark matter (CDM) was proposed in [22]. Although a more precise analysis requires studying the effects of interference of particles from interacting subsystems, here we estimate the entropy between two separated (central and surrounding) parts of radially inhomogeneous BEC dark matter of a galaxy. Omitting the normalization constants, we use the formula (in accordance with [21,22]): Here c(ξ) = n(ξ)/N is the fraction of particles contained in the central subsystem, n(ξ) being defined in (33). A typical dependence of the entanglement entropy on the specific radius is shown in Fig. 4. It is worth noting that, although the system is inhomogeneous, and c(ξ) is not a constant, the maximum influence of one subsystem of particles upon the other turns out to occur at x = 0.5. What concerns negative sign of S E (x): that is eliminated by explicit account of normalization constant, that is, by adding definite positive number (of the order 10 2 or higher, see e.g. [22]) which depends on the mass of dark matter particles and their number in the subregion.

Discussion
In the framework of the extended version of BEC dark matter model that involves two-and threeparticle interactions we explored the solutions of the system of equations (11)- (14) as well as the properties of thermodynamic functions. The analysis has led us to the conclusions that the interplay of the basic free parameters A, Q and B gives interesting consequences for the system under study including dark matter halo properties: (i) first, the effective interaction potential, see (18)- (21) and Fig. 1, and the density profiles; (ii) the equation of state which shows basic dependence on the parameter Q including existence of its "critical" value Q C that separates different regimes; (iii) persistence at Q > Q C of nontrivial phase structure (inherited from the restricted pure ψ 6 sub-model). The obtained information on the thermodynamic function such as the mass density profile, enabled us to infer in Sec. 3 the galactic rotation curves, and the bipartite entanglement entropy in Sec. 4. The former clearly demonstrated, with the example of dwarf galaxy M81, that the model can easily provide nice agreement with observational data. As follows from the treatment of rotational curves, the parameter Q responsible for the two-particle interaction plays essential role. Besides, as seen in Fig. 2, at Q > Q C the model reveals the most rich behavior: the curves possess two extrema that implies the metastability region. At Q ≥ Q C , only inflection points do survive. Without Q (without pair interaction) we could not have such features of unexpected behavior. Anyway, for both Q > 0 and Q = 0 (that brings us back to the situation explored in [18]) the two phases are present, though their identification involves differing thermodynamical functions. It is of interest to compare the model considered in this paper with some of the nonlinear (namely, non-polynomial) extensions of the BEC model of DM, e.g. those studied in [23,24,25,26] where the scalar potential V 0 [cosh (λκΦ) − 1] was employed, with κ = √ 8πG. The essential feature of these models is the presence of all-order nonlinearities, with the corresponding powers of the single free parameter λ (note that λ can be viewed as a deformation parameter since at λ → 0 the potential and thus the interaction are vanishing). On one hand, the model studied above is obviously simpler than the models involving cosh or cos: indeed, we encounter the first terms of series expansion if the parameters are restricted as A = √ B = λ. On the other hand, the ∼ ψ 4 plus ∼ ψ 6 model studied herein is richer in the sense that it operates with the free parameters A, B, Q, instead of single one. It is this property that allowed us to disclose (confirm) the existence of two phases and of the phase transition. Moreover, the influence of different values of these parameters was essential in our treatment of the rotation curves and the entanglement entropy, see Fig. 4 above.
Concerning the entanglement, an interesting question arises for the situation when dark matter particles are not elementary bosons, but composites built from two bosons or two fermions. In the both cases (i) the composites differ from pure bosons and are naturally realizable, as can be seen in [27], through deformed oscillators or deformed bosons (or quasi-bosons); (ii) bipartite internal entanglement entropy of quasibosons obtained in [28,29] turned out to depend on the deformation parameter. The question now is as follows: to which extent the microscopic intra-quasibosonic entanglement and its entropy do affect (superimpose with) the macroscopic entanglement and the entanglement entropy that was studied in Sec. 5 above. Also it is no doubt important to explore in detail the connection [21,30,31,32] between peculiarities of the behavior of entanglement and the phase transition (of the first order in our case), especially in the context of the properties of dark matter. We hope to explore these questions in one of our future works.