Lattice $\phi^4$ Field Theory on Riemann Manifolds: Numerical Tests for the 2-d Ising CFT on $\mathbb{S}^2$

We present a method for defining a lattice realization of the $\phi^4$ quantum field theory on a simplicial complex in order to enable numerical computation on a general Riemann manifold. The procedure begins with adopting methods from traditional Regge Calculus (RC) and finite element methods (FEM) plus the addition of ultraviolet counter terms required to reach the renormalized field theory in the continuum limit. The construction is tested numerically for the two-dimensional $\phi^4$ scalar field theory on the Riemann two-sphere, $\mathbb{S}^2$, in comparison with the exact solutions to the two-dimensional Ising conformal field theory (CFT). Numerical results for the Binder cumulants (up to 12th order) and the two- and four-point correlation functions are in agreement with the exact $c = 1/2$ CFT solutions.


Introduction
Lattice Field Theory (LFT) has proven to be a powerful non-perturbative approach to quantum field theory [1]. However the lattice regulator has generally been restricted to flat Euclidean space, R d , discretized on hypercubic lattices with a uniform ultraviolet (UV) cut-off Λ U V = π/a in terms of the lattice spacing a. Here we propose a new approach to enable non-perturbative studies for a range of problems on curved Riemann manifolds. There are many applications that benefit from this. In the study of Conformal Field Theory (CFT), it is useful to make a Weyl transform from flat Euclidean space R d to a compact spherical manifold S d , or in Radial Quantization [2], a Weyl transformation to the cylindrical boundary, R × S d−1 , of AdS d+1 . Other applications that could benefit from extending LFT to curved manifolds include two-dimensional condensed matter systems such as graphene sheets [3], four-dimensional gauge theories for beyond the standard model (BSM) strong dynamics [4,5], and perhaps even quantum effects in a space-time near massive systems such as black holes.
Before attempting a non-perturbative lattice construction, one should ask if a particular renormalizable field theory in flat space is even perturbatively renormalizable on a general smooth Riemann manifold. Fortunately, this question was addressed with an avalanche of important research [6,7,8,9] in the 70s and 80s. A rough summary is that any UV complete field theory in flat space is also perturbatively renormalizable on any smooth Riemann manifold with diffeomorphism invariant counter terms corresponding to those in flat space [9]. Taking this as given, in spite of our limited focus on φ 4 theory, we hope this paper is the beginning of a more general non-perturbation lattice formulation for these UV complete theories on any smooth Euclidean Riemann manifolds.
The basic challenge in constructing a lattice on a sphere-or indeed on any nontrivial Riemann manifold-is the lack of an infinite sequence of finer lattices with uniformly decreasing lattice spacing [2,10,11]. For example, unlike the hypercubic lattice with toroidal boundary condition, the largest discrete subgroup of the isometries of a sphere is the icosahedron in 2d and the hexacosichoron, or 600 cell, in 3d. This greatly complicates constructing a suitable bare lattice Lagrangian that smoothly approaches the continuum limit of the renormalizable quantum field theory when the UV cut-off is removed. Here we propose a new formulation of LFT on a sequence of simplicial lattices converging to a general smooth Riemann manifold. Our strategy is to represent the geometry of the discrete simplicial manifold using Regge Calculus [12] (RC) and the matter fields using the Finite Element Method (FEM) [13] and Discrete Exterior Calculus (DEC) [14,15,16,17]. Together these methods define a lattice Lagrangian which we conjecture is convergent in the classical (or tree) approximation. However, the convergence fails at the quantum level due to ultraviolet divergences in the continuum limit. To remove this quantum obstruction, we compute counter terms that cancel the ultraviolet defect order by order in perturbation theory. We will refer to the resultant lattice construction as the Quantum Finite Element (QFE) method and give a first numerical test in 2d on S 2 at the Wilson-Fisher CFT fixed point.
While our current development of a QFE Lagrangian and numerical tests are carried out for the simple case of a 2-d scalar φ 4 theory projected on the Riemann sphere we attempt a more general framework. Since the map to the Riemann sphere, R 2 → S 2 , is a Weyl transform, the CFT is guaranteed [18] to be exactly equivalent to the c =1/2 Ising CFT in flat space and therefore presents a convenient and rigorous test of convergence to the continuum theory. More general examples can and will be pursued mapping conformal field theories in flat Euclidean space R d+1 either to R × S d , appropriate for radial quantization, with t = log(r) or to the sphere, S d , Current tests of the QFE method for the 3-d Ising model in radial quantization, R 3 → R × S 2 , are underway, so we take the opportunity here to give a brief introduction to both geometries. By considering the sphere S d as a dimensional reduction of the cylinder R × S d by taking the length of the cylinder to zero, both cases are conveniently presented together in Appendix A.
The organization of the article is as follows. In Sec. 2 we review the basic Regge Calculus/Finite Element Method framework as a discrete form of the exterior calculus and show its failure for a quantum field theory beyond the classical limit due to UV divergences. The reader is referred to the literature [13,14,15,16] for more details and to Ref. [19] for the extension to Dirac fermions. In Sec. 3 we address the central issue of counter terms in the interacting φ 4 theory required to restore the isometries on S 2 in the continuum limit. Sec. 4 compares our Monte Carlo simulation for fourth and sixth order Binder cumulants with the exactly solvable c = 1/2 CFT. In Sec. 5 we extend this analysis to the two-point and four-point correlation functions. We fit the operator product expansion (OPE) as a test case of how to extract the central charge, OPE couplings, and operator dimensions.

Classical Limit for Simplicial Lattice Field Theory
The scalar φ 4 theory provides the simplest example. On a smooth Riemann manifold, (M, g), the action, is manifestly invariant under diffeomorphisms: x = f (x) ≡ x (x). We include the coupling to the scalar Ricci curvature ( ξ 0 Ric), where ξ 0 = (d−2)/(4(d−1)) and an external constant (scalar) field h. Ric = (d − 1)(d − 2)/R 2 on the sphere S d of radius R. The field, φ(x), is an absolute scalar or in the language of differential calculus a 0-form with a fixed value at each point P in the manifold, independent of co-ordinate system: φ (x ) = φ(x(x )) for x = f (x) . For future reference to the discussion in Sec. 2.2 on the discrete exterior calculus, we identify vol d = √ g d d x = √ gdx 1 ∧ · · · ∧ dx d , as the volume d-form. The scalar φ 4 theory in d = 2, 3 has a (super-)renormalizable UV weak coupling fixed point and a strong coupling Wilson-Fisher conformal fixed point in the infrared (IR), as illustrated in the phase diagram in Fig. 2.1. In passing, we also note its similarity to 4-d Yang Mills theory with a sufficient number of massless fermions to be in the conformal window [5], which is a central motivation for this research. Both theories are UV complete, perturbatively renormalizable at weak coupling, and have a strong coupling conformal fixed point in the IR. The mass terms (m 2 0 φ 2 or m 0ψ ψ, respectively) must be tuned either to or near the critical surface to reach the conformal or mass deformed theory.
To formulate a lattice action, we introduce a simplicial complex, or triangulation in 2d, as illustrated in Fig. 2.2, and a discrete action, where i labels all vertices and the sum ij runs over all links with proper length l ij . The technical requirement is to fix the weights, (K ij , √ g i , m 2 i , λ i ) as functions of bare couplings and the target manifold, on a sequence of increasingly fine tessellations so that in the limit of vanishing lattice spacing, a = O(l ij ), the quantum path integral converges to the continuum renormalized quantum theory with a minimal set of fine tuning parameters. For the φ 4 theory, there is one parameter to tune in the approach to the critical surface: the relevant mass parameter Figure 2.2: A 2-d simplicial complex with points (σ 0 ), edges (σ 1 ) and triangles (σ 2 ) (illustrated in yellow). At each vertex σ 0 there is a dual polytope, σ * 0 (illustrated in red), and at each link, σ 1 , there is a dual link σ * 1 and its associated hybrid cell σ 1 ∧σ * 1 (illustrated in blue).
The construction of our QFE simplicial lattice action can be broken into three steps: i. Simplicial Geometry: The smooth Riemann manifold, (M, g), is replaced by simplex complex, (M σ , g σ ) with piecewise flat cells (see Sec. 2.1). ii. Classical Lattice Action: The continuum field φ(x) is replaced by a discrete sum over FEM basis elements, φ(x) → W i (x)φ i (see Sec. 2.2). iii. QFE: Quantum Corrections: Quantum counter terms are added to the discrete lattice action to cancel defects at the UV cut-off (see Sec. 3.2).

