Holographic fundamental matter in multilayered media

We describe a strongly coupled layered system in 3+1 dimensions by means of a top-down D-brane construction. Adjoint matter is encoded in a large-$N_c$ stack of D3-branes, while fundamental matter is confined to $(2+1)$-dimensional defects introduced by a large-$N_f$ stack of smeared D5-branes. To the anisotropic Lifshitz-like background geometry, we add a single flavor D7-brane treated in the probe limit. Such bulk setup corresponds to a partially quenched approximation for the dual field theory. The holographic model sheds light on the anisotropic physics induced by the layered structure, allowing one to disentangle flavor physics along and orthogonal to the layers as well as identifying distinct scaling laws for various dynamical quantities. We study the thermodynamics and the fluctuation spectrum with varying valence quark mass or baryon chemical potential. We also focus on the density wave propagation in both the hydrodynamic and collisionless regimes where analytic methods complement the numerics, while the latter provides the only resource to address the intermediate transition regime.


Introduction
The physics of (2 + 1)-dimensional layers is at the core of two active research fronts in condensed matter: graphene and layered hetero-structures. The most notable examples of layered systems are the copper oxides and high-T c superconductors in general. Although the attention is usually focused on the physics along the layers, relevant information is also associated with the dynamics in the orthogonal directions [1]. This motivates the quest for theoretical models which can simultaneously capture the whole in-/off-plane behavior of multilayered systems.
The layers are in general expected to induce characteristic anisotropy both to the equilibrium and out-of-equilibrium properties of the system. However, it should be borne in mind that the anisotropy could also be regarded as an emergent effect due to the coarse-graining of the microscopic scales dominated by strongly-coupled fundamental degrees of freedom. Thus, to genuinely describe layered hetero-structures, the construction of a bottom-up model which just reproduces the lowenergy anisotropic physics may be unsatisfactory.
The description of strongly-coupled dynamics does not only require an alternative to perturbative approaches, but it can even radically affect fundamental schematizations on which intuition and physical models are founded. Most notably, the concept of quasi-particles is typically unsuited to account for the low-energy response of a strongly-correlated system, especially in a dense medium. In other words, the collective dynamics has to be directly addressed without any preestablished paradigm based on long-lived modes. This constitutes a fundamental motivation to resort to holographic descriptions [2,3] that, through the study of the small fluctuations of the gravity dual, provide explicit information on the low-energy collective behavior of the strongly-coupled boundary field theory. Layered systems are no exception. Layers introduce two qualitatively different "sectors" associated to the in-and off-plane physics, respectively. Accessing direct information about the collective low-energy response [4] is particularly important for layered systems. Recent experiments on copper oxides present puzzling outcomes, for instance the characterization of the modes emerging from the low-energy transport properties contrasts the featureless spectral density-density response [5].
In this paper we follow a top-down strategy, where we can, at least in principle, keep track of microscopic degrees of freedom. To describe the strongly-coupled physics of fundamental matter, we resort to AdS/CFT and consider a top-down intersecting D-brane construction. 1 The most natural intersecting D-brane configuration is the D3-D5-brane construction of [18,19]. In this setup the layer associated with directions along the (2 + 1)-dimensional D5-brane worldvolume introduces a co-dimension one defect in the ambient (3 + 1)-dimensional N = 4 super Yang-Mills theory. Previous studies of the D3-D5 systems where the flavor branes are considered in the probe approximation include [20][21][22][23][24][25][26][27].
In more detail, the large-N stack of space-filling D3-branes introduces adjoint matter living in 3 + 1 dimensions and a stack of D5-branes, partially wrapped along the internal manifold, introduce fundamental hypermultiplets on a R 1,2 ⊂ R 1,3 defect. However, we do not attempt to solve this localized configuration, but instead we introduce the defects through continuously smearing them over the orthogonal directions [28][29][30]. In addition to computational advantages, this procedure is also a physically appealing approximation. The smearing allows one to disentangle two consequences due to the layers: the breaking of rotations and that of orthogonal translations. Introducing the smeared distribution of layers, the interlayer separation is formally vanishing. This procedure then does not break orthogonal translations and leads to genuinely homogeneous but anisotropic configuration. In the dual ten-dimensional gravity side, the anisotropy is then encoded, e.g., in the resulting metric featuring a Lifshitz-like scaling and where the spatial direction orthogonal to the layers is singled out.
In this paper, we furthermore introduce flavor D7-branes [31] in the anisotropic D3-D5 background. We treat these D7-branes in the probe approximation, which on the dual field theory side corresponds to the partially quenched approximation, where the defect "sea quark" degrees of freedom stemming from D5-branes are unquenched while the "valence quarks" (D3-D7 strings) are quenched. We allow finite bare valence quark masses as well as a non-vanishing baryonic chemical potential for them, though we do not consider them simultaneously.
In addition to studying the thermodynamic phase diagram of the model we also discuss the excitation spectrum. In particular, we study the various propagating (sound) and purely dissipating (diffusion) regimes [32] of the longitudinal modes, both along and perpendicular to the defects. We especially focus on the peculiarities of the off-plane sector and comment on its analogies to Lifshitz models [33,34]. As for future directions, the characterization of the dynamical polarizability of the layered medium is a necessary step in view of considering semi-holographic generalizations [35]. Assigning a dynamical character to the boundary conditions of the bulk gauge field [36][37][38][39][40][41][42], one can address physically relevant questions, for instance: the propagation of light-waves through the layered medium [43][44][45] and the introduction of long-range (Coulomb-like) interactions relevant for plasmonic physics [46][47][48][49][50][51][52][53].
The paper is structured as follows. In Section 2 we describe the D3-D5-brane setup and briefly review the zero-temperature background and the corresponding black hole solutions. We also introduce a probe D7-brane and study its embedding. In Section 3 we start analyzing the thermodynamics of the background at vanishing chemical potential. We first recall the supersymmetric solution at zero temperature and show that it can be treated fully analytically. We then turn to a finite temperature analysis and discuss the two competing phases and the associated meson melting phase transition. The two phases correspond to Minkowski and black-hole embeddings of the D7-brane. In Section 4 we start analyzing the system at non-zero chemical potential. This allows us to approach questions that are relevant for the condensed matter systems. The Minkowski embeddings describe insulators while black hole embeddings are associated with metallic behavior. We will, however, mainly focus on the latter in this paper. We are especially interested in the fluctuation spectrum, the propagation of modes and their dissipation. In particular, we study the longitudinal zero sound both along and orthogonal to the layers. Adopting the two alternative real-momentum and real-frequency approaches, we characterize the response of the system. Section 5 contains our discussion and outlines directions for future research. The appendices contain several technical details omitted in the body of the paper.

