Causality and stability in relativistic viscous magneto-fluid dynamics

We investigate the causality and the stability of the relativistic viscous magneto-hydrodynamics in the framework of the Israel-Stewart (IS) second-order theory, and also within a modified IS theory which incorporates the effect of magnetic fields in the relaxation equations of the viscous stress. We compute the dispersion relation by perturbing the fluid variables around their equilibrium values. In the ideal magnetohydrodynamics limit, the linear dispersion relation yields the well-known propagating modes: the Alfv\'en and the magneto-sonic modes. In the presence of bulk viscous pressure, the causality bound is found to be independent of the magnitude of the magnetic field. The same bound also remains true, when we take the full non-linear form of the equation using the method of characteristics. In the presence of shear viscous pressure, the causality bound is independent of the magnitude of the magnetic field for the two magneto-sonic modes. The causality bound for the shear-Alfv\'en modes, however, depends both on the magnitude and the direction of the propagation. For modified IS theory in the presence of shear viscosity, new non-hydrodynamic modes emerge but the asymptotic causality condition is the same as that of IS. In summary, although the magnetic field does influence the wave propagation in the fluid, the study of the stability and asymptotic causality conditions in the fluid rest frame shows that the fluid remains stable and causal given that they obey certain asymptotic causality condition.


Introduction
In relativistic heavy-ion collisions experiment, two fast moving charged nuclei collide with each other and generate a deconfined state of matter known as Quark-Gluon-Plasma (QGP). In non-central collisions an extremely strong magnetic field (∼ 10 18 -10 19 Gauss) is also produced in the initial stages refs. [1][2][3][4][5] mostly due to the spectator protons.
The charge separation in Au+Au collisions are claimed to be observed by the STAR collaboration refs. [63][64][65]. However, it is still a challenge to extract the CME signals from the huge backgrounds caused by the collective flows refs. [66][67][68]. Therefore, it requires the systematic and quantitative studies of the evolution of the QGP coupled with the electromagnetic fields for the discovery of CME. It is widely accepted that the QGP produced in high energy heavy-ion collisions behaves as almost ideal fluid (i.e., possess very small shear and bulk viscosity). This conclusion was made primarily based on the success of relativistic viscous hydrodynamics simulations in explaining a multitude of experimental data with a very small specific shear viscosity (η/s) as an input refs. [69][70][71][72][73][74][75][76]. Most of these theoretical studies use IS second-order causal viscous hydrodynamics formalism or some variant of it. The fact that the QGP is composed of electrically charged quarks indicates that it should have finite electrical conductivity which is corroborated by the lattice-QCD calculations refs. [77][78][79] and perturbative QCD calculations refs. [80,81]. The electrical conductivity of the QGP and the hadronic phase was also calculated by various other groups (mostly using the Boltzmann transport equation) see refs. [82][83][84][85][86][87][88][89][90][91][92][93][94][95][96]. It is then natural to expect that the appropriate equation of motion of the high temperature QGP and low temperature hadronic phase under large magnetic fields is given by the relativistic viscous magneto-hydrodynamic framework. As mentioned earlier the IS second-order theory of causal dissipative fluid dynamics, although successful, known to allow superluminal signal propagation (and hence acausal) under certain circumstances refs. [97][98][99][100]. It is then important to know under what physical conditions the theory remains causal and stable in presence of a magnetic field which is also important for the numerical MHD studies of heavy-ion collisions.
Relativistic magnetohydrodynamics (MHD) is a self-consistent macroscopic framework which describe the evolution of any charged fluid in the presence of electromagnetic fields refs. [101][102][103][104][105][106][107]. In ref. [4], we have computed the ratio of the magnetic field energy to the fluid energy density in the transverse plane of Au-Au collisions at √ s N N = 200 GeV in the event-by-event simulations. Our results imply that the magnetic field energy is not negligible. In ref. [101], we have derived the analytic solutions of a longitudinal Bjorken boost invariant MHD with transverse electromagnetic fields in the ideal limit. We have found that the transverse magnetic fields will decay as ∼ 1/τ with τ being the proper time. Later, in ref. [102], we have studied the corrections from the magnetization effects and extended the discussion to (2 + 1)-dimensional ideal MHD refs. [108,109]. We have also investigated the effects of large magnetic fields on (2 + 1)-dimensional reduced MHD at √ s N N = 200 GeV ref. [110]. Very recently, we have derived the analytic solutions of MHD in the presence of finite electric conductivity, CME and chiral anomaly ref. [106] and extended the results to cases with the transverse and longitudinal electric conductivities ref. [107]. For numerical simulations of ideal MHD, one can see refs. [104,105]. As mentioned earlier in the ordinary relativistic hydrodynamics, the widely used framework is the second order IS theory [111]. The pioneering studies on the instabilities of first order hydrodynamics are shown in refs. [97,112]. Later, the systematic studies for the dissipative fluid dynamics have been done earlier with bulk viscous pressure [99], shear viscous stress [98] and heat currents [113], also see refs. [100,114]. There have been several recent studies on casualty and stability of ideal MHD in refs. [115][116][117] and reference therein.
We aim to study the stability and causality of the IS theory for MHD, whose form is derived by the complete moment expansion as done in refs. [118,119]. First, we analyze the propagating modes in ideal non-resistive MHD. Next, we discuss the causality and stability of the relativistic MHD with dissipative effects. To analyse the causality and stability of the relativistic viscous fluid, we linearise the relevant equations by using a small sinusoidal perturbation around the local equilibrium and study the corresponding dispersion relations in line with the studies in refs. [97][98][99]112] .
The manuscript is organized as follows: in section 2 we briefly discuss the energymomentum tensor of fluid for ideal MHD case and the modified IS theory. In section 3 we revisit the standard analysis of causality and stability of a system without magnetic fields. Then, in section 4 we show the stability and causality of an ideal MHD and carry out the analysis of characteristic velocities in section 5. In section 6 we consider the newly developed IS theory for non-resistive MHD. Finally, we conclude our work in section 7.