Geometry and Regge Calculus
The Regge Calculus approach to constructing a discrete approximation to a Riemann manifold, (M, g), proceeds as follows. First the manifold, M, is replaced by a simplicial complex, M σ , composed of elementary simplicies: triangles in 2d as illustrated in Fig. 2.2, tetrahedrons in 3d, etc. This graph defines the topology of the manifold in the language of Category Theory [20]. Next a discrete metric is introduced as a set of edge lengths on the graph: g µν (x) → g σ = {l ij }. Assuming a piecewise flat interpolation into the interior of each simplex, we now have the Regge representation of a Riemann manifold (M σ , g σ ) that is continuous but not differentiable. The curvature is given by a singular distribution at the boundary of the each simplex; in 2d, concentrated at the vertices, and in higher dimensions at the d−2-dimensional "hinges." The integrated curvature over the defect is easily computed by parallel transport of the tangent vector around each defect.
Each simplex is parameterized by using d + 1 local barycentric coordinates, 0 ≤ ξ i ≤ 1. Using the constraint, ξ 0 +ξ 1 +· · ·+ξ d = 1 to eliminate ξ 0 , our φ 4 action in Eq. (2.1) on this piecewise flat Regge manifold is given as a sum over each simplex, We may choose an isometric embedding into a sufficiently high dimensional flat Euclidean space with vertices at y = r n so a point in the interior of each simplex and its induced metric are given by This defines the volume element, √ G σ = det[G ij ] and inverse metric, G ij σ , as well.
Since we are not considering dynamical gravity, this piecewise flat metric is frozen (i.e. quenched) and chosen to conform as closely as possible to our target manifold for the quantum field theory. For example this can be achieved by computing the edge length l ij , from the proper distance between points on the original manifold x i , x j or from the Euclidean distance | r i − r j | in the isometric embedding space, r i = y(x i ) to first order relative to the local curvature of the target manifold. At this stage, the dynamical quantum field φ(x) is still a continuum function on the piecewise flat Regge manifold.

Hilbert Space and Discrete Exterior Calculus
The second step is the approximation of the matter field φ(x) as an expansion, 12 is the area of the triangle formed by the sites 1, 2, and the circumcenter, σ * 2 (123). The free scalar action on the entire simplicial complex is now found by summing over triangles, (2.7) Each link ij receives two contributions-one from each triangle that it borders-resulting in the total area for the hybrid cell, For higher dimensions, a natural generalization of the kinetic term is where V ij = |σ 1 (ij) ∧ σ * 1 (ij)| = l ij S ij /d is the product of the length of the link (l ij ) times the volume of the d−1-dimensional "surface", S ij = |σ * 1 (ij)|, of the dual polytope normal to the link i, j . A mass term has also been included weighted by the dual lattice volume √ g i = |σ * 0 (i)|. This elegant form was recommended in the seminal papers on random lattices by Christ, Friedberg and Lee [21,22,23] and subsequently in the FEM literature by the application of the simplicial Stokes' theorem for Discrete Exterior Calculus (DEC). However, we note this generalization is not equivalent to linear FEM for d > 2 (See Ref. [19]).
To appreciate the DEC approach [14,15,16,24], it is useful to expand a little on the geometry of a simplicial complex, S, and its Voronorï dual, S * . A pure simplicial complex S consists of a set of d-dimensional simplices (designated by σ d ) "glued" together at shared faces (boundaries) consisting of d−1-dimensional simplices (σ d−1 ), iteratively giving a sequence: σ d → σ d−1 → · · · σ 1 → σ 0 . This hierarchy is specified by the boundary operator, where i k means to exclude this site. Each simplex σ n (i 0 i 1 · · · i n ) is an anti-symmetric function of its arguments. The signs in Eq. (2.9) keep track of the orientation of each simplex. It is trivial to check that the boundary operator is closed: ∂ 2 σ n = 0. On a finite simplicial lattice ∂ is a matrix and its transpose, ∂ T , is the co-boundary operator. The circumcenter Voronorï dual lattice, S * , is composed of polytopes, σ * 0 ← σ * 1 ← · · · ← σ * d , where σ * n has dimension d − n as illustrated in Fig. 2.2. A crucial property of this circumcenter duality is orthogonality. Each simplicial element σ n ∈ S is orthogonal to its dual polytope σ * n ∈ S * . As a consequence, the volume, of the hybrid cell, σ n ∧ σ * n , is a simple product. Hybrid cells, constructed from simplices σ n in S and their orthogonal dual σ * n in S * , give a proper tiling of the discrete d-dimensional manifold with the special case |σ 0 (i)| = 1 and |σ * 0 (i)| = √ g i . This is a first modest step into discrete homology and De Rham cohomology on a simplicial complex.
The discrete analogue of differential forms, ω k , is a pairing or map, to the numerical value of a field from sites in the case of 0-forms, from links in the case of 1-forms, and from k-simplices in the case of k-forms. Of course this is familiar to lattice field theory, associating scalars (ω 0 ∼ φ x ), gauge fields (ω 1 ∼ A µ dx µ ) and field strengths (ω 2 ∼ F µν dx µ ∧ dx ν ) with sites, links and plaquettes respectively. This enables us to define the discrete analogue of the exterior derivative dω k of a k-form by replacing the continuum Stokes' theorem by the discrete map or DEC Stokes' theorem: The discrete exterior derivative is automatically closed (dd = 0) because the boundary operator is closed (∂∂ = 0). Applying Eq. (2.12) to the discrete exterior derivative of a scalar (0-form) field trivially gives, the standard finite difference approximation on each link.
Next we need to define the discrete analogue of the Hodge star ( * ). This is the first time an explicit dependence on the metric is introduced. In the continuum, a k-form is an anti-symmetric tensor, ω k (x) = (k!) −1 ω µ 1 ,µ 2 ···µ k dx µ 1 ∧dx µ 2 ∧· · ·∧dx µ k , in a orthogonal basis of 1-forms dual to tangent vectors: dx µ (∂ ν ) = δ µ ν . The Hodge star takes a k-dimensional basis into its orthogonal complement, where µ 1 , µ 2 , · · · , µ n is an even permutation of (1, 2, · · · , n). However the wedge product (or equivalently the Levi-Civita symbol) is not a tensor but a weight 1 tensor density [25]. The true volume k-form (or tensor), vol k = √ g k dx µ 1 ∧ · · · ∧ dx µ k , requires a factor of √ g k which under the Hodge star operation in Eq. (2.14) gives the identity, √ g n−k * (vol k ) = √ g k vol n−k , where we have used orthogonality to factor √ g = √ g k √ g d−k , that is, between the plane and its dual. Consequently on the simplicial complex, it is reasonable that the proper definition of the discrete Hodge star, replaces these factors √ g k , √ g n−k by finite volumes |σ k |, |σ * k | respectively. The Hodge star identity in Eq. (2.15) uniquely fixes the dual field values, ω * k , σ * k and the action of a discrete co-differential, δ = * d * through Stokes theorem in Eq. (2.12) on S * .
Putting this all together, we consider the DEC Laplace-Beltrami operator on scalar fields, (δ + d) 2 which corresponds to the action in Eq. (2.8). For 2d this is illustrated in Fig. 2.3, as the sum of fluxes through the boundaries ∂σ * 0 (i) with surface area, S ij /(d − 1)! = V ij /l ij = |σ 1 (ij) ∧ σ * 1 (ij)|/l ij . This is identical to linear finite elements in 2d. Besides providing an alternate approach to linear finite elements for constructing the discrete Hilbert space for scalar fields in d > 2 dimensions, it provides a useful geometric framework for fields with spin. We note however that additional considerations are still needed for non-Abelian gauge fields [22], Dirac Fermions [19] and Chern-Simons terms [26]. The best geometrization of simplicial field theories and the error estimates thereof are an active research topic. To complete the simplicial action, we need to add the potential term. This may be constructed in a variety of ways. If one follows strictly the linear finite element prescription, the expression for the local polynomial potential (e.g. mass and quartic terms) will not be local but rather point split on each simplex σ d . For example, in 2d, after expanding φ(ξ) = d i=0 ξ i φ i and evaluating the integral over the linear elements, the contribution for the quadratic term on a single triangle, σ 2 (123), is The general expression for a homogeneous polynomial over a d-simplex with volume V d = |σ d | is given as as sum over distinct partitions of n: Nonetheless in the spirit of dropping higher dimensional operators in a derivative expansion, we choose local terms approximating the potential at each vertex weighted by the volume , √ g i = |σ * 0 (i)|, of the dual simplex σ * 0 (i).
This approximation leads to our complete simplicial action, combining the DEC Laplace-Beltrami operator for the kinetic term and the local approximation for the bound- ary, where µ 2 0 = −m 2 0 /2. We will use this action for our discussion of S 2 .
Finally we should acknowledge that there are many other alternatives in the FEM literature worthy of consideration which may offer improved convergence and faster restoration of the continuum symmetries. Our goal here is to find the simplest discrete action on simplicial lattice capable of reaching the correct continuum theory with no more fine tuning than is required on the hypercubic lattice.

