Non-comoving Cosmology
J. A. R. Cembranos A. L. Maroto H. Villarrubia-Rojo




Departamento de  Física Teórica and Instituto de Física de Partículas y del Cosmos IPARCOS, Universidad Complutense de Madrid, E-28040 Madrid, Spain



    One of the fundamental assumptions of the standard ΛCDM cosmology is that,  
    on large scales,  all the matter-energy components of the Universe share a common rest frame. 
    This seems natural for the visible sector, that has been in thermal contact and tightly coupled in 
    the primeval Universe. The dark sector, on the other hand, does not have any non-gravitational 
    interaction known to date and therefore, there is no a priori reason to impose that it is comoving 
    with ordinary matter. In this work we explore the consequences of relaxing this assumption and study 
    the cosmology of non-comoving fluids. We show that it is possible to construct a homogeneous and 
    isotropic cosmology with a collection of fluids moving with non-relativistic velocities. Our model 
    extends ΛCDM with the addition of a single free parameter β_0, the initial velocity of 
    the visible sector with respect to the frame that observes a homogeneous and isotropic universe. This 
    modification gives rise to a rich phenomenology, while being consistent with current observations for 
    β_0<1.6× 10^-3 (95% CL). This work establishes the general framework to describe a 
    non-comoving cosmology and extracts its first observational consequences for large-scale structure. 
    Among the observable effects, we find sizeable modifications in the density-velocity and 
    density-lensing potential cross-correlation spectra. These corrections give rise to deviations from
    statistical isotropy with a dipolar structure. The relative motion between the different 
    fluids also couples the vector and scalar modes, the latter acting as sources for  metric vector 
    modes and vorticity for all the species. 


Non-comoving Cosmology
    J. A. R. Cembranos[cembra@ucm.es], 
        A. L. Maroto[maroto@ucm.es],
        and H. Villarrubia-Rojo[hectorvi@ucm.es]
    December 30, 2023
=============================================================================================================================





§ INTRODUCTION

    The isotropy and homogeneity of the Universe on large scales are two foundational assumptions
    of the standard cosmological model, the so-called ΛCDM model. These two assumptions are usually grouped under the
    name of Cosmological Principle. All the observational evidence, ranging from the extremely isotropic
    cosmic microwave background (CMB) <cit.> to the galaxy number counts and the measured expansion
    from SNIa <cit.>, supports the conclusion that 
    the Universe is very nearly isotropic on large scales.
    However, the notions of homogeneity and isotropy are inextricably linked with the election of a 
    privileged frame. For any observer moving with respect to this frame, the Universe would appear 
    anisotropic and inhomogeneous. This is precisely our situation on Earth. 

    
    Starting with the early CMB measurements <cit.>, a significant dipole 
    modulation, much larger than any other anisotropy, was found. This was readily interpreted as a 
    kinematical effect: a Doppler shifting effect arising from the relative motion of the Earth with 
    respect to the CMB rest frame, i.e. a frame in which the CMB looks isotropic. Recent analysis by 
    the Planck Collaboration <cit.>
    explored other kinematical effects, like the violation of statistical isotropy induced by the observer
    motion, and reported an independent measurement of our relative velocity with respect to the CMB
    frame. This measured velocity can, given the uncertainties, fully account for the observed dipole, 
    supporting its kinematical origin. Even if it is mostly kinematical, it may still contain an 
    intrinsic contribution. Some authors have proposed searches for the intrinsic dipole, e.g. using
    spectral distortions <cit.>.

    
    A different kind of dipole should appear in the distribution of galaxies, induced by our motion with
    respect to the matter frame, i.e. a frame in which the matter distribution looks isotropic. The 
    origin of the large scale structure (LSS) dipole lies in a combination of Doppler shifting and
    aberration effects in the galaxy number counts <cit.>. Unfortunately, current 
    observations can only loosely constrain its amplitude and direction, yielding a value compatible 
    with the CMB dipole <cit.>. Future surveys like Euclid 
    <cit.> and SKA <cit.> will measure it with unprecedented accuracy.

    
    To complete the picture, we only need to know the relative velocity between the matter and CMB frames.
    Concerning this point, ΛCDM contains the underlying assumption, that usually goes by 
    unnoticed, that both frames coincide. ΛCDM assumes matter and CMB to be comoving.
    As we will see in this work, it is possible to relax this condition. The homogeneous and isotropic
    Robertson-Walker (RW) metric can be sourced, at the background level, using non-comoving fluids. 
    Thus, we will show that it is possible to construct a viable cosmological model for non-comoving 
    fluids, with interesting phenomenological consequences and without any flagrant isotropy violation.

    
    Early theoretical work concerning non-comoving fluids was mostly developed under the framework of
    tilted universes. The term was coined in the groundbreaking work by King and Ellis <cit.>.
    The authors considered a class of homogeneous models sourced by a single moving fluid, i.e. models
    in which the fluid 4-velocity is tilted with respect to the homogeneous hypersurfaces. These tilted
    models produce homogeneous but anisotropic universes. In a different context, Coley and Tupper
    <cit.> analyzed two-fluid cosmological models with general imperfect and non-comoving
    fluids. In order to source homogeneous and isotropic RW metrics, only very special configurations
    with radial velocities were considered. Later on, Turner <cit.> proposed a theoretical 
    mechanism to produce a mismatch between matter and CMB velocities. In Turner's tilted universes,
    the presence of a near-horizon-sized perturbation, remnant of inflation, could introduce a spatial
    gradient, driving the velocity of matter.

    
    The analysis of non-comoving fluids has been extended to dark energy. Given our fundamental ignorance
    about the behaviour of the dark sector, it is conceivable that it has not ever been coupled to 
    ordinary matter and that it does not share the same rest frame. Following this idea, a model of moving 
    dark energy was proposed in <cit.>. In this case, for a dynamical dark energy fluid, 
    even if the matter and CMB frames coincide initially, they differ at late times. Different models of 
    moving homogeneous dark energy were analyzed in <cit.>, as well as its possible 
    impact in observables like the CMB quadrupole. The construction of a fully anisotropic model in which 
    the full dark sector, i.e dark energy and dark matter, is non-comoving with the CMB and ordinary matter 
    was carried out in <cit.>. The authors analyzed a Bianchi I universe in which dark 
    matter and dark energy had different relative velocities with respect to the frame of ordinary 
    matter and then derived some observables, like a modified luminosity-distance relation and CMB 
    quadrupole.

    
    From the observational point of view, a signal of the relative motion of the matter and CMB frames 
    will be the detection of a large-scale bulk flow. In recent years, several works have claimed 
    measurements of matter flows well in excess the ΛCDM predictions 
	on different scales and at different statistical confidence levels <cit.>. 
	Although there seems to be a broad agreement on the direction of the flow, the amplitude is
	still  subject to  controversy <cit.>. Such flows will be an indication of the existence 
	of a cosmological preferred spatial direction. On the other hand, detected anomalies in the low 
	multipoles of the CMB  temperature power spectrum <cit.>, such as the low-multipole 
	allignment and the dipolar or hemispherical anomalies, also suggests the presence of a preferred 
	cosmological direction <cit.>. This fact has triggered the search for mechanisms 
	which could break isotropy while keeping the 	predictions of the standard cosmology.
    
    
    This work builds upon these previous studies, but we will present the first complete analysis for the
    evolution of a set of non-comoving fluids, from the early to the late Universe, both at the
    background and perturbation level. As we will see later, it is reasonable to assume that any pair
    of tightly coupled fluids share the same velocity. Hence, we can expect that photons, baryons
    and neutrinos, being in thermal contact in the early Universe, shared a common rest frame, i.e. a 
    frame in which the plasma looked isotropic. However, there is no a priori reason to assume 
    the same about the dark sector. The dark sector, regardless of its composition, may very well possess 
    its own rest frame, with a given global velocity with respect to the visible sector. The only 
    reasonable assumption is that there is one frame that observes a homogeneous and isotropic 
    universe, i.e. a RW background. In this paper we take seriously this possibility and prove, at the 
    background level, that
    
        
  * If the dark and visible sectors are initially moving with non-relativistic 
            velocities, to first order in these velocities, it is possible to define a cosmic center 
            of mass frame that observes a RW background. Our model only contains one additional 
            free parameter β_0, the velocity of the visible sector with respect to the cosmic center 
            of mass frame deep inside the radiation era. As we will see later, 
            β_0 < 1.6× 10^-3 is a conservative limit to be in agreement with all observations.
        
  * The subsequent evolution gives rise to relative velocities between all the different 
            components of the visible sector, e.g. between baryons and photons, defining each of them its
            own rest frame.
    
    While at the perturbation level,
    
        
  * There appear couplings, of order β_0, between scalar and vector modes.
        
  * There is a production of vorticity for all the species and a net production of vorticity
            and metric vector modes.
        
  * The transfer function of every cosmological quantity acquires a dipolar contribution
            of order β_0. 
        
  * There is a violation of statistical isotropy of order β_0. Among the different
            observables where such an effect could be measured, the easiest to compute is the
            cross-correlation spectrum between different scalar perturbations, which acquire a dipolar
            contribution of order β_0.
    
    
    In this work we will limit ourselves to LSS observables, letting for a forthcoming work the analysis
    of CMB signatures <cit.>.
    The structure of the paper goes as follows. Section <ref> assesses the problem of
    constructing a homogeneous and isotropic universe using non-comoving fluids. First, only the 
    energy-momentum tensor for a perfect fluid with bulk velocity is considered in <ref>,
    where the conditions for isotropy and homogeneity are discussed. Then, the formalism is extended
    to imperfect fluids in <ref>, using the kinetic approach, where both background and
    perturbations are studied in <ref> and <ref>. Section 
    <ref> analyzes the dynamics of this model from the point of view of the
    Boltzmann equations. The free-streaming term is derived in <ref>, while the
    collision term for the photon-baryon plasma is computed in <ref>. The
    main evolution equations are presented in <ref>. The usual scalar-vector-tensor 
    decomposition is performed in section <ref>, where we describe our approximation
    scheme. Section <ref> contains the Einstein equations, which are not modified in our case.
    In section <ref>, we present a reduced version of the original system, under our
    approximation scheme. Sections <ref>, <ref> and
    <ref> are devoted to bulk velocities, scalar and vector modes, respectively.
    Finally, section <ref> contains the numerical solution of the aforementioned
    reduced system and discusses some observables.
    Section <ref> gathers the main conclusions and presents some prospects for future
    work.



§ MOVING FLUIDS

    

 §.§ Perfect fluid with bulk velocity

        

  §.§.§ Physical setting

            Let us consider a perfect fluid with energy-momentum tensor
            
    T^μ_ν = (ρ + P)u^μ u_ν + Pδ^μ_ν ,

            in a flat Robertson-Walker (RW) metric
            
    s^2 = a^2(τ)(-τ^2 + δ_ij x^i x^j) .

            Now we will consider the situation where the fluid  possesses a bulk velocity with respect
            to the frame in which the metric takes the form (<ref>). Parameterizing the
            four-velocity as
            
    u_μ = aγ(-1, v_i) ,

            from the normalization condition, u_μ u^μ =-1, we have
            
    γ = 1/√(1-v^iv_i) ,

            where the spatial indices in v_i are lowered and raised using δ_ij. With this
            parameterization the components of the energy-momentum tensor are
            
            
    T^0_0   = -ρ - (ρ + P)γ^2v^2 ,
    
                    T^0_i   = (ρ + P)γ^2 v_i ,
    
                    T^i_j   = Pδ^i_j + (ρ + P)γ^2v^iv_j .

            
            Since the non-diagonal components are not zero, this moving fluid cannot act as a source for
            the geometry (<ref>). Let us show how to construct a valid source for the homogeneous
            and isotropic metric (<ref>) using a collection of fluids.

            
            

  
Isotropy. For non-relativistic fluids, the only non-diagonal component, to first
                order in v, is
                
    T^0_i = (ρ + P)v_i .

                If instead of a single fluid we have several fluids in relative motion, they can act
                as a source for (<ref>) if they satisfy
                
    T^0_i = ∑_s T_s^0_i = ∑_s(ρ_s + P_s)v_s i = 0 .

                The physical content of this condition is that of a kind of center of mass frame
                condition. An isotropic source can be constructed, to first order in v, out of two
                non-relativistic fluids if the net flux of momentum of one fluid is counterbalanced by
                that of the other fluid. We will see later that this constraint is conserved in time, 
                so it can be implemented with an appropiate choice of the initial conditions. In 
                <cit.> it is discussed how to transform to this frame, starting
                from an arbitrary configuration of the fluids.

            
            

  
