On the Interaction Between Electromagnetic, Gravitational, and Plasma Related Perturbations on LRS Class II Spacetimes

We investigate electromagnetic, gravitational, and plasma related perturbations to first order on homogeneous and hypersurface orthogonal locally rotationally symmetric (LRS) class II spacetimes. Due to the anisotropic nature of the studied backgrounds, we are able to include a non-zero magnetic field to zeroth order. As a result of this inclusion, we find interesting interactions between the electromagnetic and gravitational variables already to first order in the perturbations. The equations governing these perturbations are found by using the Ricci identities, the Bianchi identities, Einstein's field equations, Maxwell's equations, particle conservation, and a form of energy-momentum conservation for the plasma components. Using a $1+1+2$ covariant split of spacetime, the studied quantities and equations are decomposed with respect to the preferred directions on the background spacetimes. After linearizing the decomposed equations around a LRS background, performing a harmonic decomposition, and imposing the cold magnetohydrodynamic (MHD) limit with a finite electrical resistivity, the system is then reduced to a set of ordinary differential equations in time and some constraints. On solving for some of the harmonic coefficients in terms of the others, the system is found to decouple into two closed and independent subsectors. Through numerical calculations, we then observe some mechanisms for generating magnetic field perturbations, showing some traits similar to previous works using FLRW backgrounds. Furthermore, beat-like patterns are observed in the short wave length limit due to interference between gravitational waves and plasmonic modes.


I. INTRODUCTION
Using the standard model of cosmology, based on the homogeneous and isotropic FLRW spacetimes, most astronomical observations on a cosmological scale can be explained. However, some unsolved mysteries remain. One of these mysteries pertain to the origin of the large scale magnetic fields that can be observed on the level of galaxies and clusters of galaxies. Although several mechanisms have been proposed for generating and amplifying these magnetic fields, the topic has seemingly been rather controversial [1]. Amidst the controversy, a common explanation appears to be that the magnetic fields were originally created very early in the history of the universe. It is then suggested that these primordial seed fields have been amplified later on through galactic dynamo mechanisms, creating the fields we observe today [1]. Although the galactic dynamo could greatly amplify the seeds, the proposed mechanisms that generate the seeds may be too weak, or, in some cases, lie uncomfortably close to the minimum strength requirement. In addition to the strength requirement, there are also stringent requirements on the coherence length of the seeds, further limiting the possible candidates 1 .
To increase the number of possible candidates, one way around the stringent requirements would be if some mechanism existed for amplifying the seeds after they have been created, but acting before the galactic dynamo starts. Some possible mechanisms of this kind have already been proposed. Examples include mechanisms based on gravitational wave interactions [2][3][4], and velocity perturbations in the cosmological plasma [4]. Although some of the mechanisms for amplifying the seed fields have been suggested to be too weak [5][6][7], the prospect of a mechanism based on the dynamics of a self-gravitating plasma and classical general relativity is still enticing, since it would not explicitly depend on any additional exotic physics. Therefore, it is still of great interest to study the interactions between a cosmological plasma and general relativistic effects.
The studies of self-gravitating plasmas in a cosmological scenario seem to have mainly used a perturbative approach with the FLRW spacetimes as zeroth order backgrounds [2][3][4][5][6][7][8]. However, due to the isotropy of the FLRW spacetimes, it is essentially impossible to have a non-zero magnetic field 3-vector on the background, since that field would otherwise define a preferred spatial direction. To circumvent this problem, several different approaches have been used. One approach, which is usually referred to as the weak-field approximation, is to assume that the background magnetic field is sufficiently small so that it does not contribute to the energy-momentum tensor, which involves squares of the magnetic field [2]. With this assumption, the magnetic field should not affect the underlying geometry, and hence, in a sense, not violate the isotropy of the background spacetime. However, this approach runs in to some formal problems, since the generated magnetic field perturbations that are explicitly studied are not gauge invariant in the strict mathematical sense of the Stewart-Walker lemma [3,9].
To better avoid interference from gauge related degrees of freedom, second order schemes have also been used. The magnetic field is then introduced as a first order perturbation rather than on the background [3,4,7]. The interaction between the gravitational effects and the electromagnetic fields then becomes a second order quantity, since this involves products of the electromagnetic fields and other first order variables. However, some care is still needed, as it is not strictly correct to integrate the second order results for the gauge invariant variables to obtain explicit expressions for the generated magnetic field, since this quantity is not gauge invariant [7]. Hence, one is restricted to interpreting various gauge invariant quantities to draw any conclusion about the generated magnetic fields [7].
Finally, a more drastic approach, which is the one that will be used here, is to use an anisotropic background model instead of the FLRW spacetimes. Then there is no longer any isotropy problem, and a non-zero magnetic field to zeroth order is allowed without causing severe problems with the background geometry. In this paper, the aim is therefore to investigate if it is possible to get interesting gauge invariant interactions already to first order in the perturbations if we use a certain class of anisotropic backgrounds instead of the FLRW spacetimes. More specifically, we will consider backgrounds that are homogeneous, hypersurface orthogonal, and belong to the locally rotationally symmetric (LRS) class II of spacetimes, which allow for a preferred spatial direction.
First order perturbations on LRS class II spacetimes have been studied previously, both in the perfect fluid case [10][11][12][13], and when including dissipative effects [14]. However, electromagnetic perturbations on general LRS class II spacetimes have, to our knowledge, mainly been considered as test fields, without any electromagnetic fields on the background [15,16]. Hence, it is still interesting to investigate electromagnetic perturbations on LRS class II spacetimes when including more interactions with gravitational and matter related perturbations.
Our approach to study these interactions will follow the same steps as in [14], but instead of specializing the equations to describe one-component dissipative fluids, we will consider self-gravitating plasmas. Thus, we will make use of the results for general energy-momentum tensors presented in [14], but we will then close the general system by imposing Maxwell's equations, particle conservation, and equations describing energy-momentum conservation for the plasma components. In doing so, this article reproduces the methods and equations from [17] and an upcoming contribution to the proceedings from the sixteenth Marcel Grossmann meeting on general relativity [18], but extends these accounts with additional numerical results and discussions.