Setup
In this section we review the background generated by intersecting D3-and D5-branes [18,19,29] to which we then add an extra D7-brane. The flat space and low-energy configuration of the set of all the branes can be summarized in the following table: The study of the complete dynamical system is rather involved, we therefore adopt an approximation in which the D7-brane is considered a probe in the geometry generated by the set of intersecting D3-and D5-branes. On the field theory side this corresponds to a partially quenched approximation where the fundamentals dual to strings attached to the D7-brane are not dynamical.
The geometry following from stacks of (coincident) D3-and (delocalized) D5-branes was obtained at zero temperature in [29], with subsequent generalization to a black hole geometry in [30]. The solutions to the field equations were obtained by the smearing approximation [28], in which the D5-branes are distributed homogeneously along the internal directions orthogonal to the D3's, as well as along one of the directions parallel to the D3's. The D5-branes can then be viewed as multiple parallel layers that create a codimension-one defect in the (3+1)-dimensional theory living on the worldvolume of the stack of D3-branes. This theory is very rich and has lots of interesting features [54]. The present analysis focuses on the extra D7-brane probing the multilayer geometry and introducing new flavor degrees of freedom in the fundamental representation of the gauge group (quarks). Similar bulk constructions [55][56][57][58][59][60] but with a completely (2+1)-dimensional field theory dual are connected with ABJ(M) Chern-Simons matter theories [61,62] coupled with fundamentals [63,64].
Let us now review in detail the D3-D5 geometry of [29,30]. The ten-dimensional Einstein frame metric at finite temperature can be written as: where ds 2 5 is the following five-dimensional metric: where φ is the dilaton. The anisotropy in the spatial directions in (2.3) is due to the fact that the D5-branes are extended along x 1 x 2 and smeared along x 3 . The warp factor h and the emblackening factor B in (2.3) are given by: The radius of curvature R is determined by the number of colors N c , i.e., the number of D3-branes: where Q c ∼ N c (see below). The dilaton field φ for this geometry is non-trivial and determines the spatial anisotropy: where Q f ∼ N f , with N f being the number of D5-branes per unit length in the x 3 direction. The precise relations between Q c , Q f , and N c , N f are where we are working in units in which g s = α = 1. The metric dŝ 2 5 of the internal part is [29]: In (2.8) χ and τ are angular coordinates taking values in the range 0 ≤ χ ≤ π/2 and 0 ≤ τ ≤ 2π, whereas ω 1 , ω 2 , and ω 3 are left-invariant SU (2) one-forms (see (A.3) in Appendix A for their explicit representation in terms of three angles). The complete type IIB supergravity background contains Ramond-Ramond five-and three-forms, whose explicit expression is not needed in this work (see [29,30]). The different fields satisfy the equations of motion of ten-dimensional type IIB supergravity with D5-brane sources. When r h = 0 (B = 1) our solution is supersymmetric and invariant under a set of Lifshitz-like anisotropic scale transformations in which the x 3 coordinate transforms with an anomalous exponent z = 3 [29]. Let us recall that the temperature T of a black hole is given in terms of the tt and rr components of the metric by the general formula, leading to: (2.10) Using (2.5) we can recast this relation in terms of Q c as: Let us now add a D7-brane probe to the D3-D5 geometry. According to the array (2.1), we extend the D7-brane along (x 0 , x 1 , x 2 , x 3 , r) and the three angular directions of the three-sphere corresponding to the one-forms ω 1 , ω 2 , and ω 3 . In addition we can consistently set τ = constant and consider the embedding χ = χ(r) . (2.12) Since we want to add charge in the fundamental representation in the dual theory, we consider a D7-brane with a potential one-form A on the worldvolume: Therefore, in order to fix completely the configuration, we have to determine the functions χ(r) and A t (r) by solving the equations of motion derived from the DBI action: where F = dA and the factors of e φ result from working in the Einstein frame. 2 The induced metric on the worldvolume for our ansatz is: where χ = dχ/dr. The integrand of the DBI reads where θ is the polar angle used to represent the one-forms ω i (A.3). The Lagrangian then is: from where one can derive the equations of motion for the embedding and the gauge potential. The equation of motion for A t can be integrated once, giving where d is a constant of integration which is proportional to the charge density. It is rather straightforward to solve (2.18) for A t , namely: We can now write the equation of motion for the embedding function χ(r) and use (2.19) to eliminate the worldvolume gauge field in favor of the density d: Notice that (2.20) is trivially satisfied by taking χ = 0. This solution corresponds to the so-called massless embedding which we consider in Section 4. Let us now discuss the asymptotic UV behavior of the solutions to (2.20), i.e., at large r. In this regime we can safely assume that χ is small, giving us the following differential equation An ansatz χ ∼ r α leads to an algebraic equation to be solved for the exponent α: which has two solutions .
(2.23) Therefore, deep in the UV, the angle χ behaves as: (2.24) The first term corresponds to the leading UV behavior. According to the holographic dictionary, its coefficient A is related to the mass m q of the quarks introduced by the D7-brane. Moreover, the coefficient B of the subleading solution determines the vacuum expectation value of the quarkantiquark bilinear operator O m =ψ ψ + · · · (see below). In the sections that follow, we analyze different particular solutions of (2.20). We begin by analyzing in Section 3.1 the case with T = d = 0, for which we find a simple analytic solution. In Appendix A we further show that this embedding is supersymmetric.

Vanishing density 3.1 SUSY configuration
Let us consider the D3-D7 background at zero temperature and suppose that the gauge potentials on the worldvolume of the probe D7-brane are vanishing. In this d = 0, B = 1 case the embedding equation (2.20) can be solved analytically. Indeed, one can directly verify that χ = χ(r) given by: where m r is a constant proportional to the mass (discussed at length below), is a solution of (2.20). Notice that (3.1) behaves in the UV as in the general equation (2.24) with B = 0. This means that the field theory dual to this solution has vanishing quark-antiquark condensate O m and that the constant m r in (3.1) is proportional to the quark mass. Interestingly, one can verify that (3.1) is the general solution of a first-order differential equation, which in turn can be obtained as the saturation condition of a BPS bound for the action. In order to prove this statement, let us write the Lagrangian density of this d = 0, B = 1 case as: where we have omitted the irrelevant prefactor. Let us now write the Lagrangian (3.2) for a general function χ(r) as the following square root of a sum of squares: where Z and Λ are given by: which, taking into account that Z is a total radial derivative, implies the following bound for the action S 0 = drL 0 : This lower bound only depends on the boundary values of the fields and is saturated when Λ = 0, i.e. when the following first-order equation holds: whose integration yields precisely (3.1). Notice that the previous argument shows that (3.1) minimizes S 0 for fixed boundary conditions at the UV. Thus, it solves the variational problem and, therefore, it must fulfill the Euler-Lagrange equations derived from L 0 . Actually, our derivation of (3.7) suggests that (3.1) is a supersymmetric worldvolume soliton. We check this fact directly in Appendix A by using the kappa symmetry of the D7-brane probe.
In this paper we focus mostly on non-supersymmetric configurations. However, the supersymmetric configuration (3.1) is quite useful in these studies since it allows us to regulate the thermodynamic functions of the probe at non-zero temperature. For this purpose, it is interesting to notice that the on-shell action for the BPS solution diverges as: as can be easily derived by evaluating the right-hand side of (3.6) on the solution (3.1) for r = r max → ∞.

