Tensor B mode and stochastic Faraday mixing

This paper investigates the Faraday effect as a different source of B mode polarization. The E mode polarization is Faraday rotated provided a stochastic large-scale magnetic field is present prior to photon decoupling. In the first part of the paper we discuss the case where the tensor modes of the geometry are absent and we argue that the B mode recently detected by the Bicep2 collaboration cannot be explained by a large-scale magnetic field rotating, through the Faraday effect, the well established E mode polarization. In this case, the observed temperature autocorrelations would be excessively distorted by the magnetic field. In the second part of the paper the formation of Faraday rotation is treated as a stationary, random and Markovian process with the aim of generalizing a set of scaling laws originally derived in the absence of the tensor modes of the geometry. We show that the scalar, vector and tensor modes of the brightness perturbations can all be Faraday rotated even if the vector and tensor parts of the effect have been neglected, so far, by focussing the attention on the scalar aspects of the problem. The mixing between the power spectra of the E mode and B mode polarizations involves a unitary transformation depending nonlinearly on the Faraday rotation rate. The present approach is suitable for a general scrutiny of the polarization observables and of their frequency dependence.


Introduction
The Bicep2 collaboration has recently reported the detection of a B mode polarization that has been attributed to the presence of gravitational waves of inflationary origin [1,2] (see also [3] for the detection of the B mode coming from lensing). The tensor fluctuations of the geometry are able to generate a B mode polarization in the Cosmic Microwave Background (CMB in what follows) provided the typical wavelengths of the relic gravitational waves were larger than the Hubble radius after matter radiation equality but before decoupling, i.e. at the moment when the initial conditions of the polarization anisotropies are set.
In principe the B mode polarization maybe the result of a more mundane process well known in the treatment of cold plasmas, namely the Faraday effect. According to the Faraday effect the polarization plane of the incoming radiation is rotated because of the presence of a magnetic field in a medium with finite density of charge carriers. The latter requirement is met by the predecoupling plasma that is globally neutral but intrinsically charged since the electrons and the ions have a common concentration that is O(10 −10 ) the concentration of the photons at the corresponding epoch. The only further assumption to get a B mode polarization is therefore the presence of a magnetic field that will be assumed to be stochastic not to conflict with the assumed isotropy of the background geometry.
The Faraday effect of the CMB polarization was analyzed almost two decades ago and, to some extent, even before (see [4] and references therein). These suggestions and have been subsequently discussed in a number of different articles (see [5] and references therein). There are, in principle, other sources of B mode polarization due to the magnetic fields but they are smaller than the contribution of the Faraday effect.
The first task of the present paper will therefore be to establish if the observed B mode polarization can be ascribed to the Faraday effect. The answer to the question will be, in short, that the observed B mode polarization cannot be attributed predominantly to a Faraday rotated E mode polarization since the magnetized temperature autocorrelations would be too distorted.
In the second and more technical part of the paper the attention will then be focussed on the interference of the Faraday effect with the B mode produced by the tensor mode of the geometry. The aim will be to derive a set of scaling rules that could be directly applied to the E mode and B mode power spectra. The Faraday rotation will be described terms of a stationary, quasi-Markovian and random process [17]. It will be shown that the evolution of the brightness perturbations obeys a set of stochastic differential equations that can be solved using the cumulant expansion [18,19], pioneered in similar contexts by Van Kampen.
The stochastic approach to the Faraday effect has been exploited in astrophysics where the source of linear polarization is provided by the properties synchrotron emission [20,21,22]. It has been recently suggested [23] that a consistent stochastic description can be successfully achieved in the case of the Cosmic Microwave Background (CMB in what follows) where the linear polarization is primarily provided by the adiabatic initial conditions of the Einstein-Boltzmann hierarchy [24].
The B mode induced by the stochastic Faraday effect thanks to the presence of the linear polarization of the CMB can be expressed, according to the results of [23] as: where C (EE) ℓ denotes the autocorrelation of the E mode polarization obtained in the absence of stochastic Faraday term and ω F is given by: ( 1.2) In Eq. (1.2) X F (τ ) denotes the Faraday rotation rate and the stochastic process has been assumed, for sake of simplicity, homogeneous in space. The derivation of Eq. (1.1) does not demand ω F ≪ 1 and shows that the stochastic Faraday rate affects not only the B mode polarization but, to some extent, also the E mode itself. Equations (1.1) and (1.2) have been derived in [23] by assuming that the sole sources of linear polarization were the scalar fluctuations of the geometry. It was anticipated in [23] that the results could be extended to the case where the initial source of polarization is not only provided by the scalar modes but also by the tensor modes that appear in one of the minimal extensions of the so-called ΛCDM paradigm, where Λ stands for the dark energy component and CDM for the cold dark matter contribution. While this analysis was in progress there have been claims of detection of a primordial B mode polarization by the Bicep2 collaboration [1] (see also [2]) complementing the results of the B mode from lensing [3]. Although these data are in tension with other data sets for various reasons, it seems timely to present a derivation of the analog of Eq. (1.1) when the sources of polarization is not only provided by the standard adiabatic mode but also by the tensor fluctuations of the geometry. The main result of the present analysis can be summarized by writing the analog of Eq. (1.1) in this extended set-up where the E modes and the B mode power spectra of the tensors are included: (1. 3) as in Eq. (1.1) C (EE) ℓ denotes the E mode power spectrum coming from the scalar modes of the geometry while C (BB) ℓ and C (EE) ℓ (both in calligraphic style) denote, respectively, the polarization observables induced by the tensor modes of the geometry. Following the standard terminology, the B-mode autocorrelations are denoted by BB. With similar logic, we talk about the TT, TE and EE angular power spectra denoting, respectively, the autocorrelations of the temperature, the autocorrelations of the E mode and their mutual cross correlations. It is appropriate to recall that the tensor modes of the geometry not only produce BB correlations but also EE and TT power spectra (see e.g. [25]). Comparing Eqs. (1.2) and (1.3) in the limit ω F → 0 we can appreciate that the B mode polarization disappears from Eq. (1.2) while it persists in Eq. (1.3) and it is solely given by the tensor B mode.
According to Eqs. (1.2) and (1.3) both the E mode and the B mode polarization are frequency dependent since ω F is proportional to the square of the rate and, ultimately, to the fourth power of the comoving wavelength. The stochastic approach to the Faraday rate represents an ideal framework for deriving a set of scaling laws only involving the measured polarization power spectra. The present findings support the consistency of the whole description and are suitable for a discussion of the Faraday effect when the predominant source of the B mode polarization is provided by relic gravitons with wavelengths comparable with the current Hubble radius.
The present paper is organized as follows. In section 2 we shall examine, in the light of the Bicep2 findings, the Faraday interpretation of the B mode polarization. In section 3 we shall corroborate the analysis with a numerical discussion. In section 4 we shall describe the phenomenon of stochastic Faraday mixing. The polarization observables and their scaling properties will be deduced in section 5 while the concluding remarks will be collected in section 6. To avoid digressions some relevant technical aspects of the discussions have been collected in the appendices A, B and C.