Causal relativistic fluid in presence of magnetic field
In this work we consider the causal relativistic second order theory for relativistic fluids by Israel-Stewart (IS) and also a modified form of the IS theory in presence of a magnetic field given in ref. [118], for later use we define it as NRMHD-IS theory (here NRMHD corresponds to non resistive magneto-hydrodynamics). The total energy momentum tensor of the fluid can be written as, where ε, P are fluid energy density, pressure , u µ is the fluid four velocity and Π, π µν are bulk viscous pressure and shear viscous tensor, respectively. The magnetic and electric four vectors are defined as is the field strength tensor. The space-time evolution of the fluid and magnetic fields are described by the energy-momentum conservation coupled with Maxwell's equations (2.4) The non-resistance limit means the electric conductivity σ e is infinite. In this limit, in order to keep the charge current j µ = σ e E µ be finite, the E µ → 0. Then, the relevant Maxwell's equations which govern the evolution of magnetic fields in the fluid is, For simplicity, we will also neglect the magnetisation of the QGP, which implies an isotropic pressure and no change in the Equation of Sate (EoS) of the fluid due to magnetic field (e.g. see ref. [102]).
In the original IS theory the viscous stresses Π, π µν are considered as an independent dynamical variables given by the following equations (e.g. see refs. [120][121][122]), where ζ and η are bulk and shear viscosity, respectively. The coefficients τ Π and τ π are the relaxation times for the bulk and shear viscosity, respectively and ω µν ≡ 1 2 ∆ µα ∆ νβ (∂ α u β − ∂ β u α ) is the vorticity tensor. The subscript NS means the Navier-Stokes values, Note that all of these coefficients are functions of baryon chemical potential (µ) and temperature (T ). Equation (2.7) can be derived from the kinetic theory via complete moment expansion, one can see refs. [123][124][125] for more details.
For further simplification, we also ignore the coupling of viscosity with other dissipative forces and concentrate on the following terms, We note that in principle the magnetic field may cause viscous tensor to be anisotropic as shown in ref. [126] but in this work we consider zero magnetisation and hence use eqs. (2.10), (2.11) for simplicity.

Dispersion relation in the absence of magnetic field
As is known, IS theory is a consistent fluid dynamical prescription which preserves causality provides that the relaxation time associated with the dissipative quantities (such as shear and bulk viscous stresses) are not too small refs. [97][98][99][100][112][113][114]. Here we aim to study the stability and causality of a relativistic viscous fluid (governed by the IS equations) in an external magnetic field by linearising the governing equations under a small perturbation. Before discussing the causality and stability of a relativistic viscous fluid in a magnetic field, for the sake of completeness, let us summaries here the findings without the magnetic field. We note that the following results are not new and most of them can be found in refs. [98,99,114].

