© ESO, 2012

1. Introduction

Several astrophysical observations of distant type Ia supernovae have revealed that the content of the universe is made of about 70% dark energy, 25% dark matter, and 5% baryonic (visible) matter (Riess et al. 1998; Perlmutter et al. 1999; de Bernardis et al. 2000; Hanany et al. 2000). Thus, the overwhelming preponderance of matter and energy in the universe is believed to be dark, i.e. unobservable by telescopes. The dark energy is responsible for the accelerated expansion of the universe. Its origin is mysterious and presumably related to the cosmological constant. Dark energy is usually interpreted as a vacuum energy, and it behaves like a fluid with negative pressure. Dark matter is also mysterious. The suggestion that dark matter may constitute a large part of the universe was raised by Zwicky in 1933. Using the virial theorem to infer the average mass of galaxies within the Coma cluster, he obtained a much higher value than the mass of luminous material. He realized therefore that some mass was “missing” to account for observations. This missing mass problem was confirmed later by accurate measurements of rotation curves of disk galaxies (Rubin et al. 1980; Persic et al. 1996). The rotation curves of neutral hydrogen clouds in spiral galaxies measured from the Doppler effect are found to be roughly flat (instead of Keplerian) with a typical rotational velocity v_∞ ~ 200   km   s^-1 up to the maximum observed radius of about 50 kpc. The corresponding mass profile is much more extended than the distribution of starlight, which typically converges within  ~10 kpc. This implies that galaxies are surrounded by an extended halo of dark matter, whose mass $\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{~}\mathit{r}{\mathit{v}}_{\mathrm{\infty }}^{\mathrm{2}}\mathit{/}\mathit{G}$ increases linearly with radius. This can be conveniently modeled by a classical isothermal gas whose density decreases as r^-2 (Chandrasekhar 1939). Although some authors like Milgrom (1983) propose a modification of Newton’s law (MOND theory) to explain the rotation curves of spiral galaxies without invoking dark matter, the dark matter hypothesis is favored by most astrophysicists.

The nature of dark matter (DM) is one of the most important puzzles in modern physics and cosmology. A wide “zoology” of exotic particles that could form dark matter has been proposed. In particular, many grand unified theories in particle physics predict the existence of various exotic bosons (e.g. axions, scalar neutrinos, neutralinos) that should be present in considerable abundance in the universe and comprise (part of) the cosmological missing mass (Primack et al. 1998; Overduin & Wesson 2004). Even if the bosonic particles have never been detected in accelerator experiments, they are considered as leading candidates of dark matter that might play a significant role in the evolution and structure of the universe.

If dark matter is made of bosons, it should have formed compact gravitating Bose-Einstein condensates (BEC) such as boson stars. Boson stars were introduced by Kaup (1968) and Ruffini & Bonazzola (1969) in the 1960s. Early works on boson stars (Thirring 1983; Breit et al. 1984; Takasugi & Yoshimura 1984; van der Bij & Gleiser 1987) were motivated by the axion field that was proposed as a possible solution to the strong CP problem in QCD. For particles with mass m ~ 1   GeV/c², the maximum mass of a boson star, called the Kaup mass, is much lower than the solar mass (M_Kaup ~ 10^-19   M_⊙!), so that these mini boson stars, such as axion black holes, are not very astrophysically relevant. They could play a role, however, if they exist in the universe in abundance or if the axion mass is extraordinary low (Schunck & Mielke 2003), leading to macroscopic objects with a mass M_Kaup comparable to the mass of the sun (or even higher). For example, axionic boson stars could account for the mass of MACHOs (between 0.3 and 0.8   M_⊙) if the axions have a mass m ~ 10^-10   eV/c² (Mielke & Schunck 2000). Some authors (Baldeschi et al. 1983; Sin 1994; Schunck 1998; Matos & Guzmán 1999; Hu et al. 2000; Arbey et al. 2001; Silverman & Mallett 2002) have also proposed that dark matter halos could be giant systems of “Bose liquid”, but in that case the mass of the bosons must be extremely low (m ~ 10^-24   eV/c²) to yield masses consistent with the mass of galactic halos. Such an ultralight scalar field is called “fuzzy cold dark matter” (FCDM) by Hu et al. (2000), who discuss its overall cosmological behavior. On the other hand, Colpi et al. (1986) have shown that if the bosons have a self-interaction, then the mass of the boson stars can increase considerably, even for a small self-interaction. For m ~ 1   GeV/c² it approaches the solar mass so that dark matter could be made of numerous boson stars. On the other hand, for m ~ 1   eV/c², the mass of the boson stars approaches the mass of the galactic halos. Therefore, some authors (Lee & Koh 1996; Peebles 2000; Goodman 2000; Arbey et al. 2003; Böhmer & Harko 2007; Fukuyama et al. 2008) have proposed that dark matter halos themselves could be in the form of gigantic self-gravitating BECs with short-range interactions described by a single wave function ψ(r,t). In the Newtonian limit, which is relevant on the galactic scale, the evolution of this wave function is governed by the Gross-Pitaevskii-Poisson (GPP) system (Böhmer & Harko 2007). With the Madelung (1927) transformation, the GP equation turns out to be equivalent to hydrodynamic (Euler) equations involving an isotropic classical pressure due to short-range interactions (scattering) and an anisotropic quantum pressure (or a quantum potential) arising from the Heisenberg uncertainty principle. For a standard BEC with quartic self-interaction, the equation of state is that of a polytrope with index n = 1. On large scales, scattering and quantum effects are negligible, and one recovers the classical hydrodynamic equations of cold dark matter (CDM) models that are remarkably successful in explaining the large-scale structure of the universe (Peebles & Ratra 2003). On small scales, gravitational collapse is prevented by the (repulsive) scattering or by the uncertainty principle. This may be a way to solve the problems of the CDM model, such as the cusp and the missing satellite problems (Hu et al. 2000).

To this day, the BEC model has two main virtues. (i) It can reproduce the rotation curves of several low surface brightness (LSB) galaxies (Sin 1994; Lee & Koh 1996; Schunck 1998; Guzmán & Matos 2000; Arbey et al. 2003; Böhmer & Harko 2007) and account for their almost flat behavior at long distances, and (ii) the quantum or scattering pressure can stabilize the system against gravitational collapse and lead to dark matter halos with flat central density profiles (Hu et al. 2000; Matos & Ureña-López 2001; Bernal et al. 2008). This is interesting because there is strong observational evidence that the distribution of dark matter presents a nearly constant density core (Burkert 1995). For example, in the case of LSB galaxies and dwarf spiral galaxies, the content of baryons is small and these objects are dominated by dark matter. The Hi nearby galaxy survey (THINGS) of dwarf galaxies (Walter et al. 2008) has revealed discrepancies between observations and simulations done in the framework of the ΛCDM model. In particular, the logarithmic inner slopes of the density profiles are on the order  −0.2 (Oh et al. 2011), as for LSB galaxies (de Blok et al. 2003), instead of the  −1 predicted by CDM simulations (Navarro et al. 1996). Therefore, the observations seem to favor flat density profiles over cuspy profiles. This is consistent with the BEC model, as shown in the recent paper of Harko (2011b). At the same time, this explains why the abundance of satellite galaxies predicted by CDM simulations is not observed around the Milky Way (Klypin et al. 1999). It is difficult to presume the relevance of the BEC model to account for all astrophysical observations, but there is a growing accumulation of recent material indicating that this model may be viable. In any case, it deserves to be studied in greater detail at the physical and astrophysical levels because it combines fundamental elements of gravity, quantum mechanics, and fluid mechanics.

In our previous papers (Chavanis 2011a; Chavanis & Delfini 2011), we performed an exhaustive study of the equilibrium configurations of a Newtonian self-gravitating BEC with short-range interactions. For a given value of the scattering length, we obtained the mass-radius relation M(R) connecting the noninteracting case (corresponding to low masses) studied by Ruffini & Bonazzola (1969) to the Thomas-Fermi (TF) limit (corresponding to large masses) investigated by Böhmer & Harko (2007). We also considered the case of attractive self-interaction. This corresponds to a negative scattering length (a_s < 0) yielding a negative pressure. In that case, we found there is a maximum mass ${\mathit{M}}_{\max }\mathrm{=}\mathrm{1.012}ħ\mathit{/}\sqrt{\mathrm{\left|}{\mathit{a}}_{\mathrm{s}}\mathrm{\right|}\mathit{Gm}}$ above which the system becomes unstable. It turns out that this mass is ridiculously low (being possibly as low as the Planck mass M_P = 2.18   10^-8   kg ! ), meaning that a self-gravitating BEC with attractive short-range interactions is extremely unstable (except if the mass of the bosons is extraordinarily low, on the order of m ~ 10^-22   eV/c², like in Hu et al. 2000). We proposed that the scattering length could be negative in the early universe and that it could help form structures¹. Then, the scattering length could become positive and help stabilize the structures against gravitational collapse. Of course, the mechanism by which the scattering length changes sign (due, for example, to a change in density) remains to be established so that this idea is highly speculative. We note that some atoms in terrestrial BEC experiments are reported to have negative scattering lengths (Dalfovo et al. 1999). Therefore, the possibility of negative pressure for a BEC can be contemplated. Furthermore, it has been experimentally demonstrated that, under certain conditions, it is possible to manipulate the sign and value of the scattering length (Fedichev et al. 1996). The extension of these ideas to cosmology remains an important problem that we shall not discuss further here. We assume that the dark matter in the universe is a BEC with repulsive (a_s > 0) or attractive (a_s < 0) self-interaction and theoretically explore the consequences of this hypothesis.

If dark matter halos are BECs, they have probably formed by Jeans instability. The gravitational instability of a scalar field equivalent to a BEC was considered by Khlopov et al. (1985); Bianchi et al. (1990); Hu et al. (2000); and Sikivie & Yang (2009). However, these authors all started their analysis from relativistic field equations and did not take the self-interaction of the particles into account. In our previous paper (Chavanis 2011a), we studied the Jeans instability of a Newtonian self-gravitating BEC with short-range interactions described by the Gross-Pitaevskii-Poisson system in a static universe. We considered both attractive and repulsive self-interaction and found that when the scattering length is negative the growth rate of the instability increases with respect to the case where the scattering length is positive or zero. The next step is to study the Jeans instability of a BEC in an expanding universe. This is the topic of the present paper.

