Electromagnetic response of strongly coupled plasmas

We present a thorough analysis of the electromagnetic response of strongly coupled neutral plasmas described by the gauge/gravity correspondence. The coupling of the external electromagnetic field with the tower of quasi-normal modes of the plasmas supports the presence of various electromagnetic modes with different properties. Among them we underline the existence of negative refraction with low dissipation for a transverse non-hydrodynamical mode. Previous hydrodynamical approaches have shown the ubiquitous character of negative refraction in charged plasmas and the absence thereof in neutral plasmas. Our results here extend the analysis for neutral plasmas beyond the hydrodynamical regime. As an application of these new insights we briefly discuss the case of the quark gluon plasma in the temperature dominated regime.


Introduction
The study of the linear response of a system to a small external perturbation is an essential tool to gain insight on the collective behavior of its constituents and to provide information about the transport properties at thermal equilibrium 1 . Media having electrically movable components that can easily respond to an external electromagnetic field are ubiquitous in physics: high-energy physics plasmas and condensed matter superconductors are instances of this sort. For such media usually the electromagnetic field interacts with thermal and mechanical modes; hence the electromagnetic energy finds various channels of propagation and dissipation through the medium. The response of the system to an external electromagnetic field is then in general non-local and the medium is said to be spatially dispersive (see for instance [2]). This rich dynamics generates an intriguing phenomenology which can feature exotic electromagnetic effects such as negative refraction (NR) [3] and additional light-waves (ALW) [4] as we will shortly review.
Exotic electromagnetic phenomena have recently attracted an intense theoretical and experimental (and even technological) interest especially in relation to artificial systems called meta-materials [5]. Inspired by the amazing developments obtained on metamaterials, in [6][7][8] it was first discovered that charged fluids which admit a hydrodynamical description present exotic electromagnetic properties as a general feature. The negative refraction phenomenon (in a certain range of frequencies) and the presence of additional light-waves are shown respectively in [6,7] and [8]. The former implies that the electromagnetic energy flux and the phase velocity of an electromagnetic wave through the medium propagate in opposite directions [3,5]; the latter means that, even if isotropic, a medium supports multiple electromagnetic waves with the same frequency but different wave-vectors [4]. Following these results, various authors investigated these topics in a wide variety of different setups, see for example [9]. In particular [10] provides an analysis of the actual signature of the abovementioned phenomena in condensed matter systems and it proposes suitable experiments to measure such effects.
Even though much of the previously stated theoretical progress about NR and ALW was obtained within the framework of strongly coupled plasmas, it turned out that the qualitative behavior of the discovered phenomena has a more general valence and it constitutes a universal property of electrically charged hydrodynamical systems [7]. The main purpose of the present paper is to go beyond the regime described by hydrodynamics in systems characterized by strong coupling physics. In particular, we will show that the coupling of the external electromagnetic field with the tower of quasi-normal modes of strongly coupled plasmas leads to the presence of various electromagnetic modes which present different properties. For instance, it is important to stress that we found negative refraction with low dissipation for a transverse non-hydrodynamical mode.
Gaining insight on the strongly coupled dynamics is crucial in order to describe some very interesting phases of matter which have recently been realized both in high energy physics experiments (e.g. the quark gluon plasma (QGP), see for example [11]) and in condensed matter systems (e.g. high-T c superconductors, see for example [12] for a recent review 2 ). Moreover, the existence of strong interactions among the various resonances in some of the artificially engineered meta-materials, the so called stereo-meta-materials [14], provides yet another interesting phenomenological stimulus to investigate the electromagnetic behavior and response of strongly coupled systems.
The idea of performing an analysis beyond the hydrodynamical regime aims at studying the electromagnetic wave modes that are not controlled by the hydrodynamical uni-versal behavior accounted by the equations of fluid dynamics 3 . Particular attention is then paid to observe whether the exotic phenomena arising already in hydrodynamics are present also beyond the long wavelength and long time approximations and what are their properties in this regime.
It is well known that the dynamics at strong coupling is usually very difficult to study and, outside the regime constrained by hydrodynamics, it is quite difficult to provide universal predictions. However, some interesting theoretical progress has been recently obtained thanks to techniques coming from string theory. These new methods exploit a conjectured duality between a strongly coupled quantum field theory and a weakly coupled gravity theory; such duality is usually referred to as gauge/gravity correspondence [17]. The gauge/gravity duality (and similar correspondences) allows us to quantitatively study the correlation functions of specific models and, at the same time, to understand better many interesting and generic qualitative behaviors associated to various (universality) classes of strongly coupled systems. Insights from the gauge/gravity correspondence appeared to be relevant both to high energy physics and condensed matter (see for instance [18] to have two paradigmatic examples) especially in relation to the characterization of the phase structure and transport properties. Within the gauge/gravity correspondence framework, quantitative computations of expectation values and correlation functions of strongly coupled quantum systems at finite temperature are reduced to classical gravity calculations in specific black hole backgrounds (which are actually finite temperature gravitational systems). The thermodynamical observables of the strongly coupled system in equilibrium are mapped to classical properties of the dual gravity model. A semi-classical study of the equations of motion for the fields of the dual gravity model unveils both the equilibrium thermodynamics and the linear response (i.e. the slightly out of equilibrium physics). The latter being obtained considering linearized equations of motion (equipped with appropriate boundary conditions) describing small fluctuations around the gravitational black hole background.
To our aim, the benefit of the gauge/gravity framework is twofold: first of all it allows us to study correlation functions for any value of the frequency and wave-vector in a well defined and complete setup (essential to go beyond hydrodynamics); secondly, it provides interesting hints on generic responses of strongly coupled systems beyond the hydrodynamical regime.
A possible phenomenological application of our results is to QGP physics. The neutral plasma that we study could indeed mimic the QGP plasma in a regime dominated by the temperature where the finite chemical potential is neglected. In particular, we stress that the electrodynamical, non-hydrodynamical modes described in this paper have smaller wavelengths when compared to the hydrodynamical modes studied in [7]; hence they probe smaller distances which could be relevant for actual QGP samples produced in experiments.
The paper is organized as follows. In Section 2 we briefly review the linear response theory needed to study the electromagnetic properties of a medium. Section 3 contains our main results obtained by studying the first transverse electromagnetic waves through the strongly coupled plasma. We show in particular the generic presence of ALW and we demonstrate the existence of a propagating electromagnetic mode with negative refraction and very small dissipation. In Section 4 we describe the actual calculations venturing into technical details of both the gauge/gravity setup and the peculiarities of the necessary renormalization procedure. Section 5 is dedicated to the description of our numerical and semi-analytical methods and the checks that our results have undergone. In Section 6 we discuss very briefly some phenomenological applications of our study to the QGP. In Section 7 we show the results of an analogous analysis regarding the longitudinal lightwave modes. Eventually, in Section 8 we conclude stressing the significance of the results and possible future prospects. The appendices contain a characterization of the quasinormal modes of the plasma (both transverse and longitudinal) and the details about the computation regarding the longitudinal light-wave sector.