Dispersion relation for bulk viscosity
We consider a perturbation around the static quantities X 0 , where we choose five independent variables X = (ε, u x , u y , u z , Π). Here, we only consider the system in the local rest frame, i.e. u µ 0 = (1, 0). Then, we linearise eq. (2.3), (2.10) in vanishing magnetic fields and shear viscous tensor limit and obtain a cubic polynomial equation of the form given in eq. (A.2) with X i 's are and the other two roots being zero. The solutions of this cubic polynomial are obtained from eq. (A.3). Here, we introduce a constant α = c 2 s , where c s is speed of sound. We adopt the following parametrisation of the bulk viscosity coefficient and the relaxation time refs. [99,114]: where s and T are the entropy density and the temperature, respectively. The parameters a 1 and b 1 characterize the magnitudes of the viscosity and the relaxation time, respectively. In the small wave-number limit, the dispersion relation is Whereas the asymptotic forms of the dispersion relation in this case for large k are Note that one of the roots is a pure imaginary which is also known as the non-hydrodynamic mode because it is independent of k in the k → 0 limit.

Dispersion relation for shear viscosity
We use the following parametrization taken from ref. [98] for the shear viscous coefficient and the corresponding relaxation time: Again we linearise eqs. (2.3), (2.11) (the magnetic field and the bulk viscous pressure are taken to be zero) and obtain a set of equations with nine independent variables. Two of the roots are non-hydrodynamic with corresponding dispersion relation is ω = i/τ π . Another four roots are where each roots are double degenerate, they are known as the shear modes. The remaining three modes are obtained from a cubic polynomial of the form given in eq. (A.2) with X i 's are (3.10) These modes called sound modes as given in ref. [98]. In the small k limit, the dispersion relation for the sound modes are  For details see ref. [98].

Dispersion relation in the presence of magnetic field
We extend our studies to explore the cases in a non-vanishing magnetic field. In this section, we will investigate the dispersion relation and the speed of sound in a viscous fluid in the presence of a homogeneous magnetic field. We will derive the physical conditions of causality and stability. To achieve this goal, we carry out a systematic study for the following cases, (i) non-resistive ideal MHD, (ii) viscous MHD with bulk viscosity only, (iii) with shear viscosity only, (iv) with both bulk and shear viscosity.

Ideal MHD
For an ideal non-resistive fluid in magnetic field the energy-momentum tensor eq. (2.1) takes the following form which is normalized to b µ b µ = −1 and orthogonal to u µ i.e, b µ u µ = 0. Again we consider the similar perturbation as eq. (3.1) around the equilibrium configuration in the local rest frame (u µ 0 = (1, 0)). Ignoring the second and higher-order terms for the perturbations in ε, P, u µ and B µ , the perturbed energy-momentum tensor can be expressed as, Next, using the above δT µν in the energy-momentum conservation equation and noting that ∂ µ δT µν = 0 we get the following four equations, Here, we define h = ε 0 + P 0 + B 2 0 , and use δP = αδε. The relevant Maxwell's equations which govern the evolution of magnetic fields in the fluid is µναβ ∂ β F να = 0, which can also be written in the following form (4.8) Linearizing the above Maxwell's equations lead to the following set of equations, The equations of motion are the energy-momentum conservation equations [Eqs. (4.4)-(4.7)] and the Maxwell's equations [eq. (4.9)-(4.12)]. However, we notice that eq. (4.9) does not include a time-derivative and it is a constraint equation for δB, δb x and δb y . This constraint is consistently propagated to the remaining system of equations of motion. After replacing δB by δb x and δb y , these equations become, where, δX = δε, δũ x , δũ y , δũ z , δb x , δb y , (4.14) and A is a 6 × 6 matrix of the following form, In deriving the above equations, we have also used the following condition δũ µ b µ +u µ δb µ = 0, for changing the variable from δb t to δũ z . Without loss of generality, we consider the magnetic magnetic field b µ along the z-axis and k µ lies in the x-z plane and making an angle θ with the magnetic field, i.e., (4.16) The dispersion relations are obtained by solving which gives us six hydrodynamic modes. Two of these modes are the called Alfvén modes whose dispersion relation are given as, where, v A is the speed of Alfvén wave. The fluid displacement is perpendicular to the background magnetic field in this case and the Alfvén modes can be thought of as the usual vibrational modes that travel down a stretched string. The rest four modes correspond to the magneto-sonic modes with the following dispersion relations where v M is the speed of the magneto-sonic waves, The ± sign before the square-root term is for the "fast" and the "slow" magneto-sonic waves, respectively. From eq. (4.20), it is clear that when the propagation of the perturbation is parallel to the background magnetic fields (θ = 0), then the fast magneto-sonic modes propagate with the speed of sound (c s = √ α) and the slow magneto-sonic modes propagates with the speed of the Alfvén waves (v A ). Whereas for θ = π 2 , the velocity of the slow magneto-sonic mode becomes zero and the velocity of the fast magneto-sonic wave is More discussions can be found in refs. [116,117].

