General Perturbations of Homogeneous and Orthogonal Locally Rotationally Symmetric Class II Cosmologies with Applications to Dissipative Fluids

First order perturbations of homogeneous and hypersurface orthogonal LRS (Locally Rotationally Symmetric) class II cosmologies with a cosmological constant are considered in the framework of the 1+1+2 covariant decomposition of spacetime. The perturbations, which are for a general energy-momentum tensor, include scalar, vector and tensor modes and extend some previous works where matter was assumed to be a perfect fluid. Through a harmonic decomposition, the system of equations is then transformed to evolution equations in time and algebraic constraints. This result is then applied to dissipative one-component fluids, and on using the simplified acausal Eckart theory the system is reduced to two closed subsystems governed by four and eight harmonic coefficients for the odd and even sectors respectively. The system is also seen to close in a simplified causal theory. It is then demonstrated, within the Eckart theory, how vorticity can be generated from viscosity.

First order perturbations of homogeneous and hypersurface orthogonal LRS (Locally Rotationally Symmetric) class II cosmologies with a cosmological constant are considered in the framework of the 1+1+2 covariant decomposition of spacetime. The perturbations, which are for a general energymomentum tensor, include scalar, vector and tensor modes and extend some previous works where matter was assumed to be a perfect fluid. Through a harmonic decomposition, the system of equations is then transformed to evolution equations in time and algebraic constraints. This result is then applied to dissipative one-component fluids, and on using the simplified acausal Eckart theory the system is reduced to two closed subsystems governed by four and eight harmonic coefficients for the odd and even sectors respectively. The system is also seen to close in a simplified causal theory. It is then demonstrated, within the Eckart theory, how vorticity can be generated from viscosity.