Bicepobservations and Faraday effect
The normalization provided by the Bicep2 results [1] should be satisfied by any plausible physical explanation of the observed B mode autocorrelation. The BB power spectrum should be at most O(10 −2 ) µK 2 for typical angular scales of the degree. More specifically 2 we can estimate that where ℓ ≃ 248 denotes the central value of a bunch of multipoles ranging from 231 to 265. For comparison with the observational data, we shall refer, when needed, to the order of magnitude suggested by Eq. (2.1). Indeed, for the lowest bunch of multipoles (centered around ℓ ≃ 45) G Bℓ ≃ 8.6 × 10 −3 µK 2 . The remaining eight bunches of multipoles give values between 1.3 × 10 −2 µK 2 (for ℓ ≃ 73) and 3.2 × 10 −2 µK 2 (for ℓ ≃ 317). The observational frequency of Bicep2 is O(150) GHz and it is larger or even much larger than the frequencies of some of the previous polarization experiments (see [6,7,8,9]). The maximalist perspective stipulates that the observed BB correlation is just the result of the Faraday rotated EE spectrum whose origin stems from the adiabatic (scalar) fluctuations of the geometry. Conversely, in the minimalist perspective the observed BB spectrum 2 The normalization of the B mode autocorrelation is often reduced to equivalent tensor to scalar ratio r T . We prefer to quote here the B mode autocorrelation in its physical units (i.e. µK 2 ) since we are going to investigate here a complementary perspective of the problem not relying necessarily on the tensor modes of the geometry. is originated by a primordial tensor mode with tensor to scalar ratio r T ≃ 0.2: according to this school of thought magnetic fields as well as plasma effects are absent. There is finally a third logical possibility namely that the BB and the EE correlations are mixed by the Faraday effect, as suggested by Eq. (1.3) and by the viewpoints conveyed in this investigation. The minimalist perspective cannot be excluded, a priori, but it is fair to say that the Bicep2 data are in a certain tension with other CMB data set. We shall therefore focus hereunder on the remaining two options giving specific criteria for their empirical exclusion.
In the maximalist perspective the tensor E and B modes disappear from Eq. (1.3)). Since the tensors are totally absent we could use as well Eq. (1.1). The observed B mode depends then on ω F whose explicit value can be roughly estimated as: showing that the actual value of X F /ǫ ′ is not necessarily much smaller than 1 and it is O(1) for comoving field strengths of a few nG (i.e. 1 nG = 10 −9 G) and frequencies O (10) GHz [23]. Assuming now that G Eℓ is well estimated by the measured E mode autocorrelation [24,10] we shall have The actual value of the maximum of G Eℓ is slightly smaller than 50 µK 2 (i.e. O(43) µK 2 [24,10] The value 50 µK 2 used in Eq. (2.4) maximizes G Eℓ but it corresponds to relatively high multipole moments, typically ℓ ∼ 10 3 ; for smaller multipoles (compatible with the Bicep2 observations) we have that G Eℓ < 5 µK 2 . Moreover the value of 1 nG for the magnetic field intensity is barely compatible with the distortions produced by a large-scale magnetic field on the temperature autocorrelation: we must bear in mind that the scalar inhomogeneities induced by an inhomogeneous magnetic field are the leading source of distortion of the TT, EE, and TE angular power spectra in comparison with vector and tensor modes. The magnetized CMB observables have been derived for the magnetized adiabatic mode and compared with the available experimental data with the aim of pinning down the properties of the magnetic field. In the second and third papers of Ref. [11] this analysis has been performed, for the first time, using the WMAP five-year data and later confirmed by subsequent analyses and different data sets (see, e.g. [12]).
The estimates of Eqs. (2.2), (2.3) and (2.4) can be further refined. Assuming that the Faraday rate is perturbative, the angular power spectrum of the Faraday rotated E mode can be computed. This approximation boils down to a sharp separation between the moment of formation of the polarization from the moment of the Faraday rotation of the produced polarization (see e.g. last two papers of Ref. [5]). The result of the discussion is conceptually similar to the one of Eq. (2.4) but mathematically more accurate as far as the scaling with the multipoles is concerned. Defining the normalized form of the Faraday rotation rate X F (n), we have that is the angular power spectrum of the normalized rate. In terms of the power spectrum of the Faraday rate the autocorrelation of the B-mode polarization can be computed as C where Z(ℓ, ℓ 1 , ℓ 2 ) contains also a Clebsch-Gordon coefficient 3 while C EE ℓ 2 is the E-mode power spectrum already discussed above; the sum of Eq. (2.6) must be conducted in compliance with the constraints stemming from the triangle inequality |ℓ 1 − ℓ 2 | ≤ ℓ ≤ ℓ 1 + ℓ 2 . The Faraday effect can be treated, within this approach, either with uniform magnetic fields or with stochastic magnetic fields and different analyses have been performed starting with the ones of Ref. [4] (see also [5] for an incomplete list of references).
Using analytic methods it is possible to estimate Eq. (2.6) at small angular scales (i.e. ℓ 1 ≫ 1, ℓ 2 ≫ 1 and ℓ ≫ 1) where the Clebsch-Gordon coefficient inside Z(ℓ, ℓ 1 , ℓ 2 ) can be evaluated in analogy with the semiclassical limit in non relativistic quantum mechanics. This approach to the asymptotics of the Clebsch-Gordon coefficients was originally studied by Ponzano and Regge [15] by exploiting the connection of the Clebsch-Gordon coefficients with the Wigner 3j and 6j symbols (see also [16]). This analytical technique has been exploited in the last paper of Ref. [5] for the explicit estimates of Z(ℓ, ℓ 1 , ℓ 2 ). The result is that the magnetic field must be larger than about 15 nG for different ranges of the spectral indices if we ought to be compatible with the orders of magnitude of Eq. (2.1). Thus, the semianalytical considerations suggest that to have a chance of observing the Bicep2 value of Eq. (2.1) we are led to magnetic field intensities larger than O(10) nG. Let us now corroborate the previous considerations with an explicit numeric evaluation of Eq. (2.6).

TT correlations and Faraday B mode
The two-point function of scalar curvature perturbations in Fourier space 1) shall be normalized at the pivot scale k p = 0.002 Mpc −1 ; using the WMAP 9yr data alone [24] in the light of the concordance scenario we have that A R = (2.41 ± 0.10) × 10 −9 and n s = 0.972 ± 0.013. The exact scale-invariant limit is realized when n s → 1.
The pivotal parameters of the ΛCDM paradigm can be determined on the basis of different data sets 4 and, for illustrative purposes, we shall consider only three examples. The first one is obtained by comparing the ΛCDM paradigm to the WMAP9 data alone (see, in particular, [24]): (Ω b0 , Ω c0 , Ω de0 , h 0 , n s , ǫ re ) ≡ (0.0463, 0.233, 0.721, 0.700, 0.972, 0.089), (3.2) with A R = 2.41 × 10 −9 . If we include the data sets pertaining to the baryon acoustic oscillations (see, e.g. [13]) the string of parameters of Eq. ( . Many other sets of parameters corresponding to the combinations of different data sets can be used. These differences are totally immaterial for the present considerations and we shall therefore adopt, for illustration, the fiducial set of parameters of Eq. (3.2). The minimal scenario where the magnetic fields can be consistently included is sometimes dubbed the magnetized ΛCDM scenario: the only two supplementary parameters in comparison with the ΛCDM parameters are the magnetic field amplitude and the magnetic spectral index (see below Eqs. As it can be easily verified from the definition of the Fourier transform, P B (k, τ ) has dimensions of an energy density and its square root has, therefore, the dimensions of a field intensity 5 . In Eq. (3.3), A B has the correct dimensions of an energy density and can be related to the regularized magnetic field intensity B L which is customarily employed to phrase the comoving values of the magnetic field intensity. In the case when n B > 1 (i.e. blue magnetic field spectra), In the case of white spectra (i.e. n B = 1) the two-point function is logarithmically divergent in real space and this is fully analogous to what happens in Eq. (3.1) when n s = 1, i.e. the Harrison-Zeldovich (scale-invariant) spectrum. Quasi-scale invariant spectra with red tilt (i.e. n B < 1) can arise when magnetic fields are produced in the context of conventional inflationary models (see the last paper of Ref. [12]) but for numerical purposes we shall mainly focus the attention on blue spectra.
In Fig. 1 with the full, dashed and dot-dashed lines we report the results for the BB spectrum produced by the Faraday effect and computed on the basis of Eq. (2.6) after having included the magnetic fields in the Einstein-Boltzmann hierarchy as in [11]. In Fig. 1 the values of the magnetic fields range between 15 and 20 nG while the magnetic spectral index has been fixed to n B = 1.5. As indicated, the other parameters of Fig. 1 have been fixed to the best fit values given in Eq. (3.2). Both plots in Fig. 1 share the same parameters but the plot on the right is focussed on the large angular scales while the plot on the left illustrates the small angular scales. Semilogarithmic scales are used in both plots. Figure 1 suggests that magnetic fields O(10) nG are unable to reproduce the observed Bicep2 amplitude. The analytical estimates of Eqs. (2.3)-(2.4) and (2.6) are quantitatively correct but excessively optimistic. As anticipated, the EE correlations have been purposely overestimated in Eqs. (2.3) and (2.4).
In Fig. 2 the magnetic spectral index has been fixed at n B = 1.5 (plots at the top) and at n B = 2 (plots at the bottom). The full, dashed and dot-dashed curves in the various plots of Fig. 2 denote, respectively magnetic field intensities of 1, 5 and 10 nG. In a frequentistic perspective, beyond the outer contours of the exclusion plots of the second paper Ref. [11] the parameters of the magnetized ΛCDM scenario are excluded to 95% confidence level (see in particular Figs This means that if we compute the TT correlations for the models of Figs. 1 or 2 we shall see that they are excessively distorted. According to Fig. 1 the model that is closer to the BICEP2 data is the one illustrated by the dod-dashed line. The TT correlations for such a model are compared, in Fig. 3 with the best fit to the WMAP9 data alone. The conclusion is that the models of Fig. 1 get close to the Bicep2 measurement but are already excluded by the analysis of the other CMB observables. Magnetic fields systematically less intense than the ones of Figs. 1 and 2 lead to B mode polarizations that are minute in the Bicep2 region. The value of B L ≃ O(nG) leads to a B mode signal O(10 −6 ) µK 2 . Let us finally check if and how we could find some sort of model that could be barely compatible with the BICEP2 data at large scales while still being compatible with the previous bounds on the B mode polarization at smaller scales. This aspect is discussed in Fig. 4 where the Bicep2 data points seem to be vaguely close to the curve but at smaller angular scales scales (i.e. ℓ ∼ 10 3 ) the B mode polarization is O(0.8) (µK) 2 which is actually an enormous value (potentially in conflict with existing upper limits at small angular scales) and leading to ridiculously large distortions of the TT correlations of the type illustrated in Fig. 2.

Brightness perturbations
The main objective of the the present section and of the following one will be to derive the results contained in Eq. (1.3). Further details on the cumulant expansion and on the technicalities involved in the calculation are collected in Appendix C. At the end of section 5 this technical effort will lead to the derivation of a set of scaling laws that can be directly applied to the angular power spectra.

General considerations
We shall consistently work in a conformally flat space-time whose line element and metric tensor are defined as: a(τ ) is the scale factor of a conformally flat geometry of Friedmann-Robertson-Walker type and τ is the conformal time coordinate. The fluctuations of the Stokes parameters in comparison to their equilibrium values can be decomposed as: where X = I, Q, U, V denotes one of the four Stokes parameters and where the superscripts refer, respectively, to the scalar, vector and tensor modes of the geometry. The brightness perturbations 6 defined in Eq. (4.2) are affected by the presence of the In Eq. (4.4) X F ( x, τ ) denotes the Faraday rotation rate: where ω = 2πν is the (comoving) angular frequency, while ω Be and ω pe denote the comoving Larmor and plasma frequencies, respectively; n e =ñ 0 a 3 is the comoving electron concentration.
If the time scale of spatial variation of the rate is comparable with the time scale of spatial variation of the gravitational fluctuations, X F can be considered only time dependent (i.e. a stochastic process). In the opposite situation the Faraday rate must be considered fully inhomogeneous (i.e. a stochastic field). Both possibilities will be considered. There is, in principle, also a third case, namely the situation where X F is just a (time-independent) random variable characterized by a given probability distribution. This is, in some sense, the simplest and most naive case and it follows from the present results when X F is considered a space-time constant with random distribution. It is finally useful to recall that the stochastic Faraday effect has been discussed when X F is a constant random variable in the framework of the synchrotron emission (see, e.g. the second paper of Ref. [21]).
Equations (4.3), (4.4) and (4.5) may also depend on the comoving frequency for a different reason. Since the magnetic field modifies the trajectories of the charged particles it also affects the collision matrix for the photons impinging on the electrons. The collision matrix will the inherit corrections O(f 2 e ) where f e (ω) is defined as: and z * is the redshift to last scattering. This effect is well understood but rather difficult to compute when scalar, vector and tensor modes of the geometry are included. Since the typical frequencies of the observational channel are normally much larger than the Larmor frequency of the electrons, the corrections O(f 2 e ) are negligible. This comment anticipates a possible objection and demonstrates that, in this context, the role of the magnetic field in the collision matrix can be ignored at least in the first approximation. This analysis has been however performed in a series of papers in connection with the problem of the circular polarizations (see [26] and references therein). To be precise the exact form of the correction O(f 2 e ) to the full equations will be reported below in the scalar and tensor cases.

Scalar brightness perturbations
The scalar equations can be written, in a compact form, as: The evolution equations of the polarization do not contain explicitly the fluctuations of the metric that appear instead in the evolution of ∆ where the corrections O(f 2 e ) have been only indicated but will be explicitly included in the full equations (see below); the term S (s) P denotes the standard combination: 14) The conventions adopted for the definition of the multipoles of the polarization and of the intensity include the factor (2ℓ + 1) in the expansion, i.e. δ (s) where P ℓ (µ) denote the Legendre polynomials. When X F (τ ) is a homogenous stochastic process the equation for the intensity is 7 ∂ τ δ (s) , (4.17) 7 We have also introduced for completeness and in analogy with the forthcoming tensor case the contribution coming from the scalar metric fluctuations and modifying the collisionless part of the Boltzmann equation. The corresponding details and the related conventions are collected in appendix B.
while the equations for the polarization are: (4.18) In Eq. (4.18) we included also the correction (parametrized by f 2 e ) arising when we take into account the magnetic field in the scattering term.
If X F (τ ) is a spatially homogenous stochastic process the relevant equation we shall be dealing with is: If X F is a spatially inhomogeneous stochastic process the discussion in mathematically slightly different but physically equivalent as far as the frequency scaling is concerned. More specifically, the evolution equations for δ ± will now contain a convolution and can be written as: (4.20) where we introduced, for convenience, the term b F (ν, τ ) = 2 e 3 n e /[(2π) 5/2 m 2 e a 2 (τ )ν 2 ] and neglected, for simplicity, the corrections O(f 2 e ). The addition of the spatial dependence is just a technical complication since the essential aspect is the stochastic evolution in time.

Tensor brightness perturbations
The tensor brightness perturbations can still be written in compact notation as: (4.22) In the scalar case the evolution of the intensity is coupled to the polalrization only through the source term S (s) P . In the tensor case something similar happens with the difference that the two polarizations of the relic tensor wave will determine also a specific azimuthal dependence of the other brightness perturbations. Recalling the results of the appendices A and B we can write, in Fourier space, where we have considered the case where X F (τ ) is a stochastic process.
Before computing the collision terms it is therefore necessary to specify the azimuthal structure of the relevant brightness perturbations. Since (n ·â) 2 − (n ·b) 2 = (1 − µ 2 ) cos 2ϕ and 2(n ·â)(n ·b) = (1 − µ 2 ) sin 2ϕ the intensity brightness must be written as (4.24) The functions Z ⊕ (k, τ ) and Z ⊗ (k, τ ) are obey the same equation since the components of the tensor polarization h ⊕ and h ⊗ obey the same equation. Keeping track of all the normalizations we can therefore write Eq. (4.24) as: 27) where D T (k,n) and D ± (n) are defined as: In Eq. (4.28) two time-independent stochastic variables F 1 ( k) = (F ⊕ − iF ⊗ )/ √ 2 and F 2 (k)(F ⊕ − iF ⊗ )/ √ 2 have been introduced. Each of the two stochastic polarizations F ⊕ and F ⊗ is related to the tensor power spectrum P T (k) as 29) but F ⊗ ( k)F ⊕ ( p) = 0. The source terms are, in the tensor case: (4.31) where S (t) P (k, τ ) is given by: and has been computed to lowest order in X F , i.e. δ (t) P (k, τ ). The evolution equations become, therefore, As in the scalar case, (see Eq. (4.20)), Eq. (4.34) can be generalized to the case where the stochastic process is not spatially homogeneous leading, in Fourier space, to an integrodifferential equation containing also a convolution.