Setting the scene
The electromagnetic response of a system in local equilibrium to an external electromagnetic perturbation is described by the 2-point retarded correlation function of the electromagnetic current [19]. If, as in the case we are considering, the system has spatial dispersion, the response function depends on the space positions where the currents are evaluated. In the particular case of a homogeneous and isotropic system, the correlation functions depend only on the distances and they can be further decomposed into a transverse and a longitudinal part with respect to the wave-vector q of the external perturbation. In this case, as shown for example in [2,3], it is possible to describe the macroscopic electromagnetic properties of the system with three fields: D, E and B. Moreover, the linear response of the medium (in Fourier space) is accounted for by a single tensorial function depending both on the frequency ω and the wave-vector q, the so-called permittivity tensor ij (ω, q), As a consequence of the isotropic assumption, the permittivity tensor is expressed in terms of only two scalar functions describing the transverse and longitudinal response with respect to the spatial momentum q of the external perturbation. The 2-point retarded correlation functions, denoted as G T (ω, q) and G L (ω, q) for the transverse and the longitudinal part respectively, characterize completely the permittivity tensor: where e is the electromagnetic coupling constant. The knowledge of G T,L (ω, q) allows us to study the propagation and dissipation of the electromagnetic waves inside the medium at the linear level. Indeed, the Maxwell equations for the electromagnetic field provide the dispersion relation between q and ω both for the transverse and the longitudinal channel 4 . The Maxwell equations for the transverse and longitudinal part of the electric field E imply the dispersion equations: which admit solutions q(ω) corresponding to possible electromagnetic waves through the medium 5 . The information about the propagation and dissipation of an electromagnetic mode through the system is encoded in the refractive index which is defined as n 2 (ω) = q 2 (ω)/ω 2 considered on the dispersion relation of the mode under study. Albeit we study both the transverse and the longitudinal electromagnetic sectors, we focus mainly on the detailed analysis of the transverse sector (whose properties are more interesting to us). An analogue detailed treatment of the longitudinal sector is reported in Appendix A.

