# 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 Schiﬀer

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 ﬁelds 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 ﬁnd 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 inﬁnite-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 ﬁrst time, and we ﬁnd 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 eﬀective theory canbe formulated [4, 5], it comes with an inﬁnite 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 ﬁxed point,making it eﬀectively renormalizable nonperturbatively,with only a ﬁnite number of parameters needing to bespeciﬁed from experiment. One would also hope that thenumber of relevant parameters needing to be speciﬁedfrom 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 ﬁxed point would appear as acontinuous phase transition, the approach to which de-ﬁnes a continuum limit. This is the ﬁrst 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 brieﬂy review the evidence so far that EDT mayrealize the asymptotic safety scenario for gravity. Ref. [9]introduced a ﬁne 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 ﬁrst 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 ﬁrst order line to a possible critical endpoint,and ensembles following this procedure were generatedat a number of diﬀerent lattice spacings. The globalHausdorﬀ dimension was measured using ﬁnite-volumescaling and shown to be close to four [9]. The spectraldimension, which is a fractal dimension deﬁned 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 ismodiﬁed by discretization eﬀects, but this type of eﬀecton long-distance behavior is typical when a symmetry ofa theory is broken by the regulator, for example by theﬁnite lattice spacing in the case of lattice regularization.Then a ﬁne-tuning is needed to approximately restore thesymmetry at ﬁnite lattice spacing. Residual symmetry-breaking eﬀects appear at ﬁnite 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 ﬁne-tuning is requiredto restore chiral symmetry, and even then, at ﬁnite 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 diﬀeo-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 ﬂat-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 eﬀects 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 eﬀects 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 matterﬁelds living on dynamical triangulations, this time focus-ing on scalar ﬁelds. We follow the work of de Bakkerand Smit [34] in our implementation of scalar ﬁelds.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 conﬁrm the ﬁnding that there is anattractive force between scalar particles, and after tak-ing the continuum, inﬁnite volume limit, we ﬁnd 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 ﬁnd thatour lattice spacings are smaller than the Planck lengthand that for our ﬁnest 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 ﬁelds 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 theﬁtter used in the analysis. II. EUCLIDEAN DYNAMICALTRIANGULATIONSA. The model

In the continuum, for ﬁxed 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 ﬁeld 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 ﬁeld.The partition function for lattice quantum gravity thatwe use in this work is that of Euclidean dynamical trian-gulations, where the path integral at ﬁnite 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 conﬁgura-tions used in this work are given in Ref. [9], and they aresummarized brieﬂy 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 ﬁnite-size eﬀects 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 ﬁve vertex labels. Suchconﬁgurations 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. Wedeﬁne a sweep as a ﬁxed 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 ﬁxed to S .The topology remains ﬁxed 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 inﬁnite lat-tice volume limit it is necessary to ﬁne-tune the bareparameter κ to its critical value. This amounts to atuning of the bare cosmological constant. Although thistuning is required to take the inﬁnite lattice volume limit,it does not guarantee that we can take the inﬁnite phys-ical volume limit. The inﬁnite 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-oﬀ 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 ﬁxed latticefour-volumes. The Pachner moves require the volume tovary in order for them to be ergotic, though in practicethe volume ﬂuctuations 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 ﬁducial value N f . Thisterm does not alter the action when N = N f , but it doeskeep the volume ﬂuctuations 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 signiﬁcant 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 ﬁrst 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 hasHausdorﬀ dimension 2, while the collapsed phase has alarge, possibly inﬁnite, 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 largeﬁnite-size eﬀects [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 ﬁne-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 globalHausdorﬀ 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 ﬁrst 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 ﬁnd that this selection procedure becomes un-necessary at the ﬁnest lattice spacing considered in thiswork, since we can simulate somewhat farther from thetransition line without losing the four-dimensional ﬁnite-volume scaling. This is fortunate, since future runs willbe at ever ﬁner lattice spacings, so will likely not suﬀerfrom 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 ﬁnest 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 inﬁnite-volume extrapolation.The determination of the relative lattice spacing fol-lows the procedure of Ref. [9], which looked at the returnprobability of a diﬀusion process on the lattice geome-tries. The return probability is dimensionless, but variesas a function of the diﬀusion time step, which is not.One can rescale the diﬀusion 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 diﬀer slightly from thoseﬁrst 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 reﬂect 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 ﬁelds on dynamicaltriangulations, then we show how to calculate the bind-ing energy between scalar particles on our lattices, andﬁnally we review the nonrelativistic limit of the bindingenergy in the continuum limit. (cid:96) rel β κ N Number of conﬁgurations1.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 ﬁrst column shows the relative lattice spacing, with theensembles at β = 0 serving as the ﬁducial lattice spacing.The quoted error is a systematic error associated with ﬁnite-volume eﬀects. The second column is the value of β , the thirdis the value of κ , the fourth is the number of four-simplices inthe simulation, and the ﬁfth is the number of conﬁgurationssampled. A. Scalar ﬁelds on dynamical triangulations

In order to map Eq. (3) on to the lattice, we use anaive discretization where we associate each scalar ﬁeldat x , φ ( x ), with a four-simplex (or equivalently a dual-lattice site). In this way each scalar ﬁeld is restrictedto only have ﬁve neighbors through the ﬁve 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 ﬁne 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 coeﬃcient 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 ﬁrst calculate the propagator for a freescalar ﬁeld 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 ﬁve), 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 ﬁrst is the average over all possibledistances a separation | x − y | = r away. In practice thisis not done in its entirety, and a ﬁxed 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 conﬁgurations.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-oﬀ 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 suﬃciently 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 diﬀerence in energy between a two-body system in its ground state and twice the mass ofa single particle in its ground state, with gravitationaleﬀects taken into account.With these deﬁnitions, taking the logarithm of F , weﬁnd 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 eﬀective 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 ﬁrst need to take the continuum, inﬁnite-volume limit of our calculation. To get an idea for whattype of discretization and ﬁnite-volume eﬀects 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, inﬁnite-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 ﬁelds on our latticessee a long distance eﬀective 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 ﬁt 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 eﬀective 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 ﬁnd 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 ﬁts to the binding energy versusthe renormalized mass. Finally, we discuss the inﬁnite-volume, continuum extrapolation of the exponent in thepower law as well as our determination of Newton’s con-stant.

