Magneto-optic effects of the Cosmic Microwave Background

Generation of magneto-optic effects by the interaction of the CMB with cosmic magnetic fields is studied. Effects which generate polarization such as the Cotton-Mouton effect, vacuum polarization and photon-pseudoscalar mixing in external magnetic field are studied. Considering the CMB linearly polarized at decoupling time, it is shown that photon-pseudoscalar mixing in external magnetic field, the Cotton-Mouton effect in plasma and the vacuum polarization in cosmic magnetic field, would generate elliptic polarization of the CMB depending on the photon frequency and magnetic field strength. Among standard magneto-optic effects, the Cotton-Mouton effect in plasma turns out to be the dominant effect in the generation of CMB elliptic polarization in the low frequency part $\nu_0 \sim 10^8-10^9$ Hz with degree of circular polarization $P_C(T_0)\simeq 10^{-10}-10^{-6}$ for magnetic field amplitude $B_{e0}\sim 100\, \text{nG}-1$ nG. The vacuum polarization in magnetic field is the dominant process in the high frequency part $\nu_0\geq 10^{10}$ Hz where the degree of circular polarization at present is $P_C(T_0)\lesssim 10^{-11}$ in the best scenario. It is shown that photon-pseudoscalar particle mixing in cosmic magnetic field generates elliptic polarization of the CMB as well and even in the case of initially unpolarized CMB. New limits/constraints on the pseudoscalar parameter space are found. By using current limit on the degree of circular polarization of the CMB, the upper limit of $|g_{\phi\gamma}|<4.29\times 10^{-19}(\text{G}/B_{e0})$ GeV$^{-1}$ for $m_\phi<1.6\times 10^{-14}$ eV in the weak mixing case is found. If $|g_{\phi\gamma}|<1.17\times 10^{-24}(\text{G}/B_{e0})$ GeV$^{-1}$, a value of the order $|g_{\phi\gamma}|\simeq 10^{-26}(\text{G}/B_{e0})$ GeV$^{-1}$ for $m_\phi\simeq 1.6\times 10^{-14}$ eV in the resonant case, from large scale temperature anisotropy is obtained.


Introduction
The interaction of light with matter and fields has been intensively studied in the literature and first quantitative studies dates back to Galileo, Newton, Faraday and Maxwell. Among the interesting effects that such interaction represents, there is one class of phenomena which includes the interaction of light (electromagnetic wave) with external electromagnetic fields. These phenomena manifest when an electromagnetic wave propagates through an external electromagnetic field that has been altered by the presence of the incident electromagnetic wave. If there is present only an external electric field, the effects that manifest are called electro-optic effects. Instead, if there is present only an external magnetic field, the effects that manifest belong to the category of the magneto-optic effects. In this work I study only the last effects.
Magneto-optic effects not only are important to the established physics but also allow to investigate new effects that have not been found yet. They are generally divided in three main categories that are related to transmission, reflection and absorption of the incident light by the magnetized medium. Depending on the initial polarization of the electromagnetic wave, there are essentially four magnetic-optic effects which belong to the transmission (and not only) category, the Cotton-Mouton (CM) effect, the Faraday effect and two more exotic effects which are the vacuum polarization and the mixing of photons with pseudoscalar (and also scalar) particles in external magnetic field. The reflection category includes essentially only the Kerr effect while the absorption category includes the so called molecular circular dichroism in gases and as will be shown in this work also the photon-pseudoscalar mixing in magnetic field.
In the transmission category, the CM effect has been extensively studied in the literature. It has been experimented mostly in gases, liquids, solids and to some degree even in plasma. The CM effect manifest when light propagates in a magnetized medium where the external magnetic field has a transversal component with respect to the direction of light propagation [1]. This effect also shares a property with more exotic phenomena such as vacuum polarization (or simply QED effect) and photon-pseudoscalar mixing in magnetic field. All three effects manifest only where there is a transversal component of the external magnetic field with respect to the direction of light propagation.
The vacuum polarization has been first proposed and studied in Ref. [2] and since then has received much attention from both theory and experimental physics. This effect would manifest as a phase shift between the two photon states perpendicular and parallel to the external (transverse) magnetic field that eventually give rise to a birefringence effect which intensity depends on the incident electromagnetic wave frequency. One of the most important achievement from the experimental side, is to measure the acquired QED ellipticity angle of the incident light propagating through the magnetic field. Indeed, this has been the quest for the PVLAS [3] experiment, BFRT [4] experiment and for new generation of experiments [5]. After a first claim of detection of vacuum birefringence by PVLAS experiment [6], there is still a long way to achieve the required apparatus sensitivity in order to measure the QED predicted ellipticity which is by more than an order of magnitude smaller than the current apparatus sensitivity. At current status, apparatus sensitivity is contaminated with not well understood background noise, must probably from the same apparatus and new methods have also been proposed [7].
The birefringence effect predicted by QED can also be mimicked by another magneto-optical effect, namely the photon-pseudoscalar/scalar mixing in magnetic field. In fact, as it will be shown in this work, mixing of photons with pseudoscalar particles gives rise to both birefringence and dichroism effects. Therefore experiment such as PVLAS, BFRT etc., can in principle find pseudoscalar particles such as axions, ALPs, scalar bosons if the induced birefringence or dichroism signal is bigger than the QED expected signal. Other important experiments that aim to find exotic pseudoscalar particles include the CAST and IAXO experiments [8], ADMX experiment [9] and ALPS-II [10].
Among all magneto-optic effects, the Faraday effect has received much attention in astronomy and cosmology. It manifest when an initial linearly polarized electromagnetic wave interacts with an external magnetic field that has a longitudinal component along the wave propagation direction. This coupling makes possible the rotation of the polarization plane of the incident electromagnetic wave and the rotation angle is proportional to B e d where B e is the strength of the external magnetic field and d is the length of the path. Consequently, the Faraday effect has been widely used in radio astronomy as a probe of cosmic magnetic fields, in galaxy clusters and also in the intergalactic space [11]. Measurements of the rotation angle of light received from galaxy clusters confirm the presence of a magnetic field inside them, with a magnitude of about few µG. In the intergalactic space, present studies would suggest a weaker large scale magnetic field with upper limit magnitude B e 3 − 1380 nG, see for example current limits by Planck collaboration [12] where limits of the order of 1380 nG are set from Faraday effect. On the other hand non observation of gamma rays emission from intergalactic medium due to injection of high energy particles by blazars would suggest a lower value on the strength of extragalactic magnetic field B e ≥ 10 −16 − 10 −15 G [13]. The origin of this field is still unknown and present studies suggest that it may have been created during structure formation or it may have a primordial origin, see Ref. [14] for a review on cosmic magnetic fields structure and current updated limits/constraints by Planck collaboration. In this work it is assumed that the magnetic field has a primordial origin and its amplitude is a slowly varying function of space-time coordinates, namely a slowly varying inhomogeneous field in space and time which can be also stochastic in nature, see sec. 7 for details.
In connection with the CMB physics, the Faraday effect has been used to probe the existence of primordial magnetic field [15] present at the decoupling time since it would rotate the polarization plane of the CMB. In fact, it is well known by now that the CMB posses a very small linear polarization that is believed to have been generated at the decoupling time due to Thomson scattering of CMB photons on electrons. Such a polarization is generated because of temperature anisotropies present at the decoupling epoch that eventually generate a position dependent photon intensity on the surface of the last scattering [16]. Consequently, Thomson scattering of an anisotropic background of photons on electrons would generate linear polarization of the CMB with non zero Stokes parameters Q and U [17], [18].
In general, the linear polarization pattern of the CMB can be decomposed in two modes with opposite parity, the so called E-modes (or gradient modes G) which are the dominant component of the linear polarization and B-modes (or curl modes C) which are the subdominant component of linear polarization, see Refs. [19]. The former are generated only by scalar density fluctuations of the cosmological plasma while the latter are generated by either vector or tensor modes. The generation of B-modes is induced by tensor perturbations (gravitational waves) [20], Faraday rotation of the CMB [15], gravitational lensing of the E-mode component [21] and due to primordial magnetic fields [22] via perturbations sourced by the magnetic field. In general, the spectrum of Bmodes is described in multipole moments l of spherical harmonics used to describe linear polarization. The location of the peak signal of B-modes as function of l would give the possibility to distinguish between signals generated by different sources of B-modes.
So far, much of attention on the CMB polarization has been focused mostly on the linear polarization. This fact, partially has been influenced by the first experimental observation of E-modes (due to primordial adiabatic scalar fluctuations) by DASI, WMAP and BOOMERANG collaborations [23] and also by the fact that many inflationary models predict an almost scale invariant spectrum of gravitational waves, which as already mentioned above, can produce B-modes which are believed to be the 'holy grail' of the inflationary theory. Moreover, since Thomson scattering is the most frequent type of scattering in the early universe and because it generates only linear polarization, other types of CMB polarization have been to some extent obscured and the V Stokes parameter has become essentially the 'lost along the way' parameter. However, it is well known that light can have two additional types of polarization, circular and elliptic which translate into a nonzero Stokes parameter V .
After this premise on the CMB linear polarization, several questions come spontaneously. Does the CMB posses only linear polarization? Does it have any degree of circular polarization? If yes, what are the generating mechanisms? Even though, there is not urgency on the study of CMB circular polarization, since the discovery of the CMB, there have been several attempts in the past and also at the present to experimentally measure it. Moreover, since CMB linear polarization has already been detected, the next step would be that of the study of circular polazation which as I will show in this paper is generated by very interesting mechanisms which are extremely important to the fundamental physics.
The first studies on the CMB circular polarization were done in connection with studies on anisotropic expansion of the universe which are characterized by some type of Bianchi models [24]. Other studies on generation of CMB circular polarization include non commutative geometry [25], electron-positron scattering in magnetized plasma at decoupling time [26], propagation of CMB photons in magnetic field of supernova remnants of the first stars [28], photon-pseudoscalar mixing in magnetic field domains [29], scattering of CMB photons with cosmic neutrino background [30]. For a recent review on other CMB circular polarization mechanisms see Ref. [31]. The first experimental attempts to measure the circular polarization of the CMB were done in Ref. [32] where no evidence for CMB circular polarization was found and only constraint on the degree of circular polarization was set. The current upper limit on the CMB circular polarization has been set by the MIPOL experiment [33], P C (T 0 ) 7 × 10 −5 − 5 × 10 −4 at the frequency 33 GHz and at angular scales between 8 • and 24 • .
In this work I study the impact of magneto-optic effects on the CMB polarization in the presence of cosmic magnetic fields where I mostly concentrate on generation of CMB circular polarization. A systematical study of the most important magneto-optic effects in the generation of a net CMB elliptic (circular and linear) polarization is done. By including all magneto-optic effects mentioned above, I derive the equations of motion for the Stoke's parameters which form a coupled system of differential equations. I use a density matrix approach to study the mixing of different magneto-optic effects and then solve the equations of motion by using perturbation theory. It turns out that among CM and vacuum polarization effects, the CM effect in plasma is the most promising effect in generation of elliptic polarization in the low frequency part of the CMB, while in the high frequency part, the vacuum polarization is the dominant one. I also will use current limit on the degree of circular polarization, to set new limits on the mass and coupling constant of pseudoscalar particles. In connection with CMB circular polarization, I calaculate its magnitude in terms of degree of circular polarization at present P C (T 0 ) and compare with experimental result(s). Generation and evolution of CMB E-mode and B-mode generated by the above mentioned effects is not studied in this work. This paper is organized as follows: In Sec. 2, I derive the equations of motion for the photon and pseudoscalar fields in an expanding universe and introduce the photon polarization tensor in magnetized medium which describes forward scattering of photons. In Sec. 3, I study the equations of motion for the density matrix in the case of open systems and establish the connection between the system Hamiltonian and the field mixing matrix. In Sec. 4, I find the equations of motion for the density matrix in an expanding universe and solve them in the case of vacuum polarization and CM effects. In Sec. 5, I present the equations of motion for the density matrix in the case when the contribution of the pseudoscalar field is included and introduce the concept of generalized Stokes parameters. Then I find perturbative solutions of the reduced Stokes vectors in transverse magnetic field. In Sec. 6, I study the generation of CMB circular polarization in the case of photon-pseudoscalar particle mixing in transverse magnetic field and set new limits on the pseudoscalar parameter space. In Sec. 7, I conclude. In this work I use the metric with signature η µν = diag(1, −1, −1, −1) and work with the natural (rationalized) Lorentz-Heaviside units (k B = = c = ε 0 = µ 0 = 1) with e 2 = 4πα.

