Linear perturbations of an anisotropic Bianchi I model with a uniform magnetic field

In this work, we study the effect of a magnetic field on the growth of cosmological perturbations. We develop a mathematical consistent treatment in which a perfect fluid and a uniform magnetic field evolve together in a Bianchi I universe. We then study the energy density perturbations on this background with particular emphasis on the effect of the background magnetic field. We develop a full relativistic solution which refines previous analysis in the relativistic limit, recovers the known ones in the Newtonian treatment with adiabatic sound speed, and it adds anisotropic effects to the relativistic ones for perturbations with wavelength within the Hubble horizon. This represents a refined approach on the perturbation theory of an isotropic universe in GR, since most of the present studies deal with fully isotropic systems.


Introduction
The formation of large scale structures across the Universe is one of the most fascinating and puzzling questions, still opened in theoretical cosmology. Among the long standing problems of this investigation area is the determination of the basic nature and dynamics of the cold dark matter [1], responsible for the gravitational skeleton on which the baryonic matter falls in, forming the radiative component of the present structures.
However, also the peculiarity of the matter distribution across the Universe, in particular the possibility for large scale filaments [2], as well as hypotheses for structure fractal dimension [3,4] call attention for a deeper comprehension.
In this respect, we observe that the Universe plasma nature, both before the Hydrogen recombination and, for a part a e-mail: federico.digioia@uniroma1.it b e-mail: giovanni.montani@enea.it in 10 5 also in the later matter dominated era [5,6], has to be taken into account.
At the recombination the Universe Debye length is of the order of 10 cm and therefore the implementation of a fluid theory, like General Relativistic Magneto-hydrodynamics is to be regarded as a valid and viable approach to treat the influence of the primordial magnetic field [7] on the evolution of perturbations [6]. Nonetheless, the smallness of such magnetic field, as constrained by the Cosmic Microwave Background Radiation (CMBR) up to 10 −9 G [8,9,10,11,12,13,14], significantly limits the impact of the plasma nature of the cosmological fluid on the evolution of perturbations. As shown in [6,15], the presence of the magnetic field is able to trigger anisotropy in the linear perturbations growth and it can be inferred that in the full non-linear regime, such anisotropy grows up to account for the formation of large scale filaments.
Apparently, a weak point in the perspective traced above consists of the small plasma component surviving when the Hydrogen recombines and in the observation that the most relevant cosmological scales enter the non-linear regime in such a neutral Universe. Instead, it can be surprisingly demonstrated [5,6,15] that the coupling between the neutral and ionized matter is very strong at spatial scale of cosmological interest (for overdensities of mass greater than 10 6 solar masses, the Ambipolar Reynold number is much greater than unity for redshift 10 < z < 1000). Thus, the dynamical features, for instance anisotropy, that we recover for the plasma component clearly concern the Universe baryonic component too. This statement is not affected by the presence of dark matter gravitational skeleton in formation, simply because the radiation pressure prevents, up to z ∼ 100 the real fall down of the baryonic fluid into the gravitational well. In fact, the large photon to baryon ratio, about 10 9 (also constant during the Universe evolution), maintains active a strong Thomson scattering process, even after the hydrogen is recombined into atoms [16,17,18,5,6].
These considerations are to underline that a single fluid General Relativistic Magneto-hydrodynamics formulation is an appropriate tool to investigate the impact of the Universe plasma features on structure formation, at least for a large range of the cosmological thermal history.
In this context many works have been developed, mainly assuming as negligible the backreaction of the magnetic field on the isotropic Universe, see [19] and references therein. However, the presence of a magnetic field rigorously violates the isotropy of the space and the (essentially) flat Robertson-Walker geometry must be replaced by a Bianchi I model. This paper faces the general question of how the linear perturbations evolve on a background Bianchi I cosmology, thought as a weak perturbation of the isotropic case, but treated in its full generality for arbitrary large magnetic fields.
We discuss in detail the structure of the perturbation equations in the synchronous gauge and the specific form of the spectrum time dependence in specific important limits, like the large scale limit, when the dependence on the wavenumber can be suppressed, and the sub horizon limit, when the dependence on the wavenumber is dominant.
Furthermore, the change of the Jeans scale, when passing from the ionized to the (essentially) recombined Universe, is determined for the small scales, shedding light on the role of the magnetic field and on the real nature of the gauge perturbations.
We recover the slowing-down of the growing mode in super-horizon scales, long known in FRW models. This effect is very small given the upper limits on the cosmological magnetic fields, of order O v 2 A 1. At sub-horizon scales, we generalise the solutions of [20] and [21], which in turn generalise the results of [19]. While they consider random (i.e. isotropic) magnetic fields to preserve the FRW model, we work in the anisotropic case and also consider a nonvanishing sound speed.
Finally, we stress that, along the whole analysis, we compare our results with previous achievements in literature, providing a significant contribution to the understanding of the different effects that the Universe anisotropy, due to the magnetic field, induces on the perturbation evolution and stability.
We notice that there is another paper about this matter [22], which was the first analytical study to address this issue. There, the authors study the model in 3 different physical limits with specific anisotropies, while we completely relate the background anisotropy to the magnetic field.
The paper is structured as follows: in section 3 we summarize the exact GRMHD equations in the 3+1 covariant formalism; in section 4 we find the solution for the background Bianchi I model, then we write the equations for the perturbations in synchronous gauge in section 5 and we find the gauge modes in section 6; finally we solve our system in some specific cases in section 7 and we compare our results with present literature.
2 General properties of the Bianchi I models As we already said, it is impossible to accommodate a magnetic field in a isotropic model. Moreover, although present observations show that the isotropic FRW model describes very well the present universe, it is only a very special description of the universe towards the initial singularity, while the general one should incorporate anisotropy [23,24].
In the first stage of the universe evolution the matter contribution is negligible, while it is necessary to have a isotropic matter field to achieve the isotropization of the model [25,26]. The general solution is constructed through the Bianchi VIII and IX models [23,24,26], but we will focus for simplicity on a single Kasner era and so we will use a Bianchi I model. The Bianchi I model is similar to the FRW one, but with three different scale factors. It is intrinsically anisotropic in vacuum, i.e. the three cosmic scale factors are never all equal; moreover, in vacuum one of the three scale factor always decreases with time, meaning that one of the spatial direction is contracting.
Near enough to the cosmological singularity, any matter source in the form of perfect fluid energy density, having equation of state p = wρ always behaves as a test fluid, i.e. it induces negligible backreaction, as far a 0 < w < 1. Since the background magnetic field energy density is a radiationlike term in the Universe and it is associated to an equation of state p = ρ/3, near enough to the singularity, we can expect a typical vacuum solution of the Kasner form [24,26].
The more general Bianchi IX model can be described as a succession of Kasner epochs, in which the different directions exchange time evolutions, alternating moments of growing and decreasing [26]. For more detailed informations regarding the Bianchi models we recommend [27].
Clearly, as soon as the Universe expands enough, the matter source can no longer be negligible and, if the pressure term is isotropic, the solution must correspondingly isotropize, i.e. the three scale factors tend to be equivalent. This process of isotropization is particularly efficient in the case of an inflationary paradigm [28,26], when a vacuum energy, having an equation of state p = −ρ is dominating the Universe dynamics.
The relevance of our study for the structure formation takes place when the isotropization process reduced the Bianchi I cosmology to a flat Robertson-Walker Universe, except for the residual intrinsic anisotropy due to the presence of a background magnetic field.
There exist already a large number of studies regarding Bianchi I models, analysing cases with different values for the barotropic index w of the matter source in addition to the magnetic field. [29] was probably the first to address their stability. [25] studies the effect of a pure magnetic matter component, [30] contains analytic solutions for dust w = 0 and radiation w = 1/3, [31] contains solutions for w = 1 and 1/3 ≤ w ≤ 1 and for the pure magnetic case, [32] analyses the case of vacuum energy w = −1. The nature of the solutions depends on the values of various constants, it can collapse isotropically or anisotropically, only in the longitudinal or in the transverse direction towards the Big Bang. In general the magnetic fields accelerates expansion (or decelerates collapse) in the transverse direction of the magnetic pressure and it decelerates expansion (or accelerates collapse) in the direction of the magnetic tension. For general properties of the solutions, see [33].
Some interesting cases are analysed in [30]: if B 2 /ρ → 0 towards the singularity then the magnetic field effects are negligible; if B 2 /ρ does not approach 0, then it is constant and both fluids determine the dynamics, or the magnetic field causes a rapid expansion in the transverse direction and this change of the dynamics causes B 2 /ρ → 0. Moreover, [32] shows that in presence of a cosmological constant the magnetic field has a strong effect at early times, decelerating the collapse in the transverse direction and accelerating it in the longitudinal one, and is negligible at later times, when the vacuum energy causes accelerated expansion in both directions; the authors also describe the shape of the singularity.
It should be noted that in general the presence of the magnetic field causes a slowing down in the process of isotropization, making the shear more important; this way the CMB gives a strong constraint on primordial homogeneous magnetic fields [9,10].

