Gravitational instability of the primordial plasma: anisotropic evolution of structure seeds

We study how the presence of a background magnetic field, of intensity compatible with current observation constraints, affects the linear evolution of cosmological density perturbations at scales below the Hubble radius. The magnetic field provides an additional pressure that can prevent the growth of a given perturbation; however, the magnetic pressure is confined only to the plane orthogonal the field. As a result, the"Jeans length"of the system not only depends on the wavelength of the fluctuation but also on its direction, and the perturbative evolution is anisotropic. We derive this result analytically and back it up with direct numerical integration of the relevant ideal magnetohydrodynamics equations during the matter-dominated era. Before recombination, the kinetic pressure dominates and the perturbations evolve in the standard way, whereas after that time magnetic pressure dominates and we observe the anisotropic evolution. We quantify this effect by estimating the eccentricity $\epsilon$ of a Gaussian perturbation in the coordinate space that was spherically symmetric at recombination. For perturbations at the sub-galactic scale, we find that $\epsilon = 0.7$ at $z=10$ taking the background magnetic field of order $10^{-9}$ gauss.


General Remarks
Our theoretical knowledge of the Universe is based on the Standard Cosmological Model, that provides a convenient framework to satisfactorily explain the majority of cosmological observations, like the anisotropy pattern of the cosmic microwave background (CMB) [1,2,3], the large-scale structure of the Universe [4,5,6], the Hubble diagram of distant type Ia supernovae [7,8,9] and the abundances of light elements [10].
The Standard Cosmological Model relies on the assumption that the Universe, at least at large scales, is highly homogeneous and isotropic, and its geometry is thus described by the Robertson-Walker metric. In fact, the distribution of luminous red galaxies shows that the present Universe is homogeneous on scales greater than ∼ 100 Mpc [11], while the isotropy of the CMB itself (which has a black-body distribution at T = 2.73 K with temperature fluctuations of order 10 −5 or less) is an indication of the isotropy of the Universe as a whole and a strong evidence for homogeneity at the time of hydrogen recombination (nearly 400.000 years after the Big Bang, corresponding to a cosmological redshift z rec = 1100). On the other hand, below the "homogeneity scale" of 100 Mpc, the distribution of matter is definitely inhomogeneous. Such a dichotomy between the smoothness in the matter-energy distribution at z = z rec and the clumpiness of the recent Universe (for z 1) below a certain scale is explained by the mechanism of gravitational instability: the structures we observe today have been formed through Email address: lattanzi@fe.infn.it (Massimiliano Lattanzi) the growth of tiny density perturbation seeds that, accordingly to the currently accepted model, were created in the early Universe during a phase of inflationary expansion.
The presence of a large scale (i.e., coherent over a Hubble length), strong magnetic field is forbidden by the observed isotropy, as it would naturally single out a preferred spatial direction. However, a background and uniform magnetic field could be present at cosmological scales provided that its intensity is small enough. In particular, upper limits on the present field intensity of order ∼ 10 −9 G have been derived from observations of the CMB temperature anisotropies [12,13] and of its temperature-polarization correlation [14,1] 1 . Smaller scale fields, on the other hand, could be as strong as 10 −6 G.
The effects of large-scale magnetic fields on the evolution of cosmological structures have been studied extensively in the literature (for a complete review, see Ref. [16] and references therein), where both Newtonian and general relativistic treatments, the latter often using covariant and gauge-invariant techniques, are present. It is known that, among others, magnetic fields slow down (and possibly prevent) the growth of perturbations and can produce vorticities and shape distortions in the density field [17,18]. In this Letter, our goal is to revisit the issue of the existence of a "magnetic Jeans length" and of its dependence on the direction along which the perturbation propagates, other than to give a realistic estimate of its value in a matter-dominated Universe. The presence of a magnetic Jeans length has been discussed in several papers [16,17,18,19,20] both in a Newtonian and in general relativistic framework (for an analysis of the standard Jeans mechanism in the presence of dissipative effect see also [21,22,23]), but some of the analyses failed to recognize its angular dependence. Here, we present a neat derivation of the relevant instability scales bases on the equations of magnetohydrodynamics (MHD) on an expanding Universe. We discuss a simple generalization of the magnetic Jeans length suited for two-fluid systems, and clarify some misunderstandings that are present in the literature. We also show the results of the numerical integration of the coupled MHD and Poisson equations. Finally, we numerically study the distortion introduced in the density field by the anisotropy in the critical length.
The Letter is organized as follows. In Sec. 2, we characterize the Universe as a plasma. In Sec. 3, we introduce the basic equations, and in Sec. 4 we carry out the linearization procedure. We derive the existence of the magnetic Jeans length from analytical considerations in Sec. 5, while in Sec. 6 we show some numerical results. Finally, we draw our conclusions in Sec. 7.

