Holography and magnetohydrodynamics with dynamical gauge fields

Within the framework of holography, the Einstein-Maxwell action with Dirichlet boundary conditions corresponds to a dual conformal field theory in presence of an external gauge field. Nevertheless, in many real-world applications, e.g., magnetohydrodynamics, plasma physics, superconductors, etc. dynamical gauge fields and Coulomb interactions are fundamental. In this work, we consider bottom-up holographic models at finite magnetic field and (free) charge density in presence of dynamical boundary gauge fields which are introduced using mixed boundary conditions. We numerically study the spectrum of the lowest quasi-normal modes and successfully compare the obtained results to magnetohydrodynamics theory in $2+1$ dimensions. Surprisingly, as far as the electromagnetic coupling is small enough, we find perfect agreement even in the large magnetic field limit. Our results prove that a holographic description of magnetohydrodynamics does not necessarily need higher-form bulk fields but can be consistently derived using mixed boundary conditions for standard gauge fields.


Introduction
The essence of the holographic duality lies in the correspondence between a gravitational theory in d + 1 dimensions and a "dual" field theory in d dimensions which is formally defined via the so-called GPKW (Gubser, Polyakov, Klebanov, Witten) master rule [1,2]. Importantly, the bulk gravitational action does not uniquely define the boundary dual field theory but, for that scope, needs to be supplemented with boundary conditions for all the bulk fields.
Let us consider the example of a bulk gauge field A µ . Formally, in the (d + 1)dimensional bulk, one usually writes a gravitational action of the form: where R is the Ricci scalar, Λ a negative cosmological constant and L (F µν ) a generic Lagrangian for the bulk field A µ which is written in terms of its field strength F := dA, as imposed by gauge invariance in the bulk. The simplest example possible is a Maxwell kinetic term, L (F µν ) = −F 2 /4 with F 2 := F µν F µν , which gives rise to the so-called Einstein-Maxwell action. The ellipsis in Eq.(1.1) indicates the presence of other possible bulk fields which are irrelevant for the present discussion and therefore not shown. Also, for simplicity, let us neglect possible couplings between the field strength F µν and other bulk fields, e.g., dilaton couplings. The asymptotic solution for a massless gauge field in an asymptotically Anti-de Sitter bulk geometry is then generically given by: where r is the radial holographic direction, r = ∞ the location of the AdS boundary. At this point, A (0),(1) µ (t, x) are two undetermined functions, usually denoted as the leading/subleading terms, which must be fixed by the choice of boundary conditions at r = ∞. The canonical (but not unique) procedure, which goes under the name of standard quantization, is to fix the value of A (0) µ at the boundary and dynamically determine the value of A (1) µ . From a dual field theory perspective, this corresponds to a CFT deformation of the type: where J µ is a U(1) current operator in the dual CFT with conformal dimension ∆ = d − 1 whose vacuum expectation value is determined by A (1) µ (see [3] for details). Gauge invariance in the bulk, A µ → A µ + ∂ µ ξ, implies that the dual U(1) current is conserved, ∂ µ J µ = 0. Additionally, the current operator two-point function JJ (whose spatial component is related to the electric conductivity σ) can be easily computed in linear response theory by looking at the ratio A (1) /A (0) . Within this picture, it is often said that a local symmetry in the bulk (in this case, a local U(1) gauge symmetry) corresponds to a global symmetry in the dual boundary theory. In other words, following the prescription just described and using standard quantization, the dual field theory does not possess a local U(1) symmetry and, as a consequence, no electromagnetism nor long-range Coulomb interactions are present therein.
Having in mind possible applications of the holographic duality to realistic systems, this outcome might appear rather disappointing and limiting. There are indeed several circumstances in which the role of dynamical gauge fields and electromagnetic interactions are fundamental for the correct physical description and cannot be neglected. Plasma physics and condensed matter systems are clear examples of this sort. Even if too often ignored, within the holographic framework, the solution to this problem is more than ten years old. Starting from the seminal works by Witten [1,4] (see also [5][6][7][8]), the existence and meaning of different boundary conditions for bulk vector fields were analyzed in detail by Marolf and Ross in 2006 [9] 1 and few years later applied for the first time in the context of holographic superconductors [11][12][13], to be contrasted with the more famous holographic superfluid model of Hartnoll, Herzog and Horowitz [14]. 2 Without going into details, the main idea of this program is to promote the b.c.s. for the gauge field at the AdS boundary to the most general form: which takes the name of mixed boundary conditions. 3 The standard quantization is recovered by setting β = 0, while the alternative quantization by setting α = 0. It was early realized (mostly for the simpler case of a bulk scalar field) [4,6,[17][18][19] that these most general b.c.s. are connected to double trace deformations in the boundary field theory and specific SL(2,R) transformations in the moduli space of the dual CFTs.
As already mentioned, one of the first concrete realizations of these modified boundary conditions in the context of applied holography has appeared in the construction of a "real" holographic superconductor model [12,13,20] which has turned out to be fundamental for the study of certain specific properties of superconductors such as vortices [11,[21][22][23][24][25][26][27][28][29] or the Meissner effect [30]. Similar types of mixed b.c.s. have been considered in the study of anyons physics in [31][32][33]. A more recent explosion of efforts to incorporate, understand and utilise the effects of dynamical electromagnetism in the boundary field theory is connected to the study of magnetohydrodynamics and plasmons physics. From one side, the latter has been initiated by Gran, Tornsö and Zingg [34] using mixed b.c.s. for the bulk gauge field and has been investigated in several directions [35][36][37][38][39][40][41][42]. The connection of this framework to the diagrammatic Random Phase Approximation (RPA) [43] and double trace deformations in the dual field theory has been explained in [44,45]. On the contrary, the aspects related to magnetohydrodynamics have been so far dominated by the usage of higher-form symmetry structures as proposed in the original work by Hofman, Grozdanov and Iqbal [46] which has been implemented within the holographic framework in [47][48][49] and discussed from a hydrodynamic perspective in [50][51][52][53]. In particular, Grozdanov and Poovuttikul [47] have demonstrated a complete match between the magnetohydrodynamic expectations [50] and the bulk higher-form picture. Moreover, in this class of theories, a bona fide photon has been identified [54] 4 and interestingly described as a Goldstone mode of the emergent higher-form symmetry [55].
Is the higher-form bulk picture really necessary to have dynamical electromagnetism at the boundary? What is its relation with the mixed boundary conditions discussed so far? Can we recover magnetohydrodynamics in the boundary field theory without using higherform symmetries? Most of the answers to these questions have been already addressed in a beautiful work by Higginbotham and DeWolfe [56] (see Fig. 1 therein for a nice summary of their results). In a nutshell, the Hodge dual operation performed in the bulk, and needed to pass from the standard Maxwell picture to the higher-form description (which in AdS 5 appears as a Maxwell action for a 3-form field strength), does not leave the boundary conditions unchanged. On the contrary, it hiddenly modifies the standard Dirichlet b.c.s. into mixed b.c.s. rendering the original global U(1) symmetry in the boundary dynamical. As a consequence, the higher-form formalism introduces non-trivial physics in the boundary field theory because it corresponds to deforming the original Dirichlet boundary conditions for the Maxwell bulk gauge field A µ . In other words, despite the bulk physics is unchanged because of the harmless Hodge dual, the field theory interpretation of the two scenarios is completely different. 5 Therefore, the naive expectation is that one could be able to obtain the same results, i.e., to obtain dynamical electromagnetism and magnetohydrodynamics in the boundary, by sticking to the maybe less elegant but more direct gauge field picture and deforming its asymptotic b.c.s. without the need of any higher-form structure. This possibility and its outcomes are the subject of this paper.
More concretely, we are asking whether a standard Einstein-Maxwell action: implemented with the "right" (and indeed not Dirichlet) boundary conditions for the gauge field A µ is able to provide the physics of a dual system exhibiting a dynamical U(1) symmetry with EM (electromagnetic) long-range interactions.
In the case of standard Dirichlet b.c.s., the dual field theory is a finite temperature CFT with a conserved U(1) current J µ in presence of an external, and not dynamical, gauge field A µ (and therefore an external magnetic field B as well). This scenario has been studied in several works [59][60][61][62][63] and the complete consistency between the holographic picture and the dual hydrodynamic framework has been recently verified in [62]. This same system has also been studied in presence of explicit and/or spontaneous breaking of translations [64][65][66][67][68][69][70] and anomalies [71]. Here, in analogy with the higher-symmetry analysis in [47], we aim at running a parallel program for the case in which the U(1) symmetry, and correspondingly the magnetic field B, in the boundary are dynamical. Our holographic results will be compared to the magnetohydrodynamics derived in Ref. [50] reduced, for simplicity, to two spatial dimensions. We will study the system at finite charge density and finite dynamical magnetic field and explore the regime of strong magnetic field.
Finally, we will discuss the more speculative possibility of modifying the nature of the dual field theory not by using boundary conditions nor by performing a Hodge duality in the bulk but rather by substituting the original Maxwell term F µν F µν in the bulk with a non-canonical higher derivative action of the form (F µν F µν ) N/2 .
Aware of the issues of "naturalness" in the effective field theory sense, we will consider this case as a toy model to understand better the implementation of symmetries and boundary conditions in bottom-up holography. The idea for an action as in (1.6) is borrowed from the holographic axions model [72] in which this type of higher order kinetic terms [73] has been employed to realize the spontaneous symmetry breaking of translations in the dual field theory [74]. Therein, this procedure turned out to be equivalent to modifying appropriately the boundary conditions for the axion fields responsible for the breaking of translations [58]. Despite its odd nature, the bulk action written in terms of non-canonical kinetic terms exactly reproduces the structure and dynamics of viscoelasticity theory [75] (see [76] for a review on the topic) proving its validity as an effective bulk description. Here, we will perform the same analysis for a bulk gauge field with non-canonical kinetic term as in Eq.(1.6). As we will explore in detail, the deformation of the bulk action as in (1.6) automatically modifies the nature of the coefficient A (0) µ (t, x) in the asymptotic expansion of the gauge field (1.2) from leading to subleading. This implies that, assuming standard quantization for the theory in (1.6), the coefficient A (0) µ (t, x) is not anymore a source for an external field A µ but rather the expectation value of the current J µ . This is exactly what would happen by considering the standard Maxwell action, N = 2, but with alternative b.c.s. for the bulk gauge field A µ . Indeed, the two frameworks will give analogous results.
Structure of the paper -In Section 2, we revisit the the magnetohydrodynamic framework of [50] in 2 + 1 dimensions and obtain the dispersion relation of the low-energy modes at finite charge density and dynamical magnetic field; in Section 3, we introduce our holographic setup and the precise boundary conditions used: in Section 4, we present all the main results of our work with modified mixed b.c.s. for the bulk gauge field A µ ; in Section 5, we discuss the features of a higher-derivative bulk model and its meaning from the point of view of the dual field theory; finally, in Section 6, we conclude and discuss a few points for future investigation. Appendix A discusses some interesting outcomes regarding the regime of validity of the hydrodynamic framework.

Hydrodynamics with dynamical gauge fields
In this section, we consider relativistic magnetohydrodynamics in (2+1) dimensions, including the effects of finite charge density and magnetic field. With the term relativistic magnetohydrodynamics we refer to the hydrodynamic description of a relativistic charged fluid in presence of long-range EM interactions mediated by a dynamical gauge field (see [77] for a recent review on the topic). This is very different from the situation (which is often loosely labelled in the same way) in which the magnetic and electric fields are external and non dynamical (see [78] for a review). From a practical perspective, the computations presented in this section are identical to those in [50] but in the simpler situation of (2+1) dimensions. The main simplification with respect to [50] arises from the fact that a magnetic field in two spatial dimensions cannot be associated to a proper vector field and therefore one cannot define an angle θ between the magnetic field B and the wave-vector k. This avoids several complications related to the anisotropy of the system in (3+1) dimensions. As a downside, the dynamics in (2 + 1) is less richer than that in (3 + 1) dimensions. For example, it does not include the so-called Alfvén waves nor the separation between fast and slow magnetosonic waves.

Setup
Let us start by considering the generating functional Z[g µν , A µ ]: where Φ denotes a set of dynamical fields, A µ an external gauge field coupled to the field theory current J µ (Φ) and g µν a fixed external metric coupled to the stress tensor T µν (Φ). Using Eq.(2.1), and the standard functional derivative prescription, the n-point functions of the corresponding conserved current operator J µ (and not only) can be computed. Within the holographic business, this would correspond to imposing the standard Dirichlet boundary conditions on the bulk gauge field (see more details below). From the generating functional in Eq.(2.1), we can define an effective action as where F denotes the free energy density. The derivative expansion of F, gives the thermodynamic pressure p at the leading order in fluctuations, i.e., at equilibrium. Here, T is the temperature, µ the chemical potential, and B the magnetic field. Moreover, we have assumed that T, µ, B are O(1) in derivatives while E is order O(∂). This is the correct assumption in case of magnetohydrodynamics (MHD) [50]. Using Eq.(2.2), one can further define the stress-energy tensor T µν and the U(1) conserved current J µ as So far, all the quantities discussed are defined in terms of the external fields g µν , A µ which are not dynamical. In order to promote the external gauge field to be dynamical, one considers the following Legendre-transformed action (2.5) In the second line of Eq.(2.5) we separate S into two pieces which correspond to the "matter contribution" S m , and the Maxwell kinetic term for the dynamical gauge field. The Maxwell kinetic term is defined using the field strength F := dA. The last term in Eq.(2.5) represents a coupling of the dynamical gauge field to an external current J µ ext . Here, λ is the square of the electromagnetic coupling. For convenience, we define the matter contribution to the stress-energy tensor T µν m and to the U(1) conserved current J µ m as In order to show the physical meaning of the action in Eq.(2.5), it is convenient to vary it with respect to A µ obtaining The vanishing of Eq.(2.7) corresponds to the standard Maxwell equations, i.e., implying that the gauge field A µ is now dynamical and coupled to the external current J µ ext as in standard electromagnetism through the coupling λ. Using the action in Eq.(2.5), the dynamical equations of motion can be summarized as where the Levi-Civita symbol is taken following the notation 012 = 1/ √ −g, and T µν EM is the stress-energy tensor of the the Maxwell kinetic term given by Likewise, we can define the contribution of the Maxwell kinetic term to the current as: With these notations, the total stress tensor and current are given by: 13) and, in absence of external sources, J ext = 0, are both conserved. In terms of the total stress tensor and total current, the EOMs in (2.9)-(2.10) become simply: as reported in [50]. Notice that the first equation in (2.10) implies the independent conservation of the external current J µ ext as well. The first two equations (2.9), can be obtained by utilising the diffeomorphism invariance (and the gauge invariance) of the action in Eq.(2.5). The other equations (2.10), are the Maxwell equation and the electromagnetic Bianchi identity, respectively. Note that the Maxwell equation determines the evolution of the dynamical gauge field A µ and the Bianchi identity is used to ensure that the electric/magnetic fields can be derived from a scalar/vector potential.
In order to solve the equations of motion (2.9)-(2.10), we need to further specify the constitutive relations for either T µν m and J µ m , or equivalently T µν and J µ . Following the standard procedure to construct hydrodynamic theories, for this purpose, we will use a gradient expansion (more details about this procedure and its validity will be provided below). In the Landau (or energy) frame, the constitutive relations at first order in derivatives are where is the energy density, ρ is the charge density, p is the pressure given in (2.3), and ∆ µν := g µν + u µ u ν the projection tensor in terms of the fluid velocity u µ . In addition, we have defined the H tensor: where M µν m is the polarization tensor (2.19). Finally, (Π µν , ν µ ) are the first order in derivatives dissipative corrections (2.25). Then the constitutive relation of matter part of the stress-energy tensor and U (1) current density are given by where m = − 1 4λ F 2 , p m = p + 1 4λ F 2 , and ρ m = ρ. All the quantities indicated with a subindex m relates to the matter part of the total action in Eq.(2.5). Using the constitutive relation for the current J µ m in (2.17), one can notice that Eq.(2.8) corresponds to the standard Maxwell equation in matter: in which J µ free := ρ m u µ + ν µ and J µ bound := ∇ ν M µν m . J µ free refers to the current of free charges while J µ bound incorporates the polarization effects. We can decompose the polarization tensor M µν m and H µν with respect to fluid velocity as and can also be identified with M µν m = 2∂p m /∂F µν , H µν = 2∂p/∂F µν . The objects D µ and H are respectively the displacement vector and the magnetic H-field, also known as the magnetic field strength [79]. In (3+1) dimensions [50], the magnetization M in Eq. (2.19) becomes the magnetic polarization vector M µ . The electric polarization vector P µ and the magnetization M are associated with the electric field E µ and magnetic field B via the susceptibilities (χ EE , χ BB ), i.e., with To understand better the physical meaning of D µ and H, it is convenient to re-write Eq. (2.18) Eq.(2.20) implies that D µ and H are also proportional to the electric and magnetic field E µ and B via the following relations in which we have defined the electric permittivity e and the magnetic permeability µ m . Using all the previous identities and definitions, we finally arrive at the following identities which connect the susceptibilities to the electric permittivity and the magnetic permeability.
Continuing with the hydrodynamic construction, the dissipative terms Π µν and ν µ are given by where η is the shear viscosity and σ the conductivity. Here, conformal symmetry has been assumed (e.g., the bulk viscosity is not appearing therein).
One can now solve the equations of motion (2.9)-(2.10) together with the constitutive relations (2.15) and obtain the low-energy excitations of the system. For this purpose, we consider the following set of fluctuations (δT, δµ, δu i=x,y , δE i=x,y , δB) around the equilibrium configuration with (T eq , µ eq , B eq ) being the equilibrium values for the corresponding thermodynamic quantities. In order to simplify the notations, we will drop the subfix eq in the following. Note that the fluctuations of the electric field (δE i ) and magnetic field (δB) appear explicitly in this case as a direct manifestation of the dynamical, rather than external, gauge field. Indeed, this is one of the major differences with the hydrodynamics with external gauge fields considered for example in [59][60][61][62].
In Fourier space, assuming the spacetime dependence of the fluctuations to be proportional to e −iωt+ikx , the linearised equations of motion can be rewritten in matrix form as where M(ω, k) is a 6 × 6 matrix and s A = {δT, δu i=x,y , δE i=x,y , δB}. The matrix M(ω, k) is a function of the thermodynamical susceptibilities 28) and the various thermodynamical parameters ( , p , T , µ , ρ , B , σ , η , e , µ m ). Then, the dispersion relations (or eigenmodes), ω = ω(k), can be obtained by solving the condition det M(ω, k) = 0 . (2.29) The complete analytic expressions of Eq.(2.29) and of the dynamical matrix M(ω, k) itself are rather lengthy and therefore not made explicit. We will show the dispersion relations of the low-energy modes only in the low ω, k expansion in the next sections. Furthermore, in the main text, to simplify our formulas and avoid clutter, we will only show the formulas in the low B expansion. The small B expansion is sometimes called the "weak field" limit and it is somehow equivalent to assuming the magnetic field to be of order O(∂) [50]. When expanding the expressions in this limit, all the thermodynamic quantities, such as , p, will be independent of B. To avoid clutter, we will still denote them in the same way, avoiding additional subscripts. Readers interested in the complete expressions are referred to the GitHub repository available here.