A. Correlation functions

The correlators deﬁned in Eqs. (12) and (13) are ob-tained from exact inversions of the matrix L xy calculatedon a given lattice conﬁguration 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 conﬁguration from one, ﬁve, 20, and 60 in orderto assess the eﬀect the number of source simplices usedhas on the statistical error. We ﬁnd that for a largenumber of sources, say 60, the systematic errors associ-ated with modeling the deviations of the data from themodel ﬁt function are dominant. These deviations couldbe due to excited states, and ﬁnite-size and discretizationeﬀects. In order to avoid this diﬃculty of modeling thesystematics, we use a single source per conﬁguration inour main analysis. All source simplices are selected ran-domly from the largest three-volume cross-section of theentire lattice. This is done by ﬁrst shelling a conﬁgura-tion starting from a source chosen at random, and thenonly selecting sources for the propagator from the largestslice in the shelling. We ﬁnd that restricting our sourcesto come from the largest three-slice minimizes ﬁnite lat-tice spacing eﬀects, 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 deﬁned in Eq. (15) is shown inFig. 3. Both plots use log-linear coordinates and show re-sults for several masses. These ﬁgures 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 ), forﬁve 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 diﬀerent 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 oﬀ of the mother universe, where thebaby universes can be quite long, although their cross-section is of order the lattice spacing. This eﬀect 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 eﬀectappears to vanish in the continuum limit. It is useful tokeep this in mind when choosing a ﬁt window to extractmasses from our correlation functions, since it sets anupper bound on how far in the Euclidean time extent wecan ﬁt and still expect our model ﬁt function to describethe data. This bend in the data that we see is most likelydue to baby-universe eﬀects at long distances. At shortdistance scales we expect the usual discretization eﬀects,as well as excited state contamination. Thus, the ﬁt win-dow to the correlation function is rather constrained inour current approach.The peak in log[ F ] is one of the ﬁrst clues on how toextract a physically motivated answer for the binding en-ergy. From the deﬁnition of E b in Eq. (17) the coeﬃcientof 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 speciﬁc 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 inﬂection pointdenotes the end of the valid ﬁtting region according toEq. (18), since after this point the long-distance eﬀectsbegin to dominate the shape of the function. Across allbare masses and ensembles, we ﬁt to a region that beginsat r = 1 and ends around the inﬂection point of log[ F ].We ﬁt this same range in the one-particle correlator, and log[ G ] p -values FIG. 5. A histogram of the p -values extracted from ﬁts tolog[ G ]. This histogram contains the p -values from the indi-vidual jackknife ﬁts for all ensembles and mass values useddownstream in the analysis, combined into a single data set. the F function.The choice of ﬁt 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 asﬁt parameters. By ﬁtting 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 ﬁts are done with non-linear least squares ﬁttingincluding the correlations of the dependent data. The er-rors are estimated using single-elimination jackknife re-sampling, including the oﬀ-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 ﬁts, the data is not blocked, but the errors areinﬂated to reﬂect the increased error due to autocorrela-tions. The ﬁts are performed under a jackknife, and thecorrelation matrix is reconstructed for each individual ﬁtunder the jackknife from the data on each jackknife sub-ensemble. By including correlations in the ﬁt, the χ per degree of freedom is expected to be a reliable mea-sure of goodness of ﬁt. We compute from the χ andthe number of degrees of freedom a conﬁdence interval(a p -value) for the ﬁt, correcting for ﬁnite sample size.We make a histogram of p -values from the ﬁts for whichthe ﬁt parameters are propagated through to the rest ofthe analysis. This includes ﬁts from all ensembles. Theresulting histogram of p -values is relatively uniform, andis shown in Figs. 5, and 6 for the correlator ﬁts, and F ﬁts, respectively. Only the lowest bin possesses a smallspike. Since the ﬁts 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 ﬁts tolog[ F ]. This histogram contains the p -values from the indi-vidual jackknife ﬁts 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 ﬁt to the N = 16 ,

