Field Theory of Disordered Elastic Interfaces at 3-Loop Order: The $\beta$-Function

We calculate the effective action for disordered elastic manifolds in the ground state (equilibrium) up to 3-loop order. This yields the renormalization-group $\beta$-function to third order in $\epsilon=4-d$, in an expansion in the dimension $d$ around the upper critical dimension $d=4$. The calculations are performed using exact RG, and several other techniques, which allow us to resolve consistently the problems associated with the cusp of the renormalized disorder.


Introduction
Disordered systems are notoriously difficult to treat, since naive perturbation theory leads to absurd results, as exemplified by the phenomenon of dimensional reduction [1]. Two main paths out of this dilemma have been pursued: Replica symmetry breaking [2], and the functional renormalization group. The latter goes back to the work by Wilson [3] and Wegner and Houghton [4]. For disordered systems these methods were first used by Daniel Fisher [5]. However it took until 1992 that Narayan and Fisher [6,7], shortly thereafter followed by Natterman, Stepanow, Tang and Leschhorn [8], recognized that the disorder correlator, which plays the role of the coupling constant in the functional renormalization group (FRG) treatment, has to assume a cuspy form. The physical origin of this cusp lies in the metastability of the zero-temperature states which dominate the partition function, as recognized by Balents, Bouchaud and Mézard [9] in 1996. Only in 2006 this property was given a precise meaning as an observable, which can be measured in a numerical simulation both for the statics [10,11], the driven dynamics [12,13], and in an experiment [14]. This has nicely been reviewed in [15]. It was important conceptually, since the very existence of the cusp had in the early days questioned the validity of the method. Once this question of principle solved, it remained the problems of feasibility and practicality: First, whether there is a controlled loop or ε-expansion, and second how to implement a method which makes sense of the cusp in this loop expansion, and more particularly of the derivatives at the cusp. The latter change sign, depending on whether the limit is taken for positive or negative argument, not to mention the additional problems arising for a higher-dimensional field [16]. While these problems were conceptually simpler to solve for depinning [17], due to the Middleton-theorem [18], for the statics the question is more delicate. A consistent solution has been given at 2-loop order, based on renormalizability, recursive construction, or consistency schemes (the "sloop-algorithm" to be discussed below) [19,20], or exact RG [21][22][23]. At 3-loop order we performed several independent calculations. Here we give the resulting β-function.
The analysis of the fixed point will be published separately [24]. There we will extract the roughness exponent ζ , obtain the fixed-point functions R to 3-loop order, give the correction-toscaling exponent ω, as well as the momentum dependent 2-point function.
Here and below . . . denote thermal averages, and (. . . ) disorder ones. In the zerotemperature limit,the partition function is dominated by the ground state, and we may drop the explicit thermal averages. The random potential can without loss of generality [19,20] be chosen Gaussian with second cumulant (2.5) R 0 (u) takes various forms: Periodic systems are described by a periodic function R 0 (u), randombond disorder by a short-ranged function, and random-field disorder of variance σ by R(u) −σ |u| at large u.
To average over disorder, we replicate the partition function n times, Z n =: e −S , which defines the effective action S, (2.6) We used the notations introduced in Eqs. (2.2) and (2.3). In presence of external sources j a , the n-times replicated action becomes from which all static observables can be obtained. a runs from 1 to n, and the limit of zero replicas n = 0 is implicit everywhere.
− π 2 9 = 0.5859768096723648... (3.4) = 1.4140231903276352... . (3.5) It has some limitations though. Indeed, if one applies this procedure to the 3-loop calculation, one finds that it renders unique all but one ambiguous diagram, namely , (4.1) which has no subdivergence. Thus there are no counter-terms which could lift the ambiguities. This diagram must therefore be computed directly and we found that it can be obtained unambiguously by the sloop elimination method. We give an explanation of this method in section 6.1; it is also well documented, see section VD of [20].

5) Reparametrization invariance:
From standard field theory, one knows that RG functions are not unique, but depend on the renormalization scheme. Only critical exponents are unique. This is reflected in the freedom to reparametrize the coupling constant g according to g −→g(g) where g(g) is a smooth function, which has to be invertible in the domain of validity of the RG β-function.
Here we have chosen a scheme, namely defining R(u) from the exact zero momentum effective action, using dimensional regularization, and a mass. One can explore the freedom in performing reparametrizations. In the functional RG framework, reparametrizations are also functional, of the form where B(R, R) is a functional of R. For consistency, one has to demand that B(R, R) has the same analyticity properties as R, at least at the fixed point R =R * , i.e. B(R, R) should as R be twice differentiable and then have a cusp. Details can be found in Section 7.
As far as applicable, all methods, who are genuinely different, give consistent results. This is strong evidence that the problem has a unique field theory, which we identify in this paper to 3-loop order. In particular, the ambiguities which arise in perturbation theory due to the cusp turn out to be superficial and are absent in our treatment. Let us now turn to actual calculations using these methods. We start with the ERG approach. We will then use renormalized field theory and a combination of the above-mentioned methods. Let us stress that each of these two calculations was done independently by one of the authors, which serves as a non-trivial check of the RG β-function such obtained.

Calculation via the exact renormalization group
In this section we derive the 3-loop flow equations by means of the exact renormalization group (ERG). In condensed matter this RG is sometimes called "functional RG" because it is based on exact flow equations formulated for thermodynamic functionals. To avoid possible confusions, we will use the term "functional RG" only in the sense of perturbative field theory, i.e. as a loop expansion.