Finite temperature
In the rest of this section we restrict ourselves to the case in which the charge density d still vanishes while keeping the temperature T finite. We want to study the thermodynamics of the probe D7-brane. For that purpose it is quite convenient to introduce a new isotropic radial variable u, similar to the one defined in [56,65]. We define u by means of the following differential relation with r: dr This equation can be readily integrated as: In this new variable the horizon is located at u = 1. Notice also that asymptotically in the UV: Eq. (3.10) can be easily inverted as: wheref (u) is defined as:f (u) = 1 + u − 10 3b = 1 + u − 15 4 . Let also define f (u) as: (3.14) The emblackening factor B in terms of f andf can be written as: Let us now use these results to write the 10d metric in the variable u. The internal metric dŝ 2 5 is the same as in (2.8). The non-compact part ds 2 5 takes the form: with the dilaton φ given by: Let us consider now an embedding ansatz for the D7-brane in which τ is constant and χ = χ(u). The induced metric is: whereχ = dχ/du. Since we are considering the case in which the charge density vanishes, the worldvolume gauge field is also zero. Therefore, the Lagrangian for the probe in the new radial variable is: from which the equation of motion can be readily obtained, Let us now study the UV behavior of the solutions of (3.20). By expanding χ in power series around u = ∞ we can solve (3.20) asymptotically. We get: where m and c are related to the quark mass and condensate, respectively. The precise relation is worked out in detail in Appendix B. It turns out that the bare quark mass m q can be written in terms of the parameter m as  where the coefficient c is a pure number in our units (B.21). As in [65], given the UV limit (3.21), there are two types of embeddings depending on how the brane behaves at the IR. At low temperature (or large mass parameter m) the D7-brane probe closes off outside the horizon and we have a so-called Minkowski embedding. On the other hand, if the temperature is large enough (small m) the probe terminates at the horizon and we have a socalled black hole embedding. For some intermediate value of m there is a phase transition between these two types of configurations. Due to the different boundary conditions they satisfy at the IR it is quite convenient to choose different variables to analyze these two types of embeddings. We do it separately in the two subsections that follow. In Appendix C we study the critical embeddings close to the transition, while in Appendix D we analyze the thermal screening of the quark-antiquark potential.

Black hole embeddings
In order to study the black hole embeddings of the probe, let us introduce the variable η, defined as: We parametrize our embeddings by a function η = η(u). Sincė the induced metric on the worldvolume of the D7-brane for this parametrization becomes: It follows that, in these variables, we have for the DBI determinant: Let us now obtain the reduced action of the probe, denoted by I bulk , defined as: where S is the action, V 4 is the (infinite) volume of the 4d Minkowski spacetime and N is: The explicit form of I bulk is given by Taking into account (3.24) and that sin χ ≈ χ − χ 3 /6 + . . ., we obtain from (3.21) the UV behavior of η: This behavior can be directly obtained by solving the equation of motion derived from (3.30) as a power expansion around the UV. Notice that the term u −3 is not present in the UV expansion of η. In order to get the embedding function η(u) in the whole range of the holographic coordinate u we must solve numerically the Euler-Lagrange equation of motion derived from the action (3.30). For a black hole embedding we must impose regularity conditions at the horizon u = 1, which are easily seen to be: where η h is an IR constant which determines the UV constants m and c. By varying η h we get a relation c = c(m), which determines the condensate as a function of the quark mass at fixed temperature or, equivalently, the condensate as a function of the temperature at fixed m q . The corresponding results are plotted in Fig. 1. We notice that c is negative for small m (or large T ), vanishes near m ≈ 1 and becomes positive for larger values of m. For high temperatures (or low m) the embedding function η remains small for all values of the holographic coordinate u and the equation for the embedding can be linearized and solved analytically. This analysis is presented in detail in Appendix E. We find that c decreases linearly with m (c ≈ −0.198m, see (E.22)). This result is in good agreement with the numerical results, as shown in Fig. 1.
We want now to address the problem of determining the thermodynamic functions of the probe. The first step in this analysis is finding the free energy density F of the D7-brane. The expression of F can be obtained from the Euclidean on-shell action of the brane, which is divergent and should be regulated appropriately by adding a suitable boundary term. Indeed, plugging the UV behavior (3.31) into the right-hand side of (3.30), we find that I bulk diverges in the UV as: It is quite instructive to compare (3.33) to the behavior (3.8) for the SUSY embedding at zero temperature. The three terms in (3.8) correspond to the first, third, and fourth terms in (3.33) (take into account that r max ∼ u 9/8 max ). The second term in (3.33) is independent of the embedding and is due to the finite temperature.
To construct a boundary action that cancels the divergences of I bulk we first consider the boundary metric at u = u max It follows that: We want the renormalized action to vanish for the SUSY embeddings at T = 0. To fulfill this requirement we notice that the on-shell action of such embeddings depends on η as cos 4 χ = (1−η 2 ) 2 evaluated at the boundary (see (3.6)). Taking this fact into account, let us write the boundary action as: where α, β 1 , and β 2 are constants to be determined by imposing the condition that I bulk + I bdy is finite as u max → ∞. We have included in (3.36) a prefactor containing powers of the functions f andf to take into account the redshift occurring at T = 0 due to the emblackening factor of the metric. Actually, this term is needed to cancel the second term on the right-hand side of (3.33). Plugging in the UV behavior of the embedding function η (eq. (3.31)), we find By requiring the cancellation of the leading term of I bulk + I bdy , we fix the value of the constant α to be: Moreover, the cancellation of the u 3 4 max term requires that: Using these values of the constants, we can write I bdy as: Notice also that taking the solution β 2 = −β 1 = 1 of (3.39), the boundary action can be written in terms of the emblackening factor B as: Let us rewrite the divergent terms on the right-hand side of (3.40) as: Using this result, we can write I bdy as the following integral: The free energy density of the probe is given by: where G(m) is the following integral: We have extended the upper limit of integration to ∞ since it is a convergent integral. For black hole embeddings, u min = 1 in (3.44) and (3.45).
Once F is known, we can obtain the other thermodynamic functions by using the standard relations between them. For example, the entropy density s of the probe can be computed as: Since N ∼ T 4 , ∂ T N = 4N /T and therefore Let us now calculate the second term in this expression where we took into account that, for fixed quark mass m q , the mass parameter m behaves as m ≈ T −b (B.8). The derivative of F/N with respect to m can be extracted from (B.20), yielding where F/N can be obtained from (3.44) for black hole embeddings or from (3.66) for Minkowski embeddings. The internal energy E is given by E = F + T s: Let us next compute the heat capacity, defined as: Proceeding as with F to compute this derivative, we get: From the behavior of F at high temperature obtained in Appendix E, (eq. (E.36)), it is immediate to find the leading behavior of s, E, and c v : We have checked this high temperature behavior numerically (see Fig. 2).

