Screening of Coulomb interactions in Holography

We introduce Coulomb interactions in the holographic description of strongly interacting systems, by performing a (current-current) double-trace deformation of the boundary theory. In the theory dual to a Reissner-Nordstr\"om background, this deformation leads to gapped plasmon modes in the density-density response, as expected from conventional RPA calculations. We further show that by introducing a $(d + 1)$-dimensional Coulomb interaction in a boundary theory in $d$ spacetime dimensions, we recover plasmon modes whose dispersion is proportional to $\sqrt{|\mathbf{k}|}$, as observed for example in graphene layers. Moreover, motivated by recent experimental results in layered cuprate high-temperature superconductors, we present a toy model for a layered system consisting of an infinite stack of (spatially) two-dimensional layers, that are coupled only by the long-range Coulomb interaction. This leads to low-energy `acoustic plasmons'. Finally, we compute the optical conductivity of the deformed theory in $d = 3 + 1$, where a logarithmic correction is present and we show how this can be related to the conductivity measured in Dirac and Weyl semimetals.


Introduction
The AdS/CFT correspondence [1][2][3] has become an important tool for studying strongly coupled quantum field theories, and in the last decade it found an ever increasing range of applications into the realm of condensed-matter physics [4][5][6][7]. Many recent experiments demonstrate that in strongly interacting systems as, for example, the cuprates hightemperature superconductors, the observed behavior can deviate quite drastically from well known condensed-matter results based on weakly coupled theories that admit a quasiparticles description. In particular, these materials exhibit a strange-metal phase for temperatures higher than the superconducting critical temperature, where quasiparticle excitations appear to be absent (see e.g. ref. [7]). The promising aspect of the gauge/gravity duality is that it allows us to study systems without quasiparticle excitations. It can be used to compute thermodynamics quantities and response functions with relative little effort, by mapping strongly interacting quantum field theories on the boundary of an asymptotically AdS spacetime, to classical gravity in the curved bulk spacetime. Trying to tune holographic models such that they reproduce the properties of laboratory condensed-matter systems as closely as possible is therefore one important goal of ongoing research in holographic applications to condensed matter. A common feature in conventional holographic calculations is that they describe neutral systems, where the low-energy hydrodynamic excitations in the longitudinal channel contain sound modes. However, it is well established that in metals the density fluctuations of the charged electrons give rise to a different type of collective excitations known as plasmons [8]. In order to have more realistic models of strongly interacting metals and superconductors, it is thus important to modify the holographic theory to include the effect of the long-range Coulomb interaction that turns the linear sound modes into plasmon modes.
In textbook condensed-matter response calculations the electron-electron interaction is introduced in the description through a self-consisted field method known as the randomphase approximation (RPA) [9][10][11][12], where the induced charge without Coulomb interaction is replaced with the screened charge. In a translational invariant system the density-density response function χ(ω, k) then becomes, in terms of the non-interacting response function χ 0 (ω, k), equal to where (ω, k) is the dynamical dielectric function, and V (ω, k) the Fourier transform of the Coulomb potential. In non-relativistic systems, where the Fermi velocity is much smaller than the speed of light, the latter is only a function of the magnitude of the momentum |k|. In particular in 2 + 1 and 3 + 1 spacetime dimensions we have that V (ω, k) = e 2 2 0 |k| , 2 + 1 dimensions V (ω, k) = e 2 0 k 2 , 3 + 1 dimensions, (1.2) with −e the charge of the electron, and 0 the dielectric constant.
In this paper, we develope a generic procedure to include the above RPA correction, and hence the effect of Coulomb interactions, into holographic models at nonzero density. A first step towards studying plasmon modes in holography was recently made in refs. [13,14], where it was shown that plasmon quasinormal-modes in the Reissner-Nordstrom model can be studied by imposing a new type of boundary conditions on the bulk linearized equations of motion. Here we explain how these boundary conditions naturally arise by introducing a (current-current) double-trace deformation of the boundary theory, and how this gives rise to a RPA-like form of the response function. Using also the example of the Reissner-Nordström metal, we numerically study the full frequency and momentum dependence of the spectral functions and show how the prescribed procedure turns the low-energy sound modes into gapped plasmon modes.
In the particular case of a (3 + 1)-dimensional theory, we also show that by introducing long-range Coulomb interactions in the holographic theory, we obtain an anomalous cutoff-dependent logarithmic behavior in the real part of the optical conductivity, analogous to the one observed in the conductivity of Dirac and Weyl semimetals [15].
In (2+1)-dimensions, we adapt our procedure to correctly describe a system where electrons are constrained to two spatial dimensions, but where the Coulomb potential exists in three spatial dimensions. This allows us to obtain a low-energy dispersion relation ω ∝ |k| as expected in charged (spatially) two-dimensional system as, for example, graphene [16,17]. Furthermore, since experiments on (2 + 1)-dimensional materials are often performed on a multi-layer system, we propose a toy model describing an infinite stack of (2+1)-dimensional layers, where the dynamics in each layer is determined holographically through a dual Reissner-Nordström theory and where the coupling between the layers is given only by the long-range Coulomb interaction. A recent example of the interest in these multi-layered systems is given in ref. [18], where the spectral function for a layered copper-oxide hightemperature superconductor has been measured. These experimental results suggests that it is indeed the Coulomb interaction that governs the coupling between different layers, validating the assumptions of our toy model, and hinting at the fact that low-energy charge fluctuations might be very relevant for the description of the dynamics of high-T c superconductors. We therefore show how our toy model qualitatively reproduces the observed spectral functions, with low-energy linear excitations, often referred to as 'acoustic plasmons', believed to possibly play an important role in the mechanism of high-temperature superconductivity [19][20][21].
The structure of the paper is as follow. In section 2 we briefly summarize the Reissner-Nordström model and its spectral functions containing the sound modes. In section 3 we explain how to deform the theory to describe a charged system with plasmon excitations and show how the density-density response function obtains the RPA-like form. In section 4 we present the numerical results for the spectral functions with plasmon excitations in a (2+1)-dimensional theory and propose a toy model for a layered system in 3+1 dimensions. Finally, section 5 contains the results for the optical conductivity for a (3 + 1)-dimensional system with Coulomb interactions.

