Tight coupling expansion and fully inhomogeneous magnetic fields

The tight coupling expansion, appropriately generalized to include large-scale magnetic fields, allows the estimate of the brightness perturbations of CMB anisotropies for typical wavelengths that are larger than the Hubble radius after matter-radiation equality. After discussing the basic features of the the pre-decoupling initial conditions in the presence of fully inhomogeneous magnetic fields, the tight coupling expansion is studied both analytically and numerically. From the requirement that the amplitudes and phases of Sakharov oscillations are (predominantly) adiabatic and from the inferred value of the plateau in the temperature autocorrelation, the effects of the magnetized contribution can be systematically investigated and constrained.


Introduction
A distinctive signature of pre-decoupling (adiabatic) initial conditions is that the temperature autocorrelation has the first (Doppler) peak for ℓ d ≃ 220 [1,2,3]. For the same set of initial conditions, the cross-correlation power spectrum between temperature and polarization will have the first (anticorrelation) peak for ℓ c ≃ 3 ℓ d /4 < ℓ d [1].
Defining as k the comoving wave-number and as c sb the sound speed of the baryonphoton system, the large-scale contribution to the temperature autocorrelation oscillates as a cos (kc sb τ dec ), while the cross-correlation oscillates as sin (2kc sb τ dec ). The first (compressional) peak of the temperature autocorrelation will then correspond to kc sb τ dec ∼ π (i.e. ℓ d ≃ 220), while the first peak of the cross-correlation will arise for kc sb τ dec ∼ 3π/4 (i.e., approximately, ℓ c ≃ 3 ℓ d /4). This result can be obtained analytically by working to firstorder in the tight-coupling expansion [4,5,6]. In this framework, the key assumption is that the initial conditions of the Einstein-Boltzmann hierarchy are characterized, prior to matter-radiation equality, by a solution that is, at least predominantly, adiabatic. A corollary to the previous statements is that the peak in the cross-correlation occurs for typical wavelengths that are larger than the Hubble radius at the the time of photon decoupling.
In a previous paper [7] it has been shown that large-scale magnetic fields, present prior to equality affect the CMB initial conditions. If only the adiabatic mode is present, the fully inhomogeneous magnetic fields contribute to the Sachs-Wolfe plateau as a subdominant non-adiabatic component that may be either correlated or uncorrelated with the adiabatic mode. If the pre-decoupling initial conditions are given by a mixture of adiabatic and nonadiabatic modes, the magnetized contribution may also mix with the isocurvature modes so that the total number of parameters defining the initial conditions increases. It is desirable to complement and extend the analysis of [7] by scrutinizing the impact of pre-decoupling magnetic fields on intermediate scales. Within the tight coupling expansion, appropriately generalized to include the effects of fully inhomogeneous magnetic fields, the fluctuations in the brightness induced by the (adiabatic and non-adiabatic) scalar modes can be estimated.
The effects of magnetic fields on the scalar brightness perturbations are typically neglected. To give an example, let us consider how the impact of large-scale magnetic fields is usually computed, for instance, in the Faraday effect. This case is interesting, since, the rotation of the CMB polarization plane indeed presupposes that the magnetic fields do not affect the process of formation of polarization. The influence of large-scale magnetic fields on CMB polarization can be separated, for illustration, into two physically distinct effects. Prior to decoupling, gravitating magnetic fields modify the initial conditions of CMB anisotropies and affect the evolution equations of the photon-baryon system. Therefore, the induced degree of polarization will be sensitive to the presence of magnetic fields whose role will be, at this stage, to modify the features of the adiabatic initial conditions. This is the first effect. Once the polarization is formed, the polarization plane of the CMB may be Faraday rotated [8,9]. From this second effect interesting bounds on the magnetic field intensity are usually derived.
This philosophy motivated various independent studies on the possible Faraday rotation of the CMB polarization degree [10,11,12,13]. In some cases the magnetic field has been assumed just uniform [10], but there are also calculations analyzing the situation where the magnetic field is fully inhomogeneous [11,12]. In spite of the relevant technical differences, the common assumption of the mentioned investigations is that the magnetic field is negligible in the process of formation of CMB polarization. The latter assumption can be rephrased by saying that the initial conditions of the Einstein-Boltzmann hierarchy are set in such a way that the magnetic fields are absent. In the absence of large-scale magnetic fields the only source of polarization arises to first-order in the tight coupling expansion. The magnitude of the obtained contribution is just a consequence of the adiabaticity of the solution. If largescale magnetic fields are consistently introduced before matter-radiation equality, the usual adiabatic (as well as non-adiabatic) initial conditions of the Einstein-Boltzmann hierarchy are modified and, consequently, also the degree of polarization of the CMB is affected.
This example motivates a systematic analysis of the impact of magnetized initial conditions on the brightness perturbations of the radiation field. As anticipated above, this calculation will be conducted within the tight-coupling approximation which is known to be rather effective in obtaining the brightness perturbations for large and intermediate scales. This analysis has never been done before, to the best of our knowledge. A related motivation is that, in the absence of magnetic fields, the tight coupling approximation is employed in numerical Boltzmann codes at early times to avoid the integration of stiff Euler equations. The results of the present investigation allow to use the same strategy also in the presence of fully inhomogenous magnetic fields.
Several aspects of CMB physics are germane to large-scale magnetization (see [14,15] for two topical review articles). The astrophysical evidence for large-scale magnetization can be summarized as follows. Large-scale magnetic fields are an important element of the physics of the interstellar medium and have been measured, through various techniques in galaxies [16,17] and in rich clusters [18,19,20,21]. The determination of large-scale magnetic fields associated with loosely gravitationally bound systems (like superclusters) is still debatable [22,23,24], but, nonetheless, extremely intriguing. As far as galactic magnetic fields are concerned, there is the hope that in the future all sky survey of Faraday rotation (achievable through the Square Kilometer Array [SKA] [17,18]) will allow a Faraday "tomography" of the Milky Way addressing the dark corners (number of reversals, pitch angle) of the morphological structure of the closest large-scale magnetic field we observe through different surveys of pulsars such as the Parkes and the southern galactic plane surveys (see [25,26,27] and references therein). There is the hope, in the present framework, that the present-day coherent field in spiral galaxies may keep the "memory" of the structure of the initial seed field (see [16] for a review of this long-standing problem). The morphological features of the magnetic fields profiles in clusters are less understood than in the case of galaxies. However, attempts have been made in constraining the radial profiles of magnetic fields in clusters [19,20]. Determinations of intracluster magnetic fields through statistical samples of rich clusters [21] seem to suggest that (ordered?) intra-cluster magnetic fields may be as large as the µG (see, for instance, [28] for a review on the features of our magnetized Universe).
Compressional amplification (taking place during the gravitational collapse of the protogalaxy) allows to connect the observed field to a protogalactic field, present prior to gravitational collapse, between 0.1 nG and 0.01 nG. It is plausible, within the dynamo hypothesis, that the protogalactic fields could even be much smaller than the nG. In recent years, a lot of progress has been made both in the context of small-and large-scale dynamos i.e. the process by which the kinetic energy of the protogalactic plasma is transferred, by differential rotation, into magnetic energy. This progress [29] (see also [30] for an introduction to astrophysical dynamos) has been driven both by the higher resolution of numerical simulations and by the improvement in the understanding of the largest magnetized system that is close to us, i.e. the sun. In contrast with the previous lore, it is now clear that the dynamo action demands a change in the topology of the flux lines. As a consequence, large-scale dynamos should also produce small-scale helical fields that quench (i.e. prematurely saturate) the α effect [29] (see also [31]). The morphology and strength of magnetic fields in clusters may be related to the way the dynamo effect saturates.
The simplest way to understand why large-scale magnetic fields may be relevant for CMB physics is to think of a protogalactic magnetic field in the nG range. If this field is naively blue-shifted at the decoupling epoch, i.e. at a redshift z dec ≃ 1100 its strength could be as large as the mG, i.e. roughly six orders of magnitude larger. In the past ten years, indeed a lot of work has been done with the purpose of constraining large-scale magnetic fields using CMB physics. Historically the first type of configurations used in this type of exercise have been weakly uniform magnetic fields that would partially break the isotropy of the background geometry [15]. It was than realized that this assumption induces a correlation of the a ℓ−1,m and a ℓ+1,m multipole coefficients arising in the angular power spectrum (see [32,33,34] and references therein). From this observation, uniform magnetic fields can be constrained.
Another (complementary) possibility is that magnetic fields are fully inhomogeneous. In this case the magnetic fields do not break the spatial isotropy of the background but they can affect virtually all observables relevant in the theory of CMB anisotropies. In particular, in a series of papers, Barrow and Subramanian [35,36,37] pointed out possible effects of magnetic fields on the induced vector CMB anisotropies as well as on the polarization power spectra induced by the same modes. Within a slightly different approach the vector modes (as well as the tensor modes) have been discussed in [38]. More recently [39] a full numerical analysis of the vector and of the tensor modes induced by fully inhomogeneous magnetic fields has been presented. In [40] the signatures of magnetic helicity have been studied always in the perspective of the vector modes. Finally, in a series of papers, the evolution of weakly inhomogeneous magnetic fields has been discussed within the covariant formalism and in the presence of a uniform magnetic background [41].
In spite of the fact that vector modes are probably one of the characteristic signatures of large-scale magnetic fields, scalar modes are also induced [7,42]. This statement can be understood by noticing that the magnetic energy density, the magnetic anisotropic stress and the magnetic pressure contribute, respectively, to the Hamiltonian and momentum constraint and to the dynamical evolution of the curvature perturbations. Furthermore, within the magnetohydrodynamical approximation, the divergence of the Lorentz force contributes to the evolution equation of the divergence-full part of the baryon peculiar velocity.
What is missing, at the moment, is the analysis of the brightness perturbations induced by magnetized initial conditions. As anticipated, one of the purposes of the present paper will exactly be to bridge fill this gap and to allow an approximate evaluation of the brightness perturbations, also at intermediate scales, as a function of the spectral properties of fully inhomogeneous magnetic fields.
The present paper is then organized as follows. In Section 2 the initial condtions of the Einstein-Boltzmann hierarchy will be discussed for typical wavelengths larger than the Hubble radius prior to equality. In Section 3 the tight-coupling equations will be derived and discussed analytically. The numerical treatment of the system is presented in Section 4. Section 5 contains our concluding remarks.
2 Adiabatic and non-adiabatic modes 2.1 Zeroth-order tight-coupling approximation By combining the Thompson cross section (defined, in the following, as σ T ) with the ionization fraction x e and with the electron density n e , the differential optical depth is given by κ ′ = n e x e σ T a/a 0 . Since prior to decoupling the photon mean free path is much smaller than the Hubble radius at the corresponding epoch, the inverse of κ ′ can be used as an expansion parameter. This is the basic idea of the tight coupling approximation. The lowest order, the tight-coupling approximation implies 1/κ ′ → 0 (or, equivalently, that σ T → ∞). In this limit the divergences of the peculiar velocities of photons (i.e. θ γ ) and baryons (i.e. θ b ) are driven to a common value i.e. θ γ ≃ θ b = θ γb . In the following, the consequences of the zeroth-order in the tight coupling approximation will be investigated in the presence of fully inhomogeneous magnetic fields. Consider, therefore, the case where large-scale magnetic fields are fully inhomogeneous and characterized, in Fourier space,by the following two-point function In Eq. (2.2) k p is the "pivot" scale at which the spectra are normalized (see also, below, Eq. (2.35)). Since the magnetic fields are fully inhomogeneous the line element of the background geometry will not be anisotropic 2 (as it happens, for instance, in magnetized Bianchi-type backgrounds) and it can be parametrized as Consequently, the evolution equations of the background geometry will be solely determined by the total pressure and energy density of the system as where H = a ′ /a and the prime denotes a derivation with respect to the conformal time coordinate τ . Denoting the longitudinal fluctuations of the metric (2.3) as 3 and selecting the conformally Newtonian gauge [53,54], the evolution equations of photons and baryons can be obtained in terms of the density contrasts (i.e. δ γ and δ b ) and in terms of the divergence of the peculiar velocities (i.e. θ γ and θ b ): The term J × B parametrizes the Lorentz force contribution as it arises in a conducting plasma which is, however, globally neutral (see, for instance, [45,46]). In this situation, typical of magnetohydrodynamics (MHD), the possible contribution of electric fields is suppressed by inverse powers of the conductivity as it will be discussed later. By subtracting the two equations for the velocities introduced in (2.6) and (2.7), the combination (θ γ − θ b ) obeys The term at the right hand side of Eq. (2.8) can be viewed as a source term for the evolution of the difference between the divergences of the photon and baryon peculiar velocities. In spite of the magnitude of the source term at the right hand side of Eq. (2.8), (θ γ − θ b ) will quickly be driven to zero in the limit σ T → ∞. Consequently, even if pre-decoupling magnetic fields are present, the lowest order tight coupling expansion still implies that θ γ → θ b as expected.
In the tightly coupled limit, the evolution equation of θ γb = θ γ = θ b can then be swiftly obtained by summing up the equations for the velocity fields in such a way that the scattering terms cancel. The result of this procedure is, in Fourier space, where R b = (3/4)ρ b /ρ γ . For numerical purposes (see Section 4) the late-time cosmological parameters will be fixed, for a spatially flat Universe, as 4 where ω X = h 2 Ω X and Ω Λ = 1 − Ω m ; the present value of the Hubble parameter H 0 will be fixed, for numerical estimates, to 73 in units of km/(sec Mpc). In terms of the fiducial set of parameters of Eq. (2.10) the baryon to photon ratio appearing in Eq. (2.9) becomes where z + 1 = a 0 /a is the redshift (for numerical purposes the scale factor will be normalized in such a way that a 0 = 1; see below, Eq. (2.25)). In Eq. (2.9) we traded the MHD current J (appearing in the second equation of (2.7)) for a combination of the magnetic energy density and of the anisotropic stress. More specifically, it can be easily shown that Equation (2.12) holds since, in MHD, the total Ohmic current J is solenoidal (i.e. ∇ · J = 0) and given by 4π J = ∇ × B [45,46] (see also the appendix of [28]). In addition, the function σ B is nothing but a parametrization of the anisotropic stress that appears in the spatial part of magnetic energy-momentum tensor, i.e. where In various context the force-free approximation is often employed [49,50] (see also [28] and references therein). Such an approximation amounts, in practice, to set ( ∇ × B) × B → 0. 4 The values of the cosmological parameters introduced in Eq. (2.10) are compatible with the ones estimated from WMAP-3 [2,3] in combination with the "Gold" sample of SNIa [47] consisting of 157 supernovae (the furthest being at redshift z = 1.75). We are aware of the fact that WMAP-3 data alone seem to favour a slightly smaller value of ω m (i.e. 0.126). Moreover, WMAP-3 data may also have slightly different implications if combined with supernovae of the SNLS project [48]. The values given in Eq. (2.10) will just be used for a realistic numerical illustration of the methods developed in the present investigation.
In this limit, from Eq. (2.12) σ B → Ω B /4. In the present framework this assumption shall not be invoked.
According to Eq. (2.12), the leading contribution of fully inhomogeneous magnetic fields to the scalar modes of the geometry is given by the magnetic energy density and by the divergence of the MHD Lorentz force. These two quantities, affect the perturbed Einstein equations and, in particular, the Hamiltonian constraint (2.15) and the total (scalar) anisotropic stress Equations (2.15) and (2.16) follow, respectively, from the (00) and (i = j) components of the perturbed Einstein equations written in the conformally Newtonian coordinate system; the quantity δ s ρ t is the (total) scalar fluctuation of the energy-density of the fluid sources. It is relevant to remark that, in Eq. (2.16), on top of the anisotropic stress of the magnetic field, the only source of the anisotropic stress of the fluid is provided by the neutrinos that are collisionless after weak interactions fall out of thermal equilibrium for typical temperatures of the order of 0.8 MeV. In analogy with Eq. (2.13), the anisotropic stress of the neutrinos has been parametrized as The quantity δ s p t is the fluctuation of the total pressure. The last constraint on the dynamical evolution us derived from the (0i) components of the Einstein equations and it is The second term at the right hand side of Eq. (2.18) is the divergence of the Poynting vector.
In MHD the Ohmic electric field is subleading and, in particular, from the MHD expression of the Ohm law Since the Universe, prior to decoupling, is a very good conductor, the the ideal MHD limit [46] can be safely adopted in the first approximation (see also [51]); thus for σ → ∞ the contribution of the Poynting vector vanishes. In any case, even if σ would be finite but large, the second term is suppressed in comparison with the first term at the right hand side of Eq. (2.18); such a term is given by the sum of the divergence of the peculiar velocity of each fluid component weighted by the specific enthalphy (i.e. (p a + ρ a ) for a generic a-fluid): 20) or, in more explicit terms, The second equality in Eq. (2.21) follows from the zeroth order in the tight-coupling expansion and from the definition of the baryon to photon ratio (2.11). Every regular solution to lowest order in the tight coupling expansion has to satisfy the momentum constraint of Eq.
(2.18). Such a requirement implies, in particular, a restriction on the form of the peculiar velocity fields for comoving wave-numbers larger than the Hubble radius (i.e. kτ ≪ 1). The evolution of the cold dark matter (CDM) component only feels the presence of the magnetic field through the modifications induced by Ω B and σ B on of the two conformally Newtonian potentials, i.e. φ and ψ: A similar situation occurs for the evolution of the neutrino component where, however, the effect of the magnetic field enters the not only the Newtonian potentials but also the anisotropic stress that is constrained from Eq. (2.16). The relevant equations for the neutrinos are then, in Fourier space, The quantity F ν 3 appearing in Eq. (2.23) is the octupole of the (perturbed) neutrino phase space distribution while σ ν (i.e. the anisotropic stress of the neutrinos introduced in Eqs. (2.16) and (2.17)) is related to the quadrupole moment of F ν as σ ν = F ν 2 /2. The rationale for the inclusion of F ν 3 in the fluid description is that neutrinos are collisionless for temperature smaller than 0.8 MeV when the Universe is still deep into the radiation-dominated regime. In this regime the neutrinos do not interact, and, consequently, they should be rather described by the collisionless Boltzmann hierarchy of the (perturbed) neutrino phase space distribution. To address this problem, the usual approach is to use an improved description where on top of the monopole and dipole of the (perturbed) neutrino phase space distribution (connected, respectively, with δ ν and θ ν ), one also takes into account the quadrupole, i.e. σ ν and sometimes (as in the case of specific non-adiabatic modes) also the octupole. With all these ingredients in mind, the evolution of the longitudinal fluctuations of the geometry is determined from the trace spatial components of the perturbed Einstein equations, and can be written as where δ s p t is the fluctuation of the total pressure introduced in Eq. (2.17). In the present situation it will be mandatory to set initial conditions for the fluctuations of the geometry and for the various species of the plasma before radiation-matter equality but, also, after where Ω m0 and Ω r0 are evaluated at the present time and the scale factor is normalized in such a way that a 0 = 1. In Eq. (2.25) τ 1 = (2/H 0 ) a eq /Ω m0 . In terms of τ 1 the equality time is Mpc, (2.26) i.e. 2τ eq ≃ τ 1 . In this framework the total optical depth from the present to the critical recombination epoch, i.e. 800 < z < 1200 can be approximated analytically, as discussed in [52]. By defining the redshift of decoupling as the one where the total optical depth is of order 1, i.e. κ(z dec , 0) ≃ 1, we will have, approximately