In the first part of the paper (Sects. 2–7), we neglect relativistic effects and use the equations of Newtonian cosmology introduced by Milne (1934) and McCrea & Milne (1934). We study the linear development of perturbations of a self-gravitating BEC with short-range interactions in an expanding Einstein-de Sitter (EdS) universe. In the TF approximation, the equation of state of a BEC with quartic self-interaction is that of a polytrope of index γ = 2 (n = 1). In the noninteracting case, we find that the BEC behaves similarly to a polytrope of index γ = 5/3 (n = 3/2) but with a quantum mechanically modified Jeans length. In the two cases, the equation for the density contrast can be analytically solved in terms of Bessel functions.

In the second part of the paper (Sect. 8), we take special relativistic effects into account and use the equations of Newtonian cosmology with pressure introduced by McCrea (1951). In that case, the evolution of the cosmic fluid depends on the equation of state. We consider BEC dark matter with equation of state p = kρ²c² and add the contribution of radiation, baryons, and dark energy (cosmological constant). For BEC dark matter with repulsive self-interaction k > 0 (positive pressure), the scale factor increases more rapidly than in the standard ΛCDM model where dark matter is pressureless (k = 0), while it increases less rapidly for BEC dark matter with attractive self-interaction k < 0 (negative pressure). We study the linear development of the perturbations in these two cases and show that the perturbations grow faster in BEC dark matter than in pressureless dark matter.

In the third part of the paper (Sect. 9), we consider a “dark fluid” with a generalized equation of state p = (αρ + kρ²)c² having a component p = kρ²c² similar to BEC dark matter and a component p = αρc² mimicking the effect of the cosmological constant (dark energy). We find optimal parameters that give good agreement with the standard ΛCDM model. Therefore, this “mixed” equation of state could be an alternative to the cosmological constant.

While our series of papers on this subject was being written, Harko (2011a) published a paper where he also considered the formation of structures in an expanding universe made of BEC dark matter. Our independent contribution is complementary to Harko’s work and confirms his main results. In addition, we address the following issues: (i) we study the effect of the quantum pressure in the first part of the paper; (ii) we use different relativistic hydrodynamic equations to model the cosmic fluid in the second part of the paper; (iii) we study a generalized equation of state in the third part of the paper, and (iv) we consider both positive and negative scattering lengths throughout the paper. The papers of Böhmer & Harko (2007) and Harko (2011a) show that a cosmic BEC can be an interesting model for the dark matter of the universe. Our series of papers, which develop this idea, go in the same direction.

Finally, in a different perspective, Widrow & Kaiser (1993) have proposed to describe a classical collisionless self-gravitating system by the Schrödinger-Poisson system. In this approach, the constant ħ is not the Planck constant, but rather an adjustable parameter that controls the spatial resolution λ_deB through a de Broglie relation λ_deB = ħ/mv. They argue that when ħ → 0, the Vlasov-Poisson system is recovered and that a finite value of ħ provides a small-scale regularization of the dynamics. In that case, the Schrödinger-Poisson system has nothing to do with quantum mechanics since it aims at describing the evolution of classical collisionless matter under the influence of gravity (in static or expanding universes). Still, from a mathematical point of view, these equations are equivalent to those describing self-gravitating BECs without self-interaction. Therefore, the results of the first part of our paper can have application in that context, independently of quantum mechanics.

2. The Gross-Pitaevskii-Poisson system

Following Böhmer & Harko (2007), we assume that dark matter is a BEC. A self-gravitating BEC with short-range interactions is described by the Gross-Pitaevskii-Poisson system $\begin{array}{ccc}\mathit{i}ħ\frac{\mathit{\partial \psi }}{\mathit{\partial t}}& \mathrm{=}& \mathrm{-}\frac{{ħ}^{\mathrm{2}}}{\mathrm{2}\mathit{m}}\mathrm{\Delta }\mathit{\psi }\mathrm{+}\mathit{m}\mathrm{(}\mathrm{\Phi }\mathrm{+}\mathit{h}\mathrm{(}\mathit{\rho }\mathrm{\left)}\mathrm{\right)}\mathit{\psi ,}\\ \mathrm{\Delta \Phi }& \mathrm{=}& \mathrm{4}\mathit{\pi GNm}\mathrm{|}\mathit{\psi }{\mathrm{|}}^{\mathrm{2}}\mathrm{-}\mathrm{\Lambda }\mathit{,}\end{array}$where ρ = Nm|ψ|² is the density, Φ the gravitational potential, and $\mathit{h}\mathrm{\left(}\mathit{\rho }\mathrm{\right)}\mathrm{=}{}^{\mathrm{\int }}\mathit{u}_{\mathrm{SR}}\mathrm{(}{r}\mathrm{-}{{r}}^{\mathrm{\prime }}\mathrm{)}\mathit{\rho }\mathrm{\left(}{{r}}^{\mathrm{\prime }}\mathit{,t}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{{r}}^{\mathrm{\prime }}$ an effective potential that takes small-scale interactions into account. For the sake of generality, we have included the cosmological constant Λ in the Poisson equation². We write the wave function in the form ψ(r,t) = A(r,t)e^iS(r,t)/ħ where A and S are real, and make the Madelung (1927) transformation $\mathit{\rho }\mathrm{=}\mathit{Nm}\mathrm{|}\mathit{\psi }{\mathrm{|}}^{\mathrm{2}}\mathrm{=}\mathit{Nm}{\mathit{A}}^{\mathrm{2}}\mathit{,} {u}\mathrm{=}\frac{\mathrm{1}}{\mathit{m}}\mathrm{\nabla }\mathit{S,}$(3)where ρ(r,t) is the density field and u(r,t) the velocity field. We note that the flow is irrotational since ∇ × u = 0. With this transformation, it can be shown that the Gross-Pitaevskii Eq. (1) is equivalent to the barotropic Euler equations with an additional term Q called the quantum potential (or quantum pressure) arising from the Heisenberg uncertainty principle. Indeed, one obtains the set of equations $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{u}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }{u}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{u}\mathrm{·}\mathrm{\nabla }\mathrm{)}{u}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }}\mathrm{\nabla }\mathit{p}\mathrm{-}\mathrm{\nabla \Phi }\mathrm{-}\frac{\mathrm{1}}{\mathit{m}}\mathrm{\nabla }\mathit{Q,}\\ & & \mathrm{\Delta \Phi }\mathrm{=}\mathrm{4}\mathit{\pi G\rho }\mathrm{-}\mathrm{\Lambda }\mathit{,}\end{array}$with $\mathit{Q}\mathrm{=}\mathrm{-}\frac{{ħ}^{\mathrm{2}}}{\mathrm{2}\mathit{m}}\frac{\mathrm{\Delta }\sqrt{\mathit{\rho }}}{\sqrt{\mathit{\rho }}}\mathrm{·}$(7)The pressure p(ρ) is determined by the effective potential h(ρ), playing the role of an enthalpy (Chavanis 2011a), through the relation p′(ρ) = ρh′(ρ). This yields p(ρ) = ρh(ρ) − H(ρ), where H is a primitive of h. For a contact pair interaction u_SR(r − r′) = gδ(r − r′), the effective potential h(ρ) = gρ where g = 4πa_sħ²/m³ is the pseudo-potential and a_s is the s-scattering length (Dalfovo et al. 1999). For the sake of generality, we allow a_s to be positive or negative. The corresponding equation of state is $\mathit{p}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }{\mathit{a}}_{\mathrm{s}}{ħ}^{\mathrm{2}}}{{\mathit{m}}^{\mathrm{3}}}{\mathit{\rho }}^{\mathrm{2}}\mathrm{·}$(8)It corresponds to a polytropic equation of state of the form $\mathit{p}\mathrm{=}\mathit{K}{\mathit{\rho }}^{\mathit{\gamma }}\mathit{,} \mathit{\gamma }\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathit{n}}\mathit{,}$(9)with polytropic index n = 1 (i.e. γ = 2) and polytropic constant K = 2πa_sħ²/m³.

The effective potential corresponding to the polytropic equation of state (9) is $\mathit{h}\mathrm{\left(}\mathit{\rho }\mathrm{\right)}\mathrm{=}\frac{\mathit{K\gamma }}{\mathit{\gamma }\mathrm{-}\mathrm{1}}{\mathit{\rho }}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathit{.}$(10)This leads to a GP equation of the form $\mathrm{i}ħ\frac{\mathit{\partial \psi }}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\frac{{ħ}^{\mathrm{2}}}{\mathrm{2}\mathit{m}}\mathrm{\Delta }\mathit{\psi }\mathrm{+}\mathit{m}\mathrm{(}\mathrm{\Phi }\mathrm{+}\mathit{\kappa }\mathrm{|}\mathit{\psi }{\mathrm{|}}^{\mathrm{2}\mathit{/}\mathit{n}}\mathrm{)}\mathit{\psi ,}$(11)where κ = K(n + 1)(Nm)^1/n. This is the usual form of the GP equation considered in the literature (Sulem & Sulem 1999). The standard BEC, with a quartic interaction p ∝ ρ² ∝ |ψ|⁴, corresponds to n = 1. We may also recall that classical and ultra relativistic fermion stars are equivalent to polytropes with index n = 3/2 and n = 3 respectively (Chandrasekhar 1939). The GPP system (11)–(2) with n = 3/2 has been studied by Bilic et al. (2001) in relation to the formation of white dwarf stars by gravitational collapse. They show that the quantum pressure (Heisenberg) can regularize the dynamics on small scales.

For an isothermal equation of state $\mathit{p}\mathrm{=}\mathit{\rho }\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{eff}}}{\mathit{m}}\mathit{,}$(12)the effective potential is $\mathit{h}\mathrm{\left(}\mathit{\rho }\mathrm{\right)}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{eff}}}{\mathit{m}}\ln \mathit{\rho ,}$(13)and the GP equation reads as $\mathrm{i}ħ\frac{\mathit{\partial \psi }}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\frac{{ħ}^{\mathrm{2}}}{\mathrm{2}\mathit{m}}\mathrm{\Delta }\mathit{\psi }\mathrm{+}\mathrm{(}\mathit{m}\mathrm{\Phi }\mathrm{+}\mathrm{2}{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{eff}}\ln \mathrm{|}\mathit{\psi }\mathrm{\left|}\mathrm{\right)}\mathit{\psi }\mathit{.}$(14)Interestingly, a nonlinear Schrödinger equation with a logarithmic potential similar to Eq. (14) was introduced long ago by Bialynicki & Mycielski (1976) as a possible generalization of the Schrödinger equation in quantum mechanics.