Equations of motions in an expanding universe
In this section we derive the equations of motion for the photon and pseudoscalar fields propagating in a magnetized medium in the framework of the Friedemann-Robertson-Walker (FRW) metric. To start with, we write the effective action of the photon and pseudoscalar fields in curved spacetime where F µν is the total electromagnetic field tensor, Π µν is the photon polarization tensor in medium, φ is the pseudoscalar field, m φ is the mass of the pseudoscalar field, g is the metric determinant and A µ is the photon vector potential. By varying the action with respect to the electromagnetic field A ν and pseudoscalar field φ, the equations of motion are where ∇ µF µν = 0, = ∇ µ ∇ µ is the d'Alambertian operator in curved space, x µ = (t, x), ∇ µ is the covariant derivative and ∇ µ φ = ∂ µ φ. In this work we consider the case of flat (κ = 0) FRW metric with line element ds 2 = dt 2 − dx 2 (t), where t is the cosmological time and x is the spatial coordinate. The only non zero components of the affine connection in the FRW metric are Γ i 0j = (ȧ/a)δ ij and Γ 0 ij =ȧaδ ij . In general the electromagnetic field tensor F µν is given by the sum of the field tensor of incident photon field and of field tensor corresponding to the external magnetic field. In most cases the electromagnetic field tensor corresponding to the external magnetic field is the dominant term. Considering the photon propagation in external magnetic field, the equations of motion (2.2) for the components of vector potential A i and pseudoscalar field φ in the Coulomb gauge 1 are We may notice that there is an extra term in the equations of motion (2.3) with respect to the Minkowski flat space-time for the photon and pseudoscalar fields, that is 3H∂ t where H =ȧ/a is the Hubble parameter and a(t) is the cosmological scale factor. This term is the so called Hubble friction that is responsible for the damping of the fields in an expanding universe. We look for single wave vector solutions of Eqs. 2.3 of the form where k is the photon wave vector, e λ is the photon polarization vector and λ is the photon polarization index. For simplicity, we consider an electromagnetic wave propagating along the observer'sẑ axis with k = (0, 0, k), k = |k|. Without any loss of generality we choose the external magnetic field in the xz plane with coordinates B e = (B e sin(Φ), 0, B e cos(Φ)) where Φ is the angle between the magnetic field direction and photon wave vector k, cos(Φ) =B e ·n withn = k/k. Given the symmetry of the problem, only the transverse part of the external magnetic induces photon-pseudoscalar mixing. Inserting the expansion (2.4) into the equations of motion (2.3) we obtain where we used the WKB approximation, namely that ∂ t |A λ | ω|A λ | and ∂ t |φ| ω |φ|. Indeed, this approximation is well satisfied since the variation in time of external potential (that is proportional to the external magnetic field amplitude) due to universe expansion is much smaller than photon/pseudoscalar frequency. We also used the fact that photons and pseudoscalar particles are assumed to be relativistic and expanded the operator (∂ 2 t − ∇) 2k(−i∂ t + k). The system of Eqs. 2.5 can be written in a matrix form as follows 2 Here A + is the photon state perpendicular to the transverse part of B e , A × is the photon state parallel to the transverse part of B e and I the identity matrix. The diagonal elements of the mixing matrix M in Eqs. (2.6) are , while the off diagonal elements are M φγ = g φγ B e sin(Φ)/2 and iM F = −Π 12 /(2k) is the term that corresponds to the Faraday effect. The elements of the photon polarization tensor 3 Π 11 , Π 22 , Π 12 and Π 21 are calculated in momentum space [35] where we took the adiabatic limit t → t. Their expressions will be given explicitly in the next sections. As far as concerns the nature of the large scale magnetic field, in this work is assumed that it has fixed direction in the sky and its amplitude B e is a slowly varying function of space coordinates, namely homogeneous or almost homogeneous in space. On the other hand due to universe expansion, the field amplitude changes in time, namely the large scale magnetic field is non stationary.

Open systems
In this section we consider the case when photons (for example the CMB) are considered to interact with a medium, which for example can be magnetic field and cosmological plasma. Our goal is to find the equation of motion for the density matrix which in general is not trivial. Here we are interested in quantities that are proportional to the amplitude square of the fields and because we want to study the mixing of CMB photons with pseudoscalar particles, the density matrix approach is the most adapted in this situation. Another fact in favour of this approach is that the CMB is almost unpolarized where the statistical mixture is maximal and the description of such a state demands the use of the density matrix. In the case when a system couples to another system, we are dealing with open systems that exchange energy and matter between each other. Therefore, in the case of photons interacting with plasma and magnetic field, the photon number is not in general conserved and the most important processes that can change their number, in the case that we treat in this work, is photon-pseudoscalar particle mixing in the cosmological plasma.
In the general case of an open quantum system, the equations of motion for the total density matrix, in the Schrödinger picture, are given by the von Neumann equation where ρ is the total density matrix of the system and H T is the total Hamiltonian (not to be confused with the Hubble parameter of the next section). The total system, is in general the sum of a quantum system S which is coupled to another quantum system B which is called the environment or bath, namely S + B. The total system considered here is assumed to be closed, following Hamiltonian dynamics. The state of the system S, which we call the photon-pseudoscalar system, will change as a consequence of its internal dynamics and because of the interaction with its surroundings. The interaction leads to system-environment correlations, such that state changes of S, can no longer be represented in terms of unitary Hamiltonian dynamics. In this context, the photon-pseudoscalar system S is also called a reduced system. Suppose that H S is the Hilbert space of the photon-pseudoscalar system S and H B is the Hilbert space of the environment. The Hilbert space of the total system S + B would be the tensor product H = H S ⊗ H B and the total Hamiltonian has the general form H T = H S ⊗ I B + I S ⊗ H B + H I (t), where H S is the free Hamiltonian of the reduced system, H B is the free Hamiltonian of the environment, H I is the interaction Hamiltonian between the two systems S and B and I B , I S are identity operators in their corresponding Hilbert spaces. If we are interested in the observables of only system S, we can define the density operator of such system by taking the partial trace on the total density operator of the system ρ as follows where ρ S is the density operator of the system S (photon-pseudoscalar system) and Tr B is the partial trace over the environment degrees of freedom. Inserting (3.2) into the von-Neumann equation we get 3), generally depends on different processes that appears in a specific problem and on type of fields that interact with the system. In our case we deal with photons that interact with different particle fields in the cosmological plasma, such as electrons, positrons, protons, light nuclei, cosmic magnetic field and in principle with other exotic particles. We refer to these fields as background fields and the calculation of the expression Tr B [H T , ρ] would be quite involved. In fact, as one may realize at this point, there are essentially two ways on writing down the equations of motion for ρ S . The first possibility would be to start from the general expression (3.3) and use the Hamiltonian of the total system and calculate the commutator with ρ by taking the partial trace over B. The second possibility would be to start with the effective action and derive the equations of motions for the fields by including the effective polarization tensor for photons and their interaction with the pseudoscalar field. In the latter case, one can derive a Schrödinger type equation, which dynamics is governed by an effective Hamiltonian that is given by the mixing matrix M and a 'damping' term due to the Hubble friction as in (2.6). Obviously, the second method is more convenient since it bypasses all the tedious procedure in calculating the r. h. s. of Eq. (3.3). Similar approach has been widely used also in neutrino physics [37] and it is still the most used approach on calculating oscillation probabilities in presence of damping. However, the second approach mentioned above is an approximation of the first method and should not be sought as the most standard procedure.
All told, we work under the approximation where M is the field mixing matrix that is already 'traced out' since it includes the effect of background fields on photons, photon-pseudoscalar interaction and D is a 'damping' matrix that is given by D = (3/2)HI where H is the Hubble parameter. On the right hand side of (3.4) instead of the total system density matrix appears only the reduced system density matrix ρ S . This is due to the fact that the coupling between S and B is weak such that the influence of S on B is very small (the so called Born approximation). In such case, at a given time t one can approximate 4 ρ(t) ≈ ρ S ⊗ ρ B [38]. Consequently, the equation of motion for the density matrix becomes where the first term in (3.5) describes an unitary evolution and the second term describes the 'damping' of fields in an expanding universe.