Minkowski embeddings
The η = η(u) parametrization is convenient for the black hole embeddings. For Minkowski embeddings, Cartesian-like coordinates (ρ, P ) are better suited. They are defined in terms of u and χ as follows: The inverse relation is: When using these variables, we parametrize our embeddings by a function P = P (ρ), then where P = dP/dρ. Let us now write down the induced metric on the worldvolume. Taking into account that: as well as: we get: (3.60) In this expressionf and f should be understood as the functions: It follows that the DBI Lagrangian density in these variables is: Let us now study the free energy F in these variables. First of all, let us write the bulk on-shell action as: This action diverges at the UV. To characterize this divergence, we use the UV behavior of the embedding function P (ρ): Following similar steps as in Section 3.2.1, it is straightforward to find the regulator: Alternatively, one can simply perform the change of variables on (3.43) to land on (3.65). It follows that the free energy F can be represented as: As a highly non-trivial check of this expression, we have verified numerically that F vanishes when m → ∞ (see Fig. (2)). The other thermodynamic functions S, E, and c v can be obtained from the free energy F by using

Non-zero density
Let us now explore the D3-D5-D7 system at non-zero charge density. To simplify our analysis we restrict ourselves to the case in which the embedding function χ is zero which, according to the holographic dictionary, corresponds to having massless quarks. Recall that χ = 0 solves trivially the equation of motion (2.20). Moreover, when χ = 0 eq. (2.19) yields the following value for the radial derivative of the worldvolume gauge field potential A t : (4.1) By the AdS/CFT dictionary, the chemical potential is given by the UV value of A t : The last equality follows for the black hole embeddings of the D7-brane, since the temporal part of the background gauge potential has to vanish at the horizon. Notice that in the current case, at finite density, the only possible embeddings are those that enter into the black hole. At zero density, but at finite chemical potential, there is the possibility of considering Minkowski embeddings in which case one needs to allow non-zero IR values of the gauge potential [66,67].
Using the value of the dilaton φ for our background (2.6), we can directly perform the integration in (4.2), yielding where γ is the following constant: For later purposes we also define a function H as follows (4.5)

Stiffness at zero temperature
Before discussing the excitation spectrum of our system it is useful to start with computing the stiffness at zero temperature. This is given by the following thermodynamic derivative at fixed entropy density where p i is the pressure along direction x i . At zero temperature this typically agrees with ∂p i ∂ T [68] that we aim to compute. With an abuse of language we follow the existing literature and henceforth call (4.6) the speed of first sound.
Let us now take T = r h = 0 in our equations and let us examine the thermodynamics of the probe at zero temperature. When r h = 0 the second term in (4.3) vanishes and µ is related to the density d as follows: Then, the on-shell action is: where V 4 is the volume of our system in the four-dimensional Minkowski space-time. Since the function H is given by: the on-shell action is divergent and must be regulated by subtracting the action at zero density. We get: (4.10) The grand potential Ω is just Ω = −S reg on−shell . Therefore: From Ω we get the density ρ as: as well as the energy density : In order to obtain the pressures along the different directions, let us put the system in a box of sides L 1 , L 2 , and L 3 in the directions of x 1 , x 2 , and x 3 respectively. Then, the pressure along the i th direction is given by: The grand potential Ω is an extensive quantity which clearly depends linearly on L 1 and L 2 . Therefore, if p refers to the pressure in the x 1 x 2 plane, we have: It follows that the corresponding in-plane speed of sound is: In section 4.2.1, by studying the fluctuations of the probe D7-brane, we will show that this value of u coincides with the value of the speed of the in-plane zero sound. The result in (4.16) matches the value corresponding to a conformal theory in 2d and not that of 3d. Notice that at finite chemical potential there are plenty of examples in holography where one can exceed the conformal value [68][69][70], as long as the conformal symmetry is broken. Here the mechanism is quite different and can be associated with the presence of the defect instead of some intrinsic energy scale. At vanishing chemical potential there are also other mechanisms that lead to stiff equations of state: brane intersections with non-AdS asymptotics [71][72][73], softly broken translational symmetry [74], having non-relativistic scaling symmetries [75,76], dynamical magnetic fields [77], and those obtained from double trace deformations [78].
The dependence of Ω on L 3 is not clear since Q f is proportional to the density N f of D5-branes smeared along x 3 . Therefore, the calculation of the thermodynamic derivative (4.14) is far from obvious in this case. Actually, we will verify in section (4.2.2) that the dispersion relation of the off-plane modes in the zero sound channel is quite different from the ones corresponding to the sound modes propagating along the plane directions.

Zero sound
Let us now turn to discussing the collective phenomena that arise from fluctuating the fields on the worldvolume of the D7-brane. In the field theory dual, these excitations are density waves which correspond to quasinormal modes on the gravity side. We are only discuss vector fields, focusing solely on s-waves. Thus, we suppress the dependence on the internal part of the geometry. We start the analysis at zero temperature, where we can obtain analytic results. In the subsequent section, we instead focus on the high temperature result, which is also amenable for analytic treatment. In between there is the regime of finite temperature where we need to resort to numerical analysis, which accurately interpolates the limiting cases.
Let us thus begin by considering our system at exactly zero temperature. Here the dominant mode is the so-called holographic zero sound [79] which is a pretty robust collective excitation mode persisting to finite temperature [80]. Many other aspects of the holographic zero sound mode have been discussed in various contexts [16, 42, 71-73, 75, 76, 81-100]. To study these modes we perturb the worldvolume gauge potential as follows where A (0) is the gauge field of (4.1) and δA is a small perturbation about the equilibrium value. We Fourier transform to momentum space: The equations of motion for the fluctuations can be obtained by perturbing the DBI action around the configuration with χ = 0 and A = A (0) . The derivation of these equations is done in Appendix F. We work in the radial gauge by setting a r = 0. We also introduce a gauge-invariant combination, the electric field: E = k a t + ω a . In what follows we consider the case in which a and k are parallel. There are then two possibilities depending on the direction of wave propagation. When k is oriented along the (x 1 , x 2 ) plane we call it "in-plane", whereas the "off-plane" waves propagate along x 3 . We distinguish between these two cases in the subsequent sections.