Results
In this Section we describe the main results of our analysis. We focus on the refractive indexes of the first electromagnetic transverse modes inside a strongly coupled plasma composed of electrically charged constituents but globally neutral. The plasma we study is the finite temperature state of a relativistic and conformally invariant particle physics theory composed of fermions and bosons 6 . The choice of this particular kind of plasma is due to the existence of very powerful theoretical tools (provided by the gauge/gravity correspondence) to study its strongly coupled dynamics. Indeed the present analysis can be read as a first example to investigate the potentiality of the gauge/gravity correspondence to explore exotic electromagnetic phenomena of strongly coupled, finite temperature systems beyond the hydrodynamic regime.
The plasma is assumed to be homogeneous and isotropic 7 . Previous studies have shown that plasmas described with the gauge/gravity correspondence have in general an infinite set of discrete quasi-normal (complex valued) self-frequencies [21][22][23], that appear as poles in the two-point retarded correlation function of the currents. In this paper we are interested to couple these quasi-normal modes of the plasma to an external electromagnetic field and study the possible electromagnetic modes supported by the system. It is important to observe that, even if the quasi-normal modes of the medium appear as poles of the correlation functions, due to high dissipation, the quasi-particle approximation is usually not valid. This implies that at any frequency the response function can not be approximated by a single resonance. Indeed, the self-frequencies of the system are in general complex and the amount of dissipation and propagation of the quasi-normal modes is usually comparable 8 .
The propagation and dissipation of an electromagnetic wave is described by the corresponding refractive index: the real part of the index accounts for the propagation of the wave while the imaginary part describes its attenuation. We focus here on the transverse modes of the electromagnetic field. The analysis of the longitudinal channel is contained in Section 7 and Appendix A. A study of the first three transverse modes (i.e. those associated to the lowest wave-numbers) reveals the rich structure of the electromagnetic response of the medium. We have analyzed the real and imaginary parts of their refractive index and of the wave-vector q as functions of the frequency, and, in order to characterize the propagation/dissipation rate of the aforementioned modes, we have studied the absolute value of the ratio Im(n)/Re(n) 9 .
The results are presented in Figures 1, 2 and 3. For reasons that will become clear in the next Sections we plot the results using the rescaled wave-vector q and the rescaled frequency w, obtained by dividing the wave-vector number and the frequency by the temperature of the system A quite interesting picture arises from our analysis of the electrodynamic response of the plasma. The system can indeed support various electromagnetic modes, some of them evanescent and other propagating, with different peculiar characteristics. The wave-vector of one of the modes vanishes for vanishing frequency while those corresponding to the other two modes reach a finite value at zero frequency (one purely real and the other purely imaginary). The presence of a finite value of q at null frequency (regardless of its real, imaginary or complex character) is due to the coupling of the electromagnetic waves to non-hydrodynamical quasi-normal modes of the plasma. The mode presenting an imaginary q at zero frequency is very dissipative in the low-frequency region. On the contrary, the one presenting a real q at zero frequency propagates with very small damping. It is very interesting to observe that the real part of the refractive index of this last propagating electromagnetic mode is negative implying that it is a negatively refractive mode, i.e. its energy flux and phase velocity are directed in opposite directions 10 .
In some previous studies [6,7] it has been proven that for charged plasmas in the hydrodynamic regime it always exists a mode with negative refraction for small enough frequency, while there is no negative refraction for globally uncharged plasmas. In the present paper we extended the investigation beyond the hydrodynamical modes and we proved that globally neutral, strongly coupled plasmas can support modes with negative refraction too. Moreover, in [8], the presence of multiple electromagnetic modes was shown to be generic for the hydrodynamical regime. Here we extended such analysis to entail also non-hydrodynamical modes and we showed that ALW still appear to be present. In this extended context it is natural to relate the presence of ALW to the coupling of the electromagnetic radiation to the infinite tower of quasi-normal modes.
Note that we studied the modes at frequencies w 1. In this regime a neat hierarchy between the modes arises from their dissipation over propagation ratio. From Figure 3 one can observe that the negative refraction mode dominates the signal propagation. In line with the collected numerical data one expects that, increasing the frequency, more electromagnetic waves propagate significantly. To have a faithful picture of the plasma response (in terms of propagating waves) at higher frequency it is therefore necessary 10 The system we are studying is a passive medium which can only dissipate and not source energy. This implies that the sign of the imaginary part of the refractive index agrees with the direction of propagation of the electromagnetic energy flux in the medium [24]. Intuitively, in a passive medium, a wave propagates in the direction in which it is damped.
to include higher modes into the analysis 11 . The negative refraction phenomenon of the plasma at hand is then likely to be significant in the low-frequency regime. Of course a careful study of the boundary conditions would be required to actually understand which modes are excited in a specific experimental setup. In this paper we concentrate on the bulk properties of the medium and we leave this problem for future investigation.
As a final comment, let us mention that in this Section we presented the results without making reference to the specific value of the temperature or other physical parameters of the system. Let us anticipate that the qualitative behavior of the system is robust and valid in a wide range of parameters. A precise treatment of the quantitative predictions of our model requires a careful analysis of the regularization/renormalization procedure which we postpone to later sections.