II. SPACETIME DYNAMICS
To determine the behavior of both the zeroth order background and the first order perturbations, we will follow [19] and make use of the Ricci identities, the Bianchi identities, and Einstein's field equations. In doing so, we will assume that the studied spacetimes have two preferred directions as specified by the vector fields u a and n a . The field u a is assumed to be timelike and normalized as u a u a = −1, whilst n a is spacelike with n a n a = 1 2 . As for the physical interpretation of these directions, u a denotes the 4-velocity of some fundamental observer, whilst n a defines the direction with respect to which the backgrounds are assumed to be locally rotationally symmetric. It is then natural to make use of the Ricci identities for these preferred vector fields, Together with the contracted Bianchi identities and Einstein's field equations these relations provide most of the necessary equations. Here G ab ≡ R ab − g ab R/2 is the usual Einstein tensor, whilst T ab is the energy-momentum tensor, which could be general at this point. However, with a general energy-momentum tensor, these equations will not be enough to get a closed system [14]. Therefore, to arrive at a closed system, we will here specialize the energy-momentum tensor to describe a plasma together with some electromagnetic fields. Hence we assume that where T (F ) ab is due to the plasma whilst T (EM ) ab comes from the macroscopic electromagnetic fields. The electromagnetic part can in turn be written as in terms of the Faraday tensor F ab , which satisfies the usual Maxwell's equations where j a is the 4-current density.
As for the plasma contribution, we will follow the formalism in [20] and assume that this contribution consists of separate contributions from each plasma component, so that where the ith component satisfies the following energy-momentum conservation equation Here I (i) a includes interactions such as collisions between the different fluid components, whilst j (i) a is the 4-current density of the ith component. The current density satisfies the relations where ρ (i) = q c(i) N (i) is the charge density, q c(i) is the charge, N (i) is the number density, and u (i) a is the the 4-velocity of the ith component. In the following, we will also write the interaction terms and the total current density as where J a and f (i) a are orthogonal to u a . Furthermore, we will assume that each fluid component can be described as a perfect fluid, so that where µ (i) and p (i) are the energy density and pressure relative to the component rest frame, as defined by the 4-velocity u (i) a . This 4-velocity can in turn be written as where We then also impose particle conservation for each fluid component by requiring that Finally, after specifying a microscopic description of I (i) a , equations of state p (i) = p (i) (µ (i) , N (i) ), and choosing a specific frame, we now have enough equations to get a closed system. To get a simple, yet still informative, system, we will later on assume that the plasma is sufficiently cold so that p (i) = 0. Then we will also assume that the interactions between the fluid components can be described with a single scalar electrical resistivity, η, in the magnetohydrodynamic (MHD) approximation. In this approximation, the plasma is described in terms of some collective variables instead of the individual component variables. Therefore, after linearizing and harmonically decomposing the equations, we will write them in terms of the harmonic coefficients of the following total fluid energy density and average 3-velocity Using Ohm's law, the total 3-current density can then be written as where E a and B a are the electric and magnetic fields defined through However, the cold MHD assumption will only be imposed after linearizing and harmonically decomposing the equations. Until then, we will continue using the more general multifluid description with non-zero p (i) and with interactions and current densities determined by Eqs. (12)- (16) .

III. COVARIANT SPLITS OF SPACETIME
As the studied spacetimes are assumed to have two preferred directions, it is natural to decompose the studied quantities with respect to those directions. Since the vector fields defining these directions are physical objects, the decompositions respect the tensorial properties of the decomposed variables. Hence they are usually referred to as covariant splits of spacetime.