Magnetized adiabatic mode
In the limit τ /τ 1 ≪ 1, the solution of the whole system of governing equations corresponding to the adiabatic mode can be obtained as an expansion in powers of kτ . When the given wavelength is larger than the Hubble radius before equality (i.e. kτ < 1 for τ < τ eq )the solution of the coupled system introduced in the previous subsection can be written as 5 where the fractional contribution of photons to the radiation plasma, i.e. R γ has been introduced and it is related to R ν , i.e. the fractional contribution of massless neutrinos, as 6 In Eq. (2.31) φ * (k) and ψ * (k) denote the values of the Fourier modes of the longitudinal fluctuations of the geometry that are (approximately) constant in time. Equations (2.28), (2.29) and (2.31) define a solution that is adiabatic since the density contrasts in matter (i.e. δ c and δ c ) are exactly 3/4 of the density contrasts in radiation (i.e. δ γ and δ ν ). In the limit Ω B → 0 (vanishing magnetic energy density) and σ B → 0 (vanishing magnetic anisotropic stress), Eqs. (2.28), (2.29) and (2.31) reproduce the standard adiabatic mode deep in the radiation dominated epoch [42,57]. From the second equation reported in (2.31) it is clear that the neutrino anisotropic stress is partially counterbalanced by the magnetic anisotropic stress. This aspect is entirely due to the presence of the magnetic anisotropic stress and it has been verified in different situations involving also entropic fluctuations [58].
In view of the forthcoming numerical applications, it is useful to characterize the adiabatic mode in fully gauge-invariant terms. In fact the quantity introduced so far are meaningful in the longitudinal gauge. The spectrum of curvature perturbations used as initial condition for the numerical evaluation of the brightness perturbations is often introduced [1, 2, 3] in terms of the curvature perturbations on comoving orthogonal hypersurfaces [55,56], customarily denoted with R. In the longitudinal gauge we have that this quantity is simply given by In spite of the fact that R is here defined in terms of the variables of the longitudinal gauge, its value is invariant under infinitesimal coordinate transformations. In terms of R the longitudinal fluctuations of the geometry entering the solution of Eqs. (2.28), (2.29) and (2.31) become: where R * (k) is the constant value of curvature perturbations prior to matter-radiation equality. Within the notations followed in this paper, the spectrum of R * (k) is defined as (2.37) The power spectra of Eq. (2.37) follow from the definitions of Ω B and σ B by using the two-point function of the magnetic field, i.e. Eq. (2.1). In fact, according to Eqs. (2.12) and (2.14), Ω B and σ B are quadratic in the magnetic field intensity, and, therefore, lead to two mode-coupling integrals that have been estimated in in the nearly scale-invariant limit (i.e. ε < 1) where the functions F (ε) and G(ε) are determined to be In the opposite limit (i.e. ε ≫ 1) the mode-coupling integrals arising in the expressions of Ω B and σ B are dominated by (small scale) diffusive effects ( see, for instance, [35,38]) and the power spectra are sensitive to the specific scale of the diffusion (Alfvén) damping. In Eq. (2.37) the quantities have been introduced. Notice that B L measures the magnetic field intensity smoothed over a comoving length-scale L. In fact, by using a Gaussian window function e −k 2 L 2 /2 for each magnetic field intensity, the (real space) magnetic autocorrelation function at two coincident spatial points is In the present paper the calculations are consistently conducted within the parametrization introduced in Eqs. (2.1) and (2.2). Different authors use slightly different conventions which can be, however, easily translated in the present parametrization. For instance, the magnetic spectral index defined as n in [14,35,38,39] correspond, in the present notations, to ε − 3. In [14,35,38,39] the factor 2π 2 /k 3 appearing in Eq. (2.1) is omitted. These definitions are clearly conventional but we regard as preferable the notation adopted here since it agrees with the usual way of assigning the power spectrum of curvature perturbations in CMB physics. The other remark is that the magnetic fields considered here are not helical (see, for instance, [40] and references therein). This is not a limitation since it can be easily appreciated that the helical contribution does not affect the scalar modes which are the main subject of the present investigation.
where k L = 2π/L. Equation (2.40) implies that the amplitude of the magnetic power spectrum appearing in Eqs. (2.1) and (2.2) can be directly characterized in terms of B 2 L as It should be mentioned, incidentally, that nearly-scale invariant magnetic energy spectra can be achieved both in some class of inflationary models [59] and in the case of pre-big bang models [60]. It is not our purpose here to endorse any specific model. The spirit of this investigation is more modest since we shall be content of understanding if large-scale magnetic fields can be consistently included in the picture of pre-decoupling physics and in the framework of the usual techniques commonly employed for the description of scalar CMB anisotropies. On top of the curvature fluctuations on comoving orthogonal hypersurfaces, another class of gauge-invariant quantities is related to the density contrasts both total and partial (i.e. pertaining to each single species of the plasma). The total density contrast on uniform density hypersurfaces can be defined as a functional of the longitudinal degrees of freedom . (2.43) Having introduced ζ, related to the total density contrast, it is also useful to introduce the following 4 variables Equations (2.44) and (2.45) are a useful tool for the discussion of the non-adiabatic solutions of the system as it will now be explicitly shown.