Photon polarization effects
In this section we focus on the case when the mixing matrix M is not stationary and look for solutions of equations of motion of the density matrix, Eq. (3.5). Indeed, it is more convenient to work with the density matrix than the wave equation, Eq. (2.6). In this section we consider the case of missing pseudoscalar field. As already mentioned, in the presence of an external magnetic field, excluding for the moment the case of photon-pseudoscalar mixing, there are essentially other three magneto-optic effects which depend on the external magnetic field direction and which are proportional to its strength.
In the presence of vacuum polarization, CM and Faraday effects the equation of motion of the density matrix in terms of the Stokes parameters 5 are given bẏ It is import to stress that this approximation does not imply that there are no excitations in the background fields. 5 For the definition of photon density matrix and its connection with the Stokes parameters, see Appendix A where we wrote the elements of the density matrix in terms of the Stokes parameters. In general, the r. h. s. of (4.1) would depend on the temperature T rather than t. Therefore, expressing the time as t = t(T ), the time derivative in the FRW metric becomes ∂ t = −HT ∂ T where H = −Ṫ /T . The system (4.1) can be written in the following matrix form S (T ) = A(T ) · S(T ) + (3/T )I · S(T ), (4.2) where S is the Stokes vector 6 defined as S = (I, Q, U, V ) T and A(T ) is a matrix defined as The system (4.2) is a first order system of linear differential equations with variable coefficients. Even though the matrix A(T ) that enters (4.2) looks very simple, generally the system (4.2) has no closed form of solutions. However, it is possible to find analytic solutions by using the perturbation theory. Indeed, as we will see in what follows, for the parameter space of the photon/pseudoscalar momentum k and magnetic field strength B e which we study in this work, one has in most cases the condition M F |∆M |. This condition on the other hand depends on Φ and for values of Φ → π/2, the Faraday term vanishes. In this case the condition M F |∆M | would not be valid anymore. Therefore, we focus on for the moment in the case when Φ = π/2 in such way that condition M F |∆M | holds and split the matrix A(T ) in the following way where we wrote ∆M (T )/(HT ) = G (T ) and 1 is a parameter that depends on momentum k, magnetic field strength B e and on the angle Φ. HereG(T ) is a function that depends only on the temperature T . The numerical factor in the product G (T ) is included in the parameter . The expression of will be given in the next sections. The second term which appears in (4.2) corresponds to the Hubble friction and its contribution to S(T ) appears as a damping factor of the form exp[−3 T (1/T )dT ] and it is common to all components of the Stokes vector. The easiest way to see it, is by observing that matrix A(T ) commutes with (3/T )I for every T . For the moment we concentrate on the solution of Eq. (4.2) without the damping term and include it in the final result.
We look for solution of the Stokes vector up to first order in as follows Inserting expansion (4.4) into Eq. (4.2) and collecting the appropriate terms we get the following matrix equations S 0 (T ) = A 0 (T )S 0 (T ), (4.5) We may observe that for different cosmological temperatures, the commutator of [A 0 (T 1 ), A 0 (T 2 )] = 0 which allows us to find the following exact solution for S 0 (T ) where T i is the initial temperature and where T < T i is the CMB temperature after the decoupling time. We may observe that the homogeneous part of Eq. (4.6) has the same solution as Eq. (4.5) with the replacement S 0 (T ) → S 1 (T ). The non homogeneous part of Eq. (4.6) can be solved with the method of the variations of constants. Performing several algebraic operations that involve matrix exponentiations, collecting all the appropriate terms together and including the term corresponding to the Hubble friction we get the following solutions for the components of the Stokes vector to the first order in : There are several interesting considerations that can be made about (4.9)-(4.11). In the first place we may notice that each solution is proportional to the initial values of the Stokes parameters Q i , U i and V i , as one would expect from a first order system of linear differential equations. This implies that if the initial conditions are all zero, as for example in the case of unpolarized light, it would remain unpolarized during the universe expansion. If this is the case, the Faraday effect, the vacuum polarization and CM effect would not have any impact on the CMB polarization at all. The only way that these effects can have an impact on the CMB polarization, would be if the CMB is initially polarized. As already mentioned in the introduction section, Thomson scattering would generate CMB linear polarization only if there are anisotropies in the CMB temperature (or intensity). If the incident light is initially unpolarized and anisotropic, Thomson scattering generates outgoing polarized light with non zero Stokes parameters I and Q while V = U = 0. This is a general property of Thomson scattering for anisotropic incident light. Since the parameters Q and U depends on the coordinate system, one can rotate the system to a common one, in such a way to have U = 0. It can be shown that in the rotated system, the temperature anisotropy of the CMB generates non zero initial Stokes parameters Q i and U i at the decoupling time [17], [39] where σ T is the Thomson scattering cross section, σ B is the cross sectional area of the scattered light and a 22 is the second multipole coefficient used in expanding the incident photon intensity in spherical harmonics Y lm . We have intentionally labelled with i the values of Q and U at the decoupling time and use them as the initial conditions in (4.9)-(4.11). However, as it has been well studied in the literature, Thomson scattering does not generate circular polarization and in the case of CMB this fact is confirmed since equation which governs evolution of V parameter due to Thomson scattering has no source term [17]. Consequently in this section we assume that at decoupling In the second place we may note from (4.9)-(4.10) that due to magneto-optic effects, light during its propagation contemporary has its polarization plane rotated and there is generation of phase shift between photon states A + and A × . This behaviour is well know in optics and the medium which induces such effects is usually referred as rotated retarder, namely it rotates the polarization plane and generates ellipticity at the same time. We may note that for V i = 0 the expressions for Q and U are the same as in the case of solely Faraday effect taken into account. Moreover, as it has been discussed in Sec. 2, in case when light is linearly polarized, ∆M (T ) = 0 the Faraday effect alone does not generate circular polarization. Also we may note the contribution of the Hubble friction term to the Stokes parameters, namely (T /T i ) 3 that for convenience reasons we putted it on the r. h. s. of (4.8)-(4.11). The effective scaling of the Stokes parameters due to universe expansion is not (T /T i ) 3 but (T /T i ) 2 since the Stokes parameters contain an intrinsic scaling of (T /T i ). This is due to the fact that in the WKB approximation the fields have an intrinsic normalization of 1/ ω(t) with ω(t) being photon/pseudoscalar energy. We did not show explicitly this factor for simplicity which gives a multiplicative factor to fields proportional to T 1/2 and to the Stokes parameters proportional to T . Since the scaling factor due to universe expansion is common to all Stokes parameters and because we are mostly interested in their ratio or expressions that contain their ratio, such as polarization degrees P L,C , this term eventually cancels out. For example, the degree of linear polarization of the CMB remains constant during universe expansion to first order in where we took for simplicity I i = 1 and V i = 0 in (4.9)-(4.10). The total rotation angle of linear polarization of the CMB due to the Faraday effect is given by ψ F (T ) = F (T )/2. The contribution of CM and vacuum polarization effects to linear polarization does not appear to the first order in . Their contribution appears only to second order in but for our purposes only expansion to first order is important in this section. Apart the fact that the polarization plane of the CMB is rotated due to the Faraday effect, another interesting effect is the generation of circular polarization with non zero Stokes parameter V (T ). Even in the case when there is not circular polarization at the decoupling time, it is generated afterwards due to vacuum polarization and CM effects. In the case of vanishing V i we have In both terms on the r. h. s. of (4.13) enters the function F (T ) which represents the effect of the Faraday effect. To have an analytic expression for F (T ) we need first the expression for M F (T ) which is given by one of the off-diagonal terms of Π ij . The Faraday effect is induced by the longitudinal component of the magnetic field with respect to k, namely by B L = B e cos(Φ). Consequently, linearly polarized electromagnetic wave propagating along the direction of the external magnetic field, has its polarization plane rotated with an angle proportional to B L . This occurs because the right and left handed indexes of refraction n R and n L are different from each other, which make possible mixing between linearly polarized states A + and A × . The expression for the Faraday term is given by where ω 2 pl = 4παn e /m e is the plasma frequency, n e is the free electron number density, m e is the electron mass and we used k ω for photons. Here ω c = eB e /m e is the cyclotron frequency with e being the electron charge. During propagation of the electromagnetic wave in a magnetized medium, the wave polarization remains unchanged for initial linearly polarized wave, but the linear polarized states A + and A × propagate with a new index of refraction ∆n F in the medium which is given by ∆n F = n R − n L .
The last thing that remains to calculate is the expression for the Hubble parameter which enters F (T ) in (4.7) and in G (T ). In general, its expression in the case of zero spatial curvature (κ = 0) is given by where H 0 is the Hubble parameter at the present epoch, , Ω Λ is the density parameter of the vacuum energy, Ω M is the matter density parameter and Ω R is the density parameter of relativistic particles. According to the Planck collaboration [40], values of density parameters of nonrelativistic matter and vacuum energy are respectively h 2 0 Ω M = 0.12 and Ω Λ = 0.68 with h 0 = 0.67. The density parameter of relativistic particles it is straightforward to calculate, Ω R = 4.15 × 10 −5 h −2 0 which includes the contribution of photons and three neutrino species assumed to be nearly massless. The contribution of the external magnetic field to the energy density budget of relativistic fields can be safely neglected since its energy density is ρ B (T 0 ) 10 −7 (B 0 /nG) 2 ρ γ (T 0 ).