The plasma features of the pre-and post-recombination Universe
In this Section, we aim at characterizing the plasma features of the cosmological fluid. Between the time of e + e − annihilation (i.e., for a temperature 2 T m e , corresponding to a redshift z ∼ 10 9 ) and the present (z = 0), the matter-energy content of the Universe is provided by electrons, protons and neutrons (the three species being collectively referred to as baryons in the cosmological jargon), photons, neutrinos, and two elusive components dubbed dark matter and dark energy, that presently account for more than 99% of the total energy budget of the Universe. However, dark energy has been subdominant for most of the past history of the Universe and can be safely neglected at redshifts z > 1. Dark matter and neutrinos interact only gravitationally with the other components and can be neglected as long as the plasma properties of the fluid are concerned.
The cosmological baryon-to-photon ratio is extremely small and equal to n b /n γ 6.1 × 10 −10 , where n γ and n b are the photon and baryon number densities, respectively. Both n b and n γ scale with redshift as (1 + z) 3 ; their present values are n b (z = 0) 2.5 × 10 −7 cm −3 and n γ (z = 0) 410 cm −3 . Most of the baryons in the Universe are in the form of 1 H nuclei (i.e., isolated protons) so that, for simplicity, in the following we assume n b = n p (n p being the proton number density). Furthermore, n p is also equal to the electron number density n e = n p , because of the charge neutrality of the Universe.
We study separately the properties of the cosmological fluid before and after the time of hydrogen recombination occurring at z = z rec = 1100 (T 0.25 eV). For z > z rec , the protons and electrons are free and thus one deals with a fully ionized 2 All throughout the paper, we use natural units with c = = k B = 1. plasma. In this regime, photons and baryons are tightly coupled due to Thomson scattering, and share a common temperature where the present photon temperature T 0 γ = 2.73 K. After recombination (z < z rec ), most of the electrons and protons exist in the form of neutral hydrogen atoms, and only a small residual ionized fraction x e = 2.5×10 −4 survives, making the fluid a weakly ionized plasma. In this regime, the photon temperature still scales as (1 + z), while that of baryons evolves in the same way only until z = 100, due to residual scatterings that keep them in thermal equilibrium with photons; after that time, their temperature decreases faster, as (1 + z) 2 . The baryonic fluid remains neutral until the time of reionization, when the UV radiation produced by the first stars ionizes again the hydrogen present in the cosmological medium. This is likely to have happened around z 10, however the precise details of the reionization history are still largely unknown, and for this reason we limit our analysis to redshifts z > 10.
In the following, we will assume the presence of a background homogeneous magnetic field B(z), whose contribution to the total energy density of the Universe can be considered negligible. We recall that the field intensity B(z) scales as (1 + z) 2 and, unless otherwise stated, we take its value at the present time to be B(z = 0) = 10 −9 G.