E mode and B mode power spectra
In the present discussion there are in fact two E mode polarizations and two B mode polarizations. All the four angular power spectra are affected by the stochastic Faraday rate.
The two E modes come from the adiabatic scalar mode and from the tensor fluctuations of the geometry: the scalar E mode induces a B mode but also the tensor B mode is modified by the stochastic Faraday effect. The total polarization power spectra are related to their scalar and tensor components by means of a unitary transformation. To clarify the discussion as much as possible we shall carefully distinguish between the scalar and the tensor contributions to the polarization power spectra.

Scalar case
The brightness perturbations ∆ where N ℓ = (ℓ − 2)!/(ℓ + 2)!. In real space the scalars ∆ (s) E (n, τ ) and ∆ (s) B (n, τ ) can be expressed in terms of the generalized ladder operators [30] raising and lowering the spinweight of a given function: (5.5) The differential operators appearing in Eqs. (5.4) and (5.5) are generalized ladder operators (see [30], first paper) whose action either raises or lowers the spin weight of a given fluctuation. They are defined within the present conventions as acting on a fluctuation of spin weight j: For instance K ( The lengthy expressions reported in Eqs. (5.8) and (5.9) simplify greatly since the scalar modes do not have azimuthal dependence: The angular power spectra for the E mode polarization and for the B mode polarization are then defined as where ... denotes the ensemble average. We can finally determine a (E) ℓm and a (B) ℓm within the set of conventions followed here: The solutions for δ (s) ± (see Eq. (C.5)) must then be used inside Eq. (5.13). The the obtained expression can be averaged over the stochastic process by using Eqs. (C.8)-(C.11). Thus, if the scalars would be the only source of polarization the angular power spectra of the polarization can then be expressed as where C

(EE) ℓ
is the E-mode autocorrelation produced by the standard adiabatic mode and in the absence of Faraday mixing: In Eq. (5.16) j ℓ (x) are spherical Bessel functions and x = k(τ − τ 1 ). As already mentioned, the results have been derived in the sudden decoupling limit but can be extended to the case where the visibility function has a finite thickness [27,28,29]. Note, finally, that sometimes it is common to separate S (s) is the power spectrum of the constant adiabatic mode defined as in section 2 (see Eq. (3.1)). The result of Eq. (5.14) has been already anticipated in [23] and will now be completed by computing the tensor contribution.

Tensor case
Also the tensor brightness perturbations ∆ (t) ± (n, τ ) can be expanded in terms of spin-±2 spherical harmonics ±2 Y ℓ m (n), i.e. (5.17) In full analogy with the previous case the E and B modes are defined as ℓ m have the property of being invariant under rotations on a plane orthogonal ton and they can be be expanded in terms of (ordinary) spherical harmonics: ℓ m Y ℓ m (n), (5.19) where N ℓ = (ℓ − 2)!/(ℓ + 2)!. From Eqs. (5.4) and (5.5) written in the tensor case we can ℓm and, in particular, we shall have B (k,n, τ ). (5.20) After some algebra ∆ (t) E (τ, k,n) and ∆ (t) B (τ, k,n) can be expressed as: (5.22) where x = k(τ − τ 1 ); the functions Y ± (k, x, τ, τ 1 ) and q ± (n) are defined as Furthermore Λ E (x) and Λ B (x) are two differential operators The EE and BB angular power spectra are then defined as (5.26) As in the scalar case, using Eqs. (C.8)-(C.11) and the results of section Appendix C, the stochastic averages can be computed and the angular power spectra of the polarization can then be expressed as 27) where C (EE) ℓ and C (B) ℓ (in calligraphic style) are the E-mode autocorrelation produced by the standard adiabatic mode and in the absence of Faraday mixing: 5.30) As in the scalar case the results have been derived in the sudden decoupling limit [27,28,29]. Furthermore, as in the scalar case it is possible to separate the polarization source as S (t) P (k, τ 1 ) = P T (k)S (t) P (k, τ 1 ) where P T (k) is the tensor power spectrum defined as in Eq. (4.29) and related to the scalar power spectrum as P T (k) = r T P R where r T is the tensor to scalar ratio measuring the ratio between the scalar and tensor contributions at the pivot scale k p (see also Eq. (3.1)). Equations (5.27) hold under the academic hypothesis that the tensors are be the only source of polarization. In the realistic case Eqs. (5.14) and (5.27) shall be modified even further.