Reissner-Nordström black brane
In this section we introduce the Reissner-Nordström model that we use throughout the paper as a concrete example of the proposed procedure, that we introduce in detail in section 3. The dynamics of this holographic model has been thoroughly studied in refs. [22][23][24]. In order to be able to compare directly with our later results, however, we briefly review its dynamical properties by computing the spectral functions. We then recover that the low-energy excitations of the theory in the longitudinal sector are linear sound modes and a diffusive mode.
Throughout the paper, we call d the spacetime dimension of the boundary theory, dual to a gravitational (d + 1)-dimensional bulk theory and we use the mostly-plus metric. We denote with r the coordinate of the extra spatial dimension of the bulk theory, with the boundary at r → ∞ unless explicitly stated otherwise. When Fourier transforming in the direction orthogonal to r, we denote with k = (ω/c, k) the d-dimensional momentum. We use, with a slight abuse of notation, greek lower-case tensor indices for both the (d + 1)dimensional bulk and the d-dimensional boundary, as the range of the indices is clear from the context. Moreover, in computing dynamical quantities, we always use a gauge in which all the r-components of the field fluctuations are zero, so that no confusion should arise.

Gravity action
The simplest holographic model with a nonzero-density boundary theory is the Reissner-Nordström model, described by a bulk Einstein-Maxwell action, that in SI units reads , i.e., the bulk field A t has the dimensions of a chemical potential. Since we want to describe a uniform boundary theory, i.e., the chemical potential and the energy-momentum tensor do not depend on the spacetime coordinates, we look for a solution of the form A µ = (0, A 0 (r), 0) and g µν = g µν (r).
For notational convenience, we first rewrite the action in terms of dimensionless fields and coordinates, by redefining the fields as with L 2 = −d(d − 1)/2Λ the AdS radius squared and r 0 denotes the position of the blackbrane (outer) horizon. Notice that this rescaling thus fixes the horizon atr 0 = 1. However, using the scaling symmetry of the actioñ we can obtain any solution from the solution withr 0 = 1. From now on we only use dimensionless variables and omit the tildes. In the end, the theory we are going to study, including boundary terms, is described by the action: , (2.4) where now every field and coordinate is dimensionless, as is the prefactor in front of the integral that for now we take to be equal to one, by an appropriate scaling of the action, but we briefly come back to this prefactor later on. The second term is the Gibbons-Hawking-York boundary term, with γ µν ≡ g µν − n µ n ν the induced metric on the boundary, n µ the unit vector normal to the boundary, and the determinant γ of γ µν is taken only in the directions orthogonal to n µ . Finally, K is the trace of the second fundamental form K = γ µν ∇ µ n ν . This boundary term is necessary to have a well-defined Dirichlet problem on the boundary [25][26][27]. In the case at hand of an asymptotically flat spacetime we have n ν = √ g rr δ r µ and K reduces to K = − √ g rr γ µν Γ r µν , with Γ r µν the bulk Christoffel symbols. The last term in the action contains the counterterms necessary to regulate the divergences in the theory, and its explicit form depends on the dimension d of the boundary spacetime. We present later the counterterms in the case of d = 2 + 1 and d = 3 + 1. For a detailed treatment see, for instance, ref. [28].
The classical equations of motion derived from the above action are solved by the Reissner-Nordström background, describing an asymptotically AdS charged black brane, that in the units we introduced takes the form for d > 2, that, according to the holographic dictionary we interpret as the chemical potential of the boundary theory. Here and below, the prime denotes derivatives with respect to the radial coordinate. From the solution of f (r), we can see that the mass M and charge Q of the black brane are expressed in terms of µ as The solution is then completely characterized by one dimensionless parameter, that we take to be T /µ, with the black-brane temperature given by The equilibrium properties of the theory are well known (see e.g. ref. [6]) and can be derived from the GKPW rule after regularizing the action at the cutoff scale r = r uv by inserting the boundary counterterm We give here for later convenience the nonzero components of the equilibrium expectation values in the chosen units, that for d > 2, are: Here S cl is a reminder that the action is evaluated at the classical solution and the superscript in δg µ means that we are deriving with respect to the leading-order term in the asymptotic expansion of the field for large r. We see here that, in order to compare holographic results to experiments, we should choose, in this bottom-up approach, the prefactor in front of the action in eq. (2.4), that we conveniently scaled away, such to tune the thermodynamics quantities to values close to the experimental ones.