Stringy setup and gauge/gravity computation
In this Section we provide some details about the actual class of systems we consider, the computations we performed and the associated dual gravitational setup. The material collected in the first two Subsections is quite standard and describes how to perform computations in the gauge/gravity correspondence at finite temperature. Nevertheless, we decided to report it in some detail since some subtleties of the computation are relevant for our analysis concerning the electromagnetic modes. Experts in gauge/gravity correspondence could however speed through the discussion below and go to Subsection 4.3.
As we have previously stated we consider the electromagnetic response of strongly coupled neutral plasmas. The gauge/gravity correspondence allows us to describe the quantum dynamics at strong coupling of a gauge field theory defined on 4-dimensional Minkowski space-time in terms of a classical gravitational model living in a 5-dimensional Anti-de Sitter space-time (AdS 5 ). Furthermore, the gauge/gravity correspondence maps the finite temperature phase of a quantum field theory into a finite temperature gravity configuration: namely a black hole solution which becomes asymptotically AdS 5 at large distances (AdS black hole).
We do not restrict ourselves to any particular model, but we study instead the generic behavior of the universal sector of a strongly coupled theory containing the energy momentum tensor T µν and the U (1) conserved current J µ .
In the context of the gauge/gravity correspondence, these two minimal ingredients are mapped respectively to a 5-dimensional metric g mn and a 5-dimensional vector field A m , whose dynamics is encoded in the 5-dimensional Maxwell-Einstein action 12 where Λ is the cosmological constant and N is, in the gravity model perspective, a normalization constant. The gravitational constant has been related to the AdS 5 radius of curvature L and the string constant α which have been both put to 1.
The ground state at finite temperature T and zero charge density of a strongly coupled quantum field theory is dually described by the neutral AdS 5 black hole solution of the action (4.1) specified by the following bulk metric written in the Poincaré patch The horizon of the black hole is located at u = 1, while the AdS 5 conformal boundary corresponds to u = 0.
In the gauge/gravity framework, the study of the fluctuations of the 5-dimensional photon A m on the background specified by (4.2) provides information on the retarded 2point correlation functions of the U (1) current J µ of the dual (3+1)-dimensional quantum field theory [25]. As already stated, the knowledge of the retarded 2-point correlation function of the currents is the required element to discuss the electromagnetic properties of such strongly coupled plasmas.

Vector transverse fluctuations
We focus on the fluctuations of the 5-dimensional bulk photon A m . We fix partially the 5-dimensional gauge invariance imposing A u = 0; note that this 5-dimensional "radial" gauge condition still leaves the usual 4-dimensional gauge freedom for the remaining components of the vector potential A µ . In other words, we study the fluctuations of the bulk gauge field along the 4-dimensional space-time directions parametrized with µ = t, x, y, z where the dual theory is actually defined. Moreover, we consider the fluctuations which are transverse to the direction of propagation, namely the fluctuations that are perpendicular to the spatial wave-vector q.
Since the system is invariant with respect to spatial rotations, we can in general take the spatial momentum of the fluctuations to be along the z direction so that the 4-momentum is k µ = (ω, 0, 0, q). The transverse fluctuations are therefore those along the x and y directions and we introduce the label α = x, y to refer to this transverse space. Note that, as we have cylindrical symmetry around the z direction, we can focus on the gauge field fluctuations along x without spoiling the generality of the treatment. Specifically, we adopt the following ansatz for the Fourier components Plugging this ansatz into the Maxwell equation derived from the action (4.1) and considered on the background (4.2), we obtain which, by cylindrical symmetry, is valid for the fluctuations along a generic transverse direction. The u, q and ω dependence of φ is understood, the prime indicates the derivative with respect to the radial coordinate u; we also adopted the dimensionless frequency w = ω/(2πT ) and momentum q = q/(2πT ) already introduced in (3.1), where T is the temperature of the dual field theory 13 or, equivalently, of the AdS black hole in (4.2). A near-horizon analysis of Equation (4.5) yields the following asymptotic behavior for The generic solution to the differential equation is given by the superposition of an ingoing and an out-going solutions. We select the in-going solution imposing the condition a + = 0 in accordance with the idea that nothing can be emitted by the black hole horizon at the classical level. This choice corresponds to compute retarded correlators in the dual field theory. As we are confronted with a second order differential equation we need to impose one further boundary condition. Relying on the linearity of Equation (4.5) we choose a − = 1; as we will shortly explain, we are eventually interested in ratios like φ /φ which are completely insensitive to this specific choice. As we will describe in detail in Section 5, once equipped with these horizon boundary conditions, we are able to solve the differential problem propagating the solution from the horizon. We now turn the attention to the near-boundary region corresponding to small u. This asymptotic UV analysis is necessary to extract from our model the physical quantities we are interested in. A near-boundary term-wise study of the equation of motion (4.5) near the AdS boundary at u = 0 shows that the asymptotic expansion of the field φ is there where the φ i andφ i coefficients are independent of the u coordinate, and satisfy 14 (4.10) Note that, since the equation of motion is a second order differential equation, we can express all the coefficients of this term-wise solution as functions of the first two coefficients φ 0 and φ 1 .