Remark: for a BEC at T = 0, the pressure arising in the Euler Eq. (5) has a meaning different from the kinetic pressure of a normal fluid at finite temperature. It is due to the self-interaction of the particles encapsulated in the effective potential h(ρ) and not to thermal motion (since T = 0). As a result, the pressure can be negative (!) contrary to the kinetic pressure. This is the case, in particular, for a BEC described by the equation of state (8) when the scattering length a_s is negative. In terrestrial BEC experiments, some atoms like ⁷Li have a negative scattering length (Fedichev et al. 1996). On the other hand, the constant T_eff appearing in Eq. (12) is just an “effective” temperature since it arises from a particular form of self-interaction and has nothing to do with the kinetic temperature (which here is T = 0). In particular, this effective temperature can be negative. In that respect, we note that linear equations of state p = αρc² with α < 0 have been introduced heuristically to account for the accelerated expansion of the universe (see the review of Peebles & Ratra 2003 and Sect. 9).

3. Newtonian cosmology

In this first part of the paper, we use the Newtonian cosmology introduced by Milne and McCrea & Milne in 1934. The great advantage of the Newtonian treatment is its simplicity³. Furthermore, the equations are identical to those derived using the theory of general relativity, provided the pressure is negligible in comparison with the energy density ρc² where c is the speed of light (we go beyond this limitation in the second part of the paper). This makes this approach attractive. Furthermore, its simplicity allows us to introduce novel ingredients such as the quantum pressure that is derived from the classical Gross-Pitaevskii equation. The validity and limitation of Newtonian cosmology have been discussed by various authors such as Layzer (1954), McCrea (1955), Callan et al. (1965) and Harrison (1965).

We recall now the basics of Newtonian cosmology. We consider a spatially homogeneous solution of Eqs. (4)–(7) of the form $\mathit{\rho }\mathrm{\left(}{r}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,} {u}\mathrm{\left(}{r}\mathit{,t}\mathrm{\right)}\mathrm{=}\frac{\mathit{ȧ}}{\mathit{a}}{r}\mathit{,}$(15)where a(t) is the scale factor and H = ȧ/a is the Hubble “constant” (in fact a function of time). The Euler Eqs. (4)–(5) reduce to $\begin{array}{ccc}& & \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{b}}}{\mathrm{d}\mathit{t}}\mathrm{+}\mathrm{3}{\mathit{\rho }}_{\mathrm{b}}\frac{\mathit{ȧ}}{\mathit{a}}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\begin{array}{}\mathrm{¨}\\ \mathit{a}\end{array}}{\mathit{a}}{r}\mathrm{=}\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{b}}\mathit{.}\end{array}$The pressure p and the quantum potential Q do not appear in the theory of the homogeneous model since they enter Eq. (5) only through their gradients. Equation (16) leads to the relation ${\mathit{\rho }}_{\mathrm{b}}{\mathit{a}}^{\mathrm{3}}\mathrm{~}\mathrm{1}\mathit{,}$(18)which corresponds to the conservation of mass. Taking the divergence of Eq. (17) and using the Poisson Eq. (6), we obtain the cosmological equation $\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{a}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\mathrm{=}\mathrm{-}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\mathit{a}\mathrm{+}\frac{\mathrm{\Lambda }}{\mathrm{3}}\mathit{a}\mathit{.}$(19)Using Eq. (18), its first integral is ${\left(\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}\right)}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{(}\mathrm{8}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\mathrm{+}\mathrm{\Lambda }\mathrm{)}{\mathit{a}}^{\mathrm{2}}\mathrm{-}\mathit{\kappa ,}$(20)where κ is a constant of integration. Equations (18)–(20) are the Newtonian equations of an isotropic and homogeneous universe. They coincide with the equations derived by Friedmann (1922, 1924) for Λ = 0 and κ =  ± 1 and by Einstein & de Sitter (1932) for Λ = 0 and κ = 0 from the theory of general relativity when the pressure is much lower than the energy density ρc². In that case, κ is the curvature constant and space is flat (κ = 0), elliptical (κ = 1), or hyperbolic (κ =  −1)⁴.

The Einstein (1917) static universe corresponds to $\mathit{a}\mathrm{=}\mathrm{1}\mathit{,} {\mathit{\rho }}_{\mathrm{b}}\mathrm{=}\frac{\mathrm{\Lambda }}{\mathrm{4}\mathit{\pi G}}\mathit{,} \mathit{\kappa }\mathrm{=}\mathrm{\Lambda }\mathit{.}$(21)However, this universe is unstable against perturbations in a (Eddington 1930, Harrison 1967). The Einstein-de Sitter (EdS) universe corresponds to Λ = 0 and κ = 0. In that case, Eqs. (19) and (20) reduce to $\begin{array}{}\mathrm{¨}\\ \mathit{a}\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\mathit{a,} {\mathit{ȧ}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{8}}{\mathrm{3}}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}{\mathit{a}}^{\mathrm{2}}\mathit{.}$(22)This yields $\mathit{a}\mathrm{\propto }{\mathit{t}}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathit{,} \mathit{H}\mathrm{=}\frac{\mathit{ȧ}}{\mathit{a}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}\mathit{t}}\mathit{,} {\mathit{\rho }}_{\mathrm{b}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{6}\mathit{\pi G}{\mathit{t}}^{\mathrm{2}}}\mathrm{·}$(23)Both inflationary theory (Guth 1981) and observations favor a flat universe (κ = 0). In that case, Eq. (20) can be written ${\mathit{H}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{8}}{\mathrm{3}}\mathit{\pi G}\mathrm{(}{\mathit{\rho }}_{\mathrm{b}}\mathrm{+}{\mathit{\rho }}_{\mathrm{\Lambda }}\mathrm{)}$, where ρ_Λ = Λ/8πG is the dark energy density, and it gives a relationship between Hubble’s constant and the total density of the universe. With the present-day value of the Hubble constant H₀ = 2.273 × 10^-18   s^-1, the Einstein-de Sitter model (κ = Λ = 0, p ≪ ρc²) leads to an age of the universe t₀ = 2/3H₀ ~ 9.3 billion years, while the ΛCDM model, which is the standard model of cosmology, predicts 13.75 billion years. The effect of dark energy is considered in Sects. 8.4 and 9.

Remark: we can obtain Eqs. (18)–(20) directly if we use the “naive” picture that the universe is a uniform sphere of radius a(t), density ρ_b(t), and mass M. The conservation of mass $\mathit{M}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{\rho }}_{\mathrm{b}}{\mathit{a}}^{\mathrm{3}}$ leads to Eq. (18). On the other hand, applying the Gauss theorem at the border of the sphere, Newton’s equation can be written d²a/dt² =  −GM/a² + Λa/3 leading to Eq. (19). This is the equation of motion of a fictive particle of position a in a potential V(a) =  −GM/a − Λa²/6. Its first integral is given by Eq. (20) where E =  −κ/2 can be regarded as the energy of the fictive particle. The case E = 0 (EdS universe) corresponds to the situation where the velocity of the fictive particle is equal to the escape velocity. In this naive picture, a may be interpreted as the “radius” of the universe. In fact, a better derivation (which does not assume a finite sphere) is to apply Newton’s equation to an arbitrary fluid particle located in r and write d²r/dt² =  −G(4πρ_br³/3)/r² + Λr/3. Setting r = a(t)x and dividing by x yields Eq. (19).

4. Quantum Euler-Poisson system in an expanding universe