Homogeneity. Homogeneity is easily implemented when the fluids are at rest, but
                we need to be cautious in our context. Consider two observers locally related by a boost
                
                    
  * (τ, x̌), 𝒪 frame in which (<ref>) is
                        satisfied and the metric takes the form (<ref>).
                    
  * (τ̃, x̌̃̌), 𝒪̃ frame moving with respect
                        to 𝒪 with velocity β̌.
                
                In the 𝒪̃ frame, the transformed coordinates are obtained applying a local
                Lorentz transformation
                
    x̃^μ = Λ^μ_ν(β) x^ν ,

                and the metric looks inhomogeneous
                
    s^2 = a^2(τ(τ̃, x̃̌̃))(-τ̃^2
                            +δ_ijx̃^ix̃^j) .

                The same applies to other time-dependent quantities like ρ and P. To provide
                a consistent source, we will require that the energy-momentum tensor of each fluid is
                homogeneous in the 𝒪 frame, i.e. the frame that observes an isotropic and
                homogeneous metric, not in the comoving frame with the fluid.

                                            
            Notice that on the tangent space at a given space-time point, we can always define an
            orthonormal basis ě_a given by
            
    ě_a = e^μ_ax^μ ,     e^μ_a=a^-1(τ)δ^μ_a .
            
            and the corresponding orthonormal basis in the 𝒪̃ frame reads
            
    ẽ̌_a = ẽ^μ_ax̃^μ ,
                            ẽ^μ_a=a^-1(τ(x̃))δ^μ_a .

            Thus, as expected, the two basis are just related by a Lorentz transformation on the tangent
            space
            
    ẽ̌_a = (Λ^-1(β))^b_aě_b .

                
            To sum up, we can source a flat RW metric (<ref>) with a collection of 
            non-relativistic moving fluids as long as
            
                
  * We are in the center of mass frame, where
                    
    ∑_s (ρ_s + P_s)v_s^i = 0 .

                
  * The energy-momentum tensor is homogeneous in that frame
                    
    ∂_iT^μ_ν = 0 .

            
            
        

  §.§.§ Evolution

            Finally, in this subsection we will analyze the evolution of a perfect fluid with bulk 
            velocity. Assuming homogeneity
            
    ∂_i T^μ_ν = 0 ,

            the conservation equation
            
    ∇_μ T^μ_ν = 0  ,

            for a flat RW metric in conformal time (<ref>) yields the equations of motion
            
    ∂_0(a^3T^0_0)    = ℋa^3T^i_i ,
                        
    ∂_0(a^4T^0_i)    = 0 ,

            where ℋ≡ȧ/a is the conformal Hubble parameter and 
            ≡∂_0≡∂/∂τ. As it was anticipated, 
            (<ref>) proves that the center of mass constraint 
            (<ref>) is conserved in time. Writing explicitly the components 
            (<ref>), the equations of motion can be rewritten in terms of v_i and
            ρ. The equation for the velocity can be expressed
            as
            
    ∂_0v_i - ∂_0log((1+w)a^4γ^2ρ)v_i = 0 ,

            where w≡ P/ρ is the equation of state of the fluid.
            In the absence of interactions, the velocity does not change its direction, so we only need
            to follow the evolution of its magnitude.
            Combining (<ref>) and (<ref>), it is possible
            to obtain
            
    ρ̇   = (v^2-3)(1+w)/1-wv^2ℋρ + ẇ/1-wv^2v^2ρ ,
                        
    v̇   = (1-v^2)(3w-1)/1-wv^2ℋv 
                        + ẇ/1+w1-v^2/1-wv^2 v ,

            assuming w≠ -1.
            It is worth particularizing these results to two kind of fluids.
            
                
  * Radiation, w = 1/3.
                    
    ρ̇   = -4ℋρ ,
    v̇   = 0 .

                    The fluid moves with constant velocity and with the usual scaling ρ∝ a^-4.
                
  * Matter, w = 0.
                    
    ρ̇   = (v^2-3)ℋρ ,
    v̇   = -(1-v^2)ℋv .

                    The equations can be solved analitically in this case
                    
    ρ   = ρ_0/a^2√(v_0^2 + a^2(1-v_0^2)) ,
    
                            v    = v_0/√(v_0^2 + a^2(1-v_0^2)) , 
    γ^2ρ   = ρ_0/a^4(1-v_0^2)√(v_0^2+a^2(1-v_0^2)) ,

                    where v_0 and ρ_0 are the velocity and density today. If the fluid starts with 
                    ultrarelativistic initial conditions, it behaves as radiation γ^2ρ∝ a^-4
                    until the velocity drops down and it enters the non-relativistic regime. In the
                    non-relativistic regime, to first order in v, the velocity slows down with the
                    expansion v∝ a^-1 and the density scales as usual 
                    γ^2ρ≃ρ∝ a^-3.
            
            Analytic expressions for a generic equation of state w(a) can be obtained in the regime
            of small velocities
            
    ρ   = ρ_0 exp(-3∫ a/a(1+w)) + v^2 ,
    
                    v    = v_0(1+w_0)/a^4(1+w)exp(3∫ a/a(1+w)) + v^2 ,

            where w_0 is the value of the equation of state today. For the particular case w=const., 
            we have <cit.>
            
    ρ   = ρ_0  a^-3(1+w) + v^2 ,
    
                    v    = v_0  a^-(1-3w) + v^2 .

        
    

 §.§ Kinetic approach

        Non-relativistic particles and the photon-baryon plasma in the early Universe can be well
        approximated as perfect fluids, and analyzed using only their energy-momentum tensor. However, 
        to describe the free-streaming of neutrinos or the details of decoupling, the fluid approximation
        is not enough and we must use a more general kinetic approach. In this section we follow closely 
        the presentation of <cit.>.

        
        The phase-space of the system is described by
        
            
  * Three positions x^i.
            
  * Three conjugate momenta P_i. These conjugate momenta are defined as the spatial
                components of the four-momentum
                
    P^μ≡ m x^μ/λ ,      P_μ P_ν g^μν = -m^2 ,

                where λ≡√(- s^2) is the proper time and the spatial index i has
                been lowered with the full metric g_μν.
        
        
        The number of particles per unit of phase-space volume is 
        
    N = g_* f(τ, x^i, P_j) ^3 x ^3 P/(2π)^3 ,

        where g_* is the number of internal degrees of freedom, e.g. the number of helicity states, and 
        f is the phase-space distribution function. In the kinetic approach, the 
        energy-momentum tensor can be defined as <cit.>
        
    T_μν =g_*/(2π)^3∫^3 P  (-g)^-1/2P_μ P_ν/P^0f(τ, 
    			    x^i, P_j) ,

		where g is the determinant of the metric g_μν. In a cosmological setting, the next step
		would be to particularize these definitions to a RW metric (<ref>) and then to 
		study perturbations over the metric and the distribution function. We will split the discussion
		in two parts. Since our main modification with respect to standard comoving cosmology concerns 
		the definition of the unperturbed distribution function, we will focus on the background in the 
		next section <ref>. Section <ref> contains the complete 
		treatment of perturbations.
        
	    

  §.§.§ Background

	        Starting again with the metric 
	        
    s^2 = a^2(τ)(-τ^2 + δ_ij x^i x^j) ,

	        and denoting the momentum as
	        
	        
    ϵ   ≡ -P_0 ,
    
    	            q_i    ≡ P_i ,

	        
	        such that q^i≡ q_i, we obtain from the mass-shell condition (<ref>)
	        
    ϵ^2 = m^2a^2 + q^2 .

	        We introduce again two related frames.
	        
	            
  * 𝒪, frame in which the metric takes the form (<ref>).
	            
  * 𝒪̃, frame moving with respect to 𝒪 with velocity β̌.
	        
	        The local Lorentz transformation (<ref>) that connects both 
	        frames yields
	        
	        
    ϵ̃   ≡Λ_βϵ = γ(ϵ - q̌·β̌) ,
    q̃^i    ≡Λ_β q^i = 𝒫^i_jq^j - γϵβ^i ,
    β̌·q̃̌̃   = γ (q̌·β̌-ϵβ) .

	        
	        We have defined
	        
	        
    𝒫^i_j    ≡δ^i_j + (γ - 1)β̂^iβ̂_j ,
    γ   ≡ (1-β^2)^-1/2 ,

	        
	        where β̂^i is a unit vector along β̌ and every spatial index has been 
	        lowered or raised with δ_ij. Next we consider a homogeneous distribution function 
	        in the 𝒪 frame
	        
    f(τ, x̌, q̌) = f_0(τ, q̌) .

	        With this distribution function, we can define the usual fluid quantities
	        
	            ρ   ≡a^-4 g_*∫^3 q/(2π)^3 ϵ  f_0 ,           P    ≡a^-4 g_*∫^3 q/(2π)^3  q^2/3ϵ  f_0 ,

	            Q^i    ≡a^-4 g_*∫^3 q/(2π)^3  q^i f_0 ,           n    ≡a^-3 g_*∫^3 q/(2π)^3  f_0 ,

	            Π^ij    ≡a^-4 g_*∫^3 q/(2π)^3(q^iq^j/ϵ
	                    -q^2/3ϵδ^ij) f_0 ,           V^i    ≡a^-3 g_*∫^3 q/(2π)^3  q^i/ϵ  f_0 ,
	        
	        that represent the usual energy, momentum, shear tensor, pressure, number and velocity
	        densities of the fluid. To relate this set of quantities with those computed in the 
	        boosted 𝒪̃ frame, we can either use their tensorial character under
	        local Lorentz transformations or the fact that f transforms as a scalar
	        
    f̃_0(τ̃, x̌̃̌, q̌̃̌) = f_0(τ(τ̃, x̃̌̃), q̌(q̌̃̌)) ,

	        where, from now on, we will denote f̃_0(τ̃, x̌̃̌, q̌̃̌) just as f̃_0(τ,q̌̃̌).
	        With this property and the Lorentz-invariant volume element we can write, for instance,
	        
    Q̃^i = a^-4g_*∫^3q̃/(2π)^3 q̃^i f̃_0(τ, q̃̌̃)
    	                = a^-4g_*∫^3q/(2π)^3(Λ_βϵ)(Λ_β q^i)/ϵ
    	                    f_0(τ, q̌) .

	        Using this procedure, one can obtain
	        
	        
    ρ̃   = ρ + γ^2(β_kβ_lΠ^kl-2β_kQ^k)
    	                    + γ^2β^2(ρ + P) ,
    Q̃^i    = γ𝒫^i_j(Q^j-β_kΠ^kj)
    	                    -γ^2β^i(ρ+P-Q^jβ_j) ,
    Π̃^ij   = (𝒫^i_k𝒫^j_l-1/3γ^2β_kβ_l
    	                        δ^ij)Π^kl
    	                    - γ(β^i𝒫^j_k+β^j𝒫^i_k
    	                        -2/3γδ^ijβ_k)Q^k
    	                    + γ^2(β^iβ^j-1/3β^2δ^ij) (ρ + P) ,
    P̃   = P + 1/3γ^2(β_kβ_lΠ^kl - 2β_kQ^k)
    	                    + 1/3γ^2β^2(ρ+P) ,
    ñ   = γ (n-V^jβ_j) ,
    Ṽ^i    = 𝒫^i_jV^j- γβ^i n  .

	        
	        It is important to stress that the preceeding relations hold as well if the quantities are
	        defined with the full distribution function f, instead of using just the background part
	        f_0, and we will make use of them when we study perturbations.

	        
	        There remains the question of how to describe the moving fluids of section <ref>
	        in terms of a distribution function. We will describe the different constituents of the universe
	        with an unperturbed distribution function that satisfies
	        
    f_0(τ, q̌) = f̃_0(τ, q̃) .

	        That is, the distribution function is homogeneous in the 𝒪 frame, i.e. the frame
	        that observes a homogeneous and isotropic universe, and isotropic in the 𝒪̃
	        frame, i.e. the frame comoving with the fluid. This parallels the discussion in section
	        <ref> and allows us to describe a fluid moving with velocity β̌. 
	        The condition (<ref>) is the main physical assumption in our work.
	        For instance, applying it to a black-body spectrum for massless particles, using 
	        (<ref>) we obtain the usual boosted distribution function
	        
    f̃_0(q̃) = 1/^ q̃/T̃-1 
    	                = 1/^ q/T(q̌)-1 = f_0(q̌) ,     
    	                T(q̌)≡γ(1-q̌·β̌/q)T̃ .

	        
	        If the distribution function satisfies (<ref>) we have
	        
    Π̃_ij=0 ,    Ṽ_i = Q̃_i = 0 ,

	        so we are indeed in the comoving frame with the (perfect) fluid. In this case, the 
	        relation between both sets of fluid variables is
	        
	        
    ρ   = ρ̃ + γ^2β^2(ρ̃ + P̃) ,
    
    	            Q^i    = γ^2β^i(ρ̃+P̃) ,
    Π^ij   = γ^2(β^iβ^j-1/3δ^ijβ^2) 
    	                (ρ̃ + P̃) ,
    
    	            P    = P̃ + 1/3γ^2β^2(ρ̃+P̃) ,
    
    	            n    = γñ ,
    
    	            V^i    = γβ^i ñ .

	        
	        Using the definition of the energy-momentum tensor (<ref>), we can write its
	        components as
	        
	        
    T^0_0   = -ρ = -ρ̃ - γ^2(ρ̃+P̃)β^2 ,
    	                
    
    	            T^0_i   = Q_i = γ^2(ρ̃+P̃)β_i ,
    	                
    
    	            T^i_j   = Pδ^i_j + Π^i_j = P̃δ^i_j +γ^2(ρ̃
    	                +P̃)β^iβ_j .

	        
	        These expressions agree with the ones obtained for a perfect fluid (<ref>). 
	        Both approaches are equivalent at this level. To first order in β we have
	        
	        
    T^0_0   = -ρ̃ ,
    
    	            T^0_i   = (ρ̃+P̃)β_i ,
    
    	            T^i_j   =  P̃δ^i_j .

	        
	        The full energy-momentum for a collection of fluids is homogeneous, as it was imposed in 
	        (<ref>), and it is also isotropic, to first order in β, if the 
	        velocities of the fluids satisfy the constraint
	        
    ∑_s (ρ̃_s + P̃_s)β_s^i = 0 ,

	        which is the same condition obtained in (<ref>). It is clear from 
	        (<ref>) that a similar constraint can always be imposed, to all orders
	        in β, to achieve T^0_i=0 but this is not enough to source a RW geometry.
	        Already to second order in β, (<ref>) contains a quadrupolar
	        anisotropy that cannot be compensated by the other fluids. In this case, we should go one
	        step further and consider a Bianchi universe. However, for the values of β that we will
	        consider, this quadrupole lies well below the observed value <cit.>. 
	        Therefore, in this work we will restrict ourselves to first order and a RW background.
        
	    

  §.§.§ Perturbations

	        Our starting point now is a perturbed flat RW metric
	        
    s^2 = a^2(τ)(-(1-A)τ^2 + 2B_iτ  x^i 
    	                + (δ_ij+H_ij) x^i x^j) .

	        Reparameterizing the momentum as
	        
	        
    P_0   ≡ -aE + δ P_0 ,
    
    	            P_i   ≡ a(δ^j_i +1/2H^j_i)p_j ,

	        
	        from the mass-shell condition (<ref>) we obtain
	        
	        
    E^2    = m^2 + p^2 ,
    δ P_0    = 1/2A aE  + ap_iB^i .

	        
	        These expressions can be regarded just as redefinitions but they have a very simple 
	        physical interpretation in terms of the vierbein (<ref>). 
	        Our choice of momentum, P_μ = e_μ^ap_a, corresponds to a choice of vierbein 
	        with components
	        
	        
    e_0^0   = a(1-1/2A) ,
    
    	            e_i^0   = 0 ,
    
    	            e_0^i   = aB^i ,
    
    	            e_i^j   = a(δ^i_j+1/2H^i_j) .

	        
	        With these definitions, p_i are the momenta measured by a locally inertial observer at
	        a fixed spatial position. It is convenient to work with a closely related set of variables
	        defined as
	        
    q_i≡ ap_i ,    ϵ≡ aE .

	        As we will see, written in terms of q_i, the Boltzmann equation in section <ref>
	        does not contain a zero-order term. 
	        We perturb the phase-space distribution accordingly
	        
    f(τ, x̌, q̌) = f_0(τ, q̌) + δ f(τ, x̌, q̌) ,

	        and define the corresponding perturbed fluid variables
	        
	            δρ   ≡a^-4 g_*∫^3 q/(2π)^3 ϵ  δf ,           δP    ≡a^-4 g_*∫^3 q/(2π)^3  q^2/3ϵ  δf ,

	            δQ^i    ≡a^-4 g_*∫^3 q/(2π)^3  q^i δf ,           δn    ≡a^-3 g_*∫^3 q/(2π)^3  δf ,

	            δΠ^ij    ≡a^-4 g_*∫^3 q/(2π)^3(q^iq^j/ϵ
	                    -q^2/3ϵδ^ij) δf ,           δV^i    ≡a^-3 g_*∫^3 q/(2π)^3  q^i/ϵ  δf .
	        
	        The components of the perturbed energy-momentum tensor are
	        
	        
    δ T^0_0   = -δρ + B^iQ_i ,
    δ T^0_i   = δ Q_i +1/2AQ_i + 1/2H^j_iQ_j ,
    δ T^i_j   = δ Pδ^i_j + δΠ^i_j -B^iQ_j 
    	                + 1/2(H^k_jΠ^i_k-H^i_kΠ^k_j) .

	        
	        Expressing the background quantities in the 𝒪̃ frame
	        
	        
    δ T^0_0   = -δρ +γ^2(ρ̃+P̃)B_iβ^i ,
    δ T^0_i   = δ Q_i +1/2γ^2(ρ̃+P̃)
    	                (Aδ^j_i + H^j_i)β_j ,
    δ T^i_j   = δ Pδ^i_j + δΠ^i_j 
    	                -γ^2(ρ̃+P̃)B^iβ_j
    	                +1/2γ^2(ρ̃+P̃)(β^iH^k_jβ_k 
    	                - β_jH^i_kβ^k) .

	        
	        It is easy to show that, to first order in β, once we apply the center of mass condition
	        (<ref>), every metric variable cancels out and does not appear in the definition 
	        of the total energy-momentum tensor.
	        
	        

  §.§.§ Non-relativistic limit

	            We end this section with an aside about the non-relativistic (NR) limit (q/ϵ≪ 1). In the next sections
	            only massless (photons and massless neutrinos) and NR particles (baryons and cold dark matter)
	            will be taken into account. In particular, we will only keep the first NR order, 
	            neglecting the pressure and sound speed of baryons.

	            
	            Although this kind of approximation is standard, we must be careful when taking the NR limit
	            in a moving frame. The proper way to account for this limit is to take it in the frame comoving
	            with the fluid, i.e. 𝒪̃. To first NR order we have
	            
    δΠ̃^ij≃ 0 ,    δP̃≃ 0 .

	            However, in the 𝒪 frame, using (<ref>), it can be
	            seen that we do have pressure and anisotropic stress
	            
	            
    δ P    ≃2/3γ^2β_kδQ̃^k 
    	                    + 1/3γ^2β^2δρ̃ ,
    δΠ^ij   ≃γ(β^i𝒫^j_k+β^j𝒫^i_k
    	                    -2/3γδ^ijβ_k)δQ̃^k
    	                    +γ^2(β^iβ^j-1/3β^2δ^ij)δρ̃ .

	            
	            To first NR order and to first order in β, we have the following useful results
	            
	                δP    ≃2/3β_kδQ^k ,           δV^i    ≃δQ^i/m ,

	                δΠ^ij    ≃(β^iδ^j_k+β^jδ^i_k 
	                        - 2/3δ^ijβ_k)δQ^k ,           δṼ^i    ≃1/m(δQ^i - β^iδρ) ,

	                δñ    ≃δρ̃/m ≃1/m(δρ-2β_kδQ^k) ,           δn    ≃1/m(δρ-β_kδQ^k) .
	            
	            Note that the preceeding results hold as well for the corresponding unperturbed quantities,
	            see (<ref>).
	            Finally, the full energy-momentum tensor for a non-relativistic species to first order in
	            β is
	            
	            
    T^0_0+δ T^0_0   ≃ -ρ̃
    	                    -δρ+ρ̃β^iB_i ,
    
    	                T^0_i+δ T^0_i   ≃ m δ V_i
    	                    +ρ̃(δ^j_i+1/2δ^j_iA+1/2H^j_i)β_j ,
    
    	                T^i_j+δ T^i_j   ≃ 
    	                    m (β^iδ V_j - β_jδ V^i) + ρ̃β_jB^i .

	            
        