Set-up of ERG equation
For each realization of the random potential V , let Z V be the partition function. By the standard replica trick we average the logarithm of Z V over disorder by introducing replicas of the field. The replicated partition function is written as a functional integral It depends on an external replica-dependent field j a (x) with a = 1 . . . n. We choose N 0 = (Z n V ) −1 such that e W [0] = 1. We denote (j, u) = a d d x j a (x)u a (x) and the replicated action is given by Correlation functions and other observables averaged over disorder can be calculated from replica averages obtained from a polynomial expansion of W [j ], see Ref. [23] for details. For example, the connected 2-point correlation function is given by The mass m 2 > 0 provides an infrared regularization, and we are interested in the limit of m 2 → 0. The ERG is set up by successively lowering the parameter m 2 , which is our RG scale. Since the action S[u] depends on m only via its quadratic part in the fields, the scale derivative of W [j ] can be expressed by a Polchinski-type equatioṅ Here q ab (p) = T −1 (p 2 + m 2 )δ ab denotes the bare inverse propagator and the derivative with respect to the scale m is denoted by a dot. Due to momentum conservation the inverse propagator in Fourier space has only diagonal elements. The inner product between two fields u and v sums over the replica index and the spatial dependence We use matrix notation such that, for example, .
The second term in S [u] is invariant under a shift with a replica-independent field. This is expressed by the so-called statistical tilt symmetry (STS) where j is a replica-independent field and g ab (q) = q ab (q) −1 . It follows that the thermal propagator lim n→0 a

Replica expansion hierarchy
We expand the effective action in the number of replica sums n!T n a 1 ,...,a n S (n) [u a 1 , . . . , u a n ] , (5.14) where we use the shorthand notation u ab (x) = u a (x) − u b (x). Due to STS the one-replica term is given by the bare inverse thermal propagator. The two-replica term is a scale-dependent functional that depends on u ab (x) only. The initial condition for R is local and given by the bare disorder distribution function Higher replica terms are not present in the bare action but they are generated by the RG flow. STS implies that S (n) [u a 1 , . . . , u a n ] = S (n) [u a 1 + v, . . . , u a n + v] (5.16) for any field v(x). It follows that the two-replica term S (2) [u a , u b ] = R[u ab ] is a functional of only one field. Because of the sum over all replica indices, we assume all n-replica terms or -cumulants to be symmetric under permutation of the replica fields. We use the following compressed notation for functional derivatives of n-replica terms to denote p 1 derivatives with respect to field u a 1 and similarly p i derivatives with respect to field u a i for i = 1, ..., n with the total number of derivatives p max = n i=1 p i and the short-hand notation S (n) [u a 1 ...a n ] = S (n) [u a 1 , . . . , u a n ]. For example, using permutation symmetry, the second functional derivative of is given by . . , u a n ](x, y)δ ab = n{S 20...0 [u a 1 , . . . , u a n ](x, y)} + 2n(n − 1){S 110...0 [u a 1 , . . . , u a n ](x, y)} .
Because we are interested in the limit of the number of replica indices n → 0, we are free to add any function that depends on k < n replicas to a n-th cumulant. This "gauge invariance" will be used later to get rid of constant terms in the cumulants. Via Legendre transformation there is a one-to-one correspondence of -cumulants to cumulants obtained from a replica expansion of W [j ], see Ref. [23]. Therefore, the -cumulants can be used to calculate observables. In particular, the exact 2-point correlation function averaged over disorder is given by u=0 . (5.20) Here, compared to leading-order perturbation theory, the second derivative of the bare function R 0 (u) is replaced by the second derivative of the renormalized functional R [u].
In order to obtain RG equations for each -cumulant, we expand the inverse in Eq. (5.12) in a geometric serieṡ insert Eq. (5.18), and count the number of replica sums. The propagators g and gqg = −ġ are diagonal in replica space. Replica sums arise from second derivatives of ˆ , their matrix products, and one additional sum from the trace. Therefore, in order to calculate the flow equation of the n-th cumulant, the geometric series in Eq. (5.21) does not contribute for l > n. On the other hand, a term in the geometric series of l-th order contributes to cumulants to all orders n ≥ l. That is, for any initial action the RG flow generates cumulants to all orders. The term l = 0 and the one-replica term in l = 1 in Eq. (5.21) are constants and can therefore be neglected due to gauge invariance. Evaluating the two-replica contributions in the l = 1 and l = 2 terms give the flow equatioṅ The evaluation at zero field arises in terms of coinciding replica indices. We note that Eq. (5.22) is a non-linear integro-differential equation for a functional. Similar equations can be obtained for higher cumulants, Ṡ (3) and Ṡ (4) ; they are given in Appendix B.1. Due to the l = 1 term in Eq. (5.21) there is a contribution from S (m+1) to Ṡ (m) . Therefore, in order to obtain exact solutions for the -cumulants, one has to consider the full, infinitely large hierarchy. Note that, formally, up to now no approximations were made; in particular, we do not encounter ambiguities when a cusp in the second derivative of the local disorder distribution function develops.