Basic equations
We will now recap the fundamental equations we'll need later; their derivation can be found in [19]. Following [19] we define the magnetic field as the spatial part of the Faraday tensor F µν in the frame comoving with the cosmological fluid; we will use the ideal MHD approximation to turn off the electric field. These equations can be easily obtained in the covariant 3+1 formalism [34,35,36,37], as done in [38,39,19,40]; we will solve them, however, in a fixed synchronous gauge. We will assume geometric units for the speed of light c and Newton's gravitational constant G in witch c = 8πG/c 4 = 1.
We describe an anisotropic system with a metric g µν with positive spatial signature (−, +, +, +) filled by a perfect fluid with energy density ρ, isotropic pressure density p, 4-velocity u µ and energy momentum tensor where h µν is the comoving spatial projector and a uniform magnetic field with Faraday tensor F µν . The time derivative of a generic tensor T ν µ iṡ its spatial projected derivative the totally antisymmetric spatial tensor where η µνρσ is the totally antisymmetric tensor with η 0123 = 1/ √ −g, and the irreducible components of the velocity derivative are It is now possible to describe the electromagnetic field in the Lorentz-Heaviside units: the electric field is E µ = F µν u ν ; the magnetic field is B µ = ε µνρ F νρ /2, with magnetic energy B 2 = B µ B µ and energy momentum tensor The equations that describe our system are the Maxwell equationṡ where J µ is the electric 4-current, and the projected Einstein equations in which R µν is the Ricci tensor. The interaction between the fluid and the magnetic field is given by It is possible to use the Maxwell equation (8a) to find the conservation law for the magnetic energẏ and to derive the fluid energy conservation law from the temporal part of the Bianchi identities u µ ∇ ν T µν = 0 from the spatial projected Bianchi identities h µ ρ ∇ ν T ρν = 0 it is possible to find the momentum conservation law which, using gives