The pre-recombination Universe
A fundamental quantity characterizing a plasma is the Debye length λ D , namely the length over which electrons screen out electric fields in a plasma. It defines the length scale over which a system can consistently considered to be a plasma. The Debye length of a hydrogen plasma at temperature T is λ D = T 4πn e e 2 (6.9 cm) where e is the proton charge. Using T = T γ and the values given above, one gets The redshift dependence of λ D implies that the comoving Debye lengthλ D ≡ λ D (1 + z) is constant during the cosmological evolution and equal toλ D 2 × 10 4 cm. Plasma effects can be important in a system when its physical dimension L is much larger than the Debye length. For the Universe, the relevant length is the Hubble radius L = l H ≡ H −1 , where H is the Hubble parameter. This length represents the maximum scale at which microphysical processes can operate in order to establish the thermodynamical equilibrium. Today, l H 10 28 cm; during the matter-dominated era l H ∝ (1 + z) −3/2 , while in the radiation-dominated era l H ∝ (1 + z) −2 . From the analysis of both these scales, it is evident that l H λ D turns out in the period considered here. Moreover, under the same hypotheses, the baryonic matter M D within a Debye sphere is also constant and given by where m p is the proton mass. This "Debye mass" clearly results to be much smaller than any other of cosmological interest and we can conclude that the cosmological fluid can be considered as neutral at all relevant scales. Another meaningful index is the so-called plasma parameter N D , i.e., the number of particles within a Debye sphere: The dependence of λ D and n p on the redshift z implies that also N D is a constant. In particular, since N D 10 7 1, the cosmological fluid results to be a weakly coupled plasma.
In order to provide a complete characterization of the cosmological plasma, we now turn our attention to the plasma dissipative properties, starting from the plasma resistivity η. For an electron-proton plasma, this is given by η = m e ν ei / n e e 2 , where m e is the electron mass and ν ei is the electron-ion collision frequency. For the case under consideration, ν ei is well approximated by the electron-electron collision frequency ν ee [24], i.e., where ln Λ C is the Coulomb logarithm, introduced to quantify the effects that small-angle-diffusion collisions have in the Coulomb scattering. A simply estimate of Λ C in a plasma is given by Λ C 12πN D , so that for the cosmological fluid, the Coulomb logarithm is 20. Substituting Eq.(5) into the expression for the resistivity given above, we get Close to recombination, the cosmological plasma has an electric resistivity equal to η(z rec ) 1.6 Ω cm, i.e., a conductivity 0.6 siemens cm −1 , a value typical of a semiconductor [in Gaussian units, η(z rec ) = 1.8 × 10 −12 s].
Let us now turn our attention to the viscous properties of the plasma. The shear viscosity coefficient of matter strongly coupled with radiation can be expressed as [25,26] where a S B 5.7×10 −8 W K −4 m −2 is the Stefan-Boltzmann constant, while τ denotes the mean collision time between particles and can be estimated as τ (n γ σ T v) −1 (here, v c and we have introduced σ T 6.6 × 10 −29 m 2 as the cross section for Thomson scattering). For a photon gas at equilibrium at temperature T , n γ 2.0 × 10 7 (T/K) 3 m −3 and then The resistivity and viscosity coefficients enter the MHD equations through the following diffusion coefficientsη ≡ η/4π andη v ≡ η v /ρ, where ρ is the density of the fluid. Taking The relative magnitude of the viscous and magnetic diffusion rates can be parameterized through the magnetic Prandtl number Pr m ≡η v /η. Using the expression above forη v andη, we obtain Pr m 5.0 × 10 13 so that Pr m 1, i.e., viscous diffusion is more important than resistive diffusion, at recombination and indeed always at the redshifts under consideration.
Having discussed the relative importance of viscosity-and resistivity-driven dissipative effects, we now analyze at which scales these effects are relevant. In a magnetized plasma, a useful parameter is the Lundquist number S ≡ Lv A /η, where L is a typical length scale and v A = (B/4πρ) 1/2 is the Alfvén velocity. The Lundquist number is basically the ratio between the resistive diffusion timescale τ r = L 2 /η and the Alfvén crossing We also compare τ A with the viscous diffusion timescale defined as τ v = L 2 /η v . The ratio S v ≡ τ v /τ A can be thought as a viscous analogous to the Lundquist number: The baryonic mass (at the average background density) contained in a sphere of radius equal to the length where S ∼ 1 can be calculated to be Both mass scales are well below the values of cosmological relevance at the redshifts of interest. The discussion above shows how the following hierarchy among the relevant time scales holds for all mass range and redshifts of interest: meaning that viscosity always dominates over resistivity, and that both dissipative effects can indeed be neglected when studying the propagation of Alfvén waves.