Spectral Fidelity on S 2
We represent S 2 embedded in R 3 parameterized by 3-d unit vectors: In order to preserve the largest available discrete subgroup of O(3), we start with the regular icosahedron illustrated on the left in Fig. 2.4 and subsequently divide each of the 20 equilateral triangles into L 2 smaller equilateral triangles. Then we project the vertices radially outwards onto the surface of the circumscribing sphere, dilating and distorting each triangle from its equilateral form, but preserving exactly the icosahedral symmetries. The total numbers of faces, edges, and vertices are F = 20L 2 , E = 30L 2 and N = 2+10L 2 , respectively, satisfying the Euler identity F − E + N = 2.
The images of the vertices on the sphere are then connected by new links consistent with the unique Delaunay triangulation [27]. In the triangulation, each vertex is connected to five or six neighboring vertices by edges x, y . The lengths are set to the secant lengths l xy = |r x −r y | in the embedding space between the vertices on the sphere, which approximates the geodesic length to O(l 2 xy ). The extension to a lattice for the R × S 2 cylinder introduces a uniform (periodic) lattice perpendicular to the spheres at t = 0, 1, · · · , L t − 1. The sequence of refinements as L → ∞ divides the total curvature into vanishingly small defects at each vertex. As an alternative formulation, we could replace the flat triangles with spherical triangles [28], introducing spherical areas, geodesic lengths, moving the curvature uniformly into the interior of each triangle. In this formulation, each triangle would have a vanishing deficit angle in the continuum limit. Both discretizations are equivalent at O(a 2 ), and as such we prefer flat triangles due to their relative simplicity.
The first test of our construction is to look at the spectrum of the free theory, where x, y = 1, ..., N enumerates the lattice sites on a finite simplicial lattice. The spectrum is computed from the generalized eigenvalue condition (2.22) where the eigenvalues at zero mass are E (0) n . Each distinct eigenvalue E n has right/left generalized eigenvectors, φ n (x)/φ * n (x), with the orthogonality and completeness relations, The spectral decomposition for the free propagator is Alternatively we may define a Hermitian form, M = g −1/4 M g −1/4 , with complex conjugate right and left eigenvectors, g 1/4 x φ n (x) = x|n and φ * n (x) = n|x , respectively. In Dirac notation the completeness and orthogonality are given by 1 = x |x x| , 1 = n |n n| and δ x,y = x|y , δ n,m = n|m respectively.
Returning to the example of the FEM sphere, S 2 , we have verified that the generalized eigenvalues which lie well below the cut-off are well fitted by the continuum spectrum, E lm = l(l + 1), with the 2l + 1 degeneracy m = −l, · · · l. Indeed, any finite eigenvalue approaches its continuum value as 1/L 2 in the limit of infinite refinement L → ∞. In addition, we note that the right eigenvectors are well approximated by the continuum spherical harmonics, Y lm (r x ), evaluated at the lattice sitesr x . This is illustrated in Fig. 2.5, where the eigenvalues are estimated by computing diagonal matrix elements, against l. On the left of Fig. 2.5, the lack of 2l + 1 degenerate multiplets on a coarse lattice with L = 8, as we approach the cut-off at large l ∼ O(L), shows the breakdown of rotational symmetry. On the right of Fig.2.5, the degeneracy of the spherical representation holds to high accuracy and the average of the 2l + 1 eigenvalues at each l level gives the correct continuum dispersion relation, l(l + 1), to O(10 −5 ) for l L = 128.
We will refer to the exact convergence of fixed eigenvalues and their associated eigenvectors to the continuum as the cut-off is removed as spectral fidelity. This is a theoretical consequence of FEM convergence theorems for shape regular linear elements as the diameter goes uniformly to zero [13]. We do not provide a proof but will assume that this property holds for our implementation if we apply the DEC [15] to the Laplace-Beltrami operator.
It is useful to compare the low spectrum of our FEM operator on S 2 , shown in Fig. 2.5, to the corresponding spectrum on the hypercubic refinement of the torus T 2 , shown in Fig. 2.6. For a finite L d toroidal lattice, the exact hypercubic spectra is given by where the discrete eigenvalues are enumerated by k µ = 2πn µ /L for integer n µ ∈ [−L/2, L/2− 1]. The eigenvectors can be found by a Fourier analysis, and are given by The zero mass spectrum E n = 4(sin 2 (k x /2) + sin 2 (k y /2)), k µ = nµπ L , n µ ∈ [−L/2, L/2 − 1] for the 2-d lattice Laplacian on a regular 32 × 32 square lattice plotted against n 2 x + n 2 y .
Converting to dimensionful variables (m, p) and holding a physical mass fixed (m ≡ m 0 /a), the dispersion relation is The Lorentz breaking term vanishes as a 2 = O(1/L 2 ). Again spectral fidelity holds for any fixed spectral value E n as the lattice converges to the continuum. However, as we approach the cut-off, it is increasingly distorted.
Clearly the spectral fidelity on the hypercubic torus is comparable to the FEM spectra on S 2 . This is not surprising; indeed the square lattice can also be viewed as a FEM realization. One simply divides each square into two right angle triangles and notices that the formula in Eq. (2.8) implies a zero contribution on the diagonal links. This generalizes to higher dimensional hypercubic lattices when using the DEC form. Of course, the major difference between the square lattice on R 2 and the simplicial sphere S 2 is the former breaks one rotational isometry of the continuum but conserves two discrete subgroups of translations, whereas the later breaks all three isometries of O(3) down to a fixed finite subgroup independent of the refinement.

Obstruction to Non-Linear Quantum Path Integral
Having demonstrated the spectral fidelity of our construction on S 2 at the Gaussian level, we turn next to the more difficult problem of an interacting quantum theory, starting with the FEM action given in Eq. (2.19). We have performed extensive Monte Carlo simulations for the path integral given the FEM action in Eq.  . Instead we see that there is an alarming instability at large L = O(100), which we believe is due to the lack of well-defined critical surface with a second order phase transition required to define a continuum limit. A more graphic indication of this failure is present when we consider the average value of φ 2 i as a function of position on the sphere as shown on the righthand side of Fig. 2.7. Due to a spatial variation in the UV cut-off, this is not uniform, with variation that can be seen by comparing the regions close to and far form the 12 exceptional five-fold vertices of the original icosahedron. As we will show, this failure is also evident in a lattice perturbation expansion. At small λ 0 a spherically asymmetric contribution to φ 2 i is given by a UV divergent one-loop diagram.
In conclusion, in the application of FEM to quantum field theory, we have encountered a fundamentally new problem. The FEM methodology has been developed to give a discretization for non-linear PDEs that converges to the correct continuum solution as the simplicial complex is refined. In the context of a Lagrangian system for a quantum field, this implies a properly implemented FEM should therefore guarantee convergence to all smooth classical solutions of the EOM as the cut-off is removed. However, quantum field theory is more demanding. The path integral for an interacting quantum field theory is sensitive to arbitrarily large fluctuations-even in perturbation theory-on all distance scales, down to the lattice spacing, due to ultraviolet divergences. This amplifies local UV cut-off effects on our FEM simplicial action.
Our solution is to introduce a new Quantum Finite Element (QFE) lattice action that includes explicit counter terms to regain the correct renormalized perturbation theory. We conjecture that, for any ultraviolet complete theory, if the QFE lattice Lagrangian is proven to converge to the continuum UV theory to all orders in perturbation theory, this is sufficient to define its non-perturbative extension to the IR. We believe this is a plausible conjecture consistent with our experience on hypercubic lattices for field theories in R d , but it is far from obvious. It needs careful theoretical and numerical support to determine its validity. In Sec. 4 and Sec. 5 we give extensive numerical test for φ 4 on S 2 on the critical surface. In particular the new Monte Carlo simulation of the Binder cumulant U 4 on the right hand side of Fig. 4.1 including the quantum counter term gives a critical value of the Binder cumulant, U 4,cr = 0.85020(58)(90), in agreement with the continuum value, U * 4 = 0.8510, to about one part in 10 3 . While this is promising, the limitations due to statistics and the restriction of our studies to the simplest scalar 2-d CFT is duly acknowledged.