ε-expansion for T = 0
Since we cannot treat an infinite hierarchy, we perform an additional expansion in ε = 4 − d.
To this aim we split the disorder distribution functional R[u] into a local and a non-local part If u(x) = u 0 is a constant field, then R [u 0 ] = 0, so that only the local part contributes, are functionals, whereas R(u) is a function. For m → ∞ the disorder-distribution function has only a local part, which we assume to be small. That is, R 0 (u) and all its derivatives are uniformly bound by a small constant. 1 We also assume that the local part R(u) of the renormalized disorder distribution function remains small. Then the ε-expansion can be obtained by expanding the replica expansion hierarchy in R(u). From now on we set the temperature to T = 0. Because the rescaled temperature T becomes small for small m, the ε-expansion can also be obtained for T > 0 by a composite expansion in T and R(u). Suppose that R(u) ∼ O(ε). Since for T = 0, Eq. (5.22) is quadratic in R , the non-local part of the renormalized disorder distribution function will be R [u] ∼ O(ε 2 ). A similar argument for the higher -cumulants gives S (n) ∼ O(ε n ) for n ≥ 3. The assumption that also all derivatives of the cumulants remain of the same order has to be checked; we will do so up to order ε 4 , that is, 3-loop order. With this method the 2-loop order was already obtained in Ref. [23].
The 1-loop equation can be obtained by an expansion to second order in ε and can directly be read off from Eq. (5.22). We use that (5.24) where the second term is already O(ε 2 ) and does not contribute to Eq. (5.22) at 1-loop order. We therefore finḋ where R (u) := R (u) − R (0). The local part is obtained by inserting the constant field u(x) = u 0 and dividing by L ḋ where after Fourier transformation I 1 = p g(p) 2 ∼ O( m −ε ε ), and so İ 1 ∼ O(1). The diagram I 1 is evaluated in Eq. (A.6). In order to have the simplest possible formulas, we will absorb a factor of εI 1 | m=1 into the renormalized disorder, see Eq. (6.43). This effectively sets I 1 to m −ε /ε. For an n-loop integral I n we will have to evaluate the ratio I n /I n 1 . We believe this to be the most convenient convention for obtaining standardized expressions.
Up to rescaling Eq. (5.26) is the standard 1-loop FRG equation. The solution of this flow equation corrects R 0 (u) ∼ O(ε) to the renormalized R(u) to order ε 2 . The non-local part in terms of this renormalized disorder-distribution function is given bŷ Superficially, the ε-expansion seems to work. However, we assumed that R [u] is of order ε 2 and likewise all derivatives of R [u]. In fact, the existence of the cusp in R (u) of the 1-loop solution appearing at a finite scale destroys our assumptions. Due to this cusp, R (u) is not differentiable at u = 0. The left-and right-sided limits of R exist but do not coincide R (0 The fourth derivative R (u) is uniformly bound for u = 0 but it is infinity at u = 0. If we would only need up to two derivatives of R, the ε-expansion would work without caveats. However, even the computation of the 2-point correlation function in Eq. (5.20) requires a second derivative of the non-local part of R [u], that is, a second derivative of Eq. (5.27) at zero field. There enters a third and fourth derivative of R(u(x)) that have to be evaluated at u(x) = 0. Furthermore, in the derivation of higher orders in ε, that is, higher orders in the expansion of the replica hierarchy in R(u), one encounters higher derivatives as well.
For the calculation of observables via analytic continuation it suffices to evaluate derivatives of R(u) at u = 0 ± , if the left-and right-sided limits give the same result for the observable.
For example, the ambiguity is avoided if odd derivatives of R(u) enter the equations only squared. However, we work with fluctuating fields u(x) and, for example, the second derivative of Eq. (5.27) contains a term R (u(x 1 ))R (u(x 2 )). While this is a square term we have to ensure that either u(x) → 0 + or u(x) → 0 − uniformly for all x. That is, ambiguities can be avoided if we restrict to non-crossing configurations.
Note that none of our methods can handle observables involving crossing configurations. However, handling those is not necessary for the present purpose. The calculation of correlation functions requires only configurations which are infinitesimally close to a uniform one (see e.g. the discussion in section V.E of [23]). Hence the two methods presented here are consistent, and consistent with each other. From now on, the limit of two fields u a (x) and u b (x) being equal in a -cumulant is understood as where u 0 is a constant field. That is, all fields are assumed to be close to a uniform configuration. We demonstrate in the next two subsections that in this weak limit the 3-loop β-function can be derived consistently. That is, it does not matter if the right limit u 0 → 0 + or left limit u 0 → 0 − are taken in Eq. (5.28).

ε-expansion to 2-loop order
As an instructive example we review the 2-loop ERG calculation done in Refs. [21,23]. In order to obtain Eq. (5.22) to order ε 3 we have to compute S (3) 1 ) to this order. We first concentrate on the second term. Note that we expand in the renormalized disorder distribution function R(u). This gives is of order ε and we have to expand the non-local part R to second order in R(u), that is, we have to insert the second derivative of Eq. (5.27) As described above, in order to calculate R [0 + ](z 1 , z 2 ) we first insert a constant field u 0 and then take the limit u 0 → 0 + to obtain It is important to note that in this "weak limit" we obtain no ambiguities since R (0 + ) 2 = R (0 − ) 2 has a straightforward analytic continuation. The feedback from S (3) is calculated by retaining only terms of order ε 3 . The calculation is relegated to appendix B.1. The result from Eq. (B.1) is To this order it integrates to S (3) [u abc ] = 3 tr gR ab gR ab gR ac − tr gR ab gR bc gR ac + O(ε 4 ) , Here the trace is over real space. The symmetrization over fields {. . .} can be written as S (3) . In a next step we take the functional derivatives of these terms with respect to u a (x) and u b (y). Then the limit b → a is performed by first replacing u b (y) by u a (x) + u 0 and then sending u 0 → 0 + . The remaining fields u a (x) only occur in the combination u a (x) − u c (x) and can directly be set to zero. When taking the functional derivatives of Eqs. (5.34), it is helpful to remember that we set b = a afterwards. Therefore A 2 does not contribute to S (3) 110 [0, 0, −v] and A 1 contributes only if the derivatives act on the first two R ab in the trace. Finally, the term A 3 is a symmetric functional in u ac and u bc and can be symmetrically expanded in u ab as outlined in the appendices of Refs. [23,41]. Setting u c (x) = u(x), we obtain In the above equation and from now on we use the shorthand notation R (x) := R (u(x)) (and similar for higher derivatives expect for x = 0 + and x = u). Inserting Eqs.
The 2-loop β-function known from FRG calculations [19,20] is the local contribution and is obtained by inserting a constant field and dividing by L d where The O( 1 ε ) term in İ A is canceled by I 1İ1 , ensuring a finite β-function. The non-local part integrates tô and