The post-recombination Universe
After the time of recombination z rec = 1100, the cosmological plasma exists in a weakly ionized state, the neutral and ionized components having densities ρ n ρ b and ρ i = x e ρ b ρ n , respectively. We take the residual ionization fraction x e constant and equal to 2.5 × 10 −4 in the range z rec > z > 10. The results of the previous Subsection can be generalized to show that the hierarchy (14) holds also in this regime. In spite of the small value of the ionization fraction, the magnetic field could still affect the dynamics of the whole system in view of the interactions between neutral and charged particles. In particular, the magnetic forces acting on the charged particles can be communicated to the neutrals through collisions. However, if the coupling is not tight enough, the neutrals feel the magnetic field but drift with respect to the ions in a process termed ambipolar diffusion. Its relevance at a given length scale L is quantified by the ambipolar Reynolds number R amb [27,28,29,30] where γ in = 1.9 × 10 −9 cm 3 s −1 [31] is the ion-neutral drag coefficient due to collisions between the two species, v 2 A = B 2 /4πρ n is the Alfvén velocity in the tightly coupled limit, v is the characteristic velocity of the fluid, and L amb is the ambipolar length, i.e., the scale where R amb = 1. The ambipolar Reynolds number is just the ratio between the ambipolar diffusion timescale If R amb 1, the neutrals are uncoupled from the plasma. On the contrary, when R amb 1 the dynamics of the two components can be described through ordinary single-fluid MHD with an additional dissipative term [29]. This term becomes progressively less important as R amb grows and can be neglected in the limit R amb 1, or L L amb . Assuming that the evolution of the fluid is driven by Alfvénic phenomena, i.e., v ∼ v A , the ambipolar length is given by and R amb = τ amb /τ A . Thus, R amb 1 if and only if the tightcoupling condition τ ni τ = τ A is satisfied. It is straightforward to check that for 1100 > z > 10, it is always R amb > 1 at scales larger than a few tens of comoving kiloparsecs, meaning that the following hierarchy holds: so that the ions and neutrals are tightly coupled and ambipolar diffusion can be safely neglected. In order to better illustrate this point, in Figure 1 we plot the mass contained in a sphere of radius L amb at the background baryon density as a function of redshift. It is evident that, in the redshift range considered, R amb 1 for all scales M 10 6 M . In this regime, we can therefore neglect the dissipative term mentioned above and use single-fluid ideal MHD. We conclude that ambipolar diffusion does not affect the dynamics of the cosmological plasma after recombination.

Basic equations
In this Section, we derive the basic equations describing the linear evolution of instabilities in the cosmological fluid, mod- eled as a magnetized plasma. In this respect, we underline that the full investigation of the perturbative dynamics of the Universe would require a general-relativistic treatment, in order to correlate the matter and geometrical fluctuations 3 . However, as long as one is interested in scales much smaller than the Hubble radius, i.e., L H −1 , a Newtonian treatment provides a consistent description of the dynamics. Nonetheless, in this scenario the expansion of the Universe can be accounted as the bulk background motion of the fluid [25,33,34].
The starting point of our treatment is the Eulerian set of equations governing the fluid motion, on which one can develop a perturbative theory by adding small fluctuations to the unperturbed cosmological background solution. The zerothorder dynamics is derived by considering a flat homogeneous and isotropic Universe whose energy density is dominated by non-relativistic matter, and correctly describes the expansion of the Universe. We assume that a background magnetic field is present, whose contribution to the total energy density of the Universe can be considered negligible.
Let us now start by briefly recalling the basic equations of non-relativistic, ideal and single fluid MHD, which govern the plasma motion. The mass conservation and the Newtonian gravitational field are described by the continuity and Poisson equations, the single-fluid dynamics is described by the Euler equation in presence of a magnetic field B and, finally, the electromagnetic interaction can be summarized by the frozen-in and the Gauss laws. Such equations read respectively, where ρ is the mass density, v is the velocity field, Φ is the gravitational potential and G is Newton constant. This system constitutes the base of our perturbative approach.
To derive the zeroth-order dynamics, we assume the usual Robertson-Walker metric, i.e., ds 2 = dt 2 − a 2 (t) d 2 , where a = a(t) represents the cosmological scale factor, and a perfect fluid energy-momentum tensor as the matter source of the gravitational field, i.e., T µ ν = diag [ ρ 0 , −P 0 , −P 0 , −P 0 ], with ρ 0 = ρ 0 (t). In this scheme, the behavior of the mass density with time is obtained from the energy-momentum conservation law T ν 0; ν = 0 and from the Friedmann equation, i.e., respectively (the dot (˙) denotes the total derivative with respect to synchronous time). Here H =ȧ/a is the Hubble parameter and K = const. is the curvature factor. Setting the matter-dominated Universe equation of state (EoS) P 0 ∼ 0 (P 0 ρ 0 ) in Eq. (19), the zeroth-order solution of the system (18) turns out to be whereρ andB 0 are dimensional constants, r (r =|r |) denotes the radial coordinate vector and, of course, a(t) satisfies Eq. (20). We observe how this non-stationary solution characterizing the background dynamics is not affected by the so-called "Jeans swindle" proper of the static solution [25,33,34].
To obtain now the explicit time dependence of the unperturbed quantities involved in the model, we restrict the analysis to the flat case, i.e., K = 0. From the Friedmann equation (20) and using the solution for ρ 0 , one readily obtains a = 6πGρ Finally, we recall that the adiabatic sound speed is defined by v s = ∂P/∂ρ. For a general specific heat ratio γ, we assume that the pressure varies as P = Kρ γ , so that the speed of sound is given by (22c)