A. 1 + 3
Starting with a 1 + 3 split with respect to the timelike vector field u a , we write the metric as where U ab = −u a u b is a projection tensor along u a whilst h ab projects onto the 3-dimensional hypersurfaces orthogonal to u a [21]. Using Eq. (25), all relevant objects can then be decomposed relative to u a and the 3-dimensional hypersurfaces. Some of the quantities from Sec. II have already been decomposed in this manner, such as the component 4-velocities, the total current density, and the Faraday tensor. However, it remains to fully decompose the total energy-momentum tensor and the Riemann tensor. Performing this decomposition we encounter the following projected variables where µ is the energy density, q a is the energy flux orthogonal to u a , p is the isotropic pressure, π ab is the trace-free anisotropic pressure, whilst E ab and H ab are the electric and magnetic parts of the Weyl tensor C abcd [21]. When defining the anisotropic pressure π ab , we have introduced a notation with angular brackets to denote the projected, symmetric, and trace-free (PSTF) part of a tensor relative to the projection tensor h ab .
We also obtain some new projected variables when considering the covariant derivative of the preferred vector field. Defining the time derivativeṪ and the spatial derivative fully projected onto the 3-dimensional hypersurfaces we can write the covariant derivative of u a as Here Θ ≡ D a u a is the expansion rate, ω ab ≡ D [a u b] is the vorticity, and σ ab ≡ D a u b is the rate of shear [21]. Instead of working directly with the vorticity tensor, we will write it in terms of the vorticity vector After having decomposed all relevant variables with respect to u a , we can go one step further and also decompose the 3-dimensional hypersurfaces with respect to n a , leading to a so-called 1 + 1 + 2 covariant split [19]. For this purpose, we write the projection tensor h ab as h ab = N ab + n a n b .
where N ab projects onto the 2-dimensional sheets orthogonal to both u a and n a . Using Eq. (36), the relevant 3-vectors can then be decomposed as Occasionally, we will use a notation with a bar over an index to denote a projection with N ab , so that ψā = N b a ψ b . As for PSTF 3-tensors ψ ab , these are decomposed as where n a Ψ a = 0 and Here we have used a notation with curly brackets to denote the PSTF-part of a tensor with respect to the projection tensor N ab . For the PSTF 3-tensors relevant here, we therefore write π ab = Π n a n b − In a similar manner as before, some new projected variables are also obtained when considering the spatial derivative of n a . Defining the spatial gradient along n a ,T ab... cd... ≡ n e D e T ab... cd... , and the spatial derivative fully projected onto the 2-sheets, the spatial derivative of n a can be written as where ab ≡ abc n c . Here φ ≡ δ c n c describes an expansion of the 2-sheets, ξ ab ≡ δ [a n b] describes a twist of these sheets, ζ ab ≡ δ {a n b} is the so called distortion, and a a ≡n a [19]. After having decomposed the relevant variables, we also need to decompose the equations from Sec. II. This involves decomposing the covariant derivatives of scalars, 2-vectors, PSTF 2-tensors, ab and N ab . After inserting all of the decompositions into the equations, they are then contracted with various combinations of u a ,n a , ab and N ab to extract the individual equations for the projections. However, since this is a rather long process, we will omit the details here and continue with the perturbative approach.

IV. BACKGROUND SPACETIMES
We begin by considering the equations to zeroth order. As such, we need to specify the properties of the background spacetimes to be studied. These are chosen to be homogeneous and hypersurface orthogonal members of LRS class II, which is characterized by that ξ, ω ab , and H ab are all identically zero [15,22]. To be able to use the same type of harmonic decomposition for all of the studied spacetimes, we will also impose the requirement that φ = 0 [12]. Additionally, it is also assumed that the fluid velocities V (i) , and the interaction terms F (i) and ε (i) all vanish identically to zeroth order. Due to all of these assumptions, the only relevant non-zero dynamical variables to zeroth order are Setting all other variables to zero in the equations from the previous section and using the spatial homogeneity of the background, the quantities in S (0) are found to satisfy the relationṡ where It should here be noted that we need an additional equation for the pressures p (i) to get a closed system. Hence, on specifying equations of state p (i) = p (i) (µ (i) , N (i) ) and a set of initial conditions, the behavior of the background spacetime is fully determined by the system above. Furthermore, on introducing coordinates, the line element for this class of spacetimes can be written as [12,22] where the scale factors a 1 and a 2 are related to the expansion rate and the shear through the relations As for the precise nature of the function f κ (ϑ), this will depend on the 2D curvature [11] R = 2κ For κ = 0 the 2-sheets are flat with either f 0 (ϑ) = 1 or f 0 (ϑ) = ϑ 2 . If instead κ = 1, the sheets are closed with f 1 (ϑ) = sin 2 ϑ. Lastly, if κ = −1 the sheets are open with f −1 (ϑ) = sinh 2 ϑ.