Fluctuations and response functions
The background solution of the gravity action in eq. (2.4) defines the equilibrium properties of the boundary theory. In order to study the response of the field theory to small perturbations, we have to consider fluctuations of the classical fields on top of the background solution: (2.13) Exploiting rotational invariance to fix the momentum along the x direction, so that k = (ω, ±|k|, 0, . . . , 0), the fluctuations decouple according to their parity under O(d − 2) acting on x 2 , . . . , x d−1 [29]. In particular, we are interested in the longitudinal channel associated with spin 0, as it is the one containing the density-density response function that in d = 3+1 and d = 2 + 1 are described by the components δΦ ≡ (a t , a x , a r , h xt , h tt , h xx , h yy = h zz , h rr , h tr , h xr ) , (2.14) and will ultimately lead to hydrodynamic diffusion and sound modes. It is convenient to work in a gauge where h rµ = 0, as well as A r = a r = 0. The set of coupled fluctuations then becomes Solving the coupled set of linearized equations of motion for these fluctuations, with infalling boundary conditions at the black-brane horizon, allows us, after eventual renormalization of boundary divergences (see ref. [30] for a detailed treatment), to compute the retarded Green's functions of the dual boundary theory. To do so, we follow the procedure in ref. [31]. A caveat is that to extract the matrix of Green's functions for a set of M coupled fields we need M independent solutions. In the case at hand, however, our choice to work in a gauge where all the radial components are zero does not completely fix the gauge freedom of a µ and h µν . This implies that we can only find two independent solutions of the linearized equations of motion, corresponding to the two physical degrees of freedom. The remainig solutions needed to compute the full Green's function matrix are pure gauge solutions that can be constructed analytically. In appendix A, we show how to construct these gauge solutions and we provide their explicit form for the theory at hand. From the retarded Green's function G R (ω, k) we can obtain the spectral function defined as The spectral function is an interesting quantity as it tells us the rate of work done on the system by a small external perturbation at a given frequency [8]. Sharp peaks in the spectral function, that correspond to poles in the Green's function, denote long-lived collective modes, with a lifetime determined by the width of the peak.
In figures 1 and 2 we present the spectral functions of a (2 + 1)-dimensional system for the density response A ρρ as well as for the energy-density response A for two different values of the temperature. Here, as in the rest of the paper, we work in the grand-canonical ensemble, hence all the plots are rescaled to a solution with unit chemical potential. Since we are describing a strongly interacting system, we expect from hydrodynamics considerations (see e.g. ref. [32]) to find two linear sound modes due to energy-momentum conservation. As the poles are shared by all Green's functions, although the residues differ, in both spectral functions the low-energy linear sound modes ω = ±v s |k| are evident. In particular, up to the numerical precision obtained, the speed of sound satisfies the relation v s = 1/ √ d − 1 in the low-energy regime [24,33]. Moreover, we can see that the height of the peak of the sound modes in the energy-density spectral function increases as the ratio T /µ is raised. In the density spectral function we observe, in addition to the sound modes, the expected diffusive mode due to charge conservation that corresponds to the broad quadratic band.  Figure 1: The density spectral function for two different values of the temperature in 2+1 dimensions. Notice that in addition to the linear sound modes with a low-energy behavior of the dispersion given by ω = ±|k|/ √ 2, we can clearly observe the diffusive mode, that becomes more dominant as we increase T /µ.
Here the linear modes become less sharply peaked when the temperature is raised as part of the spectral weight shifts into the diffusive mode. Note that the diffusive mode is absent in A , because in that mode the charge fluctuations do not induce a velocity field v(t, x) in the fluid, hence also no energy transport. A thorough analysis of these modes has been presented for (2 + 1)-dimensional systems in ref. [24]. In the next section, we turn to the problem of deforming the holographic theory in order to describe a charged quantum field theory with plasmon excitations.

Plasmon modes in holographic theories
In this section we explain how to introduce Coulomb interactions in the boundary quantum field theory, and show how this gives rise to an RPA-like form of the response functions, that turns the linear modes observed in the Reissner-Nordström spectral function into long-lived plasmon excitations.
The holographic theory presented in the previous section is dual to a charge-neutral boundary theory, that is, the Maxwell field in the bulk acts only as a source for the current operator J µ in the boundary theory, and there is no dynamical photon [34]. However, in condensed matter we are often interested in describing charged systems, i.e., systems with a long-range Coulomb interaction. We, therefore, need to modify the boundary theory to add dynamical photons coupled to the conserved current J µ , which are described by the The energy-density spectral function for two different values of the temperature in 2 + 1 dimensions. The sound modes are clearly visible and we can notice that the associated peaks becomes higher for a higher value of the temperature.