Ultraviolet Counter Terms on the Simplicial Lattice
To remove the quantum obstruction to criticality for the FEM simplicial action, we begin by asking if we can add counter terms to the action Fig. 2.19 to reproduce the renormalized perturbation expansion order by order in the continuum limit at the UV weak coupling fixed point. On a hypercubic lattice, it has been proven for φ 4 theory [29] that this can be achieved by taking the lattice UV cut-off to infinity, holding the renormalized mass and coupling fixed. Here we will suggest how this can be achieved on a simplicial lattice for 2d and 3d. While we do not attempt a proof, the similarity with the hypercubic example strongly suggests that a proof could be found at the expense of increased technical difficulty.
The φ 4 theory is super-renormalizable in 2d and 3d, with one-loop and two-loop divergent diagrams, respectively, only contributing to the two-point function. The oneloop diagram is logarithmically divergent in 2d and linearly divergent in 3d. The two-loop diagram is UV finite in 2d, but logarithmically divergent in 3d. The divergent contributions renormalize the mass via the one-particle irreducible (1PI) contribution to the self-energy, as illustrated in the first two panels in Fig. 3.1. In both 2d and 3d the three-loop diagram, and all higher-order diagrams, are UV finite. In 4d there is a logarithmically divergent contribution to the four-point function which contributes to the quartic coupling λ 0 ; however the complete non-perturbative theory is the trivial free theory with no interesting IR physics.

Lattice Perturbation Expansion
The perturbation expansion for the φ 4 theory on our FEM lattice (or indeed any lattice) starts with the partition function, by expanding in the quartic term. In our FEM representation, the action is Beltrami operator (K i,j ) and the bare mass (m 0 ).
Following the standard Feynman rules for a perturbative expansion in λ 0 , we compute the propagator, φ i φ j = Z −1 Dφ φ i φ j e −S , to second order: After amputating the external lines, we find the inverse propagator is the 1PI simplicial self-energy.