ε-expansion to 3-loop order
In 3-loop order we have to compute S to order ε 4 . The flow equation for the three-replica cumulant at T = 0, see Eq. (B.1), is given bẏ Because there is a feeding term from the fourth -cumulant S (4) we have to calculate S (4) to order ε 4 . The only contribution in Eq. (B.2) is S (4) 4 that integrates in this order to S (4) [u abcd ] = 3 4 tr gR ab gR ac gR ad gR ad + 2 tr gR ab gR ac gR cd gR ac (5.43) − 4 tr gR ab gR ac gR cd gR ad + tr gR ab gR bc gR cd gR ad where again R ab (x, y) = R [u ab ](x, y). In order to obtain S (4) 1100 [u aabc ] the equation has to be symmetrized over replica fields and two functional derivatives have to be taken. This lengthy but straightforward calculation is not reproduced here. The limit of identical replica fields in the first and second entry again has to be taken in the weak limit.
For brevity we introduce the symbol d dm g that formally denotes a scale derivative that acts only on the propagators g that were differentiated in the initial 1PI flow equations, see Eqs.
where A 1 , A 2 , and A 3 are given in Eq. (5.34). This term was easily integrated in 2-loop order since a scale derivative acting on R gives an additional order of ε, that is, Here we also need the next order, so we have to calculate and reinsert the 1-loop result for Ṙ [u] from Eq. (5.25) to obtain this expression to order ε 4 . It is sufficient to insert the local part of R into A 1 , A 2 , and A 3 . Apart from the feeding from S (4) 1100 [u aabc ], there are two more 3-loop contributions to the flow of S (3) . One arises by inserting also non-local contributions to 1-loop order from Eq. (5.27) into A 1 , A 2 , and A 3 . And, finally, there is a cross term R × S (3) from the third line of Eq. (5.42). Here we can insert the 2-loop solution S (3) with local R's into the right-hand-side of the flow equation to obtain the complete result at 3-loop order. These 3-loop contributions are easily integrated since scale derivatives acting on cumulants would introduce additional loops. The details of this calculation and the resulting functional S (3) [u abc ] to 3-loop order are given in Appendix B.2. In order to obtain S (3) 110 [0, 0, −u] it is again convenient to use a symmetric replica expansion. Setting u a = u b again requires the weak limit; in addition to potentially problematic terms ∼ R (0 + ) 2 we also encounter R (0 + )R (5) where R is the non-local contribution of the 2-loop solution Eq. (5.39). Taking two functional derivatives and taking the weak limit with non-crossing configurations produces anomalous terms 110 [0, 0, −u] in Eq. (5.22) allows us to rearrange terms such that they are total derivatives acting on the propagators g only. In summary we obtain to 3-loop ordeṙ given by Eq. (5.36) and the 3-loop contribution Here we once again introduced shorthand notations g xy := g(x, y) and R x := R(u(x)) except for x = 0 + and likewise for derivatives of R. Inserting a constant field and dividing by L d gives the 3-loop contribution to the β-function with the following integrals which are calculated in appendix A. It turns out that the combinations occurring in Eq. (5.49) are finite for ε → 0, so our counting of orders in ε is consistent, and the theory is 3-loop renormalizable. Due to gauge invariance we can add any scale-dependent function to R(u) that does not depend on the fields. In this way we can drop all constants from the β-function. The constants in Eq. (5.49) arise directly from Eq. (5.48). In the derivation of the latter we neglected gauge terms in S (3) and S (4) , so these constants are arbitrary. Assuming non-crossing configurations, that is, using Eq. (5.28) for taking limits u a − u b → 0, allows us to derive all anomalous terms in the β-function without ambiguities. With this assumption, the ε-expansion is a straightforward expansion of the exact hierarchy of flow equations for the -cumulants in powers of the effective local disorder distribution function R(u). Presumably, this will work to all orders in ε.
Crossing configurations could not be treated and are an open problem. It is doubtful that the standard ε-expansion can be applied. This is because R(u) and all its derivatives are not a small parameter suitable for an expansion if u = 0 cannot be avoided.
In order to make contact with the result obtained by an alternative method later in Section 6, we rescale where ζ is the roughness exponent. The rescaled function R still depends on the RG scale m and satisfies the RG equation to 3-loop order given in Eq. (3.1).