Total polarization power spectra
If the scalar and tensor modes of the geometry are simultaneously present the total angular power spectra are obtained from Eqs. (5.14) and (5.27) and the result is: ( 5.31) As already mentioned in the introduction C (EE) ℓ and C (BB) ℓ denote, respectively, the total E mode and B mode power spectra, while C (EE) ℓ is the E mode coming from the adiabatic scalar mode. The calligraphic power spectra C (EE) ℓ and C (BB) ℓ denote, respectively, the angular power spectra coming from the tensor modes. The explicit expressions of ω F can be rather different but the frequency dependence will always be the same: since ω F is quadratic in the rates it will always scale as 1/ν 4 ≃ λ 4 where λ denotes the wavelength of the channel. Since the scale factor is normalized in such a way that a 0 = 1, physical and comoving frequencies coincide today but not in the past. Equation (5.31) does not assume that the rotation rate is perturbative. In the limit ω F ≪ 1 Eq. (5.31) becomes: Equations (5.31) and (5.32) imply that the E mode coming from the scalar adiabatic fluctuations can be turned into a B mode. However, if a tensor B mode is present, the total angular power spectrum of the E mode is also affected. Finally, in the limit ω F → 0 there is no stochastic mixing and the standard situation is recovered: both the scalar and tensor contributions to the E mode and only the tensor contribution to the total B mode power spectrum.