Magnetized non-adiabatic modes
The inclusion of fully inhomogeneous magnetic fields not only affects the adiabatic mode but also the other (non-adiabatic) modes whose physical features have been investigated in the conventional case (see, for instance, [ (2.49) In Eq. (2.49) c 2 s i and c 2 s j are the sound speeds of each (generic) pair of fluids of the mixture. According to Eq. (2.48), the fluctuation in the total pressure may arise either because of a inhomogeneity in the energy density or thanks to the chemical inhomogeneity of the plasma. Chemically inhomogeneous means, in this context, that the plasma is constituted by various fluids with equations of state that are different and that entail necessarily a spatial variation of the sound speed. As a consequence of Eqs. (2.48) and (2.49), the evolution of ζ can be directly obtained from the evolution of δ s ρ t stemming from the first-order covariant conservation equation: Neglecting the Ohmic electric field we do obtain the wanted evolution of ζ, i.e.
(2.51) By using Eqs. (2.51) and (2.44) together with the explicit form of the evolution equations in the longitudinal gauge, the analytic form of the non-adiabatic contributions can be obtained. In the following, as an example, the magnetized CDM-radiation mode and the magnetized neutrino entropy mode will be presented. For the case of the CDM-radiation mode the solution, in the limit τ < τ 1 and kτ < 1 can be written as where In the case of this solution, the longitudinal fluctuations of the geometry are vanishing for τ < τ 1 and the relevant entropy fluctuations, i.e. S cγ and S cν are constant and proportional to S * . In the limit Ω B → 0 and σ B → 0 the solution (2.52) reduces to the conventional CDM-radiation mode. While in the standard case the octupole moment of the (perturbed) neutrino phase-space distribution vanishes, in the magnetized case it is proportional to kτ and is negligible for length-scales larger than the Hubble radius prior to equality. In analogy with Eq. (2.35), the spectrum of entropy perturbations will be parametrized as where k p is the same pivot scale introduced in Eq. (2.35) while n s is the spectral index of the entropy fluctuation. In the case of the magnetized neutrino entropy mode, the solution can instead be presented as: In the case of the neutrino-entropy mode S νγ = 0 and it is proportional to a constant that has been denoted in Eq. (2.55) asS * (k); the spectrum ofS * (k) can be defined in full analogy with Eq. (2.54). An interesting feature of the magnetized neutrino-entropy mode (common to the conventional case when magnetic fields are absent) is that the curvature perturbations on comoving orthogonal hypersurfaces (i.e. R defined in Eq. (2.33)) vanish for wavelengths larger than the Hubble radius before equality. This statement can be verified by using Eq.

Higher-order tight coupling expansion
After having taken the derivative with respect to the conformal time τ of the first relation appearing in Eq. (2.6), Eq. (2.9) can be used to eliminate the baryon-photon velocity; the result will of this manipulation is: where, according to Eq. (2.12), the Lorentz force has been traded for a combination of the magnetic energy density and of the normalized anisotropic stress. Equation (3.1) stems directly from the zeroth order in the tight-coupling expansion.
To discuss the polarization, we have to go (at least) to first-order in the tight coupling expansion. For this purpose, it is appropriate to introduce the evolution equations of the brightness perturbations of the I, Q and U Stokes parameters characterizing the radiation field. Since the Stokes parameters Q and U are not invariant under rotations about the axis of propagation the degree of polarization P = (Q 2 + U 2 ) 1/2 is customarily introduced [6,66]. The relevant brightness perturbations will then be denoted as ∆ I , ∆ P . This description, reproduces, to zeroth order in the tight coupling expansion, the fluid equations that have been presented in Section 2 to set initial conditions prior to equality. For instance, the photon density contrast and the divergence of the photon peculiar velocity are related, respectively, to the monopole and to the dipole of the brightness perturbation of the intensity field, i.e. δ γ = 4∆ I0 and θ γ = 3k∆ I1 . The evolution equations of the brightness perturbations can then be written, within the conventions set by Eqs. (2.3) and (2.5) as 3) include a dependence on P 2 (µ) = (3µ 2 − 1)/2 ( P ℓ (µ) denotes, in this framework, the ℓ-th Legendre polynomial); ; µ =k ·n is simply the projection of the Fourier wave-number on the direction of the photon momentum. In Eqs. (3.2) and (3.3) the source term S P is defined as (3.6) Equations (3.2) and (3.3) con now be expanded in powers of τ c = |1/κ ′ |. In particular we can write Equations (3.2) and (3.3) will be preliminarily phrased in terms of a Boltzmann hierarchy coupling together the different multipoles of the brightness perturbations. Then the obtained equations will be expanded according to (3.7). The result of the first step of this procedure is 14) Having derived Eqs. (3.8)-(3.11) and Eqs. (3.12)-(3.15), the system can now be expanded to the wanted order in τ c . As anticipated, to zeroth order in the tight coupling expansion Eq. (3.1) can be reobtained. In fact, to lowest order in τ c , which is the same of Eq. (3.1) if we recall that 4∆ I0 = δ γ . The presence of the magnetic field modifies the evolution of the zeroth order monopole and, consequently, according to Eq.  (3.19) which also implies Since τ c grows very fast during decoupling, to have a better quantitative estimate of the effect the procedure outlined in [6] can be followed. By taking the first derivative of the definition of S P and by using Eqs. (3.11) and (3.15) we can get, directly the following evolution equation for S P Equation (3.21) can be used to estimate the source term in Eq. (3.3) and, therefore, the polarization can be determined as where σ dec ≃ 70-80 Mpc is the width of the visibility function κ ′ e −κ which has a Gaussian form [67]. The quantity I 1 appearing in Eq. (3.22) is Let us now perform an oversimplified estimate. Let us ignore neutrinos and magnetic fields and let us compute out of Eq. (3.2) the ordinary Sachs-Wolfe contribution assuming a constant adiabatic mode that propagates smoothly through the matter-radiation transition. The constant mode will be ψ m ≃ −(3/5)R * where R * is the constant mode of curvature fluctuations. Then, from Eqs.
(3.24) which implies that the cross-correlation between temperature and polarization oscillates as sin (2kc sb τ dec ).
It is finally interesting to compute the expression of the characteristic damping scale of the fluctuations arising from the dispersion relations. From Eq. (3.4), by taking the Laplace transform, i.e. v b (k, τ ) = e i ω(τ ′ )dτ ′ v b (k, ω) (and similarly for the brightness perturbations), we will have that the following relation holds to second order in the tight-coupling expansion In the standard derivation of the dispersion relations the metric fluctuations are neglected. This procedure is also justified in the present case provided the spectrum of the magnetic energy density is not too steep in frequency. In practice this is the case for ε sufficiently smaller than 1. In this case also the magnetic contribution can be neglected. In this situation, we can Laplace transform also the other equations (i.e. for the monopole, dipole and quadrupole). Then the first-order tight coupling expressions can be used into Eqs. (3.9) and (3.25) to eliminate v b . The resulting expression will then be which implies that ω = ω 1 + iω 2 where , (3.27) In practice, for the typical set of parameters introduced around Eq. (2.10) k −1 d ≃ 3 Mpc, which sets, grossly speaking, the typical scale of diffusion damping.