MHD with bulk viscosity
Next, we consider QGP with finite bulk viscosity and a non-zero magnetic field. Usually, the bulk viscosity is proportional to the interaction measure (ε − 3P )/T 4 of the system and hence supposed to be zero for a conformal fluid. Lattice calculation as in refs. [127,128] shows that the interaction measure has a peak around the QGP to hadronic phase cross-over temperature T co . For the sake of simplicity, here we take ζ/s = constant in the following calculation. The energy-momentum tensor in this case takes the following form, As before, we can decompose the energy-momentum tensor into two parts: an equilibrium and a perturbation around the equilibrium i.e, Here, the perturbed energy-momentum tensor takes the following form, We choose the independent variables as, These conservation equations can be cast into the form AδX T = 0 and setting det A = 0, we get, The solution of eq. (4.26) gives the following dispersion relation, These two solutions of eq. (4.29) correspond to the Alfvén modes where v A is the Alfvén velocity. The rest five modes obtained from eq. (4.27) correspond to the magneto-sonic modes. Generally, quintic equations cannot be solved algebraically. Fortunately, we find solutions for some special cases discussed below. For θ = 0, we find that two modes coincides with the Alfvén modes in eq. (4.29) and the remaining three modes are obtained from a third-order polynomial of the form given in eq. (A.2), with the coefficients X 0 , X 1 , X 2 given as The solutions of this cubic polynomial can be written as where l = 1, 2, 3, ξ is the primitive cubic root of unity, i.e., ξ = −1+ For θ = π/2, the eq. (4.27) reduces to a third-order polynomial of the form eq. (A.2), where X i 's are given as, where v f is the group velocity for the fast magneto-sonic waves defined in eq.(4.21) and the other two roots are zero. Note that all three roots in eq. (4.31) are complex because the coefficients of eq. (4.30) are complex and hence the phase velocity of any perturbations may contains a damping or growing and an oscillatory component. The left panel of Fig. 1 shows the imaginary part of the normalised ω as a function of the k/T and the right panel shows the group velocity as a function of k/T for different values of magnetic fields. Note that the imaginary part of the non-propagating mode increases and imaginary part of the propagating modes decreases when the magnetic field increases. But it is clear that (ω) always lies in the upper half of the complex plane for the parameters considered here. This implies that any perturbation will always decay and the fluid is always stable. Also, for this parameter set-up the group velocity v g ≤ 1, so the wave propagation is causal.
If we take the small k limit, Eqs. (4.26) and (4.27) yield the following modes: For this case the group velocity is observed to be same as the velocity for the ideal MHD. We analyse the causality of the system by following ref. [98] where it was shown that to guarantee the causality requires that the asymptotic value of the group velocity should be less than the speed of light. Alfvén mode in eq. (4.26) remains unaffected due to the bulk viscosity and hence always remain causal. For the magneto-sonic waves in the large k limit, we take the following ansatz ω = v L k in eq. (4.27) and collect terms in the leading-order of k, this yields v 4 L − xv 2 L + y = 0, (4.34) where, The velocities v L are Here, we see that unlike the small k limit, at large k the group velocity is affected by the transport coefficients. In order to have causal propagation, one demands v 2 L ≤ 1, which yields a causal parameter-set for the two branches, which correspond to the fast or slow magneto-sonic modes, Contour plot of the various causal regions is shown in Fig. 2, where b 1 is defined in eq. (3.4).
For the fast branch, we find that, although the asymptotic velocities depend on the magnitude of the magnetic field and the direction θ, the critical value, i.e., b 1 = 1.5 (red solid line), is independent of them. The slow branch is similarly B and θ dependent but moreover is causal throughout the parameter space.