2-point correlation function
The expressions obtained above allow us to extract the non-local part, relevant to evaluate the full 2-point correlation function [24]. The latter is obtained as follows: The Fourier transform of Eq. (5.20) reads (5.52) Because of the limit n → 0 in Eq. (5.20), where n is the number of replica fields, this is an exact expression; in particular, there are no contributions of three-or higher-replica terms to the 2-point function. As in Sec. 5, the expression R [0 + ](z, z ) denotes the second functional derivative of R[u] with respect to u(z) and u(z ) that is evaluated in a weak limit u(x) ≈ const. → 0. A precise definition of the weak limit is given in Eq. (5.28). For the expansion in the renormalized local disorder function R(u) with a constant field u we use that where the non-local part R [0 + ](z, z ) has to be expanded to 2-loop order, like in the derivation of the 3-loop β-function. Taking two derivatives of Eq. (5.39), evaluated at a constant field u, and sending u → 0 + or u → 0 − we find (5.54)

Calculation using the sloop elimination method
Here we discuss a different way to do the contractions, using "excluded replicas", which will finally lead to a rather efficient algorithm for calculating the anomalous terms.
We start by a 1-loop diagram involving two disorder vertices, after having done one Wickcontraction. For simplicity of notation we are not writing space-indices and momentum integrals, which are unimportant for the following discussion.
At the next step, the following contractions are possible (restoring the integral) The second term is a 3-replica contribution (contribution to the third cumulant of the disorder), thus not of interest to us. The correction to the disorder at 1-loop order therefore consists of the first and last term, equivalent to the first and last two diagrams, This is equivalent to the result obtained in Eq. (5.26). An alternative approach consists in remarking that in Eq. (6.1) the terms a = b, and a = c could be dropped, since they are constants, thus will not be contracted in the next step. We thus start from which after one Wick-contraction leads to (6.5) The 2-replica term (the double sum) is, as expected, the same result as obtained in Eq. (6.2). While the second line contains only excluded replica sums, there can not be any ambiguity. The latter may only appear in the ensuing projection onto non-excluded replica sums. This is indeed the case for the hat diagram ∼ R (u)R (0 + ) 2 , as the reader is invited to check on his own, starting from We will therefore in the following present an improved projection method, which we have termed the "sloop-elimination" method. (The name may be thought of as "super"-partner of a normal loop, thus sloop, which cancels part of it.) The idea of the method is very simple. Let us consider the second term on the second line of Eq. (6.2). It is a three-replica term proportional to the temperature. In a T = 0 theory such a diagram should not appear, thus it can identically be set to zero: Projecting such terms to zero at any stage of further contractions is very natural in our present calculation (and also e.g. in the exact RG approach, where terms are constructed recursively and such forbidden terms must be projected out). It is valid only when (i) the summations over replicas are free (ii) the term inside the sum is non-ambiguous. These conditions are met for any diagram with sloops, provided the vertices have at most two derivatives. (One can in fact start from vertices which either have no derivative or exactly two.) Subtracting this term from Eq. (6.5) immediately yields the result (6.3). While this could also have been done directly, let us illustrate the power of the procedure on an example. We want to contract the expression (6.7) with a third vertex R where we have already dropped constant terms which will disappear after the contractions. Also note that implicitly here and in the following the vertices are at points x, y, z in that order. We will contract the third vertex twice, once with the first and once with the second, i.e. look at the term proportional to I A = x,y,z g(x − y) 2 g(x − z)g(y − z).
Performing the first contraction between points x and z yields Similarly, the second contraction yields (with the standard combinatorial factor of 1/2) This non-trivial identity tells us that the sum of all the terms (or diagrams) thus generated upon contractions must vanish. Stated differently: A sloop, as (6.7) as well as the sum of all its descendents vanishes. Note that this is not true for each single term, but only for the sum. A property that we request from a proper p-replica term is that upon one self contraction it gives a (p − 1)-replica term. It may also give T times a p-replica term (a sloop) but this is zero at T = 0, so we can continue to contract. Thus we have generated several non-trivial projection identities. The starting one is that the 2-replica part of (6.7) is zero, since (6.7) is a proper 3-replica term. This is what is meant by the symbol "≡" above and the last identity is the one we now use.
Indeed compare (6.10) with (6.6). One notices that all terms apart from the first in (6.6) appear in (6.10). They also have the same relative coefficients, apart from the third one of (6.6). Thus one can use (6.10) to simplify (6.6): (6.11) The function R (u) 2 , which appears in the last term, is continuous at u = 0. It is thus obvious how to rewrite this expression using free summations and extract the 2-replica part This coincides with the contribution of diagram A in the ERG approach, see the second term of Eq. (5.37). We can write diagrammatically the subtraction that has been performed as δ (2) A R = − , (6.13) where the loop with the dashed line represents the sub-diagram with the sloop, i.e. the term (6.10) (with in fact the same global coefficient). The idea is that subtracting sloops is allowed since they vanish. The advantage of the method is that all intermediate results are uniquely defined.
There are other possible identities, which are descendants of other sloops. For instance a triangular sloop gives, by a similar calculation: (6.14) This however does not prove useful to simplify δ (2) A R. Remains to calculate the 3-loop diagrams, shown on Fig. 1. This is achieved in appendix C. Since the above method generates a large number of identities, one can wonder whether they are all compatible. We have checked that this is indeed so, but we have not attempted a general proof.