Correlation functions and holographic renormalization
In this Section we follow the gauge/gravity prescription to compute the retarded twopoint correlator of the transverse currents of the strongly coupled plasma. Indeed, once we found the solution to Equation (4.5) with the boundary condition explained in the previous section, we can obtain the retarded correlator of the current [25]. More specifically, as it is usual in quantum field theory, one can derive the correlators through functional differentiation with respect to the sources; from the gravity model standpoint, the procedure corresponds to functionally differentiate the on-shell bulk action with respect to the boundary value of the fluctuating fields. The on-shell gauge field action for the transverse vector field (described by the solution φ of the previous Section) becomes (N T ) 2 16 This latter is to be integrated on the boundary manifold. Interestingly, the on-shell action reduces to boundary terms and only the contribution from the AdS conformal boundary at u = 0 is not vanishing. The contribution from the horizon at u = 1 actually vanishes because f (u) is there null. As a consequence, we just focus on the contribution to the primitive of (4.12) at asymptotically small u; we expand the field φ and its derivative as in (4.7) discarding the terms that vanish at u = 0 and we obtain It is important to observe that there exists a logarithmically diverging term. This signals the need to renormalize the model. The occurrence of a large-volume divergence in the gravity theory living on the asymptotically AdS geometry corresponds (through gauge/gravity correspondence) to the UV divergences of the dual quantum field theory. Let us describe the renormalization procedure [28] of the model at hand. At first we regularize the on-shell action considering a small u = cutoff, (4.14) We then add an appropriate boundary counter-term where γ represents the determinant of the metric induced by the bulk metric on the u = shell, the indexes i, j run only over the 4 boundary directions andc is a real numerical constant (we will later comment on this constant). Notice that the overall factor is chosen to express S c.t. in terms of the gothic variables w and q. The renormalized action is obtained from and, explicitly, we have where we have introduced c = 2(1 −c) for later convenience. We remind the reader that φ 0 represents the boundary value of the bulk fluctuation field which, through the gauge/gravity dictionary, is mapped to a source of the dual boundary theory. It is therefore with respect to φ 0 that we are interested in taking derivatives of the on-shell bulk action. Specifically, the renormalized 2-point correlation function for the transverse current is obtained taking the second order derivative of the renormalized on-shell action (4.17) with respect to φ 0 . In Fourier space we obtain (4.18) In the last passage we exploited the linearity assumption and the fact that for zero source φ 0 the fluctuation profile becomes trivial 15 . From (4.18) we can define the renormalized retarded correlator for the transverse current: the label (c) is a reminder of the dependence of the correlation function on the coefficient c in front of the contact term. The correlation function (4.19) is the fundamental quantity that allows us to compute the electromagnetic response of the strongly coupled plasma according to equations (2.3) and (2.4), as explained in Section 2. It is easy to check that the function defined in equation (4.19) satisfies the usual properties for a response function of a causal quantum field theory: namely it has poles in the negative imaginary part of the complex plane and its imaginary part is negatively defined when ω and q are real.

Contact term
As explained above, the computation of the current-current correlator requires the renormalization of the on-shell action. In particular, the renormalization demands to consider a counter-term to cancel the logarithmic divergence in (4.13). However, as usual in quantum field theory, the choice of the counter-term is not unique, and only its diverging part is specified by the renormalization procedure; the finite part should instead be fixed by some experimental measure. Such ambiguity is accounted by the real constant c in the definition of the correlator in (4.19), that is actually not fixed by any consistency requirement related to symmetries and Ward identities. The arbitrariness of the finite part of the counter-term introduces a term proportional to w 2 − q 2 in the correlator (4.19) which is usually referred to as a contact term 16 . Therefore, the constant c in the correlator should be in principle fixed with an experimental measure of an observable containing it.
Once the value of c is fixed, the model does provide quantitative results depending on N (a quantity which, roughly speaking, is associated to the number of degrees of freedom of the system) and the temperature T (which represents a physical scale of the model considered on the black hole solution). As we will see, it is particularly important to underline that the qualitative behavior of the electromagnetic modes in the medium, and more specifically of the refractive index, is a robust feature with respect to the actual value of c.