MHD with shear viscosity
Many theoretical studies indicate that shear viscosity over entropy η/s has a minimum near the crossover temperature T co and rises as a function of temperature on both sides of T co ref. [129]. In this section, we consider a fluid with a non-zero shear viscosity but vanishing bulk viscosity. We may expect that the present scenario is applicable for the initial state of the QGP phase where T ∼ (4 -8)T co (for top RHIC and LHC energy in Au+Au and Pb+Pb collisions) and the bulk viscosity is vanishingly small in that temperature range. The energy-momentum tensor for a fluid with zero bulk and non-zero shear viscosity in a magnetic field takes the following form, According to the IS second-order theories of relativistic dissipative fluid dynamics, the spacetime evolutions of the shear stress tensor are given by eq. (2.11). For a given perturbation in the fluid, the energy-momentum tensor and the shear stress tensor can be decomposed as, π µν = π µν 0 + δπ µν . Where the perturbed energy-momentum tensor is, As usual, to solve the set of equations eq. (2.11), the conservation of the perturbed energymomentum tensor (eq. (4.41)), and eqs. (4.9)-(4.12) for obtaining the dispersion relation we write them in a matrix form where, δX = (δε, δũ x , δũ y , δũ z , δb x , δb y , δπ xx , δπ xy , δπ xz , δπ yy , δπ yz ) and the matrix A given in eq. (B.1). The det(A) = 0 gives, where, From eq. (4.43) we get two non-propagating and stable modes, Eq. (4.44) is a third-order polynomial equation and the analytic solution for this type of equations are discussed in appendix A. Eq. (4.45) is a sixth-order polynomial equation which is impossible to solve analytically. We can still gain some insight for a few special cases which are discussed below. For θ = 0, eq. (4.44) still remains a third-order polynomial equation and the coefficients of that polynomial can easily be obtained from eq. (4.44) as On the other hand, eq. (4.45) can be factorized into two third-order polynomial equations. The coefficients of one of such the third-order polynomial equation are whereas the coefficients of the remaining other third-order polynomial equation from eq. (4.45) are same as eq. (4.48) The roots of these third order polynomial equations are discussed in appendix A with the given X i s. We checked that the dispersion relations obtained from these equations with the coefficients given in eq. (4.49) are same as the sound mode in ref. [98].
For θ = π/2, one root of eq. (4.44) vanish and other two roots are of the form, The remaining three modes from eq. (4.45), are obtained from a cubic polynomial with X i 's given as: The corresponding roots can be calculated using the formula given in appendix A. The left panel of Fig. 3 shows the dependence of the imaginary parts of ω as a function of k/T and the right panel shows the group velocity as a function of k/T for different values of magnetic field for θ = 0. Various lines corresponds to different magnetic fields: qB = 0(blue lines), qB = 5m 2 π (green lines), qB = 20m 2 π (red lines). Fig. 4 shows the same thing but for θ = π 2 (eq. (4.52)). In the small k limit the dispersion relations that we get from eqs.  Note that, the first root have a degeneracy five.
In the large k limit we use the ansatz ω = v L k and keep only the leading-order terms in k, then the velocities v L are where,  The asymptotic causality condition for the shear-Alfvén mode can readily be obtained as, where b is defined in eq. (3.8). We observe that in this case the wave velocity and the causality conditions depend on both the magnitude of the magnetic field and direction of propagation of the perturbation. To explore the inter-dependency we show various causal regions as a function of b and θ as a contour plot in Fig. 5 (top left). We notice that the critical value of b at θ = 0 is b c = 1 and this value is independent of magnitude of the magnetic field. In the other extreme, i.e. for θ = π/2, the critical value is b c = [1+B 2 /(ε+P )] −1 , i.e., b c decreases with increasing magnetic field. In the limit of vanishing magnetic field it has been found in ref. [98] that for causal propagation of the shear modes b ≥ 1 should be satisfied. In the presence of magnetic field we found that, this constraint can be relaxed to even smaller values of b, given that the waves move obliquely. The causality constraint of the fast and slow waves in eq. (4.54) can be written in the form of (4.37). The simplified expression for the magneto-sonic modes can be written as We show various causal regions as a function of b and θ as a contour plot in Fig. 5 (top  right and bottom). The critical value of b, i.e., b c = 2 (obtained from (4.54)) is independent of the angle θ and the magnitude of magnetic field for the fast magneto-sonic mode. In the absence of a magnetic field this value coincides with that obtained for the sound mode in ref. [98]. Similarly, the slow magneto-sonic mode yields the critical value of b c = 1, independent of the angle θ and the magnitude of magnetic field B. It is still interesting to see that although the critical values of the fast and slow modes are B and θ independent, the asymptotic velocities are nevertheless dependent. Increasing the magnetic field, increases the asymptotic group velocities but the causal region always remains causal no matter how large the magnetic field becomes.