Background model
We assume that our system is homogeneous and perturbed at first order by weak inhomogeneous perturbations. At the background level we have a homogeneous universe with an isotropic perfect fluid and a uniform magnetic field: such field cannot live with an isotropic metric, such as FRW, but it can be accommodated in an anisotropic model. We must use one of the Bianchi models because of the homogeneity and our model fits best in a Bianchi I universe, which is the simplest anisotropic generalization of FRW, so our metric in synchronous gauge is These type of models were widely studied in literature in different assumptions and physical limits (see for example [41,30,31,25,32]); we are here interested mainly in their behaviour after the matter-radiation equivalence, where the magnetic field can be reasonably small compared to the matter component. This regime was already studied in different works, for example by [25] in radiation dominated universe; here we will recap [10], which accounts for different type of anisotropic stresses in both radiation an matter dominated universe. We will, however, amend for their time behaviour in matter dominated universe and we will not neglect higher order corrections in the isotropic components.
We assume that the magnetic field is oriented along the 3 axis, so the system is axisymmetric and a 1 = a 2 ; for simplicity we call a = a 1 = a 2 and c = a 3 . We have It is now straightforward to write the Einstein equations (9) 2ä and the energy conservation laws for the system (12) and (11) ρ + 2ȧ We define the Alfvén velocity, which is the energy ratio between magnetic field and fluid (note the factor 1/2 which differs from the usual definition) witch is responsible for the intensity of the anisotropies, the isotropic expansion H and the anisotropy parameter S If we now assume a barotropic fluid with equation of state p = wρ and w = const the Einstein equation (17a) becomes subtracting equation (17c) from equation (17b) we get and summing 2 times equation (17b) to equation (17c) we eventually have From the definition of v 2 A (20) and from the energy conservations (18) and (19) we havė If we now assume that the magnetic field energy is small compared to the fluid energy we have v 2 A 1 and if we write A it is easy to see from equations (22) and (24) that at 0-order in v 2 A we recover FRW and we have The anisotropy is described by S and equation (23) becomes at first order in v 2 while equation (25) giveṡ The isotropic part is contained in equations (22) and (24), which form a system whose solution is We are interested only in anisotropies caused by the magnetic field so we will put to 0 the homogeneous solution of each equation, with the exception of (29).

