Generalised global symmetries and magnetohydrodynamic waves in a strongly interacting holographic plasma

We begin the exploration of magnetohydrodynamics (MHD) in strongly coupled plasmas by constructing and analysing a holographic dual to a recent, generalised global symmetry-based formulation of dissipative MHD. The simplest holographic dual to the effective theory of MHD that was proposed as a description of plasmas with any equation of state and transport coefficients contains dynamical graviton and two-form gauge field fluctuations in a magnetised black brane background. The dual field theory, which is closely related to the large-$N_c$, $\mathcal{N} = 4$ supersymmetric Yang-Mills theory at (infinitely) strong coupling, is, as we argue, in our setup coupled to a dynamical $U(1)$ gauge field with a renormalisation condition-dependent electromagnetic coupling. After constructing the holographic dictionary, we compute the dual equation of state and transport coefficients, and for the first time analyse phenomenology of MHD waves in a strongly interacting, dense plasma with a (holographic) microscopic description. From weak to extremely strong magnetic fields, several predictions for the behaviour of Alfv\'{e}n and magnetosonic waves are discussed.


I. INTRODUCTION
Magnetohydrodynamics (MHD) is a hydrodynamic theory of long-range excitations in plasmas (ionised gases) (see e.g. [1][2][3][4]), which has been applied to systems ranging from the physics of fusion reactors to astrophysical objects. In the modern language of hydrodynamics formulated as an effective field theory [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19], MHD should describe the dynamics of infrared (IR) chargeneutral states in terms of massless effective degrees of freedom. These plasma ground states are characterised by an equation of state with a finite magnetic field. On the other hand, the electric field is suppressed due to the screening of electromagnetic interactions and is only induced on shorter length scales than the (thermodynamic) size of the system. In their standard form, the equations of motion that describe the evolution of plasmas are formulated as a combination of macroscopic fluid equations (continuity equation and the non-dissipative Euler, or dissipative Navier-Stokes equation), coupled to the microscopic electromagnetic Maxwell's equations. In ideal, non-dissipative form, the set of dynamical equations is The magnetic field is constrained by Eq. (1) is the continuity equation and Eq. (2) the Euler equation in the presence of the Lorentz force J × B, with J given by the low-frequency limit of the Ampere's law (∂ t E → 0) Eq. (3) is the Faraday's induction law with the electric field fixed by the assumption of the ideal Ohm's law which is derived by taking the conductivity in the (Lorentz transformed) Ohm's law J/σ = E+ v× B to infinity, i.e. σ → ∞. The constraint equation (5) is the magnetic Gauss's law. Since the ideal Ohm's law completely fixes E, the electric Gauss's law plays no role in the equations of MHD. Eq. (4) is the adiabatic equation of state relating density and pressure. Usually, one takes γ = 5/3. Altogether, Eqs. (1)- (4) give eight dynamical equations for eight unknown functions ρ, p, v and B, subject to the magnetic field constraint (5). While the above equations are closed, solvable and have been successfully applied to a variety of phenomena in plasma physics, they are only applicable within the specific assumptions used to construct them. This means that they are only valid for electromagnetism controlled by Maxwell's equations in the limit of ideal Ohm's law (no possibility of strong-field pair production, etc.) and for the specific equation of state in Eq. (4). This equation of state encodes a separation between the fluid and the charge carrying sectors, for which the justification, beyond assuming weakly coupled Maxwell electromagnetism, also assumes very weak interactions between the fluid degrees of freedom and electromagnetism inside the plasma. Concretely, the latter statement is reflected in the equation of state permitting no dependence on the magnetic properties controlled by the charged sector. Furthermore, because of a lack of a symmetry principle behind the construction of ideal MHD, these equations are difficult to extend unambiguously to the most general, higherorder, dissipative theory in the gradient expansion (the Knudsen number expansion) [20][21][22][23]. 1 As such, the traditional formulation of MHD lacks generality and cannot be compatible with a variety of IR effective theories of plasmas that could (in principle) be derived from quantum field theory, in particular, in the presence of a strong magnetic field.
These issues were addressed in a recent work [24], where MHD was formulated by following the effective field theory philosophy behind the construction of relativistic hydrodynamics (see e.g. [23,25]). Namely, MHD was formulated by only considering global conserved operators and writing them in terms of the most general hydrodynamic gradient expansion of the IR hydrodynamic fields [24]. 2 With such an expansion in hand, conservation equations then completely determine the temporal dynamics of a plasma with any equation of state. As in hydrodynamics, all of the details of the equation of state and transport coefficients are left to be determined by the microscopics of the underlying theory.
The two relevant global symmetries describing the long-range dynamics of a plasma were argued to give the stress-energy tensor T µν and a conserved anti-symmetric two-form current J µν [24]: While T µν corresponds to conserved energy-momentum, J µν is the manifestation of a generalised global U (1) symmetry, which can be sourced (and gauged) by a two-form gauge field b µν [33]. H ν µσ is a three-form field strength that can be turned on by an external two-form gauge field, H = db ext . This generalised global symmetry is a consequence of the absence of magnetic monopoles and directly corresponds to the conserved number of magnetic flux lines crossing a co-dimension two surface (in a four dimensional plasma). Normally, it is expressed in terms of the (topological) Bianchi identity where F = dA and A is the abelian electromagnetic field. In the language of a two-form current used in Eq. (9), The power in identifying Eq. (10) as a conservation equation of a global symmetry becomes apparent when one attempts to describe a phase of matter dominated by electromagnetic interactions, but without massless photons, i.e. the particles associated with A. In fact, this is precisely the situation in a plasma in which long-range electric forces are (Debye) screened and the photons are massive. Treating J µν as a globally conserved operator without invoking a massless gauge field A can then be used directly to organise the infra-red dynamics of such states [24]. We note that in the language of generalised global symmetries, photons only become massless particles when the (one-form) symmetry is spontaneously broken-i.e. photons are the Goldstone bosons present in 2 See also [26] and Ref. [27], which includes a valuable comparison of various related past works, such as [28][29][30].
For a new treatment of charged fluids in an external electromagnetic field, see [27,31]. Of further interest is also a recently proposed field theory description of polarised fluids [32].
the broken symmetry phase. 3 From this point of view, the Maxwell action is the effective Goldstone boson action that realises this symmetry non-linearly. The equations of motion (8) and (9) give seven dynamical equations (and one constraint). To solve them, we introduce the following hydrodynamical fields: a velocity field u µ , a temperature field T , a chemical potential µ that corresponds to the density of magnetic flux lines and a vector h µ , which can be though as a hydrodynamical realisation of a fluctuating magnetic field. The vectors are normalised as u µ u µ = −1, h µ h µ = 1, u µ h µ = 0, together resulting in 10 − 3 = 7 degrees of freedom. A (directed) velocity flow of the plasma breaks the Lorentz symmetry from SO(3, 1) to SO (3), which is further broken by the additional vector (magnetic field) to SO (2). 4 The projector transverse to both u µ and h µ is defined as ∆ µν = g µν + u µ u ν − h µ h ν and has a trace ∆ µ µ = 2. The constitutive relations for the conserved tensors with a positive local entropy production [17] and charge conjugation symmetry can now be expanded to first order in derivatives as [24] T µν = (ε + p) u µ u ν + p g µν − µρ h µ h ν + δf ∆ µν + δτ h µ h ν + 2 (µ h ν) + t µν , where The thermodynamic relations between ε, p and ρ, which need to be obeyed by the equation of state p(T, µ) are dp = s dT + ρ dµ .
Furthermore, for the theory to be invariant under time-reversal, the Onsager relation implies that ζ × ≡ ζ × . Thus, first-order dissipative corrections to ideal MHD are controlled by seven transport coefficients: η ⊥ , η , ζ ⊥ , ζ , ζ × , r ⊥ and r . Each one can be computed from a set of Kubo formulae presented in [24,27] and reviewed in Appendix A. The transport coefficients should obey the following positive entropy production constraints: η ⊥ ≥ 0, η ≥ 0, r ⊥ ≥ 0, 3 We note that the order parameter that distinguishes between a broken and an unbroken magnetic one-form symmetry is an expectation value of the 't Hooft loop operator. When the symmetry is preserved, then the expectation value of the loop operator obeys the area law, WC ∼ exp {−T Area[C]}. On the other hand, in the symmetry broken phase with massless photons, the expectation value obeys the perimeter law, WC ∼ exp {−T Perimeter[C]} [33]. 4 Note that at zero temperature, in a plasma with a non-fluctuating temperature field, the symmetry is enhanced to SO(1, 1) × SO(2) [24]. r ≥ 0, ζ ⊥ ≥ 0 and ζ ⊥ ζ ≥ ζ 2 × . In absence of charge conjugation symmetry, the theory has four additional transport coefficients, resulting in total in eleven transport coefficients [27]. The precise connection between the above formalism of MHD using the concept of generalised global symmetries and MHD expressed in terms of electromagnetic fields, which match in the limit of a small magnetic field (compared to the temperature of the plasma), was established in Ref. [27].
Since the effective theory [24] makes no assumption regarding the microscopic details of the plasma, then, should such details somehow be computable from quantum field theory, or otherwise, the effective MHD can be used in solar plasma physics, fusion reactor physics, astrophysical plasma physics and even QCD quark-gluon plasma resulting from nuclear collisions. Of course, computing the microscopic properties of such systems is extremely difficult. In this work, we will resort to using holographic duality. By using standard holographic methods applicable to hydrodynamics [34][35][36], our analysis will provide us with the required microscopic data of a strongly interacting toy model plasma needed to describe the phenomenology of MHD waves.
The paper is structured as follows: first, in Section II, we review important aspects of gauge theories with a matter sector coupled to dynamical U (1), which can describe a plasma in the IR limit. In particular, we focus on the discussion of how to couple a strongly interacting field theory with a holographic dual to dynamical electromagnetism, all within a holographic setup. Then, in Section III, we explore this holographic setup in detail, develop the holographic dictionary and use it to compute the microscopic properties of the dual plasma, i.e. the equation of state and firstorder transport coefficients. In Section IV, we then use this data to analyse the phenomenology of propagating MHD modes-Alfvén and magnetosonic waves. Finally, we conclude with a discussion and a summary of the most important findings in Section V. Three appendices are devoted to a derivation of the relevant Kubo formulae (Appendix A), details regarding the derivation of horizon formulae for the transport coefficients (Appendix B) and a derivation of the magnetosonic dispersion relations (Appendix C).
Note added: We note that in addition to this paper on the holographic dual of MHD from the perspective of generalised global symmetries, a closely related work, i.e. Ref. [37], also studies various aspects of generalised global symmetries in gauge-gravity duality and holographic dual(s) of [24]. Although the two concurrent and complementary papers focus on different aspects of holography, there is overlap between our Sections II B and III and parts of Ref. [37].