We now study the instability of the homogeneous background and the growth of perturbations that ultimately lead to the large-scale structures of the universe. We work in the comoving frame (Peebles 1980). To that purpose, we set ${r}\mathrm{=}\mathit{a}\mathrm{\left(}\mathit{t}\mathrm{\right)}{x}\mathit{,} {u}\mathrm{=}\frac{\mathit{ȧ}}{\mathit{a}}{r}\mathrm{+}{v}\mathit{,}$(24)where v is the peculiar velocity. For the moment, we allow arbitrary deviations from the background flow. We first write the Poisson Eq. (6) in the comoving frame. Integrating Eq. (17) and using Eq. (19), we find that the background gravitational potential is ${\mathrm{\Phi }}_{\mathrm{b}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{\begin{array}{}\mathrm{¨}\\ \mathit{a}\end{array}}{\mathit{a}}{\mathit{r}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\begin{array}{}\mathrm{¨}\\ \mathit{a}\end{array}\mathit{a}{\mathit{x}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathit{r}}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{\Lambda }}{\mathrm{6}}{\mathit{r}}^{\mathrm{2}}\mathit{.}$(25)This result can also be obtained by integrating the Poisson equation for a homogeneous system or by using the Gauss theorem. If we introduce the new potential φ = Φ − Φ_b, i.e. $\mathit{\phi }\mathrm{=}\mathrm{\Phi }\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\begin{array}{}\mathrm{¨}\\ \mathit{a}\end{array}\mathit{a}{\mathit{x}}^{\mathrm{2}}\mathit{,}$(26)the Poisson Eq. (6) becomes $\mathrm{\Delta }\mathit{\phi }\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{a}}^{\mathrm{2}}\mathrm{(}\mathit{\rho }\mathrm{-}{\mathit{\rho }}_{\mathrm{b}}\mathrm{)}\mathit{,}$(27)where the Laplacian is taken with respect to x. Now, following Peebles (1980) and taking the quantum pressure into account, we find that the hydrodynamic Eqs. (4) and (5) can be written in the comoving frame as $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{3}\mathit{ȧ}}{\mathit{a}}\mathit{\rho }\mathrm{+}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{a}}\mathrm{(}{v}\mathrm{·}\mathrm{\nabla }\mathrm{)}{v}\mathrm{+}\frac{\mathit{ȧ}}{\mathit{a}}{v}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho a}}\mathrm{\nabla }\mathit{p}\mathrm{-}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathit{\phi }\\ & & \mathrm{+}\frac{{ħ}^{\mathrm{2}}}{\mathrm{2}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{3}}}\mathrm{\nabla }\left(\frac{\mathrm{\Delta }\sqrt{\mathit{\rho }}}{\sqrt{\mathit{\rho }}}\right)\mathrm{·}\end{array}$It is convenient to write the density in the form $\begin{array}{ccc}\mathit{\rho }\mathrm{=}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\left[\mathrm{1}\mathrm{+}\mathit{\delta }\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\right]\mathit{,}& & \end{array}$(30)where ρ_b ∝ 1/a³ and δ(x,t) is the density contrast (Peebles 1980). Substituting Eq. (30) into Eqs.  (27)–(29), we obtain the quantum barotropic Euler-Poisson system in an expanding universe $\begin{array}{ccc}& & \frac{\mathit{\partial \delta }}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{\delta }\mathrm{\left)}{v}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{a}}\mathrm{(}{v}\mathrm{·}\mathrm{\nabla }\mathrm{)}{v}\mathrm{+}\frac{\mathit{ȧ}}{\mathit{a}}{v}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho a}}\mathrm{\nabla }\mathit{p}\mathrm{-}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathit{\phi }\\ & & \mathrm{+}\frac{{ħ}^{\mathrm{2}}}{\mathrm{2}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{3}}}\mathrm{\nabla }\left(\frac{\mathrm{\Delta }\sqrt{\mathrm{1}\mathrm{+}\mathit{\delta }}}{\sqrt{\mathrm{1}\mathrm{+}\mathit{\delta }}}\right)\mathit{,}\\ & & \mathrm{\Delta }\mathit{\phi }\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}{\mathit{a}}^{\mathrm{2}}\mathit{\delta }\mathit{.}\end{array}$For ħ = p = 0, we recover the usual Euler equations of a cold gas (T = 0) in an expanding universe (Peebles 1980)⁵. For a barotropic equation of state p = p(ρ), the Euler Eq. (32) can be rewritten $\begin{array}{ccc}\frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{a}}\mathrm{(}{v}\mathrm{·}\mathrm{\nabla }\mathrm{)}{v}\mathrm{+}\frac{\mathit{ȧ}}{\mathit{a}}{v}\mathrm{=}\mathrm{-}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\delta }\mathrm{)}\mathit{a}}\mathrm{\nabla }\mathit{\delta }\mathrm{-}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathit{\phi }& & \\ \mathrm{+}\frac{{ħ}^{\mathrm{2}}}{\mathrm{2}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{3}}}\mathrm{\nabla }\left(\frac{\mathrm{\Delta }\sqrt{\mathrm{1}\mathrm{+}\mathit{\delta }}}{\sqrt{\mathrm{1}\mathrm{+}\mathit{\delta }}}\right)\mathit{,}& & \end{array}$(34)where ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\rho }\mathrm{\right)}\mathrm{=}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{\delta }\mathrm{\left)}\mathrm{\right)}$ is the square of the velocity of sound in the evolving system. It generically depends on position and time. For an isothermal equation of state, ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathit{/}\mathit{m}$ is a constant and for a polytropic equation of state ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathit{K\gamma }{\mathit{\rho }}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}$. For a standard BEC with a quartic self-interaction, described by the equation of state (8), the square of the velocity of sound is $\begin{array}{ccc}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{s}}{ħ}^{\mathrm{2}}\mathit{\rho }}{{\mathit{m}}^{\mathrm{3}}}\mathit{,}& & \end{array}$(35)and we obtain the system of equations $\begin{array}{ccc}& & \frac{\mathit{\partial \delta }}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{\delta }\mathrm{\left)}{v}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{a}}\mathrm{(}{v}\mathrm{·}\mathrm{\nabla }\mathrm{)}{v}\mathrm{+}\frac{\mathit{ȧ}}{\mathit{a}}{v}\mathrm{=}\mathrm{-}\frac{\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{s}}{ħ}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{b}}}{{\mathit{m}}^{\mathrm{3}}\mathit{a}}\mathrm{\nabla }\mathit{\delta }\\ & & \mathrm{-}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathit{\phi }\mathrm{+}\frac{{ħ}^{\mathrm{2}}}{\mathrm{2}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{3}}}\mathrm{\nabla }\left(\frac{\mathrm{\Delta }\sqrt{\mathrm{1}\mathrm{+}\mathit{\delta }}}{\sqrt{\mathrm{1}\mathrm{+}\mathit{\delta }}}\right)\mathit{,}\\ & & \mathrm{\Delta }\mathit{\phi }\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}{\mathit{a}}^{\mathrm{2}}\mathit{\delta }\mathit{.}\end{array}$

5. Linearized equations

Of course, δ = φ = 0 and v = 0 is a solution of the quantum Euler-Poisson system (31)–(33) in the comoving frame, corresponding to the pure background flow (15). However, this solution may be unstable to some perturbations. If we consider small perturbations δ ≪ 1, φ ≪ 1, |v| ≪ 1 and linearize the foregoing equations we obtain $\begin{array}{ccc}& & \frac{\mathit{\partial \delta }}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathrm{·}{v}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{ȧ}}{\mathit{a}}{v}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{a}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\nabla }\mathit{\delta }\mathrm{-}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathit{\phi }\mathrm{+}\frac{{ħ}^{\mathrm{2}}}{\mathrm{4}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{3}}}\mathrm{\nabla }\mathrm{\left(}\mathrm{\Delta }\mathit{\delta }\mathrm{\right)}\mathit{,}\\ & & \mathrm{\Delta }\mathit{\phi }\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\mathit{\delta }{\mathit{a}}^{\mathrm{2}}\mathit{,}\end{array}$where ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}$ now denotes the square of the velocity of sound in the homogeneous background flow. It is just a function of time. Taking the time derivative of Eq. (39) multiplied by a, the divergence of Eq. (40), and using Eq. (41), these equations can be combined into a single equation governing the evolution of the density contrast $\begin{array}{ccc}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\delta }}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{+}\mathrm{2}\frac{\mathit{ȧ}}{\mathit{a}}\frac{\mathit{\partial \delta }}{\mathit{\partial t}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{\Delta }\mathit{\delta }\mathrm{+}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\mathit{\delta }\mathrm{-}\frac{{ħ}^{\mathrm{2}}}{\mathrm{4}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{4}}}{\mathrm{\Delta }}^{\mathrm{2}}\mathit{\delta }\mathit{.}& & \end{array}$(42)For ħ = 0, we recover the equation first derived by Bonnor (1957). In his case, the pressure p is a kinetic pressure. Expanding the solution in Fourier modes of the form δ(x,t) = δ_k(t)e^ik·x, we obtain $\begin{array}{ccc}\begin{array}{}\mathrm{¨}\\ \mathit{\delta }\end{array}\mathrm{+}\mathrm{2}\frac{\mathit{ȧ}}{\mathit{a}}\mathit{\delta ̇}\mathrm{+}\left(\frac{{ħ}^{\mathrm{2}}{\mathit{k}}^{\mathrm{4}}}{\mathrm{4}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{4}}}\mathrm{+}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{-}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\right)\mathit{\delta }\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(43)where, for brevity, we have noted δ(t) for δ_k(t). For a standard BEC, using ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{s}}{ħ}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{b}}\mathit{/}{\mathit{m}}^{\mathrm{3}}$, the foregoing equation can be rewritten $\begin{array}{ccc}\begin{array}{}\mathrm{¨}\\ \mathit{\delta }\end{array}\mathrm{+}\mathrm{2}\frac{\mathit{ȧ}}{\mathit{a}}\mathit{\delta ̇}\mathrm{+}\left(\frac{{ħ}^{\mathrm{2}}{\mathit{k}}^{\mathrm{4}}}{\mathrm{4}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{4}}}\mathrm{+}\frac{\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{s}}{ħ}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{b}}{\mathit{k}}^{\mathrm{2}}}{{\mathit{m}}^{\mathrm{3}}{\mathit{a}}^{\mathrm{2}}}\mathrm{-}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\right)\mathit{\delta }\mathrm{=}\mathrm{0.}& & \end{array}$(44)

6. Quantum Jeans length

6.1. Expanding universe

