On the locally rotationally symmetric Einstein-Maxwell perfect fluid
OOn the locally rotationally symmetric Einstein-Maxwell perfect fluid
Daniela Pugliese , ∗ and Juan A. Valiente Kroon † School of Mathematical Sciences,Queen Mary University of LondonMile End Road, London E1 4NS, UK Institute of Physics, Faculty of Philosophy & Science,Silesian University in Opava, Bezruˇcovo n´amˇest´ı 13,CZ-74601 Opava, Czech Republic (Dated: September 28, 2018)We examine the stability of an Einstein-Maxwell perfect fluid configuration with a privilegeddirection of symmetry by means of a 1 + 1 + 2-tetrad formalism. We use this formalism to cast,in a quasi linear symmetric hyperbolic form the equations describing the evolution of the system.This hyperbolic reduction is used to discuss the stability of solutions of the linear perturbation.By restricting the analysis to isotropic fluid configurations, we made use of a constant electricalconductivity coefficient for the fluid (plasma), and the nonlinear stability for the case of an infinitelyconducting plasma is also considered. As a result of this analysis we provide a complete classificationand characterization of various stable and unstable configurations. We found in particular that inmany cases the stability conditions is strongly determined by the constitutive equations by means ofthe square of the velocity of sound and the electric conductivity, and a threshold for the emergenceof the instability appears in both contracting and expanding systems.
PACS numbers:Keywords: locally rotationally symmetric solutions, 1 + 1 + 2-formalism, perturbation, stability problem,Magnetohydrodynamics
I. INTRODUCTION
The stability problem of plasma configurations is an important issue in a variety of astrophysical scenarios involvingfor example stellar objects and accretion disks, and various phenomena of the high energy Astrophysics, related tothe accretion disks with the instability processes as the accretion or the Jet emission. In this article we consider thesituation where gravity plays a decisive role in determining both the equilibrium states of the configurations thenthe dynamical phases associated to the instability, requiring a full general relativistic analysis [1–9]. Consequentlythese systems are described by the coupled Einstein-Maxwell-Euler equations. Notable examples are the generalrelativistic (GR) magnetohydrodynamic (MHD) systems. Often very complicated, the perturbation analysis must beconducted with suitable assumptions on the symmetries for the system (for example in the toroidal accretion disks)and the dynamics, numerical approaches are often required. A major challenge in dealing with these systems is tofind an appropriate formulation of the problem: from one side to set up a formulation adapted to the configurationsymmetries and, on the other side a suitable formulation of an initial value problem is necessary for the constructionof the numerical solutions, in order to ensure the local and global existence problems and the analysis of the stabilityof certain reference solutions[10] for a general discussion see for example[6–9, 12–16].In this work we set the Einstein-Maxwell-Euler equations in a quasi-symmetric hyperbolic form, we explore thestability properties of a perfect fluid configuration with a preferred direction of symmetry coupled with the electro-magnetic field. The formalism is adapted to the description of a general locally rotationally symmetric system, aremarkable example is the simple case of a spherically symmetric configurations. For a more specific discussion on thestability of spherically symmetric plasma see for example [17–22]. The special case of an infinitely conducting plasmadescribes an adiabatic flow so that the entropy per particle is conserved along the flow lines.The plasma configuration instability, especially in the geometrically thin toroidal structures orbiting around anattractor (for example in the Shakura-Sunyaev accretion disks) is often described by the magneto-rotational instability(MRI). The dissipative (visco-resistive) effects are essential in these models as they allow the transport of angularmomentum in the configurations in accretion on the central object. In fact, in the geometrically thin configurations itis assumed that the time scales of the dynamical process (balance in pressure and gravitational and centrifugal forces) ∗ Electronic address: [email protected] † Electronic address: [email protected] a r X i v : . [ a s t r o - ph . H E ] J un are less then the thermal one (for dissipative heating and radiation) that is less then the viscous ones (dissipativestresses and consequent angular momentum transport). The magnetic field, the dissipative effects and the radialgradient of the plasma relativistic angular velocity are therefore essential for the MRI instability. However, someaspects of the theoretical framework of the MRI and accretion process are still to be clarified. An intriguing issuefor example is the so-called visco-resistive puzzle: eventually high values of and resistivity and viscosity have to beassumed. In this work we investigate the stability problem for the systems with the magnetic field but not dissipativeeffects, showing the emergence of the instability and a threshold for the occurrence of the unstable modes which isessentially regulated by the conductivity parameter σ J and the speed of sound v s , even for configurations with morespecific symmetries, considered here in classes and subclasses of solutions, stable and unstable for linear perturbation.Considering both contracting and expanding systems (according to the sign of the kinematic expansion scalar Θ), weshow a threshold for the instability of the system determining two ranges for the density of matter and especially forthe shear scalar Σ along the privileged direction of the system. More specifically we consider a barotropic equation ofstate: when the fluid entropy is a constant of both space and time, an equation of state to link the pressure p to thematter density ρ can be given in the form p = p ( ρ ) —see e.g. [10]. In the present work we restrict our attention toisotropic fluids. We consider a one species particle fluid ( simple fluid ) and, since no particle annihilation or creationprocesses is expected, we use the equation of conservation of particle number. Moreover, we assume a polytropicequation of state with a constant velocity of sound and we specify the form of the conduction current using the Ohm’slaw , so that a linear relation between the conduction current and the electric field holds. By restricting our attentionto the isotropic fluids configurations we can make use of a constant electrical conductivity coefficient for the fluid(plasma). These assumptions simplify considerably the analysis of the stability problem of linear perturbations.A central aspect of the stability analysis to be pursued in this article is the construction of a quasi linear symmetrichyperbolic evolution system for the variables of the configuration. This, in turn, ensures the well-posed Cauchyproblem for the system —in other words, the local existence and uniqueness result for the Einstein-Maxwell-Eulerequations. By prescribing suitable initial data on an initial hypersurface, a unique solution exists in a neighbourhoodof that hypersurface. This solution depends continuously on the initial data. Accordingly, we first write the evolutionequations for the independent components of n variables collected in n -dimensional vector v used to obtain a suitablesymmetric hyperbolic evolution system of the form A t ∂ t v − A j ∂ j v = Bv , (1)where A t and A j and B are smooth matrix valued functions of the coordinates ( t, x ) and the variables v with theindex j associated to some spatial coordinates x . The system is symmetric hyperbolic if the matrices A t and A j aresymmetric and if A t is a positive-definite matrix. The evolution equations is complemented by constraint and theconstitutive equations. As the purpose of this work this is the study of the properties of the evolution system and theanalysis of its linear stability, the problem of the propagation of the constraints will only be briefly discussed referringfurther details to the literature —see in particular [23, 24].Our analysis is based on an adapted 1 + 1 + 2 -tetrad formalism for the locally rotationally symmetric spacetimes (LRS)[25], the simplest example being the spherical symmetric configurations —that is, a covariant decomposition ofEinstein-Maxwell perfect fluid field equations which is particularly suitable for LRS systems. This formalism is anextension of the usual 1 + 3-formalism in which the existence of a privileged timelike vector field u a assumed —inapplications involving the description of fluids, it is natural to let u a to follow the congruence generated by the fluid[26]. In the 1 + 1 + 2-formalism the presence of a further (spatial) vector field n a is assumed. This gives rise to a asecond split of the 1 + 3-reduced equations on the plane parallel and orthogonal to n a . This type of decomposition isparticularly useful in the presence of symmetries as one can naturally fix the spatial vector on the privileged symmetrydirection (from now on, for simplicity “radial direction”) at each point of locally rotationally symmetric spacetimes[25]. In the case of spherically symmetric spacetimes for example (here indicated with LSS), it is natural to choose n a to point in the radial direction of the spherical symmetry. After the decomposition, all tensors are covariantly splitinto scalars, vectors, and transverse-traceless 2-tensors, with respect to n a . In the case considered here, as a result ofthis split, it can be shown that it is only necessary to consider radially projected tensors. Further discussion on thedecomposition of Einstein-Maxwell-perfect fluid equations can be found in [25, 27, 29]. In [20] the same formalismhas been used to analyse self gravitating spherically symmetric charged perfect fluid configurations in hydrostaticequilibrium. Details on the 1 + 3 and 1 + 1 + 2-decompositions of the Einstein-Maxwell equations in LRS spacetimescan be found in [30–32]. We consider perturbations of the metric tensor, of the matter and EM fields, all the quantitiesshare the same preferred direction of symmetry, as in locally rotationally symmetric (LRS) classes II space-times asdescribed in [25] [49].Taking into account the gravito-electromagnetic (GEM) effects, we provide a complete classification of the solutionsin terms of the scalars of the Weyl’s conformal tensor. A general discussion and classification, of the solutions of theEinstein-Maxwell equations, through the scalars of the Weyl tensor can be found in[39].Here the consider the electromagnetic fields considered through real vector functions, while in many of the LSRspacetimes in [27, 29, 30, 33, 34] was naturally used the complex variable ψ = E + iB and ψ ∗ = E − iB in order todecouple the equations with the appropriate symmetries and obtain linear equations in the fields. As pointed out in[27] there are two ways to proceed depending of how one considers the coupled fields and test fields: specifically onhow one views the Maxwell’s and Einstein’s equations, and the gravitational effects on the EM field, in other wordsthe analysis of test fields on a background. The gravitational background must, in any case, have the same symmetriesof LRS systems: a suitable example is therefore the spherically symmetric schwarzschild solution. Alternatively, onecan consider other scenarios, we refer to a general discussion of the analysis in [27, 29, 30, 33, 34]. Here we simplyobserve that such a system can be self-gravitating or not; one can consider the perturbations of the background orotherwise to fix the spacetime.In this article we address the more general case, considering also the perturbations of the gravitational part, interms of scalars of the Weyl tensor. We assumed the symmetries to be preserved by the perturbation and limitingthe analysis to constant velocity of sound and conductivity parameter. To simplify the discussion of the results, thecalculation procedures, and the complete systems before the assumption of one privileged symmetry of direction, arespecified in some Appendix Sections.The procedure of hyperbolic reduction used in the present article follows the presentation given in [40]. Thisparticular analysis is independent of geometric gauge considerations[50]. In the present work we also provide a suitablepropagation equation for the fluid radial acceleration. To this end, we introduce an auxiliary field corresponding tothe derivative of the matter density projected along the radial direction. suitable field and evolution equations canbe obtained for this quantity. The resulting evolution system is then used to analyse the stability problem for smallnonlinear perturbations of a background solution. More precisely, we perform a first order perturbation to v of theform v (cid:55)→ (cid:15) ˚ v + ˘ v , where the parameter (cid:15) sets the order of the perturbation while ˘ v describes the (linear) perturbationof the background solution. Now, assuming the background variables ˚ v to satisfy the evolutions we end up with anevolution system for the perturbations of the form˚ A t ∂ t ˘ v − ˚ A j ∂ j ˘ v = ˚ B ˘ v . The core of the stability analysis consists of the study of the term ˚ B using some relaxed stability eigenvalue conditions.The procedure to analyse stability used here is adapted from[24] —see also discussion in [42] and cited references.Under suitable circumstances it can be regarded as a first step toward de analysis of non-linear stability. In our case,the elements of matrix ˚ B are, in general, functions of the space and time coordinates. For a general discussion onthe time dependent case and the case of non constant matrix coefficient (depending on both time and space) we referto [24]. The case of a linearised system where the coefficients are constant matrices is discussed in [45]. Finally, thecase of systems with vanishing eigenvalues has been discussed in [46–48]. In our case, a fully analysis of the stabilityproperties of the system turns out extremely cumbersome because of the form of ˚ B associated to the present problem.We will proceed with the analysis of the values of the eigenvalues using indirect methods based on the inspection ofthe characteristic polynomial. In order to keep the problem manageable, we analyse a number of simplified systemsobtained by making some assumptions about the configuration. More precisely, we consider background configurationswith a vanishing radial acceleration case for the reference solution, and we explore particular models with fixed valuesof the kinematic scalars. This analysis constitutes the main result of the article. We provide a detailed classificationby considering systems with particular kinematic configurations defined by fixing the radial acceleration of the fourvelocity, and the 1 + 1 + 2-projected expansion, shear, twisting and the vorticity of the system: the stability conditionscan be strongly determined by the constitutive equations by means of the square of the velocity of sound and theelectric conductivity.Since a significant part of this work was dedicated to the formalization of the problem in symmetric hyperbolicform, a first part of the article was necessarily devoted to the presentation of the formalism and the explanation ofthe adopted notation. To simplify the discussion of the problem and the illustration of the results the article has beendeveloped into three parts: in the first ( I ) part, from Sec. (II) to Sec. (V), we introduced the 1 + 1 + 2 formalismdecomposing the equations and set system in symmetric hyperbolic form, which is the first outcome of this work.The second ( II ) part, Sec. (VI) and Sec. (VII), develops the perturbations and the system stability is analysed. Part II constitutes the main part of this article with a major discussion of the main results on the system stability in thefundamental classes of the solutions. A more specific discussion on the other subclasses can be found in the finalpart of this article. The third ( III ) part is constituted by the Appendix Sections, deepening details of the I and II parts, and explaining some important aspects of the decomposition procedure. We show the perturbed equations insymmetric hyperbolic form for the general case of non-zero radial acceleration discussing in detail the conditions forthe unstable configurations belonging to the various classes and subclasses of solutions.In details, the present article is structured as follows: the 1 + 3-formalism is briefly reviewed in Section II A. The1 + 1 + 2-decomposition that will be used in our present analysis is discussed in Section II B. In Section III we writethe 1 + 1 + 2-equations for the LRS system. Section IV provides a discussion on the thermodynamical quantities of thesystem. A summary of the evolution equations is given in Section V. The re-parametrised set of evolution equationsconsidered for the stability analysis is given in Section V A. Section VI discusses the perturbation to the first orderof the variables. Section VI A provides some general remarks on the set of perturbed equations and system stability.Section VII contains the main results concerning the nonlinear stability of the symmetric hyperbolic system. Someconcluding remarks are given in Section VIII. Finally, in Appendix A we provide an alternative symmetric hyperbolicsystem for the fluid variables. The 1 + 1 + 2-decomposition of these equations is in Appendix B. Some general noteson the evolution equations and hyperbolicity considerations can be found in Appendix C. II. PRELIMINARIES
We consider the stability problem for a configuration with a spatial direction of symmetry (privileged direction of thesystem admitting a one-dimensional isotropy group), described by the Einstein-Maxwell-Euler equations for a perfectfluid. By applying the 1 + 1 + 2 decomposition i we can take full advantage of the symmetries of the LRS system, andmoreover this procedure allows to construct in an easy, and relatively immediate manner, a quasi-linear hyperbolicsystem. This ensures also the consistency of possible numerical approaches to the problems, without requiring theintroduction of any auxiliary variable to handle both the propagation equations along the timelike direction then theconstraint part of the system. Morover, the covariant and gauge-invariant perturbation formalism turn to be especiallysuitable for dealing with spacetimes with some preferred spatial direction, not necessary spherical symmetry in thebackground, and to the application possibly to the case of gravitational wave propagation by introducing a radial unitvector, decomposing all covariant quantities with respect to this [13, 25, 27, 28]. All the equations and the quantitiesrelated to the fields and the curvature tensor will be decomposed according to the 1 + 1 + 2 procedure. For thispurpose it will be necessary to first fix the notation induced by the 1 + 3-decomposition as introduced Sec. (II A). Theprocedure of 1 + 1 + 2-decomposition, adapted to the symmetries of the system, will be discussed with some detailsin Sec. (II B).
A. The -formalism
The implementation of the 1 + 3-formalism used in the present article follows, as much as possible, the notation andconventions of [26]. We consider 4-dimensional metrics g ab with signature ( − , + , + , +). The Latin indices a, b, c... willdenote spacetime tensorial indices taking the values (0 , , ,
3) while i, j, k... will correspond to spatial frame indicesranging over (1 , , g ab will be denoted by ∇ a . Whenever convenient, weuse the semicolon notation. As usual, one has that ∇ a g bc = g bc ; a = 0.In what follows, the timelike vector field ( flow vector ) u a will describe the normalised future directed 4-velocity ofthe fluid. It satisfies u a u a = −
1. Indices are raised and lowered with g ab . The tensor h ab ≡ g ab + u a u b is the projectoronto the three dimensional subspace orthogonal to u a , thus, one has that h ab = δ ab + u a u b , h ba h ac = h bc , h ab h ba = 3 , h ab u a = 0 . Following the standard approach of 1 + 3-formalisms, we split the first covariant derivative of u b as ∇ a u b = σ ab + ω ab + 13 Θ h ab − ˙ u b u a , (2)where σ ab ≡ D (cid:104) b u a (cid:105) with σ ab u b = 0, ω ab ≡ D [ a u b ] with ω ab u b = 0, Θ ≡ D a u a and ˙ u a = u b ∇ b u a are, respectively, the shear and the vorticity tensors, the volume expansion scalar , and the vector. Moreover, we introducethe vorticity vector ω a ≡ (cid:15) abc ω bc / (cid:15) abc = u d (cid:15) dabc , (cid:15) abe (cid:15) abe = 6 and (cid:15) dabc stands for the totally antisymmetrictensor with (cid:15) = √− det g ab . Then, σ ab u a = 0 = ω ab u a = ˙ u a u a by construction. In the above expressions, theoperator D a corresponds to the 3-dimensional covariant derivative obtained from projecting the spacetime covariantderivative in the distribution orthogonal to u b . As an example of this procedure, for a generic 2-rank tensor T bc , onehas that D a T bc = h sa h tb h pc T st ; p , and T st ; p = ∇ p T st while ˙ T st ≡ u a ∇ a T st . For clarity, whenever necessary, projectedindices of a tensor will be highlighted by (cid:104)(cid:105) - brackets. For example, we write T (cid:104) ab (cid:105) = h ac h bd T cd . B. The -formalism
The implementation of a 1 + 1 + 2-formalism to be used in the present article follows the notation and conventionsof [25]. In what follows, let n a denote a spacelike normalised vector n a chosen along a preferred direction of thespacetime and which define a further projector tensor N ba n a n a = 1 , u a n a = 0 , N ba ≡ h ba − n a n b = δ ba + u a u b − n a n b , N ba u a = N ba n a = 0 , (3)on the space orthogonal to both n a and u a —a 2-dimensional sheet . In addition to this projector, we also introducethe following rank 2 totally antisymmetric tensor (cid:15) ab ≡ (cid:15) abc n c = (cid:15) dabc u d n c , (cid:15) ( ab ) = (cid:15) ab n b = 0 , (cid:15) abc = n a (cid:15) bc − n b (cid:15) ac + n c (cid:15) ab , (4) (cid:15) ab (cid:15) cd = N ca N db − N da N cb , (cid:15) ca (cid:15) bc = N ab , (cid:15) ab (cid:15) ab = 2 . (5)Given 3-vector ψ a , one can use the projector N ab to split it as ψ a = Ψ n a + Ψ ¯ a , with ψ ¯ a ≡ ψ b N ab , lying in the sheetorthogonal to n a . In order to avoid any confusion, we make use of an overbar on index to denote the projectionwith N ab . By means of the projector N ba one can decompose the kinematical quantities ˙ u a and ω a in the form˙ u a = A n a + A a , ω a = Ω n a + Ω a , with A ≡ n a u b ∇ b u a and Ω describing the components of the acceleration andthe vorticity in the direction of n a . Now, any spatially projected, symmetric, trace-free (PSTF) tensor ψ ab can besplit into a scalar Ψ = ψ b n b , a 2-tensor Ψ ab and a vector Ψ a ψ ab = ψ (cid:104) ab (cid:105) = Ψ (cid:18) n a n b − N ab (cid:19) + 2Ψ ( a n b ) + Ψ ab , (6)where Ψ ≡ ψ ab n a n b = − N ab ψ ab , Ψ a ≡ N ba n c ψ bc = ψ ¯ a , Ψ ab ≡ ψ { ab } ≡ (cid:18) N c ( a N db ) − N ab N cd (cid:19) φ cd . (7)In the above expressions, the curly brackets indicate the transverse-traceless part of the corresponding tensor. Notethat in the present context,the term “transverse” only refers to the fact that the tensor is orthogonal to n a . In thefollowing, where useful, we will denote with ( ⊥ ) the component of any tensor orthogonal to n a and with ( (cid:107) ) thecomponent parallel to it. It is convenient to further define the following two derivatives:ˆ φ c ··· da ··· b ≡ n e D e φ c ··· da ··· b , δ e φ c ··· da ··· b ≡ N ej N af · · · N bg N hc · · · N id D j φ h ··· if ··· g . (8)Accordingly, one can write D a n b = n a a b + 12 Φ N ab + ξ(cid:15) ab + ζ ab , a a ≡ n c D c n a = ˆ n a , Φ ≡ δ a n a , ξ ≡ (cid:15) ab δ a n b , ζ ab ≡ δ { a n b } . (9)The scalar Φ denotes is the expansion of the 2-sheet generated by n a , ζ ab the shear of n a and ξ is the rotation of n a (i.e. the twisting of the 2 - sheet), finally a a is the acceleration. Finally, also define the derivative˙ n a = A u a + α a , α a ≡ ˙ n ¯ a . (10)It can then be verified that˙ N ab = 2 u ( a ˙ u b ) − n ( a ˙ n b ) = 2 u ( a A b ) − n ( a α b ) , ˆ N ab = − n ( a a b ) , δ c N ab = 0 , (11)and, furthermore, that ˙ (cid:15) ab = − u [ a (cid:15) b ] c A c + 2 n [ a (cid:15) b ] c α c , ˆ (cid:15) ab = 2 n [ a (cid:15) b ] c a c , δ c (cid:15) ab = 0 . (12)On the decomposition of the kinematic quantities Applying the decomposition (6) to the symmetric trace free tensor σ ab appearing in the 1 + 3-decomposition given by equation (2) one obtains σ ab = Σ (cid:18) n a n b − N ab (cid:19) + 2Σ ( a n b ) + Σ ab , (13)where we have introduced the scalar Σ (the totally projected part of the shear of the 3-sheet ), the vector Σ a and the 2-tensor Σ ab and defined as in equations (7).The full covariant derivative of n a and u a can now be written as ∇ a n b = −A u a u b − u a α b + (Σ + Θ) n a u b + (Σ a − (cid:15) ac Ω c ) u b + n a a b + 12 φN ab + ξ(cid:15) ab + ζ ab (14) ∇ a u b = − u a ( A n b + A b ) + n a n b ( Θ + Σ) + n a (Σ b + (cid:15) bc Ω c ) + (Σ a − (cid:15) ac Ω c ) n b + N ab ( Θ − Σ) + Ω (cid:15) ab + Σ ab . The decomposition of the Weyl tensor
Let C abcd denote the Weyl curvature tensor of the metric g ab . As it iswell known, it can be 1 + 3-decomposed in terms of symmetric, traceless 2-rank tensors defined by E ab ≡ C abcd u c u d = E ( n a n b − N ab ) + 2 E ( a n b ) + E ab , (15) B ab ≡ (cid:15) ade C debc u c = B ( n a n b − N ab ) + 2 B ( a n b ) + B ab , (16)—the so-called electric E ab and magnetic B ab parts with respect to u a . Using equation (6), these rank 2 tensorsare 1 + 1 + 2-decomposed where the scalars E and B denote, respectively, the totally projected electric andmagnetic parts of the Weyl tensor. Decomposition of the matter fields
The energy momentum tensor T ab describing the matter-field content of thesystem under consideration is of the form T ab = T f ab + T em ab , T f ab = ρu a u b + ph ab , (17)where T f ab denotes the energy-momentum tensor of the perfect fluid given by the well-known expression with ρ and p denoting, respectively, the total energy density and pressure as measured by an observer moving with thefluid. Accordingly, the timelike vector field ( flow vector ) u a denotes the normalised future directed 4-velocityof the fluid. The electromagnetic energy momentum tensor T em ab is given by T em ab = (cid:18) F ac F cb − F cd F cd g ab (cid:19) , F ab = − (2 E [ a u b ] − (cid:15) abc B c ) . (18)where F ab is the electromagnetic field (Faraday) tensor , the Faraday tensor has been split in its electric part , E a ≡ F ab u b , and its magnetic part , B a ≡ (cid:15) abc F cd , with respect to the flow. Alternatively, one can rewriteequation (18) in the form T em ab = u a u b ( E + B ) + h ab ( E + B ) + P ab + 2 G ( a u b ) , (19)where we have written E ≡ E a E a and B ≡ B a B a , and P ab denotes the symmetric, trace-free tensor given by P ab = P ( ab ) ≡ h ab ( E + B ) − ( E a E b + B a B b ) . Furthermore, G a ≡ (cid:15) auv E u B v , denotes the Poynting vector . Thus we obtain T em ab n a n b = 12 ( E + B ) − ( (cid:107) E + (cid:107) B ) , E a = (cid:107) En a + ⊥ E a , B a = (cid:107) Bn a + ⊥ B a (20)the electric E a and magnetic B a are split relative to the vector n a . Now, recall that he Maxwell equations canbe written as ∇ b F ab = J a , ∇ [ a F bc ] = 0 where J a = (cid:37) C u a + (3) j a , (3) j a = (cid:107) jn a + ⊥ j a , (21)where J a is the 4-current. where (cid:37) C is the charge density , and the 3-current is (3) j a . III. -EQUATIONS FOR LRS SYSTEM
Evolution equations adapted to the 1 + 1 + 2-decomposition of spacetime discussed in the previous Sections canbe readily obtained from the 1 + 3-evolution equations using the split of equation (6) for a symmetric, trace-freetensor and projecting the original propagation equations along the longitudinal and orthogonal directions given bythe vector n a . For a generic symmetric, trace-free tensor one obtains three equations for the projected componentsΨ, Ψ a and Ψ ab . In the case of a LRS and LSS systems it is only necessary to consider the evolution equation of theradially projected component —i.e. the scalar Φ. This approach is allowed in any spacetimes with a preferred spatialdirection at each point, the so-called locally rotationally symmetric spacetimes (LRS) —see e.g. [25].In the case LSS spacetimes, n a is a vector pointing along the axis of symmetry and can therefore be thought of asbeing a radial vector. As there is no preferred directions in the 2-surface sheet all the non-radial components of thevarious tensors are assumed to vanish. Consequently, the kinematical quantities can be decomposed as˙ u a = A n a , ω a = Ω n a , σ ab = Σ (cid:0) n a n b − N ab (cid:1) . so that n ab ≡ ∇ a n b = −A u a u b + (cid:0) Σ + Θ (cid:1) n a u b + Φ N ab + ξ(cid:15) ab , (22) u ab ≡ ∇ a u b = −A u a n b + n a n b (cid:0) Σ + Θ (cid:1) + N ab (cid:0) Θ −
12 Σ (cid:1) + Ω (cid:15) ab , (23)while the decomposition of the Weyl tensor reduces to E ab = (cid:107) E (cid:0) n a n b − N ab (cid:1) , B ab = (cid:107) B (cid:0) n a n b − N ab (cid:1) . In ananalogous manner, it can be seen that the only non-vanishing components for the matter fields are given by the scalars( ρ, p, (cid:107) E, (cid:107) B, (cid:37) C , (cid:107) j ). Evolution equation for the Maxwell fields
Starting from the usual 1+3-decomposition of the Maxwell equations—see for example [26]— and applying the general procedure described at the beginning of the previous Sectionone readily obtains the following 1 + 1 + 2-evolution equations (cid:107) ˙ E = 2 ξ (cid:107) B − (cid:0) Θ − Σ (cid:1) (cid:107) E + ⊥ E a ( α a + Σ a ) + (cid:15) ab δ a ⊥ B b + (cid:15) ab (cid:0) A a ⊥ B b + Ω a ⊥ E b (cid:1) − (cid:107) j, (24) (cid:107) ˙ B = − ξ (cid:107) E − (cid:0) Θ − Σ (cid:1) (cid:107) B − (cid:15) ab δ a ⊥ E b + ⊥ B a ( α a + Σ a ) − (cid:15) ab ( A a ⊥ E b − Ω a ⊥ B b ) . (25)The above equations are supplemented by a pair of constraint equations for for the scalars (cid:107) E and (cid:107) B —seebelow. If one now assumes LRS, equations (24) and (25) reduce to (cid:107) ˙ E = 2 ξ (cid:107) B − (cid:18)
23 Θ − Σ (cid:19) (cid:107) E − (cid:107) j, (cid:107) ˙ B = − ξ (cid:107) E − (cid:18)
23 Θ − Σ (cid:19) (cid:107) B, (26)while the constraint equations are given by (cid:107) ˆ E + (cid:107) E Φ − (cid:107) B − (cid:37) C = 0 , (cid:107) ˆ B + (cid:107) B Φ + 2Ω (cid:107) E = 0 . (27)As in this LRS one has that E a = En a , B a = Bn a and J a = jn a + (cid:37) C u a , the superscript (cid:107) can be dropped fromthe above equations without giving rise to any ambiguity. Perfect fluid equations
From the conservation of the energy-momentum tensor ∇ a T ab = 0, it readily follows that u b u a ∇ a ( p + ρ ) + ( p + ρ ) (cid:0) u b ( ∇ a u a ) + u a ∇ a u b (cid:1) + ∇ b p + ( ∇ a F ac ) F cb = 0 . (28)In what follows, we will consider the projections of equation (28) along the directions parallel and orthogonalto the flow lines of the fluid. Contracting equation (28) with u b we obtain the conservation equation u a ∇ a ρ + ( p + ρ ) ∇ a u a − u b F cb ( ∇ a F ac ) = 0 . (29)Now, in the case of an ideal conducting fluid, where E a = F ab u b = 0, the last term of equation (29) is identicallyzero and the electromagnetic field does not have a direct effect on the conservation equation along the flow lines.Contracting equation (28) with h bc one obtains the so-called Euler equation ( p + ρ ) u a ∇ a u c + h bc ∇ b p + ( ∇ a F ad ) F db h bc = 0 . (30)This last equation can also be written as( p + ρ ) u a ∇ a u b + ( u b u d ∇ d + ∇ b ) p + ( ∇ a F ad ) (cid:0) F db + F ed u e u b (cid:1) = 0 . (31)The last term of equation (31) is identically zero as a consequence of the Maxwell equation so that at the endof the day one obtains the simpler expression ( p + ρ ) u a ∇ a u c + h bc ∇ b p + ( ∇ a F ad ) F cd = 0 . Now, using equation(21) in the above one obtains that n c ( n b ∇ b p ) + ( ρ + p )( (cid:107) A n c + ⊥ A c ) + N bc ∇ b p − (cid:107) E ( n c (cid:37) C + (cid:107) ju c ) − ⊥ E c (cid:37) C − ⊥ j d ⊥ E d u c − (cid:107) B(cid:15) cd ⊥ j d − ( (cid:15) cdu ⊥ j d − (cid:107) j(cid:15) cu ) ⊥ B u = 0 , . (32)Consequently, projecting equation (32) along n a we obtain a constraint equation for the pressure p . Namely,one has that ˆ p + ( p + ρ ) A − (cid:37) C (cid:107) E − (cid:15) de ⊥ j d ⊥ B e = 0 , ˙ ρ + ( ρ + p )Θ − (cid:107) E (cid:107) j − ⊥ j c ⊥ E c = 0 . (33)where (29) has been considered for the conservation equation. In the particular case of a LSS, equation (33)read ˙ ρ + ( ρ + p )Θ − (cid:107) E (cid:107) j = 0 , ˆ p + ( p + ρ ) A − (cid:37) C (cid:107) E = 0 . (34) Evolution equation for the electric part of the Weyl tensor
Using the expressions obtained in the previousSection it can be verified that the evolution equation for the electric part of the Weyl tensor takes the form (cid:107) ˙ E = −
12 ( ρ + p )Σ + (cid:0) Σ − Θ (cid:1) (cid:107) E + 3 (cid:107) B ξ + (cid:15) cd δ c B d + (cid:15) cd ⊥ B db ζ bc +Σ a E a − Σ ec ⊥ E ce + 2 A c (cid:15) cd ⊥ B d + Ω c (cid:15) cd E d + 2 E c α c + ˙F em , (35)where ˙F em = n a n b (cid:0) − σ ab ( E + B ) − ˙ P ab − Θ P ab − D (cid:104) a G b (cid:105) − A (cid:104) a G b (cid:105) − σ c (cid:104) a P b (cid:105) c + (cid:15) cd (cid:104) a ω c P db (cid:105) (cid:1) . (36)Now, following Clarkson’s procedure for the 1 + 1 + 2-decomposition of a LSS, equation (35) becomes (cid:107) ˙ E = (cid:0) Σ − Θ (cid:1) (cid:107) E + 3 ξ (cid:107) B − ( ρ + p )Σ − ( (cid:107) E + (cid:107) B ) (cid:0) Σ2 − Θ3 (cid:1) + 13 u a ∇ a ( (cid:107) E + (cid:107) B ) , while the constraint equation is given by (cid:107) (cid:98) E + 32 (cid:107) E Φ −
13 ˆ ρ − (cid:107) B − F em = 0 , F em = 12 (cid:20)(cid:0) (cid:107) E + (cid:107) B (cid:1) φ + (cid:107) E (cid:99) (cid:107) E + (cid:107) B (cid:99) (cid:107) B (cid:21) . (37)Using, respectively, the evolution equations (26) and the equation for the electric parts of the Maxwell fieldmagnetic field we obtain u a D a E + Ej + Θ (cid:2) ( B + E ) + E (cid:3) − B ξ + (cid:2) ( p + ρ ) − ( B + E ) − E (cid:3) Σ = 0 . (38)Moreover, the constraint equation is given by n a D a E − n a D a p − n a D a ρ = A ( p + ρ ) − (cid:2) ( B + E + 3 E )Φ − B Ω (cid:3) . (39)Notice that this last equation involves also the derivatives of ρ and the pressure p projected along the radialdirection. Evolution equation for the magnetic part of the Weyl tensor
The evolution equation for the magnetic partof the Weyl tensor is given by (cid:107) ˙ B = (cid:0) Σ − Θ (cid:1) (cid:107) B − (cid:107) E ξ + 2 ¯ B c α c − (cid:15) cd δ c E d + Σ a B a − Σ ac ⊥ B ac − (cid:15) cd A c E d + (cid:15) cd B d Ω c − (cid:15) cd ⊥ E bd ζ cb + ˙B em , where ˙B em = n a n b (cid:0) curl P ab + ω (cid:104) a G b (cid:105) + (cid:15) cd (cid:104) a σ cb (cid:105) G d (cid:1) . (40)Thus, in a LSS one obtains the propagation equation and the associated constraint equation (cid:107) ˙ B = (cid:0) Σ − Θ (cid:1) (cid:107) B − ξ (cid:107) E − ξ (cid:0) (cid:107) E + (cid:107) B (cid:1) , (cid:107) (cid:98) B + (cid:107) B Φ + ( ρ + p )Ω + 3Ω (cid:107) E − Ω( (cid:107) E + (cid:107) B ) = 0 . (41) Evolution equations for the kinematic quantities
In order to discuss the evolution equations it is convenient toconsider the following third rank tensors R ( n ) abc ≡ ∇ [ a ∇ b ] n c − R abcd n d = 0 , R ( u ) abc ≡ ∇ [ a ∇ b ] u c − R abcd u d = 0 . The required evolution and associated constraint equations are then obtained from the projection of these zeroquantities along the tensors ( u a , n a , (cid:15) ab , (cid:15) abc ), we are able to find out the evolution and constraint equations forthe kinematic quantities.From the constraint R ( u ) abc u a (cid:15) bc = 0 we readily deduce the following evolution equation for the vorticity projectedalong the radial direction ˙Ω − A ξ − (cid:15) bc D b A c = (cid:0) Σ − Θ (cid:1) Ω + ( ⊥ Σ b + α b ) ⊥ Ω b . In the LSS case it readily reducesto ˙Ω = A ξ + (cid:0) Σ − Θ (cid:1) Ω , n a D a Ω + (Φ − A )Ω = 0 . (42)we obtain (in a LSS) the constraint equation through the contraction R ( u ) abc (cid:15) abc = 0. Computing u a N bc R ( n ) abc = 0we find an evolution equation for Φ, namely,˙Φ − (cid:0) Θ − Σ (cid:1)(cid:0) A − Φ (cid:1) − ξ Ω = δ c α c + A d (cid:0) (cid:15) dc ⊥ Ω c − ⊥ Σ d − a d + α d (cid:1) − ζ cd ⊥ Σ cd + a c (cid:0) ⊥ Σ c − (cid:15) cu ⊥ Ω u (cid:1) + (cid:15) uv ⊥ E u ⊥ B v which, for a LSS, reduces to˙Φ − (cid:0) Θ − Σ (cid:1)(cid:0) A − Φ (cid:1) − ξ Ω = 0 , n a D a Φ +
E − Θ − ξ + ρ − ΘΣ + Σ + Φ = 0 (43)constraint equation for Φ in a LSS is found by the contraction R ( n ) abc n a N bc . Computing the contraction u a (cid:15) bc R ( n ) abc = 0 one obtains˙ ξ = B − (cid:0)
13 Θ − Σ (cid:1) ξ + (cid:15) ab δ a α b + (cid:0) A − ΦΩ (cid:1) Ω − (cid:15) ab ζ au ⊥ Σ bu + (cid:2) (cid:15) ab ( ⊥ Σ b + α b ) + ⊥ Ω a (cid:2)(cid:0) a a + A a (cid:1) , which, for a LSS reduces to˙ ξ = B − (cid:0) Θ − Σ (cid:1) ξ + (cid:0) A − Φ (cid:1) Ω , n a D a ξ + ξ Φ − (cid:0) Θ + Σ (cid:1)
Ω = 0 , (44)the associated constraint is given, for a LSS, by the constraint R ( n ) abc (cid:15) bca = 0. Through similar calculations, oneobtains from the contractions R ( u ) abc u a g ab = 0 and R ( n ) abc u a n b u c = 0 two alternative equations for the expansionscalar: u a D a Θ − n a D a A + 12 ( B + E ) + (cid:0) p + ρ (cid:1) + Θ + Σ − A − A Φ − = 0 , (45a) u a D a Θ + u a D a Σ − n a D a A = A − ( B + E ) − E − (cid:0) p + ρ (cid:1) − (cid:0) Θ + Σ (cid:1) . (45b)Finally, in the sequel it will be useful to consider the following equations obtained, respectively, from thecontractions R ( n ) abc (cid:15) ab u c = 0, R ( u ) abc u a N bc = 0, R ( u ) abc n a N bc = 0: B − ξ Σ + (cid:0) A − Φ (cid:1) Ω = 0 , (46a) u a D a Θ − u a D a Σ + p + ρ + 2 (cid:0) Θ − Σ (cid:1) − E − A Φ − = 0 , (46b) n a D a Θ − n a D a Σ − ΣΦ − ξ Ω = 0 . (46c) IV. REMARKS ON THE THERMODYNAMICAL QUANTITIES
In the present article we consider a one species particle fluid ( simple fluid ), and denote, respectively, by n, s, T the particle number density , the entropy per particle and the absolute temperature , as measured by comoving observers.We also introduce the volume per particle v and the energy per particle e via the relations v ≡ /n and e ≡ ρ/n . Interms of these variables the first law of Thermodynamics, d e = − p d v + T d s , takes the formd ρ = p + ρn d n + nT d s. (47)Assuming an equation of state of the form ρ = f ( n, s ) ≥
0, one obtains from equation (47) that p ( n, s ) = n (cid:18) ∂ρ∂n (cid:19) s − ρ ( n, s ) , T ( n, s ) = 1 ρ (cid:18) ∂ρ∂s (cid:19) n . (48)Assuming that ∂p/∂ρ > speed of sound ν s = ν s ( n, s ) by ν s ≡ (cid:18) ∂p∂ρ (cid:19) s = nρ + p ∂p∂n > . (49)Since we are not considering particle annihilation or creation processes we consider the equation of conservation ofparticle number: u a ∇ a n + n ∇ a u a = 0 . (50)Combining this equation with equations (29) and (47) we obtain u a ∇ a s = 1 nT u b F cb ∇ a F ac . (51)0Where n is subject to equation (50) and T is given by equation (48), thus,using the 1 + 1 + 2 decompositionin aLSS one obtain the couple of equations ˙ s = + 1 nT (cid:107) E (cid:107) j ˙ n + n (cid:107) Θ = 0 . (52)In the case of an infinitely conducting plasma, where the last term of equation (47) vanishes, equation (51) describesan adiabatic flow —that is, u a ∇ a s = 0, so that the entropy per particle is conserved along the flow lines. A particularcase of interest is when s is a constant of both space and time. In this case the equation of state can be given in theform p = p ( ρ ). In the following for convenience, it is assumed that ∇ a p = v s ∇ a ρ where we have defined v s ≡ ν s .Further discussion on the thermodynamical quantities and assumptions can be found in Appendix C. V. VARIABLES AND EQUATIONS
In this subsection we summarise our analysis so far. The evolution equations for the independent components ofthe vector-valued unknown (cid:126)v = ( ρ, s, n, Ω , ξ, Φ , E, B, B , E , Θ , Σ , A ) , (53)describing the evolution of a LSS in equations (34), (52), (42), (44), (43), (26),(41), (38), (45a) and (45b) are givenby ˙ ρ = − ( ρ + p )Θ + E (cid:107) j, (54a)˙ s = 1 nT E (cid:107) j, (54b)˙ n = − n Θ , (54c)˙Ω = A ξ − ΘΩ + ΣΩ , (54d)˙ ξ = B − (cid:0) Θ − Σ (cid:1) ξ + (cid:0) A − Φ (cid:1) Ω , (54e)˙Φ = (cid:0) Θ + Σ (cid:1)(cid:0) A + Φ (cid:1) − ξ Ω = 0 , (54f)˙ E = 2 ξB − (cid:0) Θ − Σ (cid:1) E − (cid:107) j, (54g)˙ B = − ξE − (cid:0) Θ − Σ (cid:1) B, (54h)˙ B = (cid:0) Σ − Θ (cid:1) B − ξ E − ξ (cid:0) E + B (cid:1) , (54i)˙ E = (cid:0)
32 Σ − Θ (cid:1) E + 3 ξ B − ( ρ + p )Σ + (cid:18) Σ2 − Θ3 (cid:19) ( E + B ) − E (cid:107) j, (54j) u a D a Θ − n a D a A = − ( B + E ) − (cid:0) p + ρ (cid:1) − Θ − Σ + A + A Φ + 2Ω , (54k) u a D a Σ + u a D a Θ − n a D a A = A − ( B + E ) − E − (cid:0) p + ρ (cid:1) − (cid:0) Θ3 + Σ (cid:1) . (54l)In the subsequent stability analysis, we will assume a vanishing radial acceleration. It is worth pointing out that thederivative along the radial direction of the acceleration A appears in the evolution equations for Θ and Σ. In SectionV A and in Appendices A,B and C we provide a suitable propagation equation for this variable so as to complete asymmetric hyperbolic evolution system.As discussed in more detail in Section IV, the system will be supplemented by an equation of state of the form p = p ( ρ ). Furthermore, it will be necessary to specify the form of the conduction current, j a . Consistent with Ohm’slaw , we assume a linear relation between the conduction current j a and the electric field. More precisely, we set j a = σ ab E b , where σ ab denotes the conductivity of the fluid (plasma). We will restrict our attention to isotropic fluidsfor which σ ab = σ J g ab , so that J a = (cid:37) C u a + σ J E a , (55)with σ J the electrical conductivity coefficient . In terms of a 1 + 1 + 2-decomposition one has that in the particularcase of a LSS (cid:107) j a = (cid:107) Eσ J n a . (56)The system (54a)-(54j) is complemented by the constraint equations for the components of the vector-valued function (cid:126)v . Some of these constraint equations were obtained in the previous Sections. Further discussion can be found in the1Appendices to this article. As the primary aim of this work is the study of the time evolution and the analysis of itsstability, the discussion of the constraint equations will remain in a secondary level. In what follows, it is assumedthat the constraint equations are satisfied at all time. This assumption can be removed by fairly general arguments—see [23, 24]. A. Re-parametrised set of evolution equations
In what follows we assume that ∂ a v s = 0 and we introduce the variables Q and T via the relations Q ≡ T + 3Σ , T ≡ Θ − Σ so that Σ = ( Q − T ) , Θ = ( T + Σ) . (57)In terms of these one obtains a vector-valued unknown with components given by v ≡ (cid:0) ρ, E, B, E , B , Q, T, ξ, Φ , Ω , A (cid:1) . (58)and equations (54a), (54g), (54h), (54j), (54i), (54f), (54e), and(54d) take the form˙ ρ − Ej + 32 ( p + ρ ) (cid:0) ( Q − T ) + T (cid:1) = 0 , (59a)˙ E − Bξ + ET + j = 0 , (59b)˙ B + BT + 2 Eξ = 0 , (59c)˙ E − B ξ + ( B + E ) T + E T + 23 Ej + ( p + ρ )( Q − T ) = 0 , (59d)˙ B + B T + 3 E ξ + ( B + E ) ξ = 0 , (59e)˙Φ − A T − ξ Ω + T Φ = 0 , (59f)˙ ξ − A Ω + ( ξT + ΦΩ − B ) = 0 , (59g)˙Ω − A ξ + T Ω = 0 . (59h)Moreover, from equation (46a) one finds that B − A ξ + (2 A −
Φ)Ω = 0 , while from equations (54l) and (54k) oneobtains that ˙ T − A Φ − E + p + ρ + 12 T − = 0 , (60a)˙ Q − A − A + B + E + 2 E + p + 13 ρ + 12 Q = 0 , (60b)4 ˙ A v s − Q + 4 E ˆ jp + ρ − E ˙ (cid:37) C v s ( p + ρ ) − v s ( p + ρ ) (cid:18) A (cid:2) Q ( v s −
1) + 2
T v s (cid:3) ( p + ρ ) − A Ej (1 + 2 v s )+ (cid:37) C (cid:2) EQ (2 + v s ) + 2( ET v s + 2 Bξ ) (cid:3) + v s ( p + ρ ) (cid:2) Φ( Q − T ) + 4 ξ Ω (cid:3) − j (cid:2) (cid:37) C (1 + v s ) + v s (2 B Ω − E Φ) (cid:3)(cid:19) = 0 . (60c)As already discussed in Sections IV and V) the fields p and (cid:37) C can be regarded as functions of ρ , and while j c is regarded as a function of the electric field E . Equations (59a)-(60b) do not contain explicit dependence on thederivative of the matter density along the radial direction. In equation (60c) the equation of state and the constraintequation have been used to obtain a simpler expression. The equation for the radial acceleration can be recovered aftersome further calculations —see the Appendix Section. In order to simplify the stability analysis of the subsequentSections and to be able to extract some information from the rather complicated perturbed equations we will focuson the null radial acceleration case. Summarizing, our analysis has lead to a symmetric hyperbolic system consistingof eleven evolution equations for eleven scalar variables of the form A t ∂ t v − A j ∂ j v = Bv , (61) VI. PERTURBATIONS
We will perform a perturbation to first order of the variables v so that v (cid:55)→ ˚ v + (cid:15) ˘ v with the parameter (cid:15) setting theperturbation order that controls the size of perturbation and ˘ v describing the perturbation of the background solution2˚ v . Thus, assuming that the background variables ˚ v satisfy the evolution equations, and defining (d v ( a ) / d (cid:15) ) (cid:15) =0 ≡ ˘ v ( a ) ,we can obtain from equations (59a)-(59h) the following equations for the linear perturbations:˙˘ ρ − (cid:18) ˘ j ˚ E + ˘ E ˚ j − (cid:2) (˘ p + ˘ ρ )( ˚ Q + 2˚ T ) + ( ˘ Q + 2 ˘ T ) (˚ p + ˚ ρ ) (cid:3)(cid:19) = 0 , ˙˘ E + ˘ j − ξ ˚ B + ˘ T ˚ E + ˘ E ˚ T − B ˚ ξ = 0 , ˙˘ B + ˘ T ˚ B + 2 ˘ ξ ˚ E + ˘ B ˚ T + 2 ˘ E ˚ ξ = 0 , ˙˘ E − ξ ˚ B + ˘ T (cid:104)
3( ˚ B + ˚ E + 3˚ E ) − (˚ p + ˚ ρ ) (cid:105) + ˘ j ˚ E + ˘ E (2˚ j + 3 ˚ E ˚ T ) + ˘ Q (˚ p + ˚ ρ ) + (˘ p + ˘ ρ )( ˚ Q − ˚ T )+ ˘ E ˚ T − B ˚ ξ + ˘ B ˚ B ˚ T = 0 , ˙˘ B + 32 ˘ T ˚ B + ˘ ξ ( ˚ B + 3˚ E ) + ˘ B ˚ T + 3 ˘ E ˚ ξ + 2 ˘ B ˚ B ˚ ξ = 0 , ˙˘Φ − ˘ T ( ˚ A − ˚Φ) − ˘ A ˚ T + ˘Φ˚ T − ξ − ξ ˚Ω = 0 , ˙˘ ξ − ˘ B − ˘Ω( ˚
A − ˚Φ) + ˘ ξ ˚ T + ˘ T ˚ ξ − ˘ A ˚Ω + ˘Φ˚Ω = 0 , ˙˘Ω − ˘ ξ ˚ A + ˘Ω˚ T − ˘ A ˚ ξ + ˘ T ˚Ω = 0 . (62)In addition one obtains from (46a) the linearised constraint ˘ B + ˘ A (2˚Ω − ξ ) − A ˘ ξ − ˘Φ˚Ω + ˘Ω(2 ˚ A − ˚Φ) = 0 , whileperturbing to first order equations (60a)-(60b)) describing the evolution of the variables T and Q we obtain˙˘ T − ˘ E + ˘ p + ˘ ρ − ˘Φ ˚ A + ˘ T ˚ T − ˘ A ˚Φ − , ˙˘ Q − A + 2 ˘ E + ˘ p + ˘ ρ − A ˚ A + 2 ˘ B ˚ B + 2 ˘ E ˚ E + ˘ Q ˚ Q = 0 . Similar calculations render a complicated linearised equation for the perturbation of acceleration A . Thus, in order toundertake a stability analysis we make some assumptions about the configuration so as to obtain a simplified system.First, we will analyse the set (59a)-(60b)) with the assumption A = 0. Thus, we do not consider any more equation(60c). Accordingly, we consider the simplified vector-valued unknown v ≡ ( ρ, E, B, E , B , Q, T, ξ, Φ , Ω) . A. Remarks on the set of equations
We consider the set of linearised evolution equations in the form˚ A t ∂ t ˘ v − ˚ A j ∂ j ˘ v = ˚ B ˘ v . (63)General theory of linear of partial differential equations shows that systems of this form can either converge (in anoscillating manner or exponentially) to constant values, have an asymptotic exponential (oscillating) decay or beunstable —see e.g. [45]. In the rest of this article, we look at solutions in which the perturbations are stabilised onconstant values or decay. In order to study the asymptotical behaviour of the characteristic solutions of perturbedsystem it is necessary the study of the eigenvalues of the matrix ˚ B as t → ∞ . Notice that in the the matrix-valuedfunction ˚ B are, in general, functions of the coordinates. Under these circumstances, the core of the stability analysisconsists on analysing whether the matrix ˚ B satisfies some appropriate relaxed stability eigenvalue conditions . Moreprecisely, we study then the eigenvalue problem of ˚ B by looking at the sign of the real parts of the correspondingeigenvalues λ i (or modes of the perturbed system ). A similar approach has been used in [24] —see also discussion in[42]. General theory concerning the case where the coefficients linearised system are constant is discussed in [45]. Thecase of systems with vanishing eigenvalues has been discussed in [46–48] while for a general discussion on the timedependent case and the case of non constant matrix coefficient see [24].Due to the form of the matrix ˚ B associated to the present problem, a fully analysis of the stability properties ofthis system is extremely cumbersome. We proceed by analysing the sign of the eigenvalues using the fact that asufficient condition for the instability of the system is the existence of al last one positive-real part eigenvalues. Fora more extended discussion concerning the requirements of the so-called of the relaxed stability eigenvalue conditionsee [46, 48]. This simple requirement will provide a immediate way to show the conditions where instability occurs—namely, if there is at least one i such that Re ( λ i ) >
0. Crucially, in the discussion of the stability one does not needexplicitly to compute the eigenvalues of ˚ B —the sign of the real part of eigenvalues can be determined for examinationof the structure of the matrix ˚ B . Given the characteristic polynomial P (˚ B )( x ) = (cid:80) ni =0 c i x i of the matrix ˚ B , we will3make use of the following well-known results on the roots sign: the trace of a matrix is equal to the sum of itseigenvalues, hence if the trace of the matrix is positive then the system is unstable . Now, recall that the determinantof a matrix is the product of the eigenvalues. Accordingly, a necessary condition for stability is that det ˚ B (cid:54) = 0and ( − n det ˚ B > n is the number of distinct eigenvalues. Other stability criteria based on inspecting thecharacteristic polynomial are the Routh-Hurwitz criterion and the
Li´enard-Chipart theorem —see [42–44]. Whereverpossible, we will also make use of the Descartes criterion to determine the maximum number of positive and negativereal roots of the polynomial P (˚ B )( x ) for c i ∈ R . In particularly simple cases one can exploit a generalisation of theDescartes rule considering the Routh-Hurwitz criterion to determine the number of roots with positive and negativereal part of P (˚ B )( x ) by constructing the associated Routh matrix .In the case under consideration the matrix ˚ B has eleven columns and the coefficients of the characteristic polynomialare, unfortunately, large expressions without an obvious structure. Therefore, we approach the analysis of the stabilityproblem making the following assumption: ˚ j = ˚ Eσ J , ˘ j = ˘ Eσ J , ˚ p = ˚ ρv s , ˘ p = ˘ ρv s . Thus, we take the squareof the sound velocity v s = ν s to be constant and assume Ohm’s law j = Eσ J with ˚ σ j ≡ ˘ σ j and restrict our attentionto the case of null radial acceleration —i. e. A = 0.Before proceeding to the analysis of particular cases, it is convenient to note here some general properties ofthe system of symmetric hyperbolic equations under consideration. The matrix ˚ B can be explicitly written as the10 × −
12 (1+ v s ) ( ˚ Q +2˚ T ) σ J ˚ E −
12 (1+ v s )˚ ρ − (1+ v s )˚ ρ − σ J − ˚ T ξ − ˚ E B − ξ − ˚ T − ˚ B − E −
16 (1+ v s ) ( ˚ Q − ˚ T ) −
13 ˚ E ( σ J +3˚ T ) − ˚ B ˚ T − T ξ −
16 (1+ v s )˚ ρ [ (1+ v s )˚ ρ − ( ˚ B +3˚ E +˚ E )] B − E ˚ ξ − B ˚ ξ − ξ − T − B − ˚ B − E− ˚ E − − v s − E − B − − ˚ Q − − v s − ˚ T − ˚ ξ − ˚ T − ˚Ω2 − ˚Φ20 0 0 0 0 0 − ˚Φ2 2˚Ω − ˚ T ξ − ˚Ω 0 0 − ˚ T . (64) VII. DISCUSSION ON THE NONLINEAR STABILITY OF THE SYMMETRIC HYPERBOLIC SYSTEM
In this Section we discuss necessary conditions for the nonlinear stability of the configuration under consideration.We focus on the reference solution described by matrix ˚ B in (64). We present the analysis of the system stabilityfor some special configurations defined by proper assumptions on the kinematic variables given by the 1 + 1 + 2-decomposition. Each case is then identified according to the assumptions that characterise it —the fields assumedto vanish are indicated in parenthesis. We discuss in details the instability condition for the complete list of cases.Assuming the radial acceleration and other kinematic variables to be zero, we explore the implications of the assump-tion on the set up configuration and its stability. A first insight into the stability of the system can be inferred fromthe trace in terms of Σ and Θ as Tr ˚ B = 6˚Σ − (7 + v s )˚Θ − σ J . (65)It then readily follows that Tr ˚ B > not stable ifΣ > (7 + v s )Θ + σ J . (66)In the remainder of this article, we proceed to a more detailed analysis in which we classify the results according tosuitable assumptions on the remaining system kinematic variables ( T, Q, Φ , ξ, Ω). It should be pointed out, for ease ofreference, that the condition T = 0 is equivalent to the relation Θ = Σ between the expansion Θ of the 3-sheets andthe radial part of the shear of the 3-sheet Σ. Similarly, the condition Q = 0 is equivalent to Θ = − Σ, while the twoconditions T = 0 and Q = 0 imply Σ = 0 and Θ = 0. We also notice from the Maxwell equations that the evolutionof the electric (respectively, magnetic) field is coupled to the magnetic (respectively, electric) field via the twisting ofthe 2-sheet. Hence, when ξ = 0 the two evolution equations decouple and evolve with the only common dependenceon T .4 A. Preliminaries on the symmetries of the system and classes of solutions
The configuration can expand or contract along the direction of symmetry, according to Θ > < A = 0. The limit situation of null radial accelerationhas been adopted to simplify the analysis of the unstable modes, in this scheme we are able to provide a completeclassification of the stable and unstable solutions for the Einstein-Maxwell-Euler system, considering appropriaterestrictions for the set of the perturbed equations. All the kinematic quantities and fields are reduced to scalars bythe projections in the radial direction and on the orthogonal plane by the projector N ab , the set of scalars define thevector variable −→ v introduced in Eq. (58). Any further restrictions on the system, as the vanishing of other dynamicalvariables: Q i ∈ −→ ν , −→ ν ≡ { T, Q, Φ , ξ, Ω } ⊂ −→ v (67)leads to particular different solutions of the Einstein-Maxwell-Euler equations with the new symmetry conditions. Theself-gravitating systems will be especially constrained by the couple ( E , B ), in the gravitating systems the backgroundgeometry is assumed to be in one of the classes of solutions including the new symmetries. The vector −→ ν is actually arestriction of the vector variable (cid:126) v , where the fluid four-velocity u a defines the metric in its 3 + 1 form and supplies theprojected components of the Weyl tensor together with the direction of symmetry ( ξ, Φ) related through Eq. (9), andthe variation of the symmetry direction n a , used to construct the metric tensor in its 2 + 1 + 1 form. The vanishing ofthe Q i elements implies further restrictions involving the annihilation of other quantities, or the constance of Q j (cid:54) = Q i during the evolution along u a . We can then refer to the general conditions C : C ( Q i = 0) : ( E , B , E , Ω) (cid:55)→ (cid:16) p + ρ (cid:17) (68)( Q j , B ) = 0 (69)˙ Q j = 0 . (70)Conditions in Eqs. (69), specified in Table (I), define five principal classes of solutions according to the assumption ofnull Q i ∈ −→ ν . Table (II) shows the sub-classes of solutions constructed by the vanishing of a couple of scalars ( Q i , Q j ).Conditions (68) are listed in Table (III), and state the relationship between the remaining field variables of the vector −→ v i.e. the couple ( E , B ), the electromagnetic fields ( E, B ), and the matter density ρ . Finally, condition (70) is madeexplicit in Table (IV). Tables (I), (II),(III) and (IV) characterize entirely the system throughout its sub-configurationsruled by the system symmetries, providing a complete classification of the solutions. An analysis and a generaldiscussion of the solutions of the Einstein-Maxwell system in terms of the magnetic and electric parts of the Weyltensor can be found in [39]. The stability analysis will be performed on the systems with A = 0, on the five classesof configurations and their sub-classes. The I -class with ( A T ) and II -class with ( A Q ) are particularly significant.Systems T = 0 are characterized by ΘΣ >
0, i.e. the sign of the expansion is concordant with that of the radial shear,then positive for expanding systems, along the direction of symmetry, or negative for contracting systems, see alsoFig. (1). On the other side systems Q = 0 correspond to the case ΘΣ <
0. The relative sign of the scalars of the radialshear and expansion is a significant element affecting the stability of the system and the equilibrium configurations:even in the general case where the only assumption on the system is the null radial acceleration A = 0, the balance ofthe contributions given by the radial shear in the systems in expansion (Θ >
0) or in contraction (Θ < σ J , v s ). The magnetic field, although it is not constant in the time alongthe direction of symmetry, has no specific role to establish the stability of this model. It is possible to show that forthe case ( A T ) there is always a threshold for the emergence of unstable phases, while in the case of the expandingsystems with negative shear, the difference in sign between the two scalar does not involve any instability threshold.The conditions for the equilibrium of these solutions involve quantities Q i ∈ (cid:126)ν exclusively related to the fluiddynamics such as A T ξ or A ξ Ω and imply a serious constraint on the background, for that instability is in some casescertainly verified assuming that B = 0. Table (III) shows conditions C in Eqs. (68): only the electric part of the Weyltensor is determined by the matter fields and the vorticity. In general the classes ( A T ) and ( A Q ) do not impliesthe vanishing of others variables by conditions C in Eqs. (69): that is, the systems do not require any additionalsymmetry as the initial data Q j = 0 for the solutions (68). Remarkably the ( A QT ) class corresponds to the caseof null shear and null radial expansion. Tables (I) and (II) show that the couples ( ξ, Ω) and ( ξ, Φ) are related. Themagnetic field must be constant along the fluid flow for the solution ( A , T, ξ ). and conversely B = 0 where ξ is zero.The solutions T = 0 have radial vorticity Ω constant during the evolution of the system as shown in Table (IV). Ifthe radial vorticity is initially zero, then the expansion Φ of the 2-sheet generated by n a is constant along u a and theelectric part of the Weyl tensor is entirely determined by the matter field ρ . If initially Φ = 0 then ξ = 0 and therefore B = 0, providing finally a no vacuum solution with null magnetic component of the Weyl tensor, Table (III).5 A A T A Q A ξ I -class A T (cid:88) (cid:88) (cid:88) II -class A Q A T Q (cid:88) (cid:88)
III -class A Φ = A ΦΩ A T Φ = T + A ΦΩ A Q Φ = Q + A ΦΩ (cid:88) A Φ ξ B T + A Φ ξ B Q + A Φ ξ B IV -class A ξ A T ξ = T + A Φ ξ B A Qξ (cid:88) V -class A Ω A T Ω A Q Ω A ξ Ω = A ξ Ω B TABLE I:
Null radial acceleration A : the five principal classes of solutions for the linear stability the problem, identified by canceling theradial acceleration and a scalar quantity of the set Q i ∈ (cid:126)ν = { T, Q, Φ , ξ, Ω } and B , the magnetic part of the Weyl tensor. The elements ofthe Table explicit conditions C in (69). Thus, for example the III − class, conditions of ( A Φ ) configurations implies a null radial vorticityΩ or the null couple ξ and B . Then we assumed that the system has a further constraint represented by a third zero quantity Q i ∈ (cid:126)ν defining the subclasses. The check marks indicate the cases already treated, the Table exhausts the five classes ( A , Q i ) and the subclasses( A , Q i , Q j ). Particularly the III − class ( A φ ) shows the symmetry between the solutions at zeros ( A Q ), ( A T ) and { ξ Ω B} . Subclasses,( A , Q i , Q j , Q k ) are in Table. (II). A T Q A T ξ A Q Φ A Φ ξ A T Q
Φ = Q + T + A Φ ξ B (cid:88) (cid:88) (cid:88) Q + T + A ΦΩ (cid:88) (cid:88) (cid:88) A T Qξ A T Q (cid:88) (cid:88) A T Q Ω A T ξ
Ω = T + A ξ Ω B A Qξ Ω = Q + Aξ Ω B A Φ ξ Ω = Φ + ξ Ω BA T Qξ
Ω = T + Q + A ξ Ω B TABLE II:
Subclasses ( A , Q i , Q j ) of the five principal classes I − V in Table (I). The Table highlights the role of the zeros couples ( T Q )and and ξ Ω B . B. Analysis of the system stability, general conditions and results
The symmetries of the system and the adapted 1 + 1 + 2 formalism highlight the emergence of states certainlyunstable for the linear (and non linear) perturbation. The unstable phases of the system with A = 0 are primarilyregulated by the expansion (or the contraction) and the shear in the radial direction. These scalars are related by thesound velocity and the conductivity through a relation F : F : (Σ , Θ) (cid:55)→ ( v s , σ J , ρ ) , (71) ρ (cid:55)→ ( v s , σ J ) . (72)Conditions C in Eqs. (68,69,70) describe the system symmetries but not its stability, conversely the relations F inEqs. (71) relate the only kinematic variables of the radial shear and radial expansion, Σ and Θ respectively, to theconstants from the state and constitutive equations: the velocity of sound v s and the conductivity σ J . We emphasizethat Eq. (71) does not involves other dynamic variables such as the Weyl scalar or the electromagnetic contribution.Analogously for the five classes of solutions and their subclasses according to Eqs. (68,69,70) the relationship in (72)provides an upper bound on the matter density as a limiting value determined by the couple ( v s , σ J ). The density,but not to its gradient, is in fact a variable in the system and, since this is and iso-entropic and barotropic fluid, onecan use this to get information on the hydrostatic pressure the system is subjected to. C C ( A , T ): E + 2Ω = p + ρ > C ( A , Q ): E + B + 2 E = − ( p + ρ/ < C ( A , ξ ): B = ΦΩ C ( A , T, Ω): E = p + ρ/ > C ( A ξ, Ω): B = 0TABLE III: Conditions C as in Eq. (68) and Eq. (69). Class ( A , T ): ˙Ω = 0( A , T, ξ ): ˙Ω = ˙Φ = ˙ B = ˙ B = 0( A , T, Ω): ˙Φ = 0TABLE IV: Condition C in Eq. (70). Classes of solutions and null evolution of the kinematic quantities.
1. Discussion on the stability of the system with zero acceleration
A relevant constraint F , for the stability of the system with null radial acceleration, is provided by condition (66)on the trace Tr ˚ B , and it can be expressed in terms of the expansion along the privileged direction as follows( u ) : Θ < Θ u ≡ − σ J κ , κ ≡ v s > , (73)this is a first condition for the instability of the system, providing an upper bound, a threshold for the occurrence ofthe instability, on the expansion or contraction with respect to the shear projected along the direction of symmetryi.e., the system is bound to be unstable if the expansion is, algebraically, smaller then Θ u , function of the radial shear,the sound velocity and the conductivity parameter. For small or zero conductivity the limit Θ u is a fraction of thepositive shear, where the system cannot contract or an instability occurs, for negative shear the contraction is limitedby the value Θ u <
0, for contractions too fast i.e. Θ ∈ ]Θ u ,
0[ the system is unstable according to Eq. (73). We cansummarize the analysis in the as following points:
Systems in contraction Θ < and in expansion Θ > expanding configurations, described the I and IV quadrant of Fig. (1), are favored in their stable phases with respect to the systems in contraction: the preferencefor the stable expansion is especially evident for negative shear, as in IV quadrant which describes also a part ofthe ( A , Q ) systems. In the fourth quadrant the stability conditions are always guaranteed within the condition(73) as it is inf σ J >
0, the condition of the trace in Eq. (73) is not sufficient to establish the emergence ofthe instability for this system. Contrasting verses of the scalars ΘΣ < I quadrant, with positive shear, describing parts of the ( A , T )systems. For slow expansions ˙Θ ≈ − Σ in Fig. (1). Considering the instability for the contracting systems, the II and III quadrants, wecan see that the threshold for the emergence of the instability Θ u decreases in general for positive shear up tothe limiting value Σ ≡ σ J v s , the contraction is canceled and for higher shear the system begins an expandingphase, then the system will necessarily be unstable and the zone of stability will be limited to positive shear inthe range Σ ∈ ]0 , Σ [ and contraction Θ ∈ ]Θ , ≡ − σ J /κ < < Θ , the system with null shearis unstable. This region, however, increases with increasing conductivity, this means that the conductivity actsto stabilize the contracting configurations with positive shear. Systems ( A , Q ) are examples of configurations inthe III quadrant: the systems are stable for high conductivity, negative shear and sufficiently low contractions,therefore the expanding ( A Q ) systems are favored for the stability. In general for systems in contraction, in the II and III quadrants, a high conductivity is required for stable systems of ( A T ) and ( A Q ) classes respectively. The role of the radial shear
In the I and II quadrants the radial shear is positive, an increases of Σ generally actsto favor the instability of the configuration. An increase in magnitude of the negative shear in the III quadrant(ΘΣ > The role of the conductivity
In general an increase of the conductivity parameter particularly in the II and III quadrants, has a stabilizing effect on the system according to the law in Eq. (73). In the
III quadrant, for examplefor systems of the I − class or ( A T ) configurations, where the radial shear sign is equal to the contraction one,the stability region increases for low contractions and also for small conductivity as the shear is high enough inmagnitude; ultimately for high conductivity the system can remain stable even for high radial contractions. Theexpansion phase is stable for negative shear, whereas at low expansion the system is unstable. The minimumthreshold for the emergency of the instability increases with the shear but decreases with the conductivity,which has then a stabilizing effect on the system, for very high conductivity. The stable regions in the Θ − Σplane increase by increasing the shear and therefore the stability is advantaged for high conductivity, for the7expansion phase and for shear negative zero or small conductive. In particular the solutions of the II -class( A Q ) are a wider class of solutions admitting difference possible symmetries according to Table (I), stable forany conductivity at negative shear, while for expanding systems as in the II quadrants the instability regions aremaximum for smaller shear and zero for greater shear, the measure of this region increases with the conductivity.However the turning point of the expansion Θ is regulated by the σ J /κ , for null radial shear Σ it is Θ u = − σ J /κ ,the contraction is maximum as provided by Σ x and Θ x . The role of the velocity of sound
Fig. (2) emphasizes the role of the sound velocity in the determination of theunstable states of the system. The shear and expansion scalars have been normalized for the finite conductivity σ J . For the systems included in the I - II and III quadrants, the velocity of sounds plays a significant role in thedetermination of the unstable phases, while there is no threshold for the emergence of the instability, crossingwith continuity the
III and IV quadrants of the Θ − Σ plane. In the limit of very large v s the configurationstops the expansion in the I quarter or the contraction phase in III one i.e. the threshold approaches Θ u = 0.For approximately null velocity of sound it is Θ u /σ J ≈ − / / /σ J = 1 /
6, then for this configuration, the contracting phases are unstable. The solutionsΘ = 0 and
Σ = 0 belongs to ( A , T , Q ) class. For Σ /σ J > / σ J , v s ). Decreasing the velocity of the sound, the region of Θ − Σ plane, for expanding systemscorrespond to unstable regions for an increase of the positive shear. As the velocity of sound decreases forpositive shear the instability region increases. The instability in the I quadrant is carried out for sufficientlyhigh expansions and shear. The increase of the speed of sound acts to stabilize the expanding system for positiveshear and destabilize the contracting systems at Σ /σ J > /
6, in the window of smaller but positive shear, theincrease of the speed of sound acts to destabilize the system while a decrease of this corresponds to an increasethe stability provided that the contraction is sufficiently small in magnitude or Θ /σ J > − / ( v s + 7).We now focus on the effects the velocity of sound in the case of negative shear: for the expanding system in the IV quadrant, no threshold exists. For contracting systems in the III quadrant, an increase of the speed of soundincreases the stability of the system increasing in magnitude the threshold for very high contractions. A decreaseof the velocity of sound, for contracting systems with negative shear, induces the system unstable even for verysmall expansions. The increase in magnitude of the radial shear increases the system stability. Ultimately thesystem’s equilibrium depends on the velocity of sound in I and III quadrants, that is for concordant shear andexpansion or in the IV quadrant for systems in contraction with positive shear smaller than 1 / σ J . In the limitof very large velocity of sound Θ u > /σ J > /
6, for Σ /σ J ∈ ] − ∞ , /
6[ it is Θ u <
2. Stability of the subclasses with restricted symmetries
Additional restrictions on the system with Σ (cid:54) = 0 and Θ (cid:54) = 0, and therefore on the perturbations modify thecondition (73) on the reduced system as followsΣ ≥ Σ | C , Σ | C ≡ a [( b + cv s )Θ + cσ J ] , a < a (cid:28) , b (cid:29) c, c ≥ , { a, b, c } ∈ N , (74)the quantities { a, b, c } change depending on the class of solutions. In terms of the expansion Θ, analogously to Eq. (73)one has Θ ≤ Θ | C , Θ | C ≡ Σ − acσ J a ( b + cv s ) . (75)then the expansion limit is null Θ | C = 0. For Σ max ≡ acσ J , for smaller shear Σ < Σ max the expanding systemsare always stables while a threshold for the instability of the contracting systems appears: for Θ < Θ | C the systemis unstable. It is noteworthy that the shear limit Σ | C = 0 is cancelled for maximal contractions equal to Θ max ≡− cσ J / ( b + cv s ). More generally, for each subclasses of Table (I) and Table (II) the unstable phase will be regulated bysome limiting extreme values (LEV) for the systems Θ = 0 or Σ = 0 providing the boundary values for the regions8 FIG. 1:
Plot of the limiting expansion Θ u κ , defined in Eq. (73) as function of the shear in the preferred radial direction where κ ≡ v s + 7.per σ J = 0 . σ J = 100 (dashed black line ) Possibly an identical plot for the expansion Θ u in terms of σ J /κ and Σ κ .The instability regions, for Θ < Θ u , are colored: then for σ J = 0 . σ J = 100, gray region is then stable for this system. In the white region the system can bestable. Configurations ( A T ) are in the I and III quadrants (ΣΘ >
0) as in the T = 0 class of solution , for Σ = (2 / A , Q )class for Σ = − (1 / II and IV quadrants (ΣΘ < A = 0) is regulated as in the I and IV quadrants. The contraction in II and III quadrants, where the radialshear is positive and negative, respectively. The classes are summarized in Table (I). The contraction limit vanishes, and the stabilitylimit is only on the radial expansion for high enough shear according to Σ = σ J / FIG. 2:
Plot of the limiting expansion Θ u /σ J , in Eq. (73), versus the radial shear Σ /σ J for v s = 0 .
01 (black line) and v s = 100 (dashedblack line). The Shaded regions marke the sections in the Θ − Σ plane where instability occurs as Θ < Θ u . It is Θ u = 0 for Σ /σ J = 1 / FIG. 3:
Ratios of the limiting extreme values (LEV) Σ (left panel) and (LEV) Θ introduced in Eq. (76), functions of the velocity of sound v s . At at v ∗ s it is Σ x = σ − J ρ s = σ J √ ≈ . σ J while Θ x = 0 . σ J , in v (cid:48) s it is Σ s = Σ x √ ρ/σ J ≈ . √ ρ analogously Θ s = 0 . √ ρ ,in the cross point in ] v (cid:48) s , v ∗ s [ it is ρ s = 0 . σ J where Σ s = 0 . √ ρ in ( LEV ) Σ while the ratio is ρ s /σ J ≈ .
44 in (
LEV ) Θ . of the plane Θ − Σ respectively for the stability of the system: (LEV) Σ : ρ s ≡ σ J (3 + v s ) (1 + 3 v s ) , Σ s ≡ (cid:115)(cid:18)
13 + v s (cid:19) ρ, Σ x ≡ σ J v s ) . (76) (LEV) Θ : ρ s = ρ s | (LEV) Σ , Θ s ≡
32 Σ s , Θ x ≡
32 Σ x . (77)see Fig. (3). A threshold exists for the stability in the systems with the density in the two regimes ρ > ρ s and ρ < ρ s ,according to ( LEV ) Σ or ( LEV ) Θ , with conditions on the radial shear or respectively the expansion or contraction.The details of the analysis are considered in Sec. (D). For the stability in general one has the two following cases: ρ < ρ s ∆ ∈ [ − ∆ x , − ∆ s ] or ] + ∆ s , + ∞ [ (78) ρ > ρ s ] − ∆ s , ∆ s [ or ] + ∆ s , + ∞ [ , ∆ = Σ for ( LEV ) Σ , ∆ = Θ for ( LEV ) Θ (79)each of the two regions are typically regulated by additional conditions on Q i ∈ (cid:126)ν ( LEV ) Σ , for example for the I − class, ( A T ) solution, are determined by a set of conditions on B , ξ and E . We note that the limiting value on thedensity ρ s , increases with the conductivity but decreases with the velocity of sound: there are three ranges of v s tobe considered as detailed in Fig. (3). However the boundary ∆ s depends on the density ρ and therefore differentlyregulated in the three regions. For low density regimes ρ < ρ s , the density function regulates the system stabilityin a reduced region of the ρ − Σ plane through the conditions imposed on shear or the expansion especially for lowvalues. This situation is more evident for low velocity of sounds v < v ∗ s . For larger values the ranges for low densityvariations are restricted to ρ < . σ J . C. Comments on the five principal classes of solutions I -Class ( A T ) : this case is constrained by the condition C ( A , T ) obtained from Eq. (60a). This relation is aconsequence of the condition T = 0, and it will occur in the other subcases characterised by this assumption.The evolution equation for the parallel vorticity Ω is ˙Ω = 0: this quantity remains constant along the fluidmotion, for the perturbed quantity as well as the unperturbed one. Furthermore, the time evolution of thematter density only involves the kinematic variable Q . We can draw some conclusions on the stability ofthe system on the basis of the matrix trace. The system is linearly unstable if with a negative shear onthe radial direction which is bounded by ˚Σ ≤ Σ x . It is worth noting that the limiting case is defined onlyby the constitutive equation, the conductivity σ J , and the equation of state by v s —fixed by the referencesolution. Furthermore, considering the coefficients c = − c = Tr˚ B < stability c <
0. Using the condition C ( A , T ), this condition can be rephrased as the following two alter-natives: (i) ˚ ρ < ˚ ρ s with ˚Σ ∈ ( − ˚Σ x , − ˚Σ s ) ∪ (˚Σ s , + ∞ ) and (ii) ˚ ρ ≥ ˚ ρ s with ˚Σ ∈ (˚Σ s , + ∞ ), where the ( LEV ) Σ hold.0 II -Class ( A Q ) : the condition Θ = − Σ implies that the matter and field variables are related by the constraint C ( A , Q ) —see equation (60b). This relation is a consequence of the assumption Q = 0, and it will occur in theother subcases characterised by this hypothesis. We observe that the assumption Q = 0 implies in particularthat the expansion of the 3-sheets and the radial component of the shear of the 3-sheet must have oppositesign. Studying the problem for the reduced system of nine variables we find that the system is linearly unstable if Tr o B > ≤ Θ | C with ( a = c = 1; b = 9). Again, as in thecase T = 0, we obtain a constraint for the radial expansion that depends on the unperturbed constants ( σ J , v s ). III − Class ( A Φ ) : the assumptions A = 0 and Φ = 0 on equations (59f) and (59g) give rise to two possible subcases:(i) Ω = 0 or (ii) ξ = 0 and B = 0. We consider therefore these two subcases separately:( A Φ Ω ) : the trace implies the following condition of instability ˚Σ ≥ Σ | C with ( a = 2 / , b = 6 , c = 1). Wenote that for the reference solution ˚Θ and ˚Σ are related by the square of the velocity of sound and theconductivity through σ J . However, the sign of the ΣΘ is not constrained by this relation.( A Φ ξ B ) : in this case T is the only kinematical variable involved in the evolution of the electromagnetic field.The criterion on the trace of the matrix ˚ B implies the following condition of instability ˚Σ ≥ Σ | C with( a = 2 / , b = 16 , c = 3).Thus, we finally note that the subcases (i) and (ii) are characterised by similar constraints on ˚Σ and ˚Θ. IV − Class ( A ξ ) : equation (59g) leads to the condition C ( A , ξ ) :. On the other hand the condition of instability onthe trace of the matrix ˚ B leads to the inequality ˚Σ ≥ Σ | C with ( a = 1 / , b = 40 , c = 6). V − Class ( A Ω ) : in this case the system is unstable if ˚Σ ≥ Σ | C with ( a = 1 / , b = 19 , c = 3).We conclude this Section pointing out that the stability of the configuration under consideration is constrained bysimilar relations between the unperturbed radial part of the shear of the 3-sheet ˚Σ and the radial expansion ˚Θ. Theseconstraints only depend on the square of the sound velocity and the conductivity. In Sec. (D) we specialise thisanalysis to consider the subclasses of Tables (I) and (II). VIII. CONCLUSIONS
In this article we have explored the stability proprieties of an ideal, LRS Einstein-Maxwell perfect fluid system.As a first step we have formulated an hyperbolic initial value problem providing a suitable (quasilinear) symmetrichyperbolic system by means of which one can address the nonlinear stability of this system. This first result allowsus to provide a suitable formulation of an initial value problem, that is a necessary issue for the construction of thenumerical solutions, ensuring the local and global existence problems. Moreover the problem of the propagation of theconstraints can be assumed satisfied at all time and then prove this result by fairly general arguments as discussed in[23, 24]. Our analysis is based on a 1+1+2-tetrad formalism. In our calculations certain choices of kinematic propertiesof the configuration were motivated by the necessity to extract some information from the rather complicated equationsgoverning the perturbations. We studied five principal classes of solutions and different subcases, considering systemswith particular kinematic configurations. The assumption of the radial symmetry simplifies the problem and takesgreat advantage from the 1 + 1 + 2-decomposition [25]. Thus, in this work we have proceeded as follows: we wrote the1+1+2-equations for the LRS system using the radial vector, pointing along the axis of symmetry. Then, a discussionon the thermodynamical quantities of the system was been provided. For the linear perturbation analysis it was usefulto introduce a re-parametrised set of evolution equations based of a suitable combination of the radially projectedshear and expansion. The resulting evolution system was then used to analyse the stability problem for small linearperturbations of the background. We presented the re-parameterized set of evolution equations collected in a propersymmetric hyperbolic form and discussed the perturbation to the first order of the variables. We have presented themain results concerning the linear stability of this system in a classification of subcases that constitutes the principalresult of this paper. In the Appendix to this article, we provide also an alterative symmetric hyperbolic system forthe fluid fields within the 1 + 1 + 2-decomposition followed by some general notes on the evolution equations andhyperbolicity considerations. A suitable propagation equation for the fluid radial acceleration is there recovered by theintroduction of a new unknown field corresponding to the derivative of the matter density projected along the radialdirection. suitable field and evolution equations can be obtained for this quantity. The set of evolution equationsis complemented by the constraint and constitutive equations. Restricting our attention to isotropic fluids (entropy1is a constant of both space and time) we considered a one species particle fluid (simple fluid) and we introduced apolytropic equation of state with a constant velocity of sound. In order to close the system of evolution equations itis necessary to specify the form of the conduction current. Accordingly, we assumed the Ohm’s law so that a linearrelation between the conduction current and the electric field, involving a constant electrical conductivity coefficient,holds. We have assumed that the pressure of the fluid p and the charge density (cid:37) C to be functions of the matterdensity ρ , and the charge current j c function of the electric field E . Although viable, the perturbed equation forthe radially projected acceleration A turned out to be a very complicated expression of the other variables and theirderivatives. Thus, in order proceed in the stability analysis, we studied a simplified form system by taking up someassumptions on the configuration. Assuming a null radial acceleration for the reference solution, our analysis isparticularly focused on some specific cases defined by fixing the expansions of the 3-sheets and 2-sheets, the radialpart of the shear of the 3-sheet, the twisting of the 2-sheet and the radial part of the vorticity of the 3-sheet. In thisway we are also able to provide results concerning the structure of the associated LSS taking into considerations allthe different subcases. This analysis constitutes the main result of the paper. In particular, we found that in manycases the stability conditions can be strongly determined by the constitutive equations by means of the square of thevelocity of sound and the electric conductivity . In particular, this is evident for the contracting and expanding LRSconfigurations: a threshold for the emergence of the instability appears in both cases. The conditions bind mainlythe expansion (viceversa contraction) along the preferred direction with respect to different regimes of the radialshears. These results provide in a quite immediate manner information regarding were certainly unstable system.The results for this type of configurations are illustrated in Sec. (VII B) and schematically in Fig. (1) and Fig. (2),in the four fundamental cases emphasizing the role of the couple of parameters ( v s , σ J ) moreover together with theseparameters the relative sign of scalars Σ and Θ plays an essential role in the determination of the unstable phasesof the systems. For expanding configurations with ΣΘ <
0, appears no threshold for the emergence of instability bymeans of condition (66), this does not mean that expanding systems with negative sher are in any condition certainlystable, but further conditions can be provided as we discussed in dealing with different classes of solutions. Theinteresting situation is in the remaining three cases, the role of the velocity of sound and the conductivity acts in adifferent way for the systems in the three regions of the Figs. (1,2), according to our analysis to favor or not certainlyunstable states of the system; in any case there is always is a specific threshold for the contraction or expansion,above which in the first case for contracting systems with a fast contraction rate | Θ | > | Θ u ( σ J , v s ) | > < Θ < Θ u ( σ J , v s ) system is certainly unstable. The other cases and subcases show similar situations where thethreshold provided by the density values, and the couple (Σ , Θ) varies in in different ranges depending on ( ρ, v s , σ J ).The magnetic field does not have a specific role in determining the stability of the system, and the Maxwell field andthe geometrical effects enclosed by the magnetic and electric part of the Weyl tensor.Concerning the methods, the core of the stability analysis is given by the study of the non principal part of thematrix ˚ B of the system using some relaxed stability eigenvalue conditions. However, a full analysis of the stabilityproperties of the system turns extremely cumbersome because of the form of the matrix. Thus, we proceeded with theanalysis of the eigenvalues using an indirect method aimed at determining to know the sign of these. Repeatedly usethe fact that a sufficient condition for the instability of the system is the existence of at last one positive eigenvalues.This simple requirement is nevertheless able to provide a immediate way to show the conditions where the linearinstability occurs. In several places we made use of the Descartes criterion to determine the maximum number ofpositive and negative real roots of the characteristic polynomial. In particularly simple cases one can make use ofthe so-called Routh-Hurwitz criterion to determine the number of roots with positive and negative real part of thepolynomial by constructing the Routh associated matrix. It may be possible to use some of these criteria for somespecific cases and we expect future work to include a development of this paper in this direction. Acknowledgments
DP gratefully acknowledges support from the Blanceflor Boncompagni-Ludovisi, n´ee Bildt and wishes to thank theAngelo Della Riccia Foundation and thank the institutional support of the Faculty of Philosophy and Science of theSilesian University of Opava.
Appendix A: A symmetric hyperbolic system for the fluid fields
In this Section we provide a brief discussion of a construction leading to a suitable propagation equation for theradial acceleration A . In addition, we also provide an alternative set of the propagation equations for the otherkinematic and field variables. A similar construction has been given in a slightly different context in [40].2For the sake of convenience let us define g ab u a u b = u u + u i u i = (cid:15) . In the present case (cid:15) = −
1. Using the identity u a ∇ b u a = 0 we readily conclude that ∇ a u = − u i u ∇ a u i . (A1)Moreover, one has that ∇ c ∇ d u + ∇ c ∇ d u i = − u ∇ c u i ∇ d u i − u i u j ( u ) u ∇ d u i ∇ c u j . From the Bianchi identity and theinhomogeneous Maxwell equation we have that ∇ a F ac = + εJ c . (A2)In this last equation we introduce ε = −
1, to match equations (A2) and (21). In addition, we have the continuityequation for the matter density ρ u a ∇ a ρ + ( p + ρ ) ∇ a u a − εu b F cb J c = 0 , (A3)and the Euler equation ( p + ρ ) u a ∇ a u c − (cid:15)h bc ∇ b p − ε(cid:15)J b F cb = 0 . (A4)This last equation can be alternatively written as( p + ρ ) u a ∇ a u c − (cid:15) ∇ c p + u c u b ∇ b p − ε(cid:15)J b F cb = 0 . (A5)In what follows, we will consider a barotropic equation of state p = p ( ρ ) and introduce here the quantity v s definedby the relation ∇ b p = (cid:16) ∂p∂ρ (cid:17) ∇ b ρ = v s ∇ b ρ. Note that v s = ν s , where ν s is the speed of sound introduced in equation(49). Thus Eq. (A5) is now ( p + ρ ) u a ∇ a u c − (cid:15)v s ∇ c ρ + v s u c u b ∇ b ρ − ε(cid:15)J b F cb = 0 . (A6)Using equations (A3) and (A6), we obtain: ∇ c ρ = (cid:15) ( ρ + p ) v s u a ∇ a u c − (cid:15) ( p + ρ )( ∇ a u a ) u c − ε J b F cb v s + ε(cid:15)u b F bd J d u c . (A7)Using again the normalisation condition we can write u as: u = (cid:15)u = β (cid:112) (cid:15) ( (cid:15) − u i u i ) , where g = η = (cid:15) and β = β ( (cid:15) ) is a sign to be fixed according to the metric signature conventions.Now, using equations (A1) in equation (A3) we find( u ∇ ρ + u i ∇ i ρ ) + ( ρ + p ) (cid:18) ∇ i u i − u i u ∇ u i (cid:19) − εu b F cb J c = 0 , (A8)and from equation (A6) ˜ E (0) ≡ u a ∇ a u − (cid:15)v s ρ + p (cid:0) ∇ c g c (cid:1) + u v s ρ + p u b ∇ b ρ − ε(cid:15)J b F b = 0 , ˜ E ( i ) ≡ u a ∇ a u i − (cid:15)v s ρ + p (cid:0) ∇ c g ci (cid:1) + u i v s ρ + p u b ∇ b ρ − ε(cid:15)J b F ib = 0 . Thus, introducing the zero-quantity: ς i ≡ u i u ˜ E (0) − ˜ E ( i ) = 0 , we obtain the equation ς i = − (cid:18) u a ∇ a u i + u i u j u u u a ∇ a u j (cid:19) + (cid:15)v s ρ + p ∇ c ρ (cid:18) g ci − g c u i u (cid:19) + (cid:15)εJ b (cid:18) F ib − u i u F b (cid:19) = 0 . (A9)Equations(A8),(A9)) constitute, after multiplication by a suitable numerical factor, an hyperbolic system for ρ and u i .3In what follows, It turns to be necessary to introduce here the fields µ a ≡ ∇ a ρ, U ab ≡ ∇ a u b , (A10)where ∇ [ a µ b ] = 0. One readily can verify that U ab u b = 0 , U a = − u i u U ai , and that ∇ a U b = − u i u ∇ a U ib − U a u U b − U ia u U bi , (A11) u b ∇ a U cb = − U cb U ba . (A12)Now, applying a covariant derivative to the equations of motion of the fluid and commuting gives Z cb ≡ (1 + v s ) u b u a ∇ a µ c + ( ρ + p ) (cid:0) u b ∇ a U ac + u a ∇ a U cb (cid:1) − (cid:15)v s ∇ c µ b + W cb = 0 , (A13)where W cb ≡ ( ρ + p ) (cid:0) R adca u d u b + R cabd u d u a (cid:1) + U cb u a µ a (1 + v s ) + u b U ac µ a (1 + v s ) + µ c (1 + v s ) (cid:0) u b U aa + u a U ab (cid:1) +( ρ + p ) (cid:0) U cb U aa + U ac U ab (cid:1) − (cid:15) ∇ c v s µ b − (cid:15) ∇ c (cid:0) J a F ab (cid:1) + u b µ c u a ∇ a v s . Using equation (A12) we obtain u b Z cb = (cid:15)u a ∇ a µ c + (cid:15) ( ρ + p ) ∇ a U ac + X c , (A14)with X c ≡ u b W cb − ( ρ + p ) u a U cb U ba , where u b W cb = (cid:15) (cid:0) ( ρ + p ) R adca u d + µ a U ac (1 + v s ) + µ c (1 + v s ) U aa − µ b u b ∇ c v s − u b ∇ c ( J a F ab ) + µ c u a ∇ a v s (cid:1) . In addition, one has that h db Z cb = ( ρ + p ) u a ∇ a U cd − (cid:15)v s ∇ c µ d + v s u d u b ∇ b µ c + Y cd = 0 , with Y cd = h db W cb + (cid:15) ( ρ + p ) u d U cb U ba u a . Equation (A14) can be written as (cid:15) (cid:2) ( u ∇ µ c + u i ∇ i µ c (cid:1) + (cid:15) ( ρ + p ) (cid:0) ∇ i U ic − u i u ∇ U ic (cid:1) + ˆ X c = 0 , (A15)with ˆ X c ≡ X c − u (cid:0) U a U ca (cid:1) . (A16)Now, the combination ( u i /u ) h b Z cb − h ib Z cb , leads to the equation − ( ρ + p ) (cid:18) u a ∇ a U ic + u i ( u ) u j u a ∇ a U jc (cid:19) + (cid:15)v s (cid:18) g di ∇ d µ c − u i u g d ∇ d µ c (cid:19) + (cid:18) u i u g d − g id (cid:19) Y cd (A17) − ( ρ + p ) u i ( u ) u a U ba U cb + v s u i u b ∇ b µ c ( (cid:15) −
1) = 0 . (A18)Viceversa, the combination ( u i /u ) h b Z cb − h ib Z cb leads to the equation( ρ + p ) (cid:18) u a ∇ a U ic + u i u j u a ∇ a U jc u u (cid:19) − (cid:15)v s (cid:18) g di ∇ d µ c − u i u ( ∇ d µ c ) g d (cid:19) + ˆ Y ic = 0 (A19)where ˆ Y ic ≡ ( ρ + p ) u a u i u u U ba U cb − (cid:18) g d u i u − g di (cid:19) Y cd . Equations (A15) and (A19) constitute a symmetric hyperbolic system for the fields µ a and U ia .4 Appendix B: -decomposition
In this Section we study the 1 + 1 + 2-decomposition of the term Z cb in equation (A13). In what follows we fix (cid:15) = −
1. Recall that up to now ( ρ, µ c , u i , U ab ) have been used as independent variables for our evolution system.We now consider ( ρ, µ c , u i ) and the quantities ( (cid:107) A , ⊥ A a , Θ , Σ , Ω , ⊥ Σ a , ⊥ Ω a , ⊥ Σ ba ) as independent variables, the lastbeing defined by the decomposition of the U ab . This decomposition is obtained from considering (cid:107) A = u b n a U ba , and U ab = − u a n b (cid:107) A + (cid:101) U ab , where (cid:101) U ab ≡ + n a n b (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) + N ab (cid:0) (cid:107) Θ − (cid:107) Σ (cid:1) + (cid:107) Ω (cid:15) ab + ⊥ (cid:101) U ab , with ⊥ (cid:101) U ab ≡ − ua ⊥ A b + n a ( ⊥ Σ b + (cid:15) bc ⊥ Ω c ) + ( ⊥ Σ a − (cid:15) ac ⊥ Ω c ) n b + ⊥ Σ ab . Thus, equation (A13) is now Z cb = (1 + v s ) u b u a ∇ a µ c − (cid:15)v s ∇ c µ b − ( ρ + p ) u c [ u b n a ∇ a (cid:107) A + n b u a ∇ a (cid:107) A ]+( ρ + p ) (cid:2) u b n c n a ∇ a (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) + u b N ac ∇ a (cid:0) (cid:107) Θ − (cid:107) Σ (cid:1) + u b (cid:15) ac ∇ a (cid:107) Ω + n c n b u a ∇ a (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) + N cb u a ∇ a (cid:0) (cid:107) Θ − (cid:107) Σ (cid:1) + (cid:15) cb u a ∇ a (cid:107) Ω (cid:3) + ` Z cb = 0 , with ` Z cb ≡ W cb + ( ρ + p ) (cid:0) [ u b ∇ a ⊥ U ac + u a ∇ a ⊥ U cb ] − (cid:107) A [ u b ∇ a ( u c n a ) + u a ∇ a ( u c n b )] + u b ∇ a ( n a n c + N ac ) (cid:107) Θ+ (cid:107) Σ u a ∇ a (cid:0) n b n c − N bc (cid:1) + u a ∇ a ( n c n b + N bc ) (cid:107) Θ + (cid:107) Σ u b ∇ a (cid:0) n a n c − N ac (cid:1) + (cid:107) Ω [ u b ∇ a (cid:15) ac + u a ∇ a (cid:15) cb ] (cid:1) . Explicitly W cb is given by W cb = − R ca u a u b (cid:0) p + ρ (cid:1) + R bacd u a u d (cid:0) p + ρ (cid:1) + u b (cid:0) v s (cid:1) µ c D a u a + u a µ c D a u b + u a v s µ c D a u b − J a u b D c E a + J a u a D c E b + E b u a D c J a − E a u b D c J a + E b J a D c u a + u b µ a D c u a + u b v s µ a D c u a + pD a u b D c u a + ρD a u b D c u a − E a J a D c u b + u a µ a D c u b + u a v s µ a D c u b + pD a u a D c u b + ρD a u a D c u b − (cid:15) baed (cid:0) J a u e D c B d + B a ( − u e D c J d + J e D c u d ) (cid:1) + u a u b µ a D c v s + µ b D c v s . Taking into account the above expressions, we obtain from the evolution equations the following projections: u b u c Z cb : u b u c Z cb = − u c u a ∇ a µ c − ( ρ + p ) n a ∇ a (cid:107) A + u b u c ` Z cb = 0 . (B1) n b n c Z cb : n b n c Z cb = v s n c n b ∇ b µ c + ( ρ + p ) u a ∇ a (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) + n b n c ` Z cb = 0 . (B2) n c u b Z cb : n c u b Z cb = − n c u a ∇ a µ c − ( ρ + p ) n a ∇ a (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) + n c u b ` Z cb = 0 . (B3) u c n b Z cb : u c n b Z cb = v s u c n b ∇ b µ c + ( ρ + p ) u a ∇ a (cid:107) A + u c n b ` Z cb = 0 . (B4) g cb Z cb : g cb Z cb = (1 + v s ) u b u a ∇ a µ b + ( ρ + p ) n a ∇ a (cid:107) A + v s ∇ c µ c + ( ρ + p ) u a ∇ a (cid:107) Θ + g cb ` Z cb = 0 . (B5) h cb Z cb : h cb Z cb = + v s h cb ∇ c µ b + ( ρ + p ) u a ∇ a (cid:107) Θ + h cb ` Z cb = 0 . (B6) N cb Z cb : N cb Z cb = + v s N cb ∇ c µ b + 2( ρ + p ) u a ∇ a (cid:0) (cid:107) Θ − (cid:107) Σ (cid:1) + N cb ` Z cb = 0 . (B7)5 (cid:15) cb Z cb : (cid:15) cb Z cb = +2( ρ + p ) u a ∇ a (cid:107) Ω + (cid:15) cb ` Z cb = 0 . (B8)The field µ b can be decomposed as µ b = ⊥ µ b + (cid:107) µn b . This expression can be explicitly written as µ a = ( N ba µ b − u a u b ˙ µ b )+ n a n b µ b , thus the previous decomposition can be further refined. Using the torsion-free condition ∇ a µ b = ∇ b µ a weobtain from n c u a ( ∇ a µ c − ∇ c µ a ) = 0, u a ∇ a (cid:107) µ = ⊥ µ c u a ∇ a n c + u a n c ∇ c ⊥ µ a − (cid:107) µn a n c ∇ c u a , (B9)which can be read as a propagation equation for (cid:107) µ .Using equation (B9) we get from equation (B1), ( u b u c Z cb = 0):( ρ + p ) n a ∇ a (cid:107) A = − u c u a ∇ a ⊥ µ c + (cid:107) µ (cid:107) A − u b u c ` Z cb . (B10)From equation (B2), ( n b n c Z cb = 0):( ρ + p ) u a ∇ a (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) + v s n b ∇ b (cid:107) µ = v s ⊥ µ c n b ∇ b n c − n b n c ` Z cb . (B11)From equation (B3), ( n c u b Z cb = 0): u a ∇ a (cid:107) µ + ( ρ + p ) n a ∇ a (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) = ⊥ µ c u a ∇ a n c + n c u b ` Z cb . (B12)Using equation (B9) in equation (B3) one obtains( ρ + p ) n a ∇ a (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) = − u a n c ∇ c ⊥ µ a + (cid:107) µn a n c ∇ c u a + n c u b ` Z cb . (B13)From equation (B4) ( u c n b Z cb = 0) one gets:( ρ + p ) u a ∇ a (cid:107) A = − v s u c n b ∇ b ⊥ µ c − v s (cid:107) µu c n b ∇ b n c − u c n b ` Z cb , (B14)and using equation (B9) in equation (B4)( ρ + p ) u a ∇ a (cid:107) A + v s u c ∇ c (cid:107) µ = v s ⊥ µ b u c ∇ c n b − u c n b ` Z cb . (B15)From equation (B5) ( g cb Z cb = 0) we get:( ρ + p ) u a ∇ a (cid:107) Θ + ( ρ + p ) n c ∇ c (cid:107) A + v s n c ∇ c (cid:107) µ = − (1 + v s ) u b u a ∇ a ⊥ µ b − v s ∇ c ⊥ µ c − µ b u b u a ∇ a v s − (1 + v s ) (cid:107) µu b u a ∇ a n b − v s (cid:107) µ ∇ c n c − g cb ` Z cb . (B16)From equation (B6), ( h cb Z cb = 0):( ρ + p ) u a ∇ a (cid:107) Θ + v s n c ∇ c (cid:107) µ = − v s h cb ∇ c ⊥ µ b − v s (cid:107) µ ∇ c n c − h cb ` Z cb , (B17)From equation (B7), ( N cb Z cb = 0):( ρ + p ) u a ∇ a (cid:0) (cid:107) Θ − (cid:107) Σ (cid:1) = − v s N cb ∇ c ⊥ µ b − v s N cb (cid:107) µ ∇ c n b − N cb ` Z cb . (B18)From equation (B11), ( n b n c Z cb = 0):( ρ + p ) u a ∇ a ( (cid:107) Θ + (cid:107)
Σ) + v s n b ∇ b (cid:107) µ = v s ⊥ µ c n b ∇ b n c − n b n c ` Z cb . (B19)From equation (B6) ( h cb Z cb = 0):( ρ + p ) u a ∇ a (cid:107) Θ + v s n c ∇ c (cid:107) µ = − v s h cb ∇ c ⊥ µ b − v s (cid:107) µ ∇ c n c − h cb ` Z cb . (B20)Finally, we consider the equation ρ + pv s ((B11) − (B20)) , which, explicitly, is given by ( ρ + p ) v s u a ∇ a (cid:107) Σ + ( ρ + p ) n c ∇ c (cid:107) µ = ( ρ + p ) v s (cid:16) v s ⊥ µ c n b ∇ b n c − n b n c ` Z cb + ( v s h cb ∇ c ⊥ µ b + v s (cid:107) µ ∇ c n c + h cb ` Z cb ) (cid:17) , (B21)6and the equation ρ + p v s (B20) or( ρ + p ) v s u a ∇ a (cid:107) Θ + ρ + p n c ∇ c (cid:107) µ = − ρ + p v s (cid:16) v s h cb ∇ c ⊥ µ b + v s (cid:107) µ ∇ c n c + h cb ` Z cb (cid:17) . (B22)For a LSS the equations Equations (B14), (B8), (B12), (B22), (B21) constitute a symmetric hyperbolic sys-tem for the unknowns ( (cid:107) A , (cid:107) Ω , (cid:107) µ, (cid:107) Θ , (cid:107) Σ). In Section (V A) we discussed an alternative set of symmetric hyper-bolic system where the evolution equation for the radial acceleration A is coupled with the evolution of equa-tion for the variable Q ≡ Θ + 2Σ. A symmetric hyperbolic system has, therefore, been given for the vari-ables v ≡ ( ρ, E, B, E , B , Q, T, ξ, Φ , Ω , A ). Equation (60c) for the variable A has been recovered from the equation(B15) − v s (B14). In a LSS equation (B14) is explicitly given by A (1 + v s ) ( Ej − A Θ( p + ρ )) + j(cid:37) C + E(cid:37) C (cid:0) Θ − Σ (cid:1) − Bξ(cid:37) C + ( p + ρ ) u a D a A + µu a D a v s + v s u a D a µ − Eu a D a (cid:37) C = 0 , and (cid:0) Θ + Σ (cid:1) [ A ( p + ρ ) + µ (1 + v s )] − (cid:37) C (cid:2) j + E (cid:0) Σ + Θ (cid:1)(cid:3) + Ej (Φ − A ) + ( p + ρ )ΣΦ+Ω [2( p + ρ ) ξ − Bj ] − En a D a j + ( ρ + p ) (cid:0) n a D a Θ + n a D a Σ (cid:1) + u a D a µ = 0 . It is possible to show, using again the evolution equations and the constraints that the scalars n c n b ` Z cb , h cb ` Z cb , u c n b ` Z cb and n c u b ` Z cb do not contain derivatives of the variables. These terms are discussed with more detail in the followingSection. Appendix C: Evolution equations and hyperbolicity considerations
The system consisting of equations (34), (52), (42), (44), (43), (26), (41), (38), (45a), and (45b) as discussed inSection (V) and equations (B14), (B12), (B21) and B22 for the variables v ( b ) ≡ ( ρ, s, n, (cid:107) Ω , ξ, Φ , (cid:107) E, (cid:107) B, (cid:107) B , (cid:107) E , (cid:107) A , (cid:107) µ, (cid:107) Σ , (cid:107) Θ)can be written explicitly as˙ ρ = − ( ρ + p )Θ + (cid:107) E (cid:107) j, (C1)˙ s = 1 nT (cid:107) E (cid:107) j, (C2)˙ n = − n (cid:107) Θ , (C3)˙ (cid:107) Ω = (cid:107) A ξ − (cid:107) Θ (cid:107) Ω + (cid:107) Σ (cid:107) Ω , (C4)˙ ξ = (cid:107) B − (cid:0) Θ − Σ (cid:1) ξ + (cid:0) (cid:107) A − Φ (cid:1) (cid:107) Ω , (C5)˙Φ = (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) (2 (cid:107) A + Φ) − ξ (cid:107) Ω , (C6) (cid:107) ˙ E = 2 ξ (cid:107) B − (cid:0) Θ − Σ (cid:1) (cid:107) E − (cid:107) j, (C7) (cid:107) ˙ B = − ξ (cid:107) E − (cid:0) Θ − Σ (cid:1) (cid:107) B, (C8) (cid:107) ˙ B = (cid:0) Σ − Θ (cid:1) (cid:107) B − ξ (cid:107) E − ξ (cid:0) (cid:107) E + (cid:107) B (cid:1) , (C9) (cid:107) ˙ E = (cid:0) Σ − Θ (cid:1) (cid:107) E + 3 ξ (cid:107) B − ( ρ + p )Σ + (cid:0) (cid:107) Σ − (cid:107) Θ (cid:1) ( (cid:107) E + (cid:107) B ) − (cid:107) E (cid:107) j, (C10)( ρ + p ) u a ∇ a (cid:107) A = − v s (cid:107) µu c n b ∇ b n c − u c n b ` Z cb , (C11) u a ∇ a (cid:107) µ + ( ρ + p ) n a ∇ a (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) = n c u b ` Z cb , (C12) ( ρ + p ) v s u a ∇ a (cid:107) Σ + ( ρ + p ) n c ∇ c (cid:107) µ = ( ρ + p ) v s (cid:16) ( v s (cid:107) µ ∇ c n c + h cb ` Z cb ) − n b n c ` Z cb (cid:17) , (C13)( ρ + p ) v s u a ∇ a (cid:107) Θ + ρ + p n c ∇ c (cid:107) µ = − ρ + p v s (cid:16) v s (cid:107) µ ∇ c n c + h cb ` Z cb (cid:17) . (C14)In the above expressions notice that ˙ u a = (cid:107) A n a , ˙ n a = (cid:107) A u a . n c n b ` Z cb = − ( ρ + p )( (cid:107) A + 2 (cid:107) Ω (cid:107) ξ ) + ( ρ + p ) (cid:2) (3 p − ρ ) − ( p − ρ − (cid:107) E − (cid:107) B ) − (cid:107) E (cid:3) + (cid:107) µ (cid:107) A (1 + v s )+( ρ + p ) (cid:0) Θ + Σ (cid:1) (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1) − (cid:15) (cid:107) µ ˆ v s + (cid:107) ˆ E(cid:37) C + (cid:107) E (cid:0) ˆ ρ C + (cid:107) j (cid:0) (cid:107) Θ + (cid:107) Σ (cid:1)(cid:1) ,h cb ` Z cb = − (cid:107) A ( ρ + p ) + ( ρ + p )( ρ + 3 p + (cid:107) E + (cid:107) B ) + (cid:107) µ (cid:107) A (1 + v s ) + ( ρ + p ) (cid:16) Θ + (cid:0) Θ + Σ (cid:1) + 2( Θ − Σ) (cid:17) − (cid:15) (cid:107) µ ˆ v s + (cid:107) ˆ E(cid:37) C + (cid:107) E (ˆ ρ C + (cid:37) C Φ + (cid:107) j Θ) ,u c n b ` Z cb = − (cid:107) A ( ρ + p ) (cid:0) (cid:107) Σ + Θ (cid:1) + (cid:107) A ( ρ + p ) (cid:0) Θ + Σ (cid:1) − (cid:15) (cid:107) µ ˙ v s + (cid:107) ˙ E(cid:37) C + (cid:107) E ˙ ρ C + (cid:107) E (cid:107) A (cid:107) j,n c u b ` Z cb = − ( ρ + p ) (cid:0) ΦΣ + (cid:107) A (cid:0) Θ + Σ (cid:1) + 2 (cid:107) Ω (cid:107) ξ (cid:1) − (cid:107) µ (1 + v s ) (cid:0) Θ + Σ (cid:1) + ˙ v s (cid:107) µ − (cid:107) ˆ E (cid:107) j − (cid:107) E (cid:16) (cid:107) ˆ j + (cid:37) C (cid:0) (cid:107) Σ + Θ (cid:1)(cid:17) . To write down the above explicit expressions for ( n c n b ` Z cb ), ( n c u b ` Z cb ) and ( h cb ` Z cb ), ( u c n b ` Z cb ) we have used, again,the evolution equation (C7) and the constraint equations (27) for (cid:107) E . With regards to the term ˆ v s in ( n c n b ` Z cb )and ( h cb ` Z cb ) and ˙ v s in ( u c n b ` Z cb ) and ( n c u b ` Z cb ), involving the derivatives of v s we note the following: the definition v s = dp/dρ has been used by setting an appropriate equation of state p = p ( ρ ). For example, in te case of apolytropic equation of state one has that v s = and therefore set aside the terms ˆ v s and ˙ v s . Otherwise, one canconsider a generic equation of state ρ = ρ ( p ) such that dv s /dρ = d p/dρ (cid:54) = 0 and ∇ a v s = (cid:0) d p/dρ (cid:1) ∇ a ρ and assume (cid:0) d p/dρ (cid:1) = constant. In the more general case we refer to the discussion in[11] —see Section 14.4. In fact, asdiscussed in Section IV, for the more general case of a fluid which is not isentropic, i. e. s (cid:54) =constant, the equation ofstate should be written in the form p = p ( ρ, s ) —see [10] and references therein. This means that when consideringderivatives of the pressure p , it should have been taken into account that ∇ a p = v s ∇ a ρ + ( ∂p/∂s ) ∇ a s , where the timeevolution of the entropy is governed by equation (C2).The evolution equations under consideration contains terms involving ∇ a v s and ∇ a s . The evolution equation for thenew variable s a ≡ ∇ a s can be obtained by covariant differentiating equation (C2), commuting covariant derivativesand using, again, the evolution equations —see [10].Concerning the terms ˙ ρ C in ( u c n b ` Z cb ) and ˆ ρ C in ( n c n b ` Z cb ) and ( h cb ` Z cb ), one can assume, for example, that thecharge density (cid:37) C is a function of the matter density (cid:37) C = (cid:37) C ( ρ ). Here we assume it to be a constant multiple of thefluid density, ρ = (cid:37)(cid:37) C . In this case one can use again equations (C1) and ( ?? ) and the definition of v s . Note that thedensity (cid:37) C appears in the constraint equations for the electric field and the pressure. The discussion of Sections (V)and (V A) shows that (cid:37) C is only involved in the evolution equation for the radial acceleration. Here (cid:37) C is inheritedfrom the term (cid:107) µ containing information from the propagation of the matter density and pressure. Finally, we canuse equation (27) and (56) for (cid:107) ˆ j in ( n c u b ` Z cb ), assuming (cid:107) j = (cid:107) j ( (cid:107) E ). In particular, we take (cid:107) j = (cid:107) Eσ J .This system (C1)-(C13) can be written matricially as A a ( a )( b ) ∂ a v ( b ) = B ( a )( b ) v ( b ) . It is convenient to write v ( b ) =( v ( i ) , v ( A ) ) where v ( i ) = ( (cid:107) E, (cid:107) B, (cid:107) B , (cid:107) E , s, n, ρ, (cid:107) Ω , ξ, Φ , (cid:107) A ) , v ( A ) = ( (cid:107) µ, (cid:107) Θ , (cid:107) Σ) . Thus it follows that A a ( i )( i ) = u a , A a ( i )( A ) = A a ( A )( i ) = 0 ,A a ( i )( j ) = A a ( j )( i ) = 0 i (cid:54) = j, while for ( (cid:107) µ, (cid:107) Θ , (cid:107) Σ): A a ( µ )( µ ) = u a , A a (Θ)(Θ) = ( ρ + p ) v s u a , A a (Σ)(Σ) = ρ + p ) v s u a ,A a ( µ )(Θ) = ( ρ + p ) n a , A a (Θ)( µ ) = ( ρ + p ) n a , A a (Σ)( µ ) = ( ρ + p ) n a ,A a ( µ )(Σ) = ( ρ + p ) n a , A a (Θ)(Σ) = 0 , A a (Σ)(Θ) = 0 . (C15)8 Appendix D: Details on the subclasses of the solutions
Subclasses of the I -class ( A , T ) In what follows we analyse in further detail the following subcases of theconfiguration( A T ):( A T Q ) : in this case the configuration is defined by the conditions A = 0, Σ = 0 and Θ = 0. The matterdensity evolution only depends on the electric field and the current density. The evolution of the electro-magnetic fields is regulated by the twist ξ of the 2-sheet. Moreover, the radial vorticity is constant duringthe motion —i.e. ˙Ω = 0. From the evolution equations for T and Q we obtain, respectively, the conditions C ( A , T ) and C ( A , Q ). From these relations we find E = p + ρ − and 4Ω = B + E + 3 p + ρ . Thiscase has at least one zero eigenvalue associated with the vorticity. This eigenvalue has multiplicity 2. Thesum of the eigenvalues of the associated 7-rank matrix is Tr˚ B = − σ J < A T Φ ) : in this case one has ˙Ω = 0, and C ( A , T ). The conditions A = 0 and Φ = 0 imply from equations(59f), (59g)) either (i) Ω = 0 or (ii) ξ = 0 and B = 0. We consider these two subcases separately:(i) For ( A T Φ Ω ) the system is characterised by the condition C ( A , T, Ω). From the criterion on the traceone can deduce that the system is unstable if ˚Θ ≤ − ( σ J ) / (3 + v s ). Now, considering the determinantof the reduced 7 × stability is that det ˚ B < p = ˚ ρv s we obtain for stable configuration the following conditions:(ia) for ˚ ρ < ˚ ρ s then one has that ˚Σ ∈ ( − ˚Σ s , +˚Σ s ) and Q ( ˚ B , ˚ ξ, ˚ E ) < or ˚Σ ∈ ( − ˚Σ x , − ˚Σ s ) ∪ (˚Σ s , ∞ )and Q ( ˚ B , ˚ ξ, ˚ E ) > ρ ≥ ρ s then one has ˚Σ ∈ ( − ˚Σ x , +˚Σ s ) and Q ( ˚ B , ˚ ξ, ˚ E ) < or ˚Σ ∈ (+˚Σ s , ∞ ) and Q ( ˚ B , ˚ ξ, ˚ E ) >
0. With ˚ ξ (cid:54) = 0, Q ( B , ξ, E ) ≡ B ξ + 2 E σ J and considering the ( LEV ) Σ We note that the reference matter density and shear ρ s and Σ x , depend on the constants ( v s , σ J ) andthe limit Σ s depends only on the matter density and the square of the sound velocity.(ii) ( A T Φ , ξ, B ): in this case the Maxwell equations simply became ˙ B = 0 and ˙ E = − j . Moreover ˙Ω = 0and the condition C ( A , T ) holds. The corresponding reduced system has six unknowns. We notice herethe eigenvalues: λ = 0 with a subspace of dimension 3 and λ = − σ J , λ ± = (cid:0) − (3 + v s )˚Θ ± (cid:113) ρ ( v s + 1) + ( v s − (˚Θ) (cid:1) . The conditions on the sign of λ ± imposes severe restrictions on the radial expansion. These depend onthe sound velocity and the background density matter: the condition λ ± < stability )is satisfied if ˚Θ > (cid:112) (1 + v s )˚ ρ/ A T ξ ) : this case implies the equations ˙ B = ˙ B = ˙Ω = ˙Φ = 0 while the electric field satisfies ˙ E = − j . Moreover,one has the condition C ( A , ξ ) with C ( A , T ). We study the associated 8 × λ ± and λ i as defined for the case ( A T Φ , ξ, B ) and a zero eigenvalue with multiplicity 5 due tothe variables ( B, B , Ω , Φ).( A T Ω ) : in this case one has ˙Φ = 0. From equation (60a) one finds C ( A , T, Ω). The time evolution of theelectric and magnetic fields are regulated by the twisting ξ of the 2-sheet —see equations (59b)-(59c)). Thetemporal evolution of ξ is fixed explicitly by the magnetic part of the Weyl tensor —see equation (59g).On the other hand, from the trace of the associated rank 7 matrix we infer that the system is unstable if˚Σ < − Σ x . In order to obtain necessary conditions for stability one requires c = − Tr˚ B >
0, and c > LEV ) Θ :(i) if ˚ ρ < ˚ ρ s then ˚Θ ∈ ( − ˚Θ x , − ˚Θ s ) ∪ (+˚Θ s , ∞ ) and ˚ ξ (cid:54) = 0, or ˚Θ ∈ ( − ˚Θ x , +˚Θ s ) and ˚ ξ ˚ B < ρ ≥ ρ s then ˚Θ ∈ ( − ˚Θ s , +˚Θ s ) and ˚ ξ ˚ B < or ˚Θ ∈ (+˚Θ s , ∞ ) and ˚ ξ (cid:54) = 0 Subclasses of the II -class ( A , Q ) In this subsection we focus on the configurations with Θ = − Σ. Taking intoaccount the results on the system ( A Q ), we consider the following subcases.( A Q Φ ) : the condition C ( A , Q ) holds. Moreover, the conditions A = 0 and Φ = 0, imply from equation (59f)two subcases, (i) Ω = 0 and (ii) ( ξ = 0 , B = 0), respectively.(i) ( A Q Φ Ω ): from the trace of the reduced rank 7 matrix the configuration is bound to be unstable if˚Σ ≥ Σ s ≡ σ J / v s ). On the other hand, studying the sign of the characteristic polynomial9coefficient c < stability : for ˚ E < ˚ E p and (a) ˚ T > − σ J / (15 + 2 v s ) with˚ ξ < − ˚ ξ p or ˚ ξ > ˚ ξ p , or (b) ˚ E < ˚ E p and ˚ T > ˚ T p . Otherwise, for ˚ E ≥ ˚ E p with ˚ T > − σ J / (15 + 2 v s ), withthe following definitions:˚ E p ≡ (cid:0) + v s + v s (cid:1) ˚ ρ + 2 (cid:0)
50 + v s (15 + 2 v s ) (cid:1) σ J (15 + 2 v s ) , ˚ T p ≡ − v s ) σ J + (cid:115) v s ) (cid:18)(cid:0) v s (5 + 2 v s ) (cid:1) ˚ ρ − E (cid:19) + 9(13 + 2 v s ) σ J
285 + 78 v s , ˚ ξ p = √ (cid:113) (cid:0) v s (5 + 2 v s ) (cid:1) ˚ ρ − E − T (95 + 26 v s ) − T (13 + 2 v s ) σ J , with T = − σ J , v s ). This stability conditions can be alternatively expressed as follows: ρ > , ˚Σ > Σ s , and ˚ ξ > ˚ ξ s ≡ (cid:18)
2( ˚ B + ˚ E )+4 (cid:0) v s (3+ v s ) (cid:1) ˚ ρ +3˚Σ (cid:0) σ J (13+2 v s ) − v s )˚Σ (cid:1)(cid:19) . (ii) ( A Q Φ ξ B ): with the condition C ( A , Q ), while the time evolution of the radial projected vorticityis entirely regulated by the radial part of the shear —see equation (59h). The system, with a rank6 matrix, is unstable if ˚Σ ≥ Σ s ≡ σ J / v s ). Otherwise, it is stable if (a) ˚ E < ˚ E p with˚ T > − (2 σ J ) / (13 + 2 v s ) and ˚Ω > ˚Ω p ∪ ˚Ω < − ˚Ω p or − ˚Ω p ≤ ˚Ω ≤ ˚Ω p with ˚ T > ˚ T p ; (b) ˚ E ≥ ˚ E p and˚ T > − (2 σ J ) / (13 + 2 v s ), where˚Ω p ≡ √ (cid:115) − E + 2 (cid:0) v s (2 + v s ) (cid:1) ˚ ρ + 6 (73 + 26 v s + 4 v s ) σ J (13 + 2 v s ) , ˚ T p ≡ − v s ) σ J
210 + 66 v s ++ √ (cid:113) v s ) (cid:0) v s (2 + v s ) (cid:1) ˚ ρ + 3(11 + 2 v s ) σ J − E (35 + 11 v s ) − v s )˚Ω
210 + 66 v s , ˚ E p ≡ (cid:0) v s (2 + v s ) (cid:1) ˚ ρ + 2 (cid:0)
73 + 26 v s + 4 v s (cid:1) σ J (13 + 2 v s ) . Alternatively, these conditions can be reexpressed as˚Σ < Σ s and ˚Ω > (cid:18) ˚ B (3 ˚ E − − E + (cid:0) v s (9+4 v s ) (cid:1) ˚ ρ +6˚Σ (cid:0) (11+2 v s ) σ J − v s )˚Σ (cid:1)(cid:19) / . ( A Q ξ ) : in this case one has Σ = − Θ / C ( A , ξ ) from equation (59g) and C ( A Q ). The radial shear is theonly kinematical variable that explicitly regulates the time evolution of the variables ( E, B, ρ, Φ , E , B ). Theradial vorticity and the shear are related by the two evolution equations (59h) and (60a)), respectively.This is a rank 8 matrix problem. The system is unstable if ˚Σ ≥ (2 σ J ) / (51 + 6 v s ). The system is stable ifeither (a) ˚ E < ˚ E p and ˚ T > ˚ T p or − σ J / (17 + 2 v s ) < ˚ T < ˚ T p and ˚Ω < − ˚Ω p ∪ ˚Ω > ˚Ω p ; or (b) ˚ E ≥ ˚ E p and˚ T > − σ J / (17 + 2 v s ). In the above conditions we used the following definitions:˚Ω p ≡ √ (cid:113) (cid:0) v s (2 + v s ) (cid:1) ˚ ρ − E − T (25 + 6 v s ) − T (15 + 2 v s ) σ J , ˚ E p ≡ (cid:0) v s (2 + v s ) (cid:1) ˚ ρ + 4 (cid:0)
65 + v s (17 + 2 v s ) (cid:1) σ J (17 + 2 v s ) , ˚ T p ≡ (cid:115) v s ) (cid:18) (cid:0) v s (2 + v s ) (cid:1) ˚ ρ − v s ) σ J − E (cid:19) + 9(15 + 2 v s ) σ J
375 + 90 v s . < σ J / (51 + 6 v s ) , and ˚Ω > (cid:18) ˚ B + ˚ E + (cid:0) v s (9+4 v s ) (cid:1) ˚ ρ +3˚Σ (cid:0) (30+4 v s ) σ J − v s )˚Σ (cid:1)(cid:19) . ( A Q Ω ) : in this case one has that T = −
3Σ and Σ = − Θ / C ( A , Q ) : B + E + 2 E = − ( p + ρ/
3) holds. The system is unstable if ˚Σ ≥ Σ s ≡ σ J / v s ). A stable system meets the followingconditions: (a) ˚ E < ˚ E p , with ˚ T > ˚ T p or − s < ˚ T ≤ ˚ T p and (˚ ξ < − ˚ ξ p , ˚ ξ > ˚ ξ p ); or (b) ˚ E ≥ ˚ E p and˚ T > − s . In the above expressions we have used the definitions˚ T p ≡ − v s ) σ J + (cid:113) v s ) (cid:0) (5 + 3 v s (5 + 2 v s ))˚ ρ − E (cid:1) + 9(7 + v s ) σ J
165 + 42 v s , ˚ ξ p ≡ (cid:113) (5 + 3 v s (5 + 2 v s ))˚ ρ − E − T (55 + 14 v s ) − T (7 + v s ) σ J √ , ˚ E p ≡ (8 + v s ) (5 + 3 v s (5 + 2 v s ))˚ ρ + 3(57 + 2 v s (8 + v s )) σ J v s ) . The conditions can be expressed, alternatively, as˚Σ < Σ s and 16˚ ξ > (cid:0) ˚ B + ˚ E + 2(1 + v s (3 + v s ))˚ ρ − v s )˚Σ + 6(7 + v s )˚Σ σ J (cid:1) . Subclasses of the IV -class ( A , ξ ) ( A ξ Ω ): the assumptions A = 0, ξ = 0 and Ω = 0, lead from equation (59g) tothe condition C ( A ξ, Ω) : B = 0. Thus, we analyse the case: ( A ξ Ω B ), with a rank 7 matrix. The system is unstable if ˚Θ ≤ − σ J + ( v s − / v s ) . The case A = 0 , T = 0 , and Q = 0 ( A T Q Φ ): in this special case we assume the radial acceleration to be zero. Inaddition, the expansion of the 2-sheet Φ vanishes and also Σ = Θ = 0. It is worth noting that the assumption A = Σ = Θ = Φ = 0 implies from equation (59f) that ξ Ω = 0 —that is, the twist of the 2-sheet ξ togetherwith the magnetic part of the Weyl tensor vanish. This system can only accelerate in the radial direction orotherwise radially projected vorticity will vanish. The basic assumptions in this case directly lead to ˙Ω = 0,˙ ξ = B / C ( A , T ) and C ( A , Q ) hold. Thus, two subcases occur:(i) ( A T Q Φ ξ B ); in this case one further has ˙ B = ˙Ω = 0, ˙ E = − j . The only non zero eigenvalue is λ = − σ J .(ii) ( A T Q Φ Ω ); in this case, it can be shown that the trace of the reduced matrix Tr˚ B = − σ J < C ( A , T, Ω) and C ( A , Q ) is B + E + 2 E = − (cid:0) p + ρ (cid:1) < C ( A , T, Ω)and C ( A , Q ) must be satisfied simultaneously.( A T Q ξ ) : in this case the problem is simplifies considerably and we find the relations C ( A , ξ ), C ( A , T ) and C ( A , Q ). Moreover, one has that ˙ E = − j . The evolution equations ˙ B = 0, ˙ B = 0, ˙Ω = 0, and ˙Φ = 0 giverise to repeated zero eigenvalues. In addition one has the eigenvalue λ = − σ J < A T Q Ω ) : in this case one readily has that ˙Φ = 0. The associated rank 5 matrix has a zero eigenvalue withmultiplicity 2. In addition, one has that c = Tr˚ B = − σ J <
0. Thus, imposing the condition c c > ξ > − ( ˚ B + ˚ E + 3˚ E ) /
26. However, in this case the relations C ( A , T, Ω) and C ( A , Q ) hold. Thesecannot be satisfied for ordinary matter density (i.e. such that ρ >
0) and the given equation of state.
The case A = 0 , T = 0 , ξ = 0 ( A T ξ Ω ): the assumptions A = 0, ξ = 0 and Ω = 0 lead using equation (59g) to B = 0. From equation (60a) one obtains the conditions C ( A , T, Ω), and ˙Φ = ˙ B = 0. Thus, we analyse thesystem ( A T ξ Ω B ). There is a zero eigenvalue with multiplicity 3. In addition, one has λ = − σ J < λ ± = ( − (3 + v s ) ˚3Σ ± (cid:113) v s ) ˚ ρ + ( v s − ( ˚3Σ) ), the condition λ ± < > v s ) (cid:113) ˚ ρ/ (3( v s + 5 + 2 √ v s + 5 − √ . The case A = 0 , Q = 0 , ξ = 0 ( A Q ξ Ω ): this case reduces to ( A Q ξ Ω B ) with C ( A , Q ). From the associated rank6 matrix trace we infer that the system is certainly unstable if ˚Σ ≥ Σ s ≡ σ J / v s ). On the other hand, ifthe system is stable , then the following conditions must be verified: (a) ˚ E > ˚ E p with ˚ T > − s , or (b) ˚ E ≤ ˚ E p and ˚ T > ˚ T p where we introduced the following notation˚ E p ≡ v s ) (2 + 3 v s (2 + v s ))˚ ρ + 3(61 + 4 v s (6 + v s )) σ J v s ) , ˚ T p ≡ − v s ) σ J + √ (cid:113) (59 + 20 v s ) (cid:0) v s (2 + v s ))˚ ρ − E (cid:1) + 6(5 + v s ) σ J
177 + 60 v s . This condition correspond to the requirements c > B < c = − Tr˚ B > c c >
0. It can be also written as: o Σ < Σ s and ( o E ) +( o B ) < − [3+ v s (9+4 v s )] o ρ +3 o Σ[3(59+20 v s ) o Σ − v s ) σ J ]. The case A = 0 , Φ = 0 and ξ = 0 The case ( A Φ ξ Ω ) reduces to ( A Φ ξ Ω B ). We can say that the rank 6 systemis unstable when ˚Θ ≤ − σ J ) / v s ). The case A = 0 , T = 0 , Q = 0 and ξ = 0 In this case one has ( A T Q ξ Ω ) which reduces to ( A T Q ξ Ω B ). Thissystem is characterised by ˙Φ = 0, ˙ B = 0, ˙ E = − j . There are the eigenvalues λ = 0 with multiplicity 4 and λ = − σ J <
0. Nevertheless, one has that C ( A , T, Ω) and C ( A , Q . These conditions cannot be satisfied by theassumptions ρ > [1] Z. B. Etienne, Y. T. Liu, & S. L. Shapiro, Relativistic magnetohydrodynamics in dynamical spacetimes: a new AMRimplementation , Phys. Rev. D , 084031 (2010).[2] J. A. Font, Numerical hydrodynamics in general relativity , Living Rev. Rel. (2003).[3] J. A. Font, Numerical hydrodynamics and magnetohydrodynamics in general relativity , Living Rev. Rel. , 7 (2007) .[4] M. Alcubierre, Introduction to numerical Relativity , Oxford University Press, 2008.[5] B. Giacommazo & L. Rezzolla,
WhiskeyMHD: a new numerical code for general relativistic MHD , Class. Quantum Grav. , S235 (2007).[6] A. Lichnerowicz, Relativistic hydrodynamics and magnetohydrodynamics , Benjamin, New York, 1967.[7] L. Anton, O.Zanotti, J.A.Miralles, J. M. Marti, J. M. Ibanez, J. A. Font and J. A. Pons,
Numerical 3+1 general relativisticmagnetohydrodynamics: A Local characteristic approach , Astrophys. J. , 296 (2006).[8] D.Radice, L.Rezzolla and F.Galeazzi,
High-Order Fully General-Relativistic Hydrodynamics: new Approaches and Tests,
Class. Quant. Grav. ,075012 (2014).[9] H.Witek, Numerical Relativity in higher-dimensional space-times,
Int. J. Mod. Phys. A , 1340017 (2013).[10] D. Pugliese and J.A. Valiente Kroon, On the evolution equations for ideal magnetohydrodynamics in curved spacetime ,Gen. Rel. Grav. , 2785 (2012).[11] Y. Choquet-Bruhat, General Relativity and the Einstein equations , Oxford University Press, 2008.[12] M. Shibata & Y. Sekiguchi,
Magnetohydrodynamics in full general relativity: formulations and tests , Phys. Rev. D ,044014 (2005).[13] T. W. Baumgarte & S. L. Shapiro, General relativistic magnetohydrodynamics for the numerical construction of dynamicalspacetimes , Astrophys. J. , 921 (2003).[14] C. Palenzuela, D. Garrett, L. Lehner, & S. Liebling,
Magnetospheres of black hole systems in force-free plasma , Phys. Rev.D , 044045 (2010).[15] A.M. Anile Relativistic fluids and magneto-fluids: With applications in astrophysics and plasma physics , CambridgeUniversityPress, Cambridge, U.K.; New York, U.S.A., 1989.[16] M. M. Disconzi,
On the well-posedness of relativistic viscous fluids,
Nonlinearity , 1915 (2014).[17] E. Horst, Symmetric Plasmas and Their Decay , Commun. Math. Phys. , 613-633 (1990).[18] S. C. Hsu, T. J. Awe, S. Brockington, A. Case, J. T. Cassibry, G. Kagan, S. J. Messer, M. Stanic, X. Tang, D. R.Welch, and F. D. Witherspoon,
Spherically Imploding Plasma Liners as a Standoff Driver for Magnetoinertial Fusion , Ieeetransactions on plasma science , 5 (2012).[19] Tai-Ho Tan and Joseph E. Borovsky, Spherically symmetric high-velocity plasma expansions into background gases , Journalof Plasma Physics , 02 239 (1986).[20] P. D. Laskyand A. W. C. Lun Gravitational collapse of spherically symmetric plasmas in Einstein-Maxwell spacetimes .Phys. Rev. D , 104010 (2007).[21] R. L. Viana, R. A. Clemente and S. R. Lopes, Spherically symmetric stationary MHD equilibria with azimuthal rotation ,Plasma Phys. Control. Fusion , 197 (1997). [22] Yan Guo and A. Shadi Tahvildar-Zadeh, Formation of singularities in relativistic fluid dynamics and in spherically sym-metric plasma dynamics , Contem. Math., , 151-161 (1999).[23] O. Reula,
Hyperbolic methods for Einstein’s equations , Living Rev. Rel. , 1 (1998).[24] O. Reula, Exponential decay for small nonlinear perturbations of expanding flat homogeneous cosmologies , Phys. Rev. D , 083507 (1999).[25] C. Clarkson, A Covariant approach for perturbations of rotationally symmetric spacetimes , Phys. Rev. D , 104034 (2007).[26] G. F. R. Ellis & H. van Elst, Cosmological models: Cargese lectures 1998 , NATO Adv. Study Inst. Ser. C. Math. Phys.Sci. , 1 (1998).[27] C. A. Clarkson, M. Marklund, G. Betschart, and P. K. S. Dunsby,
The Electromagnetic Signature of Black Hole Ring-Down , Astrophys. J. , 492 (2004).[28] M. Marklund and C. Clarkson,
The General relativistic MHD dynamo , Mon. Not. Roy. Astron. Soc. (2005) 892.[29] G. Betschart, and C. A. Clarkson,
Scalar field and electromagnetic perturbations on locally rotationally symmetric space-times , Class. Quantum Grav. , 5587 (2004).[30] R. B. Burston, , Class. Quantum Grav. Covariant perturbations of Schwarzschild black holes , Class. Quantum Grav. The covariant approach to LRS perfect fluid spacetime geometries ,Class. Quantum Grav. , 1099 (1996).[33] R. B. Burston and A. W. C. Lun, , Class. Quant. Grav. (2008) 075003.[34] R. B. Burston, , Class. Quant. Grav. , 993-1000 (2006).[36] O. Zanotti and D. Pugliese, Von Zeipel’s theorem for a magnetized circular flow around a compact object , Gen. Rel. Grav. (2015) 4, 44.[37] Bekenstein, J.D., Oron, E.: New conservation laws in general-relativistic magnetohydrodynamics. Phys. Rev. D 18,1809171819 (1978). DOI 10.1103/PhysRevD.18.1809[38] Bekenstein, J.D., Oron, E.: Interior magnetohydrodynamic structure of a rotating relativistic star. Phys. Rev. D 19,2827172837 (1979). DOI 10.1103/PhysRevD.19.2827[39] H. Stephani, D. Kramer, M. MacCallum, C. Hoenselaers, E. Herlt , Exact Solutions of Einstein’s Field Equations , CambridgeMonographs on Mathematical Physics, Paperback 2009.[40] C. Lubbe and J. A. Valiente Kroon,
A conformal approach for the analysis of the non-linear stability of pure radiationcosmologies , Annals Phys. , 1 (2013).[41] J.M. Stewart, and M. Walker, Proc. R. Soc. London
A 431 , 49 (1974).[42] A. Alho, F. C. Mena, & J. A. Valiente Kroon,
The Einstein-Friedrich-nonlinear scalar field system and the stability ofscalar field Cosmologies , In arXiv:1006.3778 , (2010).[43] Q. I. Rahman and G. Schmeisser,
Analytic Theory of Polynomials: Critical Points, Zeros and Extremal Properties . LondonMathematical Society Monographs, Clarendon Press, 2002.[44] S. Barnett,
A New Formulation of the Theorems of Hurwitz, Routh and Sturm
J. Inst. Maths Applies
240 (1971).[45] H. O. Kreiss and J. Lorenz,
Stability for time-dependent differential equations . Acta Numerica , 203 (1998).[46] H. O. Kreiss, O. E. Ortiz and O. A. Reula, Stability of quasi-linear hyperbolic dissipative systems . Journal of DifferentialEquations , 78 (1998).[47] H. O. Kreiss, G. B. Nagy, O. E. Ortiz and O. A. Reula,
Global existence and exponential decay for hyperbolic dissipativerelativistic fluid theories . Journal of Mathematical Physics , 5272 (1997).[48] O. E. Ortiz, Stability of nonconservative hyperbolic systems and relativistic dissipative fluids . Journal of MathematicalPhysics , 1426 (2001).[49] We refer to [20, 25, 27, 29, 30, 33, 34] for a general discussion concerning the configurations with a unique symmetrydirection, described by the Einstein-Maxwell-Euler system. We mention also the well known solution of toroidal magneticfield widely adopted in the axes-symmetric accretion configurations [35], and [36] for a discussion on the case of a poloidalmagnetic field where a metric representation is adapted to the direction of the field, as proved by Bekenstein & Oron[37, 38]. Moreover for a deeper and more general discussion about the MHD configurations in spherical symmetry, see forexample the work [17–22][50] As pointed out in [25] in a locally rotationally symmetric spacetime any background quantities are scalars, implying thatthe vector and tensor quantities are automatically gauge invariant, under linear perturbations as a consequence of theStewart-Walker lemma [41].[51] The non-zero coefficients of the characteristic polynomial are: c = (6˚ ξ ˚Ω) > c = ˚ ξ (8 σ J ˚ E + 18˚ B ˚ ξ + 9 σ J ˚Ω ), c = ˚ ξ (cid:0) σ J (3˚ B − B ˚ E ) + 2˚ ξ (cid:0) −
4( ˚ B + ˚ E − ξ ) (cid:1)(cid:1) , c = ˚ ξ (˚ B + 2 σ J ˚ ξ ), c = 13˚ ξ > c = − Tr ˚ B > . [52] det ˚ B = (1 + v s )˚ ξ (2 σ J ˚ E + 3˚ B ˚ ξ ) (cid:0) (1 + 3 v s )˚ ρ − Q (cid:1)(cid:1)