Zero density
We first study the hydrodynamics of a neutral plasma, at zero charge density, ρ = 0, or equivalently at zero chemical potential, µ = 0. Moreover, we may further set the following thermodynamic susceptibilities for the neutral state. This assumption will be verified a posteriori using the holographic computations in the next section. Using Eq.(2.30) together with ρ = 0, the spectrum of low-energy excitations exhibits six modes: four gapless modes and two gapped modes. Let us first discuss the simplest case, i.e., the case with a vanishing magnetic field, B = 0. For such a system, we have a pair of longitudinal sound waves together with one transverse shear diffusion mode where v 2 s = ∂p/∂ and Γ s = η/( + p). In this case, ( +p) is just the momentum susceptibility. The modes above are decoupled from the others and simply follow from the conservation of energy and momentum as in simple relativistic fluids. In addition, the remaining gapless mode can be determined from the "telegrapher equation" for electromagnetic (EM) waves [80,81]: where e is the electric permittivity, µ m magnetic permeability. In vacuum, σ = 0, Eq.(2.32) would just give rise to standard electromagnetic waves with light speed equal to c 2 ≡ 1/ ( e µ m ). Screening effects are introduced when σ = 0 and relate to the well-known skin effect (for which, differently from here, one usually assumes ω ∈ R and k ∈ C). Notice, that in our holographic setup, which does not preserve Galilean symmetry, the conductivity σ is finite even at zero charge density.
Eq.(2.32) gives rise to the following solution: where we have identified c 2 = 1/( e µ m ). The dynamics in Eq.(2.33) is sometimes labelled as k-gap [81] since the dispersion relation of the modes acquires a finite real part only above a certain critical wave-vector. This happens via a collision between a diffusive hydrodynamic mode and a relaxational non-hydrodynamic mode. As we will prove explicitly, Eq.(2.33) is an accurate description only when the k-gap is small (in this case σ/(2 c e ) T ), or also in the so-called quasihydrodynamic regime [51].
In the small wave-vector limit, Eq.(2.32) gives rise to two modes with dispersion EM waves do not propagate anymore at long distances (small k). Conversely, for σ/ e = 0, the magnetic field diffuses with a diffusion constant D B := 1/(σµ m ) while the electric field relaxes with a rate τ −1 e := σ/ e . In the language of global higher-form symmetries [46], the electric U(1) symmetry is explicitly broken while the magnetic one is preserved implying the conservation of the magnetic flux and the presence of magnetic diffusion.
The last low energy mode is a longitudinal damped charge diffusion mode Here, we have defined the charge susceptibility χ ρρ = ∂ρ/∂µ. Charge fluctuations do not diffuse anymore but they are rather relaxing with a rate equal to τ −1 e , and identical to that for the electric field E.
Next, let us turn on the magnetic field and discuss its effects on the low energy modes. We focus on the small B limit, defined as B/T 2 1 (where the eq subscripts are neglected for simplicity), and use the following identities: which are valid in that regime. The longitudinal sound waves and transverse shear mode in (2.31) are now modified into (2.37) The sound modes in (2.37) are known as magnetosonic waves. Their velocity v ms and attenuation constant Γ ms are given by , (2.38) The EM waves still follow the same dynamics as in Eq.(2.32) but their dispersion relations at small wave-vector are corrected by the presence of a finite magnetic field. In particular, the expressions in Eq.(2.34) are now modified into Finally, the dispersion relation of the damped charge diffusion mode (2.35) becomes where, once again, the limit of small magnetic field is assumed. We remind the Reader that in this limit, all the thermodynamic quantities appearing in the expressions above are not functions of the magnetic field B. This is equivalent to the "weak field" limit in [50]. Note that all the dispersion relations presented in (2.37)-(2.40) are consistent with the ones in (3+1) dimensions derived in [50]. A main qualitative difference is that in (2+1), one does not have Alfvén waves because the magnetic field is always perpendicular to the wavevector. 6 Indeed, Alfvén waves disappear even in (3+1) dimensions when the direction of propagation is perpendicular to the magnetic field [50] (θ = π/2 in their notations, see also [47]). Moreover, the distinction between fast and slow magnetosonic waves, which relies on the presence of a finite angle θ, is also absent in (2+1).