One Loop Counter Term
Since there is no analytic spectral representation of the free FEM Green's function, we compute the one loop diagram in co-ordinate space by numerical evaluation of the propagator. To be concrete, for 2d on S 2 and for 3d on R × S 2 , the Gaussian term is where the sites are now labeled by i = (x, t). The 2-d geometry, S 2 , can be viewed as a special case with a single sphere at t = 0. The integer t labels each sphere along the cylinder, and x indexes the sites on each sphere. K x,y is non-zero on only nearest neighbor links x, y . To take the continuum limit we need to define the normalization convention of our lattice constants: where in this specific context E refers to the number of edges on the simplicial graph. The extra factor of 2/3 is introduced to compensate, on average, for the six nearest neighbors per site on the sphere relative to a conventional square lattice with four nearest neighbors.
In flat 2-d space this factor of 2/3 for the triangular lattice gives the same dispersion relation, E = m 2 0 + k 2 + O(k 4 ), as the square lattice.
In lattice perturbation theory on a 2-d FEM simplicial lattice, we expect the logarithmically divergent one-loop term to give a site-dependent mass shift, The numerical computation of the one-loop diagram on S 2 does indeed have the appropriate form, The logarithmic divergence is regulated by the IR mass, m 2 . Of course on the lattice, at fixed dimensionless bare mass, m 2 0 = a 2 m 2 and bare coupling, λ 0 , there is no actual UV divergence. The renormalized perturbation expansion is defined by fixing the physical mass (m), which in the continuum limit (N → ∞) corresponds to a vanishing effective lattice spacing, a 2 m 2 = O(1/N ). On a hypercubic lattice, there is a universal cut-off (π/a x = π/a), whereas on a simplicial lattice the one-loop term separates into two terms on the right side of Eq. We have checked numerically, to high accuracy, that the first divergent term is independent of position. This is a crucial observation which we believe is a general consequence of the FEM prescription. Spatial dependence for a divergent term would imply the need for a divergent counter term to restore spherical symmetry. This departs from the standard cut-off procedure on regular lattices and would be very difficult if not impossible to implement. Moreover, we show in Fig. 3.2 that the fit coefficient is accurately determined to be √ 3/(8π), which is the exact continuum value.
A heuristic argument for the value of this constant follows from our S 2 tessellation as a nearly regular triangular lattice. In flat space, relative to a square lattice, the density of states on a triangular lattice is . The extra factor of √ 3/2 is due to the area |σ * 0 | of the dual lattice hexagon relative to the dual lattice square on a hypercubic lattice. However, since this constant defines the local charge in the 2-d Coulomb law and our linear elements do not implement local flux conservation, this is surprising result. Nonetheless in Sec. 3.3 we will prove this based on the FEM principle of spectral fidelity and the renormalization group for the mass anomalous dimension. To the best of our knowledge this is a novel extension of FEM convergence theorems.
Next consider the finite spatial dependence term in Eq. (3.9). As illustrated in Fig. 3.3, it is almost exactly a linear function of log[ Again there is a heuristic explanation for this. An almost perfect analytic expression can be found by considering the dilatation incurred by the projection of our triangular tessellation of the icosahedron onto the sphere. The icosahedron is a manifold with a flat metric except for 12 conical singularities at the vertices. Our choice of the simplicial complex began with flat equilateral triangles on each of the 20 faces of the original icosahedron. The radial projection onto the sphere is a Weyl transformation to constant curvature with conformal factor (or Jacobian of the map), where R c is the circumradius for one of the 20 icosahedral faces and x 1 , x 2 the flat coordinates on that face. This formula gives a line coinciding with the linear data in Fig. 3.3. Although this is an extremely good approximation to the counter term, it is not perfect since it neglects the conical singularities in the map near the 12 exceptional vertices of the original icosahedron. These exceptional points are visible in the upper left corner of Fig. 3.3. We find that the true numerical computation does give better results and, moreover, it is a general method not requiring our careful triangulation procedure. A more irregular procedure should also be allowed. The lesson is that the simplicial metric, g ij = {l ij }, on a given Regge manifold not only defines the intrinsic geometry but also a co-ordinate system breaking diffeomorphism explicitly. The set of lengths includes both a definition of curvature and the discrete co-ordinate system. The counter term must compensate for this arbitrary choice.
The QFE one loop counter term is introduced to cancel the finite position dependence in Eq. (3.9) that violates rotational invariance to order λ 0 . To project out the spherical component of a local scalar density, ρ x , we average over the rotation group R ∈ SO(3), Applied to ρ x = [M −1 ] xx , the subtracted lattice Green's function is This removes completely the logarithmic divergence, leaving the position dependent finite counter term which adds a contribution to the FEM action, In 2d φ 4 theory this is the only UV divergent term. We also computed the two-loop contribution to the 1PI simplicial self energy Eq. (3.5) which in co-ordinate space is given by the third power of the Green's function, [M −1 ] 3 xy , between the vertex at x and y. The position dependence, is found again by subtracting the average contribution on the sphere. The result is plotted in Fig. 3.4, against log[ √ g x ] which demonstrates that it vanishes rapidly in the continuum limit as 1/L 2 ∼ 1/N consistent with naïve power counting in continuum perturbation theory. While we do not attempt a general analysis of the simplicial lattice perturbation expansion, it is plausible based on analogous studies on the hypercubic lattice [30,31] that after canceling cut-off dependence in the one loop UV divergent term the entire expansion restores the renormalized perturbation series in the continuum limit. However our test on the 2-d Riemann sphere simplicial lattice goes further. Namely we give numerical evidence that this 2-d QFE lattice action, when tuned to the the critical surface, approaches the non-perturbative continuum theory in the IR for the Ising CFT field theory.

Universal Logarithmic Divergence
The success of our QFE action prescription depends on the coefficient of the logarithmic divergence being exact in the continuum limit. On S 2 or any maximally symmetric space, this implies the coefficient is also independent of position. At first, this observation appears surprising. After all, the UV divergence is sensitive to the local, short distance cut-off which is not uniform. We claim this is a consequence of the observed spectral fidelity of the DEC Laplace-Beltrami operator and renormalization group (RG) for a logarithmic divergence relating UV divergences to the IR regulator. Adapting the argument to a general 2-d Riemann manifold is in principle straight forward, matching the local divergence to the one loop renormalization of the continuum theory [7] at each point x.
Let us first apply this renormalization group (RG) argument to the hypercubic lattice, where we already know and understand the answer. Suppose that the one-loop diagram has a logarithmic term, with an unknown position dependent coefficient c x . To isolate this coefficient, we take the logarithmic derivative, finding (3.16) so that γ 1 (m 2 ) = 2c x + O(a 2 m 2 ). This is the one-loop contribution to the anomalous dimension of the φ 2 operator. (To clarify scaling, we have re-introduced the lattice spacing a so that p = k/a, m = m 0 /a have mass dimensions.) By power counting, the integral resulting from the logarithmic derivative is UV finite in the continuum limit. We may now introduce a new cut-off, Λ 0 , separating the IR from the UV: m Λ 0 π/a. To estimate the integral, we can take the continuum limit π/a → ∞, followed by taking the IR regulator to zero: m 2 /Λ 2 0 → 0. This proves that the coefficient c x = 1/(4π) is independent of position and identical to the continuum value on R 2 . Of course, this is not surprising for flat space. Now we can apply the same reasoning to the simplicial sphere S 2 , beginning with the spectral representation of the Green's function, On the basis of FEM spectral fidelity, assume the low eigenspectrum for l ≤ L 0 is well approximated to O(1/N ) by spherical harmonics, where R/a is the radius of the sphere in units of the lattice spacing. Matching the area of the sphere with N triangles we have 4πR 2 = a 2 N √ 3/2. For R = 1 the spectrum is l(l + 1), but with our convention of setting a = 1, we find R = √ 3N/(8π). Next, we fix the eigenvector normalization by requiring x √ g x |Y 00 | 2 = N/4π.
By splitting the spectral sum at the cut-off L 0 , and using the addition formula, 4π m Y * lm (r x )Y lm ( r y ) = (2l + 1)P l (r x · r y ) for l < L 0 , we have where we introduced µ 2 = R 2 (m 2 0 /a 2 ). The first term is indeed rotationally invariant. At x = y, the asymptotic limit of the first term is √ 3 log(L 0 )/(8π), which is the desired behavior. As with the infinite, flat lattice, we wish to convert the sum of the first term to an integral in the limit N → ∞: with Λ 2 0 = a 2 L 0 (L 0 + 1). Again, if we take a logarithmic derivative, we have a UV convergent integral that is saturated by rotationally invariant contributions up to power corrections at O(m 2 0 /Λ 2 0 ), (3.23) In the limit O(m 2 0 /Λ 2 0 ) → 0, we isolate the coefficient in Eq. (3.15) and prove that the divergence is both independent of position on the sphere with its coefficient identical to the continuum one loop diagram on the sphere. The essential assumption we have made is the spectral fidelity property of the discrete DEC Laplace-Beltrami operator. In Appendix A we provide the exact continuum propagator on the sphere, and identify the coefficient of the logarithmic divergence by relating the angle to the lattice cut-off by θ 2 ∼ 1/L 2

Comments on Structure of the Counter Terms
In spite of our very limited example of the QFE counter terms on S 2 , we hazard a general interpretation which we hope will guide extensions to any smooth Riemann manifold. This is based in part on current study of counter terms in 3-d φ 4 theory on R×S 2 . Nonetheless, we acknowledge the fact that our examples so far may have very special features not shared more generally.
The first special feature is that the S 2 manifold is an example of a maximally symmetric manifold. This property is shared by flat space, spheres, and anti-de Sitter space with constant curvature. On a d-dimensional maximally symmetric manifold, there is a full set of d(d+1)/2 isometries implying that all points have identical metric properties. A consequence is in these cases, the counter terms can be easily defined by projecting out the average over the full set of isometries as in Eq. (3.11) above. While these maximally symmetric manifolds are actually of paramount interest, we are also confident that this restriction is not necessary. On a general manifold, the procedure would be to compute the continuum UV divergence at each site and subtract it before defining the counter term.
The second special feature of the one-loop diagram in 2d in φ 4 theory is that it is both local on the lattice and logarithmically divergent in the continuum. In 3d the situation changes. First the one-loop term in the continuum is linearly divergent. However on the lattice such a divergence is not manifest because it enters as a dimensional factor of 1/a, which in the bare lattice action is chosen implicitly to be 1/a = O(1). Also all UV power divergent diagrams are finite in the IR so that there is no need for a dimensionless mass regulator m 0 = am. We have computed the one-loop diagram on the simplicial cylinder R × S 2 , and found numerically that the counter term is almost identical to the 2-d case. This is easy to understand. The spectral decomposition of the 3-d Green's function, is a natural generalization to the case on S 2 given in Eq. (3.20). We see that the breaking of rotational invariance on R × S 2 is very similar to the case of S 2 for each mode ω n and as such the overall breaking term will also be very similar to the 2-d case.
Next, in 3d, there is a new, logarithmically divergent two-loop term. Again, a comparison between 2d and 3d is interesting. The one-loop term in Hamiltonian form can be viewed simply as a problem of normal ordering: φ 4 (x) =: φ 4 (x) : + φ 2 (x) φ 2 (x) 0 . The diagram has no unitarity cuts and can be simply dropped when defining a normal ordered perturbation expansion. In 3d, the two-loop term introduces a non-local contribution with essential unitarity cuts that must be preserved in a renormalized perturbation theory. On the lattice, the continuum position space two-loop diagram is replaced by the element-wise third power of the Green's function, Again, we must subtract the rotationally invariant contribution at fixed t 1 − t 2 and |r x − r y | 2 = 2 − 2 cos(θ xy ). We have found numerically that although the rotationally invariant part has a power fall-off dictated by naïve dimensionality, the residual breaking term, δµ 2 xy,t 2 −t 1 φ x,t 1 φ y,t 2 , falls off exponentially in lattice units. Consequently, in physical units, it is exponential in (t 2 −t 1 )/a and |r x −r y |/a and therefore to leading order may be replaced by a local counter term in the QFE action. We postpone further analysis of counter terms to future publications.
We are also considering alternative methods that avoid the difficulty of computing individual UV divergent perturbative diagram. For example, as a proof of concept, we have implemented the Pauli-Villars (or Feynman-Stuekelberg) approach. It introduces an intermediate scale of Pauli-Villars mass, M P V π/a separating the IR from the UV and protects the UV diagram from the dependence of the variable cut-off of the simplicial lattice. This separation of scales is another way to exploit the spectral fidelity of the low spectrum, but this time to all orders in perturbation theory. The implementation amounts to the addition of a ghost field giving a simplicial equivalent of the continuum propagator, (3.27) Now, reaching the continuum requires a double limit: a → 0 at fixed aM P V , followed by aM P V → ∞. We have implemented this on S 2 by adding a quadratic PV term, −M −2 P V φ x K xz K zy φ y / √ g z , to the FEM action in Eq. (2.19), seeing qualitatively that this does work. However, in addition to the cost of the double limit, in the context of our current φ 4 simulations, the PV term has the technical disadvantage that it prevents the use of the very efficient cluster algorithm of Ref. [32].
We are also exploring extensions of the renormalization group approach in our demonstration of the one-loop logarithmic divergences in 2d in Sec. 3.3. For asymptotically free theories in 4d, such as the non-Abelian gauge theory, the Lagrangian, F 2 /g 2 , has a dimensionless coupling with a logarithmic divergence. We anticipate the use of the RG approach and perhaps the scale setting properties of Wilson flow [33] to correct the scheme dependence of a simplicial lattice. As demonstrated in the classic paper by A. Hasenfratz and P. Hasenfratz [34] for pure non-Abelian gauge theory, the lattice scheme dependence is a one-loop effect which on a simplicial lattice should be replaced by a local site dependent multiplicative counter term for F 2 .

Numerical Tests of UV Competition on S 2
Here we present our first test of our QFE simplicial lattice construction for the nonperturbative study of quantum field theories on curved manifolds. The Monte Carlo simulation of the 2-d scalar φ 4 theory on the Riemann sphere, S 2 , must agree within statistical and systematic uncertainties with the exact solution of the Ising or c = 1/2 minimal CFT in the continuum limit. The first test is to fine-tune the mass parameter to the critical surface illustrated in Fig. 2.1 and to compare the Monte Carlo computation of bulk Binder cumulants with analytic values. In Sec. 5, we extend our tests to the detailed form of the conformal two-point and four-point correlators.

Let us begin by defining the Euclidean correlations functions
for our scalar field theory as a derivative expansion of the partition function, . Likewise replacing the current by a constant magnetic field, J(x) = h g(x), the derivatives of the partition function give the magnetic moments, . From these moments, defining homogeneous quotients, Q 2n = M 2n / M 2 n , the first three Binder cumulants are [35] These are just the connected moment expansion of the free energy, divided by appropriate factor of M 2n . Therefore they vanish in the Gaussian limit. The overall normalization has been chosen so that the Binder cumulants are unity in the ordered phase. We designate the exact value of Binder cumulants for the critical Ising CFT by U * 2n .
On the sphere S 2 the exact value of the Binder cumulants U * 2n can be computed from the solution for the conformal n-point functions on R 2 . This a consequence of the special property for conformal correlators under a Weyl transformation of the flat metric: g µν (x) = Ω 2 (x)δ µν on S 2 . For example the CFT correlators of a primary operator φ with dimension ∆ obey the identity [18] φ(x 1 ) · · · φ(x n ) gµν = 1 Ω(x 1 ) ∆ · · · 1 Ω(x n ) ∆ φ(x 1 ) · · · φ(x n ) f lat . (4.5) As also pointed out in Ref [18], even the Weyl anomalies will cancel in homogeneous ratios of CFT correlators. Consequently as noted in Ref. [36] homogeneous moments in the Binder cumulants are computed by integration over the n-point function on the compact manifold.
In our current application, the Weyl transformation is a stereographic projection from R 2 to the Riemann sphere, S 2 , which can be explicitly parameterized by w = (r x + ir y )/(1 + r z ) = tan(θ/2)e iφ ,r = (sin θ, cos φ, sin θ sin φ, cos θ) , (4.6) where w is the complex co-ordinate w = x + iy in the R 2 plane andr is unit vector in R 3 in Eq. (2.20). The resulting metric on S 2 is with Ω 2 (θ) = cos 2 (θ/2). After the Weyl rescaling the two-point function on S 2 is where θ 12 is the angle between the two radial vectors,r 1 ,r 2 , on the surface of the sphere. Incidentally a pedestrian proof uses standard trigonometric identities to show that Ω(θ 1 )|w 1 − w 2 | 2 Ω(θ 2 ) = |r 1 −r 2 | 2 = 2(1 − cos θ 12 ). We have chosen the normalization of φ so that the numerator is unity. Just as Poincare invariance on the plane implies that correlators are a function of the length (or Euclidean distance on the plane, |w 1 − w 2 |), rotational invariance on the sphere fixes the correlator to be a function of the geodesic distance, θ 12 .
On R 2 the conformal n-point correlation functions of the 2-d Ising model can be constructed, in principle, to any order [37,38,39,40,41] and in practice have been computed up to sixth order [37,38,42]. To normalize the ratios Q 2n we use the analytic result, with ∆ = 1/8 for the Ising model. This allows us to find the relationship between the normalization of the continuum field φ(x) and the normalization of our discretized field φ x in a QFE calculation. Similarly, higher moments can be computed numerically as integrals over the n-point correlators on S 2 .
The integral of the four-point function on the sphere was performed by Deng and Blöte in Ref. [36]. The computation of m 4 is naïvely an eight-dimensional integral, but after utilizing rotational invariance the integration is reduced to a five-dimensional integral. The numerical four-point integral evaluation in Ref. [36] used 1000 independent Monte Carlo estimates, yielding an estimate of m 4 = 1.19878(2), leading to a prediction of the critical value for the fourth Binder cumulant U * 4 = 0.8510061(108) after error propagation.
We computed both the four-and six-point integrals numerically using the MonteCarlo method of Mathematica's Integrate[] function. For the four-point integral, we set AccuracyGoal→4, which yields a Monte Carlo distribution with standard deviation approximately 10 −4 . We repeated the Monte Carlo estimation of the integral 100 times. The quoted error is the standard error on the mean of this sample distribution. We find m 4 = 1.1987531(116) and the corresponding Binder cumulant is U * 4 = 0.8510207(63). We again used the Mathematica MonteCarlo integrator to compute m 6 . We set AccuracyGoal→3, which yields a Monte Carlo distribution with standard deviation approximately 10 −3 . We repeated the Monte Carlo estimation 50 times. We find m 6 = 1.632851(253), and a corresponding critical Binder cumulant value of U * 6 = 0.7731441(213).

Finite Scaling Fitting Methods
To compute these cumulants, we must obtain estimates for the even moments of the average magnetization. On the lattice we define the moments of the field on the simplicial complex by where . . . denotes the ensemble average and the spatial average is shown explicitly by a summation. Numerically, these moments are easy to compute by a weighted sum of the field over the complex on each field configuration, with the weights √ g x defined to be the areas of the Voronoï cells [43] as computed on the flat triangles and normalized to one per site on average, i.e. our convention N x=1 √ g x = N in Eq. (3.7).
In our Monte Carlo calculations, the finite size of our simplicial complex will break conformal invariance. If the bare lattice parameters are tuned sufficiently close to the critical point, finite-size scaling (FSS) relations can be used to extract CFT data by fitting the volume dependence of moments or cumulants of average magnetization [44]. We give the key details below, and in Appendix B we give further details of the FSS analysis used to parameterize our numerical data as we take the infinite volume limit.
It begins by expanding the free energy in a finite lattice volume N around the critical point (g * σ , g * , g * ω ). By Z 2 symmetry the critical point for the Z 2 odd operator is g * σ = 0. Of course, we cannot directly vary the renormalized couplings (g σ , g , g ω , · · · ) but instead can vary our bare couplings (h, µ 2 0 , λ 0 ) defined in our cut-off theory. The expansion then for the Z 2 odd operator takes the form 11) and without loss of generality, we can rescale h so that α 1 = 1. There is mixing between the two even operators so that (4.12) Since our current simulations were performed at a single fixed λ 0 , we cannot directly determine the coordinate transformation R. Instead, we define the derived quantities.
where a k1 , b k1 , and c k1 are defined in Eq. (B.11) and Eq. (B.12) in Appendix B. Future simulations varying the bare coupling near to the Wilson Fisher fixed points will improve the finite volume parameterization. The result of the finite volume and scaling expansion is to parameterize our data by where we introduce the length parameter L = √ N in accord with the FSS analysis summarized in Appendix B. The ellipses are a reminder that our expansion drops higher order terms in the fitting expressions, leading to a systematic error in the included parameters. For the d = 2 Ising model, ∆ σ = 1/8 and ∆ ω = 4, so the b 21 term scales like L −9/4 whereas the c 20 term scales like L −2 . Interestingly, the leading irrelevant correction to scaling is not due to a conformal quasi-primary operator ω, but rather a breaking of conformal symmetry due to a finite volume. We know of no other example of a CFT where this occurs.
Finally, to compare our simplicial calculations to the precise analytic results described in Sec. 4.1, we can fit our estimates of the magnetization moments and use the fitted parameters to estimate the critical Binder cumulants. Note that even though we chose to approach the critical surface along a line of constant λ 0 and were therefore unable to determine the linear coordinate transformation R of Eq. (4.12) leading to some of that dependence being absorbed into redefinitions of fit parameters, e.g. a 21 vs. a 21 , our estimates of the critical Binder cumulants U * 2n are free of this ambiguity.