Vacuum polarization in external magnetic field
Having the expressions for the Faraday term M F (T ) and H(T ), we have almost all necessary ingredients to calculate the degree of circular polarization 7 at present time, V (T 0 ). Vacuum polarization and CM effects are responsible for generation of circular polarization. They are induced by the transverse component of the external magnetic field, B T . In both effects the linear polarization indexes of refraction, n + and n × are different from each other. Contrary to the Faraday effect which has its index of refraction proportional to B L , vacuum polarization and CM effects have their indexes of refraction proportional to B 2 T . In this section we consider the contribution of vacuum polarization 8 to V (T ) separately from the CM effect which will be considered in the next section. Vacuum polarization 9 occurs (not only) in the presence of an external magnetic field due to creation of electron/positron pair from the vacuum, see Fig. 1. The expressions of elements of photon polarization tensor 10 corresponding to the states A + and A × in case of vacuum polarization, for slowly varying external magnetic field in space and time over the Compton wavelength, are respectively given by [41] where κ = (α/45π)(B e /B c ) 2 and B c = m 2 e /e is the critical magnetic field. Using (4.14) and definitions of M + and M × we get Now using the definition of κ and taking into account that the photon energy scales with the temperature as ω = ω 0 (T /T 0 ) and assuming magnetic field flux conservation in the cosmological plasma with B e (T ) = B e (T 0 )(T /T 0 ) 2 we get where T 0 is the CMB temperature today and ω 0 = 2πν 0 with ν 0 being the CMB frequency at present. Deriving (4.16) we used the fact that ω c ω in the Faraday term, expressed the free electron number density n e (T ) which enters the plasma frequency as n e (T ) = 0.76 n B (T 0 )X e (T )(T /T 0 ) 3 with X e (T ) being the ionization fraction of free electrons and n B (T 0 ) is the baryon number density at present. Moreover, we assumed for simplicity that only non relativistic matter contributes to the Hubble parameter, namely Until know we did not give any explicit expression for the parameter QED . This parameter can be extracted immediately from G(T ) in (4.16) and it is given by . 8 Vacuum polarization in external magnetic field is a non linear QED effect which Lagrangian density is given by the Euler-Heisenberg Vacuum polarization considered in this paper is due to interaction of CMB photons with an external magnetic field which is different from free CMB photon-photon scattering studied in Ref. [27]. 10 The diagonal terms of the polarization tensor include the contribution of plasma effects, vacuum polarization and CM effect. Since plasma effects are the same for A + and A × and because in this section and in the next we calculate, ∆M = M + − M × , the plasma term cancels out.
With expressions of F (T ), G(T ) and QED we have all necessary quantities to calculate V (T ) in (4.13). Let us start with the first term on the r. h. s. of (4.13) which has the following dependence on the temperature where we have defined As one would expect there is an integration in X e (T ) on the r. h. s. of expression (4.17) that complicates the situation quite a lot. Indeed, there is no known analytic expression for X e (T ) which in general satisfies a complicated differential equation, see Ref. [42] for details. At the temperature T 3000 K numerical solution of the equation satisfied by X e (T ), shows that ionization fraction is X e 0.13 and drops down to X e 2 × 10 −2 at the temperature T = 2000 K. When the temperature is about 200 K it drops down to X e 2.7 × 10 −4 and remains almost constant afterwards if no reionization epoch is assumed. In this work we use the solution for X e (T ) given in Ref. [42] and interpolate it with X e 1 for T 21.8 K which corresponds to the period of end of the reionization epoch.
The integral in (4.17) has analytic solution in terms of the incomplete Euler gamma functions, if the expression for X e (T ) is constant 11 . However, one may observe that for T ≤ 2970 K, for magnetic field strength B e0 ≤ 1 nG and frequencies ν 0 ≥ 10 10 Hz, the expression inside sine function is less than unity, namely F (T ) < 1. In this case one can use series expansion and consider only the first term. Consequently we obtain We are interested in calculating the integral in (4.18) at T = T 0 and numerical calculation gives Now it remains to calculate the second term on the r. h. s. of (4.13). Based on the same arguments as we did above for the first term, the argument of cosine function is smaller than unity and we can write The value of V at present time would be for T i = 2970 K (the CMB temperature at the redshift 1 + z = 1090 corresponding to the decoupling time) and T 0 = 2.725 K If for example we take ν 0 = 30 GHz and B e0 = 1 nG, we get We have checked that numerical values derived from expression (4.20), perfectly agree with numerical solutions in the case when one assumes H H 0 √ Ω M (T /T 0 ) 3/2 . The difference between numerical solutions with total H and semi-analytic solutions with only matter contribution to H, is that in the former case numerical solutions are in general smaller by a factor less than √ 2 with respect to the latter case. We may note that the second term in (4.21) is proportional to the frequency and for higher values of ν 0 , V 0 increases linearly with ν 0 .
So far, we have derived our results in the case when Φ = π/2, which allowed us to find perturbative solution for the Stokes vector S(T ). In the case when the magnetic field is transverse, this approximation is not valid anymore since M F (T ) → 0 for Φ → π/2. However, if Φ = π/2 it is not necessary to work with the perturbative approach since the equation for the Stokes vector simplifies significantly. Indeed, for Φ = π/2 there is only mixing between the Stokes parameters U and V . The solution of equations of motion for the Stokes parameters U and V , in transverse magnetic field are immediate and read In order to estimate V (T ) at present time for Φ = π/2, we first must calculate G(T ) in the argument of sine function. Consequently, we get In general for a wide range of the parameters ν 0 and B e0 in (4.23), the argument of sine function is much less than unity and one can replace to first order the sine with its argument. If we take for example the values ν 0 = 100 GHz and B e0 = 1 nG we get If the magnetic field is 10 nG we would get V 0 2.7 × 10 −8 U i . It is worth to note that in case of transverse magnetic field, the vacuum polarization induces also a rotation of the polarization plane. Indeed, as can be inferred from (4.22), the rotation angle of the polarization plane is given by which in general is a very small quantity for vacuum polarization. Until now we kept the dependence on Φ explicit in (4.20) but it is more convenient to average over all possible orientations of B e relative to k. We must note that (4.20) has been derived by assuming Φ = π/2 which allowed us to find perturbative solution for the Stokes vector S(T ). However, we may note that in the limit Φ → π/2, the first term in (4.20) goes to zero while the second term coincides with the argument of sine function in (4.23) which has been found exactly. This fact allows us, to find the following expression for the rms of V 0 in (4.20), which for F (T ) < 1 is given by Assuming for example, B e0 = 1 nG, Q i U i we get for ν 0 = 30 GHz and ν 0 = 700 GHz respectively V 2 We may note from (4.25), that biggest contribution comes for values of Φ → π/2 or transverse fields and for higher values of B e0 , the rms of V 0 is bigger.
Another interesting case is when the arguments of sine and cosine functions in (4.22) are equal to π/2, namely G(T ) = π/2. This condition is fulfilled when When this condition is met, we would have V (T ) = U i . However, in order for condition (4.26) to be fulfilled, the value of ν 0 must be much far beyond the present observed CMB spectrum for reasonable values of B e0 .

The Cotton-Mouton effect
As briefly mentioned in the previous sections, the CM effect is a birefringence effect which is induced in a medium in presence of transverse external magnetic field and which generates elliptic polarization. This effect has been studied and experimented in gases and liquids where limits on the CM constant C CM are set or its value is established. However, the CM effect does not only manifest in gases and liquids but also in plasma 12 . The theory of this phenomena is studied to some extend classically and also quantum mechanically and for a discussion on this mechanism see Ref. [1]. After the decoupling epoch, the ionization fraction of free electrons rapidly dropped down to an almost constant value of X e 2.7 × 10 −4 and later again it increased to X e 1 at the reionization epoch. On the other hand almost all baryons would bind together to form the light elements such as atomic hydrogen and helium etc. This state of mixed hydrogen and helium gas (plus a small fraction of other light elements) with the electron plasma continued coexisting till the start of the reionization epoch. In order to study the impact of CM effect in the generation of CMB polarization, we need the value of ∆M CM for hydrogen and helium gass and for the electron plasma.
Theory of CM effect in gases has been extensively studied in the literature and for a review on the subject see Ref. [43] and references there. In case of gases, theoretical calculations give the following expression for the difference in index of refraction ∆n gas CM [43] ∆n gas CM = πB 2 e sin 2 (Φ) n gas ∆η/(4πε 0 ), (4.27) where n gas is the gas number density and ∆η is called the hypermagnetizability anisotropy. Here it is assumed that n gas obeys the perfect gas law or its closely related ideal gas law. In general ∆η will depend on the type of gas and on the incident energy of the electromagnetic radiation. Quite often the CM constant is also defined through the relation ∆n gas where λ is the wave length of the incident electromagnetic wave (not to be confused with photon helicity state). By comparing (4.28) with (4.27) we get the following expression for C CM = π ∆η n gas /λ. In our case, we are interested in only the magnetic hypermagnetizability of hydrogen and helium gases since these elements are the most abundant ones and neglect the contribution of the other light elements. Following Ref. [43], theoretical values of ∆η in the limit of zero incident photon momentum, gas temperature T gas = 273.15 K and gas pressure P gas = 1 atm are respectively given by ∆η H = 13.33 au and ∆η He = 1.06 au where 1 au of η is 2.682 × 10 −44 (4πε 0 ) G −2 cm −3 [43] with ε 0 = 1 in the rationalized Lorentz-Heaviside system. Consequently The contribution of the electron plasma to the CM effect enters the diagonal elements of the polarization tensor in magnetized plasma, Π 11 and Π 22 . The difference with respect to the Faraday effect is that CM effect is quadratic in the amplitude of transverse magnetic field and one would expect that for typical values of the cosmic magnetic field, its magnitude would be much weaker than the Faraday effect. In general, for a magnetized plasma the contribution of the CM effect to the photon polarization tensor is given by [35] We may note that in (4.30), the CM term appears only in Π 11 (the second term) while it does not appear in Π 22 . This is due to the fact that we have choosen since the beginning the transverse part of B e along the x axis with no y component. Using the definition of ∆M , for the CM effect in magnetized plasma, we get 12 I learned only recently about the CM effect in plasma.
In case of CMB, we have that photon frequency is much bigger than cyclotron frequency and we can approximate If we compare (4.31) with (4.29), we may observe that for the parameter space of magnetic field amplitude B e0 and photon frequency ν 0 of interest, the contribution of hydrogen and helium gases to the CM effect is much smaller than the contribution of electron plasma. Therefore from now we will neglect the gas contribution to the CM effect. Now we can calculate the contribution of CM effect to V (T ) in the same way as we did for the vacuum polarization. In case when Φ = π/2 and F (T ) < 1 at post decoupling epoch, we have Ti T X e (T ) T 1/2 dT dT , where in the second term in (4.32) numerical integration has been used and defined G(T ) = CMG (T ) with CM CM = −1.21 × 10 32 The second term that enters the r. h. s. of (4.13) is given by The total expression for the degree of circular polarization V (T ) at present time is given by The case when the magnetic field is completely transverse, Φ = π/2 is treated in the same way as in the case of vacuum polarization. What we need is to calculate G(T ) in the case of CM effect, which in most cases is 1. We remind that expressions for the Stokes parameters in case of transverse magnetic field are found exactly without using perturbation theory and are given in (4.22). In case when V i = 0, from (4.22) we get where we kept only the first order term in G(T ) 1. The most interesting fact, is the relation V 0 ∝ ν −3 0 in (4.35). If we consider for example B e0 = 1 nG and ν 0 = 10 8 Hz we get V 0 1.2 × 10 −4 U i while for ν 0 = 10 9 Hz we get V 0 1.2 × 10 −7 U i . In principle one can also calculate the rms of V 0 for the CM effect, as we did in case of vacuum polarization, but it is not that easy. Indeed, expression (4.33) has been derived in the approximation when F (T ) < 1, which assuming that B e0 ≤ 1 nG, it is satisfied for ν 0 > 10 10 Hz. However, the biggest contribution in V 0 comes from the low frequency part as we saw for the case when Φ = π/2. Instead of looking for analytic solution even when F (T ) < 1 is not satisfied and after estimate the rms of V 0 , one possible way to circumvent this situation, is to note that rms of V 0 is bigger for values of Φ → π/2. Indeed, we have checked the numerical solution and found that value of V 0 for Φ = π/2, fixed B e0 and ν 0 is much smaller than that of Φ = π/2. Consequently, one can approximate to very good accuracy the rms of V 0 with its value at Φ = π/2.
The low frequency part of the CMB, ν 0 ∼ 10 8 Hz for Φ = π/2 is very interesting since significant rotation of the polarization plane occurs. Indeed, if one keeps G(T ) 1 up to second order in U (T ) in expression (4.22) and assuming that at decoupling time Q i −U i , one gets for the rotation angle the following expression where we wrote ψ(T ) ψ(T i ) + δψ(T ) with |δψ| 1 and used (4.24). If we consider B e0 = 1 nG, ν 0 10 8 Hz we get δψ(T 0 ) 1.8 × 10 −9 rad, for B 10 −8 G and same frequency we would get δψ(T 0 ) 1.8 × 10 −5 rad and for B = 10 −7 G we would get δψ(T 0 ) = 0.18 rad.