In-plane zero sound
Let us consider the in-plane fluctuations at zero temperature. For concreteness, we assume that k = k e x 1 is directed along x 1 ≡ x. The corresponding equation for E ≡ k a t + ω a x is worked out in the appendix, c.f. (F.24), and takes the following form Since in this section we are interested in the T = 0 case, we take B = 1 from now on. We first study (4.20) near the horizon r = 0, where it takes the form: The solution of this near-horizon equation with infalling boundary conditions at r = 0 can be written in terms of a Hankel function of the first type: is not an integer the Hankel function has the following expansion near the origin: Let us then expand the near-horizon electric field at low frequency as: where C is a constant: Let us now perform the two limits (near-horizon and low frequency) in opposite order. At low frequency and momentum we can neglect the last term in (4.20) and the equation can be integrated once to give: where c E is a constant. To perform a second integration, let us define the rescaled densityd as follows: where α is the quantity defined in (4.9). We also define the integrals I(r) and J(r) as: (4.28) Then, E(r) at low frequency can be written as: where E 0 is the value of E at the UV boundary: The integrals I(r) and J(r) are special cases of the integrals that we define and inspect in the Appendix G, see (G.2) and (G.5), respectively. More precisely, the relationships are: According to the general formulas (G.3) and (G.6) these integrals behave near r = 0 as: where γ is the constant defined in (4.4). Therefore, E(r) behaves near the horizon as: Let us now compare (4.24) and (4.32). By looking at the terms depending on r we get: By comparing the constant terms and using (4.33) we find: By imposing Dirichlet conditions at the UV boundary, i.e. E 0 = 0, we get the following dispersion relation: In Fig. 3, where we have used the reduced quantities defined in (4.79), we compare this to the numerical results and find that they are accurately captured at small temperature. We choose to represent the dispersion in the case where we keep frequency real and take complex momentum. Then a convenient way to display the result is to plot the ratio which readily follows from (4.35). In the limit of high frequency this ratio asymptotes to √ 3. We note that at any non-zero temperature, the zero sound is not the dominant mode at small frequency, but the relevant physics is dominated by the diffusion pole discussed in the subsequent section.
Let us then study the dispersion relation at small frequency. At leading order in ω (4.35) becomes: We note that this result equals the first sound result (4.16) obtained in the previous section. Note also that we could write the dispersion relation (4.35) in terms of the chemical potential µ by using the relation: Let us now find the next order dispersion relation by writing and plugging this ansatz in (4.35), we get at first order in δω: We can now solve for δω in powers of k . The second term on the RHS is clearly of higher order. Therefore: The imaginary part of this equation gives the attenuation of the in-plane zero sound, namely: Moreover, the higher order correction to Re ω is:

Off-plane zero sound
We now analyze the case in which k = k ⊥ e x 3 lies in the direction of x 3 ≡ z. The equation for the gauge-invariant combination E = k ⊥ a t + ωa z has been obtained in (F.28) and can be written as: We first analyze (4.44) with B = 1 near the horizon r = 0, where it takes the form: The solution of this near-horizon equation with infalling boundary conditions is: In order to expand E at low frequencies, we use the expansion of the Hankel function of index zero near the origin: where γ E is the Euler-Mascheroni constant, yielding Let us now redo the computation taking the limits in the opposite order as in the previous subsection. At low frequency we neglect the last term in (4.44), which leads to an equation that can be integrated once as: where c E is a constant. A second integration yields: where α has been defined in (4. .

(4.51)
For small r this function behaves as: Using this result, together with the expansion of J(r) written in (4.31), we get: Let us now match (4.48) and (4.53). First we identify the terms containing log r, which gives: Using this result, we identify the constant terms and write the UV value of E as: By requiring Dirichlet boundary condition E 0 = 0 we arrive at the desired dispersion relation. To simplify the final expression we define the constant D as: so that the dispersion relation can be neatly written as In Fig. 3 we compare this to the numerical results and find that they are accurately captured at small temperature. Again, we have chosen to represent the dispersion in the case where momentum is complex while frequency is real. The ratio of the real and imaginary parts of the momentum in this case reads (4.58) In the limit of high frequency this ratio asymptotes to zero. Let us also study the small frequency limit of the dispersion relation just obtained. At leading order in small ω (4.57) becomes:

Diffusion modes
We now consider the fluctuation modes at non-zero temperature in the diffusive channel, i.e., solve exactly the same fluctuation equations as before but consider the other limiting case, that of very high temperature. Thus, the corresponding equations for these modes are (4.20) and (4.44) for the in-plane and off-plane diffusion, respectively. These two cases are considered separately in the two subsections that follow. The goal is to find diffusive modes in which ω is purely imaginary and related to the momentum as ω = −i D k 2 + . . ., with D being the diffusion constant. In what follows we get an analytic expression of D for the two types of waves.

In-plane diffusion
Let us begin our analysis by expanding the in-plane fluctuation equation (4.20), with B = 1, around the horizon r = r h . The expansion of the emblackening factor B is: The functions multiplying E and E in (4.20) can be expanded as: where the constants A, c 1 , and c 2 are given by: Plugging these expansions in (4.20), we get the following near-horizon fluctuation equation: We solve this equation in Frobenius series near the horizon: Therefore, also at leading order, we have:  Notice that β ∼ O(1). Moreover, as a ∼ ω ∼ O( 2 ) we can take a = 0 in (4.64) at leading order and write: with E nh = E(r = r h ).
We now perform the limits in the opposite order. First, we take the limit of low frequency. The fluctuation equation (4.20) takes the form: which can be integrated as: where c E is an integration constant. Let us now expand (4.71) near r = r h . We have: where I is the integral: Clearly, E nh = E (0) + c E I . By imposing the Dirichlet boundary condition E (0) = 0 at the UV boundary r → ∞, we get: Moreover, by comparing the linear terms in r − r h in (4.69) and (4.72) we obtain the following relation between c E and E nh : Let us now introduce the reduced quantitiesd,ω, andk as: If we change in (4.20) to a new radial variabler = r/r h , the resulting equation expressed in terms ofd,ω, andk does not contain the constants R, α, and r h . The dispersion relation in this reduced frequency and momentum takes the form: At high T the reduced densityd becomes very small and the hypergeometric function in (4.82) is approximately equal to one. At leading order, we have: Taking into account the relation betweenD and D written in (4.82), we get that at large T the parallel diffusion constant behaves as: Let us now study the opposite regime of small T (larged) and we can approximate hypergeometric function in (4.82) as ThereforeD can be approximated as: Using now that D ∼ T −1D and thatd We have solved (4.20) numerically and checked that, indeed, for high enough temperature the diffusive mode is the dominant one. This numerical analysis allows us to extract the diffusion constant and to compare the result with the analytic formula (4.82). This comparison is performed in Figure 4, where we see that the agreement between (4.82) and the numerical results is very good. When the temperature is low enough, the zero sound mode is the dominant one and the system enters a collisionless regime. This collisionless/hydrodynamic crossover is illustrated in Figure 4, where we notice that the zero sound persists at T = 0 if T is small enough. When T is increased, the imaginary part of the zero sound grows as T 7 3 until the crossover to the hydrodynamic diffusive regime takes place. The frequency ω cr and momentum k cr at which this transition occurs depend on the temperature and chemical potential. Our numerical analysis has allowed us to determine that ω cr and k cr scale with T and µ as: (4.88) In Figure 5 we display typical dispersion relations of several excitation modes for the in-plane case. The off-plane dispersion relations are qualitatively the same. We notice that the zero sound dispersion saturates at high frequencies. This is due to interactions between complex modes, as is clear from the figure. It would be tempting to identify this in-plane zero sound mode with a surface plasmon polariton, whose dispersion relation has close resemblance. However, we cannot clearly separate the underlying physics and associate the saturation directly with surface phenomena due to the following observation: We found out that for any flavor Dq-brane configuration that we have tested, the zero sound saturates if one reaches high enough frequency. This in itself is remarkable and had gone unnoticed in all the previous works in this field. It would be very interesting to understand why this happens and in particular discern if this phenomenon is due to non-linearities of the DBI action.

Off-plane diffusion
Let us now study the off-plane fluctuation equation (4.44) in the diffusive regime. We first expand the coefficients of E and E near the horizon as: whereÃ,c 1 , andc 2 are given by: Following the same procedure as in the in-plane case, we next expand E(r) at low frequency and write E(r) near the horizon as in (4.69), where now β is given at leading order by: h α d 2 + α r 16 3 h .

(4.91)
Taking the limits in the opposite order, we arrive exactly at (4.70). Thus, we can expand E(r) as in (4.72) with I being the integral written in (4.73). The constant c E can be obtained from the consistency of both expansions. This procedure yields (4.92) and, by again imposing Dirichlet boundary conditions at the UV we get the diffusion dispersion relation ω = −i D ⊥ k 2 ⊥ + . . ., where D ⊥ is given by: 16 3 h . (4.93) Let us now reintroduce reduced quantitiesd,ω, andk ⊥ which allow to absorb R, r h , and α in (4.44). Due to different scaling in the off-plane direction one needs to be a bit more careful with the momentum. Thed andω are defined as in (4.79), whereask ⊥ must be defined as: Using these definitions the dispersion relation for the off-plane diffusion can be written asω = −iD ⊥k 2 ⊥ + . . ., whereD ⊥ is related to D ⊥ as: .

(4.95)
Moreover, from the expression of D ⊥ in (4.93) we get thatD ⊥ is given by: which means thatD ⊥ =D . The fact that the reduced, non-physical, in-plane and off-plane diffusion constants can be made equal by a temperature-dependent rescaling is a reflection of the scaling symmetry of the background metric. However, we should emphasize that the physical diffusion constants D ⊥ and D behave rather differently due to the different scalings used to define the corresponding hatted quantities. For example, at high temperatureD ⊥ =D ≈ 3 4 and thus: which means that, contrary to D , the off-plane diffusion constant grows as T

Discussion
Holography offers a possibility to investigate various phases of dense matter in the strong coupling regime. Intense effort has been put to study inhomogeneous phases, 3 starting with the early works of [102,103] with the first explicit example [80] showing striped phases for fundamental matter realized with defect flavor D7'-branes in AdS 5 × S 5 spacetime [66]. The construction and analysis of such modulated ground states has led to increasingly demanding numerics [104][105][106] due to the highly non-linear DBI action. Another important merit of holographic techniques is to yield useful toy models that keep computations simple. 4 The D3-D5 geometry [29] has provided us with a powerful framework to study anisotropic media; in the present paper, we probe this medium with fundamental matter and explore both the thermodynamics and the excitation spectra along and across the anisotropic direction. We believe that our work sets the standard and sparks many other investigations in the future. One particularly interesting case is in the context of the astrophysics of compact objects where anisotropy seems to offer clues regarding the "universal relations" for neutron stars [111] and the already established holographic modeling of isotropic matter [70,[112][113][114][115][116][117].
While the D3-D5 geometry enables a multitude of investigations of anisotropic matter, it is not at all perfect. The real life anisotropic systems in laboratories consist of a finite number of layers of perhaps different materials with finite separation. The D3-D5 system is an idealization: The anisotropic direction spans an infinite range and the separation between the layers is strictly vanishing. Interestingly, contrary to common beliefs, relaxing these approximations in a controlled fashion may not be a major hurdle and can be obtained by a judicious selection of the smearing form for the D5-branes. Reaching this milestone would be very rewarding and will teach us further lessons on surface phenomena. In particular, it would allow us to dissect the issues of the surface plasmon left open in this paper. Besides, the shear-viscosity over entropy density bound [118,119] is known to be affected by anisotropy [120][121][122][123][124][125]; it would be interesting to consider the status of the bound in the presence of a further scale introduced by a finite inter-layer spacing. Furthermore, a finite-spacing would break the translations along the orthogonal direction to the layers.

A Kappa symmetry
Let us verify that the embeddings we found in Section 3.1 are supersymmetric. We will work in the following vielbein basis for the metric presented in Section 2 for B = 1: where the warp factor h is the function written in (2.4) and g and f are given by: We will write the SU(2) left-invariant one-forms ω 1 , ω 2 , and ω 3 in terms of three angles (θ, ϕ, ψ) as follows: ω 1 = cos ψ dθ + sin ψ sin θ dϕ ω 2 = sin ψ dθ − cos ψ sin θ dϕ where 0 ≤ θ ≤ π, 0 ≤ ϕ < 2π, and 0 ≤ ψ < 4π. The supersymmetric embeddings of the D7-brane are those that satisfy the kappa symmetry condition: where is a Killing spinor of the background. For a D7-brane without any worldvolume gauge field, the matrix Γ κ is given by: where γ α 1 ···α 8 is the antisymmetrized product of induced Dirac matrices and σ 2 is the Pauli matrix. We will take the following set of worldvolume coordinates: and the embedding will be determined by two embedding equations of the type: The kappa symmetry matrix for this ansatz is: To obtain the induced γ-matrices, let us write the vielbein one-forms (A.1) in a coordinate basis as: Then, the γ's are: where the Γ's are flat constant matrices for the vielbein (A.1). Specifically, the induced γ-matrices on the worldvolume for our configuration are: cos χ e g cos ψΓ 1 + sin ψΓ 2 , γ ψ = h 1 4 2 cos χ sin χ e g Γ 3 + cos χ e f Γ 5 , 2 cos χ sin ψ sin θe g Γ 1 − cos ψ sin θe g Γ 2 + sin χ cos θe g Γ 3 + cos χ cos θe f Γ 5 . (A.11) Clearly, we can factorize the antisymmetrized product in (A.8) as: where the first factor is: The second factor can be written as: where the c i coefficients are given by: Putting everything together, we can write Γ κ as: (A. 16) It was proven in [29] that the Killing spinor of the background can be written as: where η is a doublet of constant Majorana-Weyl spinor satisfying the following projection conditions: As Γ 12 commutes with the products of matrices of the rhs of (A.16), we can rewrite the kappa symmetry condition (A.4) as: Moreover, from the last equation in (A.18) one can easily demonstrate that: Thus, Γ κ acting on η can be written as: Moreover, since from (A.18) we have: we can write Γ κ η simply as: According to (A.19), we are looking for embeddings such that Γ k acts on η as plus/minus the identity. Inspecting the rhs of (A.23) we notice that the last two terms act on η as a non-trivial matrix. Therefore, we should require that the coefficient of these terms is zero, namely: which is the same as (3.7). Notice that the general solution of (A.26) is indeed (3.1). To finish with the proof of kappa symmetry, let us compute the terms of Γ κ η which are proportional to the unit matrix. One can readily prove from (2.16) and (A.15) that, when the BPS condition (A.26) holds, we have: From these two equations we get: and, as a consequence,

B The dictionary
The bare quark mass m q is obtained from the Nambu-Goto action of a fundamental string hanging from the boundary to the horizon. We will obtain m q as the limit of the constituent quark mass M c , which contains the effects of the thermal screening of the quarks and is also obtained from the Nambu-Goto action. Let us consider a fundamental string extended in t, P at constant ρ = 0. The induced metric is: whose determinant is: The constituent quark mass M c is minus the action per unit of time of the Nambu-Goto action: where the e φ 2 factor is due the fact that our metric is written in Einstein frame (g string = e φ 2 g Einstein ). When ρ = 0, we have: where J(P 0 ) is defined as the following integral: This integral can be obtained in analytic form and is given by: In order to get the bare quark mass m q , let us write P 0 = m and consider the limit m → ∞, where M c equals m q . We get: The free energy F can be written as: where N is written in (3.29), I bulk is the integral (3.30) evaluated on-shell and I bay is given in (3.40). Since the normalization factor N does not depend on the mass parameter m, we have: In order to compute its derivative with respect to the mass parameter m, let us write the bulk action as: where J is given by: The momentum appearing on the right-hand side of this equation is: Notice also that the terms containing ∂c/∂m vanish and the VEV O m is proportional to c, as expected:

C Critical embeddings
In this appendix we perform an analysis of critical embeddings, at zero density and non-zero temperature, of the probe D7-brane along the lines of [56,65]. In these configurations the D7brane touches the horizon near χ = π/2. It is then useful to define the following new coordinates (y, z) as to zoom on the vicinity of the contact point. In (C.1) C and α are constants to be determined in order to cast the induced metric into a Rindler form. In the new coordinates (C.1), the emblackening factor B(r) near the horizon reads where T is the temperature defined in (6.111), the radial and temporal pieces of the induced metric take the desired Rindler form, namely In terms of the z and y coordinates, the metric in the vicinity of the horizon, at leading order, becomes where the ellipsis represents terms that do not contribute to the embedding of the probe. Describing the D7 embedding with the parametrization y = y(z), one has dz 2 + dy 2 = (1 +ẏ 2 )dz 2 , where the dot denotes derivation with respect to z. Substituting into the metric (C.5), one obtains the associated DBI Lagrangian density The equation of motion descending from (C.6) is The Lagrangian density (C.6) and the associated equation of motion (C.7) correspond to the case n = 3 in the generic analysis performed in [65], according to the fact that three internal directions wrapped by the D7-brane collapse in the critical embedding at y = 0. More generically, the parametrization y = y(z) is suitable to describe near-to-critical black hole embeddings where the D7-brane intersects the horizon at an angle which can be seen solving the equation order by order.
To study the near-to-critical Minkowski embeddings instead, one has to parameterize the embeddings as z = z(y), where z h = z(y = 0) represents the IR radial distance of the D7-brane from the horizon. From (C.5) one gets the DBI Lagrangian where the prime denotes differentiation with respect to y. The associated equation of motion is to be solved imposing the boundary conditions that implies z (y = 0) = 0 . (C.14) The critical solution y = √ 3z , 1 + 1 − r h r 10 3 .

(C.17)
Recalling the definition of the coordinate z in (C.1) and using (C.3), near the horizon one has Adopting the Cartesian-like coordinates (3.55) and recalling (C.1), in the near-horizon region one has On the critical solution (C.15) the coordinate P becomes which corresponds to the following in-falling angle in the (P, ρ) plane: Let us next study the near-critical black hole solutions, which are described by the following deformation of the critical solution (C.15): where ξ(z) is a small function of z. Linearizing (C.7) in the deformation ξ one gets Considering a power-law solution ξ(z) ∼ z ν , we get that ν must satisfy the quadratic equation ν 2 + 3ν + 4 = 0, which has the following two solutions: The general solution to (C.23) is then given by where A and B are numerical coefficients while the power of T has been introduced to comply with the correct physical dimensions. Notice that the perturbed solution (C.25) does not respect the boundary conditions (C.9) discussed for the near-critical differential problem. The perturbation ξ(z) of the critical solution introduced in (C.22) drives us away from criticality: the critical system is dynamically unstable. The linearized problem (C.23) features a scale invariance, namely if y(z) = f (z) for a certain function f is a solution, thenȳ where µ is a number, is a solution too. The rescaled solutionȳ corresponds to the rescaled boundary conditionsȳ Following Appendix F of [56], one can show that the coefficientsĀ andB of the rescaled solution (C.26) are related to the coefficients A and B in (C.25) as follows: where the matrix M(µ) is defined as The matrix (C.30) satisfies the property where one needs to recall (C.28). Therefore one has where v is a constant vector characterizing the whole family of near-critical black hole embeddings.
Since the matrix M(x) is periodic when x is shifted by

D Thermal screening
We wish to study the quark-antiquark potential at non-zero temperature. We consider a string hanging from the boundary and reaching a minimal radial coordinate r 0 , with its endpoints a distance from each other in the isotropic x-direction. We parametrize the string with the radial coordinate, which means that the midpoint will be found when 1 x (r) = 0. The string metric is given by The potential will be given by (minus) the action .

(D.4)
For small temperatures, be approximated as where J is the following integral: Eq. (D.5) can in turn be reverted for r 0 , Moreover, the energy of the string can be written as: whereJ is given by: Finally, inserting (D.7) and removing the zero-point thermal energy introduced by the regularization, Once again we note that the first temperature correction makes the force less attractive. Similarly, we can compute the energy with a separation in the anisotropic z-direction. This case the metric reads with a separation where we used that, at first order in the correction, P 0 ≈ m and that m is related to m q as in (B.8).

E High temperature black hole embeddings
Let us consider now black hole embeddings for d = 0, which we will describe by considering the angle χ as a function of the radial variable r. The equation for χ(r) written in (2.20) reduces to: We will analyze the high temperature limit in which χ is small for all values of the radial variable r. At first order in χ the equation of motion becomes: The horizon radius in (E.2) can be scaled out by redefining the radial variable and χ, this is equivalent to setting r h = 1. The equation for χ becomes: The general solution of this equation can be written in terms of hypergeometric functions as: , where c 1 and c 2 are integration constants. The function on the right-hand side of (E.4) is divergent at the horizon r = 1 for generic values of c 1 and c 2 . Indeed, one has: and, thus, near r = 1 the angle χ behaves as: . (E.7) Let us next look at the UV behavior r → ∞ of χ. In general, for large z the hypergeometric functions behave as Thus, at the UV: χ ≈ m r r 8 9 + c r r 28 9 , (E. 10) where m r and c r are the mass and condensate parameters in the r variable. Let us prove that, when c 1 and c 2 fulfill (E.7), m r and c r are real. Using (E.9), we get that the coefficient of the leading term of χ as r → ∞ is: . (E.13) Using (E.13) one can show that the imaginary part of (E.11) vanishes and that the real part is given by: 14) The coefficient of the subleading term is: we can easily show that the imaginary part of (E.15) vanishes and the real part is given by: Eliminating c 2 between (E.14) and (E.17), we get the following relation between the condensate parameter c r and the mass parameter m r :