Scaling relations and sum rules
In Ref. [23] some relations have been derived in the absence of the tensor B mode polarization and it was argued that these relations could be modified if the tensor B mode was included from the very beginning. We are now in a position to show that this conclusion is incorrect since the relations obtained in [23] are preserved even in the presence of the tensor B mode. The reason for this simple result can be understood, a posteriori, from the unitary transformation connecting the total angular power spectra to their scalar and tensor components.
Let us therefore show that the nonlinear combinations proposed in [23], in the absence of the tensor B mode, hold also in the present case. Defining the properly normalized total power spectra of the E mode and B mode polarizations: the first nonlinear combination guessed in Ref. [23] was: The result expressed by Eq. (5.34) holds, indeed, also when G Eℓ and G Bℓ are constructed in terms of the spectra of Eq. (5.31). To demonstrate this point it suffices to insert Eq. (5.31) into Eq. (5.33) and into the definition of L 0 (ω F ) given in the first equality of Eq. (5.34). One can easily think of two further combinations with well defined scaling properties with ω F , namely: As in the case of Eq. (5.34), also Eqs. (5.35) and (5.36) can be directly verified by direct use of Eqs. (5.31) and (5.33). Different nonlinear combinations can be invented either by combining nonlinearly L 0 (ω F ), L 1 (ω F ) and L 2 (ω F ) or by concocting new variables.