§ BOLTZMANN EQUATION

    After writing the energy-momentum tensor in terms of the distribution function, we need to compute
    its time evolution. This information is encoded in the Boltzmann equation, which in the locally
    Minkowskian frame takes the form
    
    Df/ t = 𝒞[f] .

    where t≡ e_μ^0 x^μ is the time measured by the locally Minkowskian 
    observer, which for the vierbein choice in (<ref>) reads 
    t= a(1-1/2A)τ. The left-hand side, the so-called Liouville operator, describes 
    the free streaming of particles in phase space. It is defined as
    
    Df/ t≡ft + x̣^̣ịtfx^i + q̣^̣ịtfq^i .

    This operator contains the information about the space-time geometry, through its effects on the
    geodesics of the particles. The functional on the right-hand side of (<ref>) 
    is the so-called collision term. It takes into account how the number of particles per unit of 
    phase-space volume change due to collisions, i.e. local interactions. Hence, it does not contain 
    information about the underlying geometry <cit.>. The collision term takes the same 
    form as in flat space-time when written in terms of the momenta measured by a locally inertial 
    observer p_a defined above. The final form of the Boltzmann equation in conformal time is
    
    Df/τ = a(1-1/2A)𝒞[f] .

    The next section is devoted to the left-hand side of the Boltzmann equation, i.e. the Liouville operator,
    particularizing to massless and non-relativistic particles. This analysis exhausts all the information
    needed to follow the evolution of non-interacting particles, like CDM and neutrinos. However, to 
    describe the photon-baryon plasma, we must move on to the full Boltzmann equation. Section 
    <ref> studies the interaction between photons and electrons.

    
    

 §.§ Liouville operator

        In order to compute the time derivatives appearing in (<ref>), we need the geodesics in the metric (<ref>). The whole computation
        of the geodesics can be found in detail in the appendix <ref>. Using the 
        definition of the four-momentum (<ref>) and the parameterization 
        (<ref>), the final results are
        
        
    x̣^̣ịτ   = q^i/ϵ(1-1/2A)-B^i
    	            -1/2ϵH^i_kq^k ,
    q̣^̣ịτ   = 1/2ϵ∂_iA + q^jC_ij 
                    + q^jq^k/ϵD_ijk ,

        
        where the following combinations of metric variables have been defined
        
        
    C_ij   ≡∂_iB_j - 1/2Ḣ_ij ,
    
                D_ijk   ≡1/2(∂_iH_jk-∂_kH_ij) .

              
        The left-hand side of the Boltzmann equation (<ref>) is evaluated
        in the 𝒪 frame, where the geodesics are computed,
        
    Df/τ≡fτ + x̣^̣ịτfx^i + q̣^̣ịτfq^i .

        We will assume that the distribution function takes the form
        
    f(τ, x̌,q̌)    = f_0(τ, q̌) + δ f(τ, x̌, q̌)
       = f̃_0(τ, Λ_βϵ) 
                                         + δ f(τ, x̌, q̌) ,

        where f̃_0(τ, ϵ̃) is the standard isotropic distribution in the 
        𝒪̃ frame. Since (<ref>) is already first order in perturbations, 
        the left-hand side of the Boltzmann equation can be recast as
        
    Df/τ = f_0τ + δ fτ + x̣^̣ịτδ fx^i
                    +q̣^̣ịτf_0q^i .

        We will restrict our discussion to massless and non-relativistic massive particles.
        
        

  §.§.§ Massless particles

            For massless particles the Lorentz transformations (<ref>) take a simpler form
            
            
    q̃   = γ(1-n̂·β̌)q ,
    q̃^i    = (𝒫^i_jn^j-γβ^i)q ,

            
            where we have splitted the momentum into direction and magnitude
            
    q^i≡ q  n^i ,     n_iδ^ijn_j=1 .

            Using the Boltzmann equation, it can be directly checked that the unperturbed distribution 
            function in the 𝒪̃ frame only depends on q̃, if there is no 
            zero-order collision term,
            
    f̃_0(τ,q̃) = f̃_0(q̃) .

            It is convenient to work with the reduced phase-space density, integrating out the dependence 
            on the momentum magnitude, defined as
            
    ℱ(τ, x̌, n̂) ≡1/𝒩̃∫ q^3 q δ f(η, x̌, q̌) ,

            where the constant 𝒩̃ is related to the comoving energy density
            
    𝒩̃≡∫q̃^3q̃ f̃_0(q̃)
                        = 2π^2/g_*a^4ρ̃ .

            Plugging the geodesic equation (<ref>) into the right-hand
            side of (<ref>) and integrating over the
            momentum magnitude, we obtain, to all orders in β,
            
    ∫ q^3 q Df/τ = 𝒩̃[η1/γ^4
                        (1-n̂·β̌)^4   + ℱ̇ + n^i∂_iℱ
                        - 4/γ^4(1-n̂·β̌)^5(1/2n^i∂_iA
                        +n^in^jC_ij
       -1/2β^i∂_iA - β^in^jC_ij -β^in^jn^kD_ijk)] .

            Expanding it to first order in β, we have
            
    ∫ q^3 q Df/τ = 𝒩̃[4n̂·β̇̌̇   +ℱ̇ 
                        + n^i∂_iℱ -4(1/2n^i∂_iA + n^in^jC_ij)
                        (1+5n̂·β̌)
       +4β^i(1/2∂_iA+ n^jC_ij+n^jn^kD_ijk)] .

            The first moments of the angular distribution can be obtained performing the appropiate
            integrals
            
    ∫Ω/4π∫ q^3 q Df/τ   = 𝒩̃[δ̇+4/3∂_iδ v^i-4/3C_ijδ^ij-4/3β^i∂_iA-4/3D_ijkδ^ijβ^k] ,
                        
    ∫Ω/4π n^i∫ q^3 q Df/τ   = 4/3𝒩̃[β̇^i+δv̇^i 
                        + 3/4∂_jπ^ij
                        +1/4∂^iδ - 1/2∂^iA
                        -(δ^ijβ^k+δ^jkβ^i)C_jk] ,

            where we have defined
            
            
    δ   ≡∫Ω/4π ℱ 
                        = δρ/ρ̃ ,
                        
    δv̌   ≡3/4∫Ω/4π n̂ ℱ
                        = δQ̌/ρ̃+P̃ ,
                        
    π^ij   ≡∫Ω/4π(n^in^j-1/3δ^ij)
                        ℱ
                        = δΠ^ij/ρ̃ .

            
            
        

  §.§.§ Massive particles

            The results for massive particles are much more involved. In this case, we must use the
            full expressions for the Lorentz transformations (<ref>) and the
            geodesics (<ref>). However, since we will focus on non-relativistic
            particles, we can simplify the analysis restricting it to the first moments of the 
            Boltzmann equation. Summing up, in this section we compute the evolution of the number 
            density, energy and velocity perturbations.
 
            
            The energy density contrast, equation of state and sound speed are defined as
            
    δ≡δρ/ρ̃ ,    
                    w ≡P̃/ρ̃ ,    
                    c_s^2 ≡δ P/δρ .

            Note that in these definitions the background quantities are referred to the 
            𝒪̃ frame while the perturbed quantities are defined in the 𝒪
            frame. The election of intermediate variables is a matter of choice, the only condition
            being that we write the energy-momentum tensor consistently in terms of these variables. We
            stick to this convention throughout this work. The final results for the first moments of
            the distribution are
            
                
  * Number density.
                    
    a^-3g_*∫^3q/(2π)^3Df/τ = a^-3τ(γ a^3ñ) + a^-3τ(a^3δ n)+ ∂_iδ V^i -γñ[1/2β^i∂_iA + δ^ijC_ij+D_ijkδ^ijβ^k] .

                
  * Energy.
                    
    a^-4g_*∫^3q/(2π)^3ϵDf/τ =    τ[
                                γ^2(ρ̃+β^2P̃)]
                                +3ℋ(ρ̃+P̃)(1+4/3β^2γ^2)
                                +∂_iδ Q^i
       + ρ̃(δ̇+3ℋδ (c_s^2-w))
                                +δ(ρ̇̃̇+3ℋ(ρ̃+P̃))
                                
       -γ^2(ρ̃+P̃)[β^i∂_iA + δ^ijC_ij
                                +β^iβ^jC_ij+ D_ijkδ^ijβ^k] .

                
  * Momentum.
                    
    a^-4g_*∫^3q/(2π)^3q^i Df/τ =    a^-4τ[γ^2
                                β^i a^4(ρ̃+P̃)] + a^-4τ(a^4δ Q^i )
                                + ∂_j(δΠ^ij + δ^ijδ P)
       -γ^2(ρ̃+P̃)[1/2(δ^il
                                +β^iβ^l)∂_lA+(δ^jkβ^i
                                +δ^ijβ^k)(C_jk+D_jklβ^l)] .

            
            These results are exact to all orders in β, for relativistic and non-relativistic
            particles alike. They can be shown to reproduce (<ref>) and
            (<ref>) for massless particles. Next we define
            
    δ_n≡δ n/ñ ,    δ v^i≡δ V^i/Ñ .

            Assuming that the zero-order Boltzmann equation (without collisions) is satisfied, so we keep 
            only cosmological perturbations and terms with β, and expanding to first NR order and 
            to first order in β we have
            
                
  * Number density.
                     
    a^-3g_*∫^3 q/(2π)^3Df/τ≃Ñ{δ̇_̇ṅ +∂_iδ v^i -1/2β^i∂_i A
                                 -δ^ijC_ij-D_ijkδ^ijβ^k} .

                
  * Energy.
                    
    a^-4g_*∫^3 q/(2π)^3ϵDf/τ≃ρ̃{δ̇+2ℋβ_kδ v^k +∂_iδ v^i
                                -β^i∂_iA-δ^ijC_ij-D_ijkδ^ijβ^k} .

                
  * Momentum.
                    
    a^-4g_*∫^3 q/(2π)^3q^iDf/τ≃ρ̃ {β̇^i +ℋβ^i + δv̇^i
                                   +(β^iδ^j_k+β^jδ^i_k)∂_jδ v^k
                                    +ℋδ v^i
       -1/2∂^iA-(δ^jkβ^i+δ^ijβ^k)
                                    C_jk} .

            
            
    

 §.§ Collision term

        This section is devoted to the calculation of the collision term for Compton scattering between
        electrons and photons. The notation in this section is slightly different. As we mentioned 
        before, the collision term must be written in terms of the momenta p_i measured by a 
        locally inertial observer at a fixed spatial position. They are related to the momenta q_i we 
        have been using as defined in (<ref>), i.e. q_i=p_i/a.

        
        The standard physical assumptions underlying the derivation of the collision term are 
        
            
  * When written in terms of the momenta p_i, the collision term is the same as in flat 
                space, since it takes into account local information where the curvature 
                effects are not important. In the same way, the matrix element ℳ is computed 
                using quantum field theory (QFT) in flat space.
            
  * The temperature of the plasma is low enough so that the electrons are non-relativistic.
                We keep only the first order correction in the non-relativistic expansion, so we keep the
                electron velocity but neglect its pressure and sound speed. We consider the NR limit of 
                Compton scattering, i.e. Thomson scattering.
            
  * Electrons and protons are much more tightly coupled between them than to the photons. The
                velocities of free electrons, protons and the full baryonic velocity are
                the same throughout the evolution.
            
  * The angular dependence of Thomson scattering is neglected. This angular dependence has 
                proven important for 1% accuracy and especially for polarization, but we will not take it
                into account in this work. We use the angle-averaged matrix element instead.
            
  * The number of internal degrees of freedom is not included in the definition of f, i.e.
                the equilibrium distributions correspond to the usual Bose-Einstein and Fermi-Dirac
                distributions.
            
  * The medium is diluted enough so that we can neglect the quantum statistical factors
                (1± f) responsible for the Bose enhancement and Pauli blocking effects. We do not
                take into account any plasma effect from finite temperature QFT.
        
        
        With these qualifications, our starting point is the following definition of the collision term
        for the process e(p_e)+γ (p)↔ e(p_e')+γ (p')
        
    𝒞[f(p)]=1/4p∫𝒟p_e𝒟p'𝒟p'_e 
    			    (2π)^4δ (p^μ +p_e^μ - p'^μ -p_e'^μ)
    			    [f(p̌'̌)f_e(p̌'̌_e)-f(p̌)f_e(p̌_e)]∑_spins |ℳ|^2 ,

		where 𝒟p=p̌^3/(2π)^32E is the Lorentz-invariant phase-space volume
		element and the dependence of the distribution functions on space-time coordinates has been 
		omitted since it does not play any role.	In our setting, we must implement the fact that the 
		fluids are moving. The collision term is defined in the cosmic center of mass, 𝒪 
		frame, and in this frame both photons and electrons have their own bulk velocity. We will 
		represent it schematically as
		
    𝒞[f] = 1/4p∫𝒟p_e𝒟p'𝒟p'_e(2π)^4
                    δ (p^μ +p_e^μ - p'^μ -p_e'^μ)[f̅(Λ_βp̌')f̃_e(Λ_β_ep̌_e')
                    -f̅(Λ_βp̌)f̃_e(Λ_β_ep̌_e)]∑_spins
                    |ℳ|^2 ,

		where f̅ and f̃_e are the distribution functions of photons and electrons in their
		frame, moving with bulk velocities β̌ and β̌_e, respectively, with respect
		to the 𝒪 frame.

		
		Previously, the 𝒪̃ frame was defined as the frame comoving with the fluid. In
		this case, we are facing two moving fluids. We take 𝒪̃ to be the frame
		moving with velocity β̌_e with respect to 𝒪, i.e. the frame comoving with
		the electrons. Performing the integration in this frame, we have
		
    𝒞[f] = 1/4p∫𝒟p̃_e𝒟p̃'𝒟p̃_e'(2π)^4
                    δ (p̃^μ +p̃_e^μ - p̃'^μ -p̃_e'^μ)[
                    f̅(Λ_βΛ_β_e^-1p̌̃̌')f̃_e(p̌̃̌_e')
                    -f̅(Λ_βΛ_β_e^-1p̌̃̌)f̃_e(p̌̃̌_e)]
                    ∑_spins|ℳ|^2 .

        The previous two equations may seem devoid of any additional content with respect to
        (<ref>). As they stand, without defining f and f_e, they correspond
        just to a renaming of functions and reshuffling of variables. The physical content lies in
        (<ref>), i.e. in the structure and relation of the background distribution
        function in 𝒪 and 𝒪̃. The 𝒪̃ frame is 
        comoving with the electrons and observes a standard isotropic equilibrium distribution. It is
        in this frame that we can perform the usual NR expansion <cit.> to get
        
    𝒞[f]=σ_T/4π p∫p̃'p̃'Ω̃'
                    [ñ_e^fullδ (p̃-p̃')
                    +ñ_eũ̌^full_e·
                    (p̌̃̌-p̌̃̌')δ (p̃-p̃')p̃']
                    (f̅(Λ_βΛ_β_e^-1p̌̃̌')
                    -f̅(Λ_βΛ_β_e^-1p̌̃̌)) ,

        where σ_T is the Thomson cross section and we have defined
        
    ñ_e^full≡ 2∫^3p̃_e/(2π)^3f̃_e(p̌̃̌_e) ,
                    ñ_eũ̌_e^full≡ 2∫^3p̃_e/(2π)^3p̌̃̌_e/Ẽ_p_ef̃_e(p̌̃̌_e) ,

        where f̅, f̃_e correspond to the full distribution functions. Just as we did for the left-hand
        side in section <ref>, we will split the distribution functions into a 
        background and a perturbation part and integrate out the magnitude of the photon momentum. The
        whole process, to all orders in β, is detailed in the appendix <ref>.
        Here we present only the final results. To first order in β, we get
        
    1/𝒩̃∫ q^3 q 𝒞[f] = ñ_eσ_T [   -(1-n̂·β̌_e)ℱ_γ 
                    + (1+3n̂·β̌_e)δ_γ
                    -8/3β̌_e·δv̌_γ -4n̂·Δβ̌ 
                    - 8/3δv̌_e·Δβ̌
       +4n̂·δv̌_e n̂·Δβ̌+4δv̌_e·(n̂
                    -β̌_e + 4n̂ n̂·β̌_e) 
                    -4n̂·β̌ δ n_e/ñ_e] ,

        where ñ_e is the number density 
        of free electrons in the 𝒪̃ frame, δ n_e and δv̌_e are defined in 
        (<ref>) and the difference of velocities is
        
    Δβ̌≡β̌-β̌_e .

        Since, to first order in β, the background quantities like ρ or n concide in 
        the 𝒪 and 𝒪̃ frames, we will drop the distinction.
        The first two moments of the photon collision term are
        
    1/𝒩̃∫Ω/4π∫ q^3 q 𝒞[f(p̌)]
                       = -4/3n_eσ_T[β̌_e·(δv̌_γ 
                        -δv̌_e) + δv̌_e·Δβ̌] ,
                        
    1/𝒩̃∫Ω/4π n^i∫ q^3 q 𝒞[f(p̌)]
                       = -4/3n_eσ_T[δ v_γ^i-δ v^i_e + Δβ^i
                        +β^iδ n_e/n_e-β^i_eδ_γ - 3/4β_e jπ_γ^ij] .

        
        

  §.§.§ Conserved quantities

            We have just computed the collision term for photons. However, the whole
            plasma is described by the coupled system
            
            
    Df/ t   = 𝒞[f,f_e] ,
    Df_e/ t   = 𝒞_e[f,f_e] .

            
            We ought to compute the collision term for electrons 𝒞_e as well. Not surprisingly,
            both terms are not independent. In fact, we can make use of some conservation laws derived
            from the full Boltzmann equation to save us most of the work. Following <cit.>,
            we write both collision terms with the compact notation
            
            
    c_eγ   =1/21/2E_p_e1/2p1/2E_p'_e1/2p'
    		    		    (2π)^4δ (p^μ +p_e^μ - p'^μ -p_e'^μ)
            			    [f(p̌'̌)f_e(p̌'̌_e)-f(p̌)f_e(p̌_e)]∑_spins |ℳ|^2 ,
    𝒞_e[f_e(p_e)]    ≡∫^3p/(2π)^3^3p'/(2π)^3^3p'_e/(2π)^3 c_eγ≡⟨ 
    		    						c_eγ⟩_pp'p'_e ,
    𝒞[f(p)]    ≡⟨ c_eγ⟩_p'p_e'p_e .
   
		    	
		    Integrating over all the momenta, it is easy to see that we have
		    	
    ⟨ c_eγ⟩_p_epp_e'p'   = 0 ,
    ⟨ (p+E_p_e) c_eγ⟩_p_epp_e'p'   = 0 ,
    ⟨ (p̌+p̌_e) c_eγ⟩_p_epp_e'p'   = 0 ,

		    	corresponding to the conservation of the number of particles, energy and momentum. Using
		    	these results, the following equalities hold
		    	
    ∫^3 p_e/(2π)^3 E_p_e 𝒞_e[f_e(p_e)] 
    		    	           = -∫^3 p/(2π)^3 p 𝒞[f(p)] ,
    		    	        
    ∫^3 p_e/(2π)^3 p̌_e 𝒞_e[f_e(p_e)] 
    		    	           = -∫^3 p/(2π)^3 p̌ 𝒞[f(p)] .

		    	This means that we can compute the first two moments of the Boltzmann equation for electrons,
		    	the only ones that we will need since they are non-relativistic, from the first two moments
		    	of the photons, already computed in (<ref>) and 
		    	(<ref>).
        
    

 §.§ Boltzmann equation for different components

        

  §.§.§ Photons

            The reduced Boltzmann equation for photons is obtained combining the Liouville operator 
            (<ref>) and the collision term (<ref>). To zero
            order in cosmological perturbations, it describes the evolution of the bulk velocity β
            
    β̇^i = -an_eσ_TΔβ^i .

            To first order in cosmological perturbations and β, we get the evolution of
            the reduced phase-space density
		    
    ℱ̇_γ   + n^i∂_iℱ_γ 
    		            -4(1/2n^i∂_iA + n^in^jC_ij)
    		            (1+5n̂·β̌)
    		            +4β^i(1/2∂_iA+ n^jC_ij+n^jn^kD_ijk)
    		            
       = an_eσ_T[-(1-n̂·β̌_e)ℱ_γ 
    		            + (1+3n̂·β̌_e)δ_γ-8/3β̌_e·δv̌_γ
    		            -8/3δv̌_b·Δβ̌
                 +4n̂·δv̌_b n̂·Δβ̌+4δv̌_b·(n̂-β̌_e + 4n̂ n̂·β̌_e)-4n̂·β̌ δ n_e/n_e+2n̂·Δβ̌A] .

		    Since protons and electrons form a single tightly coupled baryonic fluid,
		    we have substituted δv̌_e with δv̌_b, the velocity of baryons.
		    The evolution of the fluid variables can be obtained performing the appropiate angular
		    integrals. The equations for the density, combining 
		    (<ref>) and (<ref>), and the velocity, combining
		    (<ref>) and (<ref>), are
		    
    δ̇_γ+4/3∂_iδ v^i_γ-4/3C_ijδ^ij-4/3β^i∂_iA-4/3D_ijkδ^ijβ^k
    		           = -4/3an_eσ_T[β̌_e·(δv̌_γ 
    		                        -δv̌_b) + δv̌_b·Δβ̌] ,
    		                        
    δv̇^i_γ + 3/4∂_jπ^ij_γ
    		            +1/4∂^iδ_γ - 1/2∂^iA
    		            -(δ^ijβ^k+δ^jkβ^i)C_jk   = -an_eσ_T[δ v_γ^i-δ v^i_b +β^iδ n_e/n_e-β^i_e
    		            δ_γ - 3/4β_e jπ_γ^ij-1/2Δβ^iA] .

		    
        

  §.§.§ Baryons

            The evolution of the baryon density can be found using the left-hand side 
            (<ref>) and energy conservation (<ref>).
            For the velocity, we must use the left-hand side (<ref>) and
            momentum conservation (<ref>). As mentioned before, since they are much 
            more tightly coupled between them than to photons, electrons and protons form a single 
            baryonic fluid. We use β̌_e to denote the baryon velocity. To zero order in cosmological perturbations, we find the evolution of the 
            bulk velocity β_e
            
    β̇_e^i+ℋβ_e^i 
    		            = 4ρ_γ/3ρ_ban_eσ_TΔβ^i .

		    To first order in cosmological perturbations and β, the evolution of the first two
		    moments of the distribution is
            
    δ̇_b   +2ℋβ̌_e·δv̌_b +∂_iδ v^i_b-β_e^i∂_iA 
    		                -δ^ijC_ij-D_ijkδ^ijβ_e^k
       = 4ρ_γ/3ρ_ban_eσ_T[β̌_e·(δv̌_γ
    		                -δv̌_b) + δv̌_b·Δβ̌] ,
    		                
    δv̇_b^i     +ℋδ v_b^i
    		                +(β_e^iδ^j_k+β_e^jδ^i_k)∂_jδ v_b^k 
    		                -1/2∂^iA -(δ^jkβ_e^i
    		                +δ^ijβ_e^k)C_jk
       = 4ρ_γ/3ρ_ban_eσ_T[δ v_γ^i-δ v^i_b 
    		                +β^iδ n_e/n_e-β^i_eδ_γ - 3/4β_e jπ_γ^ij -1/2Δβ^iA] .
		    
		    
        

  §.§.§ Massless neutrinos

            Since we will neglect both the mass and coupling of neutrinos, they only free-stream
	        with the same left-hand side as photons. The equation for β_ν is
		    
    β̇^i_ν = 0 .

	        The equation for the evolution of the reduced phase-space density is
	        
    ℱ̇_ν + n^i∂_iℱ_ν 
    		            -4(1/2n^i∂_iA + n^in^jC_ij)(1+5n̂·β̌_ν)
    		            +4β_ν^i(1/2∂_iA+ n^jC_ij+n^jn^kD_ijk)
    		            =0 .
		    
		    
        

  §.§.§ Cold dark matter

            Cold dark matter behaves as collisionless non-relativistic matter, i.e. it just follows
	        the same equations as baryons without interactions. The equation for β_c is
	        
    β̇_c^i+ℋβ_c^i = 0 .

		    The relevant equations for the perturbations are
	        
    δ̇_c+2ℋβ̌_c·δv̌_c +∂_iδ v^i_c
    		            -β_c^i∂_iA -δ^ijC_ij-D_ijkδ^ijβ_c^k    = 0 ,
    		            
    δv̇_c^i +ℋδ v_c^i
    		            +(β_c^iδ^j_k+β_c^jδ^i_k)∂_jδ v_c^k
    		            -1/2∂^iA-(δ^jkβ_c^i
    		            +δ^ijβ_c^k)C_jk   = 0 .
		    
        
        

  §.§.§ Total fluid

            The total energy-momentum tensor, adding all the components, does not contain any explicit
            β contribution after enforcing the cosmic center of mass condition. The conservation
            of the total energy-momentum tensor, which is a direct consequence of the Einstein equations, 
            
    ∇_μ T^μ_ν = 0 ,

            gives us the conservation and Euler equations for the total fluid
            
    δ̇ + 3ℋ(c^2_s-w)δ + (1+w)∂_i(δ v^i-B^i)
                        + 1/2(1+w)Ḣ^i_i   = 0 ,
    δv̇_i+ ℋ(1-3w)δ v_i + ẇ/1+wδ v_i
                        + 1/1+w∂_i(c^2_sδ) 
                        + 1/1+w∂_jπ^j_i-1/2∂_i A    = 0 ,

            where the different variables are defined in the same way as for the individual components,
            but using the total energy-momentum tensor, e.g.
            
    δ   ≡1/ρ∑_sρ_sδ_s , 
    
                    c^2_sδ   ≡1/ρδ P ≡1/ρ∑_s δ P_s .




§ MULTIPOLE ANALYSIS

    The cosmological perturbations can be classified according to their transformation rules
    under the group of spatial rotations. This yields the usual splitting in scalar, vector and
    tensor perturbations, the only ones that contribute to the Einstein equations. Additionally, in the
    standard cosmological perturbation theory, the three types of perturbations are decoupled at the 
    linear level, a fact known as decomposition theorem <cit.>. Since, in this case, 
    the only angular contribution comes from factors of the form (n̂·k̂), i.e. the 
    angle between the line of sight and the direction of the Fourier mode, it is customary to write a 
    multipole expansion for the scalar part of the Boltzmann equation (<ref>) in terms 
    of Legendre polynomials <cit.>
    
    ℱ(τ, ǩ, n̂) = ∑_ℓ=0^∞ (-)^ℓ (2ℓ + 1)
                ℱ_ℓ(τ, ǩ)P_ℓ(n̂·k̂) .

    The Boltzmann equation then unfolds into a whole hierarchy of coupled differential equations for the
    coefficients ℱ_ℓ. The vector modes are usually neglected altogether since, even if
    initially present, they rapidly decay. The tensor modes are predicted in small quantities in many
    inflationary scenarios and their evolution must be followed when studying polarization effects.

    
    The key difference in our scenario is the existence of a new direction β̂, introducing
    new angular dependencies in the Boltzmann equation (<ref>). This means that we must
    resort to a full decomposition in terms of spherical harmonics of the form
    
    ℱ(τ, ǩ, n̂) = √(4π)∑_ℓ=0^∞∑^ℓ_m=-ℓ 
                (-)^ℓ+m√(2ℓ + 1) ℱ^m_ℓ(τ, ǩ)Y^m_ℓ(n̂) ,

    where the coefficients have been defined to match the previous ones for the scalar m=0 modes.
    In particular, in our case, the decomposition theorem no longer holds. As it can be checked, in 
    addition to the usual coupling between ℓ-1 and ℓ+1 modes, the term 
    (n̂·β̂) introduces new couplings between the modes m-1 and m+1. It is 
    possible to write down the new hierarchy of coupled differential equations for the modified 
    Boltzmann equation (<ref>) and it is important for a correct computation of
    CMB anisotropies <cit.>. However, in this work we are mainly concerned with LSS 
    observables and we do not need to trace the evolution of ultrarelativistic species with high
    accuracy.     
    Therefore, we will be working under an approximation scheme that allows us to consider only a subset of 
    these equations. Consequently, we present in this section a self-contained simplified derivation 
    of the system of equations that we will solve in later sections, bypassing the full 
    multipole decomposition.
    
    

 §.§ Scalar-vector-tensor decomposition

        Any spatial vector, in particular the velocity, can be decomposed into a divergence and a 
        divergenceless part
        
    δ v_i = ∂_i v^S + χ_i ,    ∂_iχ^i = 0 ,

        where v^S is the scalar part of the velocity and χ̌ is the vector part, the
        vorticity. In Fourier space, it can be written as
        
    δ v_i = -k̂^i/kθ + χ_i ,    θ≡ -k^2v^S .

        A spatial traceless tensor can be decomposed in a similar way
        
    π_ij = (∂_i∂_j - 1/3δ_ij∂^k∂_k)π^S
                    + 2∂_(iπ^V_j) + π^T_ij ,

        where again the vector part π̌^V is divergenceless and π^T_ij is the
        tensor part, satisfying
        
    ∂^iπ_ij^T=0 ,    δ^ijπ_ij^T=0 .

        Alternatively, we can write it in Fourier space as
        
    π_ij = -2(k̂_ik̂_j-1/3δ_ij)σ
                    + k(k̂_iπ^V_j+k̂_jπ^V_i) + π^T_ij ,
                    σ≡k^2/2π^S ,

        according to the notation of <cit.> for the scalar part of the shear tensor, σ. 
        Adapting the notation of <cit.> for a generic gauge, the metric perturbations 
        can be decomposed as
        
    s^2 = a^2(τ){-   (1+2ψ)τ^2 + 2(∂_i B-S_i) x^iτ
       +(δ_ij-2ϕδ_ij+2∂_i∂_jE + (∂_iF_j+∂_jF_i)
                 +h_ij) x^i x^j} ,

        where Š and F̌ are vector perturbations, i.e. divergenceless vectors, and
        h_ij is a tensor perturbation, i.e. a divergenceless and traceless tensor. Our previously
        defined variables for a general metric perturbation (<ref>) now take the form
        
        
    A    = -2ψ ,
    
    	        B_i    = ∂_i B -S_i ,
    
    	        H_ij   = -2ϕδ_ij + 2∂_i∂_jE+∂_iF_j +∂_jF_i + h_ij .

        
        From now on it will be convenient to work in Fourier space and to choose a basis adapted
        to the previous decomposition. The components of the line-of-sight vector n̂ are
        
        
    n̂   = sinθcosϕ x̂ + sinθsinϕ ŷ+cosθ ẑ ,
       = 1/√(2)^ϕsinθ ê_+ 
                            +1/√(2)^-ϕsinθ ê_-
                            +cosθ k̂ ,
       = √(4π/3)(-Y^+1_1ê_++Y^-1_1ê_-+Y^0_1k̂) ,

        
        where we have chosen the so-called helicity basis <cit.>
        
        
    k̂   ≡ẑ ,
    ê_+    ≡1/√(2)(x̂-ŷ) ,
    ê_-    ≡1/√(2)(x̂+ŷ) ,

        
        and our convention for the spherical harmonics matches those of <cit.> or
        <cit.>.
    
    

 §.§ Lower moments evolution

        The first moments of the Boltzmann equation have already been obtained in the previous sections,
        with ℓ=0 corresponding to the density (<ref>) and ℓ=1 to the
        velocity (<ref>). The next moment ℓ=2 can be obtained via direct integration 
        of the Boltzmann equation (<ref>) and corresponds to the shear tensor
        (<ref>). Performing the appropiate integral, we have
        
    π̇_ij   + ∂_k∫Ω/4πn^kn_in_jℱ_γ
    	                -4/9δ_ij∂_kδ v^k_γ 
    	                -4/3(β_(i∂_j)-1/3δ_ijβ^k∂_k)A 
    	                -8/15(C_(ij)-1/3δ_ijδ^klC_kl)
           +8/15β^m(D_m(ij)-1/3δ_ijD_mklδ^kl)
       = -an_eσ_T[π_ij-β^k_e∫Ω/4πn_kn_in_jℱ_γ 
    	                +4/9β̌_e·δv̌_γδ_ij 
    	                -32/15(δ v_b (iβ_e  j)-1/3δ_ijδv̌_b·β̌_e)
                   -8/15(δ v_b (iΔβ_j)
    	                -1/3δ_ijδv̌_b·Δβ̌)] .

	    To obtain the scalar, vector and tensor parts of this expression, we compute the components as
	    
    π_33≡k̂^ik̂^jπ_ij ,    π_3+≡k̂^iê_+^jπ_ij ,    π_++≡ê_+^iê_+^jπ_ij .

	    The projection of a vector V̌ in the helicity basis is computed in a similar way
	    
    V_3≡k̂·V̌,    
    	        V_+≡ê_+·V̌,     V_-≡ê_-·V̌ .

	    Projecting the equations of motion we obtain
	    
    π̇_33   +  k∫Ω/4π(n̂·k̂)^3ℱ_γ 
    	            - 4/9θ_γ -8/9β̌·ǩ A 
    	            - 8/15(C_33-1/3C^k_k)
    	            + 8/15(β^mD_m33-1/3β^mD_mklδ^kl)
       = -an_eσ_T[π_33 -∫Ω/4π
    	            (β̌_e·n̂)(k̂·n̂)^2ℱ_γ +4/9β̌_e·δv̌_γ - 32/15(δ v_b^3β_e^3-1/3δv̌_b·β̌_e)
                   -8/15(δ v_b^3Δβ_3 -1/3δv̌_b·Δβ̌)] ,
    π̇_3+   + k∫Ω/4π(n̂·k̂)^2(n̂·ê_+)
    	            ℱ_γ - 2/3β_+kA-4/15(C_3++C_+3)
    	            +4/15β^j(D_j3++D_j+3)
        = -an_eσ_T[π_3+-∫Ω/4π(β̌_e·n̂)
    	            (n̂·k̂)(n̂·ê_+)ℱ_γ 
    	            -16/15(δ v_b^3β_e^++δ v_b^+β_e^3)
    	            -8/15(δ v_b^3Δβ^++δ v_b^+Δβ^3)] ,
    π̇_++   + k∫Ω/4π(n̂·k̂)(n̂·ê_+)^2
                    ℱ_γ
    	            -8/15C_++ + 8/15β^mD_m++
       = -an_eσ_T[π_++ 
    	            -∫Ω/4π(β̌_e·n̂)(n̂·ê_+)^2ℱ_γ
    	            -32/15δ v_b^+β_e^+ - 8/15δ v_b^+Δβ_+] ,

	    where it is easy to check that the projections correspond to the scalar part and one of the
	    two vector and tensor helicities
	    
    π_33=-4/3σ ,    π_3+ =  kπ_+^V ,     π_++=π_++^T .
	    
	    There remains to perform a couple of angular integrals, writing down the appropiate coefficients
	    of the expansion (<ref>) and to substitute the metric variables defined in the
	    previous section (<ref>). After these simplifications, and rearranging terms, 
	    the first moments of the Boltzmann equation for photons are

	    
	    

  
*Scalar. (m=0, ℓ = 0,1,2)
	        
    δ̇_γ   +4/3θ_γ +4k^2/3(B-Ė)-4ϕ̇
    	                +8/3 (β̌·ǩ)(ψ-ϕ)+2k^2/3β̌·F̌
       = -4/3an_eσ_T[β̌_e·(χ̌_γ-χ̌_b)
    		            +χ̌_b·Δβ̌-/kk̂·(
    		            β̌_e (θ_γ-θ_b)
    		            +θ_b Δβ̌)] ,
            		        
    θ̇_γ   -k^2/4(δ_γ - 4σ_γ) -k^2ψ 
    	                - 4 (β̌·ǩ)ϕ̇ + 2 k^2(β̌·ǩ)(B-Ė)
    	                -k^2β̌·(Š+1/2F̌)
       = -an_eσ_T[θ_γ - θ_b - ǩ·(β̌_e(δ_γ
    	                -σ_γ) - β̌δ n_e/n_e) 
                        - 3 k/4β̌·π̌^V_γ + ψΔβ̌·ǩ] ,
                        
    σ̇_γ   - 4/15θ_γ + 3k/10ℱ^0_3
    	                -4/15(β̌·ǩ)(ϕ+5ψ)
    	                -4/15k^2(B-Ė)-2/15k^2β̌·F̌
       = -an_eσ_T[σ_γ-4/15kk̂·(β̌_eθ_γ 
    	                + 4β̌_eθ_b
    	                +Δβ̌ θ_b) + 3/10((β̌_e·k̂)ℱ_3^0
    	                -√(2/3)((β̌_e·ê_+)ℱ_3^-1
    	                +(β̌_e·ê_-)ℱ_3^+1))
                    -2/15(β̌·χ̌_γ + 4β̌_e·χ̌_b
    	                +Δβ̌·χ̌_b)] .

        
        

  
*Vector. (m=+1, ℓ = 1,2)
            
    χ̇^+_γ   -3k^2/4π_+^V - 4β_+ϕ̇ + β_+k^2(B-Ė)
    	                + /2β̌·ǩ Ḟ_+ + 1/2β_-ḣ_++
       = -an_eσ_T[χ^+_γ-χ^+_b + β_+δ n_e/n_e
    	                -β_e^+δ_γ - 1/2β_e^+σ_γ
    	                -3/4(β̌_e·k̂)π_+^V - 3/4β_e^-π_++
    	                +ψΔβ_+] ,
    	                
    π̇^V_+     +4/15(χ_γ^++√(3/2)ℱ^+1_3)
    	                +4/15β_+(ϕ+5ψ) + 4/15(S_++Ḟ_+)
    	                -2/15β_-h_++ + 2/15β_3kF_+
       = -an_eσ_T[π_+^V + 4/15k^2β_e^+(θ_γ+kℱ^0_3)
    	                +1/k√(2/15)β_e^-ℱ^+2_3 + 4/15kβ_e^3
    	                (χ_γ^++√(3/2)ℱ^+1_3)
                   +4/15k^2(Δβ_+θ_b+4β_e^+θ_b)
                        +4/15k(Δβ_3χ_b^+ + 4β_e^3χ_b^+)]

            
        

  
*Tensor. (m=+2, ℓ = 2)
            
    π̇_++   +k√(2/15)ℱ^+2_3 +4/15ḣ_++ 
    	                + 4/15(β̌·ǩ)h_++
    	            =-an_eσ_T[π_++ - 8/15(β_e^+χ_γ^+ 
    	                +4β_e^+χ_b^+ + Δβ_+χ_b^+)] .

        The corresponding results for the other helicity can be obtained substituting -↔ +
        in every sub and superscript.

        
        Several comments are in order now. In the first place, note that all the couplings between
        scalar, vector and tensor modes are introduced by terms proportional to β. In the 
        standard case, each mode evolves independently. In the second place, note the appearance
        of terms with ℓ = 3. What we present here are but the lowest moments of a whole 
        hierarchy of coupled differential equations. This system obeys a recurrence relation but
        it must be truncated at a finite, prefearably large, value of ℓ. The traditional 
        line-of-sight approach <cit.> was developed to allow a truncation at 
        lower ℓ, and it is used by every modern Boltzmann solver <cit.>. It will not be used here, 
        even though it can be adapted to our case <cit.>, since we are mainly interested in 
        the matter power spectra and not in the CMB. Instead, we will study a simplified version of the 
        system, under the following approximations.        
        
            
  * We will work to first order in β and to first order in cosmological perturbations, 
                keeping cross-products. We have been implicitly working under this assumption, since the 
                RW background is only correct to first order in β, but we will consistently carry 
                it through. 
                
            
  * We will assume that there are no initial vector or tensor modes to zero order in β. 
                The assumption is justified for vector modes, since most popular models of inflation do 
                not produce them at all. On the other hand, since we have not detected tensor modes so 
                far, and we have stringent limits on their amplitude, we assume that their amplitude is 
                small enough so we can neglect them.
                
                Under this assumption, the hierarchy is simplified. Since they are zero initially
                the only production occurs through their new couplings, i.e. it is proportional to 
                β. The vector modes are then β and the tensor modes are β^2.
                In general, we can neglect the backreaction of higher m modes into lower m modes.
                As we will see in the next section, the Einstein equations are not modified, so we
                can apply the same reasoning to the metric variables, i.e. Š and F̌ are 
                β and h_ij is β^2.
            
            
  * The last approximation is the so-called fluid approximation. We will truncate the
                hierarchy at ℓ = 2 for scalar and vector modes. This is a classic
                working assumption in approximate computations of the CMB, that may introduce up to
                10% errors in computations of the CMB spectrum <cit.>. Another source of
                error in our case is the naive truncation scheme we are using, i.e. setting to zero
                all higher moments. It is a well known fact that this truncation scheme, neglecting
                the damping produced by the transfer of power to higher moments, produces some spurious
                growth at small scales in the ultrarelativistic species. Since, in this work, we are 
                not interested in following with great precision the evolution of photons or neutrinos,
                and we have checked that they have a negligible impact in our final results, 
                we will nonetheless stick to this crude truncation scheme in our numerical solutions.
                
                To increase the accuracy of the results for ultrarelativistic species one would need to 
                evolve the full hierarchy <cit.>, or at least introduce a better truncation scheme 
                <cit.> or an effective viscosity in the equations of motion for 
                ultrarelativistic species <cit.>.
            
        
        


§ EINSTEIN EQUATIONS

    In this section we present the last piece of information needed to solve the system: the evolution
    of the metric perturbations. It is worth remembering that, to first order in β, the background
    quantities we are interested in, e.g. ρ and P, are equal in the 𝒪 and 
    𝒪̃ frames. Hence, as in previous sections, we will drop the distinction. 
    The full energy-momentum tensor for each component is
    
    
    T^0_0 + δ T^0_0   = -ρ-δρ-(ρ+P)B_iβ^i ,
    
            T^0_i + δ T^0_i   = δ Q_i+(ρ+P)(δ^j_i
                +1/2Aδ^j_i +1/2H^j_i)β_j ,
    
            T^i_j + δ T^i_j   = Pδ^i_j + δ P δ^i_j + δΠ^i_j 
                +(ρ + P)β_jB^i ,

    
    where β is different for each component. We must write now the Einstein equations with this
    source for the metric (<ref>). For the background evolution we obtain the
    standard Friedmann equations plus a condition for the cosmic center of mass frame
    
    ℋ^2    = 8π Ga^2/3∑_sρ_s , 
    
            0    = ∑β̌_s(ρ_s + P_s) ,
    ℋ̇+1/2ℋ^2    = -4π Ga^2∑_s P_s .

    Applying the condition (<ref>), the explicit β contributions to the
    full energy-momentum tensor vanish
    
    
    δ T^0_0   = ∑_sδ T_s^0_0 = -∑_sδρ_s ,
    δ T^0_i   = ∑_sδ T_s^0_i = ∑_sδ Q_s i ,
    δ T^i_j   = ∑_sδ T_s^i_j = ∑_s(δ P_sδ^i_j 
                +δΠ_s^i_j) .

    
    In our case the non-relativistic species have pressure and anisotropic stress of order 
    β. Splitting the sources into relativistic and non-relativistic components we have
    
    δρ   = ∑_NRδρ + ∑_Rδρ ,
    δ P    = 2/3∑_NRρ β_kδ v^k + 1/3∑_Rδρ ,
    δ Q^i    = ∑_NRρ δ v^i + 4/3∑_Rρ δ v^i ,
    δΠ^ij   = ∑_NRρ(β^iδ v^j + β^jδ v^i
                -2/3δ^ijβ_kδ v^k)+∑_RδΠ^ij .

    Finally, the Einstein equations read
    
        
  * (0,0)
            
    2k^2ϕ+6ℋ(ϕ̇+ℋψ) 
                        -2k^2ℋ(B-Ė) =-8π G a^2δρ .

        
  * (0,i)
            
    k(ϕ̇+ℋψ)    = -4π Ga^2δ Q_3 ,
    
                    k^2(S_++Ḟ_+)    = 16π Ga^2δ Q_+ .

        
  * (i,j)
            
    k^2(ϕ-ψ)-k^2(∂_τ +2ℋ)(B-Ė)
                           =-12 π Ga^2δΠ_33 , 
    
                    (∂_τ+2ℋ)(ϕ̇+ℋψ)
                        +ψ (ℋ̇-ℋ^2)
                           =4π Ga^2(δ P+δΠ_3 3) ,
     k(∂_τ +ℋ)(S_++Ḟ_+) 
                           = 16π Ga^2δΠ_+3 ,
    1/2(∂^2_τ +2ℋ∂_τ + k^2)h_++   = 8π Ga^2δΠ_++ .

    
    Again, the results for the - helicity can be obtained substituting -↔ +
    in every sub and superscript.
    With our notation, the Newtonian gauge can be obtained just setting B=E=0 and the synchronous
    gauge is defined as
    
    
    ψ   = B = 0 ,
    ϕ   = η ,
    
            E    = -1/2k^2(h+6η) .

    
    The Einstein equations in the synchronous gauge for scalar perturbations can be written as
    
    ḣ -2 k^2η/ℋ = 3 ℋδ ,
       η̇ = 3ℋ^2/2 k^2(1+w)θ ,
       ḧ + 6 η̈ + 2ℋ(ḣ +6 η̇) - 2 k^2η  
                = - 12 ℋ^2σ ,
       ḧ+ℋḣ = - 3 (1 + 3 c_s^2) ℋ^2δ ,

    where we have defined
    
    c_s^2δ≡1/ρδ P ,    θ≡1/ k(ρ + P)δ Q_3 ,
                    σ≡ -4/3ρδΠ_33 .

    Appendix <ref> contains a discussion about how the gauge transformations
    are modified for non-comoving fluids and, in particular, how to relate the Newtonian and synchronous 
    gauges.



§ REDUCED SYSTEM AND FINAL EQUATIONS

    This section contains the final, simplified equations that will be numerically integrated.
    In the first place, the relevant equations for the evolution of the bulk velocities are presented
    in section <ref>.
    The background follows the standard ΛCDM evolution to first order in β, but there is
    a first order effect on the perturbations.

    
    Once the evolution of the bulk velocities is known, we need to study the modified evolution of
    the perturbations. Working to first order in β and with the approximations made at the end 
    of section <ref>, our modifications to the scalar and vector modes decouple 
    and can be treated separately. This is the subject of the last two sections <ref>
    and <ref>. As mentioned before, the tensor modes are second order in β and
    we neglect them.
    
    

 §.§ Bulk velocities

        The evolution of the velocities of the different fluids is governed by
        
        
    β̇̌̇   = -1/τ_cΔβ̌ ,
    β̇̌̇_ν   = 0 ,
    β̇̌̇_e    = -ℋβ̌_e + 1/Rτ_cΔβ̌ ,
                     
    β̇̌̇_c    = -ℋβ̌_c .

        
        where we have defined
        
            R   ≡3ρ_b/4ρ_γ , 
                       Δβ̌    ≡β̌_γ-β̌_e ,

            τ^-1_c    ≡an_eσ_T ,        𝒜   ≡R/1+R ,
        
        and, as we mentioned before, the initial conditions are chosen according to the constraint
        
    ∑_s(ρ_s+P_s)β̌_s = 0 ,

        so the cosmic center of mass condition is maintained in the evolution. 
        Moreover, we will assume that all the bulk velocities (β̌,β̌_e, β̌_ν,β̌_c) 
        are aligned along the β̂ axis in the 𝒪 frame. As we will shortly see,
        and can be inferred from (<ref>) and (<ref>), when two species
        are tightly coupled their velocities evolve to become equal. Once a particle species decouples,
        the magnitude of its velocity evolves independently but, in the absence of additional interactions
        or other sources of anisotropy, it does not change its direction. In our scenario, we 
        assume that the whole visible sector has been in thermal equilibrium at some time so all its
        components, even if they are decoupled like the neutrinos, have velocities pointing in 
        the direction β̂. The only remaining contribution is the dark sector, with DM among its
        components. The dark sector in the 𝒪 frame counterbalance the flux of momentum of 
        the visible sector to achieve an isotropic universe, so it must point in the -β̂ direction.

        
        CDM and neutrinos are decoupled, but the photon-baryon system must be treated with
        some care. In the tight coupling limit, τ_c≪ 1, it is easy to see that the velocities 
        converge in direction and magnitude and we can look for an approximate solution of this system.
        Expanding perturbatively in the small parameter τ_c we have
        
    β̇̌̇   = -𝒜ℋβ̌ + τ_c ,
    Δβ̌   = 𝒜ℋβ̌τ_c + τ_c^2 .

        The differential equation can be solved to yield
        
    β̌ = β̌_0/1+R ,

        where β̌_0 is the initial velocity of the visible sector in the 𝒪 frame, the 
        only additional free parameter in our model. In a similar way, the neutrino and CDM equations can 
        be solved to give
        
    β̌_ν   = β̌_0 ,
    β̌_c    = β̌_c^todaya^-1 .

        Using the scaling (<ref>) and enforcing the constraint 
        (<ref>) during the tightly coupled regime, the value of 
        β_c^today is found to be
        
    β_c^today = -4/3β_0Ω_γ + Ω_ν/Ω_cdm .

        It is important to notice that according to the evolution of β_c (<ref>),
        early enough in time, the condition β_c≪ 1 could break down. However, this is only
        the case if the DM keeps the non-relativistic distribution at early times. 
        If the DM were light enough it could behave as a radiation-like fluid well before its bulk
        velocity reaches β_c=1. In this case, β_c would remain constant and small.
        On the other hand, if the DM is heavy, its bulk velocity can reach the relativistic 
        regime. In this case, it is worth mentioning that even if the anisotropies that would arise 
        at the background level are β_0 and not β_0^2, their effects can only 
        be relevant well before the matter-dominated era, with no observational consequences.


        In order to avoid choosing any particular framework for the dark sector, we will not follow
        the dark matter evolution using (<ref>). Instead, we follow
        the evolution of photons (<ref>), neutrinos (<ref>) and baryons
        (<ref>), and then the momentum of the dark sector, regardless of its composition, can
        be obtained imposing the center of mass condition (<ref>). Thus, our only 
        assumptions regarding the dark sector are that it is subdominant at early times and that it 
        behaves as cold matter plus cosmological constant at late times.
            
    

 §.§ Scalar modes
        
        Here we provide the scalar mode equations for the photon-baryon system. Neutrinos and CDM obey the
        same equations as photons and baryons, respectively, without the coupling term, i.e. σ_T=0.
        In addition to the definitions (<ref>) we define
        
    Δθ≡θ_γ - θ_b ,
        
        The evolution of the photon perturbations is described by the equations (<ref>,
        <ref>, <ref>). Under the approximation scheme discussed
        at the end of section <ref>, i.e. neglecting backreaction of vector modes and
        neglecting moments higher than ℓ=2, these equations take the form
        
        
    δ̇_γ+4/3θ_γ +4k^2/3(B-Ė)-4ϕ̇
    	                +8/3 (β̌·ǩ)(ψ-ϕ)
    	           = -4/3τ_c[-/kk̂·(
    	            β̌_eΔθ
    	            +θ_b Δβ̌)] ,
    θ̇_γ -k^2/4(δ_γ - 4σ_γ) -k^2ψ 
    	                - 4 (β̌·ǩ)ϕ̇ + 2 k^2(β̌·ǩ)(B-Ė)
                   = -1/τ_c[Δθ - ǩ·(β̌_e(δ_γ
                    -σ_γ) - β̌δ n_e/n_e-Δβ̌ψ)] ,
    σ̇_γ - 4/15θ_γ -4/15k^2(B-Ė)
    	                -4/15(β̌·ǩ)(ϕ+5ψ)
                   = -1/τ_c[σ_γ-4/15kk̂·(β̌_eθ_γ 
                    + 4β̌_eθ_b
                    +Δβ̌ θ_b)] .

        
        Neutrinos are described by the same system, without collision term.
        For baryons, we apply our approximation scheme to equations (<ref>) and 
        (<ref>) obtaining
        
        
    δ̇_b + θ_b -3ϕ̇ +k^2(B-Ė) 
                    -2ℋ(β̌_e·k̂)/kθ_b 
                    +2(β̌_e·ǩ)(ψ-ϕ)
                   = 1/Rτ_c[-/kk̂·(
    	            β̌_eΔθ
    	            +θ_b Δβ̌)] ,
    θ̇_b + ℋθ_b - k^2ψ + 2(β̌_e·ǩ)θ_b
    	            -4(β̌_e·ǩ)ϕ̇ + 2 k^2 (β̌_e·ǩ)(B-Ė)
    	           = 1/Rτ_c[Δθ - ǩ·(β̌_e(δ_γ
                    -σ_γ) - β̌δ n_e/n_e-Δβ̌ψ)] .

        
        Again, CDM equations take the same form, but without collision term.
        The evolution of the full energy-momentum tensor is described by
        
        
    δ̇ + 3 ℋ(c_s^2-w)δ + (1+w) θ - (1+w)(3ϕ̇-k^2(B-Ė))   = 0 ,
    θ̇ + (1- 3 w) ℋθ + ẇ/1+wθ -k^2/1+wc_s^2δ +4 k^2/3 (1+w)σ -k^2ψ   = 0 ,

        
        where
        
    c_s^2δ   = 1/ρ(1/3ρ_γδ_γ
                    +1/3ρ_νδ_ν - 2ρ_cβ_cθ_c/3k
                    -2ρ_bβ_eθ_b/3k) ,
    σ   = 1/ρ(ρ_γσ_γ + ρ_νσ_ν 
                    + ρ_cβ_cθ_c/k + ρ_bβ_eθ_b/k) .

        
        The system of cosmological perturbations is a system of linear differential equations. A generic
        cosmological perturbation g(τ,ǩ) can be written as a product of a primordial
        perturbation, encoding the initial condition, and a transfer function, encoding the subsequent
        evolution <cit.>
        
    g(τ,ǩ) = T_g(τ,ǩ)ℛ_ǩ .

        As long as we consider only the adiabatic mode, every cosmological perturbation is proportional
        to the primordial curvature perturbation ℛ_ǩ. Furthermore, since the system 
        is linear, it can be recast into a system for the evolution of the transfer functions with the
        substitution g→ T_g. It is common practice to abuse slightly of the notation and to denote
        T_g as g, the perturbation itself, and to solve the system as if it had initial conditions
        ℛ_ǩ=1. The information about the initial conditions is recovered later, in 
        the computation of the physical spectra, convolving the transfer function with the primordial
        spectrum. We follow this practice.

        
        An important difference in (<ref>) and (<ref>) 
        with respect to the β=0 case is the appearance of imaginary terms. Usually, even though 
        the Fourier coefficients are generally complex, the evolution equations are real. In this case,
        both real and imaginary part of T_g satisfy the same equation and with a judicious choice
        of initial global phase it can be made purely real. With the appearance of complex coefficients,
        real and imaginary parts form a coupled system with different equations of motion.

        
        We will assume that the global phase has been chosen so that the imaginary parts are initially
        zero, or at most β. Even if they are initially zero, the terms proportional to 
        β couple the imaginary to the real parts, driving them to a finite value proportional to 
        β. Then we are in the same situation as with the vector modes. The imaginary parts of the 
        scalar modes are determined by the real parts, but they do not backreact on them. The real parts 
        follow the standard cosmological evolution. Therefore, every perturbation transfer function can 
        be splitted as
        
    g(τ,ǩ) = g^R(τ,k) + (β̂·k̂)g^I(τ,k) ,

        where now g^R and g^I are purely real and do not depend on the direction of k̂. 
        Following this definition and the previous discussion, we are in the following situation.
        
            
  * The real part of the perturbations g^R only contains adiabatic perturbations, 
                as in standard ΛCDM, and follows the standard evolution.
            
  * The real parts act as external sources in the system for the imaginary parts,
                via contributions β. The imaginary parts of the perturbations 
                g^I are β.
        
        
        We will work in the synchronous gauge (<ref>). The modified 
        gauge-transformation properties are provided in the appendix <ref>.
        Written in the synchronous gauge, the final system for the photon perturbations is
        
        
    δ̇^I_γ + 4/3θ_γ^I + 2/3ḣ^I
                    -8β k/3η^R
    	           = 4/3kτ_c[β_eΔθ^R+θ_b^R Δβ] ,
    θ̇_γ^I -k^2/4(δ_γ^I - 4σ_γ^I)
                    +β k(ḣ^R+2η̇^R)
                   = -1/τ_c[Δθ^I - k(β_e(δ_γ^R
                    -σ_γ^R) - βδ n_e^R/n_e)] ,
    σ̇^I_γ - 4/15θ_γ^I
                    -2/15(ḣ^I+6η̇^I) - 4β k/15η^R
                   = -1/τ_c[σ_γ^I-4/15k(β_eθ_γ^R 
                    + 4β_eθ_b^R
                    +Δβ θ_b^R)] .

        
        The equations for neutrinos are the same, setting the collision part to zero. For the baryons 
        we have
        
        
    δ̇_b^I + θ_b^I +1/2ḣ^I 
                    -2ℋβ_e/kθ_b^R -2β_ekη^R              
                   = -1/kRτ_c[β_eΔθ^R+θ_b^RΔβ] ,
    θ̇^I_b + ℋθ_b^I + 2β_e kθ_b^R
                    +β_ek(ḣ^R+2η̇^R)
    	           = 1/Rτ_c[Δθ^I - k(β_e(δ_γ^R
                    -σ^R_γ) - βδ n_e^R/n_e)] .

                
        And finally for the full fluid
        
        
    δ̇^I + 3 ℋ(c_s^2-w)δ^I + (1+w) θ^I    = -1/2(1+w) ḣ^I ,
    θ̇^I + (1- 3 w) ℋθ^I + ẇ/1+wθ^I -k^2/1+wc_s^2δ^I +4 k^2/3 (1+w)σ^I    = 0 ,

                
        where, in the synchronous gauge,
        
    c_s^2δ^I    = 1/ρ(1/3ρ_γδ_γ^I
                    +1/3ρ_νδ_ν^I
                    -2ρ_bβ_eθ_b^R/3k) ,
    σ^I    = 1/ρ(ρ_γσ_γ^I + ρ_νσ_ν^I
                    + ρ_bβ_eθ_b^R/k) .

        To complete the system, we compute the variables η and aḣ using a combination of the 
        Einstein equations
        
    η̇^I    = 3ℋ^2/2 k^2(1+w)θ^I , 
    ḧ^I+ℋḣ^I    = - 3 (1 + 3 c_s^2) ℋ^2δ^I .

        
        Usually, one would integrate the equations for photons, baryons, neutrinos and CDM. Then, after
        adding all the components, one would compute the sources for the metric perturbations, i.e. the
        full energy-momentum tensor. As we mentioned before, in our setup we have found it more 
        convenient to follow a different route. Instead of tracking the behaviour of the dark sector, 
        thus choosing a particular DM framework, we follow the evolution of the whole fluid 
        (<ref>). Using these equations, the only underlying assumptions are
        
            
  * The dark sector is subdominant with respect to neutrinos and photons at early times,
                i.e. before the matter-domination era. 
            
  * There is a transition to a CDM behaviour at late times. 
        
        Under these assumptions, the only CDM contribution to the full fluid goes into the 
        equation of state w, since in the synchronous gauge it does not contribute to
        δ P or σ at first order in β. The evolution of the dark sector can be
        obtained afterwards subtracting the contributions of photons, baryons and neutrinos from the
        full fluid.

                
        We are in position now to numerically solve the system. Our strategy
        can be summarized as follows.
        
            
  * The background and the real part of the perturbation, labelled with R, follow the
                standard evolution and act as external sources in our system.
                We use class <cit.> to precompute these sources and then solve 
                numerically the system for the imaginary parts, labelled with I. The equations for the 
                evolution of β are (<ref>), (<ref>) and (<ref>).
            
  * The system to be solved for the imaginary parts consists of 
                (<ref>), (<ref>), 
                (<ref>), (<ref>), (<ref>) 
                and a neutrino contribution equal to (<ref>) but without 
                collisions.
            
  * The initial conditions are discussed in the appendix <ref>, 
                where we find the analytic super-Hubble behaviour of the perturbations. In addition to
                the usual adiabatic and isocurvature modes <cit.>, we find a “sourced
                mode” that depends on the initial value of η^R, i.e. the standard initial 
                curvature perturbation in the adiabatic mode. If we assume that 
                θ_γ(τ=0)=θ_ν(τ=0), since they have been thermally coupled at
                some point, the presence of these external sources introduces an isocurvature velocity
                perturbation
                
    θ^I_γ(τ=0) = θ^I_ν (τ= 0) = 2β_0 k η^R(τ=0) ,

                where β_0 is the initial value of the velocity of photons and neutrinos.
                All other perturbations are assumed to be initially zero. The whole analysis and the
                series expansions are described in the appendix <ref>.
            
  * The impact of perturbed recombination <cit.> is usually 
                neglected but it will be important in our case. This subject is covered in section 
                <ref>.
            
  * We need two approximation schemes to follow the numerical evolution, the so-called
                tight coupling approximation (TCA) and radiation streaming approximation (RSA). We use
                the same criteria and switches as class <cit.>. The appropiate
                equations for TCA are presented in section <ref>. While TCA is important at
                early times, RSA is crucial at late times, when the ultrarelativistic species are
                decoupled and start oscillating fast. In this case, it is time consuming to follow 
                every oscillation but it is not necessary since at late times the impact of 
                ultrarelativistic components on the metric potentials is negligible. We describe this 
                scheme in section <ref>. 
        
        
        

  §.§.§ Perturbed recombination

            The disturbances in the photon temperature field produce perturbations in the ionization 
            fraction of the electrons. In standard ΛCDM, this inhomogeneous recombination 
            produces second order effects in the CMB, but it has proven important at late times when 
            computing other observables like the 21 cm radiation, through its effects in the gas 
            temperature <cit.>.

            
            Boltzmann codes like camb <cit.> and class <cit.>
            have implemented perturbed recombination at late times. This implementations follow the 
            formulae of recfast <cit.>, including perturbations into the 
            recombination coefficient that effectively takes into account multilevel atom computations. 
            These codes also track the evolution of the gas temperature and its perturbations, that 
            modify significantly the baryon sound speed.

            
            The study of the dark ages in detail is beyond the scope of this work, but, in our case, 
            perturbed recombination plays a role in the photon-baryon system to first order in β, 
            where a perturbation in the number of free electrons δ n_e appears, e.g. see 
            (<ref>). The baryon sound speed has been neglected in our calculations, it 
            is only important at very small scales, and we will neglect the perturbations in the gas 
            temperature and the recombination coefficient as well. Defining the ionization fraction 
            <cit.> as
            
    x_e≡n_e/n_b ,

            where n_e and n_b are the number densities of free electrons and baryons, respectively, 
            we have 
            
    δ n_e/n_e≡δ_b + δ_x_e ,

            where δ_x_e≡δ x_e/x_e is the relative perturbation in the ionization
            fraction. Since δ n_e always appears multiplied by β, to study its evolution
            it suffices to take (<ref>) with β=0. After a few 
            manipulations we obtain
            
    τ(δ n_e/n_e)-τ(1/a^3n_e)
                        a^3δ n_e + θ_b + 1/2ḣ = 0 .

            Substituting the definition (<ref>) we get the final evolution equation
            
    δ̇_x_e + ẋ_e/x_e(δ_b + δ_x_e) = 0 .

            We solve this equation for each mode with initial conditions 
            δ_x_e(τ=τ_ini)=0 and using the full ionization history x_e(τ) 
            provided by the thermodynamics module in class, computed using recfast.
            
        

  §.§.§ Tight coupling expansion

            At early times, the time scale of Thomson scattering, τ_c≡ (an_eσ_T)^-1, is
            much shorter than the time scales of evolution of the background, ℋ^-1, or
            that of the modes, k^-1. In this regime, the photon-baryon system becomes computationally 
            hard to solve, since it involves widely different scales. However, we can find approximate
            expressions to follow the evolution, expanding perturbatively in the small parameter
            τ_c. This is known as the tight coupling approximation (TCA). In this work, as usual, 
            we will restrict ourselves to the lowest order in this expansion. Higher order terms in the 
            standard scenario, and the procedure to obtain them systematically, can be found in 
            <cit.> and <cit.>.

            
            To obtain the equations of motion to leading order in τ_c, we need to expand 
            Δθ^I and σ^I_γ to first order in τ_c, just as we did with the 
            TC expansion for β. Solving simultaneously for both quantities, and plugging in the 
            TC values for the real part of the perturbations and β, we get
            
    Δθ^I    = β k(δ_γ^R - δ_b^R - δ_x_e^R)
                        + τ_c{k^2/4𝒜δ_γ^I 
                        + 𝒜ℋθ_γ^I    
           +β k[𝒜(𝒜-1)ℋ(δ_γ^R-δ_b^R-δ_x_e^R)
                        -𝒜ℋδ_γ^R + 𝒜θ_γ^R - 4/15θ_γ^R
                        + 𝒜δ̇_x_e^R+𝒜/6ḣ^R
                        -2/15(ḣ^R+6η̇^R)]}+τ_c^2 ,
    σ_γ^I    = 4β/3kθ_γ^R 
                        + τ_c/15{4θ_γ^I + (ḣ^I+6η̇^I)
                        + β[k(𝒜-5)δ_γ^R + 4kη^R
                        + 8𝒜ℋ/kθ_γ^R]} +τ_c^2 .

            The final equations of motion of the photon-baryon plasma during the tightly coupled phase
            are
            
    δ̇_b^I+θ_γ^I+1/2ḣ^I    =
                        β[k𝒜/4δ_γ^R 
                        + k(3/4δ_γ^R-δ_b^R-δ_x_e^R)
                        +2kη^R+2𝒜ℋ/kθ_γ^R] ,
    δ̇_γ^I+4/3θ_γ^I+2/3ḣ^I    =
                        β[k𝒜/3δ_γ^R +8k/3η^R
                        +8𝒜ℋ/3kθ_γ^R] ,
    θ̇_γ^I + 𝒜ℋθ_γ^I
                        + k^2/4(𝒜-1)δ_γ^I    =
                        β k[𝒜(1-𝒜)ℋ(δ_γ^R-δ_b^R-δ^R_x_e)
                        -𝒜θ_γ^R -4/3θ_γ^R - 𝒜δ̇^R_x_e
                        -2η̇^R - (1+𝒜/6)ḣ^R] .
        
            The evolution of every mode starts in a radiation-dominated phase during which TCA is valid.
            We start evolving this set of equations for each mode until we switch the TCA off, according
            to the class switches <cit.>. Once we switch it off we evolve the full
            system, joining the solutions smoothly.
      
        

  §.§.§ Radiation streaming approximation

	        The evolution equations for ultrarelativistic species (<ref>), 
	        like neutrinos, can be combined into a single second-order differential equation
	        
    δ̈^I_ν = -k^2/3(δ_ν^I - 4σ_ν^I) 
    	                + 4β_ν k/3(ḣ^R + 4η̇^R) - 2/3ḧ^I ,

	        and the same applies to photons after decoupling. Once the perturbation is sub-Hubble, it
	        starts oscillating very fast. During matter domination, it is possible to ignore the 
	        contribution of ultrarelativistic species to the total density but their contribution to the
	        total velocity is still important. The radiation streaming approximation (RSA) consists on 
	        following only the non-oscillatory particular solution of these equations <cit.>.
	        Since in this regime |δ̈_ν|≪ k^2|δ_ν| and |σ_ν|≪ |δ_ν|, 
	        our RSA solution is
	        
    δ^I_ν   = -2/k^2ḧ^I + 4β_ν/k(ḣ^R+4η̇^R) ,
    θ^I_ν   = -1/2ḣ^I + 2β_ν kη^R ,
    σ^I_ν   = 0 .

	        We will apply the same equations to photons after decoupling, neglecting the small impact of 
	        reionization.
            
    

 §.§ Vector modes

        In this section, we will describe the evolution of the vector modes, following the same steps
        as in the previous section.
        Since the evolution equations, and the initial conditions, for both helicities are the same, 
        we can rewrite the vorticity for the species s as
        
        
    χ̌_s    = χ_s((β̂·ê_+) ê_-
                                + (β̂·ê_-) ê_+)
       = χ_s((β̂·x̂) x̂
                                + (β̂·ŷ) ŷ) .

        
        Then, we do not need to distinguish between helicities and we can just write one equation for 
        χ_s. The same applies to the vector part of the shear tensor π^V_s. Starting
        from (<ref>) and (<ref>), under our approximation
        scheme, i.e. neglecting tensor modes and moments higher than ℓ=2, the evolution of the
        photon vector modes is described by
        
        
    χ̇_γ +1/2β(ḣ-2η̇) -3k^2/4π^V_γ   = -1/τ_c[Δχ
                    +(βδ n_e/n_e-β_e(δ_γ 
                    +1/2σ_γ))] ,
                    
    π̇^V_γ + 4/15χ_γ + 4/15(S+Ḟ)
                    +4/15βη   = -1/τ_c[π^V_γ 
                    + 4/15k^2(β_eθ_γ + 4β_eθ_b + Δβθ_b)] ,

        
        where we have defined
        
    Δχ≡χ_γ - χ_b .

        Again, the behaviour of neutrinos can be obtained from these equations, setting to zero the 
        collision term. From (<ref>) and (<ref>), baryons and dark matter evolve 
        according to
        
    χ̇_b+ℋχ_b +β_eθ_b 
                    + 1/2β_e(ḣ-2η̇) 
                       = 1/τ_cR[Δχ + (βδ n_e/n_e
                    -β_e(δ_γ + 1/2σ_γ))] ,
    χ̇_c+ℋχ_c +1/2β_c(ḣ-2η̇)    = 0 .

        The evolution of the total vorticity, i.e. the vorticity of the full fluid, is
        
    χ̇ +ℋ(1-3w)χ + ẇ/1+wχ -k^2/1+wπ^V = 0 ,

        where
        
    π^V = 1/ρ(ρ_νπ^V_ν+ρ_γπ^V_γ
                    -1/k^2β_eθ_b) .

        Finally, the relevant Einstein equation is
        
    S+Ḟ = 16π Ga^2/k^2(ρ+P)χ .

        Our strategy for the integration of the vector modes can be summarized as follows.
        
            
  * The background and the scalar modes evolve according to standard ΛCDM and act
                as external sources for the vector modes. We use class to precompute the sources.
            
  * The system to be integrated consists of (<ref>), (<ref>),
                (<ref>), (<ref>) and a neutrino contribution equal
                to (<ref>) but without collisions.
            
  * The initial conditions are described in the appendix <ref>.
            
  * As in the scalar case, we have to take into account perturbed recombination and
                we need to implement the TC expansion and RSA in the vector case, as discussed in the
                next sections <ref> and <ref>.
        
        
        

  §.§.§ Tight coupling expansion

            As in the previous section, we need to find the approximate equations to follow the tightly
            coupled phase. Performing the same manipulations, and inserting the TC solutions for β 
            and the scalar modes, we get
            
    Δχ   = β(δ_γ - δ_b -δ_x_e)
                        +τ_c{𝒜ℋχ_γ 
                        +β[𝒜(𝒜-1)ℋ(δ_γ-δ_b-δ_x_e)
                        -𝒜ℋδ_γ+𝒜/3θ_γ + 2/15θ_γ
                        +𝒜δ̇_x_e
                                         
                        + 𝒜/6ḣ + 1/15(ḣ+6η̇)]}
                        +τ_c^2 ,
                        
    π^V_γ   = -4β/3k^2θ_γ -4/15τ_c{χ_γ + (S+Ḟ) - β[5/4δ_γ - η
                        -𝒜/4δ_γ-2𝒜ℋ/k^2θ_γ]} .

            The equation governing the evolution of the photon vorticity during TC is
            
    χ̇_γ = -𝒜ℋχ_γ - βθ_γ - 1/2β(ḣ-2η̇) + β𝒜[(1-𝒜)ℋ(δ_γ
                        -δ_b-δ_x_e)-1/3θ_γ- 1/6ḣ-δ̇_x_e] .

            
        

  §.§.§ Radiation streaming approximation

            The equations (<ref>) and (<ref>) for neutrinos, or photons
            after decoupling, can be combined into a second order differential equation
            
    χ̈_ν = -k^2/5χ_ν - k^2/5(S+Ḟ)-k^2/5β_νη
                        -1/2β_ν(ḧ-2η̈) .

            During the period of rapid oscillations, we are in the situation in which |χ̈_ν|≪ |k^2χ_ν|
            and |χ̇_ν|≪ |k^2π^V_ν|. From the previous equation and (<ref>),
            the approximate non-oscillating particular solution is found to be
            
    χ_ν   = -(S+Ḟ)-β_νη + 5/2k^2β_ν(ḧ-2η̈) ,
    π_ν^V   = 2/3k^2β_ν(ḣ-2η̇) .

        
        

  §.§.§ Semi-analytic solutions
            
            The equations obtained admit semi-analytic solutions in some regimes. During the TC regime,
            from (<ref>) and (<ref>),  we find, to lowest order in τ_c,
            
    (1+R)χ_γ   = -1/2β_0(h-2η) + 
                        β_0 𝒜(δ_γ -δ_b -δ_x_e)
                        -β_0/3∫(1+𝒜/4)θ_γτ +𝒞_γ ,
    
                    (1+R)χ_b    = -1/2β_0(h-2η)
                        -β_0/1+R(δ_γ -δ_b -δ_x_e)
                        -β_0/3∫(1+𝒜/4)θ_γτ +𝒞_γ ,

            where 𝒞_γ is a constant of integration, to be set with the initial
            condition, and β_0 is the initial velocity of the photon-baryon plasma.
            Another result that can be obtained, integrating (<ref>), is the evolution
            of CDM
            
    aχ_c    = -1/2β_c^today(h-2η)+𝒞_c .

            As discussed before, we do not specify the behaviour of the dark sector at early times. Hence,
            we do not use this equation. We solve the system instead using the total vorticity
            (<ref>) and then we obtain the vorticity of the dark sector subtracting the 
            other components.
        



§ RESULTS AND DISCUSSION

    

 §.§ Observables

        Once we have constructed a consistent system of equations, we must discuss which of the 
        intermediate variables correspond to physical observables. One of the main observables in
        cosmology is the distribution of temperature anisotropies in the CMB. The CMB is very nearly
        isotropic and described, at the background level, by an equilibrium Bose-Einstein distribution
        
    f_0 = 1/^ p/T-1 .

        Deviations from this background distribution are usually parameterized as temperature
        perturbations
        
    f(τ,x̌,p,n̂) = [exp(p/T[1
                    +Θ(τ,x̌,n̂)])-1]^-1 .

        The temperature perturbation can then be written in terms of the distribution function as
        
    (1+Θ)^4-1 = ∫ p^3 p(f-f_0)/∫ p^3 p  f_0≡Δ .

        With our definition for the distribution function (<ref>), the deviations
        from (<ref>) are
        
    Δ = 1/γ^4(1-β̌·n̂)^4-1 + ℱ_γ .

        To first order in cosmological perturbations we have
        
    Θ   = 1/γ(1-β̌·n̂)-1 
                        + 1/4γ^3(1-β̌·n̂)^3ℱ_γ ,
       = β̌·n̂
                        + 1/4(1-3(β̌·n̂) + β^2)ℱ_γ
                        + β^2 .

        This would be the temperature perturbation observed in the 𝒪 frame. From the Sun's 
        reference system, the observed temperature perturbation Θ_⊙ is
        
    1+Θ_⊙ = γ_⊙(1-n̂·β̌_⊙)(1+Θ) ,

        where β̌_⊙ is the velocity of the Solar System in the 𝒪 frame.
        Expanding to leading order in β and β_⊙ we get
        
    Θ_⊙ = (β̌-β̌_⊙)·n̂ 
                    + 1/4(1+n̂·(β̌-β̌_⊙) 
                    - 4(n̂·β̌))ℱ_γ + β^2 .

        The reduced distribution function can be decomposed schematically as
        
    ℱ_γ (n̂,β̌) = ℱ_γ^ΛCDM(n̂)
                    +(n̂·β̌)ℱ_γ^β(n̂) + β^2 ,

        where ℱ_γ^ΛCDM follows the standard evolution and 
        ℱ^β_γ contains our modification, i.e. the imaginary part of the scalar modes
        and the vector modes. Finally, we need to take into account the aberration effects. The direction 
        n̂_⊙ observed from the Solar System is related to the direction n̂ in the
        𝒪 frame as
        
    n^i_⊙ = 𝒫^i_⊙ jn^j-γ_⊙β^i_⊙/γ_⊙
                    (1-n̂·β̌_⊙) .

        Then, to first order in β, we can express the direction as
        
    n^i = n^i_⊙ - (δ^i_j-n^i_⊙ n_⊙ j)(β^⊙  j_CMB-β^j) +β^2 ,

        where we have defined the relative velocity between the Sun and the CMB rest frame
        
    β̌^⊙_CMB≡β̌-β̌_⊙ .

        It is customary <cit.> to express the deflection instead as 
        
    β̌-(n̂·β̌) n̂ = ∇̌(n̂·β̌) .

        Taking everything into account, the temperature perturbation that would be measured from Earth
        in the direction n̂_⊙ is
        
    Θ_⊙(n̂_⊙) = n̂_⊙·β̌^⊙_CMB   + 1/4(1+n̂_⊙·β̌^⊙_CMB-4(n̂_⊙·β̌))
                    ℱ_γ^ΛCDM(n̂_⊙ 
                        -∇̌(n̂_⊙·β̌^⊙_CMB)
                        +∇̌(n̂_⊙·β̌))
        + 1/4(n̂_⊙·β̌)ℱ_γ^β(n̂_⊙) + β^2 .

        The first term represents the usual kinematic dipole, i.e. the Doppler-shifting effect associated 
        with the relative motion of the observer with respect to the CMB. The second term contains 
        a dipolar modulation and aberration effects. Both effects produce a kinematic mixing of the 
        multipole coefficients. The third term is a purely dynamical contribution, i.e. the effect
        of a relative motion between different species during the evolution. While recovering standard
        results <cit.> for β̌=0, in our setting we observe two kinds of new effects.
        First, the directions of the dipole, the dipolar modulation and the aberration effects do not
        coincide. This effect comes from the fact that, in our scenario, the standard ΛCDM evolution
        is recovered in the 𝒪 frame and not in the CMB rest frame. In standard 
        cosmology both frames coincide and this difference does not arise. The second effect is an additional
        source of statistical anisotropy, coming from the modified evolution, with a dipolar pattern.

        
        The CMB dipole is very well measured, with the latest Planck value being
        β^⊙_CMB = (1.23357± 0.00036)× 10^-3 <cit.>. It is widely 
        accepted that its origin is mostly kinematical, 
        so it gives us a very precise measurement of our relative motion with respect to the CMB. The
        Planck Collaboration also measured our relative motion using the kinematic
        correlations induced between different multipoles and the resultant anisotropic signal 
        <cit.>.
        Even though the uncertainties are large in this case, and there seems to be some tension 
        <cit.>, the velocity inferred using this method is compatible with the dipole, supporting its
        kinematical origin.  The relative velocity with respect to the CMB frame is usually interpreted
        as the result of peculiar motions of the Sun and the Local Group <cit.>. However, in our
        scenario, the relative velocity would arise as a combination of the local motion, with respect
        to the matter frame, and the relative motion between the matter and CMB frames. This gives rise
        to distinctive phenomenological consequences.
        
        
        At the background level, the non-coincidence of the CMB and matter frames produces a global
        motion of large-scale structures with respect to the CMB. This effect could be potentially 
        observed as a bulk flow on the largest scales. Measurements of bulk flows, at different scales,
        have been carried out in peculiar velocity surveys  and using the kinetic
        Sunyaev-Zeldovich (kSZ) effect <cit.>. See, e.g., Table 5 
        of <cit.> and references therein for a collection of recent measurements. 
        Although there is a long history 
        of conflicting measurements and anomalously large flows on cosmological scales, in this work
        we adopt the reported limit of Planck <cit.> for two reasons. In the first place, it extends to 
        the largest scales, up to 2 Gpc, where a cleaner determination of our global flow is expected.
        In the second place, it sets the more conservative bound in our parameter β_0, in the 
        sense of being the more restrictive to us. From the reported Planck value 
        v<254 km/s (95% CL), and according to the
        time evolution of β in Figure <ref>, we obtain
        
    β_0 < 1.6×10^-3 (95% CL) .

        Note that, since we are performing a first-order computation, all our results scale trivially 
        and we will write them explicitly in units of β_0. The previous constraint is wholly 
        compatible with the local measurements of peculiar motions mentioned above. The peculiar velocity 
        of the Local Group, and other higher order structures, with respect to the CMB is inferred from 
        the movement of the Sun with respect to both of them. The constraint (<ref>) 
        yields a value β < 0.85× 10^-3 today, of the same order as the measured velocity 
        β^⊙_LG=(1.00± 0.05)× 10^-3 <cit.>, so it can be 
        accomodated without fine-tuning the directions of these relative velocities. Potentially,
        it could even constitute a component of unaccounted peculiar motions of the largest 
        structures <cit.>. This constraint also justifies our first order computation.
        In section <ref> it was discussed how to construct a RW background to β
        and how the β^2 terms introduce anisotropies, i.e. a Bianchi background. Using 
        (<ref>) we can see that the terms β^2 are in fact smaller than 
        than a typical cosmological perturbation. Therefore, it is completely justified to take the 
        RW metric (<ref>) as the background geometry.

        
        At the perturbation level, it can be proven that, to first order in β, our modification
        does not leave an imprint in the CMB temperature spectrum, i.e. C_ℓ's. To lowest order, 
        the first CMB signatures appear as deviations from statistical isotropy. It is very important
        to notice that our model produces a distinctive signature in the CMB. In the standard picture,
        as mentioned before, the motion of Earth produces a violation of statistical isotropy in the
        CMB. In our case, there would be an additional, purely dynamical, source of statistical 
        anisotropy, caused by the relative motion between matter and radiation during the evolution. 
        Both effects could in principle be disentangled. We leave this analysis for future work 
        <cit.>, focusing instead on LSS observables.

        
        The local motion of the Earth also leaves an imprint in the observed galaxy distribution
        <cit.>, even though the analysis is not straightforward in this 
        case. Upcoming galaxy surveys like Euclid <cit.> or SKA <cit.> 
        will measure the induced dipole with high precision. A significant difference between this dipole 
        and the CMB result would be difficult to accomodate in standard ΛCDM, but could be easily
        interpreted as the result of a relative velocity between the CMB and matter frames. Even 
        if no such difference is measured, the bulk motion can still be smaller than the local one, 
        and yet lead to observational signatures, as we will immediately see.

                 
        Until now, we have mainly discussed the effects on the CMB temperature perturbations.
        Another class of observables comes from the clustering of matter. The distribution and redshift 
        of galaxies give us information about density perturbations and peculiar velocities. Again, in 
        the standard matter power spectrum we do not have a first order effect. On general grounds, if 
        we consider a cosmological quantity g splitted as in (<ref>), we have
        
    |g(τ,ǩ)|^2 = |g^R(τ,ǩ)|^2+β^2 ,

        i.e. the standard result. However, in the cross-correlation between two cosmological perturbations
        g_1 and g_2 we get a first-order dipolar contribution
        
    g_1g_2^* = g_1^Rg_2^R * +  (β̂·k̂)
                    (g_1^Ig_2^R  *-g_1^Rg_2^I *) + β^2 .

        Every cross-correlation between cosmological quantities contains a dipole modification with
        this structure. This effect could be observed in the future in the cross-correlations
        between matter density and velocity <cit.>, as the precision of the surveys increases. It is
        conceivable that this effect could appear in cross-correlations between baryon and CDM 
        densities as well, even though a thorough analysis using lensing information would be in
        order. Finally, 
        the generation of vorticity, purely decaying in ΛCDM, is another distinctive feature of
        our model.

        
        Since the most accessible observables are related to velocity perturbations, we conclude this
        section clarifying a few points concerning our previous definitions. In particular, it is 
        important to relate the intermediate variables we have used with the physical velocities. The 
        velocity that would appear in the energy-momentum tensor of a fluid (<ref>) 
        is the velocity of a frame in which the momentum flux, i.e. the component T_0i, is zero 
        <cit.>. Using the boost-transformation properties (<ref>) 
        we can obtain an equation for the physical velocity U^i of the fluid
        
    γ̅𝒫̅^i_j(Q^j-U_kΠ^kj)
                        -γ̅^2U^i(ρ+P-Q^jU_j) = 0 ,

        where
        
    γ̅≡ (1-U^2)^-1/2 ,    𝒫̅^i_j≡δ^i_j + (γ̅-1)U^iU_j/U^2 .

        Working to first order in β and in cosmological perturbations we have
        
    U^i = β^i + 1/ρ̃+P̃(δ Q^i - β^kδΠ^i_k
                    -β^i(δρ + δ P)) .

        The physical velocity has two parts, a bulk velocity β_i plus a peculiar contribution
        δ u^i. For ultrarelativistic particles, the peculiar velocity can be expressed in terms
        of our previously defined variables (<ref>) as
        
    δ u^i = δ v^i -3/4β^kπ_k^i-β^iδ .

        It can be splitted into a scalar and a vector part
        
    δ u^i = -k̂^i/kϑ + ζ^i ,

        so that we have
        
        
    ϑ   = θ - (β̌·ǩ)(δ - σ) ,
    ζ^±   = χ^± - (β̌·ê_±)(δ+1/2σ) .

        
        For non-relativistic species, the results are identical setting σ=0. It is worth noting
        that this is not the only physically sensible definition of the velocity of a fluid. It can
        alternatively be defined as the velocity of the frame in which the flux of particles
        (<ref>) is zero <cit.>. Both definitions agree for non-relativistic
        fluids if the number of particles is conserved.
        
    

 §.§ Time evolution and transfer functions

        The time evolution of the bulk velocities for the different components is represented in Figure
        <ref>. All the components in the visible sector start their evolution with the 
        same velocity in the center of mass frame, and its momentum is counterbalanced by the dark sector. 
        The velocity of the neutrinos, since we are neglecting their masses, is always constant. 
        The velocity of the photon-baryon plasma is constant deep in the radiation-dominated era. Once 
        the baryonic contribution to the energy density becomes important, the plasma velocity drops 
        down as a^-1, see (<ref>),
        until decoupling. After decoupling, the velocity of the baryons keeps scaling as a^-1, like
        CDM, while the photons keep a constant velocity, with a slight late-time effect from reionization.
        Today, the cosmic center of mass, i.e. the 𝒪 frame, almost coincides with the matter 
        frame but photons and neutrinos possess a sizeable velocity.

        
        
        
        Finally, we present the transfer functions evaluated today for a range of k and their time 
        evolution for a fixed value k=10^-2 Mpc^-1, as a sample from the full results for 
        the evolution of the perturbations. All the results 
        concerning cosmological perturbations are computed in the synchronous gauge and then transformed
        back to the Newtonian gauge, that can be more easily interpreted in the Newtonian limit
        <cit.>. Figures <ref>, <ref>,
        and <ref> contain the density and velocity of CDM, baryons and the full fluid. 
        Figure <ref> contains the metric variables, i.e. Newtonian potentials and
        vector metric perturbations. The quantities with an R superscript follow the standard evolution 
        and are computed using class with the Planck 2018 <cit.> input 
        values. The modified contributions, with an I superscript, remain smaller than the standard 
        ones for most values of k, but not as small as could be expected. The difference at scales of 
        0.1 Mpc^-1 is just one order of magnitude, instead of three as could be naively 
        anticipated from β_0=10^-3, and the modifications could grow even larger above the 
        non-linearity scale.        
        
        
        
        
        
        	    
        
                        
	    
    

 §.§ Spectra
         
        
                
        We stick to the following conventions for the definition of the spectra. In the
        first place, we define the spectrum of a single variable as
        
    δ(z, ǩ)δ^*(z,ǩ')   ≡δ (ǩ-ǩ')P_δδ(z,k)
       = δ(ǩ-ǩ')|δ^R(z,k)|^22π^2/k^3𝒫_ℛ(k)
                        +β^2 ,

        where δ^R(z,k) is the real part of the transfer function and 𝒫_ℛ(k) is the usual
        nearly scale-invariant curvature spectrum. For the cross-correlations we define
        
    Reδ(z,ǩ)θ^*(z,ǩ')   ≡δ(ǩ-ǩ')P^R_δθ(z, k)
       = δ(ǩ-ǩ')δ^R(z,k)θ^R(z,k)2π^2/k^3𝒫_ℛ(k) + β^2,
    Imδ(z,ǩ)θ^*(z,ǩ')   ≡δ(ǩ-ǩ')(β̂·k̂)P_δθ^I(z, k)
       = δ(ǩ-ǩ')(β̂·k̂)(δ^I(z,k)θ^R(z,k)-δ^R(z,k)ϑ^I(z,k))2π^2/k^3𝒫_ℛ(k)+β^2 .

        The vorticity spectrum is defined in a similar way, according to definition 
        (<ref>),
        
    ζ̌(z, ǩ)·ζ̌^*(z,ǩ')   ≡δ (ǩ-ǩ')(1-(β̂·k̂)^2)P_ζζ(z,k)
       = δ(ǩ-ǩ')(1-(β̂·k̂)^2)|ζ(z,k)|^22π^2/k^3𝒫_ℛ(k)+β^3 ,

        The results for the velocity spectrum and the density cross-correlation for CDM are 
        represented in Figure <ref>. In addition to this information, the velocity-density
        cross-correlation induced by our modification shows a distinctive dipolar pattern. In the
        same way, the vorticity autocorrelation, even though its amplitude is very small, deviates from
        statistical isotropy, with a quadrupole term in addition to the monopole. Figure 
        <ref> shows the cross-correlation between the matter density and 
        the lensing potential, defined as Ψ≡ϕ+ψ. This combination is observable 
        using weak-lensing information <cit.>. Again, our additional 
        contribution becomes important at small scales, being just one order of magnitude below
        the standard result at scales k=0.1 Mpc^-1 instead of three as might be inferred
        from β_0 = 10^-3.
        Both vorticity and deviations from statistical isotropy are absent in standard ΛCDM. 
        Their presence, with the structure proposed, is a testable effect that could be used to confirm, 
        or disprove, the non-comoving scenario.

        
        As we have seen, in our modified setting, velocity spectra are the most easily accessible 
        LSS observables that show significant deviations. Peculiar velocity surveys provide useful 
        complementary information but currently are not competitive with other cosmological observables 
        to constrain standard cosmology <cit.>. Nonetheless, as the 
        precision of velocity surveys increase, such an effect might be seen or at least prove more 
        competitive than the dipole measurements in the CMB and galaxy distribution to 
        constrain β_0 <cit.>.
                   
        


§ CONCLUSIONS

    In this work we have developed the theoretical framework needed to analyze the cosmology of
    non-comoving fluids. We have shown that it is possible to relax one of the underlying assumptions
    of ΛCDM, comoving CMB and matter frames, while retaining an homogeneous and isotropic
    universe. Even if the background behaviour is preserved, the evolution of the perturbations is
    modified, leading to new phenomenological signatures. To first order in the relative velocities,
    i.e. β, we reproduce ΛCDM behaviour in the main cosmological observables, like
    the matter power spectrum and the CMB temperature power spectrum. We have postponed the full CMB
    analysis <cit.>, focusing instead on LSS observables in this work. As mentioned before, 
    we find that the effects on the autocorrelation spectra are β^2 but there are 
    β effects on every cross-correlation between cosmological perturbations. Additionally, 
    these corrections present a dipolar pattern, producing deviations from statistical 
    isotropy. We have observed that the additional contributions to the cross-correlation spectra
    become important for small scales.

    
    Another distinctive feature of this model is the production of vorticity, which is absent in 
    ΛCDM. The relative motion of the fluids induces β couplings between scalar and 
    vector modes. This in its turn leads to the production of vorticity and vector metric perturbations, 
    sourced by the scalar modes. This vector contribution also leaves a characteristic imprint in the 
    velocity spectrum, with a statistically anisotropic quadrupolar modulation. No tensor modes are 
    excited to β.

    
    Our only additional free parameter, β_0, is the initial velocity between the visible and the 
    dark sector, in the frame that observes a homogeneous and isotropic background. Measurements of
    large-scale bulk flows and the CMB dipole allow us to set a conservative limit
    
    β_0 < 1.6× 10^-3 (95% CL) .

    Satisfying this constraint, our model is compatible with current observations, and yet it can
    have striking phenomenological consequences. There is a work under way to study the impact on
    the CMB and to give concrete predictions about violations of statistical isotropy
    <cit.>. In particular, it may alleviate the tension that seems to arise when interpreting 
    the anisotropic signal in the CMB as a pure kinematical effect <cit.>. This tension corresponds
    to the dipolar modulation anomaly <cit.> that arises when analyzing 
    low multipoles. The modulation points in a direction different from the kinematic dipole and seems to 
    have a larger amplitude, even if it is compatible with zero within 2σ <cit.>.

    
    More importantly, the production of vorticity
    for the photon-baryon plasma opens an avenue for the creation of magnetic fields <cit.>.
    The origin of the galactic magnetic fields is a long-standing open problem <cit.>.
    The Harrison mechanism <cit.> is a cosmological production mechanism that needs vorticity in the 
    photon-baryon plasma to operate, but unfortunately it is absent in ΛCDM to first order
    in cosmological perturbation theory. Several studies of second-order cosmological perturbation
    theory have proven that in this case vorticity, and thus magnetic fields, is created but with an
    amplitude far too small to act as seed fields for the galactic dynamo amplification mechanism 
    <cit.>. Our setup is similar in some regards to a second 
    order computation, but our relevant scale β_0 is larger than a typical cosmological 
    perturbation and so it is our vorticity production. The associated magnetic fields in our case are 
    expected to be larger and their spectrum would have a different tensor structure as well, since we are 
    singling out a privileged direction <cit.>.

    
    After setting up the formalism for a non-comoving
    cosmology and proving its viability, this work paves the way for the search of these exciting
    new signatures.
    

    This work has been supported by the MINECO (Spain) projects FIS2014-52837-P, 
    FIS2016-78859-P(AEI/FEDER, UE), and Consolider-Ingenio MULTIDARK CSD2009-00064.
    Some symbolic computations have been performed using SymPy <cit.>
    and cadabra <cit.>.




§ GEODESICS

    In this appendix we will compute in detail the geodesics for a general perturbed RW metric
    
    s^2 = a^2(τ)(-(1-A)τ^2 + 2B_iτ x^i + (δ_ij+H_ij) x^i x^j) ,

    that is
    
    g_μν = a^2[      -1+A       B_i;       B_i δ_ij+H_ij ],    
            g^μν = 1/a^2[      -1-A       B^i;       B^i δ^ij-H^ij ] ,

    and with the constraint
    
    g^μνP_μ P_ν = -m^2 .

    Defining the proper-time parameter as λ≡√(- s^2), the standard definitions
    for the 4-velocity and 4-momentum are
    
    U^μ   ≡x̣^̣μ̣λ ,
    
            P_μ   ≡ m U_μ ,
    
            P_0    ≡ -ϵ + δ P_0 .

    The geodesics are given by
    
    Ụ^̣μ̣λ +Γ^μ_ν_ρU^ν U^ρ = 0 ,

    that can be conveniently rewritten as
    
    Ụ_̣μ̣λ = 1/2g_νρx^μU^ν U^ρ .

    This last form is especially useful. Since the background metric is homogeneous, we can keep only
    the zero order in U_μ to compute the spatial part.
    Writing the evolution in terms of the conformal time we have
    
    P̣_̣ịτ = m/U^0Ụ_̣ịλ = 1/2ϵ(ϵ^2∂_iA
                +2ϵ P^j∂_iB_j + P^jP^k∂_iH_jk) .

    The spatial momentum is redefined as
    
    P_i    ≡(δ^j_i+1/2H^j_i)q_j ,
    
            a^2P^i    = -ϵ B^i + (δ^ki-1/2H^ki)q_k ,

    where, from now on, every spatial index on a perturbed quantity is assumed to be raised or lowered 
    using δ_ij. From the mass-shell condition (<ref>) we obtain
    
    ϵ^2    = m^2a^2 + q^2 ,
    
            P_0    = -ϵ(1-1/2A) + q_iB^i  ,
    
            P^0    = ϵ/a^2(1+1/2A) .

    The geodesic equation can be written with this parameterization as
    
    q̣_̣ịτ = 1/2ϵ∂_iA + q^j∂_iB_j + 1/2ϵq^jq^k
                (∂_iH_jk-∂_kH_ij) - 1/2q^jḢ_ij .

    Further decomposing q_i into direction n̂ and magnitude q
    
    q_i≡ q n_i    →     n_iδ^ijn_j= 1 ,   q^2 = δ^ijq_iq_j ,

    we have
    
    
    q̣τ   = n^iq̣_̣ịτ ,
    ṇ_̣ịτ   = 1/q(δ^j_i-n^jn_i)q̣_̣j̣τ .

    
    Finally, using the following succint redefinition of metric variables
    
    C_ij   ≡∂_iB_j - 1/2Ḣ_ij ,
    
            D_ijk   ≡1/2(∂_iH_jk-∂_kH_ij) ,

    the final formulae needed to compute the geodesics, with the parameterization 
    (<ref>), are
    
    
    x̣^̣ịτ   = q^i/ϵ(1-1/2A)-B^i
                -1/2ϵH^i_kq^k ,
    q̣^̣ịτ   = 1/2ϵ∂_iA + q^jC_ij + q^jq^k/ϵD_ijk ,
    q̣τ   = 1/2ϵ  n^i∂_iA + q n^in^jC_ij ,
    ṇ_̣ịτ   = (δ^j_i-n^jn_i)(
                ϵ/2q∂_jA + n^kC_jk + q/ϵn^kn^lD_jkl) .

    
    
    Some relevant metric quantities, written in terms of the scalar-vector-tensor decomposition and in
    Fourier space, are
    
    A    = -2ψ ,
    
            B_i    =  k_i B -S_i ,
    
    	    H_ij   = -2ϕδ_ij - 2k_ik_jE+ (k_iF_j+k_jF_i) + h_ij ,
    
    	    C_ij   = -k_ik_j(B-Ė)- k_iS_j + ϕ̇δ_ij - /2(k_iḞ_j
    	        +k_jḞ_i) - 1/2ḣ_ij ,
    
    	    D_ijk   = -(k_iδ_jk-k_kδ_ij)ϕ - 1/2(k_ik_jF_k-k_jk_kF_i)
    	        +/2(k_ih_jk-k_kh_ij) .

    


§ FULL COMPUTATION OF THE COLLISION TERM

    Starting from the collision term (<ref>),
    
    𝒞[f]=σ_T/4π p∫p̃'p̃'Ω̃'
                [ñ_e^fullδ (p̃-p̃')+ñ_eũ̌^full_e·
                (p̌̃̌-p̌̃̌')δ (p̃-p̃')p̃']
                (f̅(Λ_βΛ_β_e^-1p̌̃̌')
                -f̅(Λ_βΛ_β_e^-1p̌̃̌)) ,

    we need to solve two integrals
    
    I_1    = ∫p̃'p̃'Ω̃'δ (p̃-p̃')
                (f̅(Λ_βΛ_β_e^-1p̌̃̌')
                -f̅(Λ_βΛ_β_e^-1p̌̃̌))
       = -4πp̃f̅(Λ_βΛ_β_e^-1p̌̃̌)
                    +p̃∫Ω̃' f̅(Λ_βΛ_β_e^-1p̌̃̌')|_p̃'=p̃ ,
    Ǐ_2    = ∫p̃'p̃'Ω̃' (p̌̃̌-p̌̃̌')
                δ (p̃-p̃')p̃'(f̅(Λ_βΛ_β_e^-1p̌̃̌')-f̅(Λ_βΛ_β_e^-1p̌̃̌))
       = -p̃∫Ω̃'(ñ̂-2ñ̂')
                    (f̅(Λ_βΛ_β_e^-1p̌̃̌')|_p̃'=p̃
                    -f̅(Λ_βΛ_β_e^-1p̌̃̌))
                    +p̃^2∫Ω̃'(ñ̂-ñ̂')
                    f̅(Λ_βΛ_β_e^-1p̌̃̌)p̃'|_p̃'=p̃ .

    Now, we will integrate out the dependence on the momentum p in the 𝒪 frame, like we
    did with the left-hand side of the Boltzmann equation. First, for simplicity, we will assume that
    β̌ and β̌_e point in the same direction, so we can obtain
    
    Λ_βΛ_β_e^-1 = Λ_Δβ ,    Δβ≡β-β_e/1-ββ_e .

    This assumption simplifies the derivation for arbitrary values of β, but it is not needed and
    in fact the first order results are independent of it. Since the expressions are already
    quite cumbersome, and this will be the only physical configuration of interest, we will adopt this 
    assumption throughout this appendix. Some preliminary results and definitions are
    
    
    ∫ q^3 q f̅_0(Λ_β q)    = 𝒩̃/γ^4(1-n̂·β̌)^4 ,
    ∫^3q̃/4π q̃ f̅_0(Λ_Δβq̃) 
                   = 𝒩̃γ^2_Δβ(1+Δβ^2/3) ,
    ∫^3q̃/4π q̃^i f̅_0(Λ_Δβq̃) 
                   = 4/3𝒩̃γ^2_ΔβΔβ^i ,

    
    where 𝒩̃ is defined in (<ref>) and 𝒫^i_e j, 
    γ_e correspond to (<ref>), evaluated with β_e. We can proceed now to 
    compute the integrals term by term. First for I_1
    
    
    ∫ q  q^2q̃ f̅(Λ_βΛ_β_e^-1q̌̃̌)
                   = 𝒩̃γ_e(1-n̂·β̌_e)(ℱ_γ 
                    + 1/γ^4(1-n̂·β̌)^4) ,
    ∫ q q^2q̃∫Ω̃ f̅(Λ_βΛ_β_e^-1q̌̃̌)
                   =4π𝒩̃/γ^3_e(1-n̂·β̌_e)^3[(1
                    +Δβ^2/3)γ_Δβ^2 
                    +γ_e^2∫Ω/4π(1-2n̂·β̌_e + (n̂·β̌_e)^2)ℱ_γ] .

    
    And for Ǐ_2
    
    ∫ q^2 q q̃^i f̅_0(Λ_Δβq̃)
                   = 𝒩̃𝒫^i_e jn^j-γ_eβ^i_e/γ^4(1-n̂·β̌)^4 ,
    ∫ q^2 q q̃^i ∫Ω̃ f̅_0(Λ_Δβq̃)
                   = 4π𝒩̃𝒫^i_e jn^j-γ_eβ^i_e/γ^4_e(1-n̂·β̌_e)^4γ^2_Δβ(1+Δβ^2/3) ,
    ∫ q^2 q ∫Ω̃ q̃^i f̅_0(Λ_Δβq̃)
                   = 16π/3𝒩̃Δβ^iγ_Δβ^2/γ_e^3(1-n̂·β̌_e)^3 ,
    ∫ q^2 q q̃ q̃^i∫Ω̃ f̅_0(Λ_Δβq̃)q̃   = -16π𝒩̃𝒫^i_e jn^j-γ_eβ_e^i/γ_e^4(1-n̂·β̌_e)^4γ_Δβ^2(1+Δβ^2/3) ,
    ∫ q^2 q q̃∫Ω̃ q̃^i f̅_0(Λ_Δβq̃)q̃   = -64π/3𝒩̃Δβ^iγ_Δβ^2/γ_e^3(1-n̂·β̌_e)^3 .

    Finally, for the electron quantities,
    
    
    ñ_e^full   ≡ 2∫^3p̃_e/(2π)^3f̃_e(p̌̃̌_e) = ñ_e + δñ_e ,
    ñ_eũ̌^full_e    = 2∫^3p̃_e/(2π)^3p̌̃̌_e/Ẽ_p_ef̃_e(p̌̃̌_e)
                = 2∫^3p̃_e/(2π)^3p̌̃̌_e/Ẽ_p_eδf̃_e(p̌̃̌_e)
                = ñ_eδṽ̌_e ,

    
    where ñ_e is the physical background electron number density, computed in its comoving frame as every
    other background quantity. Using the Lorentz transformation properties
    
    
    δñ_e    = γ_e(δ n_e - δ v_e^jβ_e  j) ,
    δṽ_e^i    = 𝒫^i_e jδ v^j_e - γ_eβ_e^iδ n_e .

    
    The final results, to first order in β, are
    
    
    ∫ q^2  q I_1    = -4π𝒩̃((1-n̂·β̌_e)ℱ_γ +4n̂·Δβ̌-(1+3n̂·β̌_e)∫Ω/4πℱ_γ
                    +2β̌_e·∫Ω/4πn̂ ℱ_γ) ,
    ∫ q^2  q Ǐ_2    = 4π𝒩̃(-8/3Δβ̌+4n̂ (n̂·Δβ̌)
                    +4(n̂-β̌_e+4n̂ (n̂·β̌_e))) ,
    δñ_e    = δ n_e-β_e jδ v^j_e ,
    δṽ^i_e    = δ v^i_e-β^i_eδ n_e .

    
    


§ GAUGE TRANSFORMATIONS WITH NON-COMOVING FLUIDS

    After an infinitesimal gauge transformation
    
    Δ x^μ = ϵ^μ ,     ϵ^μ = (T(τ,x̌), Ľ(τ,x̌)) ,

    a tensor changes as
    
    Δ T_μν = ϵ^ρT_μνx^ρ+ϵ^ρ_,μT_ρν+ϵ^ρ_,νT_μρ=ℒ_ϵ T_μν ,

    where ℒ_ϵ is the Lie derivative. Applying the result to the metric tensor, we have
    for the metric variables
    
    
    Δ A    = -2(Ṫ+ℋT) ,
    Δ B_i    = L̇_i-∂_iT ,
    Δ H_ij   = 2ℋTδ_ij+∂_iL_j+∂_jL_i .

    
    On the other hand, for the perturbed fluid variables we get
    
    
    Δδρ   = Tρ̇-2(ρ+P)β^i∂_iT ,
    Δδ Q_i    = -(ρ+P)∂_iT+T∂_τ(β_i(ρ+P))
                +1/2(ρ+P)β^j(∂_iL_j-∂_jL_i) ,
    Δδ P    = TṖ-2/3β^i∂_iT(ρ+P) ,
    ΔδΠ_ij   = -(ρ+P)(β_i∂_jT+β_j∂_iT-2/3δ_ijβ^k∂_kT) ,

    
    where we are adopting the definitions (<ref>). With our previous definition
    for the scalar-vector-tensor decomposition of the metric variables (<ref>), 
    we get, in Fourier space,
    
        Δψ   = Ṫ+ℋT ,            ΔS_+    = -L̇_+ ,

        ΔB    = -/kk̂·Ľ̇-T ,            ΔF_+    = L_+ ,

        Δϕ   = -ℋT ,            Δh_++    = 0 ,

        ΔE    = -/kk̂·Ľ .
    
    The results for the - helicity can be obtained substituting -↔ +
    in every sub and superscript.
    If we want to change from the synchronous to the Newtonian gauge, for scalar perturbations, the 
    following conditions must be satisfied
    
    
    ψ   = Ṫ+ℋT ,
    
            0    = -/kk̂·Ľ̇-T ,
    ϕ-η   = -ℋT ,
    1/2k^2(h+6η)    = -/kk̂·Ľ ,

    
    that can be solved to yield
    
    T    = 1/2k^2(ḣ+6η̇) ,
    ψ   = Ṫ+ℋT ,
    ϕ   = η+ℋT .

    Finally, for the fluid variables that we use in the main part of the computations, we have the rules
    
    δ (Newt)-δ (Syn)    = Tρ̇/ρ-2(1+w)(β̌·ǩ)T ,
    θ (Newt)-θ (Syn)    = k^2T+(β̇̌̇·ǩ)T +(β̌·ǩ)ρ̇+Ṗ/ρ+PT ,
    χ_+(Newt)-χ_+ (Syn)    = Tβ̇_++Tβ_+ρ̇+Ṗ/ρ+P ,
    δ P (Newt)-δ P(Syn)    = TṖ-2/3(β̌·ǩ)T(ρ+P) ,
    σ (Newt)-σ (Syn)    = (β̌·ǩ)(1+w)T ,
    π^V_+ (Newt)-π^V_+ (Syn)    = -β_+(1+w)T ,
    π^T_++ (Newt)-π^T_++ (Syn)    = 0 ,
    
    where we have neglected terms β L_+, making use of the fact that, according to 
    (<ref>), the transverse part of Ľ under our assumptions can be at most 
    order β. Again, we are omitting the results for the - helicity, that can be obtained 
    substituting -↔ + in every sub and superscript.
    


§ INITIAL CONDITIONS

    In this appendix we will find the appropiate initial conditions for the system of scalar and
    vector modes in sections <ref> and <ref>. We will consider the
    most general initial condition and then study the physical restrictions that we must impose.
    For β=0, our system reproduce the standard cosmology. This case has been extensively studied
    over the years and the relevant modes, i.e. one adiabatic and four isocurvature modes, have been
    identified <cit.>. In our setup, the presence of an external source gives rise to
    the existence of a new “mode” of the system, in the sense that we have a non-trivial evolution
    even if the usual adiabatic and isocurvature modes are absent. First, we will identify this particular
    solution, 
    setting to zero the other modes of the system. Note that the external sources only contain variables
    that evolve according to standard ΛCDM, so for these variables only adiabatic initial 
    conditions are considered. After identifying the effect of the sources, the most general perturbation
    can be constructed adding to the sourced mode the adiabatic and isocurvature modes. Finally, we must
    analyze what physical requirements constrain our choice of initial conditions. In particular, we
    impose that neutrinos and photons, being tightly coupled in the very early Universe, share a common
    initial velocity. Every other initial condition that is not fixed by this condition is set to zero.
    This programme is carried out in detail in the next sections.
    
    

 §.§ Scalar modes
    
	    To obtain our results, we have analyzed the most general type of perturbation, reproducing the
	    results of <cit.> but with a slight change of notation. In the first place, we use
	    an alternative approach where we integrate (<ref>) instead of the dark matter
	    equation. In this setup, matter isocurvature modes appear when taking non-zero initial conditions 
	    for δ_ν, δ_b or ḣ. In the second place, even though it is perfectly 
	    equivalent, we parameterize the neutrino isocurvature velocity mode with the initial value of 
	    θ_γ instead of θ_ν.

	    
	    Before presenting the results, some shorthand definitions that will be used later are
	    
    S_γν   ≡Ω_γ + Ω_ν ,
    ℛ_s    ≡Ω_s/S_γν ,     s=γ,ν,b,c ,
    ℛ_bc   ≡ℛ_b+ℛ_c ,

	    Assuming a universe composed of radiation and matter, where τ stands for conformal time, 
	    
    ℋ   = H_0ℛ_b c√(S_γν)τ/2 + 1/τ(H_0ℛ_b c√(S_γν)τ/4 + 1) ,
    
    	        a    = H_0√(S_γν)τ(H_0ℛ_b c√(S_γν)τ/4 + 1) .

	    During the radiation-dominated phase, they can be expanded as
	    
    a    = H_0√(S_γν)τ + O(H_0^2ℛ^2_bcS_γντ^2) ,
    ℋ   = 1/τ + H_0ℛ_b c√(S_γν)/4 - H_0^2ℛ_b c^2 S_γντ/16 + O(H_0^2ℛ^2_bcS_γντ^2) ,
    ℋ^2   = 1/τ^2 + H_0ℛ_b c√(S_γν)/2 τ - H_0^2ℛ_b c^2 S_γν/16 + O(H_0^2ℛ^2_bcS_γντ^2) ,
    ℋ̇   = - 1/τ^2 - H_0^2ℛ_b c^2 S_γν/16 + H_0^3ℛ_b c^3 S_γν^3/2τ/32 + O(H_0^2ℛ^2_bcS_γντ^2) .

	    Now, if we look for regular super-Hubble solutions and expand every cosmological variable as 
	    
    δ^I_γ = D^(0)_δ_γ + D^(1)_δ_γτ + D^(2)_δ_γτ^2 + D^(3)_δ_γτ^3 + …

	    The results for the sourced mode only (setting the initial conditions for η^I, ḣ^I, 
	    δ_ν^I, δ_b^I and θ_γ^I to zero), can be written in terms of the initial
	    value of ψ^R(τ=0)=Ψ or
	    
    η^R(τ=0) = 4ℛ_ν+15/10Ψ ,

	    as
	    
	        D^(0)_δ_γ    = 0 ,             D^(1)_δ_γ    = 4 Ψβ_0 k (4 ℛ_ν+ 15)/15 ,

	        D^(0)_δ_ν    = 0 ,                D^(1)_δ_ν    = 4 Ψβ_0 k (ℛ_ν- 1) (4 ℛ_ν+ 15)/15 ℛ_ν ,

	        D^(0)_δ_b    = 0 ,                  D^(1)_δ_b    = Ψβ_0 k (4 ℛ_ν+ 15)/5 ,

	        D^(0)_δ    = 0 ,                  D^(1)_δ    = 0 .
	    
	
	    
    D^(2)_δ_γ   = - H_0Ψ√(S_γν)ℛ_b β_0 k (ℛ_ν - 3) (4 ℛ_ν + 15)/20 (ℛ_ν - 1) ,
    
    	        D^(2)_δ_ν   = - H_0Ψ√(S_γν)ℛ_b β_0 k (4 ℛ_ν + 15)/20 ,
    
    	        D^(2)_δ_b   = - 3 H_0Ψ√(S_γν)ℛ_b β_0 k (ℛ_ν - 3) (4 ℛ_ν + 15)/80 (ℛ_ν - 1) ,
    
    	        D^(2)_δ   = H_0Ψ√(S_γν)ℛ_b β_0 k (4 ℛ_ν + 15)/20 .

	
	    
    D^(3)_δ_γ   = H_0^2Ψ S_γνℛ_b β_0 k (4 ℛ_ν + 15) (5 ℛ_νℛ_bc(ℛ_ν - 1) - 3 ℛ_b (ℛ_ν - 6))/300 (ℛ_ν - 1)^2 + 4 Ψβ_0 k^3(44 ℛ_ν + 65)/225 (4 ℛ_ν + 5) ,
    
    	        D^(3)_δ_ν   = - H_0^2Ψ S_γνℛ_b β_0 k (4 ℛ_ν + 15) (3 ℛ_b - 5 ℛ_bc(ℛ_ν - 1))/300 (ℛ_ν - 1) + 2 Ψβ_0 k^3(32 ℛ_ν + 45)/225 ℛ_ν ,
    
    	        D^(3)_δ   = 3 H_0^2Ψ S_γνℛ_b β_0 k (4 ℛ_ν + 15) (2 ℛ_b - 5 ℛ_bc(ℛ_ν - 1))/400 (ℛ_ν - 1) + Ψβ_0 k^3(80 ℛ_ν^2 + 392 ℛ_ν + 545)/450 (4 ℛ_ν + 5) .

	
	    
	        D^(0)_θ_γ    = 0 ,             D^(1)_θ_γ    = 0 ,

	        D^(0)_θ_ν    = Ψβ_0 k (4 ℛ_ν+ 15)/5 ℛ_ν  ,           D^(1)_θ_ν    = 0 ,

	        D^(0)_θ    = 0 ,              D^(1)_θ    = 0 .
	    
	
	    
    D^(2)_θ_γ   = - Ψβ_0 k^3/6 ,
    
    	        D^(2)_θ_ν   = - Ψβ_0 k^3(44 ℛ_ν^2 + 151 ℛ_ν + 135)/30 ℛ_ν(4 ℛ_ν + 5) ,
    
    	        D^(2)_θ   = - Ψβ_0 k^3(ℛ_ν + 2) (4 ℛ_ν + 15)/15 (4 ℛ_ν + 5) .

	
	    
    D^(0)_σ_ν   = 0 ,
    
    	        D^(1)_σ_ν   = 2 Ψβ_0 k (ℛ_ν + 2) (4 ℛ_ν + 15)/15 ℛ_ν(4 ℛ_ν + 5) ,
    
    	        D^(2)_σ_ν   = 2 H_0Ψ√(S_γν)ℛ_bcβ_0 k (ℛ_ν + 2)/5 (4 ℛ_ν + 5) ,
    
    	        D^(3)_σ_ν   = H_0^2Ψ S_γνℛ_bc^2β_0 k (ℛ_ν + 2) (4 ℛ_ν - 45)/30 (2 ℛ_ν + 15) (4 ℛ_ν + 5) - Ψβ_0 k^3(32 ℛ_ν^4 + 224 ℛ_ν^3 + 914 ℛ_ν^2 + 2097 ℛ_ν + 1620)/270 ℛ_ν(2 ℛ_ν + 15) (4 ℛ_ν + 5) .

	
	    
	        D^(0)_η    = 0 ,                     D^(1)_η    = - 2 Ψβ_0 k (ℛ_ν+ 2) (4 ℛ_ν+ 15)/15 (4 ℛ_ν+ 5) ,

	        D^(0)_h    = D^(0)_h ,              D^(1)_h    = 0 .
	    
	
	    
    D^(2)_η   = H_0Ψ√(S_γν)β_0 k (- ℛ_b (4 ℛ_ν + 5) (4 ℛ_ν + 15) + 40 ℛ_bc(ℛ_ν + 2))/80 (4 ℛ_ν + 5) ,
    
    	        D^(2)_h   = 3 H_0Ψ√(S_γν)ℛ_b β_0 k (4 ℛ_ν + 15)/40 .

	    
	    
    D^(3)_η   = - H_0^2Ψ S_γνℛ_b^2β_0 k (4 ℛ_ν + 15)/400 (ℛ_ν - 1) + H_0^2Ψ S_γνℛ_b ℛ_bcβ_0 k (4 ℛ_ν + 15)/240
            + H_0^2Ψ S_γνℛ_bc^2β_0 k (ℛ_ν + 2) (4 ℛ_ν - 45)/24 (2 ℛ_ν + 15) (4 ℛ_ν + 5) + Ψβ_0 k^3(- 80 ℛ_ν^3 + 568 ℛ_ν^2 + 4525 ℛ_ν + 4950)/1350 (2 ℛ_ν + 15) (4 ℛ_ν + 5) ,
    
    		    D^(3)_h   = H_0^2Ψ S_γνℛ_b β_0 k (4 ℛ_ν + 15) (3 ℛ_b - 5 ℛ_bc(ℛ_ν - 1))/200 (ℛ_ν - 1) - Ψβ_0 k^3(80 ℛ_ν^2 + 528 ℛ_ν + 655)/450 (4 ℛ_ν + 5) .

	    
	    Once we have the new behaviour of the system, we need to evaluate the
	    assignment of initial conditions. It seems reasonable to give zero initial values to our modification
	    but there is one further physical requirement that we must impose. As mentioned in the main text, 
	    if neutrinos and photons were in thermal contact in the primeval Universe it is physically sensible
	    to impose that they shared the same velocity
	    
    θ_ν(τ=0) = θ_γ(τ=0) = θ_γ^(0) ,

	    In the standard scenario this leads to θ_γ^(0)=0 and to the absence of neutrino
	    velocity isocurvature modes. However, in our case, if we consider a neutrino isocurvature velocity
	    mode on top of the sourced mode, upon imposing this restriction we get
	    
    θ_γ^(0) = 4ℛ_ν+15/5Ψβ_0 k .

	    In order to obtain the correct initial conditions, we must consider the combination of the sourced
	    mode with a neutrino isocurvature velocity mode with the previous initial condition. The final 
	    results are
	    
    δ_γ^I    = δ_ν^I = δ_b^I = δ^I = 0 + τ^3 ,
    θ_γ^I    = Ψβ_0 k (4 ℛ_ν + 15)/5 
    	                + 3 H_0Ψ√(S_γν)ℛ_b β_0 k (4 ℛ_ν + 15)/20 (ℛ_ν - 1)τ
         + 3 H_0^2Ψ S_γνℛ_b β_0 k (4 ℛ_ν + 15) ((ℛ_ν-1) ℛ_bc + 3 ℛ_b)/80 (ℛ_ν - 1)^2τ^2
    	                - 2 Ψβ_0 k^3(ℛ_ν + 5)/15τ^2 + τ^3 ,
    θ_ν^I    = Ψβ_0 k (4 ℛ_ν + 15)/5
    	                - Ψβ_0 k^3(8 ℛ_ν^2 + 62 ℛ_ν + 95)/15 (4 ℛ_ν + 5)τ^2 + τ^3 ,
    θ^I    = - Ψℛ_νβ_0 k^3(4 ℛ_ν + 15)/5 (4 ℛ_ν + 5)τ^2 + τ^3 ,
    σ_ν^I    = 2 Ψβ_0 k (4 ℛ_ν + 15)/5 (4 ℛ_ν + 5)τ
    	                + 6 H_0Ψ√(S_γν)ℛ_νℛ_bcβ_0 k/5 (4 ℛ_ν + 5)τ^2 + τ^3 ,
    
    	        h^I    = 0 + τ^3 ,
    η^I    = - 2 Ψℛ_νβ_0 k (4 ℛ_ν + 15)/5 (4 ℛ_ν + 5)τ
    	                + 3 H_0Ψ√(S_γν)ℛ_νℛ_bcβ_0 k/2 (4 ℛ_ν + 5)τ^2 + τ^3 .

    
    

 §.§ Vector modes

        Considering adiabatic perturbations in the scalar contributions, during TC and deep in the
        radiation era, the super-Hubble evolution is
        
            D^(0)_χ_γ    = D^(0)_χγ ,                                                 D^(1)_χ_γ    = 3 D^(0)_χγ H_0 √(S_γν) ℛ_b/4 (ℛ_ν- 1) , 

            D^(0)_χ_ν    = - Ψβ_0 (4 ℛ_ν+ 15)/10 ,          D^(1)_χ_ν    = 0 , 

            D^(0)_χ    = 0 ,                                                                            D^(1)_χ    = 0 ,
        
        
    D^(2)_χ_γ   = 3 D^(0)_χγ H_0^2 S_γνℛ_b (3 ℛ_b + ℛ_bc(ℛ_ν - 1))/16 (ℛ_ν - 1)^2 - Ψβ_0 k^2(8 ℛ_ν + 25)/60 ,
    
                D^(2)_χ_ν   = - Ψβ_0 k^2(8 ℛ_ν + 25)/60 ,
    
                D^(2)_χ   = 0 ,

        
    D^(3)_χ   = 0 ,
    
                D^(4)_χ   = 5 Ψβ_0 k^4/8 (8 ℛ_ν + 45) ,

        
    D^(0)_π_ν   = 0 ,
    
                D^(1)_π_ν   = 0 ,
    
                D^(2)_π_ν   = 0 ,
    
                D^(3)_π_ν   = Ψβ_0 k^2(32 ℛ_ν^2 + 268 ℛ_ν + 375)/270 (8 ℛ_ν + 45) ,

        
    D^(0)_S+Ḟ   = 0 ,
    
                D^(1)_S+Ḟ   = 0 ,
    
                D^(2)_S+Ḟ   = 5 Ψβ_0 k^2/8 ℛ_ν + 45 ,

        Again, imposing the physical requirement that photons and neutrinos had the same velocity 
        in the very early Universe, we are led to
        
    D^(0)_χ_γ = - Ψβ_0(4 ℛ_ν + 15)/10 .

        The initial conditions provided for the numerical integration are
        
    χ_γ   = - Ψβ_0(4 ℛ_ν + 15)/10
                     - 3 H_0Ψ√(S_γν)ℛ_b β_0(4 ℛ_ν + 15)/40 (ℛ_ν - 1)τ ,
           - 3 H_0^2Ψ S_γνℛ_b β_0(4 ℛ_ν + 15) (ℛ_νℛ_bc + 3 ℛ_b - ℛ_bc)/160 (ℛ_ν - 1)^2τ^2 - Ψβ_0 k^2(8 ℛ_ν + 25)/60τ^2 +τ^3 ,
    χ_ν   = - Ψβ_0(4 ℛ_ν + 15)/10 - Ψβ_0 k^2(8 ℛ_ν + 25)/60τ^2 + τ^3 ,
    χ   = 5 Ψβ_0 k^4/8 (8 ℛ_ν + 45)τ^4 + τ^5 ,
    π^V_ν   = Ψβ_0 k^2(4 ℛ_ν(8 ℛ_ν + 67) + 375)/270 (8 ℛ_ν + 45)τ^3 + τ^4 ,
    
                S+Ḟ   = 5 Ψβ_0 k^2/8 ℛ_ν + 45τ^2 +  τ^3 .
        
            