Finite density
Next, we study the hydrodynamics at finite density (ρ = 0). To be precise, by "finite density" we mean a finite density of free charges ρ = J t . Because of Maxwell's equations, we necessarily have J µ + J µ ext = 0 at equilibrium. This implies that we also have J t ext = −ρ and that the total charge density is zero, J t tot ≡ J t + J t ext = 0. Physically, we should think of this situation as a system with a finite density of free charges (e.g., electrons) together with an equal finite density of ions which render the total system neutral. In [50], this state is labelled as "charged state offset by background charge". For simplicity, we will continue with the simpler notation "finite density", well aware of the caveat discussed above.
Solving Eq.(2.29) in the small wave-vector limit, one finds six low energy modes corresponding this time to two gapless modes and four gapped modes. Note that the number of the gapless modes (and gapped modes) is different from the zero density case. This is related to the fact that the small ρ limit does not commute with the small k limit when the gauge field is dynamical [50].
We now show the structure of the low-energy modes at finite charge density and finite magnetic field B, focusing on the small B limit. We find one longitudinal diffusive mode and one transverse subdiffusive shear mode with dispersions The fact that the dispersion relation of the subdiffusive mode is not well defined at ρ → 0 is just a manifestation of the non-commutativity of limits. Note that the gapless modes in (2.41) appear also in the case of external gauge fields [62]. Interestingly, the dispersion of the diffusive mode in (2.41) is exactly the same (at this order in k) as the one in the case of external gauge fields (see Eq.(2.51)). This observation is also consistent with the results in the higher dimensional (3+1) theory [50]. On the contrary, the subdiffusive mode in (2.41) displays a slightly different form (see Eq.(2.51) for the subdiffusive mode in the case of external gauge field). One could worry whether the subdiffusive mode in Eq.(2.41) is a spurious mode whose dispersion is not robust under higher-order corrections. We will discuss this problem in detail in Section 2.5 and we will concretely show in appendix A that this is not the case.
The remaining four modes are all non-hydrodynamic and their dispersions can be found by solving the following equation: where Ω p is the plasma frequency . . (2.44) One interesting feature of the gaps at finite density (2.44) can be observed in the limit of B = 0. More precisely, depending on the value of the charge density ρ (entering through the plasma frequency Ω ρ ), the dispersions in Eq.(2.44) can be purely imaginary or complex, i.e., (2.45) Finally, setting all the dissipative coefficients (e.g., σ = 0) to zero, one finds (2.46) The plasma frequency Ω p gaps out both the sound waves (the first equation in (2.46)) and the electromagnetic waves (the second equation in (2.46)). 7 As the sound waves become the magnetosonic waves (2.37) at finite B, one may say that finite density gaps out the magnetosonic waves as well [50]. For later use, we summarize the small B correction to (2.45) as follows. For small density, (σ 2 / 2 e 4 Ω 2 p ), we have where the magnetic field B produces a real gap. For large density (σ 2 / 2 e 4 Ω 2 p ), we have where the magnetic field B produces both a real and an imaginary gap.