000 simplexensemble with β = − .

776 for log[ F ] at bare mass m = 0 . χ / d.o.f = 0 .

77 for the ﬁt at this mass andthis ensemble corresponding to a p -value of 0 . deviation from a ﬂat distribution. An example of the ﬁtfor 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 diﬀerent 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 diﬀerent volumes at ﬁxed 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 ﬁt to the N = 16 ,

000 simplexensemble with β = − .

776 for log[ G ] at bare mass m =0 . χ / d.o.f = 0 .

19 for the ﬁt 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 fourdiﬀerent 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 ﬁrst 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 ﬁts we re-scale all the binding energiesand renormalized masses to that of the ﬁducial latticespacing at β = 0 using the relative lattice spacings givenin table I.To study the power-law behavior, we must also deter-mine a ﬁt window for the masses m , to which we ﬁt thedata. To ﬁnd the beginning of a ﬁt range, we search forthe smallest bare mass for which the expected physicalinﬂection exists in the quantity log F ( e.g. in Fig. 3, ≥ m = 0 . m corresponding to thisbare mass is the beginning of the ﬁt window. The endof the ﬁt window is determined by a change of inﬂectionin the plots of binding energy versus renormalized mass.This point is where the nonrelativistic power-law behav-ior has been overtaken by eﬀects due to strong couplingat larger mass values.For our power-law ﬁts, we assume the functional form E b = Am α (27)where A and α are ﬁt parameters, which in the contin-uum, nonrelativistic limit are expected to be A = G / α = 5, as given in Eq. (22). We ﬁnd that this sim-ple ﬁt 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 ﬁt parameters. As before, A and α can be identiﬁed with their continuum, inﬁnite-volume counterparts in Eq. (22). For these two ensemblesthe criteria for selecting the starting mass value of theﬁt is never satisﬁed i.e. the correct inﬂection 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 ﬁt window is somewhat arbitrary on these ensembles.We choose a ﬁt range that begins in the region wherethe binding energy trends negative and ends before theinﬂection in the E b versus m plot at larger masses. Wevary this ﬁt range to include a systematic error due tothis choice.The form in Eq. (27)—even at ﬁnite 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 ﬁt line and the ﬁt range used (in black), for threediﬀerent lattice spacings. These are the ﬁner 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 ﬁt functions Eq. (27) and (28), and the data.These ﬁts 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 ﬁt 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 suﬃcientlysmall bare masses. m E b FitDataFit range

FIG. 10. The power-law ﬁt to the binding energy plottedagainst the renormalized mass for the N = 16 , β = − .

776 ensemble. The ﬁt range is shown in black, and thesolid line is the ﬁt to the data. The ﬁt 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 ofﬁt range, we vary the start and end points of a ﬁt rangeover a reasonable set of values guided by the quality of ﬁtand 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 ﬁt to give a total error.We perform this power-law ﬁt across all of our ensem-bles, extracting a power α and a coeﬃcient A . From thatcoeﬃcient A we calculate √ A , which we associate witha value for G at ﬁxed volume and lattice spacing. Whilethis association with G at ﬁnite lattice spacing or volumemay suﬀer from systematic errors, in the continuum, in-ﬁnite 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 ﬁt 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 ﬁt to the binding energy plottedagainst the renormalized mass for the N = 16 , β = 0ensemble. The ﬁt range is shown in black, and the solid line isthe ﬁt to the data. The ﬁt corresponds to a χ /d.o.f. = 0 . p -value of 0 . we are able to obtain their continuum, inﬁnite-volumevalues. C. Continuum, inﬁnite volume extrapolation

To preform the extrapolation to the continuum,inﬁnite-volume limit we take the simplest ansatz sug-gested from ﬁnite-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 ﬁt 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 ﬁt 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 ﬁt 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 ﬁt 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 ﬁnd curvature inour data when we include small volumes and coarse lat-tice spacings. In addition to the ﬁt ansatzes in Eqs. (29)and (30), we also perform ﬁts 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 ﬁt are consis-tent within one-sigma to the results of the ﬁt 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 ﬁt we ﬁnd χ / d.o.f. = 0 . p -value of 0.73, and the continuum, inﬁnitevolume 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 ﬁt we ﬁnd χ / d.o.f. = 0 .

37 corresponding to a p -value of 0.87, and thecontinuum, inﬁnite 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 inﬁnite-volume, con-tinuum limit as the ﬁts shown in Figs. 15 and 16, but witha diﬀerent 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 theinﬁnite-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 ﬁt from Fig. 15 however nowplotted as a function of the squared lattice spacing. Hereexample lines of constant physical volume are plotted alongwith the inﬁnite 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 ﬁt from Fig. 16 however nowplotted as a function of the squared lattice spacing. Hereexample lines of constant physical volume are plotted alongwith the inﬁnite volume limit as a solid black line, and thedata are represented in the same manner as Fig. 16. that the continuum, inﬁnite-volume limit is taken sepa-rately for α and G .These ﬁts are performed assuming that the data pointsare uncorrelated, which is reasonable given that thepoints are all from diﬀerent ensembles. For the α ex-trapolation we ﬁnd the χ / d.o.f = 0 .

56, correspond-ing to a p -value of 0.73. For the G extrapolation weﬁnd χ / d.o.f = 0 .

37, corresponding to a p -value of 0.87.The inﬁnite 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 ﬁrst appears because of the sensitive dependence of α on the eﬀective dimension. Given Eq. (23) for the de-3pendence of α on dimension, one ﬁnds 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 eﬀective 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 ﬁducial lattice units, givenour value of G in those units. The upper range of massesin our ﬁt window is at m ≈ .

02, so the nonrelativisticcondition is satisﬁed.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 ﬁrst time. Becauseour extrapolation procedure recovers the correct Newto-nian limit, we can have some conﬁdence that the valueof G computed via this method is reliable. We ﬁnd that √ G = (cid:96) P l = (3 . ± . (cid:96) ﬁd . Thus, our ﬁducial latticespacing is around 1/4 the Planck length. Our ﬁnest 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,inﬁnite 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 ﬁnd thatit is well-described by a power law in the non-relativisticlimit. At ﬁnite lattice spacing and ﬁnite volume the ex-ponent in the power law is close to what one expects ifthe eﬀective dimension of the geometry is near the mea-sured values for the eﬀective dimension on these samelattices in Ref. [9]. Only in the continuum, inﬁnite 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 ﬁrsttime within EDT. We ﬁnd that the Planck length is (cid:96) P l = (3 . ± . (cid:96) ﬁd , so that our ﬁducial lattice spac-ing is about 1/4 the Planck length, and our ﬁnest 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 scalarﬁeld on geometry is neglected, we still expect the bind-ing of particles to be governed by tree-level graviton ex-change. The eﬀects of quenching only appear at one-loopin the low-energy eﬀective 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, inﬁnite-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), Oﬃce of Science, Oﬃce 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, Oﬃce of Science, Oﬃce 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 theOﬃce of Science of the U.S. Department of Energy.

Appendix A: Correlated ﬁtting

When ﬁtting 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 ﬁts for the renor-malized masses and binding energies—while at diﬀerentbare input masses—used the same conﬁgurations for eachﬁt, thus introducing correlations between ﬁts. In orderto account for this, we did simultaneous correlated ﬁtsin both the x and y data sets using orthogonal distanceregression. Orthogonal distance regression attempts tominimize the radial distance between the ﬁt function andthe input data—as opposed to the vertical distance be-tween the ﬁt function and the y data that is typical innon-linear least squares ﬁtting.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 ﬁtting function, β k is thecollection of ﬁt 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 diﬀerent 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 deﬁ-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 ﬁt 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. Goroﬀ 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