Numerical analysis
The evolution equations in the tight coupling approximation will now be integrated numerically. The normalization of the numerical calculation is enforced by evaluating, analytically, the Sachs-Wolfe plateau and by deducing, for a given set of spectral indices of curvature and entropy perturbations, the amplitude of the power spectra at the pivot scale. Here is an example of this strategy. The Sachs-Wolfe (SW) plateau can be estimated analytically from the evolution equation of R (or ζ) by using the technique of the transfer matrix appropriately generalized to the case where, on top of the adiabatic and non-adiabatic contributions the magnetic fields are consistently taken into account. The result of this calculation can be expressed, for ℓ < 30, as [7] where the functions Z 1 , Z 2 and Z 3

.(4.4)
We recall the meaning of the various spectral indices appearing in Eqs. where τ 0 is the present observation time. In Eqs. (4.2), (4.3) and (4.4), n stands either for n r (i.e. the spectral index of the adiabatic mode defined in Eq. (2.35)) or for n s (i.e. the spectral index of the entropic mode defined in Eq. (2.54)) or for n rs (i.e. the spectral index of the cross-correlation). In principle one could also introduce a different spectral index for each cross-correlation. However, following, for instance [65] we will assume that the spectrum of the cross-correlation is solely determined in terms of the spectrum of the autocorrelations, i.e., for instance, n rs = (n r + n s )/2. In Eq. (4.1) γ rs , γ br and γ sb are the correlation angles. From the request that the SW plateau is dominated by the adiabatic mode it is possible to constrain Ω BL and, consequently, B L (see, indeed, Eq. (2.39)). The nature of the bound, however, depends on the parameters that enter Eq. (4.1). In the absence of magnetic and non-adiabatic contributions and for Eqs. (4.1) and Eq. (4.2) imply that for n r = 1 (Harrison-Zeldovich spectrum) ℓ(ℓ + 1)C ℓ /2π = A R /25 and WMAP data [12] would demand that A R = 2.65 × 10 −9 .
If the SW plateau is determined by an adiabatic component supplemented by a (subleading) non-adiabatic contribution both correlated with the magnetic field intensity the obtainable bound may not be so constraining (even well above the nG range) due to the proliferation of parameters. A possible strategy is therefore to fix the parameters of the adiabatic mode to the values determined by WMAP-3 and then explore the effect of a magnetized contribution which is not correlated with the adiabatic mode. This implies, in Eq. (4.1) that A S = 0 and γ br = π/2. Under this assumption, in Fig. 1 the bounds on B L are illustrated. The nature of the constraint depends, in this case, both on the amplitude of the protogalactic field (at the present epoch and smoothed over a typical comoving scale L = 2π/k L ) and upon its spectral slope, i.e. ε. In the case ε < 0.5 the magnetic energy spectrum is nearly scale-invariant. In this case, diffusivity effects are negligible (see, for instance, [14,15]). As already discussed, if ε ≫ 1 the diffusivity effects (both thermal and magnetic) dominate the mode-coupling integral that lead to the magnetic energy spectrum [14,15]. In Fig. 1 (plot at the left) the magnetic field intensity should be below the different curves if the adiabatic contribution dominates the SW plateau. Different choices of the pivot scale k p and of the smoothing scale k L , are also illustrated. In Fig. 1 (plot at the left) the scalar spectral index is fixed to n r = 0.951 [3]. In the plot at the right the two curves corresponding, respectively, to n r = 0.8 and n r = 1 are reported. If ε < 0.2 the bounds are comparatively less restrictive than in the case ε ≃ 0.9. The cause of this occurrence is that we are here just looking at the largest wavelengths of the problem. As it will become clear in a moment, intermediate scales will be more sensitive to the presence of fully inhomogeneous magnetic fields.
According to Fig. 1 for a given value of the magnetic spectral index and of the scalar spectral index the amplitude of the magnetic field has to be sufficiently small not to affect the dominant adiabatic nature of the SW plateau. Therefore Fig. 1 (as well as other similar plots) can be used to normalize the numerical calculations for the power spectra of the brightness perturbations, i.e. (4.5) Let us then assume, for consistency with the cases reported in Fig. 1, that we are dealing with the situation where the magnetic field is not correlated with the adiabatic mode. It is then possible to choose a definite value of the magnetic spectral index (for instance ǫ = 0.1) and a definite value of the adiabatic spectral index, i.e. n r (for instance n r = 0.951, in agreement with [3]). By using the SW plateau the normalization can be chosen in such a way the the adiabatic mode dominates over the magnetic contribution. In the mentioned case, Fig. 1  The power spectra of the brightness perturbations, i.e. Eq. (4.5), can be then computed by numerical integration. Clearly the calculation will depend upon the values of ω m , ω b , ω c and R ν . We will simply fix these parameters to their fiducial values reported in Eqs. (2.10) (see also (2.11)) and we will take N ν = 3 in Eq. (2.32) determining, in this way the fractional contribution of the neutrinos to the radiation plasma. The first interesting exercise, for the present purposes, is reported in Fig. 2 where, in the left column, the power spectra of the brightness perturbations are illustrated for a wave-number k = 0.1 Mpc −1 ; in the right column the power spectra of the same quantities are illustrated for k = 0.2 Mpc −1 . Concerning the results reported in Fig. 2 different comments are in order: • for ε = 0.1 and n r = 0.951, the SW plateau imposes B L < 1.14 × 10 −8 G; from Fig. 2 it follows that a magnetic field of only 30 nG (i.e. marginally incompatible with the SW bound) has a large effect on the brightness perturbations as it can be argued by comparing, in Fig. 2, the dashed curves (corresponding to 30 nG ) to the full curves which illustrate the case of vanishing magnetic fiels; • the situation where B L > nG cannot be simply summarized by saying that the amplitudes of the power spectra get larger since there is a combined effect which both increases the amplitudes and shifts slightly the phases of the oscillations; • from the qualitative point of view, it is still true that the intensity oscillates as a cosinus, the polarization as a sinus; • the phases of the corss-correlations are, comparatively, the most affected by the presence of the magnetic field.   The features arising in Figs. 1 and 2 can be easily illustrated for other values of ǫ and for different choices of the pivot or smoothing scales. The general lesson that can be drawn is that the constraint derived only by looking at the SW plateau are only a necessary condition on the strength of the magnetic field. They are, however, not sufficient to exclude observable effects at smaller scales. This aspect is illustrated in the plot at the left in Fig. 3 which captures a detail of the cross-correlation. The case when B L = 0 can be still distinguished from the case B L = 0.5 nG. Therefore, recalling that for the same choice of parameters the SW plateau implied that B L < 11.4 nG, it is apparent that the intermediate scales lead to more stringent conditions even for nearly scale-invariant spectra of magnetic energy density. For the range of parameters of Fig. 2 we will have that B L < 0.5 nG which is more stringent than the condition deduced from the SW plateau by, roughly, one order of magnitude.
If ε increases to higher values (but always with ε < 0.5) by keeping fixed B L (i.e. the strength of the magnetic field smoothed over a typical length scale L = 2π/k L ) the amplitude of the brightness perturbations gets larger in comparison with the case when the magnetic field is absent. This aspect is illustrated in the right plot of Fig. 2 where the logarithm (to base 10) of the intensity autocorrelation is evaluated at a fixed wave-number (and at τ dec ) as a function of ε. The full line (corresponding to a B L = 10 nG) is progressively divergent from the dashed line (corresponding to B L = 0) as ε increases.
In Fig. 4 the power spectra of the brightness perturbations are reported at τ dec and as a function of k. In the two plots at the top the autocorrelation of the intensity is reported for different values of B L (left plot) and for different values of ε at fixed B L (right plot). In the two plots at the bottom the polarization power spectra are reported always at τ dec and for different values of B L at fixed ε. The position of the first peak of the autocorrelation of the intensity is, approximately, k d ≃ 0.017 Mpc −1 . The position of the first peak of the cross-correlation is, approximately, 3/4 of k d . From this consideration, again, we can obtain that B L < 0.3 nG which is more constraining than the SW condition.
Up to now the adiabatic mode has been considered in detail. We could easily add, however, non-adiabatic modes that are be partially correlated with the adiabatic mode. It is rather plausible, in this situation, that by adding new parameters, also the allowed value of the magnetic field may increase. Similar results can be achieved by deviating from the assumption that the magnetic field and the curvature perturbations are uncorrelated. This aspect can be understood already from the analytical form of the SW plateau (4.1). If there is no correlation between the magnetized contribution and the adiabatic contribution, i.e. γ br = π/2, the SW plateau will be enhanced in comparison with the case when magnetic fields are absent. The same situation arises when the two components are anti-correlated (i.e. cos γ br < 0). However, if the fluctuations are positively correlated (i.e. cos γ br > 0) the cross-correlation adds negatively to the sum of the two autocorrelations of R and Ω B so that the total result may be an overall reduction of the power with respect to the case γ br = π/2.