Maxwell action
with µ 0 the magnetic permeability and −e the electron charge. Note that we also used the fact that our boundary theory is defined on flat Minkowski spacetime, so that √ −η = 1. Integrating out the A µ field from eq. (3.1), using the standard ξ procedure to account for the gauge freedom, we obtain in momentum space (see for example ref. [35]) Notice that eq. (3.2) assumes the form of a double-trace deformation of the boundary theory for the current operator.
As first described in refs. [36,37], a double-trace deformation in the large-N limit can be incorporated in the holographic dual simply as a change in the boundary conditions for the corresponding fields. In the gauge/gravity duality without multi-trace deformation, we can obtain the response to small fluctuations from the asymptotic behavior of the field fluctuations. For example in the case of the Maxwell-field fluctuations we have and according to the holographic dictionary, we interpret the leading-order coefficient a (0) µ as the source a s µ of the fluctuations δ J µ . Following the prescription in ref. [36], in order to include a multi-trace deformation of the response we still interpret the subleading term in equation (3.3) as the response to the source fluctuations, but we impose the boundary condition The reason for this boundary condition can be easily understood by looking at the boundary behavior of the holographic theory. The boundary contribution from the Maxwell-field fluctuations in the renormalized holographic action, eq. (2.4), evaluated on-shell is given, at second order in fluctuations, by By introducing a boundary contribution of the form of eq. (3.1) at infinity, that in the saddle point approximation of the large-N limit, we write in dimensionless units as with α 2 a dimensionless constant that depends on the charge of the system, its permittivity and on the prefactor of the action in eq. (2.4), the boundary term for the fluctuations now becomes Since the double-trace deformation does not modify the bulk theory, we see that this is equivalent to the Reissner-Nordström theory described above, where we simply modify the boundary conditions that are of the form of eq. (3.4). Notice that these boundary conditions are equivalent to the ones used in ref. [14] to study plasmonic quasi-normal modes in a d = 2 + 1 holographic theory, but that were obtained differently by looking for zeros of the dielectric function.
The effect of the deformation in eq. (3.7) on the response functions is easily derived, and in practice there is no need to solve the linearized equations of motion by shooting for the boundary condition in eq. (3.9). In fact, ignoring for the moment the gravitational part, as it does not influence the procedure for the density-density and current-current spectral functions we are interested in, we have from linear-response theory without the double-trace deformation Since the bulk theory is left unchanged by the deformation, the Maxwell-field fluctuations a µ satisfies the same linearized equations of motion with the same near-boundary behavior, the only difference is that now we interpret the leading-order coefficient a ν not as the source alone but according to the right-hand side of eq. (3.9). In terms of the source a s µ in the presence of a double-trace deformation, expression (3.10) now becomes Rearranging the equation in order to put it into the form δ J µ = χ µν a s ν to extract the response function χ µν we have that and we see that the response function assumes the RPA-like form Explicitly, introducing the polarization Π by means of we obtain that is of course independent of the gauge-choice parameter ξ. In particular, we can see that the density-density response function is , (3.16) as expected from RPA calculations.
From relativistic hydrodynamics we know that in the low-energy limit G 00 assumes the form [32] G 00 (ω, k) with the equilibrium quantities defined in eqs. (2.10)-(2.12), and v 2 s the speed of sound . This suggest that the RPA-like response function in eq. (3.16) contains a gapped plasmon mode, with the plasma frequency determined by In figure 3, we show side by side the density-density spectral function with and without double-trace deformation, for two different values of the parameter T /µ, where we see that the linear sound modes turn into a gapped plasmon mode. The dashed black line represent the hydrodynamic approximation of the plasmon mode from eq. (3.17). We can think of higher values of T /µ as moving towards the zero chemical potential limit, where there is no plasmon mode. By comparing the two figures we can indeed observe that for higher value of T /µ the low-energy gapped mode becomes less visible in the spectral function. This can be seen in more details by comparing with figure 4, where we show some slices of the spectral function for different fixed values of momentum |k|/µ. Here we also notice that the height of the peak of the plasmon modes is lower than the corresponding peak in the Reissner-Nordström solution. This is due to the screening effect of the charged particles that opposes density fluctuations. Moreover, we see that the plasmon peak initially increases as we move to smaller values of |k|/µ, as was the case for the linear modes, but it then starts to decrease for small value of |k|/µ, due to the k 2 factor in the low-energy limit of the densitydensity spectral function. Even though in this paper we mainly focus on the density-density response, the plasmon modes are of course also present in the current-current spectral function as we show in figure 5. Here the plasmon mode is more easily visible as there is no k 2 suppression, since from Ward's identities we know that This gapped plasmon mode is observed in (3 + 1)-dimensional metals, however in an experimental (2 + 1)-dimensional system we expect to find plasmon modes of the form ω ∝ |k|.
Notice that the expression for the RPA-like correction to the response function in eq. (3.16) has been derived in an arbitrary spacetime dimension and it is independent of the value of d. This implies that in any dimension, the qualitative behavior of the spectral function looks the same as the one in figure 3, although the position and residue of the poles of course change. In the next section we show how to modify the double-trace deformation to, from an experimental point of view, correctly describe plasmon modes in d = 2 + 1. The density spectral function for T 0.02µ (top) and T 1.18µ (bottom) without (left) and with (right) double-trace deformation. Here we have taken α 2 = 1/25. The black dashed line represent the plasmon modes computed in the hydrodynamic approximation. We see that the double-trace deformation turns the low-energy sound modes into a gapped plasmon mode, although the low-energy excitation are less visible as the height of the peak for low |k|/µ becomes smaller for higher values of T /µ.

Plasmon modes in d = 2 + 1
As already noted in ref. [14], the presence of a gapped mode is due to the fact that in the holographic theory of a (2+1)-dimensional model, we also constrain the boundary dynamical photons to live in 2 + 1 dimensions, while in layered materials studied in laboratories, the current is constrained on the (2 + 1)-dimensional layer, but the photons are free to move in all the spacetime dimensions, giving rise to an effective Coulomb potential V (k) ∝ 1/|k|.  . We see that the height of the peak of the plasmon modes are lower than the corresponding sound modes due to screening of the charged particles. Moreover we observe that for low momenta, the peak in the RPA response function decreases as we move towards smaller values of |k|/µ.
Put differently, the Coulomb potential between charges at positions x and x behaves as 1/|x − x | and not as log|x − x |.

Single Layer
In order to model a more realistic (2 + 1)-dimensional system, where photons that mediate the Coulomb interaction are free to move in the whole three-dimensional space, we consider a d = 2 + 1 boundary theory, dual to the (d + 1)-dimensional bulk theory from eq. (2.4), but we define the deformation of the boundary theory by starting from a (d + 1)-dimensional boundary term. Defining the (d + 1)-dimensional vector (x, z) and momentum (k, k z ), with x and k living in d spacetime dimensions while z and k z represent the extra spatial dimension normal to the layer and the associated momentum respectively, we restrict the current operator to a plane as (4.1) Inserting this expression into eq. (3.7), Fourier transforming and performing the Gaussian integral over the Maxwell field we obtain that can be integrated over k z to obtain the d-dimensional boundary deformation This in particular implies that the response function now takes the form: However, as we are ultimately interested in describing condensed-matter systems, where the Fermi velocity is considerably smaller than the speed of light, so v F c and we look into a regime where |k| ω in eq. (4.4). As a result, the density-density response function is well approximated by the familiar form , (4.5) that is expressed in terms of the static Coulomb potential expected in two-dimensional metals. Notice that the coupling constant for the single layer is α 2 /2, in accordance with eq. (1.2).
In figure 6 we see, for two different values of T /µ, how this form of the potential does indeed give rise to ω ∝ |k| modes observed in a (spatially) two-dimensional system such as graphene [17]. As before, the black dashed lines show the hydrodynamic plasmon modes In figure 7, we observe the momentum dependence of the peak for T 0.02µ, where again we notice that the height of the peak in the charged system is less pronounced with respect to the one in the Reissner-Nordström solution, due to the screening effect of charged particles. Furthermore, in figure 8, we study the temperature dependence of this mode for a fixed value of |k|/µ. In the Reissner-Nordström solution, the position of the peak for low enough momenta is independent of the temperature and given by the speed of sound ω = v s |k| = |k|/ √ 2. In contrast, when we introduce Coulomb interactions we can see that, while the peak still gets sharper and the height increases as we lower the temperature, the position of the peak is temperature dependent. For high temperatures, as the temperature is raised the mode shifts towards the position of the peak of the neutral solution (ω/µ 0.035 for |k|/µ = 0.05), since temperature fluctuations start to dominates over the effect of Coulomb interactions. As we decrease the temperature, however, while the peak initially shifts to higher frequencies, it reaches a maximum frequency for T /µ = 1/2π (purple line in the plot) before starting to move back on the frequency axis as we further lower the temperature. This behavior is in accordance with the hydrodynamic prediction for the plasma frequency in eq. (4.6). In figure 9 we show the temperature dependence of this hydrodynamic prediction in the grand-canonical ensemble, that is given by the temperature dependence of the thermodynamics quantities in eqs. (2.10)-(2.12). The green line includes the higher-order correction expressed in terms of the temperature-independent sound velocity v s = 1/ √ 2 The hydrodynamic plasma frequency shows a maximum at T /µ = 1/2π, given by that, for α 2 /2 = 1/10 and |k| = µ/100, gives ω p 0.0268µ, in agreement with the holo-  . We see that due to screening in the charged system the plasmon peaks shift towards higher frequencies and are lower in height than the corresponding peak of the sound mode. graphic results shown in figure 8. On the other hand, as T /µ → ∞, ρ 2 / → 0 and we recover the sound dispersion.