Electromagnetic modes
As explained in Section 2, the possible electromagnetic transverse modes supported by the plasma are obtained solving the dispersion equation 17 involving q and w: where e is the electromagnetic coupling constant as previously introduced in (2.3) 18 . Substituting the explicit expression of the correlation function obtained in (4.19), the wave equation becomes Note that, since our system is conformally invariant at zero temperature, as long as we work with the dimensionless quantities w and q, we have no explicit dependence on T in the wave equation. T actually provides only a scale with respect to which we measure the actual physical quantities. Furthermore it could be interesting to notice that in (4.20) the presence of a finite value of c can be reabsorbed by an overall rescaling of the correlation function; more precisely The rescaling (4.22) suggests that, as far as the study of the wave equation (4.20) is concerned, it is possible to trade the contact term with a normalization factor in front of 17 This is actually the wave equation for the transverse component of the electric field inside the medium. 18 In the gauge/gravity correspondence setting, J is the conserved current of a global U (1) symmetry of the quantum field theory. In order to obtain the electromagnetic response of the medium, the standard procedure is to introduce in the QFT action the interaction term eJA where A is the electromagnetic field considered as external and e is the associated electromagnetic coupling constant considered to be perturbatively small. The retarded correlator of the global current provides, at leading order in e, the exact result for the retarded correlator of the local current. the correlator, this latter being related to eN . We could therefore work without specifying the value of eN and reducing everything to a contact term to be fixed against a physical measure performed on the electromagnetic modes. One could for example think to fix c with the requirement that the value of q(w = 0) for a specified mode coincides with that measured in an experiment.
To check the soundness of our results we did a scan of the electromagnetic modes over a broad range of values of c and eN . The results show that the qualitative behaviour presented in Section 3 is the same for a very broad range of parameters inside the validity of the numerical computation. For convenience in the paper we decided to plot the results for the specific values: eN = 3 and c = 5.5.

Semi-numerical analysis and checks
In this Section we would like to explain briefly how we actually performed the computations and characterized the various electromagnetic modes supported by the strongly coupled plasma. We refer again to Figures 1, 2 and 3 in Section 3. The aim is to solve the wave equation (4.20) and find the dispersion relations q A = q A (w) (the label A distinguishes the different modes) connecting the complex wave-vector q and the real frequency w.
At first we solve the equation for the transverse vector fluctuations (4.5); once a solution is obtained, we can read the near-boundary coefficients and plug them into (4.19) to find the correlator. To solve the differential equation (4.5) we use two different methods whose results are later compared and cross-checked. The first method consists in integrating numerically the equation of motion (4.5) as explained in [26]. The second method is semi-analytical and consists in expanding the function φ(u) near the horizon (as in (4.6)) and solving Equation (4.5) order by order along the lines of [29]; this infrared solution is then matched with an analogous term-wise solution computed at the boundary, see (4.7). It turns out that, while the numerical method yields more precise answers, the semi-analytical method is instead more agile.
Once a solution to (4.5) is obtained, it is then possible to follow the gauge/gravity recipe and define a correlation function as described in Subsection 4.2. For very small values of w and q we check that the correlation function agrees with the analytical result of [30]. The correlator is then inserted in the electromagnetic wave equation (4.19) to search for the dispersion relations q A = q A (w) of the various transverse modes. As the correlator is a complicated rational function of q and w, the electromagnetic wave equation has in general many solutions. The set of possible solutions q A (w) can be studied relying on a simplifying approximation: We build the correlator using the previously introduced semi-analytical matched solution of (4.5); we then write the correlator as the ratio of two polynomials in w and q; we expand both the numerator and the denominator up to a suitable order in q and w around a chosen point in the (w, q) space (this statements will be shortly made more precise). The results we obtained with the aforementioned approximations are afterwards checked against the correlation function computed from the full numerical solution.
Let us be more specific on the computational procedure. We find solutions to the electromagnetic wave equation expressed as complex functions q A (w) where the frequency is a real quantity. We test the soundness of the solutions as follows: we first repeat the procedure with an increased order in w and q up to which we expand the numerator and the denominator of the correlation function; we than check that the solutions found previously (i.e. with shallower expansions) still solve the equation within the numerical tolerance. As an additional test we check that the solutions fall in the (w, q) region where the approximated correlation function adhere well to the correlation function obtained in the fully numerical approach. Eventually and more stringently we also check that the solutions obtained with the approximate method solve the electromagnetic wave equation written in terms of the fully numerical correlation function within the numerical tolerance.
We start looking for reliable electromagnetic wave solutions in the small w and q region (w, q 1); we then follow the modes iterating the approximate procedure and expanding around a generic point of their dispersion relation q(w). This allows us to analyze the modes at higher values of the frequency.