Perturbation scheme
In order to analyze the implications that the physics of an ideal magnetized plasma can have on the structure formation, we will follow the standard perturbation approach. In this respect, we consider small perturbations around the zeroth-order cosmological solution derived above, i.e., we write ρ = ρ 0 + ρ 1 (with ρ 1 ρ 0 ) and similarly for the other quantities P, v, Φ and B. Substituting the perturbed quantities in Eqs. (18) and keeping only terms up to first order, one gets where, as already discussed, the pressure and density perturbations have been related through the adiabatic sound speed, i.e., P 1 = v 2 s ρ 1 . We are assuming that B 2 0 /4πρ 0 = v 2 A 1, where B 0 = |B 0 |, in order to preserve the isotropy of the background flow.
In the following, we replace B 1 with the dimensionless magnetic fluctuation b 1 ≡ B 1 /B 0 . Moreover, the analysis of the system above can be simplified by Fourier-transforming the spatial dependence of the involved quantities, i.e., using perturbations in the form of plane waves, taking with φ 1 = {ρ 1 , v 1 , Φ 1 , B 1 } and k is the physical wavenumber scaling as 1/a(t). It is convenient to consider also the comoving wavenumber q = ak, that stays constant during the expansion. The evolution for a given harmonic can be obtained by the equations in real space with the substitutions φ 1 →φ 1 , ∇ → ik and ∂ t → ∂ t − iH(k · r). In the following, for the sake of simplicity, we will drop the tilde over the Fourier transformed variables. Then, the system (23) reduces to (hats denote unit vectors): where we have already eliminated Φ 1 by means of the Poisson equation in k−space, i.e., k 2 Φ 1 = −4πGρ 1 . It is understood that the constraint k · b 1 = 0 always hold. Decomposing now v 1 in its components v 1 and v ⊥ 1 parallel and orthogonal to the direction of q respectively, i.e., v 1 = v 1q + v ⊥ 1 (where v ⊥ 1 ·q = 0), and introducing the following scalar variables: we finally get a further simplified system: where we have defined and, of course, 0 µ 1. We stress that ω 2 0 is not positive definite.