In the noninteracting case a_s = c_s = 0, the equation for the density contrast reduces to $\begin{array}{ccc}\begin{array}{}\mathrm{¨}\\ \mathit{\delta }\end{array}\mathrm{+}\mathrm{2}\frac{\mathit{ȧ}}{\mathit{a}}\mathit{\delta ̇}\mathrm{+}\left(\frac{{ħ}^{\mathrm{2}}{\mathit{k}}^{\mathrm{4}}}{\mathrm{4}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{4}}}\mathrm{-}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\right)\mathit{\delta }\mathrm{=}\mathrm{0.}& & \end{array}$(45)From this relation, we can define a time-dependent quantum Jeans wavenumber $\begin{array}{ccc}{\mathit{k}}_{\mathit{Q}}\mathrm{=}{\left(\frac{\mathrm{16}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{4}}}{{ħ}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{.}& & \end{array}$(46)The proper quantum Jeans wavenumber obtained by writing δ ∝ e^{ik _∗ ·r} is ${\mathit{k}}_{\mathit{Q}}^{\mathrm{\ast }}\mathrm{=}\mathrm{(}\mathrm{16}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}{\mathit{m}}^{\mathrm{2}}\mathit{/}{ħ}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{4}}$, but the comoving quantum Jeans wavenumber k_Q is better suited to our problem. Recalling Eq. (18), we can write k_Q = κ_Qa^1/4 where κ_Q = (16πGρ_ba³m²/ħ²)^1/4 is a constant. The quantum Jeans length λ_Q = 2π/k_Q decreases with time like a ^− 1/4 so that, in the comoving frame, the system becomes unstable on smaller and smaller scales as the universe expands (in contrast, the proper quantum Jeans scale increases with time as a^3/4). In the cold, nonquantum, universe (a_s = c_s = ħ = 0), there is no Jeans length: all the scales are unstable. Therefore, quantum effects can stabilize the system on small scales, on the order of λ_Q, and avoid density cusps (Hu et al. 2000).

In the TF limit where the quantum potential can be neglected, the equation for the density contrast reduces to $\begin{array}{ccc}\begin{array}{}\mathrm{¨}\\ \mathit{\delta }\end{array}\mathrm{+}\mathrm{2}\frac{\mathit{ȧ}}{\mathit{a}}\mathit{\delta ̇}\mathrm{+}\left(\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{-}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\right)\mathit{\delta }\mathrm{=}\mathrm{0.}& & \end{array}$(47)From this relation, we can define a time-dependent classical Jeans wavenumber $\begin{array}{ccc}{\mathit{k}}_{\mathrm{J}}\mathrm{=}{\left(\frac{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}{\mathit{a}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}& & \end{array}$(48)The proper Jeans wavenumber is ${\mathit{k}}_{\mathrm{J}}^{\mathrm{\ast }}\mathrm{=}\mathrm{(}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\mathit{/}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. The evolution of k_J(t) depends on the equation of state. For an isothermal equation of state for which ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathit{/}\mathit{m}$, recalling Eq. (18), we can write k_J = κ_Ja ^− 1/2 where κ_J = (4πGρ_ba³m/k_BT)^1/2 is a constant. The classical Jeans length λ_J = 2π/k_J increases with time like a^1/2 so that, in the comoving frame, the system becomes unstable on larger and larger scales as the universe expands. For a polytropic equation of state for which ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathit{K\gamma }{\mathit{\rho }}_{\mathrm{b}}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}$, we can write k_J = κ_Ja^{(3γ − 4)/2} where κ_J =  [4πG(ρ_ba³)^2 − γ/Kγ] ^1/2 is a constant. The classical Jeans wavelength behaves like λ_J ∝ a^{(4 − 3γ)/2}. For γ < 4/3, the Jeans wavelength grows with time and for γ > 4/3, it decreases with time. For γ = 4/3, the Jeans wavenumber is constant in time, i.e. k_J = κ_J. Finally, for a standard BEC, ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{s}}{ħ}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{b}}\mathit{/}{\mathit{m}}^{\mathrm{3}}$, the Jeans wavenumber can be written $\begin{array}{ccc}{\mathit{k}}_{\mathrm{J}}\mathrm{=}{\left(\frac{\mathit{G}{\mathit{m}}^{\mathrm{3}}{\mathit{a}}^{\mathrm{2}}}{{\mathit{a}}_{\mathrm{s}}{ħ}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}& & \end{array}$(49)It is independent of the density. Since k_J = κ_Ja with κ_J = (Gm³/a_sħ²)^1/2, the Jeans length decreases like λ_J ∝ 1/a.

Remark: the study of the Jeans instability in an expanding universe exhibits a critical value of the polytropic index γ_crit = 4/3 (i.e. n_crit = 3). It is interesting to note that the same critical index arises when one studies the dynamical stability of spatially inhomogeneous polytropic spheres with respect to the barotropic Euler-Poisson system. It has been established that a polytropic star is stable if γ ≥ 4/3 and unstable otherwise (Binney & Tremaine 1987). Surprisingly, the same index appears in the stability analysis of a spatially homogeneous polytropic gas in an expanding universe. The static study of Jeans (see the following section) does not directly exhibit a critical value of the polytropic index⁶.

6.2. Static universe

If we assume that Eq. (43) remains valid in a static universe, and take a = 1, we get $\begin{array}{ccc}\begin{array}{}\mathrm{¨}\\ \mathit{\delta }\end{array}\mathrm{+}\left(\frac{{ħ}^{\mathrm{2}}{\mathit{k}}^{\mathrm{4}}}{\mathrm{4}{\mathit{m}}^{\mathrm{2}}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{-}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\right)\mathit{\delta }\mathrm{=}\mathrm{0.}& & \end{array}$(50)Writing the perturbation in the form δ(t) ∝ e ^− iωt, we obtain the dispersion relation $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\frac{{ħ}^{\mathrm{2}}{\mathit{k}}^{\mathrm{4}}}{\mathrm{4}{\mathit{m}}^{\mathrm{2}}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{-}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\mathit{.}& & \end{array}$(51)This dispersion relation has been studied in our previous paper (Chavanis 2011a). The generalized Jeans wavenumber is given by ${\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}{\mathit{m}}^{\mathrm{2}}}{{ħ}^{\mathrm{2}}}\left[\sqrt{{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}\mathrm{+}\frac{\mathrm{4}\mathit{\pi G}{ħ}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{b}}}{{\mathit{m}}^{\mathrm{2}}}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\right]\mathit{.}$(52)In the noninteracting case (a_s = c_s = 0), we recover the quantum Jeans wavenumber (46) with a = 1, and in the TF limit, we recover the classical Jeans wavenumber (48) with a = 1. The dispersion relation can be written ${\mathit{\omega }}^{\mathrm{2}}\mathit{/}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}\mathrm{=}{\mathit{k}}^{\mathrm{4}}\mathit{/}{\mathit{k}}_{\mathit{Q}}^{\mathrm{4}}\mathrm{+}{\mathit{k}}^{\mathrm{2}}\mathit{/}{\mathit{k}}_{\mathrm{J}}^{\mathrm{2}}\mathrm{-}\mathrm{1}$, and the generalized Jeans wavenumber ${\mathit{k}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{k}}_{\mathit{Q}}^{\mathrm{4}}\mathit{/}\mathrm{2}{\mathit{k}}_{\mathrm{J}}^{\mathrm{2}}\mathrm{)}\mathrm{[}\mathrm{\pm }\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{4}{\mathit{k}}_{\mathrm{J}}^{\mathrm{4}}\mathit{/}{\mathit{k}}_{\mathit{Q}}^{\mathrm{4}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{]}$ with  +  when ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\ge }\mathrm{0}$ and  −  when ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{<}\mathrm{0}$. For λ < λ_c, ω is real and the perturbation oscillates with a pulsation ω; for λ > λ_c, ω is purely imaginary and the perturbation grows exponentially rapidly with a growth rate $\mathit{\gamma }\mathrm{=}\sqrt{\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}}$. The Jeans instability criterion is therefore λ > λ_c (i.e. k < k_c). For ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\ge }\mathrm{0}$, the maximum growth rate corresponds to k = 0 (infinite wavelengths) and is given by $\mathit{\gamma }\mathrm{=}\sqrt{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}}$. On the other hand, considering a BEC with negative scattering length, for which ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{<}\mathrm{0}$, it is found (Chavanis 2011a) that the maximum growth rate corresponds to k _∗  = (8π|a_s|ρ_b/m)^1/2 and is given by $\begin{array}{ccc}{\mathit{\gamma }}_{\mathrm{\ast }}\mathrm{=}\sqrt{\frac{\mathrm{16}{\mathit{\pi }}^{\mathrm{2}}{\mathit{a}}_{\mathrm{s}}^{\mathrm{2}}{ħ}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{b}}^{\mathrm{2}}}{{\mathit{m}}^{\mathrm{4}}}\mathrm{+}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}}\mathit{.}& & \end{array}$(53)Note that ${\mathit{k}}_{\mathrm{\ast }}\mathrm{=}\mathrm{(}{\mathit{k}}_{\mathit{Q}}^{\mathrm{4}}\mathit{/}\mathrm{2}\mathrm{\left|}{\mathit{k}}_{\mathrm{J}}^{\mathrm{2}}\mathrm{\right|}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ and ${\mathit{\gamma }}_{\mathrm{\ast }}\mathrm{=}\sqrt{\mathrm{4}\mathit{\pi G\rho }}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{k}}_{\mathit{Q}}^{\mathrm{4}}\mathit{/}\mathrm{4}{\mathit{k}}_{\mathrm{J}}^{\mathrm{4}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. For |a_s| ≫ (Gm⁴/4πħ²ρ_b)^1/2, we find that ${\mathit{\gamma }}_{\mathrm{\ast }}\mathrm{~}\mathrm{4}\mathit{\pi }\mathrm{\left|}{\mathit{a}}_{\mathrm{s}}\mathrm{\right|}ħ{\mathit{\rho }}_{\mathrm{b}}\mathit{/}{\mathit{m}}^{\mathrm{2}}\mathrm{\gg }\sqrt{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}}$. Therefore, an attractive short-range interaction (a_s < 0) increases the growth rate of the Jeans instability.

Jeans’ classical analysis (Jeans 1902, 1929) suffers from the defect that, in general, there is no initial stationary state that is in a uniform nonrotating fluid. Therefore, using Eq. (50) in a static universe is referred to as the “Jeans swindle”⁷. As noted by Bonnor (1957), the Jeans procedure is only valid in the static Einstein universe⁸. However, in this last case, the background density of the universe is ρ_b = Λ/4πG. The justification does not apply to a uniform mass of gas with a different density. Furthermore, the Einstein universe is strongly unstable. Therefore, this justification is not valid, and it is necessary to develop the Jeans instability analysis in an expanding universe.

Remark: Jeans analysis could be valid if the growth rate of the instability were much greater than the rate of the expansion of the universe. According to Eq. (22), the rate of expansion of the universe is H = ȧ/a = (8πGρ/3)^1/2. On the other hand, when ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{>}\mathrm{0}$, the maximum growth rate of the Jeans instability is $\mathit{\gamma }\mathrm{=}\sqrt{\mathrm{4}\mathit{\pi G\rho }}$. They are exactly of the same order of magnitude, which leads to the usual conclusion that the Jeans instability analysis is not valid (see, e.g., Weinberg 1972). However, when ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{<}\mathrm{0}$, the maximum growth rate given by Eq. (53) can be much higher. In that case, the expansion of the universe could be ignored and the static Jeans instability analysis could be valid.

7. Solution in the Einstein-de Sitter universe

7.1. The fundamental differential equation