Three state density matrix and generalized Stokes parameters
In Sec. 4 we derived temperature (or time) evolution of the Stokes parameters in case when only Faraday effect, vacuum polarization and CM effects were included in the equations of motion of the density matrix. However, there still remain one effect left which is the photon-pseusocalar mixing in magnetic field. Including this new effect, the equations of motion of the density matrix become more involved and instead of four equations which had in Sec. 4, now we have nine of them. This can be verified by inserting the mixing matrix M and the damping matrix D in Eq. (3.5). After we get the following system of differential equationṡ where the sign ( * ) means complex conjugate of a C-number and each element of the mixing matrix M depends on the cosmological time t, B e and ν. We may observe that total intensity is diluted due to universe expansion onlẏ which means that trace of the density matrix is not constant in time.
The system of Eq. (5.1) still is not in the desired form since the photon intensity, that in this section we denote with I γ , is not a conserved quantity. As already discussed in Sec. 4, this is due to the fact that we are dealing with an open system interacting with the background. Even in the case of other magneto-optic effects which we treated in Sec. 4 there is interaction with the background, but with the difference that these effects conserve the photon number with momentum k. Since I γ is not conserved, it would be convenient to express the equations of elements of the density matrix in terms of generalized Stokes parameters 13 , that is extending the usual two state Stokes parameters to the case of three states.
The derivation of generalized Stokes parameters can be done in analogous way as one does with usual Stokes parameters. As shown in Appendix A, one can express the elements of the two dimensional density matrix in terms of the Stokes parameters and one can check from direct calculation that expression (A.3) can be written in terms of the Pauli matrices σ i as follows where we recall that S 0 = I γ , S 1 = U, S 2 = V, S 3 = Q and I 2×2 is the two dimensional identity matrix. The generalization of the usual two state Stokes parameters to the three state case can be done as follows Ŝ k := Tr(ρλ k ), (5.2) where ρ is the 3 × 3 density matrix,Ŝ k (for k ≥ 1) are the generators of SU(3) group and λ k (for k ≥ 1), (k = 0, 1, ...8), are the so called Gell-Mann matrices Inserting λ k into expression (5.2) we can get the explicit expressions for the generalized Stokes parameters. The first set of four Stokes parameters is given by (A.2) while the remaining set of parameters is given by The corresponding Stokes operators to the set (5.3) in the basis |A + , |A × , |φ are given bŷ Having defined the generalized Stokes parameters, now we are at the position to parametrize three state density matrix in terms of them as follows where I is the total intensity which is given by I = I γ + I φ . Using (5.4) we can write the system of Eqs. (5.1) as followṡ

Equations of motion in absence of the Faraday effect
The system of Eqs. (5.5) is in the final form and we can immediately see from the equations of motion governing usual Stokes parameters, the contribution of the pseudoscalar field to the linear and circular polarization. Let us stress since now that an exact closed analytic solution for (5.5) is not possible. However, here we consider some particular cases, by using some reasonable approximations, which allow us to find semi-analytic solutions for Eqs. (5.5). Indeed, the system (5.5) can be simplified by considering the case of transverse external magnetic field, namely Φ = π/2. This can be achieved by observing the CMB in the direction perpendicular to the external magnetic field and for this particular configuration, the Faraday effect would be completely absent.
In the case when the Faraday effect is absent, we get the following systems of decoupled differential equations in the variable T

Solution of first reduced Stokes vectorS 1
Let us focus first on the solution of first reduced Stokes vectorS 1 (T ). We may note that an exact solution is not possible unless one uses some approximations that allow to find the solution by using perturbation theory, in a similar way as shown in Sec. 4. Therefore, we split the matrix B(T ) in the following order, B(T ) = B 1 (T ) + B 2 (T ) where we recall that ∆M 1 ≡ M + − M φ = M QED Here M pl = −ω 2 pl /(2ω) is the term corresponding to plasma effects which is the same for A + and A × . In order to use perturbation theory, first we must establish which part of the matrix B can be treated as small perturbation.
Suppose first that matrix B 2 (T ) can be considered as perturbation matrix, namely we can write it as the product of a small temperature independent parameter, , with temperature depended functionsG(T ) andG 1 (T ). This situation would be true when either |∆M (T )| < |∆M 1 (T )| M φγ (T ) or |∆M 1 (T )| < |∆M (T )| M φγ (T ). We will find the corresponding parameter space in Sec. 5. Using the same formalism as we showed in Sec. 4, we get the following solutions (to the first order in ) for U and V components ofS 1 where we have defined F φγ (T ) and G 1 (T ) respectively as and ∆G(T ) = G(T ) − G 1 (T ). Even though does not appear explicitly in (5.8), it is implicitly included in G(T ) and G 1 (T ). In (5.8) we show only the solutions for U and V and do not show those for the other components of S 1 since we are not interested in 15 . So far we found the solution forS 1 in the case when elements of the matrix B 1 (T ) are much bigger in magnitude than elements of B 2 (T ), where the last matrix has been considered as perturbation matrix. However, for some values of the parameters we have also the situation when |∆M (T )| < M φγ (T ) |∆M 1 (T )|. Here we are mostly interested in the case when the pseudoscalar mixing term is bigger than |∆M (T )| because the opposite case is fulfilled for uninteresting small values 16 of g φγ . In the case when |∆M (T )| < M φγ (T ) |∆M 1 (T )| it is convenient to move the term ∆M (T ) from matrix B 2 (T ) to matrix B 1 (T ). In this case the former matrix has non zero entries only ∆M 1 (T ) while the latter matrix has non zero entries M φγ (T ) and ∆M (T ). Now, the matrix B 2 (T ) can be considered as the leading one while B 1 (T ) can be considered as perturbation matrix. However, since M φγ (T ) appears now in B 1 (T ), in order to see the small effects of the pseudoscalar field, it is necessary to look for solution to the second order in , namely we writeS 1 (T ) = S . After collecting all terms and tedious calculations we get the following perturbative solutions for U and V components ofS 1 to second order in 15 In this paper we are only interested in usual Stokes parameters Iγ , Q, U and V since they completely describe the polarization of light. If one is also interested in intensity of pseudoscalar field I φ which is related to S 8 or transition amplitudes of photons into pseudoscalar particles then are needed also expressions for the remaining Stokes parameters S 4 , S 5 , S 6 , S 7 , S 8 . We show their solutions in Appendix A. 16 The case M φγ (T ) |∆M (T )| essentially means that contribution of pseudoscalar field to the mixing is smaller than the sum of QED and CM effects. Since the last effects are very small in general, see Sec. 4, the case M φγ (T ) |∆M (T )| is not of particular interest because it is satisfied for extremely small values of g φγ . If indeed g φγ is so small, it would be very difficult to experimentally detect pseudoscalar particles, because their signal would be smaller than the QED effect even if perfect laboratory vacuum is achieved.
where we have defined G 1 (T ) and G φγ (T ) respectively as

Solution of second reduced Stokes vectorS 2
Now we focus on the solution of second reduced Stokes vectorS 2 which is the only one left. Even in this case we look for approximate solution and use perturbation theory in analogous way with the previous section. It is convenient to split the matrix C(T ) which enters in the second equation in (5.6) in the following order, C(T ) = C 1 (T ) + C 2 (T ) At this point we must establish which matrix in (5.10) can be considered as perturbation matrix. This can be done by comparing the elements of C 1 (T ) with C 2 . In the case when |∆M 2 (T )| M φγ (T ), the matrix C 2 (T ) can be considered as perturbation matrix and vice-versa in the case |∆M 2 (T )| M φγ (T ). In case when |∆M 2 (T )| M φγ (T ), we get the following solutions 17 for I γ and Q components ofS 2 to first order in where we defined G 2 (T ) = G 2 (T ) = ∆M 2 (T )/(HT ). As in the previous section does not explicitly appear in (5.11) but is implicitly included in G 2 (T ). The case |∆M 2 (T )| M φγ (T ), needs a special treatment because the term corresponding to the pseudoscalar field is subdominant. In order to explore the vast region of pseudoscalar particles parameter space, it is necessary to look for solution ofS 2 (T ) up to second order in . Therefore, we expand the second reduced Stokes vector as S 2 (T ) =S Collecting all terms, we get the following perturbative solutions for I γ and Q components ofS 2 to second order 18 in : where G 2 (T ) = Ti T G 2 (T )dT .

Pseudoscalar particle production and generation of CMB polarization
In Sec. 5 we solved the equation of motion for the reduced Stokes vectors in case of perpendicular propagation with respect to the external magnetic field B e . This particular configuration, allowed us to solve the equations of motion by using perturbation theory. In this section we focus on the impact of pseudoscalar particle production in generation of CMB circular polarization. In what follows, we concentrate only in generation of CMB polarization after the decoupling epoch and estimate the degree of circular at present epoch 19 . In order to make our treatment as simple as possible, here we concentrate for simplicity only in the case when M φγ (T ) is subdominant with respect to the other terms ∆M 1 (T ) and ∆M 2 (T ).