Evolution of the density contrast and conditions for collapse
The form (27) of the evolution equations has the advantage that it clearly expresses the relationship between the physical quantities involved, other that being very well suited for numerical integration. Some further analytical insight can however be gained by reducing it to an unique higher-order equation for the variable δ(t).
Considering the case of a matter-dominated Universe, and using the explicit time dependence of the quantities involved in the model in that case, i.e., Eqs. (22), with some algebra one can derive the following fourth-order differential equation for δ(t): where δ ( ) denotes the th derivative of δ with respect to time, and we have defined the following constants: We recall that γ is the specific-heat ratio (P ∼ ρ γ ) and that γ 4/3, i.e., ν 0. The most general solution of Eq.(29) for δ is found to be the superposition of four independent solutions δ i (i = 1, . . . , 4), given by: where p F q [(a 1 , ..., a p ); ( b 1 , ..., b q ); z] denotes the generalized hypergeometric function of argument z, the A i 's are arbitrary integration constants and The constant coefficients a and b depend, in general, on ν, Λ 2 , µ and we report their complete expressions in Appendix A.
We are now interested in discussing the asymptotic behavior of the hypergeometric functions in the limit of very small or very large argument, i.e., Λ 1 /4ν 2 t 2ν 1 or 1. As in the nonmagnetic case discussed in Ref. [25], we restrict the analysis to the range 0 ν 1/3, i.e., treating the standard regime 4/3 γ 5/3. From the asymptotic expansion of the F functions in the case of large argument, i.e., Λ 1 /4ν 2 t 2ν 1, the density contrast always shows a damped oscillating behavior with time. In fact, in this regime it always exists at least one asymptotic solution proportional to positive power of the argument Λ 1 /4ν 2 t 2ν , which results to be the leading term of the solution superposition. In this case, δ decreases with time.
On the other hand, in the limit Λ 1 /4ν 2 t 2ν → 0, the asymptotic expansion of the solutions (31), can be written as In order for the gravitational collapse to occur, at least one of the modes has to be growing, i.e., x i > 0. It is fairly easy to show that x 1 , x 2 and x 4 are always negative, whereas the sign of x 3 depends on µ and Λ 2 . In particular, when µ 0, we obtain x 3 > 0 irregardless of the value of Λ 2 , while when µ = 0, x 3 is positive if Λ 2 < 2/3. This means that, on the plane orthogonal to the magnetic field (µ = 0), a new stability condition arise if the magnetic field is strong enough. The threshold value related to Λ 1 that should discriminates the two regimes of growing and decreasing density contrast can be set as Λ 1 = 1 [25]. Remembering that ρ 0 = 1/6πGt 2 , such condition rewrites in terms of the wave number as which is substantially the same as the usual Jeans condition for gravitational instability. In fact, in the non-magnetic case, (to which our analysis reduces for ω A = 0) this is the only criterion that separates the growing and the decaying modes [25].
In a similar way, the new threshold Λ 2 = 2/3 yields to the following condition Summarizing, we find that the presence of a background magnetic field introduces an anisotropy in the stability criterion. While outside the plane orthogonal to B 0 , the stability of the perturbations is dictated only by the standard Jeans condition k ≷ k J , on that plane the unstable modes are those for which the conditions k < k J and k < k A both hold 4 . In other words, if k A < k J (basically equivalent to v A > v s ), there are Jeans-unstable modes (those in the window k A < k < k J ) that, in the orthogonal plane, are stabilized by the magnetic pressure. The window of stable modes gets wider for larger values of the ambient magnetic field, as expected. We underline that these results are qualitatively the same as those obtained for a static and uniform background [35]. A similar analysis was carried on by the authors of Ref. [17] obtaining a similar results. However, their derivation contained a mistake when separating the 4 It is easy to verify how the physical meaning of the condition is that the timescale for gravitational collapse τ c ∼ L/v c ∼ L 3 /GM ∼ 1/Gρ 0 is much shorter than both the acoustic and Alfvén timescales τ s ∼ L/v s ∼ 1/kv s and τ A ∼ L/v A ∼ 1/kv A . real and imaginary components of the evolution equations [25]. For this reason they find a second-order differential equation instead than the fourth-order one discussed here.

Numerical Analysis
In the previous Section, we have gained an important insight on the effect of a background magnetic field on the evolution of density perturbations. We now show some results obtained through the direct numerical integration of the differential system (27).