Dynamical vs. external gauge fields
We finish this section by comparing the dispersion relations of the low-energy modes for the distinct cases of dynamical and external gauge fields. The results are summarized in tables 1-2. Let us briefly remind the Reader of the results in the case of external gauge fields. In the neutral case, the spectrum displays a energy diffusion mode and a subdiffusive mode with dispersion: together with the so-called cyclotron mode At the finite density, the dispersion of these modes are modified into where the real gap, Ω B := Bρ/( + p), is the cyclotron frequency. By looking at the summary tables 1-2, one can notice that the dispersion relations in the two cases (dynamical vs. external gauge fields) are noticeably different. However, as described below Eq. (2.41), there is one exception: the diffusion mode at finite density (cfr. Eq.(2.51) vs. Eq.(2.41)). At least at such order in the wave-vector k, the dispersions are identical.

A note on the regime of validity of first-order magnetohydrodynamics
Before concluding this section about the magnetohydrodynamic framework, we would like to clarify its regime of validity and the role of possible higher order corrections. Everything discussed so far is valid in the approximation of first-order linearised hydrodynamics. The latter is the statement that the constitutive relations for the stress tensor T µν and the U(1) current J µ are expanded in dissipative corrections up to terms which are linear, or firstorder, in the gradients. The first-order dissipative corrections, which have been indicated respectively as Π µν and ν µ , are only the first of an infinite series expansion in gradients. In particular, both the stress tensor and the U(1) current can be expanded as: where the terms indicated with suffix (n) refers to dissipative corrections beyond equilibrium which are n-th order in gradients. Higher order terms correspond to shorter timescales and lengthscales and they expand the validity of the effective description towards the microscopic world.
In our analysis, the first terms which are neglected are second order in gradients, i.e., T µν (2) , J µ (2) . Since all the dynamical equations contain at least an extra derivative, these neglected corrections enter into the dynamical matrix M(ω, k) as terms ∼ k 3 , where k is the wave-vector. All in all, this means that the dispersion relations within the second order magnetohydrodynamic framework should be extracted from: with C a matrix of k−independent coefficients. In the worst case scenario, all the entries of the matrix C are nonzero. This is obviously a very conservative view but at this point, without knowing the precise form of C, the safest. Following this argument, we can confidently trust the results from first-order magnetohydrodynamics only up to the order in which the corresponding hydrodynamic coefficients are not affected by the possible corrections appearing in C. This is a well known problem in hydrodynamics which sometimes leads to the appearance of spurious poles as well. See for example Section 2.6 in [82]. To make this point clearer, let us make an example. Let us assume that as a solution of det [M(ω, k)] = 0 we get a mode whose dispersion relation within the first-order approximation can be written as: with a i random complex numbers. Now, let us assume, that the same dispersion relation extracted from the second-order formalism reads: Then, we can trust the dispersion relation obtained from the first-order formalism only up to the order at which a i =ã i or, in other words, up to the order at which the higher order corrections do not play any role. In the rest of the manuscript, whenever we will refer to hydrodynamic predictions, we will always have in mind this first order formalism truncated up to the terms in the k expansion which can be trusted within this approximation.
Notice that all the dispersion relations written so far, apart from that of the subdiffusive mode, ω ∼ k 4 in Eq.(2.41) are at most order k 2 . It is then immediate to verify that such expressions would not be corrected by possible higher order terms ∼ k 3 . The case of the subdiffusive mode could be potentially different. In particular, both the k 3 coefficient, which is zero in the first-order approximation, and the k 4 one shown in the text might in principle be affected by higher-order corrections. Nevertheless, we have verified with an accurate numerical analysis that this is not the case. We refer the Reader to Appendix A for an extensive discussion on this point.
Finally, in Appendix A, we will investigate further the validity of the dispersion relations obtained from solving det [M(ω, k)] = 0 without worrying about possible higher-order corrections. In particular, we will show that the dispersion relations obtained in that way, by assuming somehow that no higher-order corrections appear, significantly extend the range of agreement between the numerical data and the hydrodynamic predictions. This is an a-posteriori proof that many of the higher-order corrections are either zero or negligible for the problem at hand. Of course, one is not guaranteed that this is generally the case.

Holography with dynamical boundary gauge fields
In this section, we study the dynamics of a simple holographic system at finite charge density and finite magnetic field in which the gauge field is taken as dynamical in the boundary field theory using mixed boundary conditions.

Holographic setup
Let us consider the Einstein-Maxwell action in (3+1) dimensions, where we set 16πG = 1 and the AdS radius L = 1. We use Latin indices M, N, . . . for the 4-dimensional bulk spacetime coordinates and use Greek indices µ, ν, . . . for the 3dimensional boundary coordinates. In addition, let us consider the background dyonic black-brane ansatz where r = ∞ is the location of the AdS boundary and B the magnetic field. The action in Eq.(3.1) allows for a simply dyonic black-brane analytic solution given by where µ is the chemical potential, r h the horizon radius, and the black-brane mass m 0 is determined by the condition f (r h ) = 0. The various thermodynamic parameters associated with such a solution can be derived as [59,63,67,83]. We identify the bulk on-shell action in Eq.(3.1) with the matter controbution S m [g µν , A µ ] in Eq.(2.5). Moreover, we add the following boundary terms The latter, together with the bulk part, Eq.(3.1), constitute the full boundary action which has to be compared with Eq.(2.5). As a consequence, the thermodynamic quantities for T µν and J µ , which include the contributions from the Maxwell kinetic term, are given by where (T, ρ, s, , p) are the temperature, charge density, entropy, energy and pressure density, respectively. One can easily verify that these expressions satisfy the Smarr relation + p = s T + µ ρ. Furthermore, using (3.5), one can compute all the thermodynamic susceptibilities in (2.28) holographically. Doing so, we have verified that some of them vanish for this concrete solution at zero charge density as anticipated in Eq. (2.30). Notice that while the trace of the matter contribution to the stress tensor vanishes, the trace of the total stress tensor does not. In particular, we have: This result is not surprising since it corresponds to the statement that Maxwell theory in 3 + 1 dimensions is scale invariant but not conformal invariant [84,85]. As a consequence, the trace of its stress tensor does not vanish but it is equal to the total divergence of a virial current [84]. Moreover, in presence of the magnetic field, the mechanical pressure is not equal to the thermodynamic pressure. In particular, we have: The other transport coefficients (σ , η) can also be computed and read where the conductivity σ is given in [59,63,67,83] 8 and the shear viscosity η in (3.8) is obtained from the fact that the KSS bound [87,88] is not violated in the presence of both charge density and magnetic field in (3+1) dimensions. 9 Finally, we observed that the term + p which frequently appears in the hydrodynamic expressions is not affected by the Maxwell kinetic term. In particular, as evident from Eq.(3.5), we have that: (3.9)