Mixing of the B modes
It would be nice to have definite criteria in order to rule out or to rule in the explanation of the detected B mode in terms of Faraday rotation. In this respect there are two possible criteria. One is to look for the frequency scaling and the other is to look for the scaling with the multipoles. Both properties are well understood analytically and numerically. In this perspective we find that the existing bounds on the B mode polarization at low [6,7], intermediate [14] and high frequencies [9,10] should probably be revisited. Here we just list some simple considerations. The idea, in short, goes as follows. Let us suppose that one or more experiments measures both the E mode and the B mode polarizations in different frequency channels. If there is mixing between the B mode polarization of tensor origin and the B mode induced directly by Faraday rotation, then, according to Eq. (5.31) various combinations with definite scaling properties with the comoving frequency can be constructed and some examples have been listed in Eqs. (5.34), (5.35) and (5.36). If the E mode and the B mode autocorrelations are independently measured in each frequency channel of a given experiment, both scaleinvariant and scale-dependent combinations of the angular power spectra can be constructed frequency by frequency. So, for instance the combination is asymptotically frequency-independent while the combination is truly frequency-independent. Equations (5.34), (5.35) and (5.36) (or any other nonlinear combination of the power spectra constructed with the same criteria) represent a set of physical observables that can be used to discriminate between the frequency dependence induced by the stochastic Faraday mixing or by other concurrent forms of frequency scaling caused either by the known or by the yet unknown foregrounds.
In this respect we ought to conclude this section with an extremely relevant but simple comment. If the BB correlation comes entirely from the tensor modes of the geometry, the internal linear combination (ilc) technique can be applied indifferently for all the channels of the experiment that eventually detect the B mode. In practice, the ilc map is a weighted linear combination over the smoothed maps obtained from each of the different frequency channels. Conversely, if the signal contains a frequency-dependent part, as the Faraday mixing would predict, the ilc technique cannot be applied. In these circumstances, the scaling relations and the sum rules obtained here are crucial if we intend to disentangle the real physical effects from potential foregrounds.

Concluding remarks
This paper investigated the Faraday effect of the CMB as a different and more mundane source of the B mode polarization detected by Bicep2. In the first part of the paper we discussed a maximalist alternative to the tensor B mode where the whole Bicep2 data are explained by a Faraday rotated E mode polarization. In the second part of the paper we discussed the possibility where the tensor B mode interferes with the Faraday rotated E mode polarization.
It has been shown both analytically and numerically that the Faraday rotation alone cannot explain the Bicep2 data. If this happens other CMB observables will be excessively distorted. The first estimate can be obtained by maximizing the E mode autocorrelation and by computing the induced B mode polarization. In this case we see that, given the Bicep2 frequency (i.e. 150 GHz), the Bicep2 normalization can only be reproduced when the magnetic field is O (15) nG. This value of the magnetic field is too large since it would induce unobserved distortions in the temperature autocorrelations. Indeed much lower magnetic fields (i.e. O(1.5) nG ) already produce excessive distortions on the TT correlations if the magnetic power spectrum is nearly scale-invariant (i.e. n B → 1). An independent test, in this respect, is provided by the frequency scaling of the signal which can be separately discussed. Other signals of B mode polarization can be induced directly by the tensor and vector modes of the geometry induced by the magnetic fields. It is however well established that these signals are much smaller than the Faraday effect for two reasons. First the magnetized adiabatic mode of scalar origin dominates at the level of the initial conditions and at the level of the TT correlation: the vector and tensor modes of magnetic origin are comparatively smaller. Second the B mode signals induced by the vectors and the tensors are quadratic in the magnetic power spectrum while the Faraday B mode is linear in the magnetic power spectrum.
The realistic situation is therefore the one where there are two physically plausible sources of B mode polarization: the first is given by the tensor modes of the geometry, i.e. relic gravitons with present wavelengths comparable with the Hubble radius; the second is given the Faraday rotated E mode polarization. The second part of the present paper dealt with the possibility that two effects can interfere. The B mode polarization of tensor origin is virtually frequency independent. Conversely the Faraday rotated E mode polarization does depend on the frequency.
Elaborating on a recent suggestion the Faraday effect has been treated as a random, stationary and quasi-Markovian process. The stochastic treatment of this phenomenon bears some analogy with the case of synchrotron emission and the obtained results encompass and complement previous analyses where the formation of the Faraday effect has been customarily presented as a purely deterministic process in time. Within this approach a set of scaling laws only involving observable power spectra can be derived. These scaling laws, once applied to observational data at different frequencies, can be used to disentangle the Faraday induced B mode polarization from the tensor B mode.