Radiation dominated universe
For radiation dominated universe w = 1/3 and equation (29) gives Equation (28) then gives From equation (30) we get ρ. From the definitions (21) we can get the values of a and c. Finally we have

Matter dominated universe
For matter dominated universe w = 0 and eq. (29) gives From equation (28) we get For the isotropic part we proceed as before: eq. (31) gives From equation (30) we get ρ.
From the definitions (21) we can get the values of a and c. Finally we have

Perturbed equations
We perturb all the quantities that govern our system while keeping synchronous gauge, thus the perturbed metric is where B means the background value; we can define where the indices of γ µν are raised and lowered with the unperturbed metric g B µν . In the following we write the trace of γ µν as γ = γ k k . The fluid velocity perturbation is δ u µ , with The fluid energy perturbation is δ ρ and the fluid pressure but we keep v 2 S as an arbitrary function and possibly different from w; the reason of this choice will be clear in section 7.2.
The perturbed magnetic field must remain pure spatial at all orders, as shown in appendix Appendix A, so the condition B µ u µ = 0 holds at all perturbative orders and the perturbation to the magnetic field satisfies Accordingly to [24,16] the perturbed Christoffel symbols are and the perturbed Ricci tensor is We are now ready to perturb the exact equations of section 3. We notice that, because of the homogeneity of the background model, when applied to the perturbation of a scalar quantity the comoving time derivativeṡ is the same as the synchronous time derivative ∂ 0 s, so we make no difference between them in the following. The fluid energy conservation (12) becomeṡ and the magnetic field energy conservation (11) giveṡ The Einstein 00 equation is (we will always use Einstein equations with a lower and an upper index) while the 33 equation reads to remove ∂ 3 ∂ k γ 3 k from the last equation we need to use the derivative of the 03 equation with respect to the 3 index If we had used equations (9) we would have found the same results.
By imposing the null divergence of the magnetic field (8d) we get The last equation we need is the conservation of the momentum (15) (note that A µ has only the first order component): we define an index P ∈ {1, 2} that lies on the plane orthogonal to the background magnetic field and we write the divergence of the momentum conservation on the 12-plane and the derivative of the 3 component along the 3 axis The system (54)-(61) fully characterizes the evolution of the perturbed quantities and it is the ground of the following analysis. Compared to [22] we fully related the background anisotropy to the magnetic field, without the need of additional hypothesis.

Gauge Modes
Fixing the synchronous gauge does not end the freedom of coordinate choice: we can still make a gauge transformation preserving the synchronous gauge.
We follow the same scheme as of [26]: we make a generic coordinate transformation of the form with small ε µ and we keep terms up to O(ε).
The metric tensor becomes If we define to preserve the synchronous gauge we need ∆ g 0µ = 0 which gives ε 0 = ε 0 (x j ) and where ε 0 (x j ) andε i (x j ) are arbitrary functions of the spatial coordinates: we still have 4 unused degrees of freedom represented by the functions ε 0 andε i . If we take the functions ε 0 andε i of the same order of the perturbations then the transformation given by equation (64) can be seen both as a gauge transformation and as a transformation of the functions γ µν within fixed synchronous gauge: in the latter case equation (64) gives the value of ∆ γ µν . In the same way the stress-energy tensor transforms as and if we see these as transformations on the physical variables instead of the coordinates we obtain the gauge modes for δ T µν . Substituting the explicit expression of T µν as the sum of the fluid and the magnetic field components we see that the transformation acts separately on the two components and we get for the fluid density perturbation In section 5 we linearised the equations and so the gauge transformations solve our equations and we call them gauge perturbations or gauge modes: these solutions are not physical because they correspond to a simple change in the reference frame. We are looking for a physical solution for the time dependence of δ ρ so the most interesting gauge transformation is given by equation (67).
Having the knowledge of gauge modes it is possible to construct gauge invariant variables, in a similar way as done in [42]. We have so our main scalar variable should be It is easy to check that it is exactly the variable used in [19], expressed in synchronous gauge. We will, however, not need it because the vorticity part H∂ i δ u i decays in time with respect to ∂ i ∂ i δ ρ /ρ B and we are interested in late time dynamics. We will also not need the laplacian, because we'll use Fourier expansions so it will reduce to a multiplicative term: for late times we can assume δ ρ to be gauge invariant.
It is possible to watch this approximation from another perspective, shown in Appendix Appendix B.