The effective action up to 3-loop order
Using the sloop elimination method exposed in the preceding section, we have calculated all diagrams up to 3-loop order. They are presented graphically on Fig. 1, and given below. The expressions intervening in the sloop-projection algorithm are collected in appendix C. Here we give the final result for the effective action, before discussing how to obtain the β-function in the next section.
The 2-loop term is The 3-loop terms read (6.28)

Derivation of the RG-equation to 3-loop order
Let us now discuss in general the strategy to renormalize theories, whose interaction is not a single coupling-constant, but a whole function, here the disorder-correlator R(u). We denote by R 0 the bare disorder -this is the object in which perturbation theory is carried out -and by R the renormalized disorder, i.e. the corresponding term in the effective action .
We define the dimensionless bilinear 1-loop, trilinear 2-loop and quadrilinear 3-loop functions where if all arguments are the same, we only give this one argument, e.g. δ (1)

R, R) and δ (3) (R) = δ (3) (R, R, R, R). For different arguments we use the multilinear formulas
f (x, y) : Schematically, the renormalized disorder is calculated in the preceding section (where we had not explicitly written an index 0 to indicate the bare disorder). The inversion of relation (6.35) is Since an n-loop integral scales like m −nε the β-function is directly read off from (6.35), However, we need the β-function in terms of R, for which we replace R 0 by R, using Eq. (6.36), Using the results from Eqs. (6.16), (6.17) and (6.18), this is, printing one diagram and its counterterms (as dictated by the renormalization group R-operation) per line: On this form, one can explicitly check renormalizability. Since we kept the amplitudes of subdivergences, as for instance that of the 2-loop bubble-chain diagram, one can exactly see, where these terms come from. Actually the form given above is unique, even though several diagrams have the same functional dependence on R.
Let us now proceed to simplify the above equation. In order to do so, we have to choose a renormalization-scheme. We calculate the 3 leading terms in the ε-expansion of each diagram, i.e. up to order 1/ε for the 3-loop diagrams, up to order ε 0 for the 2-loop diagrams and up to order ε for the 1-loop diagram. In order to have the final result as simple as possible, we absorb a factor of εI 1 into R. This means that an n-loop integral has to be normalized by (εI 1 ) n . It is with this normalization that the amplitudes are given in appendix A. The advantage of this procedure is that integrals take the most simple form, and there are no spurious terms like ψ(1) or ζ (2). By this way, the 1-loop diagram is automatically subtracted completely and one never has to worry about its finite parts. However, we have a choice of how to subtract diagrams at 2-loop order. The most common choice is to subtract the divergent part only. The advantage of this procedure is that the 2-loop β-function takes the simplest form, with the combination of ε(2I A − I 2 1 ) in the second line of (6.39) replaced by 1 2 . The disadvantage is that then diagrams like (q) do not vanish, but have an amplitude proportional to (see last line of (6.39)) I q − I 1 I A (since I B = I 2 1 , and in our normalizations this is exact in any subtraction scheme). Now if at second order, we only subtract the diverging part of I A this combination becomes We therefore chose to always subtract the diagram exactly. At order 3 at which we are working here, this means that we have to keep the finite part of I A . This is sufficient, since the 1-loop integral is normalized to have no finite part, and since from the 3-loop integrals one only needs the diverging part anyway. Let us now use that to restate the β-function: and group together alike terms. This yields our final expression for the 3-loop β-function given in Eq. (3.1). The coefficients C 1 to C 4 , already given in Eqs. (3.2)-(3.5) are constructed from the diagrams via = 0.6085542725335131 (6.45) These constants are closely related to each other analytically.

Reparametrization invariance
It is known in standard field theory, that one can perform a change of variables, and thus formally change the β-function, while all observables remain unchanged. In the context of a functional RG, this reparametrization invariance is much larger. The function R(u) can be changed into an arbitrary functional of f [R]. The most useful such reparametrizations involve functionals f [R], which have the same structure as corrections to R, obtained perturbatively. Especially, when the field u has dimension ζ , and R times the 1-loop integral has dimension −4ζ , this means that on dimensional grounds for each additional power of R in f [R], there should be 4 derivatives. Also we do not want R(u) to have different analyticity properties, i.e. if R(u) has a r.h.s. Taylor-expansion with a missing linear term (absence of a super-cusp) then f [R] should have the same properties. The most suggestive such functional is the 1-loop contribution itself, which we study now.
The 2-loop RG-equation for the renamed disorder correlator R u reads Consider the following change of variables This is equivalent to stating that Solving this equation perturbatively yields the β-function for R u This equations tells us nothing more than that adding a coefficient of order ε to the secondorder term does not change universal results at 2-loop order. (The reader may want to verify this surprising result for the slope of the β-function at 2-loop order in a scalar field-theory.) Suppose now that β[R](u) = 0. Then this also holds for its derivatives and multiples thereof. Therefore, we can add to the fixed-point equation of the β-function terms of the form Note that these are the same terms, which appeared in equation (7.4).
In the following, we chose ζ = 0, since this yields the simplest relations. We will comment on the more general case later. Expression (7.6) then reads Adding −1/2 times (7.7) to the β-function (7.5) and choosing there λ = −1/2 to eliminate the additional 1-loop order term gives In this equation, we have traded the term proportional to R R 2 for a term of the form R R 2 .
Since the latter is uniquely defined, this allows again to fix the anomalous terms associated to R R 2 . It would be satisfactory, to have a similar result for the case ζ = 0. The above construction however yields terms of the form ζ uR u R u (7.9) plus the respective anomalous terms. Although one can of course solve differential equations involving these terms, and thus e.g. check the numerical solution of the fixed point equation to be discussed later, we have found no way to eliminate these terms, without generating even more "unusual" ones. Our search comprised rescalings of R u , of the field u, adding uβ[R] (u) to both the variable transformation and the β-function itself, and adding multiples of the β-function.
On the other hand, one can first write the β-function without rescaling, then do the non-trivial transformations given above, and finally perform the rescaling. This will simply give the standard rescaling terms.
Let us also comment on the power of reparametrization invariance at 3-loop order. While it proves to be a powerful tool for many diagrams, it is at least not applicable to fix all anomalous terms. This can be anticipated from the difference between diagrams (o) and (q) (see appendices C.3.8 and C.3.10), which is proportional to While (o) and (q) have the same normal terms, their difference is proportional to R (0 + ) 2 , thus the anomalous terms are different.