The evolution of the perturbations is more complicated to analyze in an expanding universe than in a static universe, and it turns out to be very different. For simplicity, we consider an Einstein-de Sitter universe. In that case, it is possible to obtain analytical solutions of the linearized Eq. (43) for the density contrast in some particular cases. Measuring the evolution in terms of a instead of t and using Eq. (22), Eq. (43) is transformed into $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{\delta }}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{\delta }}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}{\mathit{a}}^{\mathrm{2}}}\left(\frac{{ħ}^{\mathrm{2}}{\mathit{k}}^{\mathrm{4}}}{\mathrm{16}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}{\mathit{m}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{4}}}\mathrm{+}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}{\mathit{a}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right)\mathit{\delta }\mathrm{=}\mathrm{0.}& & \end{array}$(54)For ħ = c_s = 0, we recover the case of a cold classical universe (dust) in which the pressure is zero and quantum effects are neglected. Equation (54) reduces to $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{\delta }}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{\delta }}{\mathrm{d}\mathit{a}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}{\mathit{a}}^{\mathrm{2}}}\mathit{\delta }\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(55)and its solutions are $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{+}}\mathrm{\propto }\mathit{a,} {\mathit{\delta }}_{\mathrm{-}}\mathrm{\propto }{\mathit{a}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathit{.}& & \end{array}$(56)The growth of perturbations is algebraic in an expanding universe, while the “naive” Jeans instability analysis in a static universe predicts an exponential growth. Therefore, the results of the two analysis are very different. As mentioned long ago by Bonnor (1957) and others, this (slow) algebraic growth may be a problem for forming large-scale structures rapidly enough.

Assuming a polytropic equation of state p = Kρ^γ for which ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathit{K\gamma }{\mathit{\rho }}_{\mathrm{b}}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}$, and introducing the variables of Sect. 6.1, we can rewrite Eq. (54) in the form $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{\delta }}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{\delta }}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}{\mathit{a}}^{\mathrm{2}}}\left(\frac{{\mathit{k}}^{\mathrm{4}}}{{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}\mathit{a}}\mathrm{+}\frac{{\mathit{k}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{3}\mathit{\gamma }\mathrm{-}\mathrm{4}}}\mathrm{-}\mathrm{1}\right)\mathit{\delta }\mathrm{=}\mathrm{0.} & & \end{array}$(57)Making the change in variables δ(a) = f(a)/a, we find that the differential equation for f is $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{f}}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\left(\frac{\mathrm{3}{\mathit{k}}^{\mathrm{4}}}{\mathrm{2}{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}\mathit{a}}\mathrm{+}\frac{\mathrm{3}{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{3}\mathit{\gamma }\mathrm{-}\mathrm{4}}}\mathrm{-}\mathrm{1}\right)\mathit{f}\mathrm{=}\mathrm{0.}& & \end{array}$(58)We have not been able to find the general solution of this equation in terms of simple functions, so we consider particular cases. Some of the solutions presented below have already been discussed in the literature (see, e.g., Harrison 1967), but we give more details and put emphasis on the solutions describing self-gravitating BECs that were not discussed previously.

7.2. The TF approximation

In the TF approximation where the quantum potential can be neglected, Eq. (58) reduces to $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{f}}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\left(\frac{\mathrm{3}{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{3}\mathit{\gamma }\mathrm{-}\mathrm{4}}}\mathrm{-}\mathrm{1}\right)\mathit{f}\mathrm{=}\mathrm{0.}& & \end{array}$(59)This is a particular case of the equation studied by Savedoff & Vila (1962) in relation to the work of Bonnor (1957)⁹. It can be solved in terms of Bessel functions (see Appendix A). For γ ≠ 4/3, f(a) is given by Eq. (A.7) and the evolution of the density contrast is $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}}{\mathit{J}}_{\mathrm{\pm }\frac{\mathrm{5}}{\mathrm{2}\mathrm{(}\mathrm{4}\mathrm{-}\mathrm{3}\mathit{\gamma }\mathrm{)}}}\left(\frac{\sqrt{\mathrm{6}}}{\mathrm{|}\mathrm{4}\mathrm{-}\mathrm{3}\mathit{\gamma }\mathrm{|}}\frac{\mathit{k}}{{\mathit{\kappa }}_{\mathrm{J}}}{\mathit{a}}^{\frac{\mathrm{4}\mathrm{-}\mathrm{3}\mathit{\gamma }}{\mathrm{2}}}\right)\mathit{.}& & \end{array}$(60)Using the asymptotic expansions of the Bessel functions (see Appendix A), we note that the density contrast is oscillating for k ≫ k_J(a) and growing for k ≪ k_J(a) (on the timescale over which these inequalities are fulfilled). Therefore, wavelengths that are longer than λ_J(a) behave like in a cold universe, whereas the perturbations have an oscillatory behavior for the short wavelengths. In this sense, the results are similar to those obtained with the classical Jeans analysis. They are, in fact, different because (i) the Jeans length varies with time; (ii) the oscillating solutions can be growing or decaying; and (iii) the growth is algebraic instead of exponential. The results depend more crucially on the value of the polytropic index γ than on the Jeans length λ_J. We must distinguish four cases. (i) γ < 4/3: for small a, the perturbations grow like in a cold gas (δ ₊  ~ a) and, for large a, they undergo damped oscillations (they decay like 1/a^{(5 − 3γ)/4}); (ii) 4/3 < γ < 5/3: for a → 0 the perturbations diverge like 1/a^{(5 − 3γ)/4} while oscillating, implying that we must start the stability analysis at a_i > 0. For small a ≥ a_i, the perturbations undergo damped oscillations (they decay like 1/a^{(5 − 3γ)/4}) and, for large a, they grow like in a cold gas (δ ₊  ~ a); (iii) γ = 5/3: for small a, the perturbations show pure oscillation and, for large a, they grow like in a cold gas (δ ₊  ~ a); (iv) γ > 5/3: for small a, the perturbations undergo growing oscillations (they grow like a^{(3γ − 5)/4}) and, for large a, they grow like in a cold gas (δ ₊  ~ a). Therefore, considering the small a regime, the perturbations grow for γ < 4/3, decay for 4/3 < γ < 5/3, and grow for γ > 5/3. Considering the large a regime, the perturbations decay for γ < 4/3 and grow for γ > 4/3. For γ < 4/3, the perturbations start to grow but finally decay. In a sense, the system is asymptotically stable. However, the growth of perturbations in the initial stage can trigger nonlinear effects that may induce instabilities.

For γ = 4/3 (n = 3), the density contrast behaves like (see Appendix A): $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }{\mathit{a}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{4}}\mathrm{\pm }\sqrt{\frac{\mathrm{25}}{\mathrm{16}}\mathrm{-}\frac{\mathrm{3}{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}}}}\mathit{.}& & \end{array}$(61)This polytropic index corresponds to a gas of photons at temperature T or to a gas of relativistic fermions at zero temperature. In that case, the Jeans length is independent of time: k_J = κ_J. For k > (25/24)^1/2k_J the density contrast behaves like $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}}\cos \left(\sqrt{\frac{\mathrm{3}{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{25}}{\mathrm{16}}}\ln \mathit{a}\mathrm{+}\mathit{\phi }\right)\mathit{.}& & \end{array}$(62)It diverges for a → 0, implying that we must start the stability analysis at a_i > 0. For a ≥ a_i, the perturbations decay like a ^− 1/4 while making oscillations. For k_J < k < (25/24)^1/2k_J, the perturbations decay algebraically without oscillating. For k < k_J, the perturbations grow algebraically (more precisely δ ₊  grows and δ ₋  decays). For k → 0, the perturbations behave like in a cold gas (δ ₊  grows like a and δ ₋  decays like a ^− 2/3). This situation is relatively close to the classical Jeans stability analysis, but (i) the growth of perturbation is algebraic instead of exponential; (ii) the oscillatory solutions decay with time; (iii) there is a short interval k_J < k < (25/24)^1/2k_J that has no counterpart in the static Jeans analysis.

We now consider particular equations of state with a physical meaning. The index γ = 1 corresponds to an isothermal universe. The density contrast is given by $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}}{\mathit{J}}_{\mathrm{\pm }\frac{\mathrm{5}}{\mathrm{2}}}\left(\sqrt{\mathrm{6}}\frac{\mathit{k}}{{\mathit{\kappa }}_{\mathrm{J}}}{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}\right)\mathit{.}& & \end{array}$(63)For the index γ = 5/3 (n = 3/2), we obtain $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}}{\mathit{J}}_{\mathrm{\pm }\frac{\mathrm{5}}{\mathrm{2}}}\left(\sqrt{\mathrm{6}}\frac{\mathit{k}}{{\mathit{\kappa }}_{\mathrm{J}}}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\right)\mathrm{·}& & \end{array}$(64)This polytropic index corresponds to a gas of nonrelativistic fermions at zero temperature. For the polytropic index γ = 2 (n = 1), we have $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}}{\mathit{J}}_{\mathrm{\pm }\frac{\mathrm{5}}{\mathrm{4}}}\left(\sqrt{\frac{\mathrm{3}}{\mathrm{2}}}\frac{\mathit{k}}{{\mathit{\kappa }}_{\mathrm{J}}}\frac{\mathrm{1}}{\mathit{a}}\right)\mathit{.}& & \end{array}$(65)This polytropic index corresponds to a BEC with quartic self-interaction in the TF approximation (see Sect. 2). In that case, κ_J = (Gm³/a_sħ²)^1/2. The previous expression assumes that a_s > 0. In the case where a_s < 0, we obtain $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}}{\mathit{I}}_{\mathrm{\pm }\frac{\mathrm{5}}{\mathrm{4}}}\left(\sqrt{\frac{\mathrm{3}}{\mathrm{2}}}\frac{\mathit{k}}{{\mathit{\kappa }}_{\mathrm{J}}}\frac{\mathrm{1}}{\mathit{a}}\right)\mathit{,}& & \end{array}$(66)with κ_J = (Gm³/|a_s|ħ²)^1/2. For large a, the perturbations behave like in a cold classical gas. For small a, we find that $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}{\mathrm{e}}^{\sqrt{\frac{\mathrm{3}}{\mathrm{2}}}\frac{\mathit{k}}{{\mathit{k}}_{\mathrm{J}}}\frac{\mathrm{1}}{\mathit{a}}}\mathrm{\to }\mathrm{+}\mathrm{\infty }\mathit{.}& & \end{array}$(67)We can combine the solutions δ ₊ (a) and δ ₋ (a) to avoid the divergence at the origin. In particular, the perturbation δ(a) ∝ δ ₊ (a) − δ ₋ (a) starts from zero at a = 0 with a horizontal tangent and rapidly increases. By contrast, in a cold classical gas, a perturbation δ(a) ≃ δ′(a) ≃ 0 at a = 0 takes a long time to grow. This demonstrates that the formation of structures is faster in BEC dark matter with an attractive self-interaction than in a cold classical universe.