Analytical Solutions
If we write the perturbations as Fourier transforms we see that the system imposes different evolution to the perturbations that propagates along the background magnetic field, with ∂ P (. . . ) = 0, and the perturbations that propagates orthogonally to the background magnetic field, with ∂ 3 (. . . ) = 0. These different modes are however coupled by the magnetic stress-energy tensor tensorial nature.
To simplify the equations we use the barotropic state equation for the fluid, so p B = wρ B with w = const and δ p = v 2 S δ ρ, and the Fourier expansion for the spatial part of the perturbations, so the spatial dependence is of the form e ik j x j . We define the new variables Our differential equation system is not simple but we can solve it for small magnetic fields by keeping only terms up to first order in v 2 A : we shall remember that S is already at first order while ∆ , G, T and δ u i have also a 0-order (FRW) part; M has only the 0-order part because it is always multiplied In the same way, looking at our system also T is always multiplied by v 2 A : this is because it does not affect density perturbations unless some anisotropy is present. We also use eq. (50b) to discard terms proportional to w − v 2 S or to(v 2 S ), unless multiplied by k i k i or k 3 k 3 . This is because, while they are equal to 0 for w = const, we will need them in sec. 7.2.
The fluid energy conservation equation (54) in the new variables readṡ Similarly we rewrite the magnetic energy conservation (55) where we found the last equality by using the fluid energy conservation.
Combining Einstein 33 equation (57) with its derivative with respect to time and using the derivative of Einstein 03 equation (58) with respect to the 3-index in order to take care of ∂ i ∂ 3 γ i 3 terms we get an equation for T . Because T only appears in the system in terms that are multiplied by v 2 A , we will only need this equation at 0-order: We can use the fluid energy conservation equation (74) to eliminate G from the other equations. This way the Einstein 00-equation (56) reads We obtain the evolution equation for the divergence of the 4-velocity by summing eqs. (60) an (61); we then use equation (57) to remove the ∂ i ∂ 3 γ i 3 term and equation (59) to remove the divergence of the magnetic field. Doing so we We will need also equation (61) which reads, using equation (51a) to remove ∂ 3 δ B 3 , Thus we restated the dynamical system (54)-(61) in a more suitable form which is more appropriate for the following analysis.

Radiation dominated universe at large scales
In radiation dominated universe we have w = v 2 S = 1/3 and at large scales we can set k 2 ≈ k 3 k 3 ≈ 0. It is easy to check that, once we get rid of the scale dependent terms, eq. (76), (77) and (78) reduces respectively to This system, together with (75) and (79), is satisfied by a power law solution and could be reduced to a pure algebraic problem, but we found simpler to solve it for v 2 A = 0 and then look perturbatively for the corrections in v 2 A . We found It can be shown that the t 1/2 modes are related to a nonvanishing divergence of the background velocity ∂ i δ u i = ik i δ u i : strictly speaking, we should have imposed the k i ≈ 0 condition, thus finding only the t and 1/t modes: and recovering the usual FRW solution in the limit v 2 A → 0. Using (67) and (37) we find that 1/t is a gauge mode, while t 1−2v 2 A0 is the physical growing mode, with the correction due to the magnetic field.
We find our solution simpler than the one of [19], and with a clearer physical interpretation of the solutions, but our physical growing mode follows a slightly different temporal law, although this correction is small given v 2 A 1. We also find simpler the comparison of our solution with the non magnetic one of [16].
We see in (84) that the magnetic field reduces the growing rate of density perturbations, but by an amount of order O v 2 A 1. This effect has long been known, and it is due to the extra magnetic pressure. A similar behavoiur was found in [19] and [22], although with the differences stated above.
We finally note a difference between our solution (83) and the one of [19]: the non dominant mode is t 1/2 in our formalism, while t −1/2 in their. At a more careful analysis, our equations tend correctly to the ones of [16] for v 2 A → 0 and we obtain in such limit the same solutions of [17,26], including the t 1/2 mode. Such discrepancy is therefore between the synchronous and covariant formalisms, and it is besides the purposes of our paper.

Matter dominated universe at small scales
In this section we analyse the perturbations in a matter dominated universe (w = 0), in the regime in which the anisotropies are small with respect to the background. We expand in Fourier the spatial part of each quantity like e ik j x j , with k j = const, and we define k 2 = k i k i .
Being at small scales means k 2 H 2 and assuming v 2 S , v 2 A 1 we can greatly simplify our equations, keeping only terms in v 2 S or v 2 A that are multiplied by k 2 and dropping terms of order v 2 S and v 2 A . This means that the effect of the sound speed and the Alfvén speed is relevant only at very small scales, as we will see from the solutions of our equations. This approximation, although still relativistic and so comparable to other result in literature, for example [19], will give the nonrelativistic limit, as shown in section 7.2.2

