Using kinetic theory to examine a self-gravitating system composed of baryons and cold dark matter

We examine the evolution of non-relativistic cold dark matter gravitationally coupled to baryons with modes deep inside the Hubble radius (sub-horizon regime) using a kinetic theory approach within the realm of Newtonian theory. We obtain the general solution for the total density perturbation and we also show that a baryon perturbation catches up with the dark matter perturbation at late times, which in turn makes possible the formation of bound structures. We extend the linear perturbation analysis by considering the turn-around event, collapse of matter, and its virialization process.


Introduction
Current lore associated with the standard cosmological model would indicate that on scales larger than 100 Mpc the universe is essentially homogeneous and isotropic, but on smaller scales there are some deviations from the mean density in the form of galaxies, galaxy clusters, amongst other configurations [1,2]. A natural question is to ask: How do such structures grow in the universe? What is the basic mechanism to aggregate matter and make it collapse in the form of a bump? The first attempt to give an answer to that question was made by Jeans long time ago [3,4]. He focused on the necessary condition under which small perturbations of a gas cloud could grow exponentially, leading to the collapse of the cloud and therefore ending in the formation of stars. In other words, the Jeans mechanism describes the gravitational instability of a self-gravitating gas cloud. Interestingly enough, there exist plenty of ways to understand the Jeans mechanism and to derive such a criterion. To be more prea e-mail: kremer@fisica.ufpr.br b e-mail: martin@df.uba.ar cise, consider an initially stable and static cloud that can be initially perturbed by the environment such as a shock wave, passing spiral arms of the galaxy, etc. This configuration can collapse if the inwards directed gravitational force is bigger than the outwards directed pressure force. The critical maximal radius that allows stability depends essentially on the Newton constant, the speed of sound and and the matter density; in fact, it leads to the idea that the denser the clouds, the more unstable they become. During the collapse only a part of the gas ends up in stars, so many stars form out of one collapsing cloud, which means that young stars are born in clusters [2].
At the cosmological level, one way to get some insight in the structure formations is to look at a simplified treatment, which in this case corresponds to the Newtonian approach [2]. Indeed, the theory of Newtonian structure formation is sufficient to understand most of the processes which are well within the horizon. To do so, one must derive the full Newtonian hydrodynamics in an expanding universe and it turns out that the Boltzmann kinetic theory is the natural way to achieve such a goal [5][6][7].
The first step to understanding how the cooperative effects of baryons and dark matter may work in the process of structure formation is by inspecting a collisionless selfgravitating system composed of two components and then solving the coupled collisionless Boltzmann and Poisson equation together [8,9]. A system composed of baryons and dark matter leads to a total Jeans mass which is smaller than the one associated with a single component, indicating that a smaller amount of mass is needed to ignite the collapsing process. One could expect that the bumps with masses greater than the Jeans mass initiate the collapsing process, but an over-dense region in an expanding universe eventually recollapses and virialises. In the case of a single component with an expanding background it turned out that the "swindle" proposal may be avoided, while the Jeans instability is expected to arise in the limit of large wavelengths [8]. Nevertheless, one must stress that the Jeans instability is not only restricted to the Newtonian (or General Relativity) realm and it can emerge within the context of alternative gravity theories as well [10][11][12][13][14].
This paper is organized as follows. In Sect. 2, one presents the kinetic theory formalism for dealing with a selfgravitating system of two components within the framework of Newtonian cosmology and by doing so one also examines the general conditions to achieve the Jeans instability. Besides, the perturbation of baryons and dark matter are studied along with the conditions under which the total matter starts to virialize. In Sect. 3, the conclusions are stated. We will use the metric convention (+, −, −, −) and nongeometric units in which 8π G = 1 and c = 1 unless stated otherwise.