Phenomenology: wavelengths and the QGP
The analysis presented so far should be able to capture the electromagnetic response of a generic globally neutral strongly coupled plasma described by the gauge/gravity correspondence. A natural candidate for phenomenological application of the model at hand is the QGP studied in high energy physics experiments such as RHIC or LHC, where the plasma is dominated by the temperature, while the charge density is comparatively smaller. A similar phenomenological investigation has been already performed in [7] for the hydrodynamical modes of strongly coupled and charged plasmas. However in [7] it was observed that, even if the fluid dynamic analysis demonstrates the generic presence of negative refraction, the wavelength of the electromagnetic mode is larger than the dimension of the typical QGP sample produced in actual experiments.
In this Section, without pretending to provide any specific experimental or phenomenological prediction, we briefly investigate the wavelengths of the modes studied in the paper. In particular we show that the wavelengths of some non-hydrodynamical modes are smaller than the those of the hydrodynamical modes considered in [7] and that they are actually comparable with the typical dimension of the sample (namely few fm). This result supports the possibility that the coupling of the electromagnetic field with the non-hydrodynamical quasi-normal modes could in principle have phenomenological and experimental relevance and it calls for further investigation.
In Figure 4 we report the wavelengths of our electromagnetic transverse modes. The wavelength of a mode is given by λ = 2π/Re[q] = 1/T Re[q] in accordance with (3.1). To associate a physical order of magnitude to the wavelength of the analyzed modes we need to specify the temperature of interest. In relation to the QGP a typical value for the temperature is 200 MeV for experiments like RIHC or LHC 19 . Re-introducing the dimensionful physical constants, we find the following estimate for the wavelength ∼ 1 fm , (6.1) for frequency of the order of 10 24 Hz, as reported in Figure 4. In particular, note that the negative refracting mode of Figure 4 (solid line) presents a wavelength of order 1 fm for the full range of frequency. The size of the QGP samples in current experiments is actually of the order of few fm's [11,31,32], therefore (higher) non-hydrodynamical modes could actually probe the plasma and perhaps leave there their footprint. We leave a deeper exploration of the higher modes and their possible effects on the QGP for future work.

Longitudinal channel
In this Section we want to analyze briefly the longitudinal electromagnetic modes supported by the strongly coupled plasma and we proceed with a similar spirit as with the is satisfied. G L (ω, q) is the retarded correlator of the longitudinal current and, as in the transverse modes case, it is a rational function whose poles correspond to the longitudinal quasi-normal modes of the plasma. Equation (7.1) provides a set of dispersion relations q A = q A (ω) between the complex longitudinal wave-vectors q and the real frequency ω of the mode. In this Section we simply give an account of some of the results whose computational details can be found in Appendix A. In Figure 5 and 6 we plot the real and imaginary parts of the wave-vector 20 q, the real part of the longitudinal refractive index n = q/w and the ratio between the imaginary and the real parts of n for the first three electromagnetic modes supported by the strongly coupled plasma. These Figures should be compared with Figures 1, 2 and 3 in Section 3 which instead refer to the transverse sector. It is very interesting to observe that also the longitudinal channel supports various electromagnetic modes. However, we did not find any negative refractive longitudinal mode. All the modes have finite momentum at vanishing frequency. Nevertheless for two of them q(w = 0) is purely imaginary, while for the remaining mode is purely real. The modes having an imaginary momentum at vanishing frequency are highly dissipative in the IR regime even though the increase of the real part of q with the frequency indicates that they could become more propagating at higher frequencies, as it can be seen in Figure  6. The mode with purely real momentum at w = 0 is highly propagating already in the low-frequency region and it keeps this characteristic in the whole range of frequency that we considered. Notice that the sign of the real part of the momentum shows that all these modes have positive refracting index.