Sound speed and Alfvén speed
First we need some considerations regarding the sound speed. From a formal point of view, the sound speed is re-lated to the barotropic index w by (50a) and w = const implies v 2 S = w, so it should vanish. From a physical point of view we need a nonvanishing sound speed and we can also estimate its value. While formally the best solution to this problem would be using a two fluid model, with a different equation of state for perturbations, here we will simply drop the relation between v 2 S and w and assume that the perturbed fluid follows a different equation of state with respect to the background fluid. This is correct in the Newtonian approximation and it's in fact the standard way of handling things [16,6], while putting v 2 S = 0 at the end will recover the full covariant value of our calculations for studying pure magnetic effects.
We proceed as in [16]: we use an adiabatic sound speed where γ is the heat ratio. We write We can estimate more precisely the sound speed value, and it's possible to show that the adiabatic sound speed is [16,6] v 2 S z<z rec = 1 3 where z rec is the redshift value at recombination, k b is the Boltzmann constant, T b is the baryons temperature, m p is the proton mass and σ is the specific entropy, whose value is σ = 4a SB T 3 /3n b k B ≈ 1.5 · 10 9 , being a SB the Stefan-Boltzmann constant and T the gas temperature. We neglected any anisotropic effects in temperature, because they would be related to the next order corrections. The baryons temperature is the same of the photons until z ≈ 100, due to residual Thomson scattering, and decreases faster thereafter: We define two constants addressing the effect of sound speed and Alfvén speed after recombination. Taking the time dependence of k 2 depending only on the 0-order part of the background metric because it always appears multiplied by v 2 S or v 2 A , we have respectively For a more detailed discussion about the sound speed see [43].

Analytical solutions
Using the assumptions of section 7.2 we can greatly simplify our equations. The energy conservation (74) and the magnetic field energy conservation (75) retain the same form. The Einstein 00-equation (77) now reads The momentum conservation 78 becomes and its counterpart along the z-axis remains (79): We need the Einstein 33-equation only at 0-order in the magnetic field, after being multiplied by v 2 A , so equation (76) in our limit reads With some algebra it is possible to reduce this system to a single equation. Expanding the spatial part in Fourier, defining the anisotropy parameter µ of the solution as and using the constants (89) we find, after some algebra, where ∆ (i) is the i-th derivative of ∆ . This corresponds exactly to equation (29) of [6], except for a difference in the definition of v 2 A and so in Λ A . We believe interesting to analyse separately the two cases of ν = 0 and ν = 1/3, instead of studying them together as in [6].

Post recombination evolution
For 1100 > z > 100 we have ν = 0. The solution of (95) is where ∆ i are arbitrary constants and The only possible growing solution is x 3 , and the requirement is that it holds one of the conditions using (89) and (26), making explicit the presence of Newton's constant we get ρ = 1/6πGt 2 , conditions (98) become While the first one is the standard Jeans condition, the second one means that, orthogonally to the background magnetic field, there is a heavier requirement dependent on the strength of the magnetic field: some modes could grow in every direction but the one of the field. The presence of the magnetic field also imposes a slowing down of the growing mode: where the equal sign holds only in absence of a magnetic field, that is only if Λ A = 0.

Late times evolution
This is exactly the case analysed in [6]. For z < 100 we have ν > 0 and the solution of (95) is where ∆ i are arbitrary constants, 2 F 3 is a generalized hypergeometric function with constant coefficients a i j , b i j depending only on the constants ν, Λ S , Λ A (see app. Appendix C for the explicit value of the coefficients) and The solutions can grow only if the argument of the hypergeometric functions is small, i.e. if Λ 2 S /4ν 2 t 2ν 1 : this way we have Condition (103) is the standard Jeans condition [16]: using (89) and (26), eq. (103) translates in [6] The only solution in (104) that can grow is 3: The first one means that, in any direction but orthogonal to the background magnetic field, the only necessary condition is the standard one. The second one is an additional condition that must hold for perturbations propagating orthogonally to the background magnetic field, and using (89) and (26) it reads [6] k The presence of this new condition makes possible the existence of Jeans unstable modes, that orthogonally to the background magnetic field are stabilized by the magnetic pressure if k A < k J and k A < k < k J [6].
Studying the growing rate of this solution with more care, we see that x 3 satisfies orthogonally to the background magnetic field the growing rate is unchanged, while in other directions it is slowed down, depending on the field strength.