Concluding remarks
In the recent past, various attempts have been made to constrain large-scale magnetic fields from CMB physics. While interesting results have been obtained for the vector and tensor modes of the geometry, the scalar modes have been, comparatively, less studied. In this study a systematic strategy to tackle this problem has been proposed and developed. The basic inspiration of the present approach can be summarized by stating that the effects of fully inhomogeneous magnetic fields on CMB anisotropies is not univocal. In the standard lore the only relevant parameters to define the effects of large-scale magnetic fields on CMB anisotropies are the ones related to the magnetic field itself. For instance, if the magnetic field is uniform the only parameter will be its intensity; if the magnetic field is fully inhomogeneous the only two parameters will be the spectral amplitude and the spectral slope. The point of view of the present investigation is that this way of thinking is ambiguous and, to some extent, not so productive. In fact, the spectral parameters of the magnetic field are certainly a necessary ingredient of the analysis. They are, however, not sufficient to conduct the calculation of the brightness perturbations. The information on the spectral parameters must be complemented by the specific pre-decoupling initial conditions. According to the results of this investigation, a sound procedure to be followed can be summarized as follows: • include the magnetic fields in the system of scalar modes before matter radiation equality; • solve the system in the long wave-length approximation prior to decoupling; • use the obtained solutions for the numerical evaluation of the brightness perturbations, for instance, in the tight coupling approximation.
The ambiguity mentioned above arises during the first step of the procedure we just outlined. In fact, the inclusion of magnetic fields prior to equality does not lead to a single solution but to different solutions. We will then have quasi-adiabatic and quasi-entropic modes where a magnetized contribution can be accommodated.
In some sense, the present proposal is more conservative than previous attempts: the signature of magnetic fields on CMB physics is not directly accessible but it is mediated by the (supposedly rich) parameter space of the pre-decoupling initial conditions of CMB anisotropies. This approach is indeed common when mixtures of adiabatic and non-adiabatic modes are studied in the absence of magnetic fields. In the latter case the philosophy is to specify what kind of non-adiabatic mode we want to study and constrain. Here the situation is similar and we do hope, that, in the future, the considerations reported here will serve as an inspiration for analyzing the effects of fully inhomogeneous magnetic fields on the scalar CMB anisotropies.
This analysis reported in this paper allowed to extend the tight-coupling approximation to the case when large-scale magnetic fields participate to the dynamics of the plasma. The evolution equations, in the tight coupling approximation, have then been integrated numerically using, as normalization, the analytical estimate of the Sachs-Wolfe plateau modified by the presence of large-scale magnetic fields. The reported results allow not only to set constraints but also to include a magnetized contribution in the current strategies of parameter extraction.
From the combined analysis of the Sachs-Wolfe plateau and of the phases of Sakharov oscillations at intermediate scales, a fully inhomogeneous magnetic field uncorrelated with the adiabatic mode can be constrained to be smaller than 0.5 nG for a slightly blue spectral slope ε = 0.1 and for typical comoving wave-number of 1 Mpc −1 . The conditions implied by the SW plateau ssem to be necessary but not sufficient since magnetic fields compatible with the SW bound may alter Sakharov oscillations at smaller length-scales. In the most general situation, correlation angles may proliferate: there may be correlations between the magnetized contribution and the adiabatic or even non-adiabatic modes. In this case the bounds on large-scale magnetic fields become much weaker and sizable magnetic fields may be allowed.
The present analysis also shows that magnetic fields can affect the process of formation of CMB polarization. For instance, the usual approach to the Faraday rotation of CMB polarization is to assume that the polarization itself is only due to the adiabatic mode while the rotation of the polarization plane is only due to the magnetic field. This picture has been complemented by the observation that the magnetic fields affect the pre-decoupling initial conditions and, consequently, modify the evolution of the degree of polarization.
The tight coupling expansion is rather useful to discuss the integration of brightness perturbations in some simple cases. Moreover, it is an essential ingredient of Boltzmann codes since it is used, at early times, to avoid the integration of the (stiff) Euler equations. It will be interesting, for future studies, to analyze the space of magnetized initial conditions by solving directly the whole Boltzmann hierarchy.