A Description of the polarization
To avoid lengthy digressions, in this Appendix we are going to report some technical details that can be useful for the interested reader. In general terms the brightness perturbations can be arranged in a column matrix, be it I = (∆ I , ∆ Q , ∆ U ) obeying the following equation where Ω denotes the angular dependence and dΩ ′ = dϕ ′ dµ ′ with µ ′ = cos ϑ ′ ; all the matrix elements of W vanish except two: W Q U = −1 and W U Q = 1. The first term generically denotes the collisionless contribution while the collision term contains the matrix N (Ω, Ω ′ ). The matrix elements entering the collision term of the equation for ∆ I are: the matrix elements entering the collision integral of ∆ Q and ∆ U are:

B Scalar and tensor fluctuations
The scalar, vector and tensor components of the brightness perturbations are affected, respectively, by the scalar, vector and tensor inhomogeneties of the geometry and of the various sources. It must be however stressed that the fluctuations of the geometry affect directly only the evolution of ∆ I while the linear polarizations are indirectly affected by the fluctuations of the geometry through the source terms of the corresponding equations. With these caveats in mind we are first going to deduce the collisionless part of the evolution of the brightness perturbations. Assuming the conformally flat background introduced in Eq. (4.1), the fluctuations of the metric can be written, in general terms, as where δ s , δ v and δ t denote the inhomogeneity preserving, separately, the scalar, vector and tensor nature of the fluctuations. The scalar modes of the geometry are parametrized in terms of four independent functions. In the longitudinal gauge the scalar fluctuations of the metric are The vector modes are described by two independent vectors Q i ( x, τ ) and W i ( x, τ ) subjected to the conditions ∂ i Q i = 0 and ∂ i W i = 0. Also in the vector case it is possible to fix a gauge by choosing, for instance, Q i = 0. Finally the tensor modes of the geometry are parametrized in terms of a rank-two tensor in three spatial dimensions, i.e.
which is automatically invariant under infinitesimal coordinate transformations. Ignoring, for the moment, the collision terms we have that the collisionless parts of the evolution of ∆ I can be written as L (s) where q =n i q i denotes the comoving three-momentum. The scalar, vector and tensor contributions to the derivatives of the modulus of the comoving three-momentum are given, respectively, by (B.10) The notation introduced in Eqs. (B.5), (B.6), (B.7) for the fluctuations of the intensity can also be generalized to the linear and circular polarizations: (B.11) where the subscript can coincide, alternatively, with Q, U and V (i.e. X = Q, U, V ) and the superscript denotes the transformation properties of the given fluctuation (i.e. y = s, v, t). The vector and the tensor polarizations can be decomposed, respectively, as wherek denotes the direction of propagation and the two orthogonal directionsâ andb are such thatâ ×b =k. Given the direction of propagation of the relic tensor oriented alongk, the two tensor polarizations are defined in terms ofâ i andb i as: 14) The projections of the vector and of the tensor polarizations on the direction of photon propagationn are: Choosing the direction of propagation of the relic vector and of the relic tensor along theẑ axis, the unit vectorsâ andb will coincide with the remaining two Cartesian directions and the related Fourier amplitudes will satisfŷ Similarly Eq. (B.11) becomes, in Fourier space, The explicit form of the transport equations for the scalar and tensor modes of the geometry will be scrutinized in the remaining part of this section. The vectors as well as the circular polarization will be neglected since they play no role in the present considerations.

C Stochastic Faraday Rotation C.1 Solutions of the evolution equations
The solutions of the evolution equations discussed in sections 4 and 5 can be obtained with the techniques already discussed in [23]. For equal times (but for different Fourier modes) the fluctuations of the brightness perturbations are random with the power spectrum determined by the (nearly scale-invariant) spectrum of (Gaussian) curvature perturbations [24]. Thus, in the absence of Faraday mixing, δ (s) ± obeys then a deterministic evolution in time while the spatial fluctuations of the polarization are randomly distributed and fixed by the correlation properties of the adiabatic curvature perturbations. Conversely since X F (τ ) is now treated as a stochastic process, the evolution equation for the polarization becomes a stochastic differential equation [17,18,19] in time and its formal solution is obtainable by iteration: For the sake of simplicity the index referring to the scalar nature of the brightness perturbations has been dropped from δ ± but has been kept in the source term. We shall just restore the superscript at the very end. Neglecting the corrections O(f 2 e ), Eqs. (4.19) and (C.1) imply the following recurrence relations: The differential optical depth directly enters the visibility function giving the probability that a photon is emitted between τ and τ + dτ : The full solution of Eq. (4.19) is formally expressible as: The same formal solutions discussed in the scalar case can also be obtained in the tensor case. In particular we shall have that Eq. (4.33) can be solved by iteration as which is formally similar to Eq. (C.5) even if the source terms are different. In Eqs. (C.5) and (C.6) the visibility function enters the integrand. For analytic estimates the visibility function has the approximate shape of a double Gaussian whose first peak arises around last scattering (i.e. for τ ≃ τ r ) while the second (smaller) peak occurs for the reionization epoch. The way the visibility function can be analytically approximated has been the subject of different studies [27,28,29] that are relevant for a refinement of the present discussion. However, as we shall argue, the scaling properties of the nonlinear combinations of the polarization observables shall not be crucially affected by the details of the visibility function. The finite thickness of the last scattering surface does not affect the ratios between the different combinations of polarization power spectra discussed here so that the limit of sudden recombination can be safely adopted; in this limit the first and more pronounced Gaussian profile tends to a Dirac delta function.