Boundary conditions for dynamical gauge fields
In order to investigate the quasi-normal modes, we consider the fluctuations δg M N and δA M as where g M N and A M are the background bulk fields given in Eq. (3.2). For convenience, we choose the radial gauge: δg tr = δg rr = δg xr = δg yr = δA r = 0. Additionally, we write the fluctuations in Fourier form using the following notations with the wave-vector k aligned along the x direction. We then construct four gauge-invariant variables (see for example [62]) as 8 See also [64,65,86] for the development of transport property where B is no longer taken to be of order one in derivatives. 9 On the contrary, in the higher dimensional case [89][90][91][92][93], the KSS bound is violated at finite B.
where the index of the metric fluctuation h M N is raised with the background metric (3.2). The number of gauge-invariant variables is related to the structure of the equations of motion at the linearized level. In our case, one can find nine second-order equations with five first-order constraints. This implies that there are four independent fluctuations and therefore four gauge-invariant variables. Then, we can study the quasi-normal modes by employing the determinant method [94] in which the source matrix, S, is constructed with the AdS boundary (r → ∞) expansion of the variables (3.12). Note that, at the AdS boundary, the gauge-invariant variables (3.12) are expanded as where the superscripts (L, S) denote the leading/subleading term in the asymptotic expansion.
At this point, it is fundamental to understand how to construct the matrix S (3.18) appearing in the determinant method in the case of dynamical gauge field. For this purpose, we need to consider the boundary action (3.4). Then the variation of the total action S on-shell + S boundary produces the following equation at the AdS boundary, where S on-shell is the on-shell action from (3.1) and Π µ the radially conserved bulk current obtained from the Maxwell equation: ∂ r ( √ −g F rµ ) = ∂ r Π µ = 0. Notice that λ parametrizes the ratio between the bulk electromagnetic coupling and the boundary one. Since we have fixed the bulk one to unity, then λ corresponds directly to the boundary coupling as in the hydrodynamic description of the previous sections.
The first order variation of each terms is given by where, for convenience, gauge-invariant variables (3.12) are used. Using the Maxwell equa-tions, we then get: (3. 16) In the following, we will set h xy , h (L) yy to zero. As shown in [95], those terms would contribute only to finite local counterterms in the on-shell action and therefore they would not modify the structure of the poles that we are interested in. Doing that, Eq.(3.16) becomes Before continuing, one remark is in order. In principle, the electric and magnetic susceptibilities in Eq.(2.21) could be computed directly knowing the expression for the matter pressure p m . In the case of holography, we are able to easily compute χ BB since we dispose of a background solution at finite magnetic field. Nevertheless, we do not know how to compute the electric susceptibility χ EE because introducing a background electric field will inevitably make the full solution time dependent. Therefore, in order to proceed, we will assume that χ EE = 0. As we will see, this assumption will turn out to be a good approximation in the limit of small EM coupling, λ/T 1, but not in general (see Fig. 16 below).
We cannot exclude that this might be one of the reasons behind the disagreement between the hydrodynamic predictions and the holographic results in the concomitant limit of large B and large λ. Given this clarification, within this assumption, the electric permittivity ( e ) and magnetic permeability (µ m ) satisfy where 1/µ m = −2∂p/∂B 2 can be computed via (3.5). Interestingly, for the simple Reissner-Nordstrom solution considered, we find: This is tantamount to saying that the dual field theory is the avatar of a diamagnetic material as already noted in [14,96,97].

Results at finite electromagnetic coupling
By following the method just outlined, we are now ready to compute the low-energy excitations in our holographic model. The phase space of our system is defined by three scale invariant parameters (µ/T , B/T 2 , λ/T ). For the moment, we mainly focus on the case of small EM interactions, λ/T = 0.1. We will discuss in detail the effects of dialing the EM coupling λ in Section 4.5. We study the quasi-normal modes at zero density (µ/T = 0) and finite density (µ/T = 0), separately. Moreover, to avoid clutter in the presentation of the results, for the pure imaginary dispersion relations (e.g., the shear diffusion mode in Eq.(2.31)), we only display the imaginary part in all the figures.
Unless indicated otherwise, in all the figures of this manuscript solid lines will refer to the hydrodynamic predictions as explained in Section 2.5. On the contrary, colored dots will represent the numerical results obtained from the quasi-normal modes (QNMs) analysis.

Zero density
Let us start from the case of zero density ρ = µ = 0. We display the dispersion relation of the quasi-normal modes at zero magnetic field B/T 2 = 0 and small EM coupling λ/T = 0.1 in Fig. 1. We find that the numerical results are well matched with the dispersion relations from hydrodynamics in the expected range, k/T 1. In particular, Fig. 1 presents: • the sound waves with dispersion as in Eq.(2.31) (red); • the shear diffusion mode with dispersion as in Eq.(2.31) (yellow; see the inset); • the EM waves as in Eq.(2.32) (or (2.34)) (blue); • the damped charge diffusion mode as in Eq.(2.35) (green). 0.
The numerical results are still well fitted by the hydrodynamic formulae. The validity of the hydrodynamic framework and the match to the numerical data upon dialing the value of the magnetic field will be discussed in more detail in Section 4.4.

Finite density
From the hydrodynamic analysis of Section 2, at finite density, we do expect two gapless modes, Eq.(2.41), and four gapped modes, Eq.(2.44). In particular, depending on how large the density is, the gapped modes can exhibit distinct behaviors given by: 4 Ω 2 p ). In this section, as representative examples for each case, we choose µ/T = 0.5 for the small density case and µ/T = 5 for the large density case. A more detailed discussion about the role of the chemical potential and the transition between the two regimes will be presented in Section 4.3.
In Fig. 3, we display the quasi-normal modes at µ/T = 0.5 both for zero (top panel) and finite, but small, magnetic field (bottom panel). In both cases, the red data correspond to the diffusive mode in Eq.(2.41), the yellow data to the subdiffusive mode Eq.(2.41) and the green/blue data to the gapped modes Eq.(2.47). The strongest effect of the finite magnetic field appears in the gapped mode, Eq.(2.47). There, B produces a real gap which is absent for B = 0 (see the difference between top and bottom panels). In the case of small charge (top panel of Fig.3), first order hydrodynamics is able to generally predict the existence of the k-gap but not its location accurately if that is large (in units of k/T ). To be more precise, the k-gap curve appears well fitted by hydrodynamics only when the momentum gap is small, i.e. in the so-called quasihydrodynamic regime [51] (see Fig. 13 later in the discussion).   In Fig. 4, we show the quasi-normal modes at large chemical potential for both zero and finite magnetic field (respectively top and bottom panels therein). In all figures, the red data correspond to the diffusive mode in (2.41), the yellow data to the subdiffusive mode (2.41) and finally the green/blue data to the gapped modes (2.48). Note that the dispersion of the gapped modes is now given by Eq.(2.48).
At B = 0, the mode in Eq.(2.48) exhibits a real gap which was absent in the case of zero charge density (Eq.(2.47)). The value of the gap corresponds to the plasma frequency Re(ω) = ±Ω p . Additionally, at finite charge density, the magnetic field B contributes to both the real and imaginary parts in the limit of zero wave-vector.
In summary, using mixed boundary conditions as explained in the previous sections, all the quasi-normal modes and their dispersion relations are consistent with the expectations from magnetohydrodynamics at low wave-vector. This is the concrete proof that the modified boundary conditions are indeed rendering the gauge field at the boundary dynamical and that the boundary physics is accurately described by relativistic magnetohydrodynamics.

The effects of a finite chemical potential
In the previous sections, we did not explore in detail the role of the chemical potential on the dispersion relations of the low-energy modes. Here, we present additional results about the µ dependence at fixed magnetic field. For this purpose, we fix λ/T = 0.1 and B/T 2 = 0.
In order to proceed with this analysis, we first distinguish the two regimes of small     and large chemical potential as For λ/T = 0.1 and B/T 2 = 0, the same inequalities can be directly expressed as follows where µ/T ∼ 0.56 corresponds to σ 2 / 2 e ∼ 4 Ω 2 p . Let us start with the first regime of small chemical potential. In Fig. 5, we display the dispersion relations of the low-energy excitations at finite wave-vector.
At zero chemical potential (left panel of Fig. 5), we observe a gapless sound mode together with the diffusive mode which substitutes the propagating EM wave. The left panel of Fig. 5 can also be found in Fig. 2. By making the chemical potential finite, also the sound mode acquires a wave-vector gap and stops propagating at small k. The nonhydrodynamic modes (see for example the central panel in Fig. 5) are located, in the limit of small density, at: as derived in Eq.(2.45). As we increase µ further, the two values in Eq.(4.3) gets closer and the different non-hydrodynamic modes approach each other on the negative imaginary frequency axes (see right panel in Fig. 5). Exactly at the critical value, (µ/T ) * ∼ 0.56 , all the imaginary parts of the non-hydrodynamic modes are equal and given by iω = σ/(2 e ). Moreover, increasing the charge density the k-gap of EM waves becomes smaller, while its imaginary part at large wave-vector remains constant. At the critical value, the k-gap becomes exactly zero and after that a real gap appears. The dispersion relations of the modes above the critical value (µ/T ) * ∼ 0.56 are shown in Fig. 6. Indeed, as we increase µ further, beyond the critical value µ/T ∼ 0.56, all the lowest   non-hydrodynamic modes discussed before acquire a finite real gap and show the same dispersion at zero wavevector, k = 0, as already discussed in Eq. (2.45). This behavior is confirmed by the numerical data in Fig. 6.