Monte Carlo Results Near the Critical Surface
For the Monte Carlo simulation, we use the embedded dynamics algorithm of Brower and Tamayo [32], mixing a number of embedded Wolff cluster updates [45] Lattice size Wolff/local update ratio L < 15  with local updates consisting of one Rosenbluth-Teller sweep [46] and one over relaxation sweep [47,48,49,50]. Table 1 shows our empirically determined number of Wolff cluster updates per local update such that on average O(N ) spins are flipped between local updates. We can approximately locate the critical surface of the Wilson-Fischer fixed point by computing the Binder cumulants at fixed λ 0 , varying the volume N = L 2 = 2 + 10L 2 and the bare mass µ 2 0 , searching for the region of parameter space where U 4 ≈ 0.851. We can understand the approximate behavior of the Binder cumulants by substituting the FSS expressions in Sec. 4.2 and re-expanding [51], giving Since ∆ = 1, the A 2n term will cause the cumulant to diverge from its critical value if µ 2 0 is not tuned to µ 2 a as L → ∞. Importantly, the direction of the divergence will depend on the sign of (µ 2 0 − µ 2 a ). At smaller L, the C 2n term will tend to dominate since ∆ σ = 1/8 and ∆ ω = 4 as previously noted. It is our experience that Eq. (4.16) should not be used to actually fit the Binder cumulant data to accurately determine U * 2n since there tend to be delicate cancellations that occur in the expansion of the ratio. Instead, the moments should be analyzed separately using the parameterization in Eq. (4.15).
In the left panel of Fig. 4.1, we show the cumulants for a fixed size of the simplicial complex (L = 200) as we sweep through the critical region by varying µ 2 0 holding λ 0 = 1 fixed. Higher order cumulants are statistically noisier but since they pass through the critical region with a steeper slope [35] they have comparable statistical weight in determining the critical coupling µ 2 a ≈ µ 2 cr . On the right we show the cumulant U 4 vs. simplicial complex size L for fixed values of µ 2 0 in the critical region. As expected, the divergence at large L changes sign as µ 2 0 passes through the critical region.
We will now proceed to describe our data fitting and analysis procedure. To estimate quantities like U * 4 as accurately as possible, we use an iterative procedure to determine which ensembles, labeled by (µ 2 0 , λ 0 , L), to include in the fits to various moments of magnetization.   We start by getting a rough fit to the data using some initial data selection window µ 2 0 ∈ [µ 2 min , µ 2 max ] and L ∈ [L min , L max ] and then form the following quantities Subsequently we adjust the range of data to be fit for each moment of magnetization m n enforcing the data cuts: δa k1 , δb k1 , δc 20 ≤ w, where w is the width of the window, for all k ≤ n. It makes sense that lower moments of magnetization should have larger data selection regions since they vary more slowly in the critical region relative to the higher moments, as indicated in Fig. 4.1.
We then refit the new data selection and reselect the data for the same w using the new fit values and then iterate until the process converges to a stable data set for that w. We expect that fit parameters like a k0 will have systematic errors of O(w 2 ) due to higher order terms in the Taylor series expansion of the moments of the free energy which we did not include in our fit. So, we would like to take w to be as small as possible but also have enough data left to adequately constrain all of the important fit parameters, particularly the a k0 .
In Fig. 4.2 we see a subset of the data selected with w = 0.03 and curves with error bands computed from the best fit values of parameters shown in Table 2. Using the values shown in the table for w = 0.03, we get U 4,cr = 0.85020(58)(90) and U 6,cr = 0.77193 (37)(90) where the first error is statistical from the fit and the second is systematic. We see O(10 −3 ) agreement within the error estimates. Although even more stringent tests can be made, with the development of faster parallel code and a more extensive use of FSS analysis, we feel this is substantial support for the convergence of our QFE method to the exact c = 1/2 minimal CFT. In the next section we give further tests for the two-point and four-point correlation functions.