MHD with both bulk and shear viscosity
In this subsection, we investigate the stability and causality of a viscous fluid with finite shear and the bulk viscosity in a magnetic field. In heavy-ion collisions the initial magnetic field is very large and both shear and bulk viscosities are non-zero for the temperature range achieved in these collisions, hence the present case is most relevant to the actual heavy-ion collisions at top RHIC and LHC energies. The energy-momentum tensor is, (4.58) The small variation of the energy-momentum tensor due to the perturbed fields is, (4.59) Following the same procedure, as discussed in the previous two sections, we obtain the dispersion relations for the following independent variables, δX = (δε, δũ x , δũ y , δũ z , δb x , δb y , δπ xx , δπ xy , δπ xz , δπ yy , δπ yz , δΠ) T . (4.60) Following the usual procedure of linearisation we get a 12 × 12 dimensional square matrix A. By setting det A = 0 we have the following equations which subsequently give the dispersion relations, ω 7 + X 6 ω 6 + X 5 ω 5 + X 4 ω 4 + X 3 ω 3 + X 2 ω 2 + X 1 ω + X 0 = 0, (4.63) here, First, we find that eq. (4.61) gives two non-propagating modes of frequency ω = i τπ . Now, the eq. (4.62) is a third-order polynomial and can be solved analytically as discussed previously whereas the eq. (4.63) is a seventh-order polynomial equation and can not be solved analytically, therefore we lookout for the solution of these equations for some special cases discussed below.
For θ = 0, we obtain two cubic and a single quartic equations. The X i 's of the two cubic polynomials are same as in eq. (4.48). The dispersion relations for these cases are already discussed in the previous section, hence we will not repeat them here. The X i 's for the fourth-order polynomial equation are 65) and the corresponding roots can be calculated using the formula given in Appendix A.
For another case, we choose θ = π 2 , this time two of the roots turned out to be zero, and another two roots are the same as eq. (4.51). As before, we call these four modes as shear mode. The X i 's for the fourth-order polynomial equation are 66) and the corresponding roots can be calculated using the formula given in Appendix A. Note that the imaginary part of the propagating modes (obtained from eq. (4.66)) are degenerate and hence not shown separately in Fig. 6. The dash-dotted lines in the left panel of Fig. 6 correspond to the non-propagating modes generated due to the bulk viscosity, this is because in the small k limit they reduce to i τ Π , and in the same logic the dotted line corresponds to the non-propagating mode due to the shear viscosity. In general, we find that the (ω) is always positive for our set-up. So, for this parameter set, the fluid is always stable under small perturbation for non-zero bulk and shear viscosity. Also, we note another interesting point, when the magnetic field is increased the imaginary part of the propagating mode tends to zero i.e, the damping of the perturbation diminishes.
In the small k limit the dispersion relations from eqs. (4.61)-(4.63) become here also the first root have degeneracy of five. Similarly, in the large k limit using the ansatz ω = v L k, we obtain the asymptotic group velocities v L as: where, Now we are ready to explore the causality of a fluid in magnetic field. For this, we again check whether the asymptotic group velocity has super or sub luminal speed. We found that the theory as a whole is causal if the fluid satisfy the following asymptotic causality conditions for magneto-sonic waves: fast: (0 < y < 1) ∧ (2 √ y ≤ x < y + 1), From eq. (4.68) we find that a larger magnetic field gives a larger v L , but always remain sub-luminal given b and b 1 are larger than their corresponding critical values (discussed earlier). It is also clear from eq. (4.68) the asymptotic group velocity for non-zero bulk and shear viscosity is larger than the individual shear and bulk viscous cases. In Fig. 7 we show the contour plot of various causal regions as a function of b and θ. The critical line (red line) of the fast magneto-sonic mode (left panel) show that b and b 1 are inversely proportional. On the other hand, the causality condition for the slow magnetosonic waves is independent of b 1 . The critical value of b for the slow magneto-sonic mode is b c = 1.