V. PERTURBATIONS
With the zeroth order spacetime specified, we now consider first order perturbations on this background. To avoid unphysical modes due to the implicit mapping between the background spacetime and the perturbed manifold, we will describe the perturbations in terms of gauge invariant variables [21]. Making use of the Stewart-Walker lemma, the gauge invariant variables to first order are chosen to be quantities that vanish on the background [9]. Hence, quantities such as H ab and B a , which are naturally zero on the background, are gauge invariant to first order. However, not all quantities are zero on the background. Thus, to represent the perturbations of the quantities in S (0) , we will follow [11][12][13][14] and make use of the spatial gradients of these quantities, as these gradients vanish on the background due to the assumed homogeneity. We therefore introduce the following gauge invariant variables Note that we do not introduce any new gauge invariant variables using the "ˆ" gradients here, as these can be given in terms of the "δ" gradients and the vorticity after performing the intended harmonic decompositions later on [13,14]. Our complete set of gauge invariant first order variables therefore becomes To obtain linearized equations governing the variables in S (1) , we write the decomposed equations from Sec. II in terms of these quantities by using commutation relations for the projected derivatives, which can be found in appendix A. In doing so, we also omit every term that is of second order or higher. This yields a large set of equations involving derivatives of the projected variables with respect to both time and space. Since the equations originating from the Ricci and Bianchi identities can be found in [14], we will here only state the equations related to the electromagnetic fields, the individual plasma components, and the decomposition of the total energy-momentum quantities in terms of their constituents.

VI. HARMONIC DECOMPOSITION
To simplify the linearized equations, they are then harmonically decomposed using harmonics P k (z) and Q k ⊥ (ϑ, ϕ), which are suited for homogeneous and hypersurface orthogonal LRS class II backgrounds with φ = 0 [12]. These harmonics are defined on the background as functions satisfying the equationŝ where the differential operators are defined as∆ ≡ n a ∇ a n b ∇ b , and δ 2 = δ a δ a , whilst k and k ⊥ are comoving dimensionless wave numbers. Using these harmonics, scalars can be decomposed as where the coefficients Ψ S now only depend on time. To decompose 2-vectors Ψ a and PSTF 2-tensors Ψ ab , we also introduce vector and tensor harmonics following Using these harmonics, we may write Here we have used bars to denote the so called odd parts of the decompositions, whilst the even parts are written without bars. In the following, when considering the relations for the individual harmonic coefficients, we will omit their subscripts, as this should not lead to any ambiguities. Inserting the harmonic decompositions into the linearized equations and using the properties of the harmonics, described in appendix B, the system can be reduced significantly by solving for some of the harmonic coefficients in terms of the others. The harmonic decompositions of the linearized equations, which were presented in Sec. V, can be found in appendix C. However, here we encounter a difficulty which was not present in previous papers. At this point in previous works, the system effectively divides into two separate parts, one containing the odd coefficients and another containing the even parts [12][13][14]. One exception to the division into even and odd equations, which is also present in previous works, are quantities defined using a Levi-Civita psuedo-tensor, such as the vorticity and the magnetic part of the Weyl tensor. For these quantities it is generally the even coefficients that are present in the odd system and vice versa, but the systems still decouple from each other. However, when including a non-zero magnetic field to zeroth order, the division into even and odd subsystems does not seem possible in general, as the pseudoscalar B serves as a coupling that connects the systems together. This can be seen when considering the harmonically decomposed momentum equations for the fluid components Here it can be explicitly observed that V V (i) is coupled to V V (i) through the terms proportional to B. To avoid this problem and simplify matters even further, we will now apply the cold MHD approximation

VII. COLD MHD APPROXIMATION
In the cold MHD approximation, we assume that µ (i) ≈ m (i) N (i) and p (i) ≈ 0 to both zeroth and first order, and write the harmonic equations in terms of the collective variables The equations governing these variables are found by adding the equations for the individual fluid components. Using charge neutrality to zeroth order and assuming that the plasma consists of electrons and a much more massive ion counterpart, it can then be shown that the collective velocities satisfy the equations where we have written the currents as . This simplifies matters, and the system is ultimately seen to decouple into two subsectors, as shown in Sec. VIII. However, we will need some model for the electrical resistivity if we want to completely close the system. Therefore, when performing numerical calculations in Sec. IX, we will make the simplified assumption that η follows a Spitzer-like expression so that where T (e) is the electron temperature [23]. Assuming that the electron pressure, although neglected in the main equations, follows a polytropic equation of state and an ideal gas law, it follows from the electron number conservation and Eq. (99) thaṫ VIII. FINAL SYSTEM The harmonically decomposed system of equations in the cold MHD approximation can be significantly reduced by solving for some of the coefficients in terms of the others. However, before arriving at the final system, we should note that there is still some freedom in specifying the dyad {u a , n a } to first order. To zeroth order, this dyad is fixed. The vector field u a is then assumed to be orthogonal to the hypersurfaces of homogeneity, and coincides with the plasma velocities. As for the field n a , this is fixed since it defines the direction of the local rotational symmetry. However, when going to first order, the homogeneity and rotational symmetry are broken, and the vector field u a does not need to coincide with the plasma velocities. As a result, there is some freedom in defining u a and n a to first order. We will use this freedom to set a V , a V , A S , A V , and A V to zero. For details on the effect and use of infinitesimal transformations of the dyad {u a , n a }, the reader is referred to [11]. Finally, after specifying the dyad to first order, we arrive at two subsystems for the first order quantities. These systems are seen to close independently of each other and consist of the odd sector and the even sector All other coefficients can be given algebraically in terms of these sets.