I. INTRODUCTION
Even if most observations indicate that the universe on the large scale is homogeneous and isotropic, and to high accuracy described by the ΛCDM model [1][2][3], the upper limits on anisotropy from redshift measurements are rather high [4][5][6]. Hence it is of interest to see how different types of perturbations, like density fluctuations, gravitational waves, rotation of matter etc., evolve and interact on an anisotropic background. For some earlier works on perturbations in anisotropic universes see e.g. [7][8][9][10][11][12]. As a first generalization away from isotropy, without introducing the full complexity of general anisotropic models, we will in the following consider backgrounds with one direction of anisotropy. In some earlier papers [13][14][15] we studied perfect fluid perturbations of homogeneous and Locally Rotationally Symmetric (LRS) cosmologies of class II using a perturbation method based on the 1+1+2 covariant split of spacetime [16][17][18], which is an extension of the 1+3 covariant split [19][20][21][22][23].
In this work we extend the analysis to a general form of the energy-momentum tensor where u a is the 4-velocity of some timelike observer, µ is the energy density, p is the isotropic stress, π ab is the anisotropic stress, and q a is the energy flux relative u a . For their definitions see section II. Depending on the physical situation, u a can be specified in different ways. For example, for a one-component perfect fluid a natural choice would be the 4-velocity of the fluid, whereas for a one-component imperfect fluid one can choose either the restframe of the fluid or the system where the energy-flow q a vanishes. The 1+1+2 split will then be made with respect to u a and a spacelike unit vector n a , which coincides with the direction of anisotropy on the background. As variables, we use the kinematic quantities of the u a and n a vectors, the electric and magnetic parts of the Weyl tensor, and the energy-momentum tensor, where the latter, through Einstein's equations, is equivalent to the Ricci tensor and the cosmological constant Λ. These objects are then subject to the Ricci identities for u a and n a and the Bianchi identities. The covariant derivative is projected into three parts, one timelike derivative along u a , a spacelike derivative along n a , and one more along the perpendicular directions. These derivatives satisfy certain commutation relations, corresponding to the commutators of the basis vectors.
We will now use this formalism to study general first order perturbations on homogeneous, but anisotropic, cosmological models of LRS class II. For the perturbed quantities we use covariant variables which are zero on the background. In this way we avoid the gauge problem of identifying spacetime points on the background and the perturbed spacetime [24]. On applying the formalism, a system of first order ordinary differential equations in time and a system of first order partial differential equations are obtained for the background and perturbation variables respectively. After a harmonic decomposition of the perturbed system, it is transformed to a system of ordinary differential equations in time together with algebraic constraints. Some of the variables can then be solved for algebraically in terms of the others. After eliminating these, the system is reduced to a system of ordinary differential equations for some of the remaining variables, but in general some variables are still freely specifiable. The system can be made to close by imposing some specific physical theory. For example, for a barotropic perfect fluid the system is closed by imposing the equation of state p = p(µ).
After finding the general equations we will apply the formalism to one-component imperfect fluids. First we consider the simplified Eckart theory [25,26], which is acausal but should be a fair approximation when relaxations times are short compared to the typical macroscopic time scale, which in turn is set by the expansion rate of the universe. The Eckart theory also goes over into the standard theory for imperfect fluids in the non-relativistic limit. The system will close by imposing, apart from the equation of state as above, also the coefficients of bulk and dynamic viscosity and heat conductivity. In general, these are functions of two variables, say energy density µ and particle density N , which for the latter we also have to include an evolution equation, like the equation of continuity, to close the system. As an application, we then show how vorticity can be generated from the dissipative terms. For works on vorticity perturbations in isotropic universes see e.g. [27][28][29][30][31][32] , and for the generation of vorticity in N-body simulations see [33]. Finally, we show how the system can be closed also in a simplified causal theory [26].
The paper is organized as follows: In section II a short summary of the 1+3 and 1+1+2 covariant splits is given, and some basic elements of relativistic thermodynamics are collected in section III. The properties of the background metric are given in section IV, and the procedure for the perturbation theory together with the harmonic decomposition is summarized in section V. In section VI, the even and odd sectors of the reduced set of evolution equations for the perturbed quantities are presented, and in section VII the closed system obtained on imposing the Eckart theory is found. Similarly, the closed system for a simplified causal theory is given in section VIII. As an example of an application, the generation of vorticity in a dissipative fluid is discussed in section IX. In Appendix A some relations in the 1+1+2 formalism are collected. The commutation relations between the different differential operators are given in Appendix B, and properties of the harmonics in Appendix C. The linearized equations in the 1+1+2 formalism are given in Appendix D, and their harmonic expansions in Appendix E. Then, in Appendix F, 24 harmonic coefficients are solved for algebraically in terms of 11 coefficients determined from evolution equations and 9 freely specifiable coefficients.
II. THE 1 + 3 AND 1 + 1 + 2 COVARIANT SPLITS OF SPACETIME Our background will be characterized by two vector fields, one timelike 4-velocity u a with u a u a = −1 and one preferred spatial direction n a with n a n a = 1, and hence a 1+1+2 covariant split of spacetime, [16,17], will be useful. This is an extension of the 1+3 covariant split, which is with respect to a timelike vector field u a [19,20]. We here briefly summarize these formalisms and give the definitions of the quantities and derivativies used in the following text. For details the reader is referred to [16,17,19,20], and in appendix A some useful relations are given.
In the 1+3 formalism tensors are split into timelike and spatial parts with the projection operators U a b = −u a u b and h ab = g ab + u a u b respectively, where g ab is the 4metric. Similarly, the covariant derivative, ∇ a , is split into a covariant time derivative and a covariant projected spatial derivative aṡ and respectively.
Instead of the metric, we now use decompositions of the curvature tensor and the kinematic quantities of u a as variables. The curvature tensor is given in terms of the electric and magnetic parts of the Weyl tensor, E ab ≡ C acbd u c u d and H ab ≡ 1 2 ε ade C de bc u c respectively 1 , and by the Ricci tensor, which, through Einstein's equations, is given by the cosmological constant Λ and the energy-momentum tensor. In general the energymomentum tensor can be decomposed as where µ ≡ T ab u a u b is the energy density, p ≡ 1 3 T ab h ab is the isotropic stress, π ab is the anisotropic stress, satisfying π ab u b = 0 and π a a = 0, and q a is the energy flux relative u a , satisfying q a u a = 0. Due to the Bianchi identities, the energy-momentum tensor also satisfies ∇ b T ab = 0.
The covariant derivative of u a is given in terms of the kinematic quantities of u a as where we have introduced the acceleration A a ≡ u b ∇ b u a , the expansion Θ ≡ D a u a , the shear σ ab ≡ D a u b , and the vorticity ω ab ≡ D [a u b] , which also can be given in terms of the corresponding vector ω a ≡ 1 2 ε abc ω bc . Here T ab ≡ h c (a h d b) − 1 3 h ab h cd T cd has the properties T ab = T ba and h ab T ab = 0, and square brackets denote antisymmetrization, i.e. A [ab] = −A [ba] . The variables are now subject to integrability conditions given by the Ricci identity for u a , the Bianchi identities, and the commutators between the derivatives˙≡ u a ∇ a and D a , see e.g. [20], whereas Einstein's equations already are imposed by expressing the Ricci tensor in terms of the energy momentum tensor (3) and the cosmological constant Λ.
Similarly, when we have one more preferred vector field n a , which is assumed to be spacelike and normalized, we can make a further split into a covariant 1+1+2 formalism. This is then done with the projection operators n a n b and N b a = h b a − n a n b . In this way, 3-vectors are split into a scalar along n a and perpendicular 2-vectors, and 3-tensors are split into a scalar, a 2-vector and a 1 Here, the volume element on the 3-surfaces is given by symmetric and trace-free 2-tensor. The quantities from the 1+3 split then decompose in the following way as A a = An a + A a , ω a = Ωn a + Ω a , q a = Qn a + Q a and σ ab = Σ n a n b − 1 2 N ab +2Σ (a n b) +Σ ab . The electric and magnetic parts of the Weyl tensor, E ab and H ab , and the anisotropic part of the stress π ab are decomposed in the same way as the shear, in terms of the variables E, E a , E ab , H, H a , H ab and Π, Π a and Π ab respectively.
In a similar way which the covariant 4-derivative ∇ a was split, we can split D a into a derivative along n a and a projected derivative onto the 2-surfaces as and respectively.
On decomposing D a n b andṅ a analogously to (4), new kinematic quantities are introduced as andṅ where the acceleration of n a , a a ≡ n c D c n a , the 2-sheet expansion φ ≡ δ a n a , the twisting of the 2-sheets ξ ≡ 1 2 ε ab δ a n b , the shear of n a , ζ ab ≡ δ {a n b} , and α a ≡ N b aṅb . Here ε ab ≡ ε abc n c and curly brackets denote the projected, trace-free and symmetric part of a 2-tensor in an analogous way to the ab operation on 3-tensors. We also introduce the notation ψ a ≡ N b a ψ b for projection onto the 2-sheets.
The previous integrability conditions for the 1+3 formalism now have to be completed with the Ricci identities for n a . For a full 1+1+2 split of the Ricci identities and the Bianchi identities into evolution equations and constraints, see e.g. [16], where also commutation relations between the different differential operators introduced can be found.