E.1 On-shell action
Let us now obtain an approximate expression for the free energy F in the limit of high temperature (or small mass parameter m). In this regime our embedding is a black hole embedding with a small value of η(u) for all values of the holographic coordinate u. In a first, zero-order, approximation we can just take η = 0 in the bulk and boundary actions. We get: Then, the zero-order free energy F (0) is given by: Therefore, we can approximate F at leading order in T as: Notice that this means that F ∝ T 4 and, therefore it obeys a Stefan-Boltzmann law at leading order in T . To find the corrections to this law, let us consider the terms in the action that are quadratic in η. For the bulk action these terms are: To evaluate this integral when η(u) is a solution to the equations of motion at quadratic order, we apply the method used in appendix C.1 of [56]. Suppose that we have an action S which depends quadratically on η(u) andη(u) as: where F 1 and F 2 are known functions of u. Then, the action S evaluated on a solution to the equation of motion is: .
(E. 30) In our case the function F 1 (u) is: We found in (E.22) that c ∼ m in this high T regime. Therefore, the last term in (E.36) is proportional to m 2 . As m ∼ T −b for fixed quark mass m q , it follows that the second term on the right-hand side of (E.36) gives rise to a subleading contribution that corrects the dominant Stefan-Boltzmann law (F ∝ T 4 ) at large T . This subdominant contribution grows as T 4−2b = T 20 9 .

