Newtonian Binding from Lattice Quantum Gravity
NNewtonian Binding from Lattice Quantum Gravity
Mingwei Dai and Jack Laiho
Department of Physics, Syracuse University, Syracuse, NY 13244
Marc Schiffer
Institut f¨ur Theoretische Physik, Universit¨at Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany
Judah Unmuth-Yockey ∗ Department of Physics, Syracuse University, Syracuse, NY 13244 andFermi National Accelerator Laboratory, Batavia, IL 60510 (Dated: February 10, 2021)We study scalar fields propagating on Euclidean dynamical triangulations (EDT). In this work westudy the interaction of two scalar particles, and we show that in the appropriate limit we recoveran interaction compatible with Newton’s gravitational potential in four dimensions. Working inthe quenched approximation, we calculate the binding energy of a two-particle bound state, and westudy its dependence on the constituent particle mass in the non-relativistic limit. We find a bindingenergy compatible with what one expects for the ground state energy by solving the Schr¨odingerequation for Newton’s potential. Agreement with this expectation is obtained in the infinite-volume,continuum limit of the lattice calculation, providing non-trivial evidence that EDT is in fact a theoryof gravity in four dimensions. Furthermore, this result allows us to determine the lattice spacingwithin an EDT calculation for the first time, and we find that the various lattice spacings are smallerthan the Planck length, suggesting that we can achieve a separation of scales and that there is noobstacle to taking a continuum limit. This lends further support to the asymptotic safety scenariofor gravity.
I. INTRODUCTION
Quantum gravity is one of the great outstanding prob-lems of theoretical physics. One possible approach toobtaining a consistent, predictive theory is the asymp-totic safety scenario of Weinberg [1]. It is well-knownthat general relativity is not perturbatively renormaliz-able [2, 3]. Although a low-energy effective theory canbe formulated [4, 5], it comes with an infinite numberof unknown couplings and is thus limited in terms of itspredictive power. Asymptotic safety rests on the hopethat there is a strongly coupled ultra-violet fixed point,making it effectively renormalizable nonperturbatively,with only a finite number of parameters needing to bespecified from experiment. One would also hope that thenumber of relevant parameters needing to be specifiedfrom experiment is small, since a theory with few inputparameters is more predictive than one with many. Ifasymptotic safety is realized for gravity, nonperturbativemethods will be crucial to convincingly establish this andto explore the detailed properties of the theory.Evidence that the asymptotic safety scenario is real-ized for gravity comes mainly from the lattice [6–9] andfrom the functional renormalization group [10–12], whichwas pioneered for gravity in [13] (see, e.g., [14–22] forintroductions and reviews). This paper follows up onprevious work [9] using Euclidean dynamical triangu-lations (EDT), which is one of the original approachesto quantum gravity on the lattice [23–25]. EDT is a ∗ [email protected] discretization of Euclidean quantum gravity, where thepath-integral is given by a sum over geometries, weightedby the Einstein-Hilbert action plus the action of the mat-ter sector. In a lattice formulation of an asymptoticallysafe theory, the ultraviolet fixed point would appear as acontinuous phase transition, the approach to which de-fines a continuum limit. This is the first test that EDTmust pass. In addition to this, the theory must repro-duce classical general relativity in the appropriate limit.If EDT could pass these two tests, it would be a can-didate for an ultraviolet-complete theory of four dimen-sional quantum gravity.We briefly review the evidence so far that EDT mayrealize the asymptotic safety scenario for gravity. Ref. [9]introduced a fine tuning of the exponent of a local mea-sure term, showing that the tuning of this parameter isnecessary in order to recover semiclassical physics. Al-though the local measure term was first introduced sometime ago in Ref. [26], the nice features of this tuninghave only been appreciated recently. Once the tuningis done, four-dimensional geometries are recovered thatresemble Euclidean de Sitter space. There also appearsto be no obstacle to taking the continuum limit by fol-lowing the first order line to a possible critical endpoint,and ensembles following this procedure were generatedat a number of different lattice spacings. The globalHausdorff dimension was measured using finite-volumescaling and shown to be close to four [9]. The spectraldimension, which is a fractal dimension defined by a dif-fusion process, varies with distance scale and approachesa value close to four at long distances. The variationof the spectral dimension with distance had also been a r X i v : . [ h e p - l a t ] F e b found earlier in other approaches [27–29]. The averageover geometries in Ref. [9] gives a result that is close tothat of Euclidean de Sitter space, and the quantitativeagreement with the classical solution gets better as theproposed continuum limit is approached. The agreementbetween the classical solution and the lattice data is ac-tually the worst at long distances, but improves as thelattice spacing is reduced. This might seem surprising,since it is typically the short-distance behavior that ismodified by discretization effects, but this type of effecton long-distance behavior is typical when a symmetry ofa theory is broken by the regulator, for example by thefinite lattice spacing in the case of lattice regularization.Then a fine-tuning is needed to approximately restore thesymmetry at finite lattice spacing. Residual symmetry-breaking effects appear at finite lattice spacing and canmodify the physics that depends on the symmetry. Thecorrect long-distance physics is then only recovered whenthe symmetry is restored in the continuum limit. An ex-ample of this is the Wilson fermion formulation of latticequantum chromodynamics, where the lattice regulatorbreaks chiral symmetry. There a fine-tuning is requiredto restore chiral symmetry, and even then, at finite lat-tice spacing the chiral symmetry breaking leads to dis-tortions of the pion sector, which contains the lightestphysical states of the theory. Ref. [9] argued by analogyto the Wilson fermion case that the symmetry that isbroken by dynamical triangulations is continuum diffeo-morphism invariance.Further evidence that dynamical triangulations recov-ers the correct long-distance limit was given in Ref. [30],where it was shown how to incorporate K¨ahler-Diracfermions [31]. This approach generalizes the staggeredfermion formulation to the random lattices of dynamicaltriangulations without the need to introduce vielbeinsor spin connections. It is well known that the K¨ahler-Dirac action reduces to four copies of Dirac fermionsin the flat-space, continuum limit [32]. Ref. [30] foundevidence that this is the case for K¨ahler-Dirac fermionson dynamical triangulations in the large-volume (smallcurvature) limit, but only if the continuum limit is alsotaken. This is seen in the approximate four-fold degen-eracy in the low-lying eigenvalues of the K¨ahler-Diracmatrix, and in the degeneracy of scalar bound states inthe continuum limit. The four-fold degeneracy is liftedby lattice effects in a similar manner to what is foundwith staggered fermions in lattice quantum chromody-namics (QCD) [33], but as in lattice QCD, the degener-acy appears to be restored in the continuum limit. Theevidence that these discretization effects vanish suggeststhat continuum K¨ahler-Dirac fermions in four dimensionsare recovered at that point. An additional advantage ofK¨ahler-Dirac fermions is that they possess an exact U (1)symmetry, which is related to continuum chiral symme-try. A study of fermion bilinear condensates providesstrong evidence that this U (1) symmetry is not sponta-neously broken at order the Planck scale, implying thatfermion bound states do not acquire unacceptably large masses due to chiral symmetry breaking. These resultsfor K¨ahler-Dirac fermions in EDT are highly non-trivialand provide further evidence for the asymptotic safetyscenario for gravity and matter.In this work we continue our investigation of matterfields living on dynamical triangulations, this time focus-ing on scalar fields. We follow the work of de Bakkerand Smit [34] in our implementation of scalar fields.They looked at scalars propagating on EDT backgrounds(without the inclusion of a local measure term) and thebinding energy of bound states of scalar particles. Theyfound clear evidence of gravitational binding, such thatthere was an attractive force between the scalars, butthey were not able to make contact with the Newto-nian limit (see [35] for recent work attempting to under-stand the systematic errors associated with the originalde Bakker, Smit analysis). In the original reference [34],these authors outlined a prescription for recovering theNewtonian limit. In this work we follow their prescrip-tion, using our approach to taking the continuum limitand our general strategy, which has been successful forother observables. We confirm the finding that there is anattractive force between scalar particles, and after tak-ing the continuum, infinite volume limit, we find resultscompatible with Newtonian gravity in the non-relativisticlimit. We can therefore use Newton’s law to infer a valueof Newton’s constant G , which allows us to convert lat-tice units into units of the Planck length. We find thatour lattice spacings are smaller than the Planck lengthand that for our finest lattice spacings we are starting tosee a separation of scales, such that the lattice spacing isstarting to become much smaller than the Planck scale.This provides evidence that not only does the formula-tion recover the correct long-distance physics, but alsothat there is no barrier to taking the continuum limit.This paper is organized as follows: Section II reviewsthe EDT formulation and the details of the lattice sim-ulations. Section III reviews how to add scalar fields toEDT and how to compute the scalar propagator and thebinding energy of a bound state of two scalar particles.Section IV presents our numerical results for the bindingenergy and the analysis needed to recover the form ofthe interaction and Newton’s constant. We conclude inSection V. Finally, an appendix gives some details of thefitter used in the analysis. II. EUCLIDEAN DYNAMICALTRIANGULATIONSA. The model
In the continuum, for fixed global topology, the four-dimensional Euclidean quantum gravity partition func-tion is given by the path integral sum over all geometries,weighted by the Einstein-Hilbert action and the actionfor the matter sector, Z E = (cid:90) D [ g ] D [ φ ] e − S EH [ g ] − S M [ φ ] , (1)where the Euclidean Einstein-Hilbert action is S EH = − πG (cid:90) d x √ g ( R − , (2)with R the curvature scalar, Λ the cosmological constant, G Newton’s constant, and we adopt for the matter sec-tor a real, non-interacting massive scalar field minimallycoupled to gravity, S M = (cid:90) d x √ g (cid:18) g µν ∂ µ φ∂ ν φ + 12 m φ (cid:19) , (3)with bare mass m for the scalar field.The partition function for lattice quantum gravity thatwe use in this work is that of Euclidean dynamical trian-gulations, where the path integral at finite lattice spacingis approximated by the sum over all four-geometries con-structed by gluing together four-dimensional, equilateralsimplices. It is given by [23, 36] Z E = (cid:88) T C T N (cid:89) j =1 O ( t j ) β e − S ER (4)where the factor C T divides out equivalent ways of label-ing the vertices in a given geometry, the term in bracketsis a local measure term with the product over all trian-gles, and O ( t j ) is the order of triangle j , i.e. the numberof four-simplices to which it belongs. S ER is the Einstein-Regge action [37] of discretized gravity, S ER = − κ N (cid:88) j =1 V δ j + λ N (cid:88) j =1 V , (5)where κ = (8 πG ) − , λ = κ Λ, δ j = 2 π − O ( t j )arccos(1 / t j , andwhere the volume of a d -simplex of equilateral edge length a is given by V d = √ d + 1 d ! √ d a d . (6)It is standard to absorb the overall numerical factorsinto constants and to perform the sums in Eq. (5) so thatthe lattice action is given by the concise form S ER = − κ N + κ N , (7)with N the number of four simplices and N the numberof triangles. The parameters κ , κ , and β are inputs tothe simulations whose values must be adjusted in orderto approach the continuum limit. We do not include thematter action in the Boltzmann weight when generatingour ensembles; this is known as the quenched approxima-tion. Although this is an uncontrolled approximation,we still expect many features of the full theory to bepreserved in the quenched theory. The main advantageof quenching is that it allows us to reuse the ensemblesgenerated and analyzed in previous works [9]. B. Simulation details
The details of the generation of the lattice configura-tions used in this work are given in Ref. [9], and they aresummarized briefly here. The lattice geometries are madeby gluing together four-dimensional simplices along theirthree-dimensional faces. The four-simplices are equilat-eral, with constant edge length a , and the dynamics isencoded in the connectivity of the simplices. We sumover a class of triangulations known as degenerate trian-gulations [38]. This class of triangulations does not obeythe combinatorial manifold constraints, though it doeslead to a reduction of finite-size effects by approximatelya factor of 10 over combinatorial triangulations, givingit some numerical advantages [38]. When the triangula-tions are degenerate, not all of the neighbors of a givenfour-simplex are necessarily distinct. Also, distinct four-simplices may share the same five vertex labels. Suchconfigurations are not allowed when the lattices obey thecombinatorial manifold constraints. Evidence for the ex-istence of a continuum limit for degenerate triangulationswith a non-trivial measure term was given in Ref. [9]. Itis likely that this universality class would be shared bycombinatorial triangulations if the continuum limit doesin fact exist, as discussed in Ref. [9].The algorithm for evaluating the partition function inEq. (4) is now standard, and consists of a set of ergodiclocal moves, known as the Pachner moves, which are usedto update the geometries [24, 39, 40]. The proposed localchanges are then accepted or rejected using a Metropolisstep. Most of the lattices used in this work were gener-ated using a parallel variant of this standard algorithmknown as parallel rejection [9]. Parallel rejection partiallycompensates for the low acceptance of the proposed lo-cal changes by the Metropolis accept/reject step. Wedefine a sweep as a fixed number of attempted moves,though the precise number varies across ensembles. Asweep could be 10 or 10 attempted moves, with theacceptance rate depending strongly on the location inthe phase diagram. Acceptance rates range from around10 − to 10 − . Measurements are made after 10 to 50sweeps, and our longer runs were for tens of thousands,or even one hundred thousand sweeps.The global topology of the geometries is fixed to S .The topology remains fixed since the local Pachner movesare topology preserving, and we choose S by startingfrom the minimal four-sphere at the beginning of theMonte Carlo evolution. In order to take the infinite lat-tice volume limit it is necessary to fine-tune the bareparameter κ to its critical value. This amounts to atuning of the bare cosmological constant. Although thistuning is required to take the infinite lattice volume limit,it does not guarantee that we can take the infinite phys-ical volume limit. The infinite lattice volume limit canbe taken even in an unphysical phase where the extentof the lattice “universe” is only a few lattice spacings, sothat the extent of the physical volume is on the order ofthe ultraviolet cut-off scale. As shown in Ref. [9], there isa region of the phase diagram with extended geometries,so that the large lattice-volume limit also corresponds tothe limit of large physical volume.The simulations are done at approximately fixed latticefour-volumes. The Pachner moves require the volume tovary in order for them to be ergotic, though in practicethe volume fluctuations about the target four-volume areconstrained to be small. This is done by introducing avolume preserving term in the action δλ | N f − N | , whichkeeps the four-volume close to a fiducial value N f . Thisterm does not alter the action when N = N f , but it doeskeep the volume fluctuations from becoming too large tobe practical in numerical simulations. In principal oneshould take the limit in which δλ goes to zero, but inpractice setting it small enough does not appear to leadto significant systematic errors. In Ref. [9] the parameter δλ was taken to be 0.04, though taking this value smallerdid not lead to a noticeable change in the results.The phase diagram for the theory is shown in Figure 1.The solid line AB is a first order phase transition linethat separates the branched polymer phase and the col-lapsed phase. Neither of these phases looks much likesemiclassical gravity. The branched polymer phase hasHausdorff dimension 2, while the collapsed phase has alarge, possibly infinite, dimension. The crinkled regiondoes not appear to be distinct from the collapsed phase,but is connected to it by an analytic cross-over. Thecrinkled region requires larger volumes to see the charac-teristic behavior of the collapsed phase, suggesting thatit is a part of the collapsed phase with especially largefinite-size effects [41]. The results of Ref. [42] for combi-natorial triangulations are in broad agreement with thispicture of the phase diagram.The exponent β associated with the local measure termmust be fine-tuned in order to obtain ensembles withfour-dimensional, semiclassical properties. It was shownin Ref. [9] that close to the transition line AB the lat-tice geometries have semiclassical properties, with globalHausdorff dimension close to four and a spectral dimen-sion that approaches four at long distances. The prescrip-tion for recovering semiclassical geometries is to approachthe line AB from the left. Since a first order transitionis characterized by tunneling between meta-stable states,if we simulate too close to the transition line then someof the time a given run will be in the wrong state. Inorder to minimize contamination from the wrong phase,we only take those parts of a run that are in the correctphase. We find that this selection procedure becomes un-necessary at the finest lattice spacing considered in thiswork, since we can simulate somewhat farther from thetransition line without losing the four-dimensional finite-volume scaling. This is fortunate, since future runs willbe at ever finer lattice spacings, so will likely not sufferfrom this systematic error.Table I shows the parameters for the ensembles usedin this work. Most of these ensembles were generatedand discussed in Ref. [9], but some of them are new. Wehave added a larger lattice volume at our finest lattice β BranchedPolymerPhaseCollapsedPhase κ BA DC Crinkled Region
FIG. 1. Schematic of the phase diagram as a function of κ and β . spacing, as well as larger and smaller lattice volumes at β = 0 to aid with the infinite-volume extrapolation.The determination of the relative lattice spacing fol-lows the procedure of Ref. [9], which looked at the returnprobability of a diffusion process on the lattice geome-tries. The return probability is dimensionless, but variesas a function of the diffusion time step, which is not.One can rescale the diffusion time step so that the re-turn probability lies on a universal curve; the rescalingfactor then tells us the relative lattice spacing. The rela-tive lattice spacings quoted here differ slightly from thosefirst reported in Ref. [9]. In that work the overlap of thereturn probability curves was obtained without regardfor the lattice volumes used in the comparison. Here wematch in a self-consistent way ensembles with latticesas close to the same physical volumes as possible. Theshifts are within previously quoted errors, but the cen-tral values quoted here are likely more accurate. Theerrors quoted here are similar in size to those previouslyobtained and reflect the uncertainties in matching thecurves in this procedure. III. THEORETICAL BACKGROUND
In this section we give the necessary theoretical back-ground for extracting and interpreting the numerical dataassociated with the binding energy of two scalar particles.We begin with a discussion of scalar fields on dynamicaltriangulations, then we show how to calculate the bind-ing energy between scalar particles on our lattices, andfinally we review the nonrelativistic limit of the bindingenergy in the continuum limit. (cid:96) rel β κ N Number of configurations1.59(10) 1.5 0.5886 4000 3671.28(9) 0.8 1.032 4000 5241 0.0 1.605 2000 2481 0.0 1.669 4000 5751 0.0 1.7024 8000 4891 0.0 1.7325 16000 5011 0.0 1.75665 32000 12180.80(4) − . − . − .
776 3.0 16000 2341TABLE I. The parameters of the ensembles used in this work.The first column shows the relative lattice spacing, with theensembles at β = 0 serving as the fiducial lattice spacing.The quoted error is a systematic error associated with finite-volume effects. The second column is the value of β , the thirdis the value of κ , the fourth is the number of four-simplices inthe simulation, and the fifth is the number of configurationssampled. A. Scalar fields on dynamical triangulations
In order to map Eq. (3) on to the lattice, we use anaive discretization where we associate each scalar fieldat x , φ ( x ), with a four-simplex (or equivalently a dual-lattice site). In this way each scalar field is restrictedto only have five neighbors through the five surroundingtetrahedra (or the dual edges). The lattice action on thedual lattice can then be written as S lat M = 12 (cid:88) (cid:104) xy (cid:105) ( φ x − φ y ) + m (cid:88) x φ x (8)where (cid:104) xy (cid:105) is a nearest-neighbor pair between two duallattice sites (or two four-simplices) which is counted once,and m is the bare mass. Note that the lattice scalaraction has a shift symmetry in the absence of a mass term[30, 43, 44]. This is also true of the continuum actionEq. (3) for zero mass and in the absence of any additionalscalar self interaction terms. This shift symmetry ensuresthat the renormalized mass goes to zero as the bare massapproaches zero without any fine tuning. It is a non-trivial check of the calculations that this is the case.Expanding and collecting terms above—as well as not-ing that we work on a compact topology—and normal-izing the coefficient of the kinetic term to one, we canwrite Eq. (8) as S lat M = (cid:88) x,y φ x L xy φ y (9)with L the dual lattice Laplacian plus a mass term, andthe sums over x and y are unrestricted over the wholelattice. In the continuum L is the Klein-Gordon operatorand the inverse of the L matrix is the scalar propagator.In the next subsection we review the form of the bindingenergy between scalar particles in the presence of gravityon the lattice. B. The binding energy
The binding energy between two particles is given by E b = 2 m − M , where m is the mass of a single particle,and M is the mass of the two-body bound state. If E b is greater than zero, it indicates that the two-body statehas a lower energy than twice the energy of the single-particle state, and thus a bound state can form.To compute the binding energy between two scalar par-ticles, we must first calculate the propagator for a freescalar field living on a triangulation. Given the latticescalar action of Eq. (9), we need only invert the matrix L xy to obtain the propagator. For an arbitrary triangu-lation, we can write down L xy explicitly, L xy = ( D x + m ) δ xy − A xy (10)where D x is the number of neighboring four-simplices toa four-simplex (in our case always five), m is the bareinput mass, δ xy is the Kronecker delta, and A xy is theadjacency matrix, which has matrix elements A xy = (cid:40) x and y share a dual edge0 otherwise. (11)The matrix elements of L − contain the propagators be-tween two four-simplices. To compute the binding en-ergy, we need the squares of the propagators as well,( L − xy ) . The interpretation of ( L − xy ) is the two-particlepropagator between x and y .With the one- and two-particle propagators, we cancompute the average two-point correlation function as afunction of geodesic distance on the lattice. The one- andtwo-particle, two-point correlation functions are, respec-tively [34], G ( r ) = (cid:42) (cid:80) x,y L − xy δ | x − y | ,r (cid:80) x,y δ | x − y | ,r (cid:43) (12)and G (2) ( r ) = (cid:42) (cid:80) x,y ( L − xy ) δ | x − y | ,r (cid:80) x,y δ | x − y | ,r (cid:43) (13)as a function of r , the geodesic distance between duallattice points. In Eqs. (12) and (13) there are two av-erages shown. The first is the average over all possibledistances a separation | x − y | = r away. In practice thisis not done in its entirety, and a fixed number of sourcepoints are chosen, from which propagators and distancesare calculated to all other points. The second averageis indicated with the angled brackets, which denote anensemble average—an average over configurations.From G and G (2) it is possible to compute the bind-ing energy between two scalar masses. We can use theasymptotic forms of G and G (2) , which are expected tofeature an exponential fall-off for large distances, G ( r ) ∝ e − mr r p , G (2) ∝ e − Mr r q , (14)to extract the binding energy. Above, m is the renormal-ized one-particle mass, and M is the two-particle renor-malized mass. Consider the function, F ( r ) = G (2) G , (15)which for sufficiently large source-sink separations be-comes F ( r ) ∝ e − ( M − m ) r r q − p . (16)We can identify the binding energy with the exponent, E b ≡ m − M, (17)and we can identify the power-law exponent, γ ≡ q − p .Thus E b measures the difference in energy between a two-body system in its ground state and twice the mass ofa single particle in its ground state, with gravitationaleffects taken into account.With these definitions, taking the logarithm of F , wefind log F ( r ) (cid:39) E b r + Z E − γ log r, (18)where Z E is a constant. The single particle propagatorhas the same form after taking the logarithm,log G ( r ) (cid:39) − mr + Z G − p log r. (19)Besides the constant shifts present in both equations, therenormalized mass and the binding energy both appearas linear terms, while the power-law behavior appears asa logarithmic correction. In the physically relevant case,and using our conventions, E b and m should be positive,indicating an attractive gravitational force and particleswith positive energy density, respectively. In the nextsubsection we consider the non-relativistic limit of thebinding energy between two scalar masses. C. Binding energy in the non-relativistic limit
In this subsection we consider the energy levels associ-ated with gravitational bound states. We follow Ref. [34]in assuming that the mass of the constituent particles ismuch lighter than the Planck mass such that the parti-cles are weakly coupled. In this case the magnitude ofthe binding energy is well below the mass of the con-stituent particles and we can work in the nonrelativis-tic limit. In that limit the effective description for thebound states of two scalar masses by gravity is simplygiven by Schr¨odinger’s equation with Newton’s gravita-tional potential. The bound state solution is then justthe Schr¨odinger solution for positronium, but with thefundamental constants that appear in the Coulomb po-tential exchanged for those of Newton’s potential.Thus we consider the Schr¨odinger equation, −∇ ψ ( r, θ, φ ) + 2 µ ( U ( r ) − E ) ψ ( r, θ, φ ) = 0 , (20) with potential U ( r ) = − Gm r , (21)where m is the particle mass, µ is the reduced mass equalto m/ G isNewton’s constant. The solution to this equation for apotential of this form is well known. The energy eigen-values, given by the principal quantum number n , are E n = G m n (22)in units where (cid:126) = c = 1.The ground state binding energy for two identicalscalar masses held together by gravity is therefore E = m G /
4, giving us a clear prediction for the dependenceof the binding energy on the constituent mass. However,in order to make contact with this power law using latticemethods we first need to take the continuum, infinite-volume limit of our calculation. To get an idea for whattype of discretization and finite-volume effects we mightexpect, we recall that the fractal dimension of the geome-tries varies as a function of distance scale. Even at thelargest distance scales that we can probe, the spectraldimension measured on the lattices does not get muchabove three. It is only in the continuum, infinite-volumelimit that the spectral dimension at long distances ex-trapolates to four, as it must to reproduce our world.In a 2+1 dimensional world, the same derivation for thebinding energy would lead to E ∝ m , so we might ex-pect a large extrapolation of the exponent in the powerlaw, from ∼ ∼
5, if the scalar fields on our latticessee a long distance effective dimension of around three atcoarse couplings and small volumes. This expectation isborn out by the data, as we show in the following section.To further interpret our results, it is useful to modelthe binding energy as a function of the dimension d . Bydimensional analysis, we have E ∝ m in 1+1 dimen-sions. Taking a simple quadratic fit to the three valuesof the exponent quoted here for integer dimension leadsto the scaling E ∝ m α , with α = d − d + 5 . (23)The values of d that determine this dependence interpo-late over the range of effective long distance dimensionsthat were found for our lattices in Ref. [9], so this simpleestimate of the relationship between α and d should bereliable enough to determine what numerical value of d is implied by the α obtained from our simulations.The other limit that must be taken in our calculationis the non-relativistic limit. Again following Ref. [34],we can get a rough estimate of the size of relativisticcorrections to Eq. (22). Consider the Hamiltonian H = 2 (cid:112) m + p − Gm r . (24)If we replace the momentum p by 1 /r and minimize theenergy, we find E b = 2 m − m (cid:114) − G m . (25)This implies that Gm / (cid:28) IV. NUMERICAL RESULTS
In this section we discuss the details of the numeri-cal analysis of the correlators, which yields the bindingenergy and renormalized mass. We then explain the de-tails of the power-law fits to the binding energy versusthe renormalized mass. Finally, we discuss the infinite-volume, continuum extrapolation of the exponent in thepower law as well as our determination of Newton’s con-stant.
A. Correlation functions
The correlators defined in Eqs. (12) and (13) are ob-tained from exact inversions of the matrix L xy calculatedon a given lattice configuration for a given bare massvalue. We do not consider every simplex on the lattice asa possible source. Instead, we vary the number of sourceson each configuration from one, five, 20, and 60 in orderto assess the effect the number of source simplices usedhas on the statistical error. We find that for a largenumber of sources, say 60, the systematic errors associ-ated with modeling the deviations of the data from themodel fit function are dominant. These deviations couldbe due to excited states, and finite-size and discretizationeffects. In order to avoid this difficulty of modeling thesystematics, we use a single source per configuration inour main analysis. All source simplices are selected ran-domly from the largest three-volume cross-section of theentire lattice. This is done by first shelling a configura-tion starting from a source chosen at random, and thenonly selecting sources for the propagator from the largestslice in the shelling. We find that restricting our sourcesto come from the largest three-slice minimizes finite lat-tice spacing effects, and it is the same procedure that wehave used in previous work on the spectral dimension [9]and for our studies of K¨ahler-Dirac fermions [30].An example of correlator data is shown in Fig. 2, andan example of the ratio F defined in Eq. (15) is shown inFig. 3. Both plots use log-linear coordinates and show re-sults for several masses. These figures display a featurethat appears across all ensembles, to varying degrees. Wesee around r ≈
10 lattice spacings a bend in the corre-lator, and a peak where the derivative of log[ F ] changessign. We also see this feature is a function of the bare r l o g [ G ( r )] m = 0.01 m = 0.02 m = 0.03 m = 0.04 m = 0.05 FIG. 2. The logarithm of the two-point correlator, G ( r ), forfive masses on the N = 8000, β = 0 ensemble. We see abend in the data which is pushed progressively out to larger r values as the bare mass is increased. The distances displayedon the horizontal axis are in units of the distance between thecenters of adjacent four-simplices, i.e. a dual edge length. r l o g [ F ( r )] m = 0.01 m = 0.02 m = 0.03 m = 0.04 m = 0.05 FIG. 3. The logarithm of the ratio between the two-particle,two-point correlator, and the square of the one-particle, two-point correlator, F ( r ). Five different bare masses are shownon the N = 8000, β = 0 ensemble. We see a peak in thedata, separating a positively sloped region and a negativelysloped region, which is pushed to larger r values for larger baremass values. The distances displayed on the horizontal axisare in units of the distance between the centers of adjacentfour-simplices, i.e. a dual edge length. mass, and as the bare mass is increased, this bend ispushed out further to larger distances. The same thinghappens as the volume of the system is increased. InFig. 4 we see the peak is pushed to larger distances asthe system volume is increased. This indicates that theturn-over in the data is most likely due to long-distancelattice artifacts. As noted in Ref. [9], there are baby uni-verses that branch off of the mother universe, where thebaby universes can be quite long, although their cross-section is of order the lattice spacing. This effect is mostpronounced on our coarsest lattices, where it can signif- r l o g [ F ( r )] N = 2000 N = 4000 N = 8000 N = 16000 N = 32000 FIG. 4. The logarithm of the ratio of the two-particle correla-tor to the square of the one-particle correlator as a function ofdistance for multiple volumes at β = 0, and m = 0 .
02. Wesee as the volume is increased the peak is pushed to largervalues of r indicating the turn-over in the data is most likelya long-distance lattice artifact. The distances displayed onthe horizontal axis are in units of the distance between thecenters of adjacent four-simplices, i.e. a dual edge length. icantly modify long-distance physics, although the effectappears to vanish in the continuum limit. It is useful tokeep this in mind when choosing a fit window to extractmasses from our correlation functions, since it sets anupper bound on how far in the Euclidean time extent wecan fit and still expect our model fit function to describethe data. This bend in the data that we see is most likelydue to baby-universe effects at long distances. At shortdistance scales we expect the usual discretization effects,as well as excited state contamination. Thus, the fit win-dow to the correlation function is rather constrained inour current approach.The peak in log[ F ] is one of the first clues on how toextract a physically motivated answer for the binding en-ergy. From the definition of E b in Eq. (17) the coefficientof r should be positive if the two-particle state is bound.We see this is only possible between r = 0 and the peakaround r ≈
10 lattice spacings (for the specific ensemblein Figs. 2 and 3). In fact, the existence of such a regionis already encouraging, since it implies there exists anattractive force between scalar masses inside the dynam-ical triangulations framework. This was noticed alreadyin Ref. [34].Additionally, looking at Fig. 3 we can see a change inconcavity for the larger masses around a value of r ≈ i.e. the peak. This inflection pointdenotes the end of the valid fitting region according toEq. (18), since after this point the long-distance effectsbegin to dominate the shape of the function. Across allbare masses and ensembles, we fit to a region that beginsat r = 1 and ends around the inflection point of log[ F ].We fit this same range in the one-particle correlator, and log[ G ] p -values FIG. 5. A histogram of the p -values extracted from fits tolog[ G ]. This histogram contains the p -values from the indi-vidual jackknife fits for all ensembles and mass values useddownstream in the analysis, combined into a single data set. the F function.The choice of fit function is decided by Eqs. (18)and (19). Thus, we use a function of the form, f ( r ) = Xr + Y + Z log r (26)for both the log[ F ] and log[ G ] data, with X , Y , and Z asfit parameters. By fitting the log[ F ] and log[ G ] data tothe functional form in Eq. (26) we can extract the bindingenergy, the renormalized mass, and the exponents β and γ as a function of the bare mass.The fits are done with non-linear least squares fittingincluding the correlations of the dependent data. The er-rors are estimated using single-elimination jackknife re-sampling, including the off-diagonal terms in the corre-lation matrix. The size of autocorrelation errors is esti-mated using a blocking procedure; the data is blocked un-til the errors no longer increase. In order to retain enoughinformation to resolve the correlation matrix when per-forming fits, the data is not blocked, but the errors areinflated to reflect the increased error due to autocorrela-tions. The fits are performed under a jackknife, and thecorrelation matrix is reconstructed for each individual fitunder the jackknife from the data on each jackknife sub-ensemble. By including correlations in the fit, the χ per degree of freedom is expected to be a reliable mea-sure of goodness of fit. We compute from the χ andthe number of degrees of freedom a confidence interval(a p -value) for the fit, correcting for finite sample size.We make a histogram of p -values from the fits for whichthe fit parameters are propagated through to the rest ofthe analysis. This includes fits from all ensembles. Theresulting histogram of p -values is relatively uniform, andis shown in Figs. 5, and 6 for the correlator fits, and F fits, respectively. Only the lowest bin possesses a smallspike. Since the fits in this bin are scattered throughoutthe parameter values of the analysis more or less at ran-dom, we do not ascribe an additional error to this slight log[ F ] p -values FIG. 6. A histogram of the p -values extracted from fits tolog[ F ]. This histogram contains the p -values from the indi-vidual jackknife fits for all ensembles and mass values useddownstream in the analysis, combined into a single data set. r l o g [ F ( r )] FitFit rangeData
FIG. 7. An example of a fit to the N = 16 ,
000 simplexensemble with β = − .
776 for log[ F ] at bare mass m = 0 . χ / d.o.f = 0 .
77 for the fit at this mass andthis ensemble corresponding to a p -value of 0 . deviation from a flat distribution. An example of the fitfor the N = 16 ,
000 simplex ensemble with β = − . F ] and log[ G ] data can be seen in Figs. 7 and 8.Given the results for the binding energy and the renor-malized mass for a wide range of bare mass values onmany different ensembles, we are able to test the the-ory presented in Sect. III C. This is done in the followingsubsections. B. Mass dependence of the binding energy
The dependence of the renormalized mass on bare massis shown in Fig. 9 for four different volumes at fixed lat-tice spacing ( β = 0). It is clear from this plot that therenormalized mass goes to zero as the bare mass also r l o g [ G ( r )] FitFit rangeData
FIG. 8. An example of a fit to the N = 16 ,
000 simplexensemble with β = − .
776 for log[ G ] at bare mass m =0 . χ / d.o.f = 0 .
19 for the fit at this massand this ensemble corresponding to a p -value of 0 . approaches zero, which is a consequence of the shift sym-metry of the lattice action. This provides a useful checkof our calculation.The dependence of the binding energy on the renor-malized mass is shown in Figs. 10, 11, 12 and 13 for fourdifferent ensembles. In order to make contact with New-tonian gravity, we must look for a power-law dependencefor the binding energy as a function of renormalized mass,as given in Eq. (22). As a first step, in order to eventu-ally be able to compare results across lattice spacings, weput the results in the same lattice spacing units. Beforeperforming the fits we re-scale all the binding energiesand renormalized masses to that of the fiducial latticespacing at β = 0 using the relative lattice spacings givenin table I.To study the power-law behavior, we must also deter-mine a fit window for the masses m , to which we fit thedata. To find the beginning of a fit range, we search forthe smallest bare mass for which the expected physicalinflection exists in the quantity log F ( e.g. in Fig. 3, ≥ m = 0 . m corresponding to thisbare mass is the beginning of the fit window. The endof the fit window is determined by a change of inflectionin the plots of binding energy versus renormalized mass.This point is where the nonrelativistic power-law behav-ior has been overtaken by effects due to strong couplingat larger mass values.For our power-law fits, we assume the functional form E b = Am α (27)where A and α are fit parameters, which in the contin-uum, nonrelativistic limit are expected to be A = G / α = 5, as given in Eq. (22). We find that this sim-ple fit function is a good description of the data on all0of our ensembles, with two exceptions: our two coarsestensembles ( β = 1 . β = 0 . E b = A | x − B | α + C (28)with A , B , α , and C the fit parameters. As before, A and α can be identified with their continuum, infinite-volume counterparts in Eq. (22). For these two ensemblesthe criteria for selecting the starting mass value of thefit is never satisfied i.e. the correct inflection in log[ F ]is not observed for any bare mass. This is most likelydue to large discretization errors on these coarse latticesmasking the physical behavior. Therefore the start ofthe fit window is somewhat arbitrary on these ensembles.We choose a fit range that begins in the region wherethe binding energy trends negative and ends before theinflection in the E b versus m plot at larger masses. Wevary this fit range to include a systematic error due tothis choice.The form in Eq. (27)—even at finite lattice spacing—isre-enforced by the existence of the shift symmetry, whichensures that the bare mass is only multiplicatively renor-malized, and hence, the binding energy is strictly pro-portional to some power of the renormalized mass. InFigs. 10, 11, and 12 we show examples of the binding en-ergy plotted against the renormalized mass, along witha best fit line and the fit range used (in black), for threedifferent lattice spacings. These are the finer lattices at (cid:96) rel = 0 . (cid:96) rel = 0 .
8, and one of the coarser ensem-bles at (cid:96) rel = 1, respectively. In Figs. 13 and 14 we showthe same quantities for the extra coarse, (cid:96) rel = 1 .
59 en-semble, and (cid:96) rel = 1 .
28 ensemble. We see good agreementbetween the fit functions Eq. (27) and (28), and the data.These fits are done by taking the correlations in thedata into account. We use weighted orthogonal dis-tance regression [45, 46] to incorporate correlations inthe renormalized mass and in the binding energy datasets simultaneously. For the weights we use the inversecovariances in both data sets to obtain the best fit tothe data points using χ minimization. A detailed dis- m m N = 4000 N = 8000 N = 16000 N = 32000 FIG. 9. The renormalized mass plotted against the baremass for four volumes with β = 0. Here we see the renor-malized mass is multiplicatively renormalized for sufficientlysmall bare masses. m E b FitDataFit range
FIG. 10. The power-law fit to the binding energy plottedagainst the renormalized mass for the N = 16 , β = − .
776 ensemble. The fit range is shown in black, and thesolid line is the fit to the data. The fit corresponds to a χ /d.o.f. = 0 .
59, with a p -value of 0.62. cussion of this procedure can be found in appendix A.To assess a systematic error associated with the choice offit range, we vary the start and end points of a fit rangeover a reasonable set of values guided by the quality of fitand tabulate the results. We then calculate the standarddeviation of those results and include it as a systematicerror, adding it in quadrature to the statistical error ofthe result from the central fit to give a total error.We perform this power-law fit across all of our ensem-bles, extracting a power α and a coefficient A . From thatcoefficient A we calculate √ A , which we associate witha value for G at fixed volume and lattice spacing. Whilethis association with G at finite lattice spacing or volumemay suffer from systematic errors, in the continuum, in-finite volume limit the quantity √ A should extrapolateto G . With these results for α and G across ensembles,1 m E b FitDataFit range
FIG. 11. The power-law fit to the binding energy plottedagainst the renormalized mass for the N = 4 , β = − . χ /d.o.f. = 0 . p -value of 0.79. m E b FitDataFit range
FIG. 12. The power-law fit to the binding energy plottedagainst the renormalized mass for the N = 16 , β = 0ensemble. The fit range is shown in black, and the solid line isthe fit to the data. The fit corresponds to a χ /d.o.f. = 0 . p -value of 0 . we are able to obtain their continuum, infinite-volumevalues. C. Continuum, infinite volume extrapolation
To preform the extrapolation to the continuum,infinite-volume limit we take the simplest ansatz sug-gested from finite-size scaling, and the discretization de-pendence suggested by the symmetries of the theory.This approach is similar to what has been used in ourprevious analyses on dynamical triangulations, with thechoice of fit functions given by α = H α V + I α (cid:96) + J α V + K α (cid:96) + L α (29) m E b FitDataFit range
FIG. 13. The power-law fit to the binding energy plottedagainst the renormalized mass for the N = 4 , β = 1 . χ /d.o.f. = 1 . p -value of 0 . m E b FitDataFit range
FIG. 14. The power-law fit to the binding energy plottedagainst the renormalized mass for the N = 4 , β = 0 . χ /d.o.f. = 1 . p -value of 0 . and G = H G V + I G (cid:96) + J G V + K G (cid:96) + L G , (30)where H i , I i , J i , K i , and L i are fit parameters fortheir respective quantities. Here V is the physical vol-ume, and (cid:96) rel is the relative lattice spacing. We includequadratic corrections in the inverse physical volume, andthe squared lattice spacing, since we find curvature inour data when we include small volumes and coarse lat-tice spacings. In addition to the fit ansatzes in Eqs. (29)and (30), we also perform fits dropping the ∼ (cid:96) term,which we are able to do when we also drop the two coars-est lattice spacings. The results from this fit are consis-tent within one-sigma to the results of the fit to the fulldata set including the (cid:96) term.2 V rel = 0 rel = 0.7 rel = 0.8 rel = 1 rel = 1.28 rel = 1.59 FIG. 15. The power α as a function of the inverse physicalvolume (expressed in units of 1000 simplices) for all of theensembles (colored), as well as the continuum limit (in black).Here quadratic corrections in 1 /V as well as (cid:96) were used tomodel the extrapolation. For this fit we find χ / d.o.f. = 0 . p -value of 0.73, and the continuum, infinitevolume value is α = 4 . V G rel = 0 rel = 0.7 rel = 0.8 rel = 1 rel = 1.28 rel = 1.59 FIG. 16. Newton’s constant G as a function of the inversephysical volume (expressed in units of 1000 simplices) for allof the ensembles (colored), as well as the continuum limit(in black). Here quadratic corrections in 1 /V as well as (cid:96) were used to model the extrapolation. For this fit we find χ / d.o.f. = 0 .
37 corresponding to a p -value of 0.87, and thecontinuum, infinite volume value is G = 15(5). The extrapolations for α and for G are shown inFigs. 15 and 16, respectively, where they are plottedagainst the inverse physical volume. There, lines ofconstant lattice spacing are drawn, and the zero-lattice-spacing limit is shown in black. In Figs. 17 and 18 weshow the same extrapolations to the infinite-volume, con-tinuum limit as the fits shown in Figs. 15 and 16, but witha different cross-section through the two-dimensional pa-rameter space. In these plots we show the values for α and G plotted versus the squared lattice spacing. There,lines of constant physical volume are drawn with theinfinite-volume limit shown as a solid black line. Note rel V = 2 000 V = 8 000 V = 32 000 FIG. 17. The same data and fit from Fig. 15 however nowplotted as a function of the squared lattice spacing. Hereexample lines of constant physical volume are plotted alongwith the infinite volume limit as a solid black line, and thedata are represented in the same manner as Fig. 15. rel G V = 2 000 V = 8 000 V = 32 000 FIG. 18. The same data and fit from Fig. 16 however nowplotted as a function of the squared lattice spacing. Hereexample lines of constant physical volume are plotted alongwith the infinite volume limit as a solid black line, and thedata are represented in the same manner as Fig. 16. that the continuum, infinite-volume limit is taken sepa-rately for α and G .These fits are performed assuming that the data pointsare uncorrelated, which is reasonable given that thepoints are all from different ensembles. For the α ex-trapolation we find the χ / d.o.f = 0 .
56, correspond-ing to a p -value of 0.73. For the G extrapolation wefind χ / d.o.f = 0 .
37, corresponding to a p -value of 0.87.The infinite volume, continuum extrapolated values are α = 4 . G = 15(5). The errors on these quanti-ties are fairly large, around 20%-30% but the value of α is consistent with what is needed to recover the Newto-nian limit. In fact, this result is somewhat better thanit at first appears because of the sensitive dependence of α on the effective dimension. Given Eq. (23) for the de-3pendence of α on dimension, one finds that our 1 σ errorband on α implies a 1 σ error range for the dimension ofbetween 3 . .
1. Thus our result implies Newtonianbinding in a dimension quite close to 4 at long distancescales in the continuum limit. Note that at coarser cou-pling and smaller lattice volumes where the effective di-mension of the lattice geometries is around 3 or lower,the value of α drops to 2 or lower, which is what we ex-pect given the dependence of α on dimension. Thus ourresults are consistent with expectations.The value of G that we should expect is not known apriori , since it sets the lattice spacing in physical units,but it should satisfy certain consistency checks, as dis-cussed in Sect. III C. The constraint that the binding benonrelativistic implies that Gm / (cid:28)
1, which trans-lates into m (cid:28) .
13 in our fiducial lattice units, givenour value of G in those units. The upper range of massesin our fit window is at m ≈ .
02, so the nonrelativisticcondition is satisfied.One particularly nice feature of this calculation is thatour value of G allows us to determine the lattice spacingin units of the Planck length for the first time. Becauseour extrapolation procedure recovers the correct Newto-nian limit, we can have some confidence that the valueof G computed via this method is reliable. We find that √ G = (cid:96) P l = (3 . ± . (cid:96) fid . Thus, our fiducial latticespacing is around 1/4 the Planck length. Our finest lat-tice spacing is around 1/6 the Planck length, so we seethat the lattice spacing can indeed be made smaller thanthe Planck length. V. DISCUSSION & CONCLUSION
One of the tests that any formulation of lattice gravitymust pass is that it must have the correct classical limit.In this work we have shown that the interaction of scalarparticles is well-described by Newton’s potential in theappropriate non-relativistic, classical limit. This conclu-sion comes from our study of the binding energy of thetwo particle bound state as a function of the constituentscalar particle mass. The analysis makes use of a numberof ensembles with multiple volumes and lattice spacings,allowing us to extrapolate our results to the continuum,infinite volume limit.Our calculation passes a number of non-trivial cross-checks. We show numerically that the renormalizedscalar mass approaches zero as the bare mass approacheszero, which is expected given the shift symmetry of thelattice action. We study the two-particle binding energyas a function of its constituent mass, and we find thatit is well-described by a power law in the non-relativisticlimit. At finite lattice spacing and finite volume the ex-ponent in the power law is close to what one expects ifthe effective dimension of the geometry is near the mea-sured values for the effective dimension on these samelattices in Ref. [9]. Only in the continuum, infinite vol-ume limit does the power law correspond to that of the Newtonian potential in four dimensions. The consistencyof the binding with Newtonian gravity allows us to ex-tract a value of G from the calculation, and knowing G allows us to determine the Planck scale for the firsttime within EDT. We find that the Planck length is (cid:96) P l = (3 . ± . (cid:96) fid , so that our fiducial lattice spac-ing is about 1/4 the Planck length, and our finest latticespacing is around 1/6 the Planck length. This shows thatas the continuum limit is approached the lattice spacingbecomes a decreasing fraction of the Planck length, sug-gesting that there is no barrier to taking a continuumlimit.Although the calculation is done in the quenched ap-proximation, such that the back-reaction of the scalarfield on geometry is neglected, we still expect the bind-ing of particles to be governed by tree-level graviton ex-change. The effects of quenching only appear at one-loopin the low-energy effective theory, so they should not af-fect the recovery of the classical limit, as seen here.Our recovery of the Newtonian potential within EDTis a strong cross-check of all of the ingredients that gointo the formulation of lattice gravity used here, fromour determination of the relative lattice spacing to ourapproach to taking the continuum, infinite-volume limit.This non-trivial result provides powerful evidence thatEDT is in fact a theory of gravity, and that investiga-tions of the formulation as a possible realization of theasymptotic safety scenario should continue. ACKNOWLEDGMENTS
The authors thank Claude Bernard and Simon Catter-all for valuable discussions, and we thank Simon Catteralland Fleur Versteegen for comments on the manuscript.JL was supported by the U.S. Department of Energy(DOE), Office of Science, Office of High Energy Physicsunder Award Number DE-SC0009998. JUY was sup-ported by the U.S. De- partment of Energy grant DE-SC0019139, and by Fermi Research Alliance, LLC underContract No. DE-AC02- 07CH11359 with the U.S. De-partment of Energy, Office of Science, Office of High En-ergy Physics. MS is supported by the German AcademicScholarship Foundation and gratefully acknowledges hos-pitality at Syracuse University and at CP3-Origins, Uni-versity of Southern Denmark, during various stages ofthis project. Computations were performed in part onthe Syracuse University HTC Campus Grid and weresupported by NSF award ACI-1341006. Computationsfor this work were also carried out in part on facilitiesof the USQCD Collaboration, which are funded by theOffice of Science of the U.S. Department of Energy.
Appendix A: Correlated fitting
When fitting the binding energy versus the renormal-ized mass, both quantities come with errors, and the dif-4ferent data points are actually correlated, both in the x and y data-sets. This is because the fits for the renor-malized masses and binding energies—while at differentbare input masses—used the same configurations for eachfit, thus introducing correlations between fits. In orderto account for this, we did simultaneous correlated fitsin both the x and y data sets using orthogonal distanceregression. Orthogonal distance regression attempts tominimize the radial distance between the fit function andthe input data—as opposed to the vertical distance be-tween the fit function and the y data that is typical innon-linear least squares fitting.The idea is the following: One seeks to minimize thefunction, Q ( { ¯ x } , { β } ) = (A1) (cid:88) i,j ( f (¯ x i , β k ) − y i ) w − ij ( f (¯ x j , β k ) − y j ) + (A2)(¯ x i − x i ) (cid:15) − ij (¯ x j − x j ) (A3)where ¯ x i is the collection of ideal x values which mini-mizes the expression, f is the fitting function, β k is thecollection of fit parameters, y i is the collection of y data, x i is the collection of x data, w − is the inverse of thecovariance matrix for the y data, and (cid:15) − is the inverseof the covariance matrix for the x data. The problemcan be recast, with a change of variables, into a different expression to be minimized, Q ( { δ } , { β } ) = (cid:88) i,j (A4)( f ( x i + δ i , β k ) − y i ) w − ij ( f ( x j + δ j , β k ) − y j ) + (A5) δ i (cid:15) − ij δ j . (A6)Now, there are two sets of parameters to minimize, the δ i , and the β i , where δ i = ¯ x i − x i is the residual in the x data.The above expression can be recast into a very simpleform for numerical purposes. The w − and (cid:15) − can becombined into a single block-diagonal matrix, C − = w − ⊕ (cid:15) − , (A7)and the residuals can be combined into a single vector, v = ( f ( x + δ, β k ) − y ) ⊕ δ. (A8)This allows the above function to be written as, Q ( { δ } , { β } ) = (cid:88) i,j v i C − ij v j . (A9)For this work, the covariance matrices are positive defi-nite, so we can take the above expression one step furtherand perform a Cholesky decomposition, C − ij = (cid:88) k L ik L Tkj (A10)which in turn allows us to express Q as a simple dot prod-uct between vectors, Q = (cid:80) i r i r i , with r i = (cid:80) k v k L ki .To minimize Q we use the Levenberg-Marquardt algo-rithm , and minimize with respect to the δ i and the β i .The resulting β values are the desired fit parameters forthe input function, f . [1] S. Weinberg, “ULTRAVIOLET DIVERGENCES INQUANTUM THEORIES OF GRAVITATION,” in Gen-eral Relativity: An Einstein Centenary Survey (1980) pp.790–831.[2] G. ’t Hooft and M. J. G. Veltman, Ann. Inst. H. PoincarePhys. Theor. A , 69 (1974).[3] M. H. Goroff and A. Sagnotti, Nucl. Phys. B , 709(1986).[4] J. F. Donoghue, Phys. Rev. Lett. , 2996 (1994),arXiv:gr-qc/9310024.[5] J. F. Donoghue, Phys. Rev. D , 3874 (1994), arXiv:gr-qc/9405057.[6] J. Ambjorn, J. Jurkiewicz, and R. Loll, Phys. Rev. D , 064014 (2005), arXiv:hep-th/0505154.[7] J. Ambjorn and R. Loll, Nucl. Phys. B , 407 (1998),arXiv:hep-th/9805108.[8] J. Ambjorn, A. Gorlich, J. Jurkiewicz, and R. Loll, Phys.Rev. D , 063544 (2008), arXiv:0807.4481 [hep-th].[9] J. Laiho, S. Bassler, D. Coumbe, D. Du, and J. T. Neelakanta, Phys. Rev. D , 064015 (2017),arXiv:1604.02745 [hep-th].[10] C. Wetterich, Phys. Lett. B , 90 (1993),arXiv:1710.05815 [hep-th].[11] U. Ellwanger, , 206 (1993), arXiv:hep-ph/9308260.[12] T. R. Morris, Int. J. Mod. Phys. A , 2411 (1994),arXiv:hep-ph/9308265.[13] M. Reuter, Phys. Rev. D , 971 (1998), arXiv:hep-th/9605030.[14] M. Niedermaier, Class. Quant. Grav. , R171 (2007),arXiv:gr-qc/0610018.[15] R. Percacci, , 111 (2007), arXiv:0709.3851 [hep-th].[16] D. F. Litim, PoS QG-Ph , 024 (2007), arXiv:0810.3675[hep-th].[17] M. Reuter and F. Saueressig, New J. Phys. , 055022(2012), arXiv:1202.2274 [hep-th].[18] A. Eichhorn, Front. Astron. Space Sci. , 47 (2019),arXiv:1810.07615 [hep-th].[19] A. D. Pereira, in Progress and Visions in Quantum The- ory in View of Gravity: Bridging foundations of physicsand mathematics (2019) arXiv:1904.07042 [gr-qc].[20] M. Reichert, PoS , 005 (2020).[21] J. M. Pawlowski and M. Reichert, (2020),arXiv:2007.10353 [hep-th].[22] A. Bonanno, A. Eichhorn, H. Gies, J. M. Pawlowski,R. Percacci, M. Reuter, F. Saueressig, and G. P. Vacca,Front. in Phys. , 269 (2020), arXiv:2004.06810 [gr-qc].[23] J. Ambjorn and J. Jurkiewicz, Phys. Lett. B , 42(1992).[24] M. E. Agishtein and A. A. Migdal, Mod. Phys. Lett. A , 1039 (1992).[25] S. Catterall, J. B. Kogut, and R. Renken, Phys. Lett. B , 277 (1994), arXiv:hep-lat/9401026.[26] B. Bruegmann and E. Marinari, Phys. Rev. Lett. ,1908 (1993), arXiv:hep-lat/9210002.[27] J. Ambjorn, J. Jurkiewicz, and R. Loll, Phys. Rev. Lett. , 171301 (2005), arXiv:hep-th/0505113.[28] O. Lauscher and M. Reuter, JHEP , 050 (2005),arXiv:hep-th/0508202.[29] S. Carlip, Class. Quant. Grav. , 193001 (2017),arXiv:1705.05417 [gr-qc].[30] S. Catterall, J. Laiho, and J. Unmuth-Yockey, Phys.Rev. D , 114503 (2018), arXiv:1810.10626 [hep-lat].[31] E. Kahler, Rend. Math , 425 (1962).[32] T. Banks, Y. Dothan, and D. Horn, Physics Letters B , 413 (1982).[33] A. Bazavov et al. (MILC), Rev. Mod. Phys. , 1349 (2010), arXiv:0903.3598 [hep-lat].[34] B. V. de Bakker and J. Smit, Nucl. Phys. B , 476(1997), arXiv:hep-lat/9604023.[35] J. Smit, (2021), arXiv:2101.01028 [gr-qc].[36] S. Bilke, Z. Burda, A. Krzywicki, B. Petersson,J. Tabaczek, and G. Thorleifsson, Phys. Lett. B ,279 (1998), arXiv:hep-lat/9804011.[37] T. Regge, Nuovo Cim. , 558 (1961).[38] S. Bilke and G. Thorleifsson, Phys. Rev. D , 124008(1999), arXiv:hep-lat/9810049.[39] M. Gross and S. Varsted, Nucl. Phys. B , 367 (1992).[40] S. Catterall, Comput. Phys. Commun. , 409 (1995),arXiv:hep-lat/9405026.[41] D. Coumbe and J. Laiho, JHEP , 028 (2015),arXiv:1401.3299 [hep-th].[42] J. Ambjorn, L. Glaser, A. Goerlich, and J. Jurkiewicz,JHEP , 100 (2013), arXiv:1307.2270 [hep-lat].[43] R. G. Jha, J. Laiho, and J. Unmuth-Yockey, PoS LAT-TICE2018 , 043 (2018), arXiv:1810.09946 [hep-lat].[44] S. Catterall, J. Laiho, and J. Unmuth-Yockey, JHEP ,013 (2018), arXiv:1806.07845 [hep-lat].[45] P. T. Boggs and J. E. Rogers, Contemporary Mathemat-ics112