Characteristic velocities for bulk viscosity
The characteristic curves can be seen as the lines along which any information is transported in the fluid, for example small perturbations, discontinuities, defects or shocks etc travel along one of these characteristic curves refs. [130][131][132]. Here we take the effect of nonlinearity in the propagation speed which is ignored in the linearisation procedure discussed earlier. Without the loss of generality we consider the (2 + 1)-dimensional case with only here, Q n = (ε, u x , u y , b x , B, Π) and R m = (0, 0, 0, 0, 0, Π). We parametrize the fluid velocity as u µ = (cosh θ, sinh θ cos φ, sinh θ sin φ, 0) and the b µ = (sinh θ, cosh θ cos φ, cosh θ sin φ, 0). The matrix elements of P t mn , P x mn , P y mn are given in Appendix B. We find the characteristic velocities (v ch x , v ch y ) by solving the following equations, For simplicity, here we take fluid in the LRF i.e, u µ = (1, 0, 0, 0) and the magnetic filed along the y-axis b µ = (0, 0, 1, 0). Then the characteristic velocities are, v ch , ± α + ζ τ Π (ε+P +Π) , (5.5) here h = ε + P + B 2 and the other roots are zero. The characteristic velocities obtain in eqs. (5.4), (5.5) are same with the eq. (4.36) for θ = π 2 and θ = 0, respectively provide Π = 0. So we conclude that, the asymptotic group velocity obtained by linearizing the MHD-IS equations is same as the characteristic velocities.

Results from the modified IS theory
So far all the results we discussed were obtained for viscous fluid in a magnetic field within the frame-work of the IS theory. However, in a recent work ref. [118] a modified form of the IS theory due to the magnetic field was derived. The modified theory which we call as NRMHD-IS from now on shows that the relaxation equation for the shear-stress tensor contains additional terms, here we neglected most of the terms and only keep the term which couples magnetic field and the shear viscosity. The simplified NRMHD-IS equation takes the following form τ π d dτ π <µν> + π µν = 2ησ µν − δ πB Bb αβ ∆ µν ακ g λβ π κλ . (6.1) Where δ πB is a new coefficient appearing only due to the magnetic field and b αβ = − αβγδ u γ b δ is an anti-symmetric tensor which satisfy b µν u ν = b µν b ν = 0. The rank-four traceless and symmetric projection operator is defined as ∆ µν Before proceeding further, a few comments on the NRMHD-IS equations are in order. It is well known that in the presence of a magnetic field, the transport coefficients split into several components, namely three bulk components and five shear components refs. [50,[116][117][118]. The information of these anisotropic transport coefficients are hidden inside the new coupling terms of the modified IS theory eq. (6.1). Note that the first-order terms on the right-hand sides are proportional to the usual shear-viscosity. These terms can be combined with the first-order terms on the left-hand side and, after inversion of the respective coefficient matrices, will lead to the various anisotropic transport coefficients. On the other hand, when solving the full second-order equations of the modified IS theory, one does not need to replace the standard viscosity with the anisotropic transport coefficients, since the effect of the magnetic field, is already taken into account by the new terms in these equations. Regarding modified second-order theory with finite bulk viscosity, we would like to mention that, there is still no existing theory that yields three distinct bulk components in Navier-Stokes limit (for details see ref. [118]) and it is still an open issue.
The last term of eq. (6.1) is the only non-trivial term added to the conventional IS theory for which we already discussed the results in previous sections. So, here we only consider the last term of eq. (6.1) and calculate the corresponding correction to the old results.
First, we add a perturbation to the new term which contributes to δπ µν , While calculating eq. (6.2) we use the fact that in the local rest frame the unperturbed shear stress tensor vanishes i.e., π µν 0 = 0, and as a consequence δB, δb µν terms are absent in eq. (6.2). For later use we define the projection of a four-vector A µ as A <µ> = ∆ µ ν A ν , which is orthogonal to u µ .
Using these new definitions we write the eq. (6.2) in a more simplified form as, In the LRF, the following components of the b µν are found to be non- where b µ taken as (0, 0, 0, 1). For the (3 + 1) dimensional case there are five independent equations for the shear stress according to the IS theory. For each five equations there are corresponding components of the δĨ µν which for our case are δĨ xx = −δ πB B 0 δπ xy , δĨ xy = 1 2 δ πB B 0 (δπ xx − δπ yy ), δĨ xz = − 1 2 δ πB B 0 δπ yz , δĨ yy = δ πB B 0 δπ yx and δĨ yz = 1 2 δ πB B 0 δπ xz . We include these new terms to the corresponding IS equations that we previously derived in section (4.3). Here also we get a 11 × 11 matrix. As usual, we derive the dispersion relations from det(A) = 0 which is a eleventh-order polynomial equation. Since finding the analytic solution of this polynomial equation is not possible, here we investigate some special cases.
In the hydrodynamical-limit i.e, in the small k limit we get the following modes Note that the frequency of a few non-hydrodynamic modes are changed due to the new coupling terms appearing in the NRMHD-IS theory.
For the large k limit we use the ansatz ω = v L k and take only the leading order terms in k which yields the following velocities here, 6) and the remaining roots are zero. Since the causality of the fluid depends on the asymptotic causality condition which here is given in eq. (6.5) and turned out to be the same as eq. (4.54). So it is clear that the causality condition remains same as eq. (4.56) whereas the dispersion relations gets modified.