Full relativistic case
If we put v 2 S = 0 we recover the exact relativistic solution. As we can see from the previous solutions, the growing condition is Moreover, the solution is with x i given by (97) with Λ S = 0, or equivalently by (102).
If we compare our result with [19], we identify the anisotropic behaviour and we obtain the correct Newtonian limit of [6]. However, our solutions are different and we are unable to explain such discrepancy: we can argue they may have found some sort of average effect, however this is not clear, given the strong anisotropy of the model: the magnetic Jeans wavenumber is present only in one direction, the one with µ = 0.
In case Λ 2 A 1 we have and setting µ 2 = 1/3 the solutions x 3 and x 4 recover eq. (31) of [20] and eq. (31) of [21], so our small scales solution of sec. 7.2 is a generalization of their work, while including a nonvanishing sound speed and pressure.

Numerical integration
To better show our results, we numerically integrated the system (54)-(61), using estimates from [44] to set the numerical values for the background functions. We followed the same procedure of [6] to determine the initial conditions: we started the integration from a very early time and we verified that the initial perturbations were outside the Hubble horizon and we used the large scale solution to match the initial conditions to the growing mode; in our case such conditions come from eq. (84). We assumed to perturb only the baryon component of the universe, while leaving the CDM component unperturbed; a rigorous treatment should rely on a multifluid model, but we ague that we can still extract meaningful information within our approximation. Practically speaking, this assumption means that every quantity present in our equations at perturbative level must be replaced by its baryonic component, while the background model still depends on CDM. Our equations are still correct, because the background interaction is only due to energy density, while at perturbative level every dependence on CDM disappears, except from background quantities.
We choosed to study the same scales of [6], i.e. k ≈ (17, 1.7, 0.37) Mpc −1 normalized at present time, corresponding to baryonic masses of M ≈ (1.5 × 10 8 , 1.5 × 10 11 , 1.5 × 10 13 ) M and roughly equivalent respectively to a dwarf galaxy, a galaxy and a galaxy cluster. The results of the numerical integration are shown in figure 2.
Our results must be compared to the ones of [6]. Until equivalence (z ≈ 3400) we are in radiation dominated universe and the comparison is obvious: our solutions grow, while theirs decay; this is because in [6] the authors always consider matter dominated universe.
After equivalence, in both cases we are subject to a decaying period, followed by a new growth after recombination, but in our case this happens for a shorter time; most of the anisotropic effects comes in this era, because before equivalence the thermal pressure is much stronger than the magnetic one and most of the anisotropy is suppressed, so they are less relevant in our simulations. This is clear in fig. 2b, where we see almost no anisotropy. As a further confirmation, it can be shown that ∆ (z ≈ 10)/∆ (z ≈ 1100) has the same value in both the analysed cases, so the main anisotropic contribution comes from the region 3400 z 1100.
After recombination we have a behaviour similar to [6], because here we are at scales were the Newtonian approximation is correct. The apparent discrepancy in the oscillating behaviour of fig. 2a is mainly due to the (small) difference in the numerical values of the background functions, because the oscillating behaviour is very sensible to such   numbers; however, the qualitative evolution is the same, with the µ = 0 case beginning to decay because of the magnetic Jeans length [6] (eq. (99b) and (107)). In this region our solution has a slightly faster growth than [6], we argue this may be caused by some residual relativistic effects, but it has to be investigated with more care. Moreover, in the last region we should be outside of the linear regime, so we would need a full nonlinear treatment.