Layered system
In a similar fashion to what we did in the last section, we present here a toy model for an infinite stack of (2 + 1)-dimensional layers. We show that we recover the form of the Coulomb potential for a layered electron gas [38,39], and we present the density spectral function for a double-trace deformation of this form. The latter shows a dispersion relation for the We see that the peaks become less pronounced as we raise the temperature due to the effect of temperature fluctuations. Moreover, for high temperatures, the position of the plasmon mode shifts towards the position of the sound mode in the neutral system, as temperature fluctuations starts dominating over the Coulomb interaction. As we lower the temperature we see that the plasma frequency reaches a maximum at T /µ = 1/2π (purple), before moving to lower values as the temperature is decreased further. : Temperature dependence of the hydrodynamic plasma frequency for |k|/µ = 0.01 and α 2 /2 = 1/10, with (green) and without (blue) a higher-order momentum correction v 2 s k 2 . We see that the plasma frequency has a maximum at T /µ = 1/2π, corresponding to ω p 0.0268µ.
low-energy plasmon excitations that changes from linear to a gapped mode as a function of Bloch momentum p in the direction perpendicular to the layers. This behavior of the plasmon modes has very recently been observed in high-temperature cuprate superconductors consisting of stacked conducting copper-oxide layers [18]. For our toy model, we make the approximation that there are both strong in-plane shortrange interactions, and a long-range Coulomb interaction that couples the different layers.
We, therefore, model each single layer as an independent d = 2 + 1 boundary system dual to the 3 + 1 Reissner-Nordström bulk theory, but we define a (3 + 1)-dimensional current as with n ∈ Z the layer index. The layers are stacked along the z axis with the distance between each layer, and J µ (x, n ) is the boundary operator computed from holographic calculations with sources a (0) n not necessarily equal to each other.
Using eq. (4.9) into the boundary deformation in eq. (3.7) and choosing the gauge ξ = 1 for convenience we obtain the action with µ = 0, . . . , d − 1, and we have integrated out A z (x, z). Fourier transforming and integrating out the Maxwell field, we obtain We can next perform the integral over the momentum k z and further Fourier transform J µ (k, n ) = 2π π/ −π/ dp J µ (k, p)e ipn , so that we can then sum over the layer indices to obtain the desired double-trace deformation that gives the Coulomb potential for a layered electron gas. Notice that when the layer sources are in phase, i.e., cos(p ) = 1, we recover the potential for a (3 + 1)-dimensional system proportional to 1/k 2 for |k| 1 and such that |k| The spectral function then contains, for p = 0, a gapped plasmon mode, as it is shown in figure 10a, where we plot the density spectral function at T 0.02µ for a Reissner-Nordström solution with a double-trace deformation of the form of eq. (4.12). However, as the only coupling between the layers is due to the Coulomb interaction, the Green's functions of each single layer without the double-trace deformation contains linear modes typical of a (2 + 1)-dimensional system and with a dispersion ω = |k|/ √ 2 for low energies, and the resulting gapped mode is therefore different from the one in a (3 + 1)-dimensional material with the same characteristics. On the other hand for any cos(p ) = 0 we can distinguish three different regimes. For |k| |cos(p )| we obtain a constant potential 2|k| sinh (|k| ) cosh (|k| ) − cos(p ) (4.14) In the Reissner-Nordström system with a linear low-energy dispersion relation, this has the effect of renormalizing the speed of these sound modes. However, for |cos(p )| |k| 1, we recover the 1/k 2 potential. Finally, for |k| 1 we obtain a /2|k| potential. As for the (2 + 1)-dimensional case, in condensed-matter systems we are interested in the limit where ω in the potential satisfies ω |k| as it is suppressed by a factor of v f /c 1, with v f the Fermi velocity. In this limit, the latter potential gives rise to a dispersion relation ω ∝ |k|. This behavior can be seen in figures 10b-10d. In particular, by looking at the hydrodynamic approximation we can see that, by introducing a double-trace deformation as in eq. (4.12), we obtain a response function whose poles are shown in figure 11 for several value of p , from p = 0, that gives the gapped plasmon mode, to p = π that leads to the lowest speed of sound, since in the hydrodynamic limit the renormalized speed of sound is given bỹ In figures 10a-10d we see that the low-energy behavior of the acoustic plasmon modes is well described by the hydrodynamics approximation (shown as a black dashed line). As we move to higher frequencies, however, we observe a discrepancy between the holographic result and the hydrodynamic approximation. This is not unexpected as hydrodynamics is only a long-wavelength theory.  Figure 10: The density spectral function for a layered system at T 0.02µ. When the sources are in phase (a), i.e., cos(p ) = 1, we recover the gapped mode of (3+1)-dimensional materials. However, if the sources are out of phase, for k |cos(p )| the dispersion relation is linear. We can also observe that for |cos(p )| k 1, we recover the quadratic modes (b). For k 1 we instead recover the dispersion relation ω ∝ |k|, as more noticeable in (c)-(d). The black dashed line represent the plasmon modes computed in the hydrodynamic approximation. In all the plots we used α 2 = 1 and = 10.