Conclusions
The current work goes beyond the previous results of refs. [116,117] which used first order viscous MHD. As is well known the first order gradient terms in the energy-momentum tensor breaks causality, which is reflected from the existence of the superluminal mode. This prohibits the application of viscous MHD in relativistic systems and it is necessary to have rigorous treatment which the present work aims. The remedy is to go beyond the first viscous corrections in hydrodynamics, and to include second order terms as well.
We studied here the stability and causality of the relativistic dissipative fluid dynamics within the framework of the standard and modified IS theories in the magnetic field. By linearising the energy-momentum conservation equation, relaxation equations for viscous stresses, and the Maxwell's equations we obtain the dispersion relations for various cases.
In the absence of viscous stresses, the dispersion relation yields the well-known collective modes namely the Alfvén, slow and fast magnetosonic modes. For the bulk viscous case the Alfvén mode turned out to be independent of the bulk viscosity. The asymptotic causality constraint for the magneto-sonic modes is independent of the magnetic field and the angle of propagation. For the fast mode, the causality condition is the same as that previously derived in ref. [99] in the absence of magnetic field. The slow mode, on the other hand, remains causal throughout the parameter space. We also derived the causality bound with finite bulk viscosity using the full non-linear set of the equation using the method of characteristics and found that it agrees with the result obtained using small perturbations.
In the presence of shear viscosity, the causality constraint for the two magnetosonic modes is found to be independent of the magnetic field and the angle of propagation. Shear-Alfvén modes, on the other hand, do depend on them. We found that the causality constraint is changed in presence of a magnetic field. For the modified IS theory in the presence of shear viscosity, new non-hydrodynamic modes emerge but the causality constraint remains unaltered. Finally, in the presence of both shear and bulk viscosity, we have deduced the causal region of parameter space.
There are many possible directions for future work, namely, the study of causality bounds:(i) in resistive, second-order dissipative MHD where the electric field is non-zero and contributes in the equations of motion. ref. [119] (ii) theories which have spin degrees of freedom allows to include effects of polarization and magnetization ref. [133]. These and other interesting questions will be addressed in the future.

A Solutions of dispersion relations
In general, the hydrodynamic dispersion relations arise as solutions to P n (X 0 , X 1 , ..., X n−1 ) = 0, (A.1) where P = det A, is a n th order polynomial obtained from the determinant of matrix A after linearising the MHD equations. In this appendix, we enlist the roots of certain polynomials P n that we will encounter throughout this work. For n = 3, the polynomial P 3 is of the form and the corresponding roots are given as Here k = 1, 2, 3, ξ is the primitive cubic root of unity, i.e., ξ = −1+ √ −3 2 and the other variables are defined Similarly, for n = 4, the polynomial P 4 is of the form ω 4 + X 3 ω 3 + X 2 ω 2 + X 1 ω + X 0 = 0, (A. 5) and the corresponding roots are given as

B Details of matrix A defined in section 4.3 and the characteristic velocities
By linearising the energy-momentum conservation equations, Maxwell's equations and IS equation for shear viscosity, we write these in the matrix form as eq. (4.42). Here the form of matrix A is where f = 1 + iωτ π . Similarly we can write the matrix A for the modified IS theory, also for both the bulk and shear viscosity case.