Subdominant pseudoscalar contribution: M φγ
In the case when the term M φγ (T ) is smaller than other terms in matrices B(T ) and C(T ), which essentially corresponds to the weak mixing case, namely M φγ (T ) |∆M 1,2 (T )|. In this case the expressions for the Stokes parameters I γ and Q are given by (5.12) while for U and V are given by (5.9). Let us concentrate at the post decoupling epoch and assume that at T = T i the CMB is very weakly polarized due to Thomson scattering and consider the generation and evolution of polarization for T ≤ T i . In what follows, we assume that the cosmological plasma is not populated by other relic pseudoscalar particles at T ≤ T i , the CMB is not circularly polarized at T = T i and significant pseudoscalar particle production starts at T = T i in already existing cosmic magnetic field. In this case we have, V i = S 4i = S 5i = S 6i = S 7i = 0 and conservation of particle number gives I γ (T i ) = √ 3S 8i . We use these values as initial conditions in expressions (5.9) and (5.12). In what follows we are not interested in the evolution of other Stokes parameters and will not be considered.
It is important at this stage to find the pseudoscalar parameter space that satisfy the condition of weak mixing. M pl − M φ where here we are also assuming that |M φ | is bigger than QED and CM terms. Therefore, it remains to confront the term M φ with the plasma term M pl , where the former depends on the pseudoscalar mass m φ . Therefore we have either |M φ (T )| > |M pl (T )| or |M φ (T )| < |M pl (T )|. Consequently, depending on the pseudoscalar mass, we have respectively either |∆M 1 (T )| |M φ (T )| or |∆M 1 (T )| |M pl (T )|. Being the plasma term much bigger than QED and CM terms, the previous conditions imply that in most cases we have |∆M 1 (T )| > |∆M (T )|. Based on the same arguments one can easily show that also ∆M 2 M pl − M φ .
In the weak mixing case we have |∆M (T )| < |∆M 1 (T )| since |M × (T )| = |M φ (T )| and therefore it remains to find the parameter space only for M φγ (T ) |∆M 1,2 (T )| where in most practical cases we have |∆M 1 | |∆M 2 | for |M φ | = |M pl |. We find that conditions |M φ | > |M pl | and M φγ |∆M 1,2 | are satisfied for all T 0 ≤ T ≤ T i at the post decoupling if 20 2 × 10 −10 eV < m φ , g φγ 9.57 × 10 15 Hz ν 0 In the case when |M φ | < |M pl | we have M φγ |∆M 1,2 | |M pl | which are satisfied for all T at post decoupling epoch if whereX e is the average value of X e (T ) at the post decoupling epoch. The reason of having chosen the average value will be clear below. Now that we have established limits of validity for I γ and Q (M φγ (T ) |∆M 2 (T )|) and for U and V (M φγ (T ) |∆M 1 (T )|) in the weak mixing, we can focus on the generation of CMB circular polarization. In the Stokes parameters I γ and Q do appear trigonometric functions that have as argument G 2 (T ) while in U and V have as argument G 1 (T ). Since these functions are given respectively by the integral of ∆M 2 (T ) and ∆M 1 (T ), it may be convenient to separate if either the plasma term M pl or the mass term M φ dominates in ∆M 1,2 . As already mentioned, the QED and CM terms are much smaller than the plasma term. Consider first the case when the plasma term M pl dominates in ∆M 1,2 . Considering only the matter contribution to the Hubble parameter H(T ), we get We may note that G 1,2 is given as the integral of the inverse square root of temperature times the ionization fraction X e . As already discussed in Sec. 4, there is not an analytic function for X e which satisfies a complicated differential equation. In Sec. 4 we were able to find semi-analytic solutions for integrals involving trigonometric functions which have as argument integrals of X e . For most practical cases, the argument of those trigonometric functions was much smaller than unity, but here we may note that G 1,2 is never less than unity for realistic values of ν 0 . So, in this section we cannot approximate the cosine or sine of G 1,2 with unity or G 1,2 to first order.
It is desirable to have analytic or at least semi-analytic expression for the degree of circular polarization as we did in Sec. 4. Since there is no known analytic expression for X e , it is convenient to replace it with its average value in G 1,2 at post decoupling epoch, namelyX e 0.023. PuttingX e into G 1,2 we obtain G 1,2 (T ) . Now with G 2 (T ) given, we can calculate the integrals which appear in I γ in expression (5.12). By integrating, we obtain 3) where we have definedã ≡ 1.19 × 10 18 (Hz/ν 0 ) and b ≡ 6.47 × 10 43 g φγ /GeV −1 2 (B e0 /G) 2 . Defining y = I c + I s , we get the following expression for the intensity at present 84ã]) and I i = 1. Let us concentrate on V parameter in (5.9) and on the first term proportional to U i , since other terms are absent with our choice of initial conditions. We may note the first term within parenthesis which corresponds to the QED and CM effects while other terms correspond to mixing of pseudoscalar term with QED and CM terms. The first thing to point out, is that the term corresponding to the QED and CM effects is smaller than other terms because we are in the situation when |∆M | is smaller than M φγ . The second thing to note is that appear double integrals which involve sine and cosine functions in the same integral. Let I cs be the double integral in the order cosine and sine functions and I sc be the integral for the opposite order. In case when the plasma term M pl dominates M φ in G 1 , is possible to find analytic expressions for I cs and I sc which are respectively given by Now putting expressions for I cs and I sc in the first term in V (T ), we get Using (6.5) and (6.4) together with expression for y, we get the following expression for the degree of circular polarization |V |/I γ at T = T 0 to second order in perturbation theory An important thing to note is that photon intensity must decrease and never increase in the case of pseudoscalar particle production. This happens because in our case we are not assuming photon injection in the medium by some external source that would eventually increase the photon intensity. Since I γ must decrease or at least remain constant, this implies that the quantity −y + y Q i in (6.4) must be negative. Indeed, this is true and can be easily verified by evaluating the second and third terms in (6.4) for a given frequency. All told, is a necessary condition but not sufficient. We must also require that 1−y +yQ i > 0 since the intensity is a positive quantity. If we consider for example the working frequency of MIPOL experiment, ν 0 = 33 GHz, for 1 − y + yQ i > 0 to be satisfied 21 we must have |g φγ | < 1.66 × 10 −15 (G/B e0 ) GeV −1 . Another consideration to be made is related with the parameter V in (6.5). As we can see in (6.5) the first term on the r. h. s. corresponds to the QED and CM effects while the second and third corresponds essentially to mixed terms. However, since perturbative expansion to second order does not reveal the asymptotic behaviour of the series and because we expect that pesudoscalar contribution to V to be small or (T i /T 0 ) 3 V (T 0 ) U i , we require the additional condition that the dominant term 211.38ã −1 b 1. In this case for ν 0 = 33 GHz we get |g φγ | < 5.13 × 10 −20 (G/B e0 ) GeV −1 . After some algebraic operations in (6.6) and requiring that P C (T 0 ) < 7 × 10 −5 (MIPOL upper limit), we get the following constraint on g φγ from upper limit on the degree of circular polarization which is within the constraint (6.2) and we took U i −Q i with Q i 10 −6 . The MIPOL upper limit (6.7) is also satisfied if |g φγ | 5.13 × 10 −20 (G/B e0 ) GeV −1 for 211.38ã −1 b 1, which is a more conservative limit and satisfies both MIPOL upper limit and the constraint of perturbation theory.
In the domain of circular polarization, now it remains to study the last case when the term M φ dominates the plasma term M pl in ∆M 1,2 . Proceeding in the same way as above, we calculate the functions G 1,2 which enter the trigonometric functions in I γ and V . An important difference now is that we do not have to worry about the ionization fraction since it does not appear in M φ . Inserting all necessary quantities into G 1,2 we get The next step is to calculate the double integrals I cs and I sc . However, in this case there are no known analytic solutions for both type of integrals so we must evaluate them numerically together with y for some specific values of the parameters. It would be more convenient first to write y = b f 1 (m φ , ν 0 ) and I cs −I sc = bf 2 (m φ , ν 0 ) and after calculate f 1 and f 2 numerically for given values of m φ and ν 0 . Let us recall that now we are in the situation in which the constraints of perturbation theory are given by (6.1). We may consider for example m φ = 10 −8 eV and ν 0 = 33 GHz which corresponds to the working frequency of MIPOL experiment. Using numerical integration we obtain f 1 = 9.37 × 10 −24 and f 2 = 6.74 × 10 −11 . Second, using the relation P C (T 0 ) = | − (G(T 0 ) + bf 2 )U i |/(1 − y(1 − Q i )) we get the following constraint |g φγ | < 1.26 × 10 −16 (G/B e0 ) GeV −1 . (6.8) If we consider for example m φ = 10 −6 eV for the same working frequency we would obtain f 1 = 9.37 × 10 −32 , f 2 = 6.74 × 10 −15 and the above limit would be two orders of magnitude weaker.