The strong magnetic field limit
Now, we extend our analysis to the regime of strong magnetic field, B/T 2 1. Our main scope is to better understand the regime of validity of the magnetohydrodynamic description beyond the small B limit. Here, we present the results at fixed λ/T = 0.1.
For this purpose, we first investigate the thermodynamic parameters ( , p, s) and the transport coefficients (σ, η, diffusion constant, etc) appearing in the magnetohydrodynamic description. We recall that p is the thermodynamic pressure and does not correspond to the mechanical pressure.

Thermodynamics and transport coefficients
In Fig. 7, we display the thermodynamic parameters ( , p, s) together with the electric conductivity σ at fixed density as a function of the magnetic field. The shear viscosity η is trivially related to the entropy s via the KSS relation (3.8) and it is therefore not shown. In the weak B field regime (B/T 2 1), all the observables are µ-dependent constants. On the other hand, in the strong B field regime (B/T 2 1), we find the following asymptotic behaviors: where all the coefficients are independent of the ratio µ/T . Interestingly, increasing the value of B, the thermodynamic pressure p becomes negative  Figure 7. Thermodynamic parameters ( , p, s) and conductivity σ from weak to strong magnetic field at µ/T = 0, 5 (black, red). at a critical value B * . This critical value can be obtained analytically from Eq.(3.5). Its behavior as a function of µ/T is showed in Fig. 8. Using the Smarr relation (p = − + sT + µρ), one can also understand the behavior of both p and σ in the strong B regime. In that limit, ∼ B 3/2 scales faster than the entropy, s ∼ B. As a result, at large B one has p ∼ − and therefore a negatively divergent thermodynamic pressure. In addition, the electric conductivity can also be rewritten as σ = sT sT +µρ 2 which approaches unity in the strong B limit. Similar results have been obtained in AdS 5 using a higher-form language [47].  Table 3. Approximate asymptotic behavior for the transport coefficients appearing in the dispersion relations of the hydrodynamic modes in weak-and strong-field limits at zero density.

Transport coefficients and magnetohydrodynamics
We now turn to the analysis of the dispersion relations of the low-energy modes (lowest QNMs) and in particular of the coefficients appearing up to order k 2 . Our task is to verify the validity of the hydrodynamic description at large magnetic field. Let us start with the simplest case of a neutral plasma. Here we do expect four gapless modes: two magnetosonic waves together with the two (shear/magnetic) diffusive modes The specific expression for the coefficients above is lengthy and is provided in the GitHub repository available here. In Fig. 9, we display (v ms , Γ ms ) as a function of B. We find that the speed of magnetosonic waves interpolates between the conformal sound speed v 2 ms = 1/2 at weak B field and the speed of light v 2 ms = 1 at strong B field. The normalized attenuation constant Γ ms T displays a non-monotonic behavior. It vanishes at large B and it asymptotes a constant at zero B. Notice that, using the asymptotic forms presented in table 3, a negative value of λ would result in a superluminal magnetosonic waves at strong magnetic field. This is not surprising since λ must be taken as positive.
We also show the B-dependence of the diffusion constants for shear and magnetic diffusion in Fig. 10. We find that all the diffusion constants vanish in the strong B regime, which is consistent with the previous literature [47,62,98]. We find that all the transport coefficients of the gapless modes, obtained by fitting the numerical dispersion relations, are in perfect agreement with the magnetohydrodynamic formulae even in the large B limit. This is somehow surprising but, as we will see, strongly dependent on the value of the EM coupling λ. We will comment more on this point in the next sections and in the conclusions. This agreement implies that the formula for the conductivity σ given in Eq. (3.8) works well for all values of B, even beyond the small B regime. To complete this section, we can also analytically derive the asymptotic behavior of all these coefficients. For the zero density case, these are given in table 3.
In Fig. 10 we also show the value of the damping of the first non-hydrodynamic mode in function of the dimensionless magnetic field. In this case, the prediction from magnetohydrodynamics are not in well agreement with the numerical data for B/T 2 1. This is not surprising. It is simply due to the fact that the imaginary part of this non-hydrodynamic mode becomes large, i.e., Im(ω) ∼ T , and the mode moves away from the regime of validity of linearised hydrodynamics.
Next, we perform the similar analysis at finite density. In this case, as demonstrated in Section 2.3, we have two gapless modes, the longitudinal diffusive mode and subdiffusive mode, whose dispersions are given by The concrete form of the diffusive parameters is cumbersome but can be found in the GitHub repository available here. For the caveats related to the validity of the subdiffusive dispersion within first order hydrodynamics see the discussions in Section 2.5 and appendix A. In Fig. 11, we display the B-dependence of the two diffusive parameters. D long vanishes in the strong B limit independently of the value of the charge density. On the contrary, D subdiff reaches a constant value at B/T 2 → ∞. More precisely we find that in the strong B limit: In addition, we also find that both transport coefficients are suppressed at larger density. Once again, the results obtained by the fitting method are in good agreement with mag- netohydrodynamic predictions at finite density even in the strong B regime. The apparent fluctuations in the numerical data visible in Fig. 11 are just due to numerical precision. Finally, we discuss also the dynamics of the non-hydrodynamic modes at finite charge density by dialing the strength of the magnetic field B.
In Fig. 12, we display the B-dependence of the imaginary and real gaps at zero wavevector, k = 0. In the weak B regime, as demonstrated in Section 2.3, the behavior of the gaps depend on the density. The precise expressions are provided in Eq.(2.47) for small density and in Eq.(2.48) for large density.
Interestingly, in the strong B limit, we find that one of the two pair of modes approaches the origin of the complex frequency plane. More specifically, both its real and imaginary parts at zero wave-vector go to zero in the limit B/T 2 → ∞. This mode in the strong B limit becomes an emergent propagating magnetosonic wave with speed v 2 ms = 1 and vanishing attenuation constant. It is very tempting to describe this mode as an emergent photon in the strong B limit. We are not aware of similar observations in the previous literature. This point needs further investigation in the future. On the other hand, the remaining pair is pushed away from the hydrodynamic limit since both its real and imaginary part diverge. Notice that whenever its imaginary part becomes too large, magnetohydrodynamics breaks down and its predictions are not anymore in good agreement with the numerical data.
The asymptotic behavior for these non-hydrodynamic modes in strong B regime can be obtained analytically and it reads where (4.9) corresponds to the pair approaching the hydrodynamic limit in Fig. 12, while (4.10) to the pair which is gapped away in the large B limit.
In summary, we find that as long as the EM coupling λ is small, the magnetohydrodynamic predictions for the hydrodynamic modes are in perfect agreement with the numerical results even in the strong B regime. At this point, we are not able to provide a solid derivation of why this is the case and how general this is. We will comment further on this point in the conclusion section. In the next section, we discuss the effect of λ.