Plasmons in 3 + 1 dimensions and conductivity
In this section we apply the double-trace deformation to the holographic calculation of the (3 + 1)-dimensional optical conductivity in the Reissner-Nordström metal. Although it would seem straightforward to extend the results of the previous section to a (3 + 1)-  Figure 11: Effect of the double-trace deformation for a layered system on a response function with low-energy sound modes. For p = 0 we obtain a gapped plasmon mode, while for any value of p ∈ (0, π] the low-energy dispersion relation is linear, with p = π having the lowest speed. The plot shows the plasmon dispersion for p = {0, π/50, π/25, nπ/8}, with n ∈ [1,8] an integer. dimensional model, there is an important subtlety as, in all even boundary dimensions, the theory contains a logarithmic divergence that introduces a scale into the theory that needs to be determined experimentally. Without a double-trace deformation this scale only enters the real part of the Green's functions, and it, therefore, does not modify the spectral functions. However, when introducing the double-trace deformation the scale affects the imaginary part of the Green's functions as well, and can be observed in the spectral function, or equivalently, in the real part of the conductivity. We show that this gives rise to a form of the conductivity that resembles the one measured in Weyl and Dirac semimetals.

Renormalization and anomaly
In order to study the properties of the boundary theory we need to regularize the divergences of the boundary action. The necessary counterterms depend on the spacetime dimension of the theory, and for the gravitational part of the action they are treated in detail in ref. [28] up to d = 6.
In the case of d = 2 + 1 treated above, the Maxwell term is finite as r → ∞, and the only UV divergence comes from the gravitational part that we regularized by adding the counterterm in eq. (2.9). More interesting, however, is the case of even boundary spacetime dimension, as the 3 + 1-dimensional system we are interested in. The asymptotic expansion of the fields as r → ∞ is where the expansion coefficients are functions of k, and the logarithmic terms only appear in even dimensions. All the higher-order coefficients as well as the coefficient of the logarithmic term are local functions of the leading-order coefficients a in the present case, that require the full solution of the linearized equations of motion. 1 In d = 3 + 1, we see that in addition to the divergence regulated by the counterterm in (2.9), both the Maxwell and the gravitational field then present a logarithmic divergence, related to the conformal anomaly of the boundary theory [30].
As we show below, the regularization of the logarithm introduce a renormalization scale, k uv , that cannot be determined from the theory but enters as an experimental parameter. Since we are mostly interested in studying the density-density response function, here we focus on the logarithmic term coming from the expansion of the Maxwell field. The case of the metric field is analogous. The boundary term relevant for the J µ J ν correlation functions is Inserting the asymptotic expansion of eq. (5.1) in d = 3 + 1, and the expressions for the induced boundary metric and the vector normal to the boundary, we obtain So, in order to take the limit we need to regularize the logarithmic divergence. This can be done by inserting a scale-dependent counterterm that gives us the regularized boundary term δS (2) A,bdy = where we defined k uv =k uv √ e. We therefore see that the response of the current operator to small fluctuations depends on the choice of scale: The coefficients a ν can be expressed in terms of the source term by matching coefficients 1 In the asymptotic expansion in eq. (5.1) it is usually convention to include all the momentum dependence in the coefficients, hence to adsorb the log(1/|k|) term into the coefficient with the same power of r. However here we keep the momentum dependence in the term log(r/|k|), as it makes the following argument more clear.
in a near-boundary expansion, and we obtain This implies that the dependence of the Maxwell Green's functions on the renormalization scale k uv appears in a purely real term were we denoted with G µν the part of the Green's functions independent of k uv . The logarithmic correction is thus not visible in the spectral function, nor in the real part of the conductivity. However, due to the form of the RPA-like response function in eq. (3.15), this is not the case anymore when we introduce the double-trace deformation, and the scale-dependent term can also be observed in the real part of the optical conductivity In particular for k = 0, the a x fluctuations decouples from all the others, and the currentcurrent Green's function takes the form and, as mentioned previously, for the neutral system the scale-dependent logarithm only affects the imaginary part of the optical conductivity When we include the double-trace deformation we see from eq. (3.11) that the optical conductivity for a system with Coulomb interactions becomes so that the real part of the conductivity now contains a signature of the scale-dependent logarithm . (5.13) This form of the conductivity, with a logarithmic term depending on a scale ω uv to be determined experimentally, resembles the optical conductivity of Dirac and Weyl semimetals [15]. In the next sections, we give some explicit examples of conductivities in 3+1 dimensions with the double-trace correction. In particular, in section 5.2 we compute the optical conductivity in pure AdS 5 , without backreaction, that can be found analytically. In section 5.3 we instead study the effects of the double-trace deformation on the optical conductivity of the Reissner-Nordström metal.