II. MATTER COUPLED TO ELECTROMAGNETIC INTERACTIONS
A microscopic theory from which an effective description of a plasma can arise comprises of a matter sector that interacts through an electromagnetic U (1) gauge field. The simplest example of such a theory is quantum electrodynamics. In other theories, the matter sector may itself exhibit complicated physics with additional gauge interactions, such as in QCD. In this work, the theory that we will study contains an infinitely strongly coupled holographic matter sector (closely related to N = 4 supersymmetric SU (N c ) Yang-Mills) with infinite N c . Because of the coupling between matter and dynamical electromagnetism, the holographic setup and the interpretation of results is somewhat subtle. For this reason, we begin our discussion by reviewing some relevant aspects of quantum field theory in a line of arguments similar to [38].

A. Quantum electrodynamics
The simplest example of a theory coupling matter to electromagnetism is quantum electrodynamics (QED). QED is a U (1) gauge theory that contains a (massive) Dirac fermion ψ (describing electrons and positrons) and a massless photon field A µ : 5 D µ is the gauge covariant derivative that couples A µ to the fermion current (with the coupling e scaled out from the interaction). For a detailed discussion of various properties of QED, see e.g. [39][40][41].
The stress-energy tensor of the theory is In the massless limit (m = 0), the theory is classically scale invariant, which is reflected in the vanishing trace of the stress-energy tensor, T µ µ = 0. Quantum mechanically, the theory does not remain scale invariant. The trace receives a correction proportional to the beta function of the electromagnetic coupling, This is the anomalous breaking of scale invariance-the so-called trace anomaly. The running electromagnetic coupling e(µ) depends on the renormalisation group scale µ. To first order in perturbation theory, the beta function is which, integrated on the interval µ ∈ [M, Λ], gives the running coupling Here, M is some IR RG scale at which the electric charge takes the renormalised physical value, e r = e(M ), and Λ is the UV cut-off. Note that at the Landau pole, Λ = Λ EM , the left-hand-side of (26) vanishes. On the other hand, the expectation value of the stress-energy tensor is a physical quantity and therefore cannot depend on µ. This statement is encoded in the following identity, which leads to the Callan-Symanzik equation: Since we are interested in neutral IR plasma states in QED that can be described by an effective theory of MHD, we can consider the (ground state) expectation value of the photon field to produce a non-zero magnetic field and a vanishing (screened) electric field, B is the magnitude of the "background" magnetic field pointing in the x 3 = z direction. The IR spectrum of the theory has a gapped-out photon, i.e. long-range charge neutrality, which allows us to neglect quantum fluctuations of A µ . For such a plasma state, Eq. (24) yields Furthermore, the expectation value of the stress-energy tensor can be conveniently split into the matter (containing matter-light interactions) and the purely electromagnetic parts, where in the second line, we chose to evaluate the expectation value at the UV cut-off µ = Λ. Note that because T µν is µ-independent (cf. Eq. (27)), this choice does not influence the final value of T µν .