Conformal Correlator on S 2
With our results for the Binder cumulants, we are confident that the QFE lattice action on S 2 has the correct critical limit. Our analysis of magnetization moments and cumulants involves global averages over n-point correlators. Here we turn to correlators to examine more closely other consequences of conformal symmetry and to provide more stringent tests for our QFE lattice action. We begin with the two-point functions S 2 , and explain how, with rotational invariance, the exact analytic result can best be tested against our numerical result via a Legendre expansion. We next move to the exact four-point correlator, which takes on a relatively simple form as a sum of Virasoro blocks. Finally we use this as a toy example to learn how to compute CFT data (dimensions and couplings) from the conformal block operator product expansion (OPE). In particular, we extract the central charge c from the energy momentum tensor contribution to the scalar four-point function.
In Fig. 5.1, we compare the analytic value to the moments of our correlation data. The fit is excellent showing only slight deviations at high l due to cut-off effects. Our Legendre coefficients for various values of L, with the naïve dependence on l scaled out for clarity of presentation. On the right, we show the relative difference between the simplicial and continuum Legendre coefficients as we take the continuum limit for fixed l.
simplicial calculation of the two-point function is made at a fixed value of µ 2 0 closest to the pseudocritical values µ 2 a in Table 2. In Fig. 5.1, the normalization of the simplicial Legendre coefficients c l (L) are chosen so that c 0 (L) = c cont 0 = 8/7. In the left panel, we see that for any fixed value of L the difference between the simplicial and continuum values is small at small l but increases with increasing l. Yet as L increases, all simplicial coefficients converge towards their continuum values. The right panel shows how the simplicial coefficients behave at fixed l as we approach the continuum limit. The scaling curves shown assume a naïve 1/L scaling of simplicial artifacts, which turns out to be a good description of the data.
More generally in a CFT on the sphere, we will fit the two-point function to determine the operator dimension, ∆ σ . This can be done directly in co-ordinate space by fitting to Eq. (5.1) up to an overall normalization, or numerically to the Legendre coefficients, c cont l (∆ σ ) = 1 −1 dz (2/(1 − z)) ∆σ P l (z), which obey the closed form recursion relation, The overal normalization c cont 0 is of course proportional to m 2 defined in Eq. (4.2). With our current simulations of φ 4 theory on S 2 , both methods agree with the exact value ∆ σ = 0.125 at the percent level.
From our analysis of the two-point function, we have a direct demonstration of how rotational symmetry is restored in the continuum limit. Because the two-point function can be represented as an expansion in Legendre functions, it is a function of only the angle between any two points on the sphere and therefore rotationally invariant. For any fixed finite l, no matter how large, the simplicial coefficient c l (L) converges to the continuum one as L → ∞.

Four Point Functions
The c = 1/2 minimal CFT has only three Virasoro primaries 1, σ, , with an OPE expansion, In our earlier paper [19], we constructed a simplicial Wilson representation of Dirac Fermions S 2 . The map to the holomorphic, ψ(z), and anti-holomorphic,ψ(z), components of free Majorana fields on all 2-d Riemann surfaces [52] allowed us to define =ψψ and therefore compute the four-point functions: 1 2 3 4 and σ 1 1 σ 2 1 without Monte Carlo simulations. Here we test our QFE φ 4 construction on S 2 to compute the σ 1 σ 2 σ 3 σ 4 correlator. In the continuum limit the leading behavior of the lattice field is the primary operator: is a function invariant under any of the 24 permutations of the four positions. This permutation invariance is easily enforced even on a finite lattice. In the continuum, the eight real coordinates can be reduced to five real coordinates using rotational invariance alone.
Conformal symmetry allows for a further reduction to four real coordinates, the two angular variables θ 13 and θ 24 defined as in Eq. (5.1), plus two real conformal cross ratios u and v u ≡ |r 1 −r 2 | 2 |r 3 −r 4 | 2 5) or equivalently a single complex variable z.
where w is the complex variable for the R 2 prior to the stereographic project to S 2 in Eq. (4.6). In these variables, the four-point function can be written where the dependence on θ 13 and θ 24 cancels in the ratio. While G(z) has the exchange symmetries in the t-channel (1 ↔ 3, 2 ↔ 4), combining this properly with the Möbius map z → 1 − z, z → 1/z gives 24 copies, transforming the blue fundamental zone on the left-hand side of Fig. 5.2 to the entire complex z plane. In the specific case of the c = 1/2 minimal model, the explicit form of G(z) is known in polar coordinates z = re iθ as The first form is a sum of the two Virasoro blocks and the second form, given in Ref. [36], is the explicitly crossing symmetric form. For visualization, we plot G(z) in Fig. 5.2. Given that independent field configurations on our simplicial complexes can be generated in nearly O(L 2 ) time, the dominant cost of computing the four-point function would naïvely be computing the product of four fields which scales as O(L 8 ) and would be intractable for large L close to the continuum limit. We can reduce the cost by only sampling a subset of values of the four-point function on each independent configuration we generate. This is achieved by drawing O(L 2 ) quartets of random points on the simplicial complex. This balances the computational cost equally between generating independent field configurations and sampling the four-point function.
We then determine the product of the four simplicial fields, the unique coordinates (θ 13 , θ 24 , z) under the permutations of the positions, and finally an estimate of G(z) by dividing our product of simplicial fields by the continuum two-point function, Eq. (4.8), which is justified by our success of the previous section in showing that the simplicial two-point function converges correctly in the continuum limit, For convenience when studying the four-point function, we bin the data by coordinate on the complex plane. We partition the unit disk |z| ≤ 1 on the complex plane into a radial grid of N 2 0 bins of equal area. The angular width of each bin is ∆θ = 2π/N 0 and we choose N 0 mod 4 = 0 so that the bins are centered at θ n = 0, π/2, π, 3π/2. This arrangement is ideal for using standard discrete cosine transforms (DCT-I) for computing Fourier coefficients. For this study we chose N 0 = 64.
We present a subset of our results for the conformal portion of the four-point function, Eq.