Some curves displaying the evolution of the density contrast δ(a) in the different cases listed above are represented in Figs. 1–6 for illustration. The last figure shows the effect of negative scattering on the growth of perturbations.

Fig. 1 Evolution of the perturbation δ + (a) in an isothermal universe (γ = 1). This corresponds to the case γ < 4/3. We have taken k/κJ = 10.

Fig. 2 Evolution of the perturbation δ + (a) in a photonic or in a relativistic fermionic universe (γ = 4/3).

Fig. 3 Evolution of the perturbation δ + (a) in a polytropic universe with γ = 1.55. This corresponds to the case 4/3 < γ < 5/3. We have taken k/κJ = 0.637.

Fig. 4 Evolution of the perturbation δ + (a) in a nonrelativistic fermionic universe or in a BEC universe without self-interaction (γ = 5/3). We have taken k/κJ = 0.707.

Fig. 5 Evolution of the perturbation δ + (a) in a BEC universe with repulsive quartic self-interaction in the TF limit (γ = 2, as > 0). This corresponds to the case γ > 5/3. We have taken k/κJ = 0.5.

Fig. 6 Evolution of the perturbation δ(a) in a BEC universe with attractive quartic self-interaction in the TF limit (γ = 2, as < 0). We have numerically solved Eq. (57) with κQ →  + ∞ and κJ → iκJ, with the initial condition ai = 0.1, δ(ai) = 10-5 and δ′(ai) = 0. From bottom to top: k/κJ = 0.2,0.4,0.6,0.8,1. The lowest curve corresponds to a cold classical gas.

7.3. The noninteracting BEC case

In the noninteracting case a_s = c_s = 0, corresponding to a BEC universe without self-interaction, Eq. (58) reduces to $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{f}}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\left(\frac{\mathrm{3}{\mathit{k}}^{\mathrm{4}}}{\mathrm{2}{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}\mathit{a}}\mathrm{-}\mathrm{1}\right)\mathit{f}\mathrm{=}\mathrm{0.}& & \end{array}$(68)Comparing this equation with Eq. (59), we see that the dependence on a is the same as for a polytrope γ = 5/3. Therefore, f(a) is given by Eq. (A.7) with $\mathit{\lambda }\mathrm{=}\mathrm{3}{\mathit{k}}^{\mathrm{4}}\mathit{/}\mathrm{2}{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}$ and α =  −1/2 yielding $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}}{\mathit{J}}_{\mathrm{\pm }\frac{\mathrm{5}}{\mathrm{2}}}\left(\sqrt{\mathrm{6}}\frac{{\mathit{k}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{2}}}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\right)\mathit{.}& & \end{array}$(69)For small a, the perturbations oscillate, and, for large a, they behave like in a cold gas (δ ₊  grows like a). In a sense, a cold bosonic universe behaves similarly to a cold nonrelativistic fermionic universe (compare Eqs. (69) and (64)). However, the dependence on k is different (k² instead of k), and the expression of the Jeans length is also different (compare Eqs. (46) and (48)).

7.4. The nonrelativistic fermionic case γ = 5/3

A nonrelativistic gas of fermions at T = 0 is usually described by a polytropic equation of state p = Kρ^γ with index γ = 5/3 (n = 3/2) and polytropic constant K = (1/5)(3/8π)^2/3h²/m^8/3 (Chandrasekhar 1939). This equation of state arises from the Pauli exclusion principle. Such a semi-classical description is valid in the TF approximation, which is exact when N →  + ∞. In some cases, it can be relevant to go beyond the TF approximation and take the quantum pressure arising from the Heisenberg principle into account. This can regularize the dynamics on small scales, as shown by Bilic et al. (2001) in their study of the formation of fermion stars by gravitational collapse. In our problem, this general situation is described by Eq. (57) with γ = 5/3. This yields $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{\delta }}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{\delta }}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}{\mathit{a}}^{\mathrm{2}}}\left(\frac{{\mathit{k}}^{\mathrm{4}}}{{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}\mathit{a}}\mathrm{+}\frac{{\mathit{k}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}\mathit{a}}\mathrm{-}\mathrm{1}\right)\mathit{\delta }\mathrm{=}\mathrm{0.}& & \end{array}$(70)The equation for f(a) is $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{f}}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\left(\frac{\mathrm{3}{\mathit{k}}^{\mathrm{4}}}{\mathrm{2}{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}\mathit{a}}\mathrm{+}\frac{\mathrm{3}{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}\mathit{a}}\mathrm{-}\mathrm{1}\right)\mathit{f}\mathrm{=}\mathrm{0.}& & \end{array}$(71)As noted in Sect. 7.3, the two terms have the same dependence on a. Therefore, f(a) is given by Eq. (A.7) with $\begin{array}{ccc}\mathit{\lambda }\mathrm{=}\frac{\mathrm{3}{\mathit{k}}^{\mathrm{4}}}{\mathrm{2}{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}}\mathrm{+}\frac{\mathrm{3}{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}}\mathit{,} \mathit{\alpha }\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{·}& & \end{array}$(72)We obtain $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}}{\mathit{J}}_{\mathrm{\pm }\frac{\mathrm{5}}{\mathrm{2}}}\left(\sqrt{\frac{\mathrm{6}{\mathit{k}}^{\mathrm{4}}}{{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}}\mathrm{+}\frac{\mathrm{6}{\mathit{k}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}}}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\right)\mathrm{·}& & \end{array}$(73)For small a, the perturbations oscillate and, for large a, they behave like in a cold gas.

7.5. The relativistic fermionic case γ = 4/3

A relativistic gas of fermions at T = 0 is usually described by a polytropic equation of state p = Kρ^γ with index γ = 4/3 (n = 3) and polytropic constant K = (1/4)(3/8π)^1/3hc/m^4/3 (Chandrasekhar 1939). If we go beyond the TF approximation and take the quantum pressure (Heisenberg) into account, the evolution of the density contrast is described by Eq. (57) with γ = 4/3. This yields $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{\delta }}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{\delta }}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}{\mathit{a}}^{\mathrm{2}}}\left(\frac{{\mathit{k}}^{\mathrm{4}}}{{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}\mathit{a}}\mathrm{+}\frac{{\mathit{k}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right)\mathit{\delta }\mathrm{=}\mathrm{0.} & & \end{array}$(74)The equation for f(a) is $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{f}}{\mathrm{d}{\mathit{a}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{a}}\frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\left(\frac{\mathrm{3}{\mathit{k}}^{\mathrm{4}}}{\mathrm{2}{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{4}}\mathit{a}}\mathrm{+}\frac{\mathrm{3}{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right)\mathit{f}\mathrm{=}\mathrm{0.} & & \end{array}$(75)Its solution is given by Eq. (A.17), and the evolution of the density contrast is $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{\pm }}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{4}}}{\mathit{J}}_{\mathrm{\pm }\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\mathrm{25}\mathrm{-}\mathrm{24}\frac{{\mathit{k}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathrm{J}}^{\mathrm{2}}}}}\left(\sqrt{\mathrm{6}}\frac{{\mathit{k}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathit{Q}}^{\mathrm{2}}}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\right)\mathrm{·}& & \end{array}$(76)

We note that the classical Jeans scale appears in the index of the Bessel function, while the quantum Jeans scale appears in the argument of the Bessel function.

7.6. The general case (asymptotics)

We now treat the general case, taking quantum pressure and polytropic pressure (we assume ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{>}\mathrm{0}$) into account. Since we cannot solve Eq. (57) analytically, we just give asymptotic results.

We first consider the small a regime. (i) If γ < 5/3, the polytropic pressure is negligible in front of the quantum pressure and the perturbations oscillate (see Sect. 7.3). (ii) If γ = 5/3, the polytropic pressure and the quantum pressure are of the same order of magnitude, and the perturbations oscillate (see Sect. 7.4). (iii) If γ > 5/3, the quantum pressure is negligible in front of the polytropic pressure and the perturbations undergo growing oscillations δ ∝ a^{(3γ − 5)/4} (see Sect. 7.2). Therefore, the perturbations oscillate for γ ≤ 5/3 and grow for γ > 5/3. The quantum pressure plays a stabilizing role for γ ≤ 4/3 (see Sect. 7.2).

We now consider the large a regime. (i) If γ < 4/3, the quantum pressure is negligible in front of the polytropic pressure, and the perturbations undergo damped oscillations and decay like 1/a^{(5 − 3γ)/4} (see Sect. 7.2). (ii) If γ = 4/3, the quantum pressure is negligible in front of the polytropic pressure, and we are led back to the critical case of Sect. 7.2: for k > (25/24)^1/2k_J the perturbations decay like a ^− 1/4 while making oscillations, for k_J < k < (25/24)^1/2k_J the perturbations decay algebraically without oscillating and for k < k_J the perturbations grow algebraically. (iii) For γ > 4/3, the perturbations behave like in a cold gas and grow like δ ₊  ~ a. Therefore, the perturbations decay for γ < 4/3 and grow for γ > 4/3. In conclusion, the system is asymptotically stable for γ < 4/3 and asymptotically unstable for γ > 4/3. For the critical index γ = 4/3, the system is stable for k > k_J and unstable for k < k_J.

7.7. Formation of BEC dark matter halos

The results of this stability analysis lead to the following scenario. For a gas with polytropic index γ > 4/3, including the BEC dark matter model corresponding to γ = 2, an expanding homogeneous distribution of matter is gravitationally unstable to small perturbations and a process of fragmentation follows. In the linear regime, the perturbations grow and overdense regions appear. At first, these regions experience a cold collapse and intensify. Then, when pressure effects become important at high densities and small scales, localized clusters in virial equilibrium (DM halos) form. These clusters, which correspond to complete polytropic spheres, are stable according to the general stability criterion γ > 4/3 recalled at the end of Sect. 6.1. Therefore, this scenario explains the formation of BEC dark matter halos by Jeans instability in an expanding spatially homogeneous gas, and their stabilization by pressure effects when they become dense enough. The weakly nonlinear gravitational clustering (intermediate regime) could be described by resorting to the Zeldovich approximation, incorporating the properties of the BEC (see Chavanis 2011b).