Conclusions
In this article, we have obtained the functional renormalization-group flow equations for the equilibrium properties of elastic manifolds in quenched disorder up to 3-loop order. The analysis of these findings will be given in a separate publication [24]: There we will extract the roughness exponent ζ , obtain the fixed-point functions R to 3-loop order, and give the correction-to-scaling exponent ω.
An interesting question is how the formalism derived here can be extended to N > 1 components. It had been shown in Ref. [16] that there is an ambiguity in the 2-point function already at 1-loop order. While this allowed the authors of [16] to still conclude on the β-function at 2-loop order, the problem becomes more severe at 3-loop order, and despite considerable efforts in this direction we have not been able to lift the ambiguities in some of the graphs.

A.1. General formulae, strategy of calculation, and conventions
We make use of the Schwinger parameterization In order to avoid cumbersome appearances of factors like 1 (4π) d/2 , we will write explicitly the last integral, and will only calculate ratios compared to the leading 1-loop diagram I 1 , given in the next section.
We will frequently use the decomposition trick which works well for dimension d ≤ 4. The reason for the utility of this decomposition is that it allows one to replace the massive propagator by a massless one, which is easier to integrate over, and a term converging faster for large k, which finally renders the integration finite. Special functions which appear are

A.2. The 1-loop integral I 1
The integral I 1 is defined as and is calculated as follows: We will also denote the dimensionless integral This gives us the normalization-constant for higher-loop calculations

A.3. 2-loop diagram I A
The non-trivial 2-loop integral can be written as The integrals J 1 and I 1 were solved by expanding the integrand in ε to order ε The result agrees with the one obtained by the subtraction method.

A.4. 2-loop integral I B
The trivial 2-loop diagram is A.6. I j I j = (A.20) Finally, All in all where two PolyGamma-identities were used B.1. Functional RG equations for S (3) and S (4) The flow equation of the third -cumulant in the ERG hierarchy is given bẏ We split the flow equation for the fourth -cumulanṫ

B.2. Third -cumulant S (3) to 3-loop order
In total there are four contributions to the flow of S (3) in 3-loop ordeṙ where the first contribution is known from the 2-loop calculation and reads where only the local part of R[v] is inserted, so u 1 is of order ε 3 . The second contribution comes from inserting the non-local part of R[v] to second order, that is Eq. (5.27), into 1 2 (A 1 +A 2 +A 3 ).
where we split the contributions according to different types of integrals. This is not the shortest way to write but better comprehensible. The same is done in the contributions from the RS (3) term where u 1 was used for S (3) on the right-hand side. Finally is the feeding term from S (4) , where we insert Eq. (5.43).

Appendix C. Systematic treatment of diagrams up to 3 loops: sloops and recursive construction
We present a systematic procedure to obtain (relatively) simple results for diagrams at up to 3 loops. The idea is to write the diagram, and then to consider all possible sloops which lead to the same diagram. Subtracting them with the right weight leads to results which are much simpler than those obtained by trying to reduce expressions term by term. The notation used throughout this section is We will write rather indistinguishably, in a little abuse of notation, R(u a − u b ) ≡ R u ≡ R ab , whatever is more convenient or suggestive. Below, we will give all diagrams.
There is always an additional combinatorial factor. At n-loop order, denote the number of propagators between points i and j as n i,j . Further denote the number of symmetries S as N S . Then the combinatorial factor for the contribution to R is at n-loop order, written apart from the diagram. We will give this factor at the beginning of each diagram with the same conventions as above.

C.1. 1 loop
Here we give the 1-loop diagram. A (closed) dashed line represents a sloop. Comb The final combination is Note that each sloop comes with a factor of (−1) and furthermore one has taken into account the proper combinatorial factor. This result is confirmed by the recursive-construction algorithm.