Conclusion
We developed above a self-consistent scheme for the analysis of cosmological perturbations in the presence of a magnetic field. We set up in the synchronous gauge a dynamical scheme which accounts for the effects induced by the magnetic field both on the background and the first order formulation. To this end, we considered a Bianchi I model, whose anisotropy with respect to the flat Robertson-Walker geometry is due to the privileged direction defined by the magnetic field.
We first solve in detail the equations describing the anisotropic background and then we analyse the perturbation dynamics, having awareness of the gauge contribution analytical form.
We amended for the previous analysis in [19] in the case of a super-horizon wavelength of the perturbation. In particular, our solution has a clearer comparison with the non magnetic one. We recovered the slowing-down of the growing mode caused by the magnetic pressure, and so of order O v 2 A 1. This effect has long been known in FRW models and has been analysed in Bianchi I models with particular anisotropies by [22], while we worked always relating the background anisotropy to the magnetic field without additional assumptions.
We refined the results of [19] for the sub-horizon wavelength of the perturbations, showing that an anisotropic treatment is required. We also generalised the results of [20] and [21], while including a nonvanishing sound speed and considering the anisotropic case.
We finally enforced the Newtonian limit obtained in [6], completing it with the relativistic analysis, also facing a numerical treatment. We showed that the relativistic regime limits the anisotropy induced by the magnetic field.
Overall, despite the assumption of a Bianchi I background, most of our solutions reproduce those obtained on an FRW background. At a closer look, the Bianchi I anisotropy enters the system via the S function defined in (21). At small scales the relevant terms are the ones with k 2 , and none of those are related to such anisotropy. However, when the condition H 2 k 2 does not hold, such terms become important; unfortunately, in this case the system would be much more complicated that the one of sec. 7.2. On the other hand, at large scales the background anisotropy survives, and we argue that it is mainly related to the perturbed fluid velocity. In particular, it can be shown that the solutions proportional to t 1/2 in (83) are related to δ u i , and more precisely in ∆ 2 t 1/2+4v 2 A0 we have both k i δ u i = 0 and k 3 δ u 3 = 0, while in ∆ 1 t 1/2−2v 2 A0 it holds k 3 δ u 3 = 0; the solutions ∆ grow t 1−2v 2 A0 and ∆ gauge /t, on the other hand, both have k i δ u i = k 3 δ u 3 = 0.
We stress that, in order to solve the equations, we assumed a small magnetic field and so all the effects we studied are related to v 2 A 1, and they become relevant only at small scales, due to the large wavenumber k 2 H 2 and to the also small sound speed v 2 S . This is clear by looking at fig. 1.

Appendix A: Magnetic field at perturbative level
In literature there are different definitions of the magnetic field at a perturbative level, but it is easy to recognize that not all of them satisfy the required properties. After a careful analysis we concluded that the correct one, at least with respect to the physical phenomenon we study here, it the one of [19] made through the 3+1 formalism. This way, the magnetic field is defined as the spatial projected part of the Faraday tensor F µν , while the electric field as the temporal one at all orders. There are two important reasons for this requirement. The first one is that the electromagnetic field is decomposed in electric and magnetic components by the observer and we are interested in its interaction with the cosmological fluid, so the natural observer is the fluid itself. Beside that, we force a vanishing electric field E µ = 0 through the assumption of infinite conductivity of the medium, thus we work in the limit of ideal MHD. To do this we need these fields to be defined with respect to the fluid. Using this definition there are no induced fields, reflecting the fact that the covariant form of Maxwell's formulae and of the electric and magnetic field definitions already incorporates the effects of relative motion [19].
The second reason is that with different definitions we would have a nonvanishing trace for the perturbed magnetic stress energy tensor, while this way all goes well and it is traceless. This is easy to check using the definition of perturbations from section 5.

Appendix B: Gauge behaviour in late times
We will analyse here the FRW case, to clarify the meaning of δ ρ becoming gauge invariant for late times. Following [16] and using the Newtonian approximation we see that the solutions after recombination are where γ = ν + 4/3 > 4/3 is the heat ratio of the fluid (after recombination γ 5/3), δ = δ ρ /ρ, is a constant, v 2 S is the squared sound speed and k the wavenumber. The functions J a (z) are the Bessel functions: when their argument is large they oscillate, but when the argument is small they behave like δ ± ∝ t (−1±5)/6 . (B.5) The growing mode is the physical solution we are looking for, while the other one decays to zero.
We cannot speak of gauge modes in Newtonian theory, but the decaying mode corresponds exactly to the relativistic gauge mode, and as expected it decays in time with respect to the growing one. This means that, for large times, gauge modes naturally decay to zero and we can neglect them as long as we are looking only for the growing ones.
It should be noted that in our calculations we are in the same situation: we cannot have a relativistic sound speed different from w in a single fluid model, but we make this approximation in section 7 because from a physical point of view we need a nonvanishing sound speed. This way we "break" the gauge invariance, but the gauge modes manifest themselves in one of the decaying solutions. We are only looking for growing modes, so we can safely neglect them.

Appendix C: Late times solution coefficients
We report here the values of the coefficients of the hypergeometric function appearing in (101), using δ ± defined in eq. (102c):