A. Odd Sector
The coefficients in the odd sector satisfy the following evolution equationṡ where the auxiliary variables are defined as We get the remaining relevant harmonic coefficients as algebraic expressions in terms of The remaining scalar coefficients are The remaining vector coefficients are in turn given by ik Finally, the algebraic tensor coefficients are

B. Even Sector
The evolution equations for the coefficients in the even sector arė where the additional auxiliary variables are defined as The remaining relavant coefficients can be written in the following way in terms of The scalar coefficients are ik The vector coefficients are in turn given by Finally, the only remaining tensorial coefficient is

IX. ANALYZING THE FINAL SYSTEM
Since the final equations from the previous section are still rather difficult to deal with analytically, we will mainly consider numerical examples. In doing so, we will focus on analyzing how the perturbations depend on the chosen length scales. However, before moving on to the numerical calculations, we should make some notes about the inherent wave length dependencies introduced when defining the gauge invariant variables and the harmonics.

A. Rescaling the Harmonics
When investigating the dependence on length scale, it is important to note that, due to our definition of the vector and tensor harmonics as gradients of the scalar harmonics, these will naturally contain factors of k ⊥ . Hence, the vector and tensor coefficients will contain inherent factors of k ⊥ relative to the scalar coefficients. Defining a new set of harmonics, where the inherent factors of k ⊥ have been removed, we note that where the Ψ * are the coefficients relative to the new harmonics. These harmonics have the benefit that the scalar, vector and tensor coefficients do not differ with factors of k ⊥ simply due to the definition of the harmonics. However, there are still some inherent factors of k ⊥ that need to be addressed. These factors occur for the gauge invariant variables that are defined as gradients of quantities that are non-zero on the background spacetime. If we let Φ a = δ a Φ, then the δ gradient implies that the coefficients { Φ V * , Φ V * } will naturally be one order higher in the factor k ⊥ /a 2 than the coefficients for the corresponding scalar perturbation of Φ. Hence, when investigating the length scale dependencies, we will focus on coefficients of the form V * } to avoid the aforementioned inherent factors of k ⊥ due to our definitions.

B. Dependence on Perturbation Length Scale
Equipped with the full system of ordinary differential equations and a rescaled set of harmonic coefficients, we now consider specific numerical examples. In these calculations, we will follow a similar convention as in [10] and set Λ = 1, which will fixate the remaining length dimension. As such, we will in the following treat all quantities as being dimensionless. In our calculations we will use the background spacetime specified by the following initial conditions and parameter values where r B is the fraction between the initial magnetic energy density and the total energy density. The initial conditions for µ, Σ, and Θ are identical to those for a dust background presented in [10], which ends at an expanding de Sitter solution. When not stated otherwise, we will use r B = 0.1 and the initial resistivity η(t 0 ) = 10 −3 . We will also choose the k-space direction k = 2k ⊥ = 2k(t 0 )/ √ 5, where When not stated otherwise, the equations are solved using ode45 in MATLAB with a relative tolerance of 10 −9 and an absolute tolerance of 10 −34 .

Tensor Perturbations
As a first example, we consider the generated magnetic fields due to perturbations in some tensorial coefficients. Following a similar line as in [4] we will let the initial tensorial shear perturbation represent some initial gravitational wave content. The aim is then to investigate how the generated magnetic fields depend on the chosen length scale. Defining the initial Hubble length L H = 3/Θ(t 0 ) and a characteristic length scale for the perturbations L = 2π/k(t 0 ), we solve the final system of equations for values of L/L H in the range 10 −3 , 10 3 . As for the initial perturbations, we choose or equivalently where k m is the maximum value of k(t 0 ) in the studied range. Requiring B V (t 0 ) = 0 and using Eqs. (129) and (149), we also let in particular, the perturbations of the electromagnetic fields and the plasma velocities are all zero at t 0 . Solving the system, we then get the results shown in Fig. 1a-c. For super-horizon scales, L L H , the curves in the logarithmic plot are linear with a slope close to two. Decreasing L below L H , the values extracted at the fixed time t = 8 begin to oscillate. In the region around L/L H = 10 −1 , the curves displaying the maximum values are almost linear with a slope of approximately one. Eventually when decreasing L even further, it can be seen that the curves approach a constant value. Finally, it should be remarked that the noise-like spikes at values of L/L H slightly larger than unity, as can be seen in Fig. 1c, are most likely due to numerical problems encountered when the auxiliary variable B, defined in Eq. (110), passes through zero for t > t 0 . As the final system of ordinary differential equations have been derived under the assumption that B is non-zero, a zero-crossing for B clearly poses a problem, as the derivatives involving factors of B −1 are not well defined when B = 0. In addition to the derivatives not being well defined, we also encounter problems with the coefficients that are given as algebraic expressions in terms of the coefficients in the final system, as their corresponding expressions may involve factors of B −1 . Therefore, we also have to be careful when specifying initial conditions for the perturbations when B(t 0 ) ∼ 0, since the coefficients with algebraic expressions, which should remain small relative to the background, may become too large even if the initial conditions for the variables in the ODE system appear small.