III. THERMODYNAMICS
To get an accurate description of relativistic nonequilibrium or dissipative thermodynamics the model should be based on a detailed kinetic theory. For some covariant and causal fluid descriptions valid close to equilibrium see, e.g. [34][35][36][37]. Even these are quite complicated and we will later turn to some simplifications, see e.g. [26,38].
Recall that the energy momentum tensor is given by equation (3) as In the following we will take u a to be the 4-velocity of matter (particle frame) when dealing with thermodynamics. Another possible choice is the so called energy frame, which makes q a = 0. For a dissipative fluid p =p + p ζ B , wherep is the equilibrium pressure of the fluid and p ζ B represents the bulk viscosity. Similarly π ab is the shear viscosity and the heat flow is given by q a . Now introduce also the particle flux N a and entropy flux S a , given by and respectively, where N is the particle number density, S the specific entropy and R a a dissipative term satisfying R a u a = 0. In line with [36] we define the temperature T from the Gibbs relation which is subject to the integrability condition The time evolution of the temperature T is then given by [26] From (11) we also see that in general the equilibrium scalars µ,p, N , S and T can be given in terms of just two of them, say µ and N , asp(µ, N ), S(µ, N ) and T (µ, N ). Hence we will also need an evolution equation for N , given by requiring particle conservation The evolution equation for the energy density is contained in the general set of equations obtained from the Ricci and Bianchi identities. To get the theory consistent with the second law of thermodynamics one should also have which will put restrictions on the dissipative quantities p ζ B , π ab and q a .

A. The Eckart Theory
A simple covariant, but acausal, theory was given by Eckart [25,26]. It can give a good description when relaxation times are short compared to relevant time scales, such as the expansion rate of the universe.
With R a = q a , the heat flux, the expressions below are consistent with the second law of thermodynamics as stated in equation (15). The bulk viscosity is given by in terms of the coefficient of bulk viscosity, ζ B = ζ B (µ, N ). Similarly π ab , the shear viscosity, is given by in terms of the coefficient of dynamic viscosity, η = η(µ, N ). Finally, the heat flow is given by where κ = κ(µ, N ) is the coefficient of thermal conductivity. These are natural relativistic generalizations of the corresponding non-relativistic laws. The main difference from them is that the acceleration appears in the heat equation.

B. Simplified Causal Theory
A simplified, but causal, version of the theory developed in [36,37] can be obtained by neglecting couplings between viscosity and heat, as well as neglecting some terms which will be of higher order, see [26]. The equations (16)-(18) then become evolution equations, given by in terms of the relaxation times τ ζ , τ π and τ q , which in general are assumed to be functions of µ and N .

IV. BACKGROUND SPACETIMES
To see the effects of anisotropy without going to the most general case, we consider for the backgrounds a class of cosmological models with local rotational symmetry (LRS), which has the flat Friedmann model as a special case. These are characterized by two vector fields, the 4-velocity u a and the symmetry axis n a , and hence are suitable for a 1+1+2 split. The models are homogeneous and hypersurface orthogonal and belong to the LRS class II, see e.g. [14,39,40], which is characterized by a vanishing magnetic part of the Weyl tensor H ab , vanishing vorticity ω ab , and vanishing twist ξ of the 2-sheets perpendicular to n a . These sheets are maximally symmetric with 2-curvature scalar R = 2K/a 2 A closed system is obtained if explicit expressions for p and Π are given. In the context of a dissipative onecomponent fluid, this is equivalent to specifying the equilibrium pressurep, the bulk viscosity p ζ B , and shear viscosity Π, as functions of, e.g., the energy density µ and the number density N . The evolution equation (14) for N then has to be added to the system. The time evolution of the corresponding background temperature T is obtained from (13) aṡ For the case of a dissipative fluid, we should also note that since we are allowing for non-zero p ζ B and Π, the background spacetime will not describe a thermodynamic equilibrium. Since the Eckart theory and its generalizations can be seen as mainly being valid close to thermodynamic equilibrium [37], we should therefore keep in mind that we in principle ought to require that p ζ B and Π, although possibly non-zero, are small. On applying the Eckart theory on an expanding and anisotropic background, this is most naturally enforced by choosing suitably small coefficients of viscosity.