Operator Product Expansion
In Sec. 5.2, we calculated the simplicial approximation to the scalar four-point function which, after taking a ratio with a product of two-point functions respecting t-channel symmetry, is described by a conformal function G(z) invariant under z → 1 − z and z →z. We computed the function everywhere in the unit circle, which can then be extended to the whole complex plane utilizing the permutation symmetry of the full fourpoint function. Our representation of the function G(z) is a tabulation of values on an equal area grid in polar coordinates (r, θ).
In the conformal bootstrap the expansion around z = 0 is used based on the operator product expansion (OPE). It is interesting to ask how well our Monte Carlo simulation on S 2 can determine the data in an OPE expansion, where G s (z) is s-channel symmetric, i.e. symmetric under interchange 1 ↔ 2 and 3 ↔ 4, ∆ O are the scaling dimensions of conformal primary operators O, and l labels the spin. The functions g ∆ O ,l (z) are called conformal blocks whose functional form is completely determined by conformal symmetry. In d = 2, there is an explicit representation of the conformal blocks in terms of hypergeometric functions, The expansion in Fourier modes is given by With a conserved energy momentum tensor T , there are special terms in the OPE expansion of the scalar four-point function which have integer powers: 1 + λ 2 T g T,2 (r, θ). The identity operator term is normalized to unity by convention. The conformal block for the energy momentum tensor has a closed form expression, (5.13) = 1 2 r 2 cos(2θ) + 1 2 r 3 cos(3θ) + 9 20 r 4 cos(4θ) + · · · , and the OPE coefficient λ 2 T is related to the central charge of the CFT, where C T = 2c = 2(1/2) = 1 for the d = 2, c = 1/2 minimal model. We can extract an estimate of the central charge c from our simplicial calculation of the scalar four-point function by modeling the result with the first few terms in the OPE expansion G s (r, θ) ∝ 1 + λ 2 g ,0 (r, θ) + λ 2 T g T,2 (r, θ) We used a DCT-I to Fourier transform the data, simultaneously fitting m = 0, 2 to four parameters: the overall normalization, λ 2 , λ 2 T and ∆ . We choose a single value of µ 2 = 1.822410 which is close to the pseudocritical value, and we restrict the fitting range in r around 0.5 where the data shows the least discretization error. A fit for L = 36 and 0.25 ≤ r ≤ 0.75 is shown in Fig. 5.4. A summary of several fits is shown in Table 3. As expected, the fit values converge towards the continuum CFT results (∆ = 1, λ 2 = 1/4, c = 1/2) as L increases towards the continuum limit and as the fit range in r is confined to a narrower region around r = 0.5 where the discretization artifacts appear smallest.

Conclusion
We have set up the basic formalism for defining a scalar field theory on non-trivial Riemann manifolds and tested it numerically in the limited context of the 2-d φ 4 theory by comparing it against the c = 1/2 Ising CFT at the Wilson-Fisher fixed point. This test is at a sufficient accuracy compared with the exact result to encourage us that we are able to define this strongly coupled quantum field theory on a curved manifold, in this case S 2 . Faster parallel code could push these tests to much higher precision. While this may be pursued in future works, our goal here was to sketch a framework for lattice quantum field theories on a nontrivial Riemann manifold as a quantum extension of FEM which we refer to as Quantum  Our construction first relies on theoretical issues in the application of FEM that, in spite of the vast literature, are not to our knowledge proven in enough generality. Nonetheless, the documented experience with FEM for PDEs does argue that this framework likely extends when properly formulated to convergence to the classical limit of any renormalizable quantum field theory on a smooth Euclidean Riemann manifold [19]. A rather non-trivial extension of the classical finite element/DEC method for the lattice Dirac Fermion and non-Abelian gauge fields on a simplicial Riemannian manifold is provided in Ref. [19] and Ref. [22] respectively. The local spin and gauge invariance on each site is achieved by compact spin connections and compact gauge links. This we believe provides the basic tools for classical simplicial constructions.
Second, and most importantly, we show that reaching the continuum for a quantum FEM Lagrangian requires a modification of the action by counter terms to accommodate UV effects. We conjecture that for a super-renormalizable theory with a finite number of divergent diagrams, a new QFE lattice action can be constructed by a natural generalization of φ 4 theory example presented here that removes the scheme dependence of the FEM simplicial lattice UV cut-off. The detailed generalization of this conjecture needs to be investigated with many more examples and most likely alternative approaches.
Third, once we have successfully constructed the QFE action with UV counter terms consistent with the continuum limit to all orders in perturbation theory, we conjecture that it represents a correct formulation of the non-perturbative quantum field theory on a curved Riemann manifold. This parallels the conventional wisdom for LFTs on hypercubic lattices that when lattice symmetries are protected (Lorentz, chiral, SUSY) only a finite set of fine tuning parameters are need to reach the continuum theory. Further research is need to establish the range of applicability of our QFE proposal. More examples are obviously need to support this proposal.
There are many more direct extensions of this formulation that we are currently entertaining. The first is to apply QFE to the 3-d φ 4 theory in both radial quantization and the 3-d projective sphere, S 3 . Our current software relies on a serial code and the very efficient Brower-Tamayo [32] modification of the cluster algorithm. We are in the process of developing parallel code that will allow for large lattices and higher precision, which will also be a necessity in going beyond 2d and scalar fields. Although new software requires substantial effort, we believe that all of the fundamental data parallel concepts and advanced algorithms utilized in lattice QCD will still be applicable here. We note there are advantages to studying CFTs on spherical S d simplicial lattices. There are no finite volume approximations. The entire R d is mapped to the sphere. For radial quantization the only finite volume effect is the one radial dimension along the R × S d−1 cylinder. With periodic boundary conditions, the finite extent of the cylinder is also an interesting parameter for the study of CFTs at finite temperature. Radial quantization on R × S d−1 allows direct access to the Dilatation operator to compute operator dimensions and the OPE expansion. Even our 2-d example is worth further investigation at higher precision, and with an emphasis on new studies enabled by a fully non-perturbative formulation of the φ 4 theory on S 2 . For example, we can consider new ways to compute the renormalization flow of the central charge from the UV to the IR. Further research is on going to hopefully justify this optimism.
during the completion of this work. R.C.B. and E.S.W. acknowledge support by DOE grant DE-SC0015845 and T.R. and C.-I T. in part by the Department of Energy under contact DE-SC0010010-Task-A. T.R. is also supported by the University of Kansas Foundation Distinguished Professor Grant of Christophe Royon.

B Finite Size Scaling for Cumulants and Moments
This follows the renormalization group arguments first presented in [44] and reviewed in [51]. We begin with the free energy Eq. expressed not as a function of the bare couplings h, µ 2 0 , λ 0 but rather the relevant couplings g σ , g λ and the irrelevant couplings {g ω , . . . } of the Wilson-Fisher fixed point. For the lattice spacing, we assume a ∼ 1/L. Under a change of scale by a factor l, the free energy renormalizes F (g σ , g , {g ω , . . . }, a) = l −d F (l yσ g σ , l y g , {l yω g ω , . . . }, la) + G (g σ , g ) (B.2) where y O ≡ d − ∆ O , F is the singular, or scaling, part of the free energy and G is the regular part. In what follows, we will continue to use y O instead of ∆ O for compactness of notation.
Before we continue, we point out an important physical concept. At finite L, the difference between the lattice scalar operator φ and the conformal primary operator σ is indicated by the presence of higher derivative terms α 2n+1 in the cumulant expressions. But, in the scaling limit L → ∞ all those terms vanish and the cumulants are dominated by the moments of the free energy with respect to the coupling g σ of the conformal primary operator σ. In this sense, the simplicial operator φ(x) becomes the primary operator σ(x) in the scaling limit.
We can now apply the important physical principle that close to the critical point we can expand the derivatives of the free energy in a Taylor series F (k) = a k0 + a k1 (g − g * ) L y + · · · + b k1 (g ω − g * ω ) L yω + · · · . (B.11) G (k) = c k0 + c k1 (g − g * ) + · · · (B.12) Substituting these expressions into those for the cumulants defined in Eqns. (B.7) to (B.10) nearly gives us an expression we can fit to computed data. Of course, we cannot directly vary the renormalized couplings (g σ , g , g ω , · · · ) but instead can vary our bare couplings (h, µ 2 0 , λ 0 ) defined in our cutoff theory. Near the Wilson-Fisher fixed point we can relate the two sets of couplings by expanding So, substituting Eq. (4.12) first into Eqns. (B.11) to (B.12) leads to expressions that can be used to model cumulant data computed in a cutoff theory close to the Wilson-Fisher fixed point.