B. Strongly interacting holographic matter coupled to dynamical electromagnetism
We now turn our attention to the holographic strongly interacting theory that will be investigated in the remainder of this paper. Throughout our discussion, it will prove useful to think of the matter sector as that of the best understood holographic example-the conformal N = 4 supersymmetric Yang-Mills theory (SYM) with an infinite number of colours N c and an infinite 't Hooft coupling λ. However, as will become clear below, the theory dual to our holographic setup will not be precisely the N = 4 SYM coupled to a U (1) gauge field, but rather its deformation, of which the microscopic definition will not be investigated in detail. Instead, the model studied here should be considered as a bottom-up construction-the simplest dual of a strongly coupled plasma, which can be described with magnetohydrodynamics in the infrared limit.
The field content of N = 4 SYM are four Weyl fermions, three complex scalars and a vector field, all transforming under the adjoint representation of SU (N c ). The theory also has an SU (4) R R-symmetry owing to its extended supersymmetry. The adjoint fields together represent the matter content of a hypothetical plasma, which further requires the fields to be (minimally) coupled to an electromagnetic U (1) gauge group (with e the electromagnetic coupling). In N = 4 SYM, this can be achieved by gauging the U (1) R subgroup of SU (4) R . Under U (1) R , the Weyl fermions transform with the charges {+3, −1, −1, −1}/ √ 3 and the complex scalars all have charge +2/ √ 3 (for details regarding the choice of the normalisation, see [38]). Such a system can be considered as a strongly coupled toy model for a QCD plasma in which the quarks interact with photons as well as with the SU (3) vector gluons. A crucial fact about N = 4 SYM is that the R-current of N = 4 becomes anomalous in the presence of electromagnetism. For this reason, the U (1) R , which is gauged, is also anomalous and thus the theory has to be deformed in some way to reestablish its self-consistency. As pointed out in [38], one way to do this is by adding a set of spectator fermions that only interact electromagnetically and "absorb" the anomaly. We will assume that the gauge anomaly can be cancelled by some deformation of the theory so that the quantum expectation value of the U (1) R R-current J µ R remains conserved, ∇ µ J µ R = 0. We can then write the total bare action of the SU (N c ) × U (1) gauge theory as where A µ is the dynamical electromagnetic gauge field and F = dA. The expectation value of the conserved operator J µ R contains a trace over the colour index of the adjoint matter field and therefore scales as N 2 c . Since it is coupled to a single photon, the Maxwell part of the total plasma action S plasma contains no powers of N c .
As in the QED plasma, we will consider the photons to be gapped out from the IR spectrum so that A µ will only produce a (classical) magnetic field In order to maintain the neutrality of the plasma, we will set the electric U (1) R chemical potential to zero, µ R = A 0 = 0. 6 For this reason, the electric one-form (or vector) conserved U (1) R R-current will play no role in the hydrodynamic IR limit of the theory, so J µ R = 0. The plasma has a conserved stress-energy tensor to which both the matter (along with its interaction with the electromagnetic field) and the purely electromagnetic sectors contribute, The trace of the superconformal theory again experiences an anomaly proportional to the beta function of the electromagnetic coupling (cf. Eq. (29)), which in N = 4 theory turns out to be one-loop exact in the presence of a background electromagnetic field and follows from a special case of the NSVZ beta function due to the fact that the U (1) R sector has a remaining N = 1 supersymmetry (see Refs. [38,44,45] The beta function for the inverse electromagnetic coupling is then with the fermionic and the scalar R-charges being q α f = {+3, −1, −1, −1}/ √ 3 and q α s = {2, 2, 2}/ √ 3, respectively. In analogy with Eq. (26) in QED, by integrating the beta function equation, we find It is essential to stress that even though our holographic theory will not be exactly dual to the N = 4 SYM theory, it will give us the same trace anomaly and thus the same electromagnetic beta function. Since the NSVZ beta function (35) is only sensitive to the matter content, this match can be interpreted as our working with a theory with the U (1)-gauged matter content and R-charges of N = 4 but with a deformed Lagrangian and possibly additional matter that is ungauged under the U (1). Beyond the stress-energy tensor of the theory discussed thus far, the only other (generalised) global symmetry of interest to describing a plasma state is the higher-form U (1) symmetry that corresponds to the conserved number of magnetic flux lines crossing a two-surface. The symmetry results in a conserved two-form current J µν = 0 and was discussed in Section I. The generating functional of the field theory that can be used to study MHD of a magnetised plasma in which the two globally conserved operators are T µν and J µν is The simplest holographic dual of such a state is one that contains a five-dimensional bulk with a dynamical graviton (metric tensor G ab ) described by the Einstein-Hilbert action, a negative cosmological constant and a two-form bulk gauge field B ab : 8 In standard (Dirichlet) quantisation, the two fields asymptote to g µν and b µν at the boundary and source T µν and J µν . Furthermore, H is the three-form defined as H = dB. In component notation, B = 1 2 B ab dx a ∧dx b and H = 1 6 H abc dx a ∧dx b ∧dx c . The two-form gauge field action is the bulk Maxwell Lagrangian F ∧ F written in terms of the five-dimensional Hodge dual three-form H = F , giving the Lagrangian term H ∧ H. In most of our work, we will set e H = L = 1. Because the two bulk theories are related by dualisation, the background solution to the equations of motion derived from (38) give rise to the same magnetised black brane solution known from the Einstein-Maxwell theory [46].
In the absence of the two-form term, the action (38) arises from a consistent truncation of type IIB string theory on S 5 and is upon identification of the Newton's constant κ 5 = 2π/N c dual to pure N = 4 SYM at infinite N c and infinite 't Hooft coupling λ. For reasons discussed above, the full dual of the action (38) is unknown and we are not aware of a mechanism for deriving this action from a consistent truncation of ten-dimensional type IIB supergravity. Nevertheless, for purposes of comparing the sizes of matter and electromagnetic contributions to the total operator expectation values, it will prove useful to keep the definition of κ 5 in terms of the number of colours N c of the hypothetical dual deformed N = 4 SYM coupled to dynamical electromagnetism.
To show further evidence that the action (38) is a sensible dual of a strongly coupled MHD plasma, it is useful to elucidate the connection between Eq. (38) and the Einstein-Maxwell theory. To put an uncharged holographic theory in an external magnetic field, one normally adds the Maxwell action F 2 with F = dA to the Einstein-Hilbert bulk action. If one imposes the Dirichlet boundary conditions on the bulk one-form A a , then A a sources the R-current J µ R at the boundary, d 4 x J µ R δA µ , and thus the electromagnetic field A µ is external and non-dynamical. The investigation of the physics of such a setup was initiated in [46] and studied in numerous subsequent works, including recent [29,30,38,47,48]. Instead, one can work in alternative quantisation and impose Neumann boundary conditions on A a . Such a choice exchanges the interpretation of the normalisable and the non-normalisable mode in A a . From the dual field theory point of view, this can be interpreted as the Legendre transform of the boundary coupling, leading to the variation d 4 x A µ δJ µ R . Physically, this means that in alternative quantisation, an external current sources a dynamical (boundary) vector field (see e.g. [49,50]). The two boundary theories, one with Dirichlet and one with Neumann boundary conditions, are related by a double-trace deformed RG flow. As we will see in Section III, the correct way to interpret our theory will be in terms of a double-trace deformed CFT (with a broken scale invariance), which will require us to impose mixed boundary conditions [51] on the two-form field. 9 Since J µ R is conserved, one can express it through an anti-symmetric b µν as µνρσ ∂ ν b ρσ , which, upon integration by parts, yields a dualised d 4 x J µν δb µν , where J µν is the anti-symmetric current from Eq. (11). From the point of view of the bulk, as in a lower-dimensional theory [56], the Einstein-Maxwell bulk (quantum) path integral runs over the metric and the Maxwell field A a . Alternatively, one can write the path integral over the fields strength F ab , but at the expense of ensuring the Bianchi identity dF = 0 by introducing a Lagrange multiplier B ab : Since the second (Bianchi identity) term vanishes for any classical field solution, it has no influence on the saddle point of the path integral. However, it does generate a non-zero contribution to the boundary action, i.e. N 2 c /8π 2 e H d 4 x µνρσ F ρσ B µν , which is precisely the source term d 4 x J µν b µν once we identify B µν ∼ b µν and J µν ∼ µνρσ F ρσ . The precise dictionary between the bulk and boundary quantities will be discussed in Section III B. By varying the action with respect to F ab , one obtains the equation of motion Then, the field strength F ab can be integrated out in the saddle point approximation which gives the two-form gauge field Lagrangian term from Eq. (38). Furthermore, in the language of the Einstein-Maxwell theory, by using Eq. (40), one finds the relation between the one-form R-current 9 See Refs. [52][53][54][55] and references therein.
J µ R and the B ab field: where u is the radial coordinate and u = 0 the boundary of the bulk spacetime. Thus, imposing the Dirichlet boundary condition on B ab (in this case, with an additional double-trace deformation) corresponds to treating J µ R as a source, which is the same as performing alternative quantisation discussed above. This is consistent with the interpretation that the dual field theory of (38) contains dynamical photons. Furthermore, as we will see from a detailed holographic renormalisation in Section III B, the (double-trace) boundary counter-terms, which are required to keep the on-shell action finite, will give us precisely the Maxwell theory for A µ (dual of b µν ) on the boundary, including a renormalised electromagnetic coupling e r , as in QED. 10 All further details of this holographic setup will be presented in Section III.

III. HOLOGRAPHIC ANALYSIS: EQUATION OF STATE AND TRANSPORT COEFFI-CIENTS
In this section, we study the relevant details of the simplest holographic theory with Einstein gravity coupled to a two-form bulk field, cf. (38), which can source a two-form current associated with the U (1) generalised global symmetry in the boundary theory. As our goal is to study the phenomenology of MHD waves in a strongly coupled plasma using the dispersion relations of [24], we will use holography only to provide us with the necessary microscopic data: the equation of state and the transport coefficients.
In Section III A, we will begin by discussing details of the magnetic brane solution [46,58] supported by the bulk action introduced in Section II B. In Section III B, we will consider holographic renormalisation of the theory in question and show how the bulk gives rise to a dual theory coupled to dynamical electromagnetism (as in Section II). In particular, we will derive the expectation values of the stress-energy tensor T µν and the two-form J µν and show that they satisfy the Ward identities (8) and (9). We will also recover and match all expected renormalisation group properties, such as the beta function of the electromagnetic coupling, from the point of view of the bulk calculation. In Section III C, we will then compute and analyse thermal and magnetic properties of the equation of state of the dual plasma. Finally, in Section III D, we will derive the membrane paradigm formulae for the seven transport coefficients required to describe first-order dissipative MHD [24] and compute them. 11 Further details regarding the horizon formulae for the transport coefficients can be found in Appendix B. 10 We note that the way the Maxwell Lagrangian arises on the boundary is equivalent to the way holographic matter can be coupled to dynamical gravity on a cut-off brane [57]. There too, a holographic counter-term gives rise to the Einstein-Hilbert action at the cut-off brane (the boundary) of a more intricately foliated bulk. As shown by Gubser in [57], such a theory can result in a radiation (CFT)-dominated FRW universe at the boundary with the stress-energy tensor of the N = 4 SYM driving the expansion. 11 These horizon formulae are analogous to the expression for shear viscosity in N = 4 theory [59]. For more recent discussions of other transport coefficients that can be computed directly from horizon data, see e.g. [60][61][62][63][64][65][66].

A. Holographic action and the magnetic brane
A holographic action dual to a plasma state with a low-energy limit that can be described by MHD was stated in Eq. (38). Including the boundary Gibbons-Hawking term and the (relevant) holographic counter-term, the full action is where trK is the trace of the extrinsic curvature of the boundary (∂M ) defined by an outward normal vector n a . For convenience, we now set e H = 1. The two-form H µν is defined as a projection of the three-form field strength in the direction normal the boundary, H µν = n a H aµν . C is a dimensionless number that needs to be adjusted to fix the renormalisation condition, of which the details will be discussed below. The equations of motion that follow are Since the theory (42) is S-dual to the Einstein-Maxwell theory, we can express the magnetised black brane solution of [46] by dualising the Maxwell terms and writting The equations of motion (43) for this ansatz reduce to three second-order ordinary differential equations (ODE's) for {F, V, W} and one additional first-order constraint. The equation of motion derived from the variation of the two-form (44) is automatically satisfied. The equations are equivalent to those derived from the Einstein-Maxwell theory [46] upon identification of the Maxwell bulk two-form The undetermined functions F , V and W are can be obtained numerically by using the shooting method. We first expand the background fields near the horizon as 12 The metric ansatz is chosen to have the form used in [47]. It can be obtained from the ansatz ds 2 = −U (r)dt 2 + e 2V (r) (dx 2 + dy 2 ) + e 2W (r) dz 2 + dr 2 /U (r) used in [46] by performing a coordinate transformation r = r h / √ u and shifting V and W by constant − ln v and − ln w, respectively, which are chosen so that the near-boundary expansion has the form after solving the equations of motion order-by-order near the horizon. The scaling symmetry of our background ansatz then allows us to rescale dx and dy so that v h 0 = w h 0 = 0. Next, we match the numerical solutions generated by shooting from the horizon towards the boundary, where the analytical near-boundary expansions of the metric functions are As before, one can solve for the coefficients can be removed by residual diffeomorphism freedom of the metric ansatz [47]. For a given value of B = vB/r 2 h , we can therefore generate a numerical background by shooting from the initial conditions of the functions set by the near-horizon expansion with The numerical value off is chosen so that the near-boundary expansion has f b 1 = 0. The near-boundary behaviour of this function then determines the properties of the dual field theory. Note that the theory is governed by a one-parameter family of such numerical solutions characterised by the is the Hawking temperature (see Eq. (84)). In practice, this ratio can be tuned by changing the parameter B of the background ansatz. The numerical solver encounters stiffness problems when B ≈ √ 3, i.e. where the temperature is close to zero. All of our numerical results will therefore stop near T / √ B = 0. In this work, we do not attempt an independent analysis of the theory at T = 0.

B. Holographic renormalisation and the bulk/boundary dictionary
The next step in analysing the dual of (42) is a systematic holographic renormalisation. In this section, we derive the one-point functions of the stress-energy tensor T µν and the two-form current J µν , and show that they satisfy the Ward identities of magnetohydrodynamics (8) and (9) [24], which in terms of operator expectation values take the form whereH = db is the field strength of the background gauge field b in field theory. The precise definition of these quantities will become clear below. Since we are only interested in the expansion of MHD to first order in the gradient expansion around a flat (boundary) background, it will be sufficient to only work with terms that contain no more than two derivatives along the boundary directions. The procedure for obtaining holographic renormalisation will closely follow Refs. [67,68]. 13 We begin by writing the bulk metric in the Fefferman-Graham coordinates [67] so that near the boundary, ρ ≈ 0, the metric g µν can be expanded as Note that Greek (boundary) indices in a tensor A µν are raised with the metric g µν (0) , which satisfies g There are two types of covariant derivative that we will use: ∇ µ are defined with respect to the metric g µν (ρ, x), while ∇ µ and ∇ µ ≡ g µν (0) ∇ µ are defined through the metric g µν (x). The Ricci tensors of g µν and g The components of bulk two-form gauge field B ab in the boundary field theory directions can similarly be expanded near the boundary as In the boundary directions, the three-form field strength is defined as in the same way at each order. Note that both quantities B (0) µν and B (1) µν are related to the two-form gauge field source of the boundary theory, d 4 x √ −g J µν δb µν . The variation of the regularised bulk on-shell contribution from the H 2 term, evaluated at the boundary cutoff ρ = ρ Λ , is The expectation value of the operator sourced by B In holography, such boundary conditions are knows as mixed boundary conditions. They arise in the presence of double-trace deformations [51], which is precisely how the logarithmically running H 2 term in the renormalised on-shell action should be interpreted. From the point of view of the boundary field theory, as we will see below, this term is a consequence of dynamical boundary electromagnetism-it is the boundary Maxwell action. Now, since the source b µν is a physical quantity, it cannot depend on the cut-off scale ρ Λ . Hence, the renormalisation group equation prompts us to set the value of C ∼ 1/ √ ρ Λ , which makes the on-shell action finite in the limit of ρ Λ → 0. As we will see below, the proportionality constant in the relation between C 2 and ρ Λ sets the value of the renormalised electromagnetic coupling (or the Landau pole scale), and corresponds to the choice of the renormalisation group condition. With these boundary conditions in hand, the expectation value of J µν can then be obtained by taking a variational derivative of the on-shell action with respect to the source b µν .
The Ward identities (48) can be obtained by solving the equations of motion (43) and (44) [68]. In Fefferman-Graham coordinates (49), these equations (together with the trace of (43)) become where g −1 denotes the matrix inverse of g (in components, this is g µν ) and where Expanding equations (55)-(59) around small ρ, we find that Since g (1) µν is proportional to second derivatives of the boundary metric, and we are only keeping track of terms up to second orders in boundary derivaties, we can ignore terms with g (1) µν . The remaining equations of motion can thus be written as The expectation values of the stress-energy tensor and the two-form current follow from the generating functional (37): In holography, W is computed from the (on-shell) action (42), giving us 14 Note that the expectation value of T µν scales as N 2 c , the expectation value of J µν is of order O(1). By using Eq. (62) and the fact that H µν = n ρ H ρµν = −2B (1) µν + O(ρ), we find that the boundary two-form current is conserved: Using the definition (11), which gives J µν = 1 2 µνρσ F ρσ and connects Eq. (68) with the Bianchi identity, we find that B (1) sets the expectation value of the Maxwell field strength F . Furthermore, the (regularised) stress-energy tensor (66) becomes where Λ is the UV cutoff of the theory. As discussed in Section II, the choice of the dimensionful constant C must now be made in order to fix the renormalisation condition, which will render the renormalised expectation value T µν physical and finite in the limit of Λ → ∞. To see how the constant C in Eq. (69) is related to our discussion in Section II, we write the last term by introducing a mass scale M : What can be seen from Eq. (70) is that this splitting precisely reproduces the way the logarithmic divergence enters into the stress-energy tensor from two different pieces of the Lagrangian: the matter content (with its coupling to the photons) and the electromagnetic (Maxwell) part: with the two terms being Finally, we note that the electromagnetic T EM µν would follow precisely from the Maxwell boundary action, which induces a double-trace deformation into the boundary field theory (see discussion below Eq. (53)) upon using Eq. (64) and the fact that the bulk B (1) determines F µν : where the last equality follows from the fact that quantum fluctuations of the photon field are suppressed in the boundary QFT. Our holographic calculation thus fully reproduces Eq. (33), which followed from the field theory discussion in Section II B. Furthermore, the running electromagnetic coupling constant matches the one found from field theory (cf. Eq. (36)) [38]. Hence, our holographic setup appears to contain the U (1)-gauged matter content of the N = 4 SYM theory.
In terms of bulk quantities, the renormalised stress-energy tensor and the two-form current are where, as in Section II, e r is the renormalised coupling which needs to be set by experimental input-the renormalisation condition. In practice, this constant is fixed by choosing the value of C in (66). For the same reasons as in QFT, there is therefore an inherent ambiguity in holographic results, which has to be fixed by external physically-motivated input.
We conclude this section by noting that the relation (63) and a relation betweenH, H (0) and H (1) implies that the Ward identity for the stress-energy tensor satisfies Eq. (8), or in terms of our holographic notation, ∇ ν T µν =H µ λσ J λσ , as in Eq. (48).

C. The equation of state
To find the equation of state of our theory, we can use the renormalised stress-energy tensor (76) and the two-form current (77) computed in the previous section. Expressed in terms of the near-boundary expansions (47), which can be read off from the numerical background, we find where we have used the (renormalised) fine-structure constant α = e 2 r /4π of the electromagnetic coupling in the plasma, which, as e r , is fixed by the choice of C in (66), and rescaled it by N 2 c /2π 2 (or |β(1/e 2 )|) asᾱ The couplingᾱ (or alternatively, the Landau pole) has to be fixed by experimental observations as in any other quantum field theory, which is not easy in an unrealistic toy model. In studying strongly coupled MHD, it is phenomenologically relevant to not only consider the matter and light-matter interactions, but to also include large electromagnetic self-interactions encoded in the Maxwell action. However, since we are working with a holographic large-N c matter sector and a single photon, it is unnatural to expect a Maxwell term of the same order. The choice that we make here is to set the rescaled constantᾱ to the physically motivatedᾱ = 1/137. There are several ways to think about this choice: one is imagining that our plasma contains magnetic properties, which have non-trivial scalings with N c , while another interpretation may assume that the bulk studied here could remain a valid dual of a theory with a reasonably small N c . Of course, by considering only a classical bulk theory, we are restricting the strict validity of any computed observable to the limit of N c → ∞. As soon as one moves towards finite N c , it becomes crucial to estimate the size of subleading 1/N 2 c corrections (topological expansion in the string coupling g s )-an endeavour in holography (and string theory) which to date has been largely neglected and will continue being neglected in this work. 15 A less problematic limit is that of the infinite 't Hooft coupling, which is also implied by the choice of our action. 16 Perhaps the best interpretation is one of an "agnostic choice" led by our having to fix a free parameter to some value. We will return to a more careful investigation of the dependence of our results on this choice in Section IV C.
The expectation values of the stress-energy tensor expressed in (78)-(80) are related to the MHD stress-energy tensor in Eq. (12) by We note that, as required in a conformal field theory with a trace anomaly induced by electromagnetic interactions, the trace of stress-energy tensor is non-zero. The holographic two-form current, is related to the equilibrium magnetic flux line density appearing in the MHD equation (13) as J tz = ρ. Temperature and entropy can be expressed in terms of the background geometry as and are therefore independent of the renormalised electromagnetic charge. The chemical potential, which corresponds to the density of magnetic flux lines, can be computed by using the thermodynamic identity ε + p = sT + µρ (cf. (20)): 15 For some discussions of 1/N 2 c corrections to the thermodynamic free energy (the equilibrium partition function) and hydrodynamic long-time tails, see [70][71][72][73][74]. 16 For recent discussions of coupling-dependent holography, see [75][76][77][78][79] and references therein. c ). Using the above relations, we can perform two consistency checks on our holographic setup and numerical calculations of the background. First, the value of the pressure computed from the stress-energy tensor component T xx = p can be compared with the value of the Euclidean on-shell action, p = −i(βV 3 ) −1 S on−shell , where β = 1/T and V 3 is the spatial volume of the theory. Secondly, we can compute ε + p − µρ from the stress-energy tensor evaluated near the boundary and by using the thermodynamic relation (20), check whether its values agrees with sT computed purely from horizon quantities. Both calculations show consistency of our setup in that we find T xx = −i(βV 3 ) −1 S on−shell and T tt + T zz = sT , within numerical precision.
We can now plot various thermodynamic quantities in a dimensionless manner by dividing them by appropriate powers of B. The natural dimensionless parameter with respect to which we present our numerical results is T / √ B. The results for the energy density, pressure, entropy density and chemical potential are shown in Figure 1. The theory has two distinct regimes: the lowand the high-temperature regimes, or alternatively, the strong and weak magnetic field regimes, respectively. The high-temperature regime T / √ B 1 is one to which MHD has been historically applied and to which the formulation of MHD, which assumes a weak-field separation between fluid and charge degrees of freedom can be applied. The claim presented in the Ref. [24] is that within the dual formulation, however, MHD applies for all values of T / √ B provided that the state remains in the hydrodynamic regime. The profiles of the thermodynamic functions in Figure 1 show a smooth crossover between the two regimes, which occurs around By using numerical fits, the equation of state in the two limits behaves as expected on dimensional grounds [24]. We present our numerical results in Table I.
We also note that the value of the pressure at low temperature strongly depends on the renormalised (re-scaled) fine structure constantᾱ, which we set toᾱ = 1/137.

D. Transport coefficients
Next, we compute the seven transport coefficients, η ⊥ , η , r ⊥ , r , ζ ⊥ , ζ and ζ × , by using the Kubo formulae derived in [24,27] and reviewed in Appendix A. The procedure only requires us to turn on time-dependent fluctuations of the background fields without any spatial dependence, G ab → G ab + δG ab (t) and B ab → B ab + δB ab (t). The perturbations asymptote to the boundary sources δg µν of the dual stress-energy tensor and the two-form current. In the absence of spatial dependence, the fluctuations decouple into five separate channels, from which the seven transport coefficients are computed, with each channel containing one independent dynamical second-order equation. The sets of decoupled fluctuations responsible for their respective transport coefficients are η ⊥ : δG xy , η : δG xz , δB tx , δB xu , ζ ⊥ , ζ , ζ × : δG tt , δG xx , δG yy , δG zz , δB tz , δG tu , δB zu , δG uu , r ⊥ : δB xz , δG tx , δG xu , r : δB xy , with only one of the three bulk viscosities being independent. Each one of the transport coefficients can then be related to a membrane paradigm formula and computed entirely in terms of the horizon quantities. We summarise these relations here and discuss their derivation below: where b (−) and Z (−) are the time-independent solutions of the fluctuations δB xz and Z s = δG x x + δG y y − (2V /W )δG z z , respectively. The arguments denote that the functions are evaluated either at the horizon, u = 1, or the boundary, u = 0. Note that the value at the boundary is set by the Dirichlet boundary conditions.
What we see is that the ratio of the transverse shear viscosity (w.r.t. the background magnetic field) to entropy density is universal, resulting in η ⊥ /s = 1/4π. Furthermore, the expressions for η and r only depend on the background quantities v and w, while ζ ⊥ , ζ and r ⊥ also depend on the fluctuations of the fields. 17 In order to derive the horizon formulae, we use the Wronskian method (see e.g. [64]). Here, we will only explicitly show the derivation of the transverse resistivity r ⊥ . The other formulae from Eq. (88) are derived in Appendix B. First, we combine the equations of motion for the relevant fluctuations, δB xz , δG tx , and δG xu , into a single second-order differential equation by eliminating the metric fluctuations, Since we are only computing first-order transport coefficients, it is sufficient to solve Eq. (89) to linear order in ω. To find the solution, we assume that there exists a time-independent solution b (−) xz (u), which asymptotes to a constant both at the boundary and the horizon. At the boundary, this asymptotic value is related to the source of the two-form background gauge field, i.e. b xz . It can be expressed as an integral over the Wronskian W R of (89): 17 For a holographic derivation of bulk viscosity in neutral relativistic hydrodynamics, see [60]. where The near-boundary and the near-horizon expansions of b Finally, δB xz (ω, u) is then the following linear combination of the two solutions: The coefficient α(ω) can be determined by imposing regular ingoing boundary conditions at the horizon, which corresponds to computing a retarded dual correlator [80,81]: The functionB xz (u) is regular at the horizon. This choice of the boundary condition implies that near the horizon, δB xz behaves as By substituting this expression into the expectation value (67) of the two-form current J µν , we obtain The expression on the right-hand-side of the second equation is obtained by using the relation between δB xz and δb xz in (53). Note that the dependence on the electromagnetic coupling enters into the one-point function at order ω 2 and, thus,ᾱ plays no role in the holographic formula for the resistivity; r ⊥ and other first-order transport coefficients are independent of the renormalised electromagnetic coupling. Finally, using the Kubo formula for r ⊥ , which is derived and presented in Eq. (A5) of Appendix A, we recover the expression presented in Eq. (88). All of the six remaining transport coefficients can be obtained by following the same procedure. We refer the reader to Appendix B for their detailed derivations.
The plots of the (dimensionless) transport coefficients η , ζ , r ⊥ and r as a function of T / √ B are presented in Figure 2. The remaining three viscosities can easily be inferred from Eq. (88). In particular, η ⊥ /s = 1/(4π), ζ ⊥ = ζ /4 and ζ × = −ζ /2. We note that all transport coefficients satisfy the positive entropy production bounds discussed in Section I. It is interesting that the bulk viscosity inequality ζ ⊥ ζ ≥ ζ 2 × is saturated, i.e. ζ ⊥ ζ = ζ 2 × in the plasma studied here for all parameters of the theory.
The plots of (dimensionless) first-order transport coefficients as a function of T / √ B.
We can now investigate the behaviour of the transport coefficients in the two extreme limits of T / √ B → 0 and T / √ B → ∞, i.e. the strong-and the weak-field regimes, respectively. The leading-order power-law scaling (which we assume) and the coefficient follow from numerical fits. The results are presented in Table II. Since the entropy density s vanishes in the limit of zero temperature, all first-order transport coefficients vanish in the strong-field limit of T → 0. Furthermore, as we will see, all (first-order) dissipative effects also vanish in the T → 0 limit. These observations are consistent with predictions of [24] based on symmetry arguments.
In the regime of a weak magnetic field, T √ B, we find that both shear viscosities η ⊥ and η converge to η ⊥ = η = s/(4π) as B/T 2 → 0. On the other hand, the longitudinal bulk viscosity limits to ζ → 4η/3, which is consistent with the fact that as B/T 2 → 0, the evolution of the plasma should be governed by uncharged relativistic conformal hydrodynamics (see e.g. [25] or Appendix B). Indeed, both resistivities, r ⊥ and r , also tend to zero in the limit. We also note that the weak-field behaviour of r ⊥ and r is consistent with the assumption used to construct standard (ideal) MHD, whereby conductivity is taken to infinity, σ ≈ 1/r → ∞, and whereby one adds corrections proportional to 1/σ. 18 In other words, small weak-field resistivities are compatible with the assumption of ideal Ohm's law, which gives rise to Eq. (7) (see also our discussion around this equation in Section I.). Furthermore, note that in standard MHD, only one resistivity (conductivity) is typically added to include dissipative corrections. What we see is that in our theory, the two resistivities take similar values in the weak-field limit in which standard MHD applies. However, in the strong-field limit, they assume drastically different values, including a different scaling with T / √ B. This observation therefore further points to the important role of anisotropic effects in MHD [24] and the necessity for using the formulation of [24,27] as one moves from the weak-to the strong-field regime.
The fact that r ⊥ and r tend to zero both in the limits of T / √ B → 0 and T √ B → ∞, along with the positivity of the entropy production bounds r ⊥ ≥ 0 and r ≥ 0 [24], implies that there always exists a maximum value of the resistivities at some intermediate T / √ B. It would be interesting to find the sizes of these maxima in experimentally realisable systems and probe the regimes of the "least conductive" plasmas. Finally, it would be interesting to further investigate the connection between maximal r and various discussions of lower bounds on conductivities, e.g. [82][83][84].

IV. MAGNETOHYDRODYNAMIC WAVES IN A STRONGLY COUPLED PLASMA
We are now ready to use the information obtained from the holographic analysis of Section III to study dissipative dispersion relations of magnetohydrodynamic waves in a toy model of a strongly coupled plasma. We will use the theory of MHD [24], which is a phenomenological effective theory, and supplement it with microscopic details-the equation of state and transport coefficients-of the holographic setup investigated above. We will be particularly interested in the dependence of the MHD modes on the angle between momentum and magnetic field, as well as the ratio between temperature and the strength of the magnetic field. The 't Hooft coupling of interactions in the matter sector is not tuneable in our model, however, the electromagnetic coupling is. In all sections, except in Section IV C, it will be set to α = 2π 2 /137N 2 c . Before presenting the numerical result, we review the relevant facts about MHD modes. For a detailed derivation of these results, see Ref. [24] and for a discussion of the general procedure, see Refs. [25,85]. First, we write the hydrodynamic variables u µ , h µ , T and µ in terms of oscillating modes perturbed around their near-equilibrium values, e.g. u µ → (1, 0, 0, 0) + δu µ e −iωt+ikx sin θ+ikz cos θ , so that θ ∈ [0, π/2] measures the angle between the equilibrium magnetic field pointing in the z-direction and the wave momentum k in the x-z plane. The dispersion relations ω(k) are then derived from the equations of MHD, i.e. Eqs. (8) and (9), with the external H µνρ = 0. The solutions depend on the angle θ, temperature T and the strength of the magnetic field (or the chemical potential of the magnetic flux number density), parametrised in our solutions by B. Any dimensionless quantity will only depend on the single dimensionless ratio T / √ B. The resulting modes can be decomposed into two channels-odd and even under the reflection of y → −y. The first channel is the transverse Alfvén channel. The second is the magnetosonic channel with two branches of solutions: slow and fast magnetosonic waves.
The linearised MHD equations of motion (8) and (9) need to be expanded in the hydrodynamic regime in powers of small ω/Λ h 1 and k/Λ h 1, where Λ h is the UV cut-off of the effective theory. In standard MHD, where T √ B, then Λ h ≈ T , whereas in the strong-field regime of T √ B, the cut-off can be set by the magnetic field, then Λ h ≈ √ B. As argued in [24], hydrodynamics may exist all the way to T → 0, even when δT = 0. Such an expansion, performed to some order, gives rise to a polynomial equation in ω and k. For example, in the Alfvén channel, within first-order dissipative MHD, r ⊥ cos 2 θ + 2r sin 2 θ η ⊥ sin 2 θ + η cos 2 θ k 4 = 0 . (98) The two solutions of the quadratic equation for ω are given by where D A,+ and D A,− are One can now series expand ω(k) = D 0 k + D 1 k 2 , or alternatively, plug this ansatz in Eq. (98) and solve it order-by-oder in k. What we find is the Alfvén wave dispersion relation [24]: where the speed is given by V 2 A = µρ/(ε + p). The dispersion relation appears to be well-defined for any angle θ ∈ [0, π/2] between momentum and equilibrium magnetic field. In particular, if we were to take the θ → π/2 limit, (101) would yield two diffusive modes, both with dispersion relation which are, however, unphysical and only result from an incorrect order of limits of k and θ.
As can be seen from the structure of the square-root in Eq. (99), the expansion in small k is only sensible so long as k 2 V 2 A cos 2 θ/(D A,− ) 2 . Hence, even for a small finite k, this expansion is inapplicable for angles θ near θ = π/2 where cos θ becomes very small. In fact, for the propagating modes cease to exist altogether and the two modes become purely imaginary (diffusive to O(k 2 )). The transmutation of two propagating Alfvén modes into two non-propagating modes occurs when the inequality in (103) is saturated, i.e. at the critical angle θ c when Re[ω] = 0: In other words, the plasma exhibits propagating (sound) modes for 0 ≤ θ < θ c and non-propagating (diffusive) modes for θ c < θ ≤ π/2. We plot the dependence of the critical angle θ c on k/ √ B and T / √ B for the Alfvén waves in our model in Figure 3. What we see is that for small k/ √ B and small T / √ B, the transition to diffusive modes occurs closer to θ c ≈ π/2. For any fixed and finite T / √ B, Eq. (104) indeed implies that θ c → π/2 as k → 0.
We note that as already pointed out in [27], the limits of k → 0 and θ → π/2 do not commute and we obtain different results depending on which expansion (k ≈ 0 or θ ≈ π/2) is performed first. If one first takes the limit θ → π/2, then Eq. (98) becomes which instead of Eq. (102) results in two non-degenerate diffusive modes The dispersion relation (101) is therefore only sensible at a finite T / √ B and infinitesimally small k/Λ h .
In the magnetosonic channel, the story is entirely analogous to the one described for the Alfvén waves. By expanding around k ≈ 0 first, we obtain the dispersion relation of [24]: where the speed of magnetosonic wave is given by The functions V A , V 0 , V S and V appearing in (108) are The susceptibilities are 19 The two types of magnetosonic waves, corresponding to ± solutions in (108), are known as the fast (with +) and the slow (with −) magnetosonic waves. We refer the reader to Appendix C for further details regarding the derivation of the magnetosonic modes. Each pair of the propagating slow magnetosonic modes also splits, in analogy with the Alfvén waves, into two non-propagating diffusive modes for θ ≥ θ c . The critical angle θ c for magnetosonic modes is also defined as in the Alfvén channel: the angle at which Re[ω] = 0. We plot the numerically-computed dependence of the magnetosonic θ c on k/ √ B and T / √ B in Fig. 3. As can be seen from the plot, the critical angles for the two types of waves are independent. However, they show similar qualitative dependence on the parameters that characterise the waves. We summarise the θ-dependent characteristics of MHD modes in Fig. 4. We observe the pattern of a transmutation of sound modes into diffusion to be different in the weak-and strongfield regimes. Namely, the two magnetosonic waves interchange their dispersion relations at small θ. Since the complicated expressions for dispersion relations greatly simplify at θ = 0 and θ = π/2, we state them below. The sound mode dispersion relations, denoted by S, are S1 : and the diffusive modes, denoted by D, are   D1 : In the regime of a large T / √ B, the results agree with those of [27]. Furthermore, using the asymptotic form of the thermodynamics quantities and transport coefficients in the T / √ B → ∞ limit, one can show that these modes reduce to sound and diffusive modes of uncharged relativistic hydrodynamics.
In the strong-field regime, which cannot be described within standard MHD, the speeds of S1 and S3 become large and approach the speed of light in the limit of T → 0. It is clear that in the strong-field regime, MHD sound waves can easily violate any causal upper bound on the speed of sound [86][87][88][89]. Furthermore, as discussed above, all diffusion constants vanish and the system becomes controlled by second-order MHD [24], which we do not investigate in this work. All details regarding angle-dependent wave propagation are presented in Section IV B.

A. Speeds and attenuations of MHD waves
Here, we plot the speeds (phase velocities) and first-order attenuation coefficients of the three types of MHD sound waves: the Alfvén and the fast and slow magnetosonic waves for the holographic strongly coupled plasma discussed above. These results assume an infinitesimally small value of momentum k, and follow from first expanding the polynomial equation of the type of (98) around k ≈ 0 and writing each dispersion relation as ω = ±vk − iDk 2 . The speeds v (presented in Fig. 5) and attenuation coefficients D (presented in Fig. 6) are then plotted for all 0 ≤ θ ≤ π/2, which, as discussed above, is only physically sensible when θ c → π/2, i.e. as k → 0. The angular profiles of the speeds and the dissipative attenuation coefficients show distinct behaviour in the strong-, the crossover (cf. Eq. (86)) and the weak-field regimes. In particular, the speeds of sound enter the weak-field regime, where they reduce to well-known standard MHD results, rapidly after the temperature exceeds T / √ B ≈ 0.7. There, Alfvén and slow magnetosonic waves travel with very similar speeds for all θ and their speeds coincide at θ = 0 and θ = π/2. The situation is different in the strong-field regime where the profiles of speeds qualitatively match the strong-field predictions of [24]. There, slow magnetosonic and Alfvén waves can travel faster at small θ, with speeds comparable to those of fast magnetosonic waves. At θ = 0, the Alfvén speed equals that of fast, instead of slow, magnetosonic waves (cf. Fig. 4). It should also be noted that there exists a value of T / √ B in the crossover regimes where all three speeds are equal at θ = 0.
The attenuation coefficients, computed with all seven transport coefficients [24,27], are computed for the first time for a concrete microscopically (holographically) realisable plasma and therefore difficult to compare with other past results. What we observe is that the Alfvén waves experience the strongest damping for all values of T / √ B. Beyond that, the qualitative behaviour again displays distinct angle-dependent features in the three regimes, which are apparent from Fig. 6. A noteworthy, but not a surprising fact is that the strength of attenuation appears to be much more strongly dependent on the angle between momentum and magnetic field in the regime of small T / √ B. Furthermore, in the crossover regime, we find that the strengths of fast and slow magnetosonic mode attenuations interchange roles as T / √ B increases. In plots at T / √ B = 0.5 and T / √ B = 0.66, there exists an angle θ at which the two attenuation strengths coincide.

B. MHD modes on a complex frequency plane
By assuming a finite value of momentum k, a full analysis of the spectrum requires us to take into account the transmutation of sound modes into non-propagating diffusive modes. The pattern of this behaviour, as a function of the angle between momentum and the direction of the equilibrium magnetic field θ, was summarised in Fig. 4. Motivated by holographic quasinormal mode (poles of two-point correlators) analyses, we plot the motion of the MHD modes on the complex frequency plane-here, as a function of θ and T / √ B. One should consider these plots as a prediction of how the first-order approximation to the hydrodynamic sector of the full quasinormal spectrum computed from the theory (42) is expected to behave. In Fig. 7, we plot the typical θ-dependent trajectories of ω(θ) for Alfvén and magnetosonic modes in distinctly strong-and weak-field regimes. At all temperatures (except at T = 0 where D = 0), the behaviour is consistent with our previous discussions, including the fact that the transmutation of Alfvén and slow magnetosonic waves into diffusive modes occurs at lower θ c as k/ √ B increases. In the crossover temperature regime (around T / √ B ≈ 0.6), we can observe in more detail the interplay between fast and slow magnetosonic modes, which was noted in Section IV A. While the speed of fast magnetosonic waves always exceeds that of slow waves, their attenuation strengths exchange roles around T / √ B ≈ 0.675, which manifests in a characteristically distinct behaviour for θ < θ c , presented in Fig. 8 (see also Fig. 6). The θ-dependence of Alfvén waves remains qualitatively similar to those depicted in Fig. 7.
For a fixed θ < θ c , where θ c depends on k and T / √ B, we plot the typical behaviour of ω(k) as a function of T / √ B in Fig. 9. At T = 0, all poles start from the non-dissipative regime (the real w axis), with the speed of fast magnetosonic waves given by v = 1. As they move towards larger T / √ B, the Alfvén and the slow magnetosonic modes again asymptote to each other, eventually transforming into diffusive modes, while the speed of the fast magnetosonic modes gradually converges towards that of neutral conformal sound with v = 1/ √ 3.
In the high temperature limit, the "collision" of the Alfvén and, independently, the slow magnetosonic poles on the imaginary axis occurs close to the real axis, which follows from the fact that for both types of waves, as T / √ B → ∞. The Alfvén waves then become the diffusive modes of uncharged conformal hydrodynamics with ω = −iηk 2 /(2sT ). As for our final plot, in Fig. 10, we present the dependence of the four diffusion constants and one sound attenuation coefficient on the temperature at θ = π/2 (cf. Fig. 4 and Eqs. (111)-(112)). The modes D1, D3 and S1 reduce to dispersion relations of uncharged relativistic hydrodynamics. D2 and D4 are new.  10. Plots of the four diffusion constants (D1, D2, D3, D4) and the sound attenuation (S1) as a function of T / √ B at θ = π/2. Black, red and blue curves depict dissipative coefficients that originate from the Alfvén, slow magnetosonic and fast magnetosonic waves, respectively.

C. Electric charge dependence
We end our discussion of MHD dispersion relations by investigating their dependence on the choice of the U (1) coupling constant, or equivalently, the position of the Landau pole, which has so far been set to the (N c -rescaled)ᾱ = 1/137. All dependence onᾱ enters into the expectation value of the stress-energy tensor through the term proportional to H µν H µν ln C (cf. Eq. (66)), which contributes no terms linear in ω. For this reason, while the equation of state strongly depends on α, the first-order transport coefficients do not. Hence, all speeds of sound and attenuation (and diffusive) coefficients depend on the choice ofᾱ through the equation of state and susceptibilities.
What we observe is that the speeds of waves and attenuation coefficients strongly depend on the renormalised electromagnetic coupling, so, unsurprisingly, the strength of electromagnetic interactions plays an important role in the phenomenology of MHD. For concreteness, we only present the detailed behaviour of the Alfvén waves (with speed V A cos θ), which reduce to the neutral hydrodynamic diffusive mode D3 (and D4) at θ = π/2. Both V A and the diffusion constant of D3, D D3 , strongly depend onᾱ. For a small variation in the values ofᾱ, we plot the results in Fig. 11. 20 To show the importance of a sensible choice of the renormalisation condition, we also vary the coupling over a larger range (toᾱ = 80/137), where we see that the system develops unphysical behaviour with instabilities. As is apparent from Fig. 12, Alfvén waves become unstable at low T / √ B.  Fig. 4), i.e. V 2 0 , which is plotted for comparison.
In all to us known literature, the unavoidable choice of the constant C, which setsᾱ, is made in a different way. C is either chosen so that the logarithmic terms vanish altogether, or so that it sets the UV scale to that of the magnetic field, which is convenient when studying strong magnetic fields as e.g. in [38,47]. Here, we wish to point out some of the consequences of setting C to either of the two standard options. The first option, which eliminates the logarithmic terms, results in the following thermodynamics quantities: 20 We remind the reader that in the boundary Lagrangian, the electromagnetic coupling is scaled out from the covariant derivatives. Thus, only the Maxwell term depends on er. As we vary er, we keep the strength of the electromagnetic field fixed. 12. The plot of the Alfvén V 2 A at a varyingᾱ ranging fromᾱ =ᾱ 0 toᾱ = 80ᾱ 0 , whereᾱ 0 = 1/137. We see that asᾱ increases, the waves develop an instability in the strong-field regime.
The second choice results in While these two renormalisation conditions are suitable for studying certain physical setups involving static electromagnetic fields, we claim that they lead to unphysical results when the boundary U (1) gauge field is dynamical. By comparing the renormalised stress-energy tensor (78)-(80) to expressions in (114) and (115), we find that the two choices correspond to the renormalised coupling being e 2 r → ∞ and e 2 r ∼ ln B, respectively. An infinite U (1) coupling is unphysical in a plasma state. The problem with the second choice is that if extrapolated to the weak-field regime, ln B/M , where M is some scale, can become negative and e r imaginary, which is again unphysical. Thus, these choices may lead to instabilities and superluminal propagation, which were absent from our results withᾱ near 1/137. We plot the Alfvén speed parameter V A for the two couplings from (114) and (115) in Fig. 13.