The role of the electromagnetic coupling
In the previous sections, we have fixed λ/T = 0.1 and considered only the limit of small electromagnetic coupling. We now investigate the role of the EM coupling λ and the interpolation between the two limits λ → 0 and λ → ∞. Similar analyses in the context of plasmons can be found in [38,40,41].
For this purpose, we first consider the dynamics of EM waves which is given within the magnetohydrodynamic framework by Eq.(2.32). In Fig. 13, we compare the numerical QNMs data from the holographic model with the magnetohydrodynamic preditions. As expected, EM waves are screened by Coulomb interactions and they start propagating only above a certain cutoff wave-vector k , the k-gap. The magnetohydrodynamic framework, in the first-order approximation, gives: From Fig. 13, we can indeed observe that at low wave-vector the EM waves are not propagating but they rather split into a diffusive mode and a non-hydrodynamic one, as predicted by magnetohydrodynamics. Moreover, for small values of the EM coupling, λ/T 1, the magnetohydrodynamic formula, Eq.(2.32), is in very good agreement with the numerical data and accurately predicts the value of the cutoff wave-vector k . The onset of propagation moves to larger wave-vectors by increasing the EM coupling λ or equivalently the magnetic field B (see for example the red data in Fig. 13) and the accuracy of the magnetohydrodynamic description decreases. This is simply the sign that higher order corrections in Eq.(2.32) become important at k/T ∼ O(1). Interestingly, the imaginary gap of the nonhydrodynamic mode can be derived analytically (see Eq.(4.14) below) and it is in perfect agreement with the numerical data even when Im[ω]/T 1. In order to derive the inverse relaxation time of the non-hydrodynamic mode analytically, we notice that the equations of motion for the gauge field fluctuations (δa i=x,y ) decouple in the limit of k = 0 and µ = B = 0. Then, the equations can be solved analytically and give the following leading/subleading coefficients near the AdS boundary (4.12) Continuing, we use the boundary conditions described in Section 3.2 which at k = 0 are given by By combining the two last equations, we finally arrive at two independent solutions (4.14) The first mode corresponds to the hydrodynamic diffusive mode at k = 0 while the second mode is the non-hydrodynamic one visible in Fig. 13 (top panel). Taking into account that in the limit µ = B = 0 we have σ = 1, this analytically confirms the identification of the parameter λ in the b.c.s. as λ = 1/ e , at least in the regime of small λ coupling. 10 These results are in agreement with those found in [41]. It is interesting to notice that the inverse relaxation time of the non-hydrodynamic mode can be analytically derived for arbitrary values of the EM coupling λ, in perfect agreement with the numerical data. Following the trend in Fig. 13, one could anticipate the appearance of a propagating free photon for λ → 0. This is indeed the case as explicitly shown in Fig. 14 for λ = 0. This outcome is maybe not surprising since the mixed b.c.s. imposed reduce in the limit of λ → 0 to the decoupled boundary Maxwell equations in vacuum: which clearly displays a freely propagating photon. In other words, we can understand this limit as the one in which the EM interactions are vanishing and therefore all the effects of polarization and screening disappear. The dynamics of the Maxwell field decouples from the current. From a technical perspective, this comes from the fact that the b.c.s. used do not reduce to the standard Dirichlet ones in the limit of λ → 0. On the contrary, they reduce to the Dirichlet b.c.s. times an independent and external factor (ω 2 − k 2 ) = 0. The spectrum of the theory in this limit is then the same as the one for a CFT with nondynamical U(1) symmetry times a freely propagating (and infinitely living) photon. Our results are consistent with those found using higher-form symmetries in [54]. What about the opposite limit of λ → ∞? As already hinted in Section 3.2, our boundary conditions in this limit do not boil down to the ones usually defined as alternative quantization. In particular, we do not fix directly the subleading term of the bulk gauge field, as for example done in [23]. Before getting there, let us first ask a different question. Is the hydrodynamic framework of Section 2 still reliable in the large λ limit? From Fig.  13, one can already notice that, in the small wave-vector regime, the gapless modes can be still well described by the hydrodynamic predictions even at large λ. In Fig. 15, we show the speed and attenuation constant of the magnetosonic waves together with the diffusive parameters of shear and magnetic diffusion for different values of B and λ. We do observe that the numerical data are not matching well the predictions from magnetohydrodynamics in the regime of large magnetic field and concomitant large EM coupling. Interestingly, the dynamics of shear diffusion is still perfectly described by hydrodynamics.
In addition to the hydrodynamic modes, one can further discuss the λ dependence of the gap of the non-hydrodynamic mode. For this purpose, we focus on the behavior of T Figure 14. The emergent propagating photon at zero density, zero magnetic field and zero EM coupling λ/T = 0.    Eq.(2.44). Fitting their dispersion relation, we obtain the value of the plasma frequency Ω 2 p and the damping parameter σ/ e numerically. In Fig. 16, we find that at small EM coupling their values are in good agreement with our expectations from magnetohydrodynamics: Away from the small λ limit, both the plasma frequency and the inverse relaxation time approaches a constant which is not anymore well approximated by the hydrodynamic predictions. Two comments regarding this discrepancy are in order. First, this might imply that, for large values of λ, the non-hydrodynamic mode are already too far away from the hydrodynamic window and, not surprising, the predictions from hydrodynamics, Eq.(2.44), are not reliable anymore. Second, note that for the solid lines in Fig. 16, we set χ EE to zero and assume that χ EE does not depend on λ. This is probably not the case. In particular, we do not expect χ EE to be generically zero in our holographic model. It would be interesting to find an independent way to calculate χ EE holographically. The main difficulty is that, by switching on a background electric field, time dependence is unavoidable.

Alternative quantization and a bulk experiment with non-canonical kinetic term
In this section, we discuss the possibility of modifying the nature of the dual field theory not by using boundary conditions nor by a Hodge duality in the bulk but rather by substituting the original Maxwell term F 2 with a higher derivative action of the form (1.6). Moreover, we compare the results of this experiment with the results obtained in the λ → ∞ limit.

Higher derivative bulk action
Let us consider a higher derivative bulk action as where the AdS asymptotic behavior of the gauge field reads Note that the higher derivative bulk action in Eq.   the gauge field (5.2) can be leading or subleading, i.e., This implies that, for the bulk action (5.1) with N > 3, the coefficient A is not anymore a source for the external field A µ but rather the expectation value of the conjugated current J µ . In other words, the standard quantization for the bulk theories with N > 3 appears equivalent to the alternative scheme for those with N < 3 (see Fig. 17 for a graphic summary). Fixing the value of A (0) µ (t, x) in a theory with N > 3 does not correspond to setting the value of a non-dynamical external gauge field (the source in common jargon). The scope of this section is to try to make sense of a bulk theory with a non-canonical kinetic term with N > 3 and in particular to understand which is the nature of its dual field theory. Finally, we would like to ask whether this bulk theory is equivalent, and in which sense, to using the standard Maxwell kinetic term (N = 2) but with alternative boundary conditions. For brevity, we will often refer to the N = 4 theory with Dirichlet b.c.s. as the "F 4 theory" and to the Maxwell action N = 2 with alternative b.c.s. as the "F 2 theory".

Low-energy modes and magnetohydrodynamics
As a concrete example, we will focus on the neutral state and compute the dispersion relation of the lowest quasi-normal modes in the F 2 model (N = 2) with alternative quantization (Neumann b.c.s) and in a high-derivative F 4 model (N = 4) with standard quantization (Dirichlet b.c.s). Note that in both cases the boundary condition corresponds to fix the value of A (1) µ (t, x) at the boundary. Notice also that, for theories with N = 2, one cannot safely take the limit of zero charge and zero magnetic field since the bulk fluctuations would then suffer of a strong coupling problem. One can explicitly check this fact by looking at the generalized bulk Maxwell equation: (5.4) in which the effective EM coupling in the bulk would be given by: such that g eff −→ ∞ when F 2 → 0. This is totally analogous to the case in which the bulk action is a higher-derivative theory for massless scalar, see discussion in [74].
Let us now consider the neutral state with a finite magnetic field for which the background bulk solution is given by with the corresponding thermodynamic variables For convenience, we have definedP =: T xx . Note that,P is not equal to the thermodynamic pressure in presence of a magnetic field [59][60][61][62][63]. 11 We then consider the fluctuations defined in Eq. (3.11) to study numerically the quasi-normal modes of the system. We impose the Neumann/Dirichlet b.c.s for the gauge fields, while we keep the Dirichlet b.c.s for the metric fluctuations. In what follows, Neumann/Dirichlet b.c.s denote the boundary conditions for the gauge fields only. We find that the quasi-normal modes for both the F 2 model (with Neumann b.c.s) and F 4 model (with Dirichlet b.c.s) exhibit four gapless modes: a pair of sound waves, a shear diffusion mode and a magnetic diffusion mode. Moreover, we empirically observe that their dispersion relations at finite magnetic field are well approximated by the following formulae: The meaning of the diffusion constant D is associated with the magnetic diffusion constant D mag at small B explained below. The numerical results for the dispersion of sound waves and shear diffusion (first two set of modes in Eq.(5.8)) are shown in Fig. 18 and they are in perfect agreement with the formulae above. Let us emphasize that in the limit of small magnetic field, the formulae presented in Eq.(5.8) can be consistently derived from   hydrodynamics by taking the λ → ∞ limit. In particular, in that regime, we find from hydrodynamics v 2 ms = which agree with Eq.(5.8) at small B. Similarly, from hydrodynamics, the diffusion constants of shear and magnetic diffusion, in the limit of λ → ∞ and small magnetic field, are  given by Notice that, using Eq.(2.24) in the λ = ∞ limit, the magnetic susceptibility is given by: and it is negative. Then, D mag > 0. For the magnetic diffusion mode, we find agreement between the hydrodynamic predictions and the numerical data only in the low-B limit and for the F 2 model (see left panel in Fig. 20). We do not believe that the failure of the magnetohydrodynamic theory for the F 4 theory in the small B limit is meaningful. On the contrary, that is just a signal of our failure in correctly identifying the hydrodynamic transport coefficients, such as the conductivity σ, in the F 4 theory. We plan to revisit these results and the transport dynamics of the F 4 theory in more detail in the future. In summary, this analysis once more shows that the magnetohydrodynamic theory only fails in the EM sector and only in the concomitant limit of large (in this case infinite) EM coupling and large magnetic field. Moreover, it shows that modifying the bulk action with a higher-derivative kinetic term is equivalent to considering the standard kinetic term with alternative boundary conditions. This is in close analogy with the case of holographic models with broken translations [72,76].