Weak CMB polarization at decoupling
So far, in our treatment of generation of CMB polarization, we have assumed a priori that the CMB acquired a small polarization due to Thomson scattering at decoupling time, with non zero Stokes parameters Q i and U i . However, even though this assumption may seems reasonable, it is not accurate because of the fact that measurements of the CMB properties are done at present epoch and not at decoupling. In fact, if cosmological magnetic fields were present at decoupling epoch, production of pseudoscalar particles would generate CMB polarization independently on Thomson scattering and consequently there is no unambigous way to tell, if the linear polarization experimentally observed at present is due to either Thomson scattering or photon-pseudoscalar particle mixing 22 .
Another important assumption about generation of CMB polarization, is essentially generated at decoupling time and it remains invariant during subsequent evolution of the universe. Even this assumption is not completely end of the story, since for example the CMB may acquire a very small additional polarization during the reionization epoch due to Thomson scattering. So based on these two basic assumptions one would have that degree of linear polarization is generated by Thomson scattering and it remains almost invariant after the decoupling epoch.
In the previous sections we found upper limits/constraints on g φγ from current limit on the degree of circular polarization which is not directly generated by Thomson scattering. However, one may note that these upper limits/constraints found for g φγ in case of circular polarization, we have to high accuracy P L (T 0 ) P L (T i ). These considerations would suggest that based on current limit on the degree of circular polarization and consequently on derived limits/constraints on g φγ from it, in principle the observed CMB linear polarization could be generated by a combination of Thomson scattering and magnetic-optic effects such as photon-pseudoscalar particle mixing where the latter is the subdominant component.
In this section, we consider another possible situation, the other way round, where the contribution of Thomson scattering to the linear polarization is supposed to be very small, so it can be completely neglected and generation of linear polarization at post decoupling time (or large angular scales polarization) is mostly due to photonpseudoscalar particle mixing. Obviously this assumption does not mean that there is no generation of polarization due to Thomson scattering at decoupling time. Based on this assumption, we would have that the CMB at decoupling is almost unpolarized with Q i 0, U i 0, V i 0. Let us concentrate for the moment on the degree of linear polarization which is given by P L (T ) = (Q 2 (T ) + U 2 (T )) 1/2 /I γ (T ). In the weak mixing case, from (5.12) and (5.9) and considering for example the case when |M pl | > |M φ |, we get the following expression for P L (T ) at present P L (T 0 ) = y 1 − y (weak mixing and |M pl | > |M φ |), (6.9) where we assumed the CMB unpolarized at decoupling and U (T ) = 0 for unpolarized light. Since P L depends on y = 4ã −2 b(1−cos[52.84ã]), it explicitly depends on the photon frequency ν 0 . Even though one can easily calculate P L at a given frequency, is more convenient to calculate its average value on a given interval. If we consider for example 10 8 Hz ≤ ν 0 ≤ 10 11 Hz, the average value ofã −2 (1 − cos[52.84ã]) 2.35 × 10 −15 . Considering that the degree of linear polarization of CMB at present is P L (T 0 ) 10 −6 , we get the following value for |g φγ | 1.28 × 10 −18 (G/B e0 ).
We may note that in case when the argument of sine function in (6.10), F φγ (T 0 ) is less than unity (or g φγ B e0 < 1.17 × 10 −24 G GeV −1 ) and because in general P L 1, from (6.10) one would get sin 2 [F φγ (T 0 )] 2P L (T 0 ). Considering that P L (T 0 ) 10 −6 , we get in the resonant case the following frequency independent value |g φγ | 1.66 × 10 −27 (G/B e0 ) GeV −1 if g φγ B e0 < 1.17 × 10 −24 G GeV −1 (6.14) The solution (6.14) corresponds to the last set in (6.13) for n = 0. Another important thing to note in the case of linear polarization is U (T ) = 0 for initially unpolarized CMB and consequently the angle of the polarization ellipse is tan[2ψ(T )] = 0, which implies a horizontal linear polarization and no rotation of the polarization plane 24 . In case of circular polarization, one may observe from (5.8) and (5.9) that V (T ) = 0 for initially unpolarized CMB, which means no generation of circular polarization for transverse magnetic field. Mixing of CMB photons with pseudoscalar particles would also generate secondary CMB temperature anisotropy from an almost initially thermalized state 25 . In order to prove this statement, consider the CMB at decoupling completely in almost thermal equilibrium and consequently almost unpolarized, Q i 0, U i 0, V i 0. Indeed, this assumption is well motivated since |Q i |, |U i |, |V i | I i in both weak and strong mixing regime, see Eqs. (5.11)-(5.12). Consider the evolution of intensity I γ as a function of T for two specific observation directions: parallel and perpendicular to B e . In the direction parallel to B e , as we already have seen in Sec. 4, it is induced only the Faraday effect and other magneto-optic effects are absent. As we saw there, the intensity for parallel propagation changes only due universe expansion, (T i /T ) 3 I γ (T ) = I γ (T i ). The intensity observed parallel to B e is the same as that observed without the presence of the external field or unperturbed universe. On the other hand, the observed intensity for perpendicular propagation with respect to B e , for initially unpolarized CMB, is given by (5.11) which in case of resonant mixing is ( In case of thermalized CMB at decoupling we have that To linear order in δT , we find the following relation from the black body intensity where we defined x ≡ 2πν 0 /T 0 and T 0 is the average value of the temperature at present, averaged over all directions in the sky. In case when x < 1 or ν 0 < 5.63 × 10 10 Hz, namely Rayleigh-Jeans regime, we have essentially δI γ /I γ δT /T 0 , while for x > 1 (Wien regime) one must use the whole expression (6.15). We can use (6.15) to find the value of g φγ in the resonant case 26 . Therefore, we have The temperature anisotropy of the CMB depends on the angular separation between two points in the sky. Since we are comparing the intensity between parallel and perpendicular observations, this means an angular separation scale of 90 • . According to WMPA9 collaboration [44], the temperature anisotropy at 90 • or multipole moment l = 2, is δT /T 3 × 10 −5 . From (6.15) and (6.16) we get the following value of g φγ |g φγ | 9.12 × 10 −27 x e x e x − 1 where we considered for simplicity only the case when when F φγ (T 0 ) < 1 and assumed that the contribution of large scale magnetic field to temperature anisotropy is subdominant to photon-pseudoscalar mixing contribution. In general (6.16) has multiple solutions similar to (6.13) when F (T 0 ) > 1. In the Rayleigh-Jeans part of the spectrum we get |g φγ | 9.12 × 10 −27 (G/B e0 ) GeV −1 .
24 This conclusion applies for unpolarized CMB at decoupling and for transverse magnetic field. 25 In addition to generation of CMB temperature anisotropy by photon-pseudoscalar mixing, there is also generation of temperature anisotropy by the large scale magnetic field itself. Consequently the total temperature anisotropy is given by the sum of photonpseudoscalar mixing and magnetic field temperature anisotropy contributions. 26 Here we consider for simplicity only the resonant case, however expression (6.16) is also valid in the strong mixing case. In general given the value of temperature anisotropy, the trigonometric Eq. (6.16) has multiple solutions in both strong and resonant mixing regimes.
In the weak mixing case and for |M pl | > |M φ | we get (6.17) where we used the expression for I γ in (5.12) with Q i = 0 and used the definition 27 y for |M pl | > |M φ |. As we can see from (6.17), δI γ is proportional to the fast varying function, 1 − cos[52.84ã], which forã = 2nπ/52.84 is equal to zero, with n ∈ Z, independently on the value of g φγ . In case when |M pl | < |M φ |, the value of y = bf 1 (m φ , ν 0 ) can be calculated numerically as we did for the case of circular polarization for given values of m φ and ν 0 .

Discussion and conclusions
In this work we have studied the most important magneto-optic effects and their impact in the generation of CMB polarization. We presented a systematic study of each of them where we mostly focused on the generation of CMB circular polarization. In this work we found the equations of motion for photon and pseudoscalar fields in an external magnetic field in the WKB approximation, and then found the equations of motion for the Stokes parameters by using density matrix approach as shown in Sec. 3. The resulting equations describe the mixing of different magneto-optic effects which obviously complicate the situation but on the other hand give richer scenarios. In Sec. 4 we studied the vacuum polarization and CM effects separately, in order to isolate the contribution of each of them to CMB polarization. They are second order magneto-optic effects on magnetic field amplitude B e and are responsible for generation of phase shifts between the states A + and A × . These effects generate CMB elliptical polarization only in the case when the CMB is initially polarized. In this work we concentrated in the post decoupling epoch and worked under the hypothesis that the CMB acquired a small polarization at decoupling time due to Thomson scattering. We used perturbation theory and found the evolution as a function of T of the Stokes parameters. We studied in particular the generation of circular polarization which is represented by the Stokes parameter V , in cases of observation angles Φ = π/2 and Φ = π/2.
The contribution of vacuum polarization and CM effects to V depends essentially on Φ, B e0 , ν 0 and on the magnitude of the Stokes parameters at decoupling which, on the other hand, depend on the temperature anisotropy. In this work we assumed that V i = 0 at decoupling while the other parameters are non zero. The magnitude of the parameters Q i and U i obviously are smaller than temperature anisotropy and observations of CMB linear polarization give an order of magnitude of Q i ∼ U i ∼ 10 −6 .
In the case of vacuum polarization and Φ = π/2, the degree of circular polarization is proportional to Q i and U i , as shown in (4.20) and in most cases is the term proportional to U i which dominates. This term on the other hand is proportional to ν 0 and B 2 e0 . Consequently, significant generation of circular polarization would occur in the high frequency part of the CMB and for higher values of B e0 . In this work we used in our estimates a canonical value of B e0 ∼ nG but in principle higher values are possible. If for example one observes the CMB in the Wien region, say at ν 0 700 GHz and the magnetic field is of the order of 100 nG, the degree of circular polarization would be of the order P C ∼ 10 −11 while for B e0 ∼ nG is four orders of magnitude smaller.
Also for the CM effect, the degree of circular polarization is proportional to the initial values of Stokes parameters at decoupling and to B e0 , ν 0 and Φ. One distinguishing feature of the CM effect is the relation between V 0 and ν 0 which is V 0 ∝ ν −3 0 . This relation makes the CM effect quite appealing in regard to generation of circular polarization in the Rayleigh-Jeans part of the spectrum. For Φ = π/2 the degree of circular polarization is given in (4.33) where the first term is proportional to Q i and the second term is proportional to U i . The term proportional to Q i shares a common feature with the vacuum polarization by the fact it gets contribution from the Faraday effect which is encoded in ρ. Under the approximations used in Sec. 4, the term proportional to Q i is smaller than that proportional to U i . The latter coincides with the solution found for V in case when Φ = π/2 for G(T ) 1 and F (T ) < 1. This means that the contribution of the CM effect to circular polarization is bigger in the limit Φ → π/2. The degree of circular polarization is substantive in the frequency region ν 0 ∼ 10 8 − 10 9 Hz while for higher frequencies ν 0 ∼ 10 11 Hz, the contribution of CM effect to V 0 is subdominant to vacuum polarization. For example, if ν 0 ∼ 10 8 Hz and B e0 ∼ nG, the degree of circular polarization for the CM effect would be P C ∼ 10 −10 while if B e0 ∼ 100 nG, P C ∼ 10 −6 .
In this work, we also studied the generation of elliptic polarization due to photon-pseudoscalar particle mixing in cosmic magnetic field, with emphasis on the degree of circular polarization. Differently from the vacuum polarization and CM effects, photon-pseudoscalar mixing has in addition two more independent parameters which are m φ and g φγ . We studied this mechanism in case of only transverse magnetic field and used perturbation theory to find the evolution in T of the Stokes parameters. We used perturbation theory in two mixing regimes, namely weak and strong mixing and estimated the degree of circular polarization at present epoch.
Since the parameters g φγ and m φ are free and in general span a wide range of values, we used the present upper limit on the degree of circular polarization in order to constrain g φγ and m φ . These parameters on the other hand are constrained by the mixing regimes, therefore the limits that we presented are valid in these regimes. In the strong mixing regime, in general one has to solve trigonometric equations or inequations which have as independent variable g φγ B e0 . The solutions generally, depend on an integer number n and consequently they are not unique. The interval of values of g φγ B e0 can in principle be narrowed by complementary constraints on g φγ from other methods. On the other hand, in the weak mixing case there is not such a dependence on n. In this case by using the upper limit on P C obtained from MIPOL experiment, we got the constraint |g φγ | < 4.29 × 10 −19 (G/B e0 ) GeV −1 for m φ < 1.6 × 10 −14 eV.
Other limits have been obtained late time generation of degree of linear polarization and by considering the case of almost unpolarized CMB at decoupling. In the weak mixing case, we obtained the average value over frequency of |g φγ | ∼ 10 −18 (G/B e0 ) for m φ < 1.6 × 10 −14 eV and P L (T 0 ) 10 −6 . In the strong mixing case, again one obtains values of g φγ B e0 that depends on n and therefore there is no unique solution. The same thing happens even in the resonant case with the particular case that, if, g φγ B e0 < 1.17 × 10 −24 , then from P L 10 −6 we get the value |g φγ | 1.66 × 10 −27 (G/B e0 ) for m φ 1.6 × 10 −14 eV. As in the case of vacuum polarization and CM effects, photon-pseudoscalar particle mixing generates non uniform degree of circular polarization and rotation of the polarization plane across the sky. This fact can be used in order to understand if the observed linear polarization at present has a non uniform component across the sky. If this would be true, it might be due to photon-pseudoscalar particle mixing if it is the dominant mechanism of generation of linear polarization among magneto-optic effects studied.
From the experimental side, it turns out that among CM and QED effects, the CM effect is the most promising effect on generating circular polarization in the low frequency part of the CMB due to the dependence V 0 ∝ ν −3 0 , while the vacuum polarization is the dominant one in the high frequency part due to V 0 ∝ ν 0 . The degree of circular polarization due to the CM effect in the low frequency part, in general, is bigger than that generated by vacuum polarization at high frequencies. Moreover, vacuum polarization and CM effects generate a rotation of the polarization plane of the CMB and this rotation together with the degree of circular polarization are not uniform across the sky because they depend on the observation angle Φ. These facts would suggest that observation of CMB circular polarization is more likely to happen in the low frequency part of the CMB, mostly due to the CM effect and if, the observation frequency range is not a big detection issue. On the other hand, if one is interested in the measurement of the rotation angle of the polarization plane, the non uniformity of the rotation across the sky might be an issue.
In order to detect CMB circular polarization, probably the most convenient frequency range would be for ν 0 ∼ 10 8 − 10 9 Hz where the degree of circular polarization would be in the interval P C 10 −13 − 10 −10 for B e0 1 nG due to CM effect where the higher value of P C corresponds to the lower value of ν 0 . In this frequency range the contribution of vacuum polarization is completely negligible with respect to CM effect. The vacuum polarization is dominant to the CM effect in the Wien regime and for ν 0 ∼ 700 GHz we found the interval P C 10 −15 − 10 −11 for the interval B e0 1 − 100 nG. Assuming for the moment that frequency observation range is not an issue and it can be fixed based on experiment characteristic, one main problem on the detection of CMB circular polarization is related to the magnetic field amplitude which is poorly known. In the case of intergalactic magnetic field, usually upper limits are found by different methods and its amplitude is expected to be less than 1 nG (canonical value) up to a value of less than 100 nG. Using a canonical value of B e0 1 nG, one would expect that the degree of circular polarization to be P C 10 −10 in the Rayleigh-Jeans regime and P C 10 −15 in the Wien region. These values of the degree of circular polarization corresponds essentially to the case when Φ = π/2 and for different values of Φ, P C is usually smaller and not uniform across the sky.
The photon-pseudoscalar mixing contributes to circular polarization as well and based on values of g φγ and m φ found from current upper limit on circular polarization from MIPOL experiment, its contribution might be bigger than CM and vacuum polarization effects. So, let us assume for example that B e0 1 nG and circular polarization would be detected with degree of circular polarization with value in the range, 10 −10 P C 10 −6 . This would mean there is a contribution to P C due to photon-pseudoscalar mixing which is much bigger than other magneto-optic effects or the amplitude of the magnetic field might be higher than assumed or circular polarization is generated by another effect not considered in this work. If one would detect CMB circular polarization with average value of P C 10 −10 , it is more likely that CM and vacuum polarization effects are the source of this polarization.
In all studied effects, we have assumed that the large scale magnetic field was present at the decoupling epoch therefore the field has been assumed to have primordial origin and a function of spacetime coordinates, B(x, t), namely the field is non homogeneous in space and time. In the case of vacuum polarization, the expressions for the photon polarization tensor and derived quantities such as the index of refraction have been derived under the assumption that the electromagnetic field tensor satisfies the condition |∂ µ F σρ | m e |F σρ |, see Ref. [45] for details. This condition on the electromagnetic field tensor translates into conditions on the magnetic field |∂B i e (x, t)/∂t| m e |B i e (x, t)| and |∂B i e (x, t)/∂x| m e |B i e (x, t)|. Obviously both the last conditions on B i e (x, t) are satisfied in an expanding universe for a large class of magnetic fields where the former can be written as H −1 (t) 2/m e = 7.74 × 10 −11 cm for the Hubble radius as a function of time and the latter condition can be written as l B (t) 2/m e = 7.74 × 10 −11 cm where l B is the variation scale in space of the external magnetic field. In the case of the CM effect, the elements of photon polarization tensor are usually derived for constant magnetic fields but since this effect is similar to the QED effect, namely is of the second order in B e , one can extend the results for constant fields also to the case when |∂ µ F σρ | m e |F σρ | in complete analogy with the vacuum polarization effect. In the case of photon-pseudoscalar mixing, the magnetic field can be either homogeneous or non homogeneous as far as the photon wavelength λ l B (t) in the WKB approximation, namely the magnetic field is a slowly varying function in space and time with respect to the photon wavelength or frequency.
It is worth to mention also what has not been studied in this work. The first thing is related to Thomson scattering and scattering of pseudoscalar particles at post decoupling epoch, namely for T < 2970 K and their absence in our density matrix formalism. In general, scattering is a mechanism of coherence breaking for mixing/oscillation processes which results in damping of the fields. In the density matrix formalism, the structure of the damping operator can be calculated by using field theory for scattering which is essentially the calculation of the commutator [H T , ρ] on the r. h. s. of (3.3) where H T includes the Hamiltonian for the Thomson scattering and that of scattering of pseudoscalar particles. However, quite often the damping term due to scattering, in case of non degenerate and non relativistic electron gas can be approximated 28 by, −i{Γ, ρ}, where Γ is the scattering rate matrix of photons and pseudoscalar particles which is diagonal in the basis |A + , |A × , |φ . Consequently, the damping term due to scattering, would have the same structure as the damping term due to Hubble friction. Therefore, the Stokes parameters would be affected by scattering but not their ratio because it cancels out exactly as the damping term due to Hubble friction. However, one must always keep in mind that this is an approximation.
The second thing is related to the case Φ = π/2 for the photon-pseudoscalar particle mixing. In Sec. 5, we found the equations of motion for the reduced Stokes vectors in case of transverse external magnetic field. In this case it was possible to find two sets of decoupled differential equations for the reduced Stokes vectors and solved the equations by using perturbation theory. If the field is not transverse, namely Φ = π/2, in general one has to solve simultaneously, a system of nine linear differential equations of the first order which can be problematic to solve even numerically because quite often they are stiff. We shall treat this problem in more details elsewhere but even at this stage we can outline very important conclusions about the nature of the solutions and the impact on the CMB polarization.
The importance of solutions of the equations of motion in the case Φ = π/2 (for photon-pseudoscalar mixing) relies in the fact, that being the system of equations linear, see Eqs. (5.5), the solutions will be proportional to initial values at a given temperature T i which does necessarily coincides with decoupling temperature. Therefore, each Stokes parameter would be proportional to I γ (T i ), Q i , U i , V i , S 4i etc., and for T < T i usual Stokes parameters (those which in general interest us) would be different from zero even in case of initially unpolarized CMB at T = T i . We saw similar situation in Sec. 6.2, where we studied the case of unpolarized CMB at decoupling for transverse magnetic field. Consequently, the CMB would acquire polarization independently on Thomson scattering, even in case when it is initially unpolarized.
This situation would be very important in order to investigate prior decoupling CMB polarization due to photonpseudoscalar mixing in external magnetic field. According to standard cosmology, generation of CMB polarization occurs at or very close to decoupling time due to Thomson scattering when the condition of tight coupling between photons and electro-baryon plasma breaks down. Indeed, for most models of generation of CMB polarization which include scalar perturbations, magnetic fields, gravitational waves etc., at the end is always the Thomson scattering which generates CMB polarization [47]. The tight coupling condition would imply that, if there is any degree of polarization prior to decoupling, generated at temperature T , it would be damped very fast due to scattering of photons with electrons and baryons. However, as we have seen and discussed in this work, photon-pseudoscalar mixing apart from generating temperature anisotropy as shown in Sec. 6.2 and spectral distortions of the CMB [48], it generates also polarization, independently on Thomson scattering. Consequently, here we advance the hypothesis that photon-pseudoscalar mixing might generate non uniform CMB polarization across the sky, even before decoupling epoch, if the rate of photon-pseudoscalar oscillation is faster than photon scattering rate with electro-baryon plasma. Obviously, all said about this hypothesis would depend on pseudoscalar field parameters m φ and g φγ . The suggested hypothesis needs further attentive study and it would be too premature to conclude that it is indeed the case.
AKNOWLEDGMENTS: This work is supported by the Russian Science Foundation Grant Nr. 16-12-10037. I would like to thank LNGS for the support received through the fellowship POR 2007-2013 'Sapere e Crescita' where part of this research was conducted.

A Photon density operator and Stokes parameters
In most cases which are of interest in physics one has to deal with quantities that are proportional to the amplitude square of fields rather than the amplitude itself and that are connected with photon polarization. Such quantities are the Stokes parameters which give a complete description of photon polarization state. Use of Stokes parameters to describe photon polarization is very convenient because allow us to deduct important information about photon polarization by using four measurable quantities associated with the photon field. The first quantity or observable expresses the intensity of the photon field while the remaining three quantities completely describe its polarization state. Stokes parameters can be applied to unpolarized, partially polarized and completely polarized light and have the mathematical convenience of not being expressed in terms of the photon amplitude, which is in general not observable. However, an observable quantity is the photon field intensity which is derived by taking the time average of the square of the amplitude.
Consider a plane wave (not necessarily monochromatic) propagating along the z direction in a given cartesian coordinate system and consider the wave at z = 0. The wave electric field vector can be decomposed along the x and y components as follows E x (t) = E x0 (t) cos[ωt + δ x ], E y (t) = E y0 (t) cos[ωt + δ y ], where δ x , δ y are respectively the instantaneous wave phases for each field component, E x0 and E y0 are respectively the instantaneous wave amplitudes and ω is the instantaneous wave angular frequency. Here we consider the hypothesis that electric fields amplitudes E x0 , E y0 and field phases δ x , δ y slowly fluctuate in time in comparison with the rapid vibration of cosine functions. In case of nearly monochromatic wave, the Stoke's parameters are defined as follows The quantities Q and U in general depends on the orientation of the coordinate system used for measurements but the quantities Q 2 + U 2 , I and V are invariant under such orientation. While the quantities Q 2 + U 2 , I and V are invariant under coordinate system rotation, Q 2 + U 2 and V are not necessarily invariant under simultaneous coordinate system rotation and non constant phase generation δ(t). It is straightforward to show that the orientation angle ψ of the polarization ellipse can be expressed in terms of Q and U as tan(2ψ) = U/Q, where 0 ≤ ψ < π physically represents the polar angle of the polarization ellipse. One can also express the total degree of polarization of the wave in terms of the Stokes parameters as follows where 0 ≤ P ≤ 1. If P = 1 the wave is completely polarized, if P = 0 the wave is completely unpolarized (natural light, Q = U = V = 0) and if P < 1 the wave is partially polarized. The degree of linear polarization is given by P L = Q 2 + U 2 /I and the degree of circular polarization is given by P C = |V |/I. The description of polarized light in quantum optics is different from classical optics and in general their connection is not obvious. In classical optics, polarization is described in terms of amplitudes and polarization ellipse while in quantum optics it is described in terms of the density matrix. Having defined the Stoke's parameters in (A.1), it is desirable to connect them with the polarization density matrix of a quantum optical system. As shown in Ref. [49], Stokes parameters are very important tool for treating polarization problems in both quantum and classical optical systems. In order to outline their connection with the polarization density matrix, let |A be an arbitrary photon state which is a linear superposition of the quantum polarization states |A + and |A × |A = c 1 |A + + c 2 |A × , where c 1 , c 2 are complex amplitudes. Their absolute value square represents the probability to find a photon respectively in the state |A + or |A × . In both classical and quantum optics, the polarization state of the wave is completely described in terms of complex amplitudes c 1 , c 2 and one can define the elements of the density matrix ρ as ρ ij = c * i c j , (i, j = 1, 2). If F is any observable of the system, its expectation value on an arbitrary state is given by F = Tr(F ij ρ ij ) where summation over repeated indices is used.
Following Ref. [49], one can associate to the Stokes parameters their corresponding quantum mechanical operators asÎ where the expectation value of each operator is given by One can find after some trivial algebra the polarization density matrix in the basis spanned by the photon states |A + , |A × ρ = 1 2 is an important representation of the density matrix in terms of the Stokes parameters and is very useful in many contexts. It is important to note that representation (A.3) is in the basis of the vector potential states and not in the electric field basis which is the most common used case. Consequently, the physical dimensions of the density matrix (and also Stokes parameters) in the vector potential basis are different from those in the electric field basis.