Kinetic theory and self-gravitating components
The spacetime evolution of the one-particle distribution function f (r, v, t) in the phase space spanned by the space and velocity coordinates (r, v) is ruled by the Boltzmann equation. The one-particle distribution function f (r, v, t) is defined by assuming that f (r, v, t)d 3 rd 3 v gives the number of particles in the volume element d 3 r about the position r and with velocity in the range d 3 v about v at time t. The Boltzmann equation in the absence of collisions between the particles but in the presence of a gravitational potential reads (see e.g. [8][9][10]) (1) On large scales the universe is spatially homogeneous and isotropic and then one can portray such a geometry using the Friedmann-Lemaître-Robertson-Walker metric [1,2]. In particular, for a spatially flat universe the line element is where a(t) is the so-called cosmic scale factor. If the universe is dominated by a perfect fluid then Einstein's field equations are reduced to some coupled differential equations known as the Friedmann and acceleration equations, respectively: Here the dot denotes the derivative with respect to time, G is the gravitational constant, while ρ and p are the mass density and the pressure of the source that generates the gravitational field.
In the case of a universe filled with dark matter and baryons one can treat such components as pressureless ones ( p 0); then the mass density as a function of the cosmic scale factor is ρ = ρ 0 (a 0 /a) 3 . Replacing this density in (3) leads to a power-law scale factor in terms of the cosmic time, where ρ 0 and a 0 describe the values of the mass density and cosmic scale factor at t = 0, respectively. In order to describe a system composed of dark and baryon matter subject to a gravitational field one can make use of a Boltzmann equation for each constituent plus a Poisson equation where both components are coupled gravitationally. From now on, one will consider two distributions functions, one for baryons,  (1): Here the index j = {b, d} denoted baryons and dark matter, respectively. The gravitational field must fulfill the Poisson equation: Provided these components are in equilibrium their corresponding distribution functions must be the well-known Maxwellians written in a comoving frame: where σ j is the dispersion velocities associated with the baryons and dark matter, respectively. Indeed, the dispersion velocities σ j do not depend on the space coordinates and are proportional to 1/a(t) for an epoch of post-recombination, namely One uses the possibility to relate the observationally measured coordinates r, also known as the physical or proper coordinates, to the comoving coordinate x by r(t) = a(t)x. One then has the proper velocity of a given observerṙ(t) = a(t)x + aẋ; thus the second term is the so-called peculiar velocity of objects. In the absence of peculiar velocities one has the Hubble-Lamaître's law v H = H (t)r, H (t) =ȧ/a being Hubble's parameter. Thus the expansion of the universe leads to objects moving away from us at a speed proportional to the distance. Further, the Hubble parameter defines a time scale and can be used to define a length scale over which physical processes act coherently. The comoving Hubble radius is defined as d H = c/H . This is also the scale at which general relativistic effects become important. For l d H Newtonian gravity is often adequate [2].
Let us justify properly why the assumption of Newtonian gravity is good enough for our purpose of studying the evolution of dark matter and baryons with the help of kinetic theory. In general the process of structure formation should be studied in a fully covariant manner within the theory of General Relativity [1]. However, Newtonian gravity is generally valid at most scales of cosmological interest. It is valid in regions which are small compared to the Hubble radius l d H and large compared to the Schwarzschild radius of most massive black holes. Newton equations of motion are recovered in the weak field limit of General Relativity, i.e., when c 2 . Conversely if one assumes that particle motions are non-relativistic, v/c 1, one can reduce an arbitrary metric to the following form in the Newtonian limit [2]: ds 2 = (c 2 + 2 )dt 2 − dr 2 , with c 2 . Here is the effective Newtonian potential that should be used in the equations of motion. Consider the following coordinate transformation law for a spatially flat FLRW metric: In the above transformation the spatial coordinates have been changed from comoving to the proper coordinates and the time has been corrected for the gravitational redshift. The transformation changes the metric to Here one has ignored cubic and higher order terms in r/d H 1 because one is only interested in regions that are small compared to the Hubble radius. The coefficient in front of dr 2 can be set to unity if v/c 1 and r/d H 1. The line element (12) then reduces to ds 2 = 1 −ä a r 2 c 2 c 2 dt 2 − dr 2 − r 2 d 2 . After comparing the latter expression with the metric in the Newtonian limit one can identify the effective Newtonian potential due to the homogeneous and isotropic background universe: Equation (13) must coincide with the gravitational potential obtained from the Poisson equation, Eq. (7). The previous finding shows that the potential and the Newtonian approximation are valid only when the particle motions are nonrelativistic, and the scales of interest are much smaller than the Hubble radius. Gravitational clustering occurs at scales which are much smaller than the scale of homogeneity, which current studies estimate to be 60 h −1 Mpc. Such a scale is indeed much smaller than the Hubble scale, 3000 h −1 Mpc. One can therefore study gravitational clustering in the Newtonian framework [2]. Summarizing the above argument, one corroborated that, deep inside the Hubble radius, Newtonian perturbation theory is in agreement with relativistic perturbation theory because the typical velocities are small compared to the speed of light. Then Newtonian and relativistic perturbation theory have to agree on the relation between single-and multi-species evolution on sub-Hubble scales [15,16]. Such a statement would confirm the results obtained from the Nbody numerical simulations within the Newtonian gravity approach being almost the same (on a sub-horizon scale) as those which emerge from solving the full non-linear Einstein equation numerically [17].
One would like to examine the basic theory for structure growth in the expanding universe within the framework of kinetic theory. To do so, one introduces small perturbations in the distributions functions along with some perturbation in the gravitational potential, where the perturbations are assumed to be small compared with the zeroth order quantities. Replacing representations (14) and (15) into the Boltzmann and Poisson equations (5) and (7) leads to the following system of equations for h b , h d and 1 : As is customary one neglected the products of and ∂ v d h d provided these terms are of second order in the perturbed variables. Furthermore, the perturbed variables will expand in a plane wave base where the physical wavenumber vector is q/a(t), while the comoving one is named q. Thus, the factor 1/a(t) in the wavenumber takes into account that the wavelength is stretched out in an expanding universe: where the amplitudes φ and h 1 j depend on time. In fact, the aforesaid amplitudes h 1 j can be written as a linear combination of the collision invariants 1, v j − H r, v j − H r 2 of the Boltzmann equations: Here A j , B j , D j are time dependent amplitudes. Now replacing (18)- (20) in (16) and (17) yields At this point one comment is in order. So far one has dealt with the Boltzmann equation for baryons and dark matter plus the Poisson equation for a fixed background cosmology. Nevertheless, one did not use any balance equations. Then a rather natural question to ask is: How will one recover a second order master equation for the total contrast density and for the partial contrast density of baryons and dark matter? The answer to that question is intrinsically connected to the kinetic theory in itself and how to obtain some derived physical quantities from the Boltzmann equation. In other words, one needs to calculate the perturbed balance equations and the way to obtain such first order balance equations is by multiplying Eq. (21) with the collision invariants 1, v j − H r, v j − H r 2 , and then integrate using some well-known results of Gaussian integrals. It turns out that one arrives at the following system of equations: From (23) and (25) it follows that d A j /dt = 0, so one can choose A j = 1 for simplicity. If one introduces the density contrasts for each component as δ ρ j = mh j d 3 v/ρ j = A j + 3σ 2 j D j , and then combines (23) and (24) along with the Poisson equation, Eq. (22), one gets a coupled system of differential equations for the density contrasts which reads where it was useful to recast Eq. (26) in a more canonical way by introducing a shift in the contrast density as follows: To identify the Jeans wavenumber for the total system (baryon plus dark matter) it requires one to introduce the total mass density ρ t = j ρ j along with the total density contrast ρ t δ t = j ρ j δ ρ j . As a result of multiplying (26) by ρ j /ρ t and then summing over the two species we find that the time evolution of the total density contrast can be written as As is well known the total adiabatic sound speed can be recast as v s = 5σ 2 /3, while for a mixture of two components it Hence for the third term on the left-hand side of (27) In order to see that the relations ρ t δ t = j ρ j δ ρ j and j 5ρ j σ 2 Here the source term corresponds to the total density and does not involve the perturbed density contrast. Equation (28) tells us that a dissipative effect enters through the usual friction term proportional to 2H . One must emphasize that the physical Jeans scale (length or wavelength) is obtained by demanding that the term δ t vanishes, namely, the Jeans wavenumber is q J = ( √ 3/2)Ha/v s . The Jeans length can heuristically be derived by balancing the sound crossing time, t s ∝ a/v s q J , with the gravitational free-fall time, t ff ∝ 1/ √ Gρ t , which yields the same result as mentioned above in an intuitive manner. Moreover, this scale seems to be sensitive to the thermal dispersion velocity of the particles through v s and the total material content ρ t ; notice that one considers the universe after decoupling, then one can focus on dark matter and baryons only. One should mention also that the Jeans wavenumber for each component in principle is different provided the propagation speeds are not the same. Here one is looking at the Jeans length of the composed system. In this context, the comoving Jeans wavenumber can be written in terms of the scale factor as q J = (4π Gρ 0 3 while the comoving Jeans length is simply λ J = 2π/q J . One can prove that perturbations can grow for q ≤ q J ; otherwise they just oscillate ( q ≥ q J ). The latter result is consistent with the behavior of the Jeans scale in a universe dominated by matter after decoupling time, q J ∝ a 1/2 [1]. One final comment: the fact that σ j = σ 0 j a −1 (which is equivalent to having a pressure p j = σ 2 j ρ j ∝ a −5 ) implies that the collisionless fluid may be treated as pressureless as long as q < q J . As one is trying to arrive at the general solution without demanding the latter condition the system of equations will be much harder to solve than the usual case. One will examine the situation where the condition q > q J holds as well.
In order to carry on it is mandatory to find a way to solve the above equation. One route is to make explicit the Hubble factor: A more useful expression than the one given by (29) can be obtained by performing a change of variables and introducing a new time, η = (2/3) ln(H 0 t), and a constant Jeans wavenumber, q 0 J = (4π Gρ 0 t ) 1/2 a 0 /v 0 s ; in this way one is able to get rid of the H 2 factor, yielding Let us analyze what kinds of limiting regimes follow from (30). For large values of the wavelength with respect to the Jeans wavelength, q/q 0 J = λ 0 J /λ 1, Eq. (30) can be approximated by whose solution as a function of the usual cosmic time is given by where C 1 and C 2 are integration constants. This solution represents the well-known growth of the total density contrast for a matter-dominated universe proportional to t and the solution of Eq. (33) in terms of the usual cosmic time admits a closed form, Here A 1 and A 2 are integration constants, while Ci(x) = − ∞ x dt (cos t/t) and Si(x) = x 0 dt (sin t/t) represent the cosine and sine integrals, respectively. We infer from Eq. (34) that it represents oscillations of the total density contrast with respect to time t for values of 4 J not too small. When the time increases, the factor 4 becomes much smaller than unity and the total density contrast will grow provided that the cosine integral term will increase for large values of the cosmic time.
In order to gain more insight in the behavior of the total density contrast with respect to cosmic time we introduce another dimensionless time, called τ = (H 0 t), and write (29) as follows: where the prime refers to differentiation with respect to τ . We solved (35) for the initial conditions δ t (10 −3 ) = 0.1 and δ t (10 −3 ) = 0 (say) in the limit cases of small and large wavelengths. The aforesaid numerical simulations are displayed in Figs. 1 and 2. Figure 1 tells us that for small wavelengths with respect to the Jeans wavelength, λ 0 J /λ = 50, the total density contrast exhibits some small oscillations for small values of time τ but eventually as τ increases the total density contrast begins to grow. In the opposite limit, with λ 0 J /λ = 0.5 we find that there is only exponential growth of the total density contrast as can be seen from Fig. 2.
The next task to tackle is to solve the master equation, Eq. (26), for each species, at least under some reasonable conditions. Equation (26) can be written as where τ = H 0 t, δ ρ j ≡ δ j and we used the 3 . For some given ratios ρ 0 j /ρ 0 t , v 0 s j /v 0 s , and q/q J 0 one can solve the coupled system of differential equations (37) numerically. As we have done in our previous analysis we will select the same values for the ratio q/q 0 J = λ 0 J /λ in the case of small wavelengths in relation to the Jeans scale (λ 0 J /λ = 50) and in the case of large wavelength λ 0 J /λ = 0.5. Accordingly, we begin by recalling that the mass density ratio ρ 0 d /ρ 0 b could be associated with the density parameter ratio 0 d / 0 b , which is about 5.5 [18], i.e., ρ 0 j /ρ 0 t = 0 d / 0 b ≈ 5.5. Perhaps not surprisingly, the ratio of the sound speeds v 0 s j /v 0 s is not usually fixed in the literature. However, we will follow [8,9] and use the value v 0 sd /v 0 sb = 1.83, which was inferred from the sim-  Figure 3 shows that, for small wavelengths with respect to the Jeans wavelength, λ 0 J /λ = 50, both contrast densities exhibit some oscillations for small values of τ , but as time proceeds they start to grow. For large wavelengths with respect to the Jeans wavelength, λ 0 J /λ = 0.5, we find that both density contrasts obey exponential growth, implying the development of the Jeans instability [cf. Fig. 4]. Indeed, we note that due to the imposed initial conditions the baryons density contrast starts from zero but at later times it becomes equal to the dark matter density contrast, since the slope of the baryon density contrast is more accentuated than the dark matter ones.
In characterizing the properties of the composed system made of dark matter and baryons and how they evolve with time, we are certainly interested in double-checking the previous numerical analysis by solving analytically the set of equations in terms of the cosmic time variable η = (2/3) ln(H 0 t). As such, it is important to notice that the abundance of the two components is encoded in the master equation (29), which reads We need to corroborate that Eq. (37) leads to the correct physics; otherwise the current formalism has no merit at all. In particular, one has to ensure that baryons can catch up with dark matter perturbations after recombination [2]. Of course, baryons only contribute a small amount to the total matter density, but all the visible structures in the universe are made of them. One explores the simplest case in which the modes satisfy the condition q q J for both baryons and dark matter. Although both components act as sources to the gravitational potential perturbation, one will assume that dark matter dominates the background density so the perturbed Poisson equation will have only as a source the perturbed dark matter density. In fact, one also will neglect in Eq. (37) the term proportional to δ j provided it appears multiplied by the factor (q/q J ) 2 , which is considerably small in this regime. The aforesaid approximations yield the following system of equations for dark matter and baryons: where f d0 = ρ 0 d /ρ 0 t . The general solution for the evolution of the dark matter perturbations is given by the superposition of three terms, namely, a constant mode, a decay mode and a growing mode: where the exponents are λ ± = − 1 4 ± 1 16 + 3 2 f d0 while δ d0 and δ ± are some constants. Given the fact that the growing mode is the relevant one at late times it will be enough to consider that δ d δ + e λ + η . Note that the growing mode of the dark matter perturbations provides the driving term in the master equation for baryons. Armed with this new fact, one can return to the task of finding the behavior of baryons which it involves to provide some initial conditions. One will assume that at early times, nearly the recombination era, Then the general solution for baryons reduces to Equation ( 41) tells that at very early times δ b 2λ + δ d0 e −η/2 , while in the opposite limit, at late times, δ b δ d0 e λ + η = δ d . So the baryon perturbation catches up with the dark matter perturbation at late times. This is another reason why we need dark matter. Without dark matter to set up gravitational potential wells for the baryons to fall into, it would be hard to explain how we can have bound baryon structures today, given that δ b (η i ) is of the order of 10 −5 . Thus, baryon perturbations would still be in the linear regime today without the help of dark matter perturbations, making bound structures like galaxies and clusters of galaxies impossible. Nevertheless, for a better understanding of the clustering properties of a system composed of dark matter plus baryons we will need to extend the analysis to a scenario where the Boltzmann equation is solved in curved spacetime for the composed system and the fully relativistic corrections of General Relativity are taken into account as well [20]. The next step of consistency is to explore the non-linear regime of total density perturbation. Suppose now, as a toy model for the formation of non-linear (gravitationally bound) structures, one has a spherical over-density of radius R and mass M embedded into the otherwise homogeneous spatially flat and matter-dominated universe (EdS). Given the fact that it is over-dense, this configuration will reach a maximum radius and subsequently contract until collapse. Such a toy model is a reasonable approximation provided the distribution of the dark matter in the universe can be considered as composed of individual so-called halos, approximately spherical over-dense clouds of dark matter which can reach highly non-linear densities in their centers. This means that one can work with a Newtonian equation of motion for the radius, where the mass of the halo can be expressed in terms of the turn-around radius and the turn-around density as M = 4π 3 R 3 ta ρ ta . The density at the turn-around event involves the critical density of the background and the over-density ζ of the halo with respect to the background at turn-around, thus ρ ta = ζρ c (a ta ) = (3H 2 ta /8π G)ζ . It is useful to introduce some more convenient variables by parameterizing all the physical quantities in terms of variables evaluated at the turnaround event, that is to say, R = R ta y, a = xa ta and τ = H ta t with H ta = H 0 a 3/2 ta . Now, the Newton equation for the spherical halo reads In order to solve (48) one could consider the case where initially the halo has zero radius and then reaches its maximum radius at the turn-around event. These assumptions imply that the boundary conditions for the shell must be: y (x = 1) = 0 and y(x = 0) = 0. After integrating once and imposing the boundary condition y (x = 1) = 0 one arrives at the velocity of the shell in terms of its radius, y = ± √ ζ (y −1 − 1) 1/2 ; the minus stands for the evolution of the shell before the turn-around event and the plus sign is applied after the turnaround. Using the Friedmann equation one can get a parametric solution τ = τ [y]: The shell must be able to collapse and by doing so it must form a structure, so the expansion of the shell must come to a halt and the radial velocity must vanish at a finite time. This stage of zero radial velocity is the so-called turn-around event where the perturbation reaches its maximum radius. For the case under study such a situation occurs at x = 1 = y and τ ta = 2/3, which in turn implies that the over-density of the halo is ζ = (3π/4) 2 . Furthermore, the collapse time is reached at τ coll = 2τ ta = 4/3 for x coll = 3 √ 4. To obtain the parameters that characterize the collapsing phase, one begins with the behavior of τ at early times. Making an expansion of τ [y, ζ ] at the lowest order possible in y around y = 0 yields τ (8/9π)y 2/3 (1 + 3y/10). Furthermore, one defines the over-density inside the halo in relation to the background density as ≡ ρ halo ζ /ρ bg = ζ(x/y) 3 . Replacing the expansion for τ into the definition of overdensity inside the halo leads to = 1 + 3y/5. Then the linear density contrast inside the halo is δ = − 1 = 3y/5. By extrapolating this formula linearly until the turn-around event one gets δ ta = δ(y)/x 3y/5x but x[τ (y)] ζ −1/3 y so the contrast linear density at the turn-around is δ ta (3/5)ζ −1/3 1.06. In the subsequent phase after the turnaround, when the collapse is reached, the inner contrast linear density is δ coll x coll (3/5)ζ −1/3 1.69. In other words, the halo of dark matter has already collapsed when its expected linear density reaches the value δ coll 1.69. Notice that the previous result is independent of the mass M, the initial over-density, and the epoch of virialization. The next step in the evolution is to consider what happens when the halo reached the virial equilibrium: essentially the potential energy of the halo must be twice that at the turn-around event and its radius must decrease at y v = 1/2, implying that v = (2x coll /) 3 ζ 178. Therefore, a halo in virial equilibrium is expected to have a mean density nearly 178 times higher than the background although in numerical simulations the density contrast is fixed at the value 200, furnishing in this way a natural definition of the virial radius of a virialized object [2]. The lesson from this simple analysis is that density perturbations can form bound structures generated by gravitational collapse after they become 200 times as dense as the background. Such a result seems to be consistent with the full results from N-body simulations where galaxies and clusters of galaxies separate out as distinct gravitationally bound structures when their densities are at least 100 times greater than the background density [2]. However, this must be thought of just as a heuristic rule provided one is using the fact that the linear theory is valid until coll ≥ 1, where the process of virialization cannot be stopped.
Having mentioned the ideal picture of halos, one must also say that the halo is not necessarily isolated from the background and there will take place a constant inflow of material into the halo, or the halo might even merge with another halo. Thus, the evolution of a halo after its formation is quite nontrivial so it cannot be easily analyzed within this simplistic model. However, if one assumes that the continuous inflow of material and the merging with other halos produce new halos which are still characterized by the same density ratio v = (2x coll /) 3 ζ 178, then one would expect that the mean density of a typical halo with a given mass M scales with redshift like ρ v = v ρ bg ∝ (1 + z) 3 , whereas its physical radius should go as r v ∝ ρ −1/3 v . Notice that the overdensity does not depend on the mass of the perturbation, on the initial over-density, nor on the epoch of virialization t v . Thus, whenever one observes an over-density of the order of δ v , one positively assumes that the corresponding structure is virialized (or close to virialization) irrespective of its mass or formation history.
Last but not least, it is physically meaningful to devote some efforts to understand how the virialization process takes places in more detail. The main goal will be to estimate the redshift for which the virialization occurs within the context of an EdS expanding universe. As mentioned before, the virialised density (inside the halo) must be at least 100 bigger than the background density (criterion): If the particles inside the virialized region have a velocity dispersion σ d , the virial theorem states that V = −2 K , V being the gravitational potential, K the kinetic energy, and the brackets indicate the average value. From the latter fact, one gets M v σ 2 d = G M 2 v /R v and therefore the radius is Inserting the latter result into Eq. (45) one obtains Taking as fiducial value 0 bg h 2 = 0.11 one gets a more or less heuristic result, Equation (47) tells us that for a typical galaxy as the Milky Way with σ d = 300 kms −1 and a total virialized mass of M v = 10 12 M , the virialized redshift should be around z v ≤ 7, so clusters of galaxies must have formed not too long ago in terms of the total evolution of the universe. Given the fact that our universe is currently accelerating it seems interesting to explore how our analysis is modified once a dark energy component (a cosmological constant) is taken into account. To do so, we proceed as before by assuming that we have a sphere (halo) of radius R enclosing a mass m and its dynamic is governed by Newton equation: where is a cosmological constant term. Notice that the cosmological constant appears provided the background cosmology has changed and the Friedmann equation now reads Equation (48) has as a first integral the energy conservation equation (per unit mass), The total energy contained in the shell of radius R can be obtained by integrating Eq. (49). In doing so, we consider that the differential mass is dm = ρ4πr 2 dr ; then the total mass becomes M = 4π R 3 ρ/3 when the total density is constant.
We would like to determine how the cosmological constant affects the virialization radius. In fact, we could expect that the inclusion of the cosmological constant would potentially affect the turn-around radius and therefore the virialization radius as well. To prove this, it is convenient to compute the total potential energy contained in a sphere of constant density using Eq. (49). It reads The virial theorem states that K = j=g,cc n j 2 V j for V j ∝ R n j , where n g = −1 for the purely gravitational potential and n cc = 2 corresponds to the potential energy associated with the cosmological constant. The interesting point as regards the theorem is that it allows us to replaced the average kinetic energy by the average potential energy in the energy conservation law. We need to compare the turn-around event with zero kinetic energy with the state of the sphere at the virialized radius. At the turn-around event the total energy is E ta = − 3G M 2 5R ta − 10 M R 2 ta , while at the virialized radius we obtained E v = 2 V cc + V g /2. We define the two dimensionless variables x = R v /R ta and ϑ = R 3 ta /3G M. The energy conservation equation (E ta = E v ) allows us to obtain a cubic equation for x = R v /R ta in order to estimate the effects introduced by : In the case without cosmological constant we basically have x = 1/2; however, we expect that the ratio R v /R ta should be less than 1/2 provided the cosmological constant contributes the total energy in the form of V cc = − 10 M R 2 . To confirm our reasoning, it is convenient to write the ratio R v /R ta as a small deviation of 1/2, namely x = 1/2 + δ with δ 1. Replacing x = 1/2 + δ in (51) and expanding at first order in δ we are able to write δ in terms of ϑ: (52) Substituting (52) in the definition of x we find that the ratio R v /R ta is given by which gives x < 1/2 for 0 < ϑ < 1 as a general result. In fact, we can go further by considering that ϑ 1 and ending with x 1/2(1 − ϑ/4) < 1/2. We can conclude that the virialized radius is smaller if the dark energy component of the universe is taken into account. In a way, we can say that the cosmological constant term in the total potential energy helps the system to reach an equilibrium configuration much faster. N -body numerical simulations show that R v 0.483R ta for a non-vanishing [21].
In order to guarantee that the system has virialized we must check that the average density inside the halo must be much bigger than the critical density of the background, with 0 = 1 − 0 bg and ζ 100. We would like to obtain a linear relation between σ 2 d and R v ; however, the inclusion of the cosmological constant spoils that possibility. Indeed, the virial theorem leads to a relation between σ 2 d and R v that is non-linear: Instead of determining z v we can consider that z v is fixed and derive the typical order of magnitude for the virial radius with the help of (54), Notice that the maximum virial radius allowed with nonvanishing is smaller than the maximum virial radius allowed without cosmological constant (57). The maximum virial radius reads We can estimate the virial radius as R v μR max v where the μ parameter is selected in such a way that the inequality (57) is satisfied; namely, we choose μ 10 −2 . Besides, we can verify that the second term in (55) is considerably small if R v μR max v . We know the observational value for the cosmological constant, = 1.105 × 10 −6 Mpc −2 , and we take as the virialized mass M v = 10 12 h −1 M , so the ratio where h = 0.7 and G being the Newton constant. The latter fact confirmed that idea that the cosmological constant introduces a small deviation in the dispersion velocity obtained from the virial theorem.

Summary
In this work, we analyzed the dynamics and the collapse of a collisionless self-gravitating system composed of dark matter and baryonic matter. This system is described by two Boltzmann equations, one for each component, gravitationally coupled through the Poisson equation. First, we derived the general solution for the total density perturbation and then we solved numerically the coupled master equation for both components in different cases, showing that for modes deep inside the Hubble horizon and under the condition q q J 0 the density of matter tends to grow. In fact, we showed that a baryon perturbation catches up with the dark matter perturbation at late times, so dark matter is the driving force which provides the gravitational potential wells for the baryons to fall into, allowing them to create bound structures. Besides, we also examined in broad terms what happens with this toy model in the non-linear regime by working within the realm of Newtonian theory. We considered the formation of non-linear structures in the form of a spherical over-density of radius R and mass M embedded into the homogeneous spatially flat and matter-dominated universe. We explored the linear regime until the turn-around event followed by the collapse and the subsequently virialization process. The viral theorem along with the fact that the virialization criterion implies that galaxies were formed at low redshift, say less than ten. We also discussed the case with non-zero cosmological constant.
It is important to mention the main elements that we should introduce in the kinetic theory to study a composed system with three components (dark matter, baryons, and dark energy) within the Newtonian approach without going to the full relativistic formulation. We have basically three distribution functions with their corresponding Boltzmann equations along with one Poisson equation that couples gravitationally all the aforesaid components. We may consider the case of dark energy decoupled from the baryon plus dark matter system. The possibility of an extra-coupling between dark matter and baryons has been suggested by the Experiment to Detect the Global Epoch of Reionization Signature (EDGES) measuring the 21-cm absorption signal from primordial neutral hydrogen at redshift z 17 [22][23][24]. This signal is considerably stronger than what is expected from the vanilla cosmic model and contains potentially new information about the true nature of dark matter. If that is the case under study, the interaction between dark matter and baryons must be specified through the collision operator. Such an interaction will affect the structure formation due to the exchange of momentum, energy and flux of heat; however, the specific result will strongly depend on the collision operator and therefore the specific interactions considered (say decay, annihilation, etc.). In fact, an interaction could potentially suppress the growth of structure in the early universe. Of course, we could expect that the structure formation will be totally ineffective at late times due to the overall expansion of the universe. Besides, one possible extension of the current work is to consider also the kinetic theory with a collision operator that takes into account the interaction between dark matter and dark energy in order to explain the current excess of brightness in the 21 cm line as proposed by several authors in the literature [25][26][27][28]. We will address the latter possibility within our formalism in the near future.