Velocity Perturbations With Vanishing Initial Ohmic Current
As a second example, we now consider how the generated magnetic field is affected by perturbations in the plasma velocity. To avoid aforementioned problems with the initial conditions when B(t 0 ) ∼ 0, we write the initial velocity perturbation as 3 Since the electrical resistivity is rather small, we also introduce the following initial perturbations in the electric field as to avoid large initial ohmic currents. Finally, requiring that B V (t 0 ) = 0, we use Eq. (149) and choose The other perturbations, that are explicitly present in the ODE system are set to zero initially. Performing the calculations over the same range of length scales as in the previous example, the results in Fig. 1d-f are obtained. There it can be seen that for superhorizon scales, the maximum and end values at t = 8 follow each other closely with a slope of approximately unity in the logarithmic scale. However, for subhorizon scales, the curves displaying the curves seemingly diverge. The maximum curve continues with a slope of roughly −1.7 whilst the curve representing the values at t = 8 follow a slope closer to −1. Finally, it should be noted that the values at t = 8 begin oscillating when decreasing the length scale below the Hubble scale. However, in contrast to the behavior in the odd sector shown in Fig. 1f, the oscillations of the quantities in the even sector, displayed in Fig. 1d-e, have become unnoticeable for length scales below 10 −2 L H .
Finally, as we are considering initial velocity perturbations, it is also interesting to investigate the generation of perturbations in the tensorial coefficients. In Fig. 2a-

Velocity Perturbations With Vanishing Initial Energy Flux
In the previous example, both velocity and electric field perturbations were introduced to ensure that the ohmic currents and magnetic fields vanished initially. However, in doing so we introduce an initial energy flux that will affect many gravitationally related variables. For instance, this flux appears in both the differential equations and algebraic relations pertaining to the electric and magnetic parts of the Weyl tensor. Hence, the previously chosen initial values may implicitly contain gravitational effects. Therefore, to better isolate the plasma related effects, we also investigate the case where the initial conditions are chosen so that the energy flux vanishes initially. For this purpose we introduce where we seemingly no longer need the factor B(t 0 ) in the velocity perturbations. It should be noted that these initial conditions for the electric field coefficients are incompatible with the ideal MHD limit, as they would in general imply an infinite initial ohmic current for η → 0. To avoid similar large initial ohmic currents in our calculations, we increase the resistivity to η(t 0 ) = 1. Finally, we set Then, on performing the calculations in the same range of length scales as in previous examples, we obtain the results shown in Figs. 1g-i. There it can be seen that, on superhorizon scales, the curves for the maximum value and the value at t = 8 decay with a slope close to −1. For subhorizon scales, the maximum curve and the general trend of the oscillating t = 8 curve both appear rather independent of the length scale. In Figs. 1h-i the values at t = 8 are notably smaller than the maximum values, which is also true for subhorizon scales in Fig. 1g. For superhorizon scales, the curves in Fig. 1g seemingly coincide.
With these initial conditions, it is also interesting to again examine the effect of the velocity perturbations on the tensorial quantities. In Fig. 2c-d we display the generated perturbations in the coefficients Σ T * and H T * . There it can be see that the generated tensorial perturbations increase with decreasing L/L H , with different slopes on super and subhorizon scales. On superhorizon scales, the curves for Σ T * have a slope of −3 while H T * have a slope of about −2 . On subhorizon scales, the curves for Σ T * and H T * both have a slope of −1. Hence, on attributing generated fields in Fig. 1g-i and Fig. 2 to plasma related effects, it can be seen that these effects become more important on subhorizon scales.