V. PERTURBATIONS
Perturbations on the backgrounds of section IV were studied in [13][14][15]. In these works the energy-momentum tensor was assumed to be that of a barotropic perfect fluid, both for the background and the perturbations.
Here these studies will be extended to the most general first order perturbations, so at this stage we do not assume anything about the quantities p, π ab and q a , but first later on specialize to dissipative fluids.
The covariant split of spacetime is advantageous for perturbative calculations. In this approach one uses covariant quantities which are zero on the background as the perturbed variables. Hence, according to the Stewart-Walker lemma [24], one avoids the gauge problem, i.e. of how to identify points on the perturbed spacetime with points on the background, see e.g. [41] for another method to construct gauge invariant objects and [27,[42][43][44] for more perturbation methods in cosmology.
The set of covariant variables we will use are the scalars which specify the background, G = {µ, p, E, Θ, Σ, Π}, and the other 1+1+2 covariant quantities defined in section II. The latter are zero on the background and hence are gauge-invariant and suitable to describe the perturbations. Later on, when specializing to dissipative fluids, we will also introduce the scalars N and T , which as well as the elements in G are nonzero also to zeroth order. When going to first order, the elements in G can be expanded as G = G (0) + G (1) , where G (0) are the quantities to zeroth order. For example µ = µ (0) + µ (1) etc. A convenient and covariant way to represent the G (1) 's is by using the gradients of the G's, i.e. by δ a G andĜ, which are zero on the homogeneous background. On performing a harmonic decomposition of the first order quantities, it is seen that theĜ derivatives may be expressed in terms of the G a ≡ δ a G [13][14][15]. Hence, the perturbations of the background quantities will be represented by the G a The complete set of first order quantities which vanish on the background is then given by When specializing to dissipative fluids, these quantities are supplemented with the first order quantities So far there is a freedom in the choice of frame in the perturbed spacetime. For a one component fluid one could, for instance, let u a be the 4-velocity of matter. Nor is the direction n a of anisotropy well-defined in the perturbed spacetime, and we will later restrict the frame by requiring a a = 0. See [13] for a discussion on how to fix the frame.
On imposing the Ricci identities for the vectors u a and n a and the Bianchi identities together with the commutator relations (see appendix B), evolution equations along u a , propagation equations along n a , and constraints are obtained for the covariant quantities. The exact equations are found in [16]. The linearized equations for the first order quantities are given in Appendix D.
The linearized partial differential equations are then transformed into ordinary differential equations by a harmonic decomposition in terms of eigenfunctions P k = P k (z) and Q k ⊥ = Q k ⊥ (ϑ, ϕ), where k are dimensionless comoving wave numbers in the n a direction and k ⊥ are the dimensionless comoving wavenumbers along the 2sheets. For their definitions and properties see Appendix C. Scalars are decomposed as where Ψ S k k ⊥ = Ψ S k k ⊥ (t). From the scalar harmonics Q k ⊥ even and odd vector and tensor harmonics are defined, see Appendix C. In terms of them, vectors are expanded as and tensors similarly as where all Ψ V,T k k ⊥ and Ψ V,T k k ⊥ are functions of time. We can now see how the hat derivatives,Ĝ, of the elements G may be expressed in terms of the G a . Expand G and G a as and respectively. From the even part of the commutation relation (B3), the coefficients ofĜ are then given in terms of the coefficients of G a as Some other simplifying relations are also easy to find. From the commutation relation (B4) and the properties of the harmonics in appendix C it follows that the odd coefficients of G a satisfy and the degrees of freedom of the vorticity can also be decreased by substitution of the harmonic expansions of Ω and Ω a in propagation equation (D25), which gives In the following sections and the appendices we will drop the indices k k ⊥ on harmonic coefficients, since their meaning will be obvious due to the superscripts S, V and T . The harmonically decomposed system obtained from the first order system in Appendix D can be found in Appendix E. It is given in terms of evolution equations and constraints, and separates into an even and an odd sector. By solving the constraints for some of the harmonic coefficients in terms of the others and then plugging these back into the differential equations, new constraints or identities are obtained. The procedure is continued until no more constraints are obtained. In this way each sector can be reduced to evolution equations for a minimal set of harmonic coefficients, in terms of which the rest of the coefficients, together with some freely specifiable coefficients, are given algebraically, see Appendix F.

VI. EVOLUTION EQUATIONS
In this section we give the reduced systems, for the even and odd sectors respectively, before imposing the thermodynamic conditions from section III. The systems obtained on imposing the Eckart theory are given in the following section VII and the systems for a simplified causal thermodynamic theory are found in section VIII. The frame has been partially locked by requiring a a = 0. Later, when specializing to an imperfect fluid, we will choose u a to coincide with the 4-velocity of the fluid.

A. Even parity sector
Here the system is reduced to evolution equations for V V , X V and φ S are expressible in terms of the other coefficients, see section F 1.
where the abbreviations have been used.

B. Odd parity sector
In the odd sector we obtain evolution equations for freely specifiable functions and the remaining quantities Ω V , p V and µ V are expressible in terms of the others, see section F 2.

VII. EVOLUTION EQUATIONS USING THE ECKART THEORY
In this section we assume the simplified Eckart theory, as described in section III A. Hence from equations (16), (17) and (18) we have in terms of the temperature T = T (µ, N ) and the coefficients of bulk viscosity ζ B = ζ B (µ, N ), dynamic viscosity, η = η(µ, N ), and conductivity κ = κ(µ, N ). The isotropic stress is given by p =p + p ζ B where the equilibrium pressurep =p(µ, N ). These equations also have to be supplemented with the equation of continuity (14) for the particle density N We saw in section VI that the harmonic coefficients p V , A S , A V , Π V , Π T and Y V in the even sector and A V , Π V and Π T in the odd sector, corresponding to the quantities p a =p a + p ζa , A, A a , Π a , Π ab and Y a ≡ δ a Π, were freely specifiable. By imposing equations (16), (17) and (18), we can solve for these quantities in terms of the others as where we have introduced the first order variables and the coefficients We also have the zeroth order equality From equation (14) and the commutator between the dot derivative and δ a , one obtains the projected time evolu- The corresponding even parity harmonic coefficients of the quantities p a , A and A a , Π a and Π ab , Y a and T a are then given by (67), (68), (69), (70) and (71) respectively as whereas the evolution equation for the even component of Z a is given bẏ Similarly, the odd parity coefficients A V , Π V and Π T are given by The closed systems of evolution equations for the even and odd sectors respectively are given below.

Even sector
For the even sector, the following closed system of evolution equations for the first order coefficients

Odd sector
Similarly, for the odd sector we obtain a closed system for the coefficients Ω S , where

VIII. CAUSAL THEORIES
In the previous section we saw how the two systems closed in the simplified acausal Eckart theory. We will here see how the system is changed from the point of view of closure when introducing a causal theory. Note that the main difference is that the Eckart equations (16)- (18) are transformed into evolution equations, as (19)-(21) for the simplified theory given in III B. Since this structure of the equations is kept in the slightly more elaborate theories in [26,38], it is sufficient to consider the simplified theory given by (19)-(21) for the purpose of deciding if the systems close, provided all introduced thermodynamical coefficients are given as functions of, e.g., µ and N .
On making a 1+2 split of the equations (19)-(21) we first note that they give two evolution equations for the zeroth order scalars p ζ B and Π respectively. Hence there are two more dynamical variables already to zeroth order. Similarly as with the other zeroth order variables we now introduce the 2-gradient of p ζ B (Y a ≡ δ a Π was already introduced in section V) as a new first order variable From p =p + p ζ B one then has that where 0 and 1 are defined in (72). By acting with δ a on equations (105) and (106) and using the commutation relation (B2), the evolution equations for p ζa and Y a are then obtained as where (111) From (20) also the following evolution equations for Π a and Π ab are obtained. The heat equation (21), splits as where T a is defined from (71) and (72). Finally, note that the evolution equations (14) and (74) for N and Z a ≡ δ a N still apply. Since we already have evolution equations for Q and Q a , or equivalently for the corresponding harmonic coefficients Q S , Q V and Q V in the general systems from section VI, equations (114) and (115) will be used to solve for A and A a . The corresponding harmonic coefficients A S , A V and A V are then obtained from and

Even sector
The even sector is now supplemented with the following evolution equations for Π and from which p V is given algebraically by

Odd sector
In the odd sector we obtain the following two evolution equations for Π V and Π and respectively. Y V , Z V and p V ζ are as before given by (41) as and corresponding to whereṗ =ṗ +ṗ ζ B . The corresponding differential equations for p V ζ and Y V , obtained from (109) and (110), will be identically satisfied. By adding the evolution equations for Π V , Π T , Y V , Z V , p V ζ and Π V , Π T to the even and odd general systems in section VI respectively, these are now seen to close.

IX. VORTICITY
From the general evolution equations (43) and (59) for the vorticity we see that the acceleration coefficients act as source terms. For a perfect fluid with a barotropic equation of state, p = p(µ), these terms can be rewritten in terms of the vorticity. Hence, vorticity cannot be generated, a result which holds to all orders in perturbation theory, see e.g. [20,28,45]. However, on applying the Eckart theory for a dissipative fluid, it can be seen from (84) and (97) that vorticity could be generated by heat flows and, as will be seen below, indirectly also by viscosity.
Thus, as an application of our equations, we now more thoroughly investigate the possibility of vorticity generation due to heat flows and viscous stresses. For this purpose, we will focus on the Eckart system. Although it has been argued that the Eckart theory could be valid for normal materials within the limits of applicability of the hydrodynamic approximation, at least in certain cases, it is in general found to be plagued by severe instabilities when considering perturbations away from thermal equilibrium [37,46,47].
Since our main goal is to show that the inclusion of viscous stresses and heat flows could indeed provide a mechanism for generating fluid vorticity, we will not worry too much about the possible instabilities of the Eckart theory at this point. A more thorough investigation of the stability conditions of the Eckart theory on a homogeneous and anisotropic cosmological background could be interesting, but is outside the scope of this paper. Instead, when performing numerical calculations, we will specify initial conditions and parameter values as to not induce any noticeable instabilities. Although the numerical values may be unphysical, the solutions could still highlight the sought mechanisms.

A. The Importance of Dissipative Effects
Before turning to numerics, we investigate the importance of the dissipative effects analytically. As a first step, an additional time derivative is applied to equation (84), which results in the following relation, where and Here we have defined the gradient of the shear viscosity and the square of the norm of the physical wave vector while A V is determined by equations (77) and (80). Looking at the source term S in (134), we can note the importance of the shear viscosity, as this source will vanish on taking the limit of vanishing shear viscosity, i.e. η, ς 0 , ς 1 → 0. In this limit, Ω V = 0 is a solution to equation (131) with the initial conditions Ω V (t 0 ) = Ω V (t 0 ) = 0. Provided that these initial conditions are satisfied, for example by specifying Q S (t 0 ) = Q V (t 0 ) = Ω V (t 0 ) = 0 as can be seen from (84), we would therefore expect that it is possible to have no vorticity generation. This absence of vorticity generation when setting η = 0 is due to that the curl-like combination forms a closed subsystem together with Ω V . This can be seen on studying the evolution equation for curlQ for η = 0. Similarly, the evolution equation for Ω V can be written aṡ for all values of η. From the closed system constituted by equations (139) and (140), it can then be seen that Hence, it is interesting to note that the generation of vorticity to first order with the Eckart theory, as illustrated through the source term in (131), is crucially dependent on the fact that the background spacetime is both anisotropic and not in thermal equilibrium, so that the combination ηΣ = −Π/2 is a non-vanishing quantity to zeroth order. However, it should also be noted that, even if η = 0, vorticity can be generated by a non-zero curlQ at t 0 , as curlQ(t 0 ) = 0 implies thatΩ In a similar manner as for Ω V , we can derive a homogeneous second order differential equation for the scalar part Ω S in the odd sector where and Hence, with the initial conditions Q V (t 0 ) = Ω S (t 0 ) = 0, which via equation (97) imply thatΩ S (t 0 ) = 0, a nonzero component Ω S cannot be created even if η = 0. For a discussion of a similar connection between vorticity and heatflow in the context of axially symmetric systems with dissipation, see [48,49].

B. Numerical results
To get even more tangible results, we now perform some simple numerical calculations. In doing so, we fixate the remaining length and temperature scales by setting Λ = 1, and k B = 1, where k B is Boltzmann's constant. Hence, all quantities are treated as being dimensionless in the following.
The numerical calculations are performed using ode45 in MATLAB to solve the even and odd sectors, given by (84)-(91) and (97)-(100) respectively, together with the background equations from section IV. To avoid some numerical issues due to subtraction of very similar values, we determine the evolution of the vorticity Ω V numerically by writing equation (84) as in (140), adding (138) to the system as an additional equation. In principle, one could use (137) to eliminate either Q V or Q S , rendering either (85) or (89) unnecessary. However, for simplicity, we solve all three equations (85), (89), and (138), using (137) to set consistent initial conditions. Any deviations between the left and right hand sides in (137) should then be due to numerical errors.
In addition to the aforementioned equations, it is also necessary to specify suitable equations of state to get a closed and numerically solvable system, For this purpose, we use that the temperature, T = T (µ, N ), and equilibrium pressure,p =p(µ, N ), should satisfy the integrability condition (12). A possible choice is then to let the pressure satisfy the ideal gas law, i.e.
and the energy density to be given by a rest-energy term N m and a thermal energy f N T /2, so that where m is the particle mass and f is the number of degrees of freedom. From these equations we can solve for T = T (µ, N ) andp =p(µ, N ) as and where γ − 1 = 2/f . With the basic thermodynamic quantities specified, we also need models for the dissipative coefficients of conductivity and viscosity. For this purpose, we assume that which holds quite generally for normal gases [50]. As for the bulk viscosity, we assume ζ B = 0 for simplicity, as the results from the previous section suggest that ζ B is less important than η when considering vorticity generation. On performing some preliminary numerical calculations using the system specified above, we note that the instabilities of the Eckart theory also seem to emerge here. These instabilities are encountered when increasing the wave numbers for a fixed set of initial conditions and fluid properties, or when the fixed set of parameters is not specified carefully to begin with. A possible explanation for some of the latter type of instabilities can be seen when considering the expressions for the accelerations, which couple into many of the evolution equations. In the Eckart theory, these expressions involve terms behaving as ∼ Q/(κT ). Therefore, a problem seems to arise when κT becomes very small. This specific problem could possibly be solved when using the simplified causal theory, since κT does not appear alone in the denominator when solving equations (116)-(118) for the accelerations.
To avoid the instability problems we will, as mentioned previously, specify parameters and initial conditions as to not induce any noticeable instabilities, even though the values may not be physical. In doing so, we assume that the order of magnitude of the variables on the background manifold do not vary drastically. The background spacetime will thus, for our purposes, be taken to be the one determined by the following parameters and initial conditions which ends at an expanding de Sitter solution.
Equipped with a background spacetime, we move on to investigate the effect of the shear viscosity when considering the first order perturbations. For this purpose, we focus on the even sector, as the results from the previous section suggest that the viscosity is more important for the vorticity generation in this sector. On assuming that all perturbations except Σ T vanish initially, we get the results shown in Fig. 1. In Fig. 1, it can be seen that the vorticity then reaches a maximum amplitude that is approximately three orders of magnitude larger than the initial shear perturbation. Thus, it would seem as though we have some new mechanism for creating vorticity. However, the precise quantitative results are questionable due to the unphysical nature of the chosen parameters and initial conditions.
It can be argued that the initial vorticity perturbation in Fig. 1 is actually generated due to numerical errors in the calculations, which were performed with a relative and absolute tolerance of 10 −9 and 10 −34 respectively. Although this might be the case, the dissipative effects can still be seen to play an important role. We can illustrate this importance by performing the same calculations as before, but setting η(t 0 ) = 0, which implies η = 0. The numerical calculations then show that both curlQ and Ω V remain identically zero throughout the whole simulation. This agrees with the differential equations from the previous section, as our choice of initial conditions imply that both Ω V andΩ V are set to zero initially through equation (84). With η put to zero, equation (131) becomes homogeneous, so with the given initial conditions Ω V should remain identically zero.
Finally, it should be noted that a non-zero vorticity Ω V is obtained numerically when specifying curlQ(t 0 ) = 0, even if η = 0. This is in agreement with the discussion below (140).

C. Perfect fluids with non barotropic equations of state
In this section we close the general first order system in section VI by requiring a perfect fluid form of the energy momentum tensor, but with a non barotropic equation of statep so that The perfect fluid form is obtained by setting the dissipative terms Q V , Q S , Π, Π V , Π T and Y V in the even sector and Q V , Π V and Π T in the odd sector, respectively, all to zero. From equations (48) and (49) in the even sector and equation (58) in the odd sector we can now solve for the components of the acceleration as and respectively. Substitution into the evaluation equations (43) and (59) for Ω V and Ω S now giveṡ Hence there are no source terms in the equations, so that vorticity cannot be generated to first order in perturbation theory. This result was also obtained in [29,30] for perturbations on isotropic Friedmann backgrounds. They also found that vorticity can be generated to second order if the density and entropy gradients are not parallel. Similarly, starting from the exact evolution equation for the vorticity in the 1+3 covariant split, [20], an evolution equation for ω a which holds to all ordered can be derived. Together with the twice contracted Bianchi identities for a perfect fluid and the commutator for a scalar Ψ, one then obtains Hence vorticity can be generated if the last term is nonzero, which can happen if 1 = 0 and the density and number density gradients are non parallel. This term is of second order if both gradients are considered small. For an example of how second order perturbation theory can be done using the covariant and gauge invariant approach, see e.g. [52].

X. CONCLUSIONS
In this paper we have presented results for general first order perturbations on homogeneous and hypersurface orthogonal LRS class II cosmologies. In doing so, we have extended previous works to general energy-momentum tensors, including effects from the energy flows q a and anisotropic pressures π ab . Using the 1 + 1 + 2 covariant formalism and harmonically decomposing the variables, we have shown that, after specifying the details of the energy-momentum tensor, it is possible to obtain closed systems of ordinary differential equations in time. Due to the generality of the results from section VI, these can be combined with any detailed description of the energy momentum tensor. As an example, we have shown how the equations can be used to describe one-component dissipative fluids after imposing a suitable thermodynamic theory to close the system. In an upcoming article, we will also describe how the results from section VI can be used to study a cosmological plasma and its associated electromagnetic fields [53].
When specializing the theory to dissipative onecomponent fluids, we have observed possible mechanisms for generating fluid vorticity to first order in the perturbations, which was seen to be a result of the inclusion of viscosity, heat flows, and a background anisotropy. This is in contrast to the barotropic perfect fluid case, where it is well known that fluid vorticity cannot be generated. However, when studying the standard Eckart theory numerically, we encountered some problems with stability. Thus, the commonly observed problems of the Eckart theory also seem to plague the theory in this context, at least with our choice of equations of state and models for the coefficients of thermal conductivity and viscosity. A more thorough stability analysis could be interesting, but our numerical observations seem to favor the idea that the Eckart theory ought to be neglected in favor of causal and stable descriptions. Here we state some relations in the 1+3 and 1+1+2 covariant splits of spacetime which are useful when deriving the equations. When deriving the basic equations from the Ricci and Bianchi identities, one often encounters the following covariant derivatives for the preferred vector fields n a and u a , the projection tensor N ab , the 2volume element ab , general scalars Ψ, 2-vectors Ψ a , and 2-tensors Ψ ab .

Appendix B: Commutation Relations
For a zeroth-order scalar field Ψ on orthogonal and homogenous LRS class II backgrounds, with φ = 0, the following first order commutation relations hold: For a first-order 2-vector Ψ a the following commutation relations, to first order, hold: and for a first-order trace-free and symmetric 2-tensor Ψ ab it holds that: Appendix C: Harmonics The properties of the harmonics associated with the background metric (22) where the scale factors a 1 ≡ a 1 (t) and a 2 ≡ a 2 (t), are here briefly summarized. For more details the reader is referred to [14,16,17,54,55].
First order scalars are decomposed as where the harmonic coefficients Ψ S k k ⊥ are functions of time. The eigenfunctions P k , which may be represented with e ik z , satisfy in terms of the dimensionless comoving wave numbers k along the n a direction, whence the physical wavenumbers are given by k /a 1 . The eigenfunctions Q k ⊥ satisfy: where a 2 is the scale factor of the 2-sheets, and k ⊥ are the dimensionless comoving wavenumbers along the 2-sheets. When R, the 2-curvature, is positive the 2-sheets are spheres and the eigenfunctions can be represented by the spherical harmonics Y m l with k 2 ⊥ = l(l + 1), l = 0, 1, 2, .... Due to the background symmetry the m-values will not appear explicitly in the equations. When R ≤ 0, and the 2-sheets are open, the k ⊥ take continuous values. For continuous values the sums go over into integrals.
Vectors Ψ a are expanded in terms of the even and odd vector harmonics [17,54,55] as Similarly, a tensor Ψ ab can be expanded in terms of the even and odd tensor harmonics Some useful relations between the harmonics are listed below. The vector harmonics satisfy the following orthogonality relation and the even and odd harmonics are related through The harmonics also satisfy the following differential rela-tionsQ For the tensor harmonics the orthogonality relation is the relation between even and odd harmonics are and the differential relations arė The vector and tensor harmonics, Q k ⊥ a , Q k ⊥ a and Q k ⊥ ab , Q k ⊥ ab defined in (C4) and (C6) respectively are defined with factors a 2 and a 2 2 . However, since the derivative δ a picks out factors of order k ⊥ /a 2 , an alternative could be to define vector and tensor harmonics as and to keep their norm to the order of unity. Correspondingly the vector and tensor harmonic coefficients change by factors k ⊥ and k 2 ⊥ respectively as and analogous for the odd coefficients. Hence, in equations with a mixture of different types of objects (scalars, vectors, tensors), different factors of k ⊥ have to be introduced to get the right relative magnitudes in powers of k between the terms. This can be of importance, for example if one want to study the equations in the limit of large wavenumbers k, the so called optical limit. The harmonic expansion of the linearized equations given in appendix E is easily transfered to the one in the basis (C21) and (C22) by noting that each factor a 2 should be associated with a factor 1/k ⊥ . For example, in equation (E1) the term k 2 ⊥ a2 α V should be replaced by k ⊥ a2 α V (on dropping the stars * ).

Appendix D: Linearized Equations
The evolution equations arė The equations containing a mixture of evolution and propagation contributions are The equations containing only propagation contributions are (D30) Lastly, the constraints are

Appendix E: Harmonic Expansion of Linearized Equations
Here the harmonic expansion of the first order equations given in section D is presented.

Even Parity
The evolution equations for the scalars arė for the 2-vectors areḢ and for the 2-tensors arė The constraint for the scalar is for the 2-vectors are and for the 2-tensors are

Odd Parity
The evolution equations for the scalars are for the 2-vectors arė and for the 2-tensors arė The constraints for the scalars are ik a 1 ξ S = Σ + Θ 3 and for the 2-tensors are

Harmonic decomposition of equations for imperfect fluids
Harmonic decomposition of Eckart's generalized equation τ qq a + q a = −κ (D a T + T A a ) , Similarly τ ππ ab + π ab = −2ησ ab (E62) gives and From (E62) one also obtains the zeroth order equation for Π. The harmonic coefficients of its gradient Y a ≡ δ a Π satisfies and respectively. From the bulk viscosity equation follows that the harmonic coefficients of its gradient p ζa ≡ δ a p ζ B satisfies and respectively. Similarly, from the continuity equation the harmonic coefficients of Z a ≡ δ a N are obtained froṁ and respectively. Finally the harmonic coefficients of the gradient of T = T (µ, N ), T a ≡ δ a T , are given by and respectively. The coefficients above are given by and τ ζµ ≡ ∂τ ζ ∂µ , τ ζN ≡ ∂τ ζ ∂N , τ πµ ≡ ∂τ π ∂µ , τ πN ≡ ∂τ π ∂N . (E77)

Appendix F: Harmonic Coefficients
From the equations in appendix E we can now solve for the even parity harmonics ζ T , E V , H V , Σ V , α V , W V , V V , X V and φ S in terms of the seven even parity independent coefficients and the freely specifiable function p V , A V , A S , Π V , Π T and Y V . Similarly we solve for the odd parity harmonics ζ T , p V and µ V in terms of the four odd parity independent coefficients and the freely specifiable functions A V , Π V and Π T .

Even parity
The even harmonic coefficients can be expressed as ik a 1 a 2