Remark: an expanding homogeneous gas is stable if γ < 4/3 and unstable if γ > 4/3, leading to fragmentation and structure formation. Inversely, an isolated dense cluster (dark matter halo or star) is stable if γ > 4/3 and unstable if γ < 4/3. Finally, a collapsing homogeneous gas is stable if γ > 4/3 and unstable if γ < 4/3, leading to fragmentation (Chavanis 2011c).

8. Bose-Einstein condensate universe

8.1. Newtonian cosmology with pressure

We now take pressure effects (coming from special relativity) into account in the evolution of the cosmic fluid. We base our study on the set of hydrodynamic equations $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{u}\mathrm{\right)}\mathrm{+}\frac{\mathit{p}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{\nabla }\mathrm{·}{u}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }{u}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{u}\mathrm{·}\mathrm{\nabla }\mathrm{)}{u}\mathrm{=}\mathrm{-}\frac{\mathrm{\nabla }\mathit{p}}{\mathit{\rho }\mathrm{+}\frac{\mathit{p}}{{\mathit{c}}^{\mathrm{2}}}}\mathrm{-}\mathrm{\nabla \Phi }\mathit{,}\\ & & \mathrm{\Delta \Phi }\mathrm{=}\mathrm{4}\mathit{\pi G}\left(\mathit{\rho }\mathrm{+}\frac{\mathrm{3}\mathit{p}}{{\mathit{c}}^{\mathrm{2}}}\right)\mathrm{-}\mathrm{\Lambda }\mathit{,}\end{array}$where c is the velocity of light. For c →  + ∞, we recover the classical Euler-Poisson system (Binney & Tremaine 1987). These equations were introduced by McCrea (1951). His initial continuity equation, which contained an incorrect pressure gradient term, was corrected by Lima et al. (1997), and we considered this correction in writing Eq. (77). These equations lead to the correct relativistic equations for the cosmic evolution. Indeed, assuming a homogeneous and isotropic solution of the form ρ(r,t) = ρ_b(t), p(r,t) = p_b(t), and u(r,t) = (ȧ/a)r, Eqs. (77)–(79) reduce to $\begin{array}{ccc}& & \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{b}}}{\mathrm{d}\mathit{t}}\mathrm{+}\mathrm{3}\frac{\mathit{ȧ}}{\mathit{a}}\left({\mathit{\rho }}_{\mathrm{b}}\mathrm{+}\frac{{\mathit{p}}_{\mathrm{b}}}{{\mathit{c}}^{\mathrm{2}}}\right)\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\begin{array}{}\mathrm{¨}\\ \mathit{a}\end{array}}{\mathit{a}}\mathrm{=}\mathrm{-}\frac{\mathrm{4}\mathit{\pi G}}{\mathrm{3}}\left({\mathit{\rho }}_{\mathrm{b}}\mathrm{+}\frac{\mathrm{3}{\mathit{p}}_{\mathrm{b}}}{{\mathit{c}}^{\mathrm{2}}}\right)\mathrm{+}\frac{\mathrm{\Lambda }}{\mathrm{3}}\mathit{,}\\ & & {\left(\frac{\mathit{ȧ}}{\mathit{a}}\right)}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{8}\mathit{\pi G}}{\mathrm{3}}{\mathit{\rho }}_{\mathrm{b}}\mathrm{-}\frac{\mathit{\kappa }}{{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{\Lambda }}{\mathrm{3}}\mathit{.}\end{array}$These are precisely the Friedmann equations that can be derived from the theory of general relativity (Weinberg 1972). The perturbation theory based on Eqs. (77)–(79) has been discussed by Reis (2003), who show that, depending on the equation of state, the results may agree or disagree with the general relativistic approach. We leave this problem open since it is not our purpose here to investigate the validity of these equations in detail. Furthermore, our main results concern the evolution of the cosmic fluid governed by the Friedmann Eqs. (80)–(82). Since the Friedmann equations can be derived from the theory of general relativity, the validity of these equations is not questioned. Equations different from Eqs. (77)–(79) have been considered by Pace et al. (2010) and Harko (2011a). Their equations involve an additional term ṗu/c² in the Euler Eq. (78). However, this term apparently leads to equations for the cosmic evolution that are different from the Friedmann equations, so this term is not considered here.

8.2. The cosmic evolution of a BEC universe

We assume that dark matter is a BEC with quartic self-interaction¹⁰ described by the barotropic equation of state (8). For convenience, we write this equation in the form $\mathit{p}\mathrm{=}\mathit{k}{\mathit{c}}^{\mathrm{2}}{\mathit{\rho }}^{\mathrm{2}}\mathit{,}$(83)where k = 2πa_sħ²/m³c² is a constant that can be positive (a_s > 0) or negative (a_s < 0). In a first step, we concentrate on the dark matter component and neglect radiation, baryonic matter, and dark energy (in particular we take Λ = 0). These additional terms are considered later. For the equation of state (83), the Friedmann Eq. (80) becomes $\frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{b}}}{\mathrm{d}\mathit{t}}\mathrm{+}\mathrm{3}\frac{\mathit{ȧ}}{\mathit{a}}{\mathit{\rho }}_{\mathrm{b}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{k}{\mathit{\rho }}_{\mathrm{b}}\mathrm{)}\mathrm{=}\mathrm{0.}$(84)This equation can be integrated into ${\mathit{\rho }}_{\mathrm{b}}\mathrm{=}\frac{\mathit{A}}{{\mathit{a}}^{\mathrm{3}}\mathrm{-}\mathit{kA}}\mathit{,}$(85)where A is a constant. If we make the physical requirement that Eq. (85) has solutions for large a, we must impose A > 0. In that case, ρ_b ~ A/a³ for a →  + ∞ which returns Eq. (18). When k > 0, which is the case previously considered by Harko (2011a), the density exists only for a > a_∗ = (kA)^1/3. For a → a_∗, ρ_b →  + ∞. When k < 0, the density is defined for all a and has a finite value ρ_b = 1/|k| when a → 0.

Combining Eqs. (81) and (83), we obtain $\frac{\begin{array}{}\mathrm{¨}\\ \mathit{a}\end{array}}{\mathit{a}}\mathrm{=}\mathrm{-}\frac{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{b}}}{\mathrm{3}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{3}\mathit{k}{\mathit{\rho }}_{\mathrm{b}}\mathrm{)}\mathit{.}$(86)When k > 0, the universe is always decelerating ($\begin{array}{}{}_{\mathrm{¨}}\\ \mathit{a}\end{array}\mathit{<}\mathrm{0}$). When k < 0, the universe is accelerating ($\begin{array}{}{}_{\mathrm{¨}}\\ \mathit{a}\end{array}\mathit{>}\mathrm{0}$) for ρ_b > ρ_c ≡ 1/3|k| and decelerating ($\begin{array}{}{}_{\mathrm{¨}}\\ \mathit{a}\end{array}\mathit{<}\mathrm{0}$) for ρ_b < ρ_c. Using Eq. (85), this can be expressed in terms of the radius¹¹: the universe is accelerating for a < a_c ≡ (2A|k|)^1/3 and decelerating for a > a_c.

To determine the temporal evolution of a(t), we assume κ = 0 (flat space) as in the Einstein-de Sitter universe. Combining Eqs. (82) and (85), we get $\mathit{ȧ}\mathrm{=}{\left(\frac{\mathrm{8}\mathit{\pi GA}}{\mathrm{3}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{\mathit{a}}{\sqrt{{\mathit{a}}^{\mathrm{3}}\mathrm{-}\mathit{kA}}}\mathrm{·}$(87)We note that, for a →  + ∞, this equation reduces to Eq. (22) with Eq. (18) so we recover the Einstein-de Sitter solution (23). It is convenient to define a_∗ = (|k|A)^1/3 and introduce R = a/a_∗. In that case, the density is given by ${\mathit{\rho }}_{\mathrm{b}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\left|}\mathit{k}\mathrm{\right|}}\frac{\mathrm{1}}{{\mathit{R}}^{\mathrm{3}}\mathrm{\mp }\mathrm{1}}\mathit{,}$(88)and Eq. (87) can be rewritten as $\mathit{Ṙ}\mathrm{=}\frac{\mathit{KR}}{\sqrt{{\mathit{R}}^{\mathrm{3}}\mathrm{\mp }\mathrm{1}}}\mathit{,}$(89)where the upper sign  −  corresponds to k > 0 and the lower sign  +  corresponds to k < 0. We have defined the constant $\mathit{K}\mathrm{=}{\left(\frac{\mathrm{8}\mathit{\pi G}}{\mathrm{3}\mathrm{\left|}\mathit{k}\mathrm{\right|}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}{\left(\frac{\mathrm{4}\mathit{G}{\mathit{m}}^{\mathrm{3}}{\mathit{c}}^{\mathrm{2}}}{\mathrm{3}\mathrm{\left|}{\mathit{a}}_{\mathrm{s}}\mathrm{\right|}{ħ}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(90)which depends on the mass m of the bosons and on their scattering length a_s. Introducing the dimensionless parameter λ/8π = a_s/λ_c = a_smc/ħ (Chavanis 2011a) measuring the strength of the short-range interactions (λ_c is the Compton wavelength of the bosons), we can rewrite Eq. (90) in the form $\mathit{K}\mathrm{=}\frac{\mathrm{2}}{\sqrt{\mathrm{3}}}{\left(\frac{\mathit{m}}{{\mathit{M}}_{\mathrm{P}}}\right)}^{\mathrm{2}}\sqrt{\frac{\mathrm{8}\mathit{\pi }}{\mathrm{\left|}\mathit{\lambda }\mathrm{\right|}}}{\mathit{t}}_{\mathrm{P}}^{-1}\mathit{,}$(91)where M_P = (ħc/G)^1/2 is the Planck mass, l_P = (ħG/c³)^1/2 the Planck length, and t_P = l_P/c the Planck time.