Exact solution in AdS 5 background
In the simple case of a Maxwell action on a AdS 5 background without backreaction, which is the probe limit where the coupling constant λ → ∞, we can compute the optical conductivity analytically. We use this example to show the importance of the logarithmic term in the optical conductivity when adding a double-trace deformation. In particular, we find that the optical conductivity in the neutral theory matches the results for a noninteracting Weyl or Dirac semimetal obtained from condensed-matter calculations [40,41]. Inserting the deformation then gives us the optical conductivity expected from Weyl and Dirac semimetals with only long-range Coulomb interactions in the RPA approximation.
Since we consider the limit of no backreaction on the metric, we are only interested in the Maxwell action with the AdS background metric (we use the variable ρ ≡ 1/r for notational convenience) The linearized equations of motion of the theory are solved, after fixing the momentum in the x direction, by withk ≡ −(−ω 2 + k 2 ) ∈ R, and H (1) and H (2) are the Hankel functions of the first and second kind, respectively. From now on we focus on the case −ω 2 + k 2 < 0 as we are ultimately interested in the conductivity for k = 0. The solution withC 2 = 0 correspond to the infalling boundary condition and can be used to compute the retarded Green's function. We thus setC 2 = 0 in the following.
For d = 3+1, we can then study the near boundary (ρ → 0) behavior of the field fluctuations that shows the logarithmic term that gives rise to a divergence on the boundary, with γ e 0.577 the Euler's gamma constant. After regularizing the boundary action with a counterterm of the form we can extract the optical conductivity that exactly reproduces the zero-temperature optical conductivity for a non-interacting Weyl or Dirac semimetal, up to a prefactor [40]. In particular, the holographic result in (5.19) describes a system with strong interactions, and the prefactor that, for convenience we have set to one by a rescaling of the action, depends on the coupling constant, contrary to the non-interacting Dirac semimetal calculation. The frequency dependence, however, is the same as it is set by the conformal symmetry of the theory.
When introducing a double-trace deformation the new response function takes the form , (5.20) so that the real part of the conductivity in the modified boundary theory with Coulomb interaction becomes: (5.21) The logarithmic term in the denominator of the conductivity is reminiscent of the conductivity of Weyl semimetals [15]. From (5.21), we can see that for ω → 0 and ω → ∞, σ ∼ ω/ log 2 (ω), and in both limits we have Re[σ] < Re[σ 0 ]. However, there is a range of values for the coupling constant α 2 where we can have Re[σ] > Re[σ 0 ] (considering only ω > 0, where the real part of the conductivity is always positive). Explicitly: This behavior can be observed in figure 12, were we plot the conductivity with and without double-trace deformation for different values of the coupling constant. In all the plots we fixed the renormalization scale ω uv = 100. Note however that, being the only scale in the theory, a change in the value of ω uv simply amounts to a change of scale. Notice also, that there is no plasmon peak in the conductivity. This is an effect of the probe limit, where the Maxwell field decouples from the gravitational background and does not backreact on the geometry. This in turns suppress the sound modes in a neutral theory, and hence the plasmon modes in a charged system.

Conductivity with Coulomb interactions
Here we compute the optical conductivity in d = 3 + 1 boundary spacetime dimensions for the Reissner-Nordström metal with and without double-trace deformation.
In figure 13, we show the optical conductivity for a neutral system for different values of T /µ, and compare it with the AdS 5 analytical solution obtained in the previous section.
In particular, the short-wavelength limit of the theory is determined by the geometry of the bulk spacetime for large r, and we therefore expect to recover the AdS 5 result for high frequencies. This is indeed what we observe in figure 13, and we can notice the logarithmic behavior at large frequency in the imaginary part of the conductivity, were we chose the value of the cutoff-scale ω uv to be the same for different values of T /µ, since the asymptotic behavior depends on the choice of scale. Moreover, we can see that at small frequencies the imaginary part of the optical conductivity of the Reissner-Nordström solution diverges as 1/ω, signaling the presence of a delta function at the origin for Re[σ(ω)]. This delta function is expected since we are dealing with a system with translation invariance, but it is suppressed in the probe limit for the AdS 5 solution as the background is kept fixed [42].
By introducing the double-trace deformation the delta function in the origin disappears, as it is turned into a peak at the plasma frequency by Coulomb interactions. In figure  14, we can clearly observe this peak, that becomes sharper and higher as we lower the temperature. In figure 15 we plot the imaginary part of the conductivity, where we see that  Figure 14: Real part of the optical conductivity for the Reissner-Nordström metal in d = 3 + 1 with a double-trace deformation (α 2 = 1/10). On the right we show an enlarged version of the figure on the left. We clearly see a peak at the plasma frequency, that becomes increasingly higher (left) and sharper (right) as we lower the temperature. Moreover, we notice that the conductivity goes to zero as ω → 0, contrary to the neutral system. the pole in the origin present in the optical conductivity of a neutral system is shifted to the nonzero plasma frequency. Moreover, with the double-trace deformation the real part of the conductivity goes to zero as ω → 0 as is expected from eq. (5.13).

Conclusion and Outlook
In summary, we proposed a general procedure to introduce screening effects of the Coulomb interactions in the holographic description of strongly interacting system. This allows us to study spectral functions of systems with charged particles that contain plasmon excitations, as it is the case in many condensed-matter system of interest. In particular, we numerically studied the effect of this procedure in a Reissner-Nordström theory, where we observed properties expected from traditional condensed-matter calculations.
In d = 3 + 1 we obtained a gapped plasmon mode in the density-density response function, and we showed that the optical conductivity with Coulomb interactions contains scaledependent logarithms, resembling the conductivities predicted in Dirac and Weyl semimetals.
In d = 2 + 1, our main result is a toy model of a system composed of a stack of (spatially) two-dimensional layers, with strong in plane interaction and with coupling between layers governed by the Coulomb interaction. In this model we see that the behavior of the lowenergy dynamics depends on the out-of-plane Bloch momentum p. When p = 0, with the inter-plane distance the density-density spectral function contains a gapped plasmon mode. However, for p = 0, we obtain an 'acoustic plasmon', with a linear low-energy dispersion relation with a renormalized speed of sound. These results shows that the model qualitatively reproduces recent experimental results [18], that suggests that the Coulomb interaction between layers might play a key role in high-temperature cuprate superconductors, proving the necessity of incorporating this interaction in holographic models if we want to study properties of these layered high-T c materials. In a single layer, we instead showed that the dispersion relation assume the form ω ∝ |k|, as observed in graphene.
In conclusion, we have seen that, even with a relatively simple holographic model as the Reissner-Nordström background, the addition of long-range Coulomb interactions presents interesting features that can more closely reproduce experimental results in strongly interacting materials. Therefore, it would be very interesting in future work to apply the double-trace deformation introduced here to different holographic backgrounds such as the holographic superconductor model and Lifshitz solutions, in order to study the behavior of their longitudinal low-energy excitation in the presence of Coulomb interactions. These gravitational theories, contrary to the Reissner-Nordström model, allow for the description of systems with zero entropy at zero temperature. As the ultimate goal is to describe laboratory condensed-matter system, a next important step is to relax the assumption of momentum conservation, as impurities in experimental materials necessarily break momentum conservation. Finally, the addition of a hyperscalying-violation factor to the metric of the Lifshitz solution can be used to study the effect of Coulomb interactions on quantum phases with hyperscaling-violation.