Further comments on magnetic diffusion
As already mentioned, we have not been able to match the magnetic diffusion constant for the F 4 model using our magnetohydrodynamic theory because we could not robustly derive the transport coefficients needed. In particular, in order to achieve this, one would need to understand how to extract the electric conductivity σ and the magnetic susceptibility χ BB in the higher-derivative F 4 model. Nevertheless, one can gain further insights on the diffusion constant D by performing a perturbative bulk analysis. From the equation for the fluctuations, one can check that such a diffusive mode originates from the gauge fluctuation sector δa y which couples to the metric fluctuation sector at finite B. In the limit of a vanishing B, one can find that the gauge sector decouples so that one can study the dynamics of the gauge fields on a fixed Schwarzschild background.
In this decoupling limit the equation of motion for δa y reads where f (r) is given by Eq.(5.6) in the limit of B = 0. Implementing standard perturbative techniques, we are able to solve the above equation analytically and obtain the Green's function for the operator dual to δa y . By looking at the poles structure of the latter, we can identify the presence of a mode whose dispersion is given by As shown in Fig. 20, the analytic result above is consistent with the numerical results in the limit of B/T 2 → 0.

On the existence of a free boundary photon in the alternative quantization scheme
For the case of alternative boundary conditions, a propagating photon ω = ±k was identified in the neutral AdS 3 case [23]. More precisely, the emergent photon was found by imposing the vanishing of the subleading term of the gauge field in the gauge-invariant way as Z (S) In order to understand the photon dispersion from (5.14), it is useful to re-express (5.14) as Using the AdS boundary expansion, one can also find the following relation which is nothing else that the conservation of the current. Putting Eqs.(5.15)-(5.16) together, one immediately obtains which has a trivial solution at ω = ±k. This is exactly the propagating photon observed in [23]. However, if one considers as boundary conditions the vanishing of the external current, as we do, the situation is different. Instead of Eq.(5.14), one has to impose The frequency dependent pre-factor cancels out and the emergent photon does not appear anymore.
In order to justify our findings, let us have another look at the standard Maxwell equation for electromagnetic waves in matter, given by: Solving this equation together with (3.19) in the λ → ∞ limit gives two modes with dispersion: where the diffusive mode is the magnetic diffusion given in (5.10) at B = 0. This means that in such a limit the photon disappears. The only way that a photon could emerge in the limit of λ → ∞ would be if σ/ e → 0. As proved numerically in Fig. 16, this is certainly not the case. In summary, in the limit of infinite EM coupling, λ → ∞, we do not find any propagating photon.

Conclusions
In this work, we have studied the low-energy dynamics of bottom-up holographic models at finite (free) charge density and magnetic field in presence of dynamical electromagnetism at the boundary. We have achieved the presence of a local U(1) symmetry in the boundary field theory by appropriately modifying the boundary conditions for the bulk gauge fields. We have then compared the numerical results from the holographic models with the predictions of magnetohydrodynamic theory in 2 + 1 dimensions. We have found perfect agreement between the two results. This proves that modified mixed boundary conditions for the bulk gauge fields provide the correct magnetohydrodynamic phenomenology in the dual field theory. Importantly, our work proves that the dual higher-form bulk description (e.g., [47]) is not necessary to obtain dynamical electromagnetism in the boundary field theory of bottom-up holographic models. This is somehow not surprising given that one could derive a precise duality between higher-form models and standard Maxwell model using different mixed boundary conditions [56].
Interestingly, we numerically observe the breaking down of magnetohydrodynamics only in the concomitant limit of large EM coupling, λ/T 1, and large magnetic field, B/T 2 1. On the contrary, we find that, as far as the electromagnetic coupling is small, the predictions from magnetohydrodynamics at small frequencies and wave-vectors are in good agreement with the numerical data even in the limit of large magnetic field. Despite the magnetic field is treated in the "strong field" limit, where B ∼ O(1), this is somehow surprising. A few possible explanations arise. (I) This is a pure coincidence valid only for the model considered. (II) We have not been able to probe very large values for the magnetic field B where maybe the predictions from magnetohydrodynamics would fail. (III) We are witnessing another case in favor of "unreasonable effectiveness" of hydrodymamics. (IV) A solid argument behind this observation exists but we have not found it yet. We find this aspect particularly interesting and we leave further investigation of this open question for the near future.
More in general, our results provide a good playground to describe holographic models with finite electromagnetic interactions in view of possible applications to plasma physics, astrophysical objects and condensed matter systems. A set of additional open questions is left for future studies.
• What is the emergent physics at infinite electromagnetic coupling in 2 + 1 dimensions and how can that be described (see [23] for earlier discussions on this point)?
• Can we find a way to compute the electric susceptibility χ EE from holography and improve our understanding at large EM coupling? It would be interesting to understand whether the models and analyses of [99][100][101] could shed light on this point.
• Can we understand better the large B limit and in particular test the recent claims made in [102] about magnetic diffusion?
• Can the numbers of gapless (hydrodynamic) modes be understood in terms of symmetries? In this sense, is the plasma frequency related to the explicit breaking of any symmetry? If the photon can be identified as a Goldstone mode, is the plasmon the manifestation of a pseudo-Goldstone mode?
• Is there an emergent photon in the strong B regime? And, why?
• What is the correct dual field theory interpretation of the higher-derivative F 2N bulk model? Which are the corresponding transport properties?
• Are the modified boundary conditions giving the correct phenomenology of superconductors once the U(1) symmetry is spontaneously broken [103] (see for example [30])?
Finally, it would be instructive to re-do our computations in a four dimensional boundary theory in which magnetohydrodynamics displays a richer, and angle dependent, spectrum with for example Alfvén waves and fast/slow magnetosonic waves [47]. Also, it would be interesting to extend our results in presence of a chiral anomaly, as done in [71] for the case of external gauge fields. We plan to report on some of these issues in the near future.

A The unreasonable effectiveness of first-order magnetohydrodynamics
In this Appendix, we provide a few more details about the discussion of Section 2.5. First, we analyze in more detail the dispersion relation of the subdiffusive mode, ω = −iD subdiff k 4 and the validity of the predictions from first-order magnetohydrodynamics presented in the main text. As already argued in Section 2.5, taking a overly pessimistic attitude, one could expect that the corrections from second-order hydrodynamics could modify the dispersion relation of the subdiffusive mode in Eq.(2.41) at order k 3 . Fortunately, this is not the case. In Fig. 21, we present an accurate analysis of the dispersion relation of the subdiffusive mode at low wave-vector. As evident from there, the dispersion displays a k 4 scaling up to low wave-vector, indicating the k 3 correction is not present. Additionally, the prediction from first-order hydrodynamics of the k 4 coefficient matches perfectly the data up to k/T ≈ 0.04. This indicates that the dispersion relation extracted from first order hydrodynamics is reliable and, at least up to k 4 order, no corrections appear.
Finally, we discuss the validity of the first order hydrodynamic formalism used in the main text. As explained in detail in Section 2.5, apart from the subdiffusive mode in Eq.(2.41), in the main text we take a conservative attitude and we consider the results from the first-order formalism only to the order at which we are sure they cannot be affected by second-order corrections, ∼ k 3 . Here, we want to relax this attitude and consider the dispersion relations from first-order hydrodynamics without expanding the solutions of det(M(ω, k)) = 0 at small wave-vector. Since M is a 6 × 6 square matrix and every entry is at most order k 2 , we do expect the final polynomial to be order six in frequency and order twelve in wave-vector. To perform this analysis we consider only the QNMs shown in the bottom panel of Fig. 3 in the main text. We show the results in Fig. 22: we also provide the data for all other cases (corresponding to Fig. 1-Fig. 6 in the main text) in the GitHub repository available here. First, in the left panel, we show that this unjustified relaxed attitude significantly enlarges the validity of the hydrodynamic predictions. The latter are now in perfect agreement with the numerical data up to k/T ∼ 0.05. This has to be contrasted with the results shown in the bottom panel of Fig. 3 in which the first-order hydrodynamic predictions fail already around k/T ∼ 0.013. In order to make this more evident, in the right panel of Fig. 22, we show both predictions for a single mode. In dashed red line we display the conservative predictions used in the main text while in solid red line the enlarged attitude described in this Appendix. The difference is evident. What is this suggesting us? These results are telling us that, at least for the system at hand, the higher order corrections which come from expanding the constitutive relations at higher-order are subleading for a quite large range of wave-vector. This might certainly not be the case in general, but it is nevertheless a nice and interesting observation. Let us conclude this appendix saying that, even considering the ∼ k 12 solution from det(M(ω, k)) = 0, the hydrodynamic predictions will not match the data at arbitrarily large values of k. Increasing the range of k further, one would see deviations as well.