Coherence effects in disordered geometries with a field-theory dual

We investigate the holographic dual of a probe scalar in an asymptotically Anti-de-Sitter (AdS) disordered background which is an exact solution of Einstein's equations in three bulk dimensions. Unlike other approaches to model disorder in holography, we are able to explore quantum wave-like interference effects between an oscillating or random source and the geometry. In the weak-disorder limit, we compute analytically and numerically the one-point correlation function of the dual field theory for different choices of sources and backgrounds. The most interesting feature is the suppression of the one-point function in the presence of an oscillating source and weak random background. We have also computed analytically and numerically the two-point function in the weak disorder limit. We have found that, in general, the perturbative contribution induces an additional power-law decay whose exponent depends on the distribution of disorder. For certain choices of the gravity background, this contribution becomes dominant for large separations which indicates breaking of perturbation theory and the possible existence of a phase transition induced by disorder.


Introduction
Holographic dualities, that relate classical theories of gravity to strongly-coupled quantum field theories, are now a forefront research area not only in high energy physics but also in quantum information and condensed matter physics. In the latter, it is emerging as a powerful tool to describe universal properties of strongly-correlated quantum systems. One of the main challenges for the application of holographic techniques in this context is the description of disorder, which is ubiquitous in realistic systems and directly responsible for a broad variety of phenomena ranging from momentum relaxation to quantum interference leading to different forms of localization [1][2][3][4][5]. The introduction of disorder in gravity backgrounds with a negative cosmological constant relevant for holography requires the solution of spatially inhomogeneous Einstein's equations, in general a difficult task.
Different approximation schemes have been proposed to make the problem technically tractable while keeping some of the expected phenomenology related to the introduction of disorder. For instance, momentum relaxation, a rather general consequence of any form of disorder, can be achieved by adding a massless scalar [6][7][8][9] that depends linearly on the boundary coordinates. Since the scalar field only couples to gravity through its derivatives, translation invariance is broken by the background but the equations of motion are still independent of the spatial coordinates which facilitates substantially the calculation, in many cases analytical, of transport properties.
As was expected, the electrical conductivity is always finite and depends directly on the strength of the translational symmetry breaking characterized by the slope, with respect to the spatial coordinates, of the scalar field in the boundary. For weak momentum relaxation, the electrical conductivity reproduces the expected phenomenology of Drude's model, which includes a peak at low frequencies, remnant of the broken translational symmetry, followed by a decay for higher frequencies. For stronger relaxation the Drude peak is suppressed, leading to an incoherent 'bad-metal' behaviour [10][11][12][13]. However even in the limit of infinite relaxation no insulating behaviour is observed. Recently, models that consider the coupling of the scalar field to the gravity and the Maxwell term managed to reproduce a vanishing conductivity in the limit of infinite relaxation [14][15][16]. Yet, in this limit the effective charge in these models vanish, so strictly speaking they cannot be considered insulators. Other effective models of momentum relaxation in holography include the memory matrix formalism [17][18][19][20], helical and Q-lattices [21][22][23] and massive gravity [24][25][26][27]. Similar results [28][29][30][31][32][33][34] have been obtained even for Maxwell fields with a spatially oscillating chemical potentials in the boundary leading to inhomogeneous Einstein's equations. We note that effects such as localization are precluded by design since a random but homogeneously distributed, and therefore delocalized, chemical potential is an input in this approach.
Disorder has also been introduced at the level of the action by a random source coupled to the dual conformal field theory operator [35,36]. By using the replica trick, it is possible to integrate out the random coupling resulting in a double trace deformation of the non-disordered theory. For marginal perturbations, renormalization group techniques suggest the existence of logarithmic corrections in the two-point correlation function that spoils conformal symmetry. Interestingly, this is in agreement with the expected behavior of certain two dimensional conformal field theories perturbed by disorder (see [37] and references therein).
Another popular approach is to consider a random field in the boundary, in most cases a random chemical potential (source) for gauge (scalar) fields, but neglecting the backreaction in the gravity background (probe limit) [38][39][40][41][42][43][44][45][46]. Analytical attempts to go beyond this probe limit have found logarithmic infra-red (IR) divergences in the gravity background, signaling the breakdown of perturbation theory [47,48]. A resummation scheme for the divergent expansion was proposed in [49][50][51][52] resulting in an averaged gravity background with an emerging Lifshitz scaling symmetry, which was shown to be a generic feature of marginal disordered deformations [53,54]. For geometries with a horizon, the most general result in this context is that of Gauntlett, Donos and collaborators [29,[55][56][57][58][59] who found closed expressions for the dc conductivity and other averaged transport coefficients in generic inhomogeneous backgrounds. As a direct consequence, bounds on the electric and thermal conductivity were proposed [60][61][62][63].
The conclusion of this more direct approach is similar to the one from the phenomenological models of momentum relaxation discussed previously, namely, it is not possible to reach an insulating state which is believed to be a distinctive feature of strong disorder in condensed matter systems. Even simpler coherence effects such as weak-localization [64], which are precursors of a metal-insulator transition, have not yet been clearly identified.
In this manuscript we propose a new approach to model disorder in holography which has the potential to reproduce some of these coherent effects. We switch perspectives and consider the effect of a disordered geometry with no horizon on a probe scalar field. This is accomplished by considering a family of three dimensional random geometries that solve Einstein's equations exactly 1 and neglecting the backreaction of the scalar in the geometry. This family is indexed by a parameter which we take to be a random function of the boundary coordinates. The scalar field feels the geometry as an effective inhomogeneous coefficient in the equations of motion. This is reminiscent, though we cannot establish a precise mapping, of a one-dimensional wave equation with a random refractive index [4,5] plus additional terms that control the evolution in the radial direction which are related to interactions in the boundary.
According to the holographic dictionary, the geometry is dual to a strongly-coupled disordered plasma living at the boundary, and the scalar field sources a boundary dual operator. We investigate numerically and analytically the properties of this disordered plasma by looking at one and two-point functions of this scalar operator. For the onepoint function our main result is the observation of coherence effects, due to interference between an oscillating source and the random geometry, that, in some cases, leads to the strong suppression of oscillations even for weak disorder. The contribution to the two-point function for a weak random Gaussian geometry is still a power-law decay for large distances with an exponent that depends on both the scalar mass and disorder correlations modeled by a non-trivial power spectrum. In some cases, this correction becomes dominant for large distances which suggests the breaking of perturbation theory and the possible transit of the system to a new disorder-driven fixed point.
The manuscript is organized as follows. In the next section we introduce the geometry we will be studying and discuss the equation of motion for the probe scalar. In section 3 we solve these equations analytically in the limit of a weak-disordered background. The one and two-point functions of the dual operator are computed for different choices of both sources and inhomogeneous geometries. In section 4, we solve the equation of motion 1 The fact that we can find an exact solution of Einstein's equations involving a free function is due to the fact that all vacuum solutions are pure diffeomorphisms in three dimensions. The generalization to higher dimensions would involve solving the equations numerically.
by numerical techniques and compute in certain cases the one and two-point correlation function of dual scalar. Section 5 is devoted to a comparison of our results with previous approaches to disorder in holography. We conclude in section 6 with a summary of results and ideas for future work. The appendices offer a wider discussion of technical points which are used throughout the main body of the paper.

Setup
In this section we introduce our objects of study. First, we introduce a family of geometries in d + 1 = 3 spacetime dimensions. This family is characterised by an arbitrary function which we can take to be random. We next introduce a minimally coupled scalar field and discuss its equation of motion and associated boundary conditions. The aim is to use this field to probe the properties of the inhomogeneous geometry, which holographically can be interpreted as a strongly-coupled disordered field theory.

Geometry
Solutions of the vacuum d+1-dimensional Einstein's Equations with a negative cosmological constant have been classified in the pioneering work of Fefferham and Graham [65,66]. For d > 2, we can integrate Einstein's Equations in a neighborhood of the boundary by requiring that the Weyl tensor vanish. This condition constrains the boundary metric to be conformally flat. However, in d = 2 the Weyl tensor vanishes exactly, leaving the conformal class of the boundary metric arbitrary [67].
In this manuscript we will be mainly interested in the following family of metrics defined by a global coordinate patch x a = (ρ, t, x) as where g tx (x) is an arbitrary function of the boundary coordinate x. In line with our discussion above, it is easy to check that this family satisfies Einstein's Equations in d = 2 dimension with a negative cosmological constant for any g tx . In these coordinates, the conformal boundary is located at ρ = 0 and the induced conformal metric g (0) is given by From a holographic perspective, this space-time encodes the degrees of freedom of a strongly-coupled field theory living on the boundary metric described above. However, the dual stress tensor vanishes since this change in the boundary metric does not induce subleading terms in the bulk metric. This resembles the situation we encounter in the axion model of [7], where the marginal, spatially dependent, scalar sources do not excite a vev. In the next sections we will be interested in studying the family of geometries in Eq.(2.1) for different choices of g tx . In particular, we will be interested in the case where g tx (x) is a random Gaussian process, as introduced in Appendix A. We will study the dynamics of a minimally coupled scalar field as a way of probing the effects of the disordered geometry. The aim is to get an insight into the nature of the strongly-coupled disordered dual field theory.

Scalar Field and Equations of Motion
We consider a probe scalar field ψ of mass m minimally coupled to the geometry in Eq.(2.1). The equation of motion is given by (∆ g − m 2 )ψ = 0, where the curved Laplacian can be written in a chart x a as ∆ g = 1 Since ∂ t is a killing vector for Eq.(2.1), we can restrict our attention to static configurations ψ = ψ(ρ, x). The equation of motion thus reads We are interested in solutions satisfying the following boundary conditions, where we have defined ν = √ 1 + m 2 . The first boundary condition assures regularity of the solution at the Poincaré horizon ρ = ∞. As discussed in Appendix B, the second boundary condition defines a field s(x) living in the boundary ρ = 0 with conformal dimension ∆ − = 1 − ν. This field sources a dual boundary operator O(x) ∼ ρ − ν+1 2 ψ of conformal dimension ∆ + = 1 + ν.
Even for simple choices of g tx , Eq.(2.2) remains largely intractable analytically. However when g tx is small we can compute perturbative corrections to the plain AdS 3 result analytically, helping us to build an intuition of the effects of the weakly disordered geometry. We will recur to numerical methods to test the analytical prediction and also to explore the region of stronger disorder not accessible to an analytical treatment.

Perturbative analytical calculation of one-point and two-point scalar correlation functions
In this section we study Eq.(2.2) perturbatively for different choices of g tx . For convenience, we will set m 2 = − 3 4 throughout this section, which is equivalent to choosing ν = 1 2 . Note this choice respects the Breitenlohner-Freedman bound, m 2 ≥ −1 in d = 2, [68,69]. We will divide our discussion by order in perturbation theory, and focus in two observables: the expectation value of the dual boundary operator (one-point function) and the two-point function. We can obtain both these quantities by solving the perturbative equations with the appropriate boundary conditions.

Zeroth Order
To set up the perturbative analysis, let g tx (x) = g(x) for a parameter 1 measuring the amplitude of g tx fluctuations. For = 0, Eq.(2.2) reduces to the equation of a massive scalar field in AdS 3 given by 2π e ikx f k (ρ), the equation above reduces to which has general solution Regularity at the Poincaré horizon Eq.(2.3) requires that a k = 0 for k < 0 and b k = 0 for k > 0, which can be written compactly as f k (ρ) = a k ρ 1/4 e −|k| √ ρ for a different constant a k . Now letting s(x) = R dk 2π e ikx s k , boundary condition Eq.(2.4) can be imposed coefficient wise to yield a k = s k . The full solution therefore reads Note that this depends directly on the Fourier components of the source. We are interested in a few particular cases that we outline below. s n cos(q n x + γ n ), where s n , q n and γ n can be freely chosen.
Noting that the finite sum can be exchanged with the integral and following the same steps as above mode-wise we find ψ (0) (ρ, x) = ρ 1/4 N n=1 s n e −|qn| √ ρ cos(q n x + γ n ). The one-point function thus reads O(x) = − N n=1 s n |q n | cos(q n x + γ n ).
Delta source and two-point function As discussed in Appendix C, the boundary-tobulk propagator is given by solving the equations of motion with a delta source, We thus have, which follows the shape of a Lorentzian (or Cauchy) distribution. Following the discussion in Appendix B, the boundary two-point function can be obtained by (3.6) Note this is not the expected result for an operator of conformal dimension 2∆ + = 3 but rather for an operator of conformal dimension 2∆ + = 2. This is because, by considering a static field, we effectively reduce the conformal dimension of the problem. Since we will be interested in static inhomogeneous configurations in what follows, this is the object we will be computing corrections for.

Second Order
Note that since g tx appears only quadratically in Eq.(2.2), there are no non-trivial order one corrections to ψ (0) . It is thus sufficient to consider ψ = ψ (0) + 2 ψ (2) + O( 4 ). Inserting into Eq.(2.2) and expanding up to second order leads to Note that the zeroth order solution act as a source for the perturbative correction. Since the full solution has to satisfy the boundary conditions, we have to apply them order by order. For instance every term in the -expansion has to be regular at the Poincaré horizon. However since we have enforced boundary condition Eq.(2.4) at zeroth order, we have to set lim ρ→0 ρ ν−1 2 ψ (2) = 0 for the inhomogeneous geometry not to correct the fixed boundary source s(x). Summarizing, we have to solve Eq.(3.7) subjected to If g tx is of Schwarz class, we can attempt to solve Eq.(3.7) in Fourier space as we did for the zeroth-order result. Letting ψ (2) 2π e ikx f k (ρ) and g(x) = R dk 2π e ikx g k , we can rewrite Eq.(3.7) in Fourier space as where the right-hand side has been evaluated applying the convolution theorem with the zeroth order solution Eq.(3.4) 2 . This integral-differential equation cannot be solved in a closed form. In the following subsections we discuss solutions for specific inhomogeneous configurations. For each choice of g, we can consider different choices of source and study the interplay between the source, the geometry and the resulting expectation value and two-point function of the dual boundary operator.

Constant geometry
The simplest example is given by taking g(x) = 2π g to be a constant. In this case g k = gδ(k) and for generic source Note that the linear differential operator on the left-hand side is exactly the same as for the zeroth-order equation. This is generic and hold at all orders in perturbation theory. The only difference is the source term on the right-hand side. By linearity, the general solution will be a linear combination of the solution for the homogeneous equation plus a particular solution. As before, the homogeneous solution has one exponentially diverging piece which should be set to zero by regularity at the Poincaré horizon. This leads to where a k is an integration constant. Close to the boundary ρ = 0, the solution behaves as and therefore boundary condition Eq.(3.9) imposes a k = − g 2 4 s k . Finally, we can write As in the zeroth-order solution, we discuss below a few cases of interest.
Constant source Let s(x) = s be a constant. Inserting in the above leads trivially to ψ (2) (ρ, x) = 0. Therefore the off-diagonal constant metric does not affect the zeroth order boundary one-point function.
Oscillating source Let s(x) = s cos qx. Inserting in the above leads to At the boundary ρ = 0 this induces a correction to the zeroth-order one-point function which is proportional to g, s k cos (q n x + γ n ). This case is similar to the above, since we can integrate term by term in the sum to give, The correction to the boundary one-point reads Again, it represents just a renormalization of the amplitude of the zeroth one-point function.
Delta source and two-point function Recall that, as discussed in the previous section and in the Appendix C that the boundary-to-bulk propagator K(ρ; x, y) can be computed by solving the equations of motion with s(x) = δ(x). We thus have Evaluating at the boundary leads to a correction to the first zeroth-order two-point function, Note that for all the cases above the constant off-diagonal metric element perturbatively decreases the amplitude of the boundary one-point function and two-point function. Since in perturbation theory the equations are linear this case can be interpreted as the mean result for a inhomogeneous geometry. The amplitude damping raises the question on whether in a non-perturbative setup the geometry can effectively suppress the boundary correlation functions.

Oscillating geometries
We now consider a generic superposition of N oscillating modes, g(x) = N n=1 A n e iωnx for ω n ∈ R. For example, an interesting particular case is ω 1 = −ω 2 = ω, A 1 = A 2 = g 2 and A n = 0 for n > 2 which correspond to g(x) = g cos ωx. The Fourier modes are given by a Dirac comb g k = N n=1 A n δ(k − ω n ) and for a generic source the equations of motion in Fourier space read, By linearity, we can solve the equation above term by term. The regular solution at the Poincaré horizon is given by the homogeneous solution Eq.(3.3) plus the sum of each individual particular solution Taking the Fourier transform, we arrive at an implicit solution for the generic source As before, we now analyse some interesting particular cases where the above integral can be done explicitly.
Constant source Let s(x) = s be a constant, i.e. s k = δ(k). As before, the above integral trivially gives ψ (2) = 0. There are no corrections to the boundary dual operator.
Plane wave source Let s(x) = se iqx , i.e. s k = sδ(k − q). The Fourier transform is given by Close to the boundary, this gives a correction to the dual one-point function There are a couple of particular cases of special interest, Consider q, ω ≥ 0 and define the relative correction ∆( , where O (0) is the zeroth order one-point result. We can identify three distinct regimes. For q > 2ω, the relative correction to the one-point function is positive and simply proportional to the geometry ∆(x) q≥2ω = − 2 2 g(x) 2 = − 1 2 g tx (x) 2 , while for q < 2ω the result depends explicitly on the relative strength of the modes, For q < ω the relative correction will be positive, while in the window ω < q < 2ω it becomes negative. Therefore the inhomogeneous geometry can either suppress (q > ω) or enhance (q < ω) the expectation value of the dual operator, depending on the coherence between the modes. The discussion is summarized in Fig.1. • g(x) = g cos ωx (A 1 = A 2 = 1 2 g, ω 1 = −ω 2 = ω, all other zero): By linearity, this case can be obtained by summing the above with the reversed q → −q. The relative correction ∆( is given by Note that, different from the above in general it is not possible to rewrite the relative correction as an explicit function of the geometry. First, lets consider the case q, ω ≥ 0. Then q + 2ω ≥ 0 and we have two subcases: q ≥ 2ω and q < 2ω. These are given by The two modes e 2iω and e −2iω can interfere constructively or destructively depending on the relative sign of q and ω. Note that the correction is symmetric under ω → −ω as expected. The q, ω < 0 case is similar up to a change of sign. Oscillating Source We now build on the previous results to analyze more intricate cases. Let s(x) = s cos qx. Linearity of Eq.(3.20) together with the example above can be used to get the following correction, Consider the particular subcases: The correction can be conveniently rearranged to give which is, as expected, a real result. However note that the source modes are now coupled with the geometry. We can still write the relative correction ∆(x) by dividing by the source, and being careful to take into account that when the source vanishes the result is zero and not divergent.
This can be analyzed as before, giving Interestingly, as with all the examples we considered above, the relative correction for the case q > 2ω factors into a contribution that only depends on the underlying geometry. This can be interpreted intuitively as follows. If the wavelength of the source (∝ q −1 ) is much smaller than the wavelength of the geometry, the dual operator will not 'feel' the inhomogeneity, leading to a coupling similar to the constant geometry case discussed in Section 3.2.1. In Fig.2 we plot the correction for different configurations of (ω, q) and inhomogeneity strength .
A n cos ω n x. This polychromatic case is a direct extension of the above. It can be obtained by taking N → 2N and taking A n → 1 2 A n for n = 1, . . . , 2N , ω n → ω n for n = 1, . . . , N and finally ω n → −ω n for n = N + 1, . . . , 2N . By carefully splitting the sum and rearranging the terms, we (3.32) Note that since both terms are symmetric under the exchange n ↔ m, we can split the sum into a diagonal and an upper diagonal part, The expressions can suggest that the diagonal part in ∆ < diverges. However this is not the case, since the denominators also vanish, and we have to take the limit carefully.
This polychromatic case is of particular interest since it can be used to simulate a discrete representation of disorder, as discussed in Appendix A. Defining ∆ω = π N a for some constant lattice spacing a > 0 and letting A n = √ ∆ω, ω n = n∆ω for n = 1, . . . , N and adding random phases i.i.d. uniformly γ ∈ [0, 2π) in the cosine arguments, g(x) become a discrete representation of the Gaussian random process for large N 1. As we will discuss in the next section, continuous disorder is more subtle, and we have to take the average at an early stage to make progress. Moreover, it is harder to implement it numerically since it can usually have discontinuous derivatives. For these reasons, this implementation is a useful representation and is specially suited for holography calculations, having been used in many previous works in the literature [42,49,54]. In this case all modes ω n are positive, and we can assume without loss of generality that q > 0. For convenience, we separate the sum in the correction in a diagonal and a non-diagonal part, ∆ = ∆ d + ∆ nd with where for convenience we abbreviated ω ± nm = ω n ± ω m . Note that only the constant term in the diagonal part survives averaging over the i.i.d. phases γ n as discussed following Eq.(A.5) in Appendix A. As expected, we reproduce the results for the constant geometry in Section 3.2.1 on average. In Fig.3 we plot the expectation value of the dual operator for N = 10 modes and fixed source q = √ 2 for a fixed realization of phases. Delta source and two-point function We now consider the case s k = 1 that leads to the corrections to the boundary-to-bulk propagator K. The integral Eq.(3.20) can be done explicitly, Recall that the two-point function of the dual boundary operator is obtained by evaluating the boundary-to-bulk propagator at the boundary. Recalling that the zeroth order result is given by 1 πx 2 , this yields to a relative correction to two-point function given by Note that both terms are symmetric over n ↔ m. Thus we can split the sum into a diagonal term and a term where n < m, Consider the following interesting examples, • g(x) = g cos ωx.
In this case we simply have which is proportional to a sinc function. Recall that we are looking at correlations of the point x with the origin. Thus the correction induced by the oscillating geometry suppress correlations of points closer to the origin. Note that for ω 1, the peak is sharper, while for ω 1 it is broader. Thus the bigger the wavelength of the inhomogeneous geometry, the less is the correction localized at x = 0, see Fig.4. We can check that the limit ω → 0 reduces to the constant case discussed in Section 3.2.1. • g(x) = N n=1 A n cos (ω n x + γ n ).
The result above can be easily generalized for many modes. Following the discussion in for the single cosine source, we have: Interestingly this discrete random process has mean zero.

Disordered Geometries
We now study the case where g(x) is a centered random Gaussian field. As discussed in Appendix A, we can decompose g in its spectral modes g(x) = R dk 2π g k e ikx , with g k also given by a centered Gaussian process with E[g k ] = 0 and E[g k g l ] = g 2 δ(k − l). Equation (3.10) becomes a stochastic integro-differential equation, and is still intractable. However since it depends on the square of the geometry we can consider its non-trivial average. Definingf k = E[f k ] and averaging over g yield As before, the constant source case give a trivial result. We now analyze non-trivial settings.
Plane wave source Let s(x) = se iqx with q > 0. Integrating Eq.(3.40) gives The solution satisfying boundary conditions Eq.(3.8) and (3.9) is given bȳ Note that the case k = q can be treated by taking the limit of the above. Fourier transforming back, (3.43) The correction to the one-point function is given by

44)
Oscillating source Let s(x) = s cos qx. As before, by linearity of Eq.(3.40), we can build the oscillating case by summing two plane waves. Using that E 1 (z) = E 1 (z) we can simplify where Ci(x) is the cosine exponential function 3 . This leads to ∆(x) = − 2 g 2 sin qx qx + Ci(qx) .
The regular solution at ρ = ∞ for the above is given by where Ei(x) = ∞ −x t −1 e −t dt is the exponential integral functions. Close to the boundary ρ = 0 we have Therefore in order to preserve the zeroth order boundary condition we need to set a k = 0. Note the appearance of two divergences: one proportional to ρ −1/4 and the other proportional to log k √ ρ. As we will see next, they recombine when taking the Fourier transform and lead to a finite result.
Using the following Fourier transform that can be calculated from the definition of Ei(x), Letting a = √ ρ > 0 and taking into account the constant term leads to (3.52) Note that the first term in Eq.(3.47) leads to a delta function with opposite sign that exactly cancels the one coming from the derivative of the sign function. As discussed in Appendix C, this result is finite close to the boundary Therefore the correction to the two-point function can be written as We shall see that the exponent for which the correction blows up as x → 0 depends on both the mass of the scalar field (here fixed to m 2 = −3/4) and on the type of disorder. In the following subsection, we will relax these conditions and explore the dependence of our results on these parameters.

Comments on other masses and correlated disorder
In this subsection we discuss how the results for the averaged two-point function generalize to different masses and correlated disorder.
We start by considering different masses. It is easy to check that for a generic mass parametrized by ν 2 = m 2 + 1 the zeroth order solution of Eq.(3.2) that satisfies the boundary conditions Eqs.(2.3),(2.4) is given by, The boundary-to-bulk propagator is, as before, obtained by setting s k = 1 and computing the Fourier transform, giving (3.56) Taking the near boundary limit, we obtain the expected dual two-point function, which the the expected result for the two-point function of a conformal operator with mass dimension ∆ + . The minus one factor is, as before, due to the time dimensional reduction. The calculation of corrections induced by the disordered geometry are exactly as before, with the only difference that we use the general solution above as a source in Eq.(3.7). Averaging the right-hand side, we get the general equation which reduces to Eq.(3.46) for ν = 1/2. The general solution for the above is a combination of hypergeometric functions. So next we consider other values of ν for which the calculations are less cumbersome. Take for instance ν = 3/2. In this case the solution satisfying the boundary conditions is given by The Fourier transform of the above can be computed exactly in the same way as in Section 3.2.3 and is given by Noting that for ν = 3/2 we have ∆ + = 1 + ν = 5/2, the dual two-point function is given by This result follow exactly the one for ν = 1/2, with the only difference that now the zeroth order result has a different mass dimension. It is not hard to check that the same is true for ν = 5/2, where the averaged two-point function is given by 3π |x| 1 |x| 6 , ν = 5/2, (3.62) and now 2∆ + − 1 = 6. Although we have not manage to prove the general result, it seems that white noise always induce the same relative corrections in the two-point function. As we will discuss next, this is not completely surprising. For higher masses the dual operator is a more relevant deformation, but we are not changing the mass dimension of disorder. To see this explicitly, note that the metric element g tx should be dimensionless. Since we are imposing g tx = Note that since σ 2 k is a variance, it needs to be a positive definite function. Generalising the previous discussion, this gives This is precisely what we found in the previous section: corrections to the two-point function decay faster in the IR limit |x| → ∞ than the leading order result. Choosing α > 0 will only make disorder more irrelevant. It is easy to check that for α > 0, we get subleading powers of ρ which do not contribute to the two-point function.
We now explore the case α = −2, when disorder becomes relevant. For simplicity, lets consider again ν = 1/2. In this case, the integral in the right-hand side of Eq.(3.40) can be written as where E α (z) = z α−1 ∞ z dt t −α e −t is the generalised exponential integral function, which in the last line we related to the exponential integral through the following recursive relation (3.65) together with E 1 (z) = −Ei(−z). The regular solution at the Poincaré horizon which does not modify the source will be given by where in the last equality we defined ω = k √ ρ. We can compute the Fourier transform of the above in following the same recipes as discussed above. The boundary-to-bulk propagator is given by 68) and the correction to the boundary two-point function thus given by (3.69) Note that different from the results in the previous section, the correction induced by the correlated Gaussian disorder dominates the decay of the propagator at the IR, |x| → ∞.
For any fixed amplitude > 0 and for |x| > 1 2 π the correction becomes more important than the zero order result, which signals a breakdown of perturbation theory.

Numerical Analysis
We have also performed a numerical analysis which allows us to confirm and extend our previous results beyond perturbation theory. Our first step is to compactify the interval of the radial coordinate. Following [49], we define y via In this new coordinate, the boundary is located at y = 0 and the Poincaré horizon at y = 1. We redefine the scalar as We can check that the near boundary behavior can be expressed in terms of χ as Here χ (0) and O correspond to the source and vev of the dual operator. Comparing with the pure AdS solution given in (3.3), we see that, at least for g tx = 0 modes with non-zero momentum decay or grow exponentially at the horizon. After turning on g tx , the near horizon asymptotics become harder to analyze, but by continuity with the pure AdS case we demand that the field vanishes there. We solve the equation of motion by discretizing it on a homogeneous grid along both the radial coordinate y and the boundary coordinate x, and solving the resulting matrix equation in Mathematica. Since most of the variability occurs along x, we use N x = 450 grid points along this direction, and N y = 50 for the radial direction. After solving the wave equation, we extract the vev by taking a radial derivative of the solution at the boundary, as seen in (4.3).
We compute approximations to the one and two-point functions, concentrating on the case of disordered geometries. We do so by employing the spectral representation of disorder discussed in Appendix A.3. More specifically, we take g tx to be given Gaussian by using (A.4) with σ(k n ) = 1 with suitable modifications described in more detail below.

One-point correlation function
We consider one-point function of the scalar in the presence of the source (4.4) As explained above, the vev corresponds to the derivative of the field at the boundary, as given by (4.3). Since the periodicities of the source and geometry need to fit in the same computational domain, the expansion of the disordered geometry will contain terms cos(2nk 0 x + γ) with integer n, see Eq. (A.4). As seen in Section 3.2.2, these are potentially problematic since the commensurability of the source with the geometry induces extra near horizon divergences. To avoid this undesirable behaviour, we only consider cosines of odd momentum in the sum Eq. (A.4). We shall see that with this modification we can still obtain some generic features of disorder which match well with the perturbation theory results for small disorder amplitude, and extend them to higher values. We extract the one-point function for k 0 = 5, and N = 20 for varying values of the disorder amplitudeV as defined in Appendix A.3. For everyV , we generate a random geometry by providing random phases in Eq.(A.4). Once we have obtained a large number of them, we take the arithmetic average at each point x. We write the average vev as a suppression factor η(V , x) times the translational invariant result, given by O 0 := O |V =0 = −k 0 cos(k 0 x). Hence, we define η by For small disorder amplitudes the x-dependence of η(V , x) is very mild. However, at larger amplitudes the x-dependence of η becomes important. We show this in Fig. 6 where  In order to estimate the overall suppression, we track the value at the peak η(V , x = 0). We show our results in Fig 7. For smallV , we observe that the averaged one-point function displays the quadratic behavior obtained in perturbation theory, although with different proportionality constant. At larger amplitudes, the power-law behavior becomes milder. Moreover, we are able to fit the distribution of values of the vev at the peak with a Gaussian centred at the average value, see Fig 8.

Two-point correlation function
We now obtain an approximation to the two-point function in the presence of a disordered geometry. In principle, this entails a highly expensive calculation which requires inserting arbitrary sources at different points and taking the variation of the action with respect to them in the presence of a spatially dependent geometry. In order to gain some insight on the behavior, we consider the more tractable calculation corresponding to the two-point function G(x, 0), which as explained above can be obtained by inserting a delta function source at y = 0. In order to regularize the delta function, we follow the strategy of [70], i.e. we take as a source the boundary-to-bulk propagator evaluated at a small cutoff. We stress that we need to take into account the fact that in our numerics the x coordinate is periodic, which changes the form of the boundary-to-bulk propagator even in the absence of disorder V = 0. In fact, it is easy to show that for a box of length 2π/k 0 the boundary-to-bulk propagator is given by Note that here y refers to the radial variable introduced in Eq.(4.1). Therefore, we approximate the two-point function G(x, 0) by the one-point function obtained in the presence of the source at small δ. To test our approximation scheme, we first derive the results for pure AdS 3 , V = 0. As expected, this approximation fails for x ≈ 0. In particular, the so-obtained vevs become very large and negative for small enough x. However, the results near the edge of the computational domain x = π are well-behaved, and match well with the analytic result, see Fig 9. Therefore, we will be able to extract meaningful results away from the cores in this region 4 . In order to capture the behavior away from the core, we take the spatial average in the interval (2π/3, π), which we denote bỹ This gives an estimate of the suppression of the two-point function in the presence of disorder. We plot our results for this quantity with k 0 = 5, N = 20, δ = 0.01 and varyingV in Fig 10. Once again, we obtain a quadratic dependence of the suppression with the disorder amplitude.

Comparison with previous results in the literature
In the introduction, we reviewed different approaches to introduce disorder or spatial inhomogeneities in the holography literature. Now we discuss similarities and differences with the one introduced in this paper. A direct comparison is in general not possible in most cases. For instance the prediction for transport coefficients of [29,[55][56][57][58][59] requires the existence of a horizon and therefore are not applicable to our case. Even with no horizon, a comparison could be problematic because the assumption of a random chemical potential prevents, at least for a non-random background, any coherence effect. The addition of a random source, investigated in [35,36], is, to the best of our knowledge, not clearly connected to our approach. Instead of turning a source for the charge density, in our model we introduce a source for the stress energy tensor. From the gravitational perspective, we expect this to have a more significant effect because all fields couple to gravity. The equivalent field theory statement is that all operators propagate on the fixed boundary geometry. Moreover the geometry in this approach is not random, so no coherence effects are expected to be observed. It is an open question whether the observation of non-perturbative logarithmic corrections for marginal disorder in the two-point function, reported in [35,36], could occur in our setting. In order clarify these issues it would be necessary to carry out a full renormalization group analysis, beyond the scope of the paper, for the parameters for which the perturbative contribution from the disordered background becomes marginal. The result of such calculation would not only shed light on the existence of logarithmic corrections but also on the possible existence of a metal-insulator transition.
The approach closer to the one studied in the paper is maybe that of [47,49,53,54] where a spatially random chemical potential, or scalar, in the boundary, backreacts in the gravity background that becomes inhomogeneous as well. However there are still important differences. At least perturbatively, interesting coherence effects are strongly suppressed, even if no horizon is present, because the only independent source of randomness comes from the scalar or chemical potential whose profile in the boundary is fixed by boundary conditions.
- 25 -In this manuscript we have proposed a new approach to study disordered holographic field theories. We have computed numerically and analytically corrections due to a weakly disordered gravity background in the one and two-point function of the scalar dual boundary operator for different choices of source and random component of the geometry g tx (constant, a superposition of plane waves and finally a Gaussian random function). The main results can be summarized as follows: • A constant geometry induces a negative constant correction that decreases the overall amplitude of the scalar one and two-point functions.
• The corrections induced by an oscillating geometry have a richer behaviour, and depend on the interaction with the source. In the simple case where the source is also an oscillating function, the perturbative correction to the one-point function can be positive or negative, depending on the relative sign of the geometry and source frequencies. For instance, if the source frequency is much bigger than the geometry's, the correction is similar in form to the constant case outlined above. However if the source frequency is smaller than the geometry's, the correction flips sign, adding constructively to the zeroth order one-point function.
• An oscillating sinusoidal geometry induces an oscillating but decaying correction to the two-point function, with an envelope ∼ |x| −α that depends on the scalar mass.
• The case in which the geometry is a superposition of oscillating modes is given, by linearity, by the sum of the single frequency results mentioned above.
• We have identified coherence effects between a sinusoidal source and the weakly random geometry introduced by a spectral decomposition. In certain region of parameters, the one-point function, which is sinusoidal, in the absence of disorder, becomes completely random even in the limit in which perturbation theory applies.
• The averaged corrections to the two-point function in the presence of a delta source and a weakly random, Gaussian distributed, gravity background, is negative and, for large distances, decays as a power-law with an exponent that depends on the type of disorder and the scalar mass. We have identified a range or parameters for which perturbation theory breaks down, as the power-law decay is slower than in the nonrandom case. This suggest an instability to a novel disorder driven fixed point which could eventually lead to a metal-insulator transition in the system.
Finally we mention some ideas for future research. First, we have only studied one instance of random geometry. Asymptotic AdS 3 solutions of Einstein's Equations have been classified in [67], and there are other families in which disorder could be introduced in a similar fashion. It would be interesting investigate whether introducing disorder in other metric components would lead to a similar universal behavior, and if not so, to identify the physics behind these differences. Second, since temperature tends to suppress coherence effects in disordered systems, in this work we have only studied holographic field theories at strictly zero temperature. But finite temperature solutions in three dimensions have also been classified [71], and therefore our work could be also generalized in this direction. In particular, it would be interesting to study the low temperature regime and compare it with our results. Third, we have only considered Gaussian distributed disorder with deltalike correlations. However our formalism is easily generalizable to more general forms of disorder where correlations between points becomes important. Fourth, a renormalization group treatment, feasible for marginal random perturbations, would shed light on the existence, or not, of a novel non-trivial disordered driven fixed point which signals an instability toward a metal-insulator transition in the system.
is commonly known as the power spectrum of the random field. For the simple case where the power spectrum is a constant σ 2 which means the distribution of f ( x) in real space is also Gaussian. Conversely, considering a non-trivial power spectrum lead to non-trivial correlations for the random field. Therefore the spectral representation of a random field is a convenient way of generating non-gaussian distributions while still working with Gaussian objects in Fourier space.

A.2 Cutoffs
However this construction is not very useful for actual applications. For instance, note that as x approaches y we get an ultra-violet divergence. This for instance can be very inconvenient in the case we are treating, since in our equations we have a lot of terms that go as 'disorder squared' at the same point. To remediate this problem, we will resolve this long wavelength divergence by introducing an ultra-violet (UV) cutoff λ and integrate only over modes | k| < λ. A pictorial way to interpret this cutoff is to say that λ introduce a length scale a = π/λ that corresponds to an underlying lattice. Modes with frequencies below this scale are then ignored. This would imply for instance that E[f ( x) 2 ] =V 2 /a. As we take the lattice spacing a → 0 we recover the expected UV divergence.
While UV divergences are a consequence of the way we introduce disorder, there can be IR divergences that are emergent in the problem, and indicate a change of behaviour in the system. Or in terms of the renormalisation group: the system flows towards a disordered fixed point. To resolve these divergences one usually introduce a box of size L, and in the end of calculation one aims to study how the system behaves as L is increased. IR singularities in the thermodynamic limit L → ∞ indicate a flow towards a new phase.

A.3 Discrete
Although the continuum implementation simplifies analytical calculations, numerically one needs a discretisation that takes into account the aforementioned observations. From now on we restrict ourselves to the case of interest, namely d = 1. Note that according to Eq. (A.2) the norm of |f ( k)| is drawn from a Gaussian distribution, while the phase is drawn from the uniform distribution on [0, 2π]. This remark leads to a useful discrete representation of the continuous spectral decomposition. We start by discretising uniformly our box of size L in N intervals of size a, i.e. L = N a. This implies the quantisation of the modes k n = n∆k with ∆k = π N a = λ N . The thermodynamic limit then becomes N → ∞. The spectral decomposition Eq. (A.1) becomes, A n cos(k n x + γ n ), (A.4) where A n ∈ R >0 corresponds to the amplitude of the Fourier modes and γ n ∈ [0, 2π] the phase. As we remarked above, in this representation each A n is a random variable taking values in a Gaussian distribution, while each γ n is a random variable taking values uniformly in [0, 2π). However a useful simplification is to make A n deterministic and only keep the phases random. One can then check that taking A n =V σ(k n )∆k and averaging uniformly in [0, 2π), dγ n 2π (· · · ), (A.5) imply that for σ(k n ) = 1 in the thermodynamic limit E[f (x)f (y)] =V 2 δ(x − y). In other words, we reproduce the Gaussian behaviour with a simpler setup in which only the phases fluctuate. Similarly, we can obtain non-Gaussianity by choosing a non-trivial function σ(k n ).

B Holographic Renormalisation
In this Appendix we discuss the details of holographic renormalisation in our inhomogeneous geometry. Let M be the underlying manifold defined by the geometry in Eq.(2.1). The action for the probe scalar is given by where n is the normal unit vector pointing outwards of the boundary ∂M and γ is the respective induced metric. Note that strictly at the conformal boundary ρ = 0 Eq.(B.2) is divergent. To regularise this divergence, we evaluate the on-shell action at a slice ρ = λ 1 which we later take to zero. The normal unit vector pointing outward the fixed ρ = λ surface is then given by n = −2ρ ∂ ρ | ρ=λ , while the induced metric is and thus √ −γ = λ −1 −g (0) with −g (0) = 1 + g 2 tx . Note that g (0) is interpreted holographically as the metric where the dual boundary field theory lives. Inserting in the on-shell action, Note that the highest diverging is in the term ∝ s 2 . To remediate this particular divergence, we add the following boundary counter-term to the action

(B.7)
In most of the manuscript we work with ν = 1/2, for which other divergences are not present, and therefore the only counter term required is the above. However, for larger values of ν there is a finite tower of higher divergences between the leading term ∝ s(x) 2 and the term ∝ s(x)A(x) which we have omitted in the ellipsis, and which should also be taken into account in order to obtain a finite one-point function. For instance, the next order divergence would be of order O λ −(ν−1) , which can be remediated by adding a derivative term S (1) ct = 1 2ν−2 ∂M d 2 √ −γφ ∆ γ φ| ρ=λ . Higher order terms can be treated in the same way, adding higher derivative terms S (k) ct accordingly (for a complete discussion, see [72]). However note these derivative terms do not contribute to the coefficient of the one-point function ∝ s(x)A(x). Accounting for all the divergences, the renormalised action reads This justifies the terminology boundary-to-bulk propagator, since it propagates the source s(x µ ) living in the boundary into a scalar field satisfying the Bulk equations of motion. From this relation it is also possible to see that the bulk-to-bulk and boundary-to-bulk propagators are related to the boundary tree-level boundary two-point function. Since we have O(x µ ) = (2ν) lim ρ→0 ρ −∆ + /2 ψ(ρ, x µ ) (C. 9) we see that if we define the boundary tree-level two-point function between the source and the expected value.