V. DISCUSSION
This work should be considered as the first holographic step in a long road towards a better understanding of magnetohydrodynamics in plasmas outside of the regime of validity of standard MHD, be it in the presence of strong magnetic fields or in a strongly interacting (or dense) plasma with a complicated equation of state and transport coefficients-all claimed to be describable within the recent (generalised global) symmetry-based formulation of MHD of Ref. [24]. In order to supply a hydrodynamical theory of MHD with the necessary microscopic information of a strongly coupled plasma, we resorted to the simplest, albeit experimentally inaccessible option: holography. Nevertheless, our hope is that in analogy with the myriad of works on holographic conformal hydrodynamics, which have led to important new insights into strongly interacting realistic fluids, holography can also help us understand observable MHD states in the presence of strong fields, high density and of strongly interacting gauge theories, such as QCD.
With this view, we constructed the simplest theory dual to the operator structure and Ward identities used in MHD of [24], investigated the relevant aspects of the holographic dictionary and used it to compute the equation of state and transport coefficients of the dual plasma state. This information was then used to analyse the dependence of MHD waves-Alfvén and magnetosonic waves-on tuneable parameters specifying the state: the strength of the magnetic field, temperature, the angle between momentum of propagation and the equilibrium magnetic field direction, as well as the strength of the U (1) electromagnetic gauge coupling. We believe that the latter feature of our model-dynamical electromagnetism on the boundary-which in the (dual) language of two-form gauge fields in the bulk allows for standard (Dirichet) quantisation, could in its own right be used for holographic studies of U (1)-gauged systems, unrelated to MHD.
Our results have revealed several new qualitative features of MHD waves, particularly in the regime of a strong magnetic field, which is inaccessible to standard MHD methods. Various properties of the equation of state, transport coefficients and dispersion relations found here, may now be compared to those in experimentally realisable plasmas, or at the least, used as a toy model for future studies of MHD. Approximate scalings in the limiting regimes of large and small T / √ B are collected in Tables I and II. Here, we summarise some of the most interesting observations: • The equation of state and transport coefficients strongly depend on the strength of the magnetic field, i.e. on whether the plasma is in the weak-field, the crossover, or the strong-field regime.
• In the weak-field regime with T / √ B 1, the system is well-described by standard MHD (see [27] for a full description) with small resistivities (large conductivity regime, which is assumed by ideal Ohm's law) and small effects of anisotropy. As T / √ B → ∞, the plasma becomes an uncharged, conformal fluid with a single independent transport coefficient, η = s/4π. In the strong-field limit of T / √ B → 0, the plasma limits to a non-dissipative regime with all first-order transport coefficients (along with sound attenuations and diffusion constants) tending to zero. Effects of anisotropy are large.
• Resistivities have a global maximum in the intermediate T / √ B regime, which indicates a regime of least conductive plasma. If the assumptions of standard MHD are correct at T / √ B 1 and the symmetry-based predictions of [24] are correct at T / √ B 1, such a regime should be generically exhibited by any plasma.
• Out of the three bulk viscosities, ζ ⊥ , ζ and ζ × , only one is independent and they saturate the positivity of the entropy production inequality, i.e. they are related by ζ ⊥ ζ = ζ 2 × . One may speculate on how general this result is and whether it is related to the suppression of entropy production at strong coupling [90,91] or perhaps some form of (holographic) universality at infinite (or strong) coupling.
• Various qualitative features of slow and fast magnetosonic modes are exchanged in the weakand strong-field regimes (usually at small angles, θ, between momentum and the equilibrium magnetic field direction), such as their asymptotic tendency to the speed of Alfvén waves and the strength of sound attenuation.
• For a finite momentum, propagating Alfvén and slow magnetosonic modes (sound modes to O(k 2 )) transmute into pairs of non-propagating, diffusive (to O(k 2 )) modes. This occurs at large angles between the direction of momentum propagation and the equilibrium magnetic field, θ c < θ ≤ π/2, where θ c is some momentum-and T / √ B-dependent critical angle (cf. Eq. (104) for Alfvén waves).
• The phenomenology of MHD modes strongly depends on the strength of the electromagnetic coupling (or the position of the Landau pole) and can, for large ranges of the coupling, lead to unstable or superluminal propagation.
Beyond the types of waves studied in this work, it would be particularly interesting to better understand the role of finite charge density, as studied in [27], within the formalism of [24]. The important question then is how the phenomenology of such MHD waves, which typically experience gapped propagation and instabilities (e.g. the infamous Weibel instability), becomes altered by strong interactions, strong fields and for more 'exotic' field content.
Finally, the holographic setup studied here will need to undergo extensive further tests and analyses in order to unambiguously establish its connection to plasma physics and MHD. In particular, it is essential to study the quasinormal spectrum of the theory to verify that the hydrodynamic modes indeed describe the small-ω and small-k expansion of the leading infrared poles. Furthermore, it will be interesting to understand the role of higher-frequency spectrum and its interplay with MHD modes. We leave all these and many other interesting questions to the future. publication. S. G. is supported in part by a VICI grant of the Netherlands Organisation for Scientific Research (NWO), and by the Netherlands Organisation for Scientific Research/Ministry of Science and Education (NWO/OCW). The work of N. P. is supported by the DPST scholarship from the Thai government and by Leiden University.
Here, we show the details of the derivation of horizon formulae of all remaining transport coefficient: η ⊥ , η , ζ ⊥ , ζ , ζ × and r . The computation are analogous to the calculation of r ⊥ in Section III D. The solution to leading order in the frequency ω can be found analytically and its near-boundary expansion gives δG y x = δh xy 1 + where δh xy is the Dirichlet background condition and the boundary theory source. If we plug this solution into to the stress-energy tensor, we find that Using Eq. (A5), we find that as stated in Eq. (88).

(ii) Shear viscosity η
Similarly to the computation of r ⊥ , the xu-component of the two-form gauge field fluctuation equation can be used to reduced the two coupled second order differential equations coupling δG xz and δB tx to a single equation: The solution to linear order in ω can again be found analytically and in the near-boundary region yields The relevant component of the stress-energy tensor is then which gives as stated in Eq. (88).
(iii) Resistivity r The only equation of motion in this channel is which leads to the near-boundary solution δB xy = δB (0) The two-form current can then be written as which yields as stated in Eq. (88).
(iii) Bulk viscosities ζ ⊥ , ζ and ζ × By counting the number of the relevant degrees of freedom, it turns out that there is only one dynamical mode in this decoupled system coming from 4 × (2 nd -order ODE's for δg tt , δg aa , δg zz , δb tz ) − 3 × (1 st -order ODE's for δg tu , δg uu , δb zu ). To find the dynamical mode, we start by solving the algebraic equations for δg tu , δg uu and δb zu from the tu and uu components of Einstein's equation combined with the zu component of Maxwell's equations. Plugging these solutions into the four second-order equations involving δg tt , δg aa , δg zz and δb tz , we find that the remaining two nontrivial equations involve only δg aa and δg zz . The single resulting equation of motion can then be expressed in terms of the gauge-invariant variable Z s (u) defined as where δg aa = δg xx + δg yy . The equation of motion for Z s can be written Z s (u) + C 1 (ω, u)Z s (u) + C 2 (ω, u)Z s (u) = 0 , where Now, suppose that the time-independent solution for Z s is Z (−) , so that Z (−) (u → 0) = Z (0) ≡ δh aa − 2δh zz (note that V /W → 1 and u → 0). The second solution, denoted as Z (+) , contains the time-dependent information and can be found from the Wronskian We then find that the near-boundary and near-horizon expansion for Z (+) are The full solution is a linear combination, Z s (u) = Z (−) +αZ (+) , and the ingoing boundary condition sets the frequency-dependent function α(ω) to be which allows us to write the solution for Z s near the boundary as and ζ The coefficients A ij are A 11 = 1 2T 2 χ 11 (µ + T χ 12 )(µ − T µ 21 ), A 12 = 1 2T ρχ 11 (µ + T χ 12 )(µ cos 2 θ + ρχ 22 sin 2 θ), µ cos 2 θ + ρχ 22 sin 2 θ , A 31 = s + ρχ 21 ε + p , A 32 = 2ρ (ε + p) sin θ A 22 , A 33 = η cos 2 θ + (η ⊥ + ζ ⊥ ) sin 2 θ ε + p , A 44 = 2ζ cos 2 θ + η sin 2 θ sT . (C3) By computing the determinant in (C1), the resulting quartic equation is where c i are functions of thermodynamics quantities, transport coefficients, k and θ. The expressions for c i in terms of A ij in (C3) are 11 s cos 2 θ + T χ 11  (C5) In principle, Eq. (C4) gives a closed-form solution for the four ω(k). In practice, the explicit solutions are extremely lengthy so it is often more convenient to find the roots of (C4) numerically