F Fluctuations
In this appendix we obtain the equations of motion for the fluctuations around the χ = 0 massless embedding of the D7-brane at non-zero temperature and density. Accordingly, let us allow the worldvolume gauge field A to fluctuate as: where A (0) is the one-form written in (4.1). Let us follow the methodology of [72] and split the total field strength as F = F (0) + f and the DBI matrix as: where X is defined as: Then, we expand the DBI determinant in series as: − det(g 8 + e − φ 2 F ) = − det(g 8 + e − φ 2 F (0) ) 1 + (F.4) Let us split the inverse of the matrix g 8 + F (0) as: where G −1 is the symmetric part and J is the antisymmetric part (G is the so-called open string metric). It follows that: It is easy to prove that the first-order term in X is a total derivative and, therefore, it does not contribute to the equations of motion of the gauge field and can be neglected. Up to second order, the Lagrangian density for the fluctuations takes the form: where H is the function written in (4.9). If we now define the prefactor L * as: then L is given by: L ∼ L * G ac G bd − J ac J bd + 1 2 J cd J ab f cd f ab . (F.10) The corresponding equation of motion for gauge field component a d is: Before going further, let us write the non-vanishing components of the open string metric for our system: The non-vanishing elements of the antisymmetric tensor are: Let us write explicitly the equations for the fluctuations. We choose the gauge in which: a r = 0 . (F.14) The equation of motion of a r in this gauge is the following first-order constraint: G tt ∂ t f tr + G xx ∂ x f xr + G xx ∂ y f yr + G zz ∂ z f zr = 0 . (F.15) Moreover, the equation for a t is: 16) and the equations for the components of a along x ≡ x 1 and y ≡ x 2 are: Finally, it remains to write the equation for the gauge field along the anisotropic direction z = x 3 , which is: ∂ r L * G zz G rr f rz + L * G zz G tt ∂ t f tz + G xx ∂ x f xz + G xx ∂ y f yz = 0 . (F.18)

F.1 In-plane propagation
Let us consider a wave propagating in the x direction. This means that all the gauge field components a ν only depend on r, t, and x. Let us Fourier transform the gauge field to momentum space as: a ν (r, t, x) = dω dk (2π) 2 a ν (r, ω, k) e −iω t + ikx . (F. 19) From now on in this section all the equations are written in momentum space. The first-order constraint (F.15) takes the form: G tt ω a t − G xx k a x = 0 , (F.20) where the prime denotes derivative with respect to r. We now define the electric field E as the gauge-invariant combination: E = k a t + ω a x .

(F.21)
This gauge-invariant combination E appears in the equation for a t and a x , which are given by: Actually, we can use (F.20) and the radial derivative of the definition of E to write a t and a x in terms of E : a t = G xx k G xx k 2 + G tt ω 2 E , a x = G tt ω G xx k 2 + G tt ω 2 E .
(F. 23) Plugging this result in the two equations in (F. 22), we arrive at a unique second-order equation for E: Moreover, the equations for a y and a z are decoupled and given by: ∂ r L * G xx G rr a y − L * G xx G xx k 2 + G tt ω 2 a y = 0 ∂ r L * G zz G rr a z − L * G zz G xx k 2 + G tt ω 2 a z = 0 .