Conclusion
In this paper we studied the electromagnetic linear response of strongly coupled neutral plasmas described by the gauge/gravity correspondence and characterized the electromagnetic modes with the lowest wave-vectors. The salient features of the present analysis are the presence of multiple electromagnetic waves with different refractive indexes and a propagating negative refracting mode with very small dissipation. Our study has been performed without adopting any hydrodynamical approximation; hence we extended some previous hydrodynamical studies of strongly coupled plasmas beyond the regime of small frequencies and wave-vectors. Our simple model highlights the potential richness of the electrodynamic response of strongly coupled plasmas and this calls for further investigations. In particular, on a phenomenological level, we showed that the characteristic wavelengths of the electromagnetic modes in the plasma could be comparable to the typical size of the QGP samples produced in high-energy physics experiments. This fact provides some arguments supporting the possible relevance of the presented exotic phenomena for actual physical systems.
Our results suggest many future lines of investigation. Similar systematic analyses of the electromagnetic properties can be indeed performed using different kinds of gravitational backgrounds featuring appealing characteristics such as finite charge density, spontaneous symmetry breaking, non-relativistic or non-isotropic setups and the presence of magnetic fields. These extensions could find application in various physical systems such as the QGP but also in condensed matter and astrophysics. We know that strong spatial dispersion is a crucial ingredient in producing the exotic electromagnetic response we are concerned with, however it would be desirable to understand whether some more precise connection between the modes of the plasma and its associated electromagnetic modes can be clarified. In particular it would be nice to have a direct understanding of the presence or not of negative refraction for the various electromagnetic modes based on the QNM structure.
An interesting and ambitious aim would be to provide a systematic and complete analysis of all the non-hydrodynamical modes.
A in the dual gravity setup: A z . However in this case A z mixes with the time component A t [29]. We consider the following ansatz and define the gauge invariant combination: which represents the electric field in the z direction. The A t and A z equations then leads to a single equation for ψ, namely From a near-boundary study we find the asymptotic behavior of the field which is identical to that of φ written in (4.7). Recall indeed that in the limit of vanishing q the longitudinal and the transverse sector coincide. Moving to finite q does not affect the general form of the UV asymptotic expansions. The same UV asymptotic behavior for the longitudinal and traverse modes implies that the holographic renormalization procedure is analogous in both channels.We can indeed use the same counter-term c(w 2 − q 2 ) for both polarizations. The on-shell action for the combination (A.2) is: This is the longitudinal version of Equation (4.11). Hence the A z A z term of the on-shell action is (N T ) 2 16 We understand that, in a similar way as for the transverse channels, the zz currentcurrent correlator related to the second functional differentiation of the on-shell action with respect to A z can be expressed as follows These last passages are done in line with [29] to which we refer for further details. Note that the tt and tz correlators are obtained in a similar way leading to It is important to notice that the set of zz, tt and tz correlators satisfy the Ward identity k µ G (c) µν (w, q) = 0 . (A.10) This fact is related to the structure of the frequency and momentum dependent factors in front of the expressions of the correlators. As a consequence, the contact term proportional to c is not fixed by the Ward identities or, in other words, it is not constrained by symmetry requirements. All the arguments about the contact terms that we have described in relation to the transverse sector can be repeated for the longitudinal sector.

B Quasi-normal modes
In order to characterize better the strongly coupled plasma under study in the main text, we report here an analysis of its internal modes. The retarded correlation function accounting for the electromagnetic response of the plasma presents poles at specific values for the momentum and frequency of the external perturbation. Such poles correspond in the dual gravitational picture to quasi-normal modes of the black hole solution (see for instance [22,23,25,29,[33][34][35][36][37]). To find the dispersion relations between q and w of the above-mentioned quasi-normal modes it is enough to look at the zeros of the inverse of the correlation function. These solutions can be represented (for instance) as complex functions q(w) of the real frequency 21 . In this Appendix we perform a characterization of the first quasi-normal modes of the plasma both in the longitudinal and transverse sectors.
The results for the transverse sector are summarized in Figure 7. All the modes we found are qualitatively similar: they all show a finite imaginary part and null real part for q as the frequency goes to zero. Such circumstance corresponds to having a highly 21 It is important to underline that these dispersion relations between q and w are the ones associated to the quasi-normal modes of the plasma and they are different from the dispersion relations for the electromagnetic modes supported by the plasma and discussed in the main text. Indeed the first ones come from solving the equations: (G T,L (ω, q)) −1 = 0, while the second ones solve equations which look approximately as q 2 = ω 2 − G T (ω, q) and ω 2 − G L (ω, q) = 0. dissipative (and non-propagating) set of modes at low frequency. At higher values of the frequency, however, the real part of q increases until it becomes significantly bigger than its imaginary part. The modes are then propagating for higher values of the dimensionless frequency w. The longitudinal sector features different kinds of modes (see Figure 8): we have a diffusive mode whose complex momentum vanishes as the frequency goes to zero. This is the hydrodynamical longitudinal mode discussed 22 in [30]. The remaining modes are closely analogous to those we found in the transverse sector. They have a purely imaginary momentum q as the frequency vanishes. The numerical investigation leads us to imagine that we have an infinite tower of such modes which again are strongly dissipating at low frequency.
As a final comment we want to underline that whenever possible we checked our results with [23].