Preliminaries
We will focus on the period of the cosmological evolution that goes from the onset of matter domination (z 3000) to the time of reionization (z 10). We start from matter domination because, before that time, the growth of density perturbation was slowed down and practically frozen by the rapid expansion of the Universe. In the matter-dominated Universe, a ∝ t 2/3 and H = 2/3t. We can ignore the presence of a dark energy component since this is sub-dominant until very recent times. The time period that we consider can be divided into two distinct phases, i.e., before and after the recombination of hydrogen occurring at z rec = 1100. Before recombination, the baryons are completely ionized and they are tightly coupled to photons, at least at scales larger than the comoving photon mean free path λ γ 1.8[(1 + z)/(1 + z rec )] −2 . At these scales, the total pressure of the fluid is given by radiation pressure. After recombination, most of the protons and electrons are in the form of neutral hydrogen atoms, leaving a small ionized fraction x e 2.5 × 10 −4 . At the scales of interest, the neutral and ionized components are tightly coupled by collisions (see Sec. 2.2) and can be treated as a single fluid. However, photons are now free streaming so that the baryon pressure is given just by kinetic pressure, dropping down by several orders of magnitude with respect to its pre-recombination value.
In view of this, we take the speed of sound of the cosmological fluid before and after recombination to be [25] v 2 s | z>z rec = 1 3 respectively, where σ = 4a S B T 3 /3n b k B 1.5 × 10 9 is the specific entropy. We recall that T b = T γ = T γ | z=0 (1 + z) for z > 100, while afterwards T b ∝ (1 + z) 2 . The expression above for the sound speed is rigorously valid only for scales corresponding to baryonic masses 10 11 M , that stay above the photon mean free path until the time of recombination. At smaller scales, baryons lose the radiation support before recombination, roughly when the given scale goes below the photon mean free path, and it is at this time that the switch between the two expressions in Eq.(36) should take place. In the following, we shall consider masses nearly as small as 10 8 M , for which the photon decoupling effectively takes place at z 3000 (very close to the time of matter-radiation equality). However, we choose to switch between the two expressions for the sound speed at z = z rec irregardless of scale and shall comment later on how we expect this choice to affect our results. Following the discussion in Sec. 2.2, the Alfvén velocity is taken to be where ρ b should always be intended as the total baryon density, both before and after recombination. A plot of the Alfvén and sound speeds as functions of redshift is shown in the upper panel of Figure 2. The sudden drop in the sound speed at recombination is due to the sharp decrease of baryon pressure after photon decoupling.
In a detailed model, the presence of different uncoupled components making up the matter content of the Universe should be taken into account. In fact, most of the matter (∼ 80%) is in the form of cold dark matter (CDM), interacting with the baryon-photon fluid only through the gravitational force. Thus, a proper treatment should rely on a two-fluid description. In the following, we shall ignore perturbations in the CDM component, however we argue that we can still draw meaningful conclusions about the perturbations in the baryonic component.
In fact, in the pre-recombination era, the large radiation pressure prevents baryons to fall into the potential wells created by CDM; this is known to be true in the non-magnetic case but we expect it to hold also in the case under consideration since, as seen in the last Section, the magnetic field only acts to in- crease stability. Thus in this regime the CDM and baryon perturbations are effectively decoupled. After recombination, the baryon density perturbations will take some time to catch up with those in the CDM component and we expect our treatment will rigorously remain valid for some time.
Before showing the results of the numerical integration, we illustrate in the lower panel Figure 2 the evolution of the standard Jeans wavenumber (34) and of its "magnetic" counterpart (35). In order to take out the change in k J and k A due to the expansion, we follow the convention to express the results in terms of the mass contained inside the corresponding length scales 1/k J,A . In particular, we consider the total baryonic mass (irrespective of the ionization state), contained inside a sphere of radius 2π/k J,A . It can be seen that the window of modes that are made stable by the magnetic field, i.e., those between the red dashed and the black solid line, spans, right after recombination, five orders of magnitude in mass.
As noted in the introduction, the existence of a magnetic Jeans length has been studied previously and all expressions for the critical wavenumber agree, apart from numerical factors, with expression (35). However, the numerical estimates of this and associated quantities that are found in the literature sometimes differ from our results. The reason seems to be that often the density ρ 0 that appears in Eq. (35) is taken to be the present critical density ρ c 9 × 10 −27 kg m −3 . This yields at the present time a magnetic Jeans length λ A ∼ 1/k A ∼ 10 kpc for B 0 (z = 0) = 10 −9 G and λ A ∼ 1 Mpc for 10 −7 G 5 . Using instead the baryon density for ρ 0 will yield values of λ A a factor ρ c /ρ b = Ω −1 b 20 larger, i.e., λ A ∼ 0.2 (20) Mpc for B 0 (z = 0) = 10 −9 (10 −7 ) G. This amounts to a factor 20 3 10 4 difference in the corresponding mass scale 6 . 5 The latter value has also been said to be of the order of the scale of a galaxy cluster, while in effect it is closer to the scale of a galaxy, as it can be seen by the fact that the mass enclosed inside a sphere of 1 Mpc radius at the critical density is ∼ 10 11 M . The reason why a galaxy is much smaller than 1 Mpc is that it has detached from the Hubble flow and undergone non-linear evolution, so that its density is much larger than the cosmological average [33]. 6 The discussion so far has ignored the gravitational action of dark matter;