A Gauge Solutions
In the theory considered in this paper, we have two gauge fields, the one-form A µ and the metric tensor g µν . In order to compute the Green's function, we then introduce fluctuations of these fields δA µ ≡ a µ and δg µν ≡ h µν . As explained in section 2, using rotational invariance to fix the momentum along the x direction k µ = (ω, ±|k|, 0, . . . , 0), the remaining fluctuations decouple according to their parity under O(d − 2) acting on x 2 , . . . , x d−1 [29].
In the longitudinal channel, in both d = 3 + 1 and d Using the gauge freedom to choose a gauge where h rµ = 0, as well as A r = a r = 0, the set of coupled fluctuations reduces to This set of fluctuations contains only two physical degrees of freedom, therefore we can only find two independent solutions to the coupled system of linearized equations of motion. However, after setting the r components to zero, we still have some left-over gauge freedom, and the remaining solutions necessary to compute the Green's functions are pure gauge solutions that can be extracted once we know the residual gauge freedom we are left with, as we explain below.

The fields are invariant under the gauge transformation
and diffeomorphisms. The latter gives In order to set the r-components of the metric fluctuations to zero we have to choose a vectorξ ν such that This defines a set of partial differential equations. Working in Fourier space and using rotational invariance to fix the momentum in the x direction, we obtain a set of coupled ordinary differential equations of the form ξ (r, ω, k) + A(r)ξ(r, ω, k) = B(r, ω, k) , (A.6) that has the general solution ξ(r, ω, k) = dr B(r, ω, k)e drA(r) + C(ω, k) e − drA(r) . (A.7) We therefore see that ξ µ is determined up to an arbitrary factor C µ (ω, k)e − drAµ(r) , that corresponds to the left over gauge freedom. In particular from eq. (A.5) we find that the residual gauge transformations are given by ξ r =C r (ω, k)e 2 drΓ r rr , (A.8) ξ α =C α (ω, k)e 2 drΓ α rα , (A.9) where α = t, x, y and there is no summation on the two repeated indices in Γ α rα . In the same way, we want to choose the scalar Λ in (A.3) in order to set the r component of the gauge field to zero. Under the diffeomorphism used to set h rµ = 0, we have that A r + a r transformed as A r + a r → A r + a r + 2g ttξ t Γ t rt A t − g ttξ t A t , (A.10) so we need to choose Λ to be Λ(r, ω, k) = dr A r + a r + 2g ttξ t Γ t rt A t − g ttξ t A t + λ(ω, k) (A.11) where λ(ω, k) is an arbitrary constant in r. Notice that this expression is also invariant under the residual gauge transformationξ t →ξ t + c t e 2 drΓ t rt . Since we have some gauge freedom left, and we know that the linearized equations of motion are invariant under all allowed gauge transformation, we have, given a solutionΦ of the linearized equations of motion, thatΦ transformed under all residual gauge transformation is also a solution. Moreover, we know that Φ = 0 is a solution. We can thus construct pure gauge solutions by choosing n independent residual gauge transformation, i.e., by fixing n linearly independent values of the vector C = (λ, C α ).
For example, in d = 2 + 1 we have from eqs. (A.5) and (A.7) that the gauge where there are no r-components of the fluctuations is defined bȳ ξ r = 1 2 dr h rr e − drΓ r rr + C r e drΓ r rr , ξ t = dr (h rt + iωξ r )e −2 drΓ t rt + C t e 2 drΓ t rt , ξ x = dr (h rx − ikξ r )e −2 drΓ x rx + C x e 2 drΓ x rx , ξ y = dr h ry + e −2 drΓ y yt + C y e 2 drΓ y ry , (A. 12) and we are left with the gauge transformations defined by Λ =λ(ω, k) , Remembering that the behavior of the considered field fluctuations under a gauge transformation is given, in Fourier space, by a x →a x − ikΛ , a t →a t + iωΛ , (A.14) and a x →a x + ikA t /f ξ t , −π/ dp 2π J µ (−k, −p) η µν 2|k| sinh (|k| ) cosh (|k| ) − cos(p ) J ν (k, p) . (B.9) For the boundary theory coming from holography, we use a set of d-dimensional solutions dual to a d + 1 bulk theory, whose boundary terms read, up to a factor of 1/2 so that the (d + 1)-dimensional stack of layers can be modelled as −π/ dp 2π A µ (k, p) J µ (−k, −p) .

(B.12)
Studying the second-order variation to extract the Green's function, we have where a s µ (k, p) is the (d+1)-dimensional source. In addition, G(k), which is independent of the Bloch momentum p, is the Green's function coming from a holographic calculation for a d-dimensional boundary theory, so that G(k)/ has the dimensions of a (d + 1)-dimensional response. When we introduce the double-trace deformation that couples the layers, the total boundary action, second order in fluctuations, then reads Notice that in the limit → ∞ the above boundary action, eq. (B.14), becomes giving the effective two-dimensional response for a single layer.