C.3.2. Diagram (i)
For a given order of the contractions, we have: = 16 6 g 4 ab + 16 g 3 ab g ac + 3 g 2 ab g 2 ac + 6 g 2 ab g ac g ad + g ab g ac g ad g ae +12 g 2 ab g ac g bc (C.32) = 16 g 2 ab g 2 ac + 2 g 2 ab g ac g ad + g ab g ac g ad g ae (C.33) = 16 4 g 3 ab g ac + 3 g 2 ab g ac g ad + g ab g ac g ad g ae + 3 g 2 ab g ac g bc . (C.34) The simplest combination is We can study another configuration, which we do not know how to draw, so call it S where we have already dropped the term a = c, which will disappear after the next contraction. Contracting (34) and then (24)  This is important since there is no counter-term in the theory for the divergence between the two leftmost vertices.

C.3.5. Diagram (l)
Diagram (l) has Comb = 1 2 3 × 1 × 1 2 . We use the notation (C.50) Diagram (l) is = 16 4f ab g 2 ab h ab + 8f ab g 2 ab h bc − 3f ab g ab g ac h bc + f bc g ab g ac h bc − f ab g 2 ac h bc +f bc g 2 ac h bc − 3f ab g ab g bc h bc − 5f bc g ab g bc h bc + f ab g 2 bc h bc + f bc g 2 ac h cd −f ac g ab g bc h cd − 2f bc g ab g bc h cd + 4f ac g ac g bc h cd + f bc g ab g bd h cd +f cd g ad g bd h cd − f bc g ac g cd h cd + f cd g ad g bd h de . (C.51) The 2-sloop is = 16 2f ab g 2 ab h bc − f ab g ab g ac h bc + f bc g ab g ac h bc − f ab g 2 ac h bc + f bc g 2 ac h bc −f ab g ab g bc h bc − 3f bc g ab g bc h bc + f ab g 2 bc h bc + f bc g 2 ac h cd − f ac g ab g bc h cd −2f bc g ab g bc h cd + 2f ac g ac g bc h cd + f bc g ab g bd h cd + f cd g ad g bd h cd −f bc g ac g cd h cd + f cd g ad g bd h de . (C.52) There is the special sloop configuration, which is obtained by starting from 16g ab (1)g ab (2) × h ac (2)R de (4). It is denoted and reads = 32 − f ab g ab g ac h bc − f ab g ab g bc h bc + f bc g ac g bc h bc + f ac g ac g bc h cd . There are of course much more possible sloops, involving three or four vertices. However, we did not use them here, and thus do not display them. The recursive-construction algorithm gives, consistent with the above It is not independent of the path of contractions. We number The simplest result is obtained by using contractions (12)(12)(34)(34)(13)(24) = 16 4 g 4 ab + 8 g 3 ab g ac + 8 g 2 ab g 2 ac + 8 g 2 ab g ac g ad + g ab g ac g ad g ae −4 g 2 ab g ac g bc − 2 g 2 ab g ac g bd − g ab g ad g bc g cd .
Another result is obtained using (12) We check that this projects to 0. Now the sloops give = 16 4 g 3 ab g ac + 2 g 2 ab g 2 ac + 6 g 2 ab g ac g ad + g ab g ac g ad g ae − 2 g 2 ab g ac g bc −2 g 2 ab g ac g bd − g ab g ad g bc g cd (C.60) = 16 4 g 2 ab g ac g ad + g ab g ac g ad g ae − 2 g 2 ab g ac g bd − g ab g ad g bc g cd . (C.61) The following combination is simple = 16(4f ab g ab 2 h ab + 4f ab g ab g ac h ab + 6f ab g bc 2 h ab + 4f bc g bc 2 h ab +3f bd g bc 2 h ab + 2f ab g bc g bd h ab + 2f bc g bc g bd h ab + f bd g bc g be h ab +3f ac g cd 2 h ab + f ac g cd g ce h ab + 2f bc g bc g bd h ac − 2f cd g bc g bd h ac −2f bc g ab g ac h bc ) .
(C.91) Sloop (13)(13), then (24)(24)(23)(34) gives = 16(4h 0 f ab g ab 2 + 4h 0 f ab g ab g ac − 2h 0 f bc g ab g ac + 2h 0 f ac g ab g ad +6h 0 f ab g bc 2 + 4f bc g bc 2 h ab + 3f bd g bc 2 h ab + 2f bc g bc g bd h ab +f bd g bc g be h ab + 3f ac g cd 2 h ab + f ac g cd g ce h ab + 2f bc g bc g bd h ac −2f cd g bc g bd h ac ) . (C.92) Sloop (24)(24), then (13)(13)(34)(24) yields = 16(4f ab g ab g ac h ab + 2f ab g bc g bd h ab + 2f bc g bc g bd h ab + f bd g bc g be h ab +f ad g cd 2 h ab + f ac g cd g ce h ab + 2f ac g bc 2 h ac + 2f bc g bc g bd h ac −2f cd g bc g bd h ac + f cd g bd 2 h ad − 2f bc g ab g ac h bc ) . (C.93) Sloops (13)(13) and (24)(24), then (23)(34) gives = 16(4h 0 f ab g ab g ac − 2h 0 f bc g ab g ac + 2h 0 f ac g ab g ad + 2h 0 f ac g bc 2 +2f bc g bc g bd h ab + f bd g bc g be h ab + f ad g cd 2 h ab + f ac g cd g ce h ab +2f bc g bc g bd h ac − 2f cd g bc g bd h ac + f cd g bd 2 h ad ) .