Results
We now discuss the results of the direct numerical integration of Eqs. (27). The initial conditions for the integration have been chosen using the fact that power-law solutions for δ can be found in the limit t → 0. There are four distinct solutions of this kind, but only one corresponds to a growing mode. We have matched the initial conditions to the asymptotic growing solution at the initial time of integration. The latter has been chosen so that all the modes of interest were outside the horizon at that time. Even if the initial time falls in what would be the radiation-dominated era, nevertheless we always consider a matter-dominated Universe. All results have been normalized to the initial value of the density contrast.
In Figure 3, we show the evolution of the density contrast for three different wavenumbers k = (17, 1.7, 0.36) Mpc −1 (normalized at the present time), i.e., for the following baryonic masses M = (1.7 × 10 8 , 1.7 × 10 11 , 1.7 × 10 13 ) M . These masses roughly correspond to the scale of a dwarf galaxy, of a galaxy and of a galaxy cluster respectively. For each mode, we show the evolution in both the direction parallel to the background magnetic field (µ = 1) and orthogonally to that direction (µ = 0). In all cases, the perturbation is initially growing but then starts to oscillate once the Jeans mass (that is growing) becomes larger than the mass of the perturbation. This happens earlier for smaller scales. In this phase the magnetic pressure does not play any role, as the much larger radiation pressure is actually providing the force that prevents the collapse. In fact, there is no difference in the evolution parallel and orthogonal to the field, as the radiation pressure is isotropic. The situation changes dramatically after recombination, when the baryon pressure drops and only the magnetic pressure can possibly oppose the growth, at least in the plane orthogonal to the field. Thus, perturbations in the direction of the field can grow unhindered, while the perturbations that are orthogonal can be stathis can be roughly taken into account by using the total matter density ρ m at the numerator of Eq. (35), while keeping the same expression for v A . This makes the value of λ A roughly twice smaller than in the case in which only baryons are considered, i.e., λ A ∼ 0.1 (10) Mpc for B 0 (z = 0) = 10 −9 (10 −7 ) G. bilized. As it can be seen from Figure 3, this is what happens for perturbations at the dwarf galaxy scale: at z = 10, the relative growth of parallel perturbations with respect to orthogonal ones is of order 100. For perturbations at the galactic scale and larger, instead, the evolution is basically the same in all directions. This can be understood by looking at the lower panel of Figure (2), from which it is clear that the pressure induced by a magnetic field of 10 −9 G can only stabilize perturbations with mass 10 10 M . We recall that we have neglected the fact that, at the dwarf galaxy scale, the support of radiation pressure is not lost at recombination but some time before (see previous Subsection). From the discussion above, it is clear that, had we taken into account this fact, the evolution of orthogonal and parallel perturbations would have begun to differentiate earlier. This goes in the direction of enhancing the anisotropic growth of perturbation and the "squeezing" effect studied in the following.
The fact that, after recombination, the evolution of the density contrast in the presence of a magnetic field changes for different directions leads to the reasonable expectation that some degree of anisotropy will be generated even in initially symmetric structures. In order to show this, we consider a Gaussian density fluctuation with standard deviation σ in coordinate space at recombination: where the x are comoving coordinates centered at the maximum of the perturbation. After Fourier-transforming, we separately evolve the different harmonics in momentum space using Eqs. (27) and we finally transform back to obtain the perturbation in coordinate space at a later time. In Figure 4 we show, for a background magnetic field directed along the y axis and with a present intensity of 1 nG, the evolution of a perturbation with σ = 0.05 kpc at recombination (so that the 3σ region encloses a mass M 1.5 × 10 8 M at the mean baryonic density, i.e., roughly the mass of a dwarf galaxy). In particular, we show equal density contours at z = 1000, 100 and 10. It is evident from the figure how the perturbation becomes progressively squeezed along the direction orthogonal to the magnetic field.
In order to quantify the anisotropy in the perturbation, we consider the isodensity contour corresponding to half the value at the peak and calculate its eccentricity = 1 − b 2 /a 2 , where a and b are the lengths of the semi-major and semi-minor axis of the contour, respectively. In Figure 5, we show how the eccentricity changes with redshift; for the parameters used above, we get 0.7 at z = 10.

Conclusions
We have studied the effect of a background magnetic field on the linear evolution of cosmological density perturbations at scales well below the Hubble length, where a Newtonian treatment can be used, focusing on the matter-dominated era. The conditions that allow for the growth of small density perturbations have been clearly stated. In particular, we have found that in the plane orthogonal to the ambient magnetic field, a new critical length appears, related to the presence of the magnetic pressure, while everywhere else outside that plane the stability is dictated by the standard Jeans criterion. This is also confirmed through a direct numerical integration of the relevant MHD equations during the matter-dominated era, and this effect is shown to be possibly important after recombination, when the magnetic pressure of baryons is much larger than their the kinetic pressure. Finally, it has been shown how the dependence of the critical scale on the angle between the perturbation wavevector and the magnetic field could lead to a sizable anisotropy in the perturbations at sub-galactic scales at the onset of non-linearity.
Our analysis has relied on some approximations: in particular, we have ignored the gravitational effects of dark matter perturbations. We have argued that this approximation limit the validity of our treatment to some time after recombination. We defer a more detailed and fully general relativistic analysis, also taking into account the different fluid components, to a future work.