C. Beat Waves
To further investigate the plasma related effects on subhorizon scales, we perform similar calculations as in Sec. IX B 2, but varying the magnitude of the background magnetic field rather than the length scale of the perturbations. Hence, we introduce the initial conditions but set V S (F ) to zero. To highlight the sought effects, we decrease the initial resistivity drastically to η(t 0 ) = 10 −12 and choose the length scale to be L = 10 −5 L H . Since the smallness of the chosen electrical resistivity appears to make the system quite stiff, we use ode15s instead of ode45 in this section, although without specifying an analytical expression for the Jacobian. As for the initial magnitude of the background magnetic field, this is, for reasons made more clear in the next section, specified through the Alfvén velocity On using that in Eq. (163) and calculating for v A (t 0 ) = 0.95, the results in Fig. 3 are obtained for the tensorial coefficients of the magnetic part of the Weyl tensor. There it can be seen that a clear beat pattern emerges.
To get an analytical explanation of this beat pattern, we now seek wave equations for the plasma velocities and the magnetic part of the Weyl tensor in the limit of large wave numbers. Looking at Eq. (134), we see that V S (F ) will not give rise to any wave equation. Using Eq. (62) and (134), we find that the evolution of V S (F ) is simply determined by the scale factors on the backgrounḋ Thus, to obtain any interesting wave equations, we instead look at V V (F ) * and V V (F ) * . In doing so, we will assume that the plasma is ideal (η → 0), and that the wave lengths considered are much smaller than the scales on the background, so that We will then assume that the frequencies for the variables V T * are of order k in this limit, so that, when applying a time derivative on one of these coefficients, the result is one order higher in k. The idea is then to truncate the equations, keeping only the terms with the highest orders in k. However, to be able to compare the orders of various terms in the equations, we first have to consider the inherent difference in order between the harmonic coefficients.
With the chosen initial conditions in this section, Ω S , Ω V * and V S (F ) are identically zero, and will hence be omitted in the following. Motivated by numerical results using the aforementioned initial conditions, we will treat the magnitudes of H as being of the same order in k, up to factors of order 10. The coefficients µ V (F ) * and B V * are in turn treated as being one order higher in k whilst Σ T * is seen to be one order lower in k than V V (F ) * . With a similar numerical motivation, we assume that the magnitudes of E are all of the same order in k, whilst the magnitude of E S is significantly smaller. The smallness of E S is consistent with the ideal limit η → 0, as E S should then tend to zero.
On using these relative orders in k and keeping only the terms of highest order in k, the following wave equations can be derivedV where the previously defined Alfvén velocity v A appears. When deriving Eqs. (180) and (181), we used Eqs. (107) and (137) to write the currents in Eqs. (105) and (135) in terms of derivatives of the electric field. The electric field coefficients were then eliminated by noting that, in the ideal limit η → 0, we must have for the currents in Eqs. (97)-(98) and (D.8) to remain finite. The wave equations for the plasma velocities can be compared with non-relativistic results for ideal MHD waves in magnetized plasmas [23]. It can then be seen that Eq. (180) is consistent with a fast magnetosonic mode in the limit of vanishing speed of sound. The fact that V S (F ) shows no wave-like behavior is also consistent with a slow magnetosonic mode in this limit. This is reasonable since we have neglected the fluid pressures, and hence also the acoustic modes. As for V V (F ) , Eq. (181) is seen to be consistent with a shear Alfvén mode. Continuing with the tensor coefficients of the magnetic part of the Weyl tensor, these can, to highest order in k, be found to satisfy the equationsḦ Assuming that V ( Here k p = |k /a 1 |, and all scale factors and background quantities should be evaluated at t 0 = 0. It should here be noted that we have set the first order derivatives to zero at t 0 . This is consistent with the first order differential equations, as they imply that the derivatives of these harmonic coefficients are of order k 0 at t 0 . As the amplitudes of the derivatives should be of order k, their initial values of order k 0 are neglected to highest order. Studying the solutions for H T * and H T * we clearly see the cause of the beat pattern in Fig. 3, which emerges due to the interaction between gravitational wave modes and the magnetized MHD waves. It should also be noted that a resonance can occur when the Alfvén velocity approaches unity. For H T * , the resonance is independent of the precise relation between k and k ⊥ , whereas the factor k p /k limits the resonance for H T * . Comparing the analytical solution for H T * with the numerical data, we see a good agreement initially in Fig. 4a. However, in Fig. 4b it can be seen that, as time increases, the numerical results become shifted relative to the analytical expression, indicating a change of frequency for the numerical results. The amplitude of the numerical results also seems to decrease slightly in comparison to the analytical expression, but this effect is less noticeable than the frequency shift for the chosen parameters. These differences between the numerical and analytical results are expected, as we, when deriving the analytical expression, have neglected frequency redshifts and damping effects due to the evolution of the background spacetime.

X. DISCUSSION
Some of the numerical results presented here agree rather well with analytical expressions presented in [4], which were found using second order gauge invariant perturbation theory on spatially flat FLRW backgrounds in the ideal and cold MHD limit. In accordance with Eq. (50) in [4] it was there found that, in the superhorizon limit, the dominant late-time contribution to the generated magnetic field perturbation due to gravitational waves depended quadratically on the ratio between the length scale of the magnetized region to first order and the Hubble scale. This  contribution was also found to be proportional to the initial shear perturbation [4] . Similar quadratic dependencies on the characteristic length scale of the magnetic field have also been reported in [2,3]. Replacing the length scale of the magnetized region in [4] with the perturbation length scale L used here, we find, despite the differences in approach and background spacetimes, a similar quadratic dependence in Figs. 1a-c, where we have assumed that the initial tensorial shear perturbations are independent of the length scale. It should also be noted that the final result in [4] for the generated magnetic field perturbation, which was obtained by integrating the equations, has been argued to not be gauge invariant [7]. Since our generated gauge invariant magnetic field variables B V * , B V * , and a 2 k −1 ⊥ B V * show similar quadratic dependencies on the characteristic length scale, despite the differences, this result seems rather robust and independent of the precise details. However, focusing on the subhorizon scale, we note some differences between our results and Eq. (50) in [4]. Setting the initial velocity perturbation to zero in Eq. (50) from [4], it predicts that the dominant late time contribution in the subhorizon limit should depend linearly on the characteristic length scale of the magnetic field. When including velocity perturbations, the dominant contribution is instead predicted to be proportional to the initial velocity perturbation and inversely proportional to the length scale of the magnetic field [4]. However, with the initial conditions, integration times, and wave lengths considered here, the generated magnetic fields in Figs.1a-c and g-i are found to be rather independent of the length scale of the perturbations in the subhorizon limit. Although the values at t = 8 oscillate in this limit, the general trend of these values and the calculated maximum amplitudes seem to be independent of L/L H . When considering Figs.1d-f, it must be noted that B(t 0 ) introduced in the initial conditions for this case contributes with a factor L −2 through the initial velocity perturbation in the subhorizon limit. Therefore, if we assume that the generated field perturbation is proportional to the initial velocity perturbation and remove the additional factors of L due to B(t 0 ) by multiplying with L 2 , we would expect the maximum values to be almost independent of the length scale in this limit whilst the values at t = 8 show a slope of unity. These subhorizon results are in contrast to the predictions from [4], from which we would expect Figs.1a-c to show a slope of unity and Figs.1g-i a slope of negative unity. However, looking at superhorizon scales in Figs. 1g-i we see the inverse dependence on length scale that [4] predicts for subhorizon scales. This may be due to that we, in some sense, have removed some gravitational effects from the initial conditions for Figs. 1g-i, so that the contributions from these effects no longer overshadow the velocity contribution in the superhorizon limit.
The exact explanation behind the difference between our subhorizon results and those in [4] is however hard to pinpoint without detailed analytical solutions. Possible, but speculative, explanations could be differences in the studied backgrounds, the initial conditions for the perturbations, the chosen model for the electrical resistivity, or that the wave lengths studied here were simply too large. It may also be due to some decaying non-dominant terms, similar to the implicit terms in [4], that have not yet completely vanished at the integration times that our results were extracted. Due to all the possible differences, it is not surprising that our results differ from those in [4]. It is, however, interesting that they qualitatively still agree rather well in the superhorizon limit.
Finally, it should be noted that the interactions we observe are crucially dependent on the background magnetic field. As an example, it can be seen in Eq. (188) that the magnitude of the generated perturbations in the magnetic part of the Weyl tensor increases with increasing B(t 0 ), especially as v A → 1. Furthermore, looking at the equations in Sec. VIII, we note that if we were to set B = 0, this would imply that the equations for the electromagnetic fields decouple from the others, at least to first order. Thus, the interactions we see here to first order are fundamentally dependent on the anisotropic nature of the background spacetime, as this is what allows us to have a non-zero background magnetic field.

XI. CONCLUSIONS
In this paper, we have studied linear electromagnetic, gravitational, and plasma related perturbations on homogeneous and orthogonal LRS class II spacetimes, where we have allowed for a non-zero magnetic field on the background. Using a 1 + 1 + 2 covariant formalism and ultimately applying the cold MHD limit, we have derived a closed set of equations governing the harmonic coefficients of the perturbations. On analyzing the system, we have observed the generation of magnetic field perturbations due to initial perturbations in tensorial quantities related to gravitational waves. The generated fields were then found to depend quadratically on the characteristic length scale in the superhorizon limit, in agreement with previous works using FLRW backgrounds. However, in the subhorizon limit, the similarities were not as apparent.
In addition to the generation of magnetic field perturbations, we also observed beat waves arising due to interference between gravitational waves and the magnetized MHD modes. Thus, replacing the background FLRW spacetime with an anisotropic LRS class II spacetime and allowing for a zeroth order magnetic field, we were indeed able to see interesting interactions already to first order in the perturbations.

A. COMMUTATION RELATIONS
Here we state the commutation relations between the projected derivatives to first order. These relations are crucial when writing the first order equations in terms of the gauge invariant variables. For a general scalar Ψ that is non-zero to zeroth order, we haveΨ For 2-vectors Ψ a we similarly haveΨā Finally, for general PSTF 2-tensors Ψ ab , it holds thaṫ

B. PROPERTIES OF THE HARMONICS
We here state some relations for the harmonic functions that are useful when performing the harmonic decomposition of the linearized equations. The even and odd vector harmonics satisfy the relations whereas for the tensor harmonics, we have (B.13)

C. HARMONIC COEFFICIENTS IN THE MULTIFLUID SYSTEM
Here we state the harmonic decompositions of the linearized equations. Since the equations originating from the Ricci and Bianchi identities can be found in [14], we here only show the additional equations pertaining to the decomposition of the total energy-momentum tensor, energy-momentum and particle conservation for the individual plasma components, and Maxwell's equations. (C.14) (C.16) for a scalar Ψ that is non-zero to zeroth order, it should also be noted that it follows directly that ik a 1 Ψ V = 2Ω VΨ . (C.36)

D. HARMONIC COEFFICIENTS IN THE COLD MHD SYSTEM
Here we state the equations for the harmonic coefficients after applying the cold MHD approximation. Only equations that differ from the multifluid case are presented.
The differential equations pertaining to the collective energy densities and velocities arė where the currents in the momentum equations have been written as (D.9) , (D.10) using Ohm's law. As for the harmonic coefficients related to the decomposition of the total energy-momentum tensor, these can in turn be written as (D.14) It should be noted that the Maxwell equation involving the divergence of the electric field as well as the equations for the particle density gradients are no longer needed to close the system. The former can instead be used to determine the charge density from some specified vorticity and electric fields.