C.2 Cumulant expansion
The randomness implies that X F (τ ) is not a deterministic variable but rather a stochastic process which is stationary insofar as the autocorrelation function Γ(τ 1 , τ 2 ) = X F (τ 1 )X F (τ 2 ) only depends on time differences, i.e. Γ(τ 1 , τ 2 ) = Γ(|τ 1 − τ 2 |); furthermore we shall also assume that the process has zero mean, even if this is not strictly necessary for the consistency of the whole approach. If τ b defines the time-scale of variation of the brightness perturbations of the polarization observables, the physical situation investigated here corresponds to τ b ≫ τ c where τ c is the correlation time-scale of X F . In the simplest case of a Gaussiancorrelated process the autocorrelation function Γ(τ 1 − τ 2 ) = F (τ 1 )τ c δ(τ 1 − τ 2 ). If the time scale of spatial variation of the rate is comparable with the time scale of spatial variation of the gravitational fluctuations, X F can be considered only time dependent (i.e. a stochastic process). In the opposite situation the Faraday rate must be considered fully inhomogeneous (i.e. a stochastic field). The statistical properties of A ± follow directly from the correlation properties of X F (τ ). If, for instance, X F (τ ) obeys a stationary and Gaussian process, for any set of n Faraday rates (characterized by different conformal times) the correlator X F (τ 1 ) X F (τ 2 ) . . . X F (τ n ) vanishes if n is odd; if n is even the same correlator equals: where the sum is performed over all the (n−1)! pairings. In the Gaussian case, the evaluation of the averages can be performed by first doing the standard moment expansion and then resumming the obtained result. As an example, from the explicit expression of A ± , it follows that where ω F is given by It follows from Eq. (C.9) that even if X F ≤ 1, ω F is not bound to be smaller than 1. The result of Eq. (C.8) holds in an approximate sense when the stationary process is only approximately Markovian, While the standard moment expansion can be formally adopted in specific cases (like the Gaussian one) it cannot be used to provide successive approximations. The reason is that any finite number of terms constitutes a bad representation of the function defined by the whole series. This difficulty is overcome with the use of the cumulants that are certain combinations of the moments. Dropping the functions and keeping only their corresponding arguments we have that the relations between the ordinary moments and the cumulants (denoted by ... ) is 10) and so on and so forth for the other moments of the cluster expansion. Substituting the naive moment expansion with the cumulant expansion we have that the average of Eq. (C.8) is given by A ± (τ, τ r ) A ± (τ, τ r ) = exp ∞ m=1 (±4i) m m! τ τr dτ m X F (τ 1 ) X F (τ 2 ) . . . X F (τ m ) . (C.11) As firstly suggested by Van Kampen (see Ref. [17,19]) in the approximately Markovian case the averages of certain stochastic processes will be given by an exponential whose exponent is a series of successive cumulants of X F . All the cumulants beyond the second are zero in the case of an exactly Gaussian process and the result reported in Eq. (C.8) is recovered. Since each integrand in (C.11) virtually vanishes unless τ 1 , τ 2 ,..., τ m are close together, the only contribution to the integral comes from a tube of diameter of order τ c along the diagonal in the m-dimensional integration space. More generally, the m-th cumulant vanishes as soon as the sequence of times τ 1 , τ 2 ,...,τ m contains a gap large compared to τ c . This is the reason why, in a nutshell, a cumulant is also rather practical in our case. If the stochastic process is not homogeneous the iterative solution of Eqs. (C.1) and (C.2)-(C.3) can be written as: (1 − µ 2 )ǫ ′ S P ( k, τ ), (C.12) ± + (ikµ + ǫ ′ )δ (1) ± + (ikµ + ǫ ′ )δ (2) ± = ∓i b F (ν, τ ) d 3 p ′ δ P ( k + p ′ , τ ) n i B i ( p, τ ) n j B j ( p ′ , τ ), (C.14) To compute the averages we must therefore specify the correlation properties of the Faraday rate. Even if the spatial dependence may reside in all the terms contributing to the Faraday rate, it is reasonable to presume that the leading effect may come from the magnetic field whose correlation function will then be parametrized as: B i ( q, τ 1 ) B j ( p, τ 2 ) = 2π 2 p 3 P ij (p) P B (p) Γ(|τ 1 − τ 2 |) δ (3) ( q + p), (C. 15) where Γ(|τ 1 − τ 2 |) = τ c δ(τ 1 − τ 2 ) in the delta-correlated case. In the same approximation exploited before and using Eq. (C.15), ω F becomes now ω F = 8b 2 F 3ν 4 dp p P B (p) τ τr dτ 1 τ τr dτ 2 Γ(|τ 1 − τ 2 |) a 2 (τ 1 )a 2 (τ 2 ) , (C. 16) where the constant b F = b F (ν, τ )a 2 (τ )ν 2 has been introduced in order to draw special attention to the frequency scaling which is the most relevant aspect of Eq. (C. 16), at least in the present approach.