# The Operator Product Expansion for Radial Lattice Quantization of 3D
$\phi^{4}$ Theory

Venkitesh Ayyar  Department of Physics and Center for Computational Science,
Boston University, Boston, MA 02215, USA  Richard C. Brower  Department of
Physics and Center for Computational Science, Boston University, Boston, MA
02215, USA  George T. Fleming<EMAIL_ADDRESS>Fermi National Accelerator
Laboratory, Batavia, IL 60510, USA  Department of Physics, Sloane Laboratory,
Yale University, New Haven, CT 06511, USA  Anna-Maria E. Glück  anna-
<EMAIL_ADDRESS>Kirchhoff-Institut für Physik, Ruprecht-
Karls-Universität Heidelberg, 69120 Heidelberg, Germany  Department of
Physics, Sloane Laboratory, Yale University, New Haven, CT 06511, USA  Evan K.
Owen  Department of Physics and Center for Computational Science, Boston
University, Boston, MA 02215, USA  Timothy G. Raben  Department of Physics and
Astronomy, Michigan State University, East Lansing, MI 48824, USA  Chung-I Tan
Brown Theoretical Physics Center and Department of Physics, Brown University,
Providence, RI 02912, USA

###### Abstract

At its critical point, the three-dimensional lattice Ising model is described
by a conformal field theory (CFT), the 3d Ising CFT. Instead of carrying out
simulations on Euclidean lattices, we use the Quantum Finite Elements method
to implement radially quantized critical $\phi^{4}$ theory on simplicial
lattices approaching $\mathbb{R}\times S^{2}$. Computing the four-point
function of identical scalars, we demonstrate the power of radial quantization
by the accurate determination of the scaling dimensions $\Delta_{\epsilon}$
and $\Delta_{T}$ as well as ratios of the operator product expansion (OPE)
coefficients $f_{\sigma\sigma\epsilon}$ and $f_{\sigma\sigma T}$ of the first
spin-0 and spin-2 primary operators $\epsilon$ and $T$ of the 3d Ising CFT.

Conformal Field Theory, Lattice Field Theory, Critical Phenomena,
Renormalization Group

††preprint: FERMILAB-PUB-23-631-T

##  I  Introduction

Since the first discovery of critical opaqueness made by Andrews in the 19th
century [1] and especially since Wilson’s breakthrough formulation of the
Renormalization Group framework for critical phenomena [2], the study of
systems at criticality has fascinated physicists throughout the decades. The
prototype of such a system with a continuous phase transition and the
workhorse in studying critical phenomena is the Ising model. Despite its
simplicity, it is of great physical relevance because universality guarantees
that the critical Ising model shows the same behavior as encountered for
critical transitions in e.g. uniaxial magnetic systems, fluids [3], and
micellar systems [4].

Motivated by experimental measurements of the scaling behaviour of physical
quantities near criticality, past theoretical efforts to study the critical
point of the Ising model had the main goal of calculating its critical
exponents. As there is no analytic solution for the Ising model in $d\geq 3$,
such calculations have to be performed using approximate and computational
methods, the most prominent of those including high-temperature expansions,
field-theoretic methods such as the $\varepsilon$-expansion or non-
perturbative methods based on approximate solutions of the RG-equations, as
well as Monte Carlo simulations. For a review, see [5]. After the development
of spin cluster algorithms [6, 7], especially the latter method combined with
finite-size scaling arguments yielded precise results for the critical
exponents [5, 8].

For comparison with experiment, it is sufficient to determine the critical
exponents of the 3d Ising model. But from a theoretical point of view, the
critical Ising model contains much more information. Following Polyakov’s
hypothesis [9], we expect the scale invariance of a system at criticality
combined with Poincaré invariance to lead this critical system to be invariant
under the conformal group $O(d+1,1)$. Consequently, the physics of the
critical 3d Ising model can be described by a conformal field theory (CFT),
the 3d Ising CFT. The local operator content of such a CFT is spanned by the
so-called primary operators and their descendants, which are created by acting
on the primaries with translation generators. The primary operators can be
assigned spins and parity and, in the case of the 3d Ising CFT, include the
leading $0^{-}$-operator $\sigma$, the leading $0^{+}$-operator $\epsilon$ and
the energy-momentum tensor $T$ ($2^{+}$), as well as infinitely many
subleading and higher-spin operators. As will be further explained in Sec. II,
a CFT is characterized by two sets of quantities - the scaling dimensions
$\Delta_{\mathcal{O}}$ and the operator product expansion (OPE) coefficients
$f_{\mathcal{O}_{1}\mathcal{O}_{2}\mathcal{O}_{3}}$ for different primary
operators $\mathcal{O}$ of the CFT.

For (at least some of) the scaling dimensions, there is a mapping to the
critical exponents, which have been thoroughly investigated with the
traditional methods reviewed in [5]

$\displaystyle\eta$ $\displaystyle=2\Delta_{\sigma}-d+2,$ $\displaystyle\nu$
$\displaystyle=1/(d-\Delta_{\epsilon}),$ $\displaystyle\omega$
$\displaystyle=\Delta_{\epsilon^{\prime}}-d.$ (1)

The OPE coefficients, however, are largely uninvestigated with these methods.
With Functional Renormalization Group methods, it was only recently possible
to determine $f_{\sigma\sigma\epsilon}$ [10]. Moreover,
$f_{\sigma\sigma\epsilon}$ and $f_{\epsilon\epsilon\epsilon}$ have been
measured with traditional Monte Carlo methods on Euclidean lattices from
three-point correlators [11, 12] and from the scaling of 2-point correlators
[13, 14] of the operators $\sigma$ and $\epsilon$.

However, as suggested by Cardy [15] in 1985, there is, in principle, a huge
advantage by replacing the traditional toroidal Euclidean lattice
approximation to ${\mathbb{R}}^{d}$ with simulations on simplicial lattices
convergent to the cylindrical ${\mathbb{R}}\times{\mathbb{S}}^{d-1}$. This is
an exact Weyl map of all CFT data to what is called radial quantization [16]
with translations generated by the dilatation operator. The challenge [15] for
$d>2$ is to define the correct lattice action on the curved manifold.
Recently, the Quantum Finite Elements (QFE) [17] project has begun to address
this by including lattice counter terms to improve the UV cut-off.

QFE tests began with $\phi^{4}$ theory at the critical point in comparison
with the minimal 2d Ising CFT on ${\mathbb{S}^{2}}$, accurately reproducing
the exact results for the $\mathbb{Z}_{2}$-odd scalar propagator and the
scalar four-point amplitude [18] as well as examples of the
$\mathbb{Z}_{2}$-even sector in the fermionic representation [19] on
${\mathbb{S}^{2}}$. The next step was to apply the QFE to critical $\phi^{4}$
theory on $\mathbb{R}\times{\mathbb{S}^{2}}$ in order to compare to the 3d
Ising CFT. The scalar propagator [20] was well represented in the continuum
limit with the scaling dimension $\Delta_{\sigma}$ determined to better than
$10^{-2}$ agreement with the conformal bootstrap [21, 22, 23, 24]. Here, we
extend the QFE method to the scalar four-point amplitude of the 3d Ising CFT.
The present work is intended as a proof of concept to show that by fitting our
lattice data to the expected form from the continuum operator product
expansion (OPE) for CFTs in radial quantization, we can extract scaling
dimensions $\Delta_{\mathcal{O}}$ and OPE coefficients
$f_{\sigma\sigma\mathcal{O}}$ of $\mathbb{Z}_{2}$-even primary operators in
the continuum limit. Because an exact OPE is an intrinsic property of CFTs,
this provides a stringent test of whether we actually reach the 3d Ising CFT
in the continuum with our QFE radial quantization. Preliminary results
presented at LATTICE 2022 [25] are extended and finalized here to higher
fidelity to assess the capabilities of this method.

In the conclusion, we discuss further improvements of the QFE method and make
comparisons with other approaches. In particular, the conformal bootstrap
spectacularly constrains the four-point functions, giving rigorous bounds not
only on the scaling dimensions of whole sets of conformal primary operators
but also on their OPE coefficients [23, 24]. Even for non-integral theories,
these offer definitive tests for Monte Carlo simulations. In addition, there
are recent successes in determining OPE coefficients using Hamiltonian methods
on the fuzzy sphere [26, 27]. While each method is complementary, comparisons
may suggest improvement for them all.

##  II  $n$-point functions in radially quantized CFTs

###  II.1 CFTs and Radial Quantization

At the Wilson-Fisher fixed point, 3d $\phi^{4}$ theory is described by the 3d
Ising conformal field theory (CFT). $d$-dimensional CFTs are quantum field
theories that are invariant under the conformal group $O(d+1,1)$ of
transformations that leave the metric invariant up to a Weyl rescaling [28]

$g_{\mu\nu}(\bm{x})\rightarrow
g_{\mu\nu}^{\prime}(\bm{x^{\prime}})=\Lambda^{2}(\bm{x})g_{\mu\nu}(\bm{x}).$
(2)

Thus, instead of considering the 3d Ising CFT on three-dimensional Euclidean
space with metric

$ds_{flat}^{2}=dr^{2}+r^{2}d\Omega^{2}=\frac{r^{2}}{R^{2}}\left[dt^{2}+R^{2}d\Omega^{2}\right]$
(3)

where $t=R\log(r/R)$, one can quantize the theory on $\mathbb{R}\times S^{2}$
with Weyl-rescaled metric

$ds_{cyl}^{2}=dt^{2}+R^{2}d\Omega^{2}=\frac{R^{2}}{r^{2}}ds_{flat}^{2}.$ (4)

$n$-point functions of scalar fields in flat space are then related to those
on the cylinder via

$\displaystyle\langle\sigma(r_{1},\Omega_{1})\sigma(r_{2},\Omega_{2})\cdots\rangle_{flat}$
$\displaystyle=\left(\frac{R}{r_{1}}\right)^{\Delta_{\sigma}}\left(\frac{R}{r_{2}}\right)^{\Delta_{\sigma}}\cdots\langle\sigma(t_{1},\Omega_{1})\sigma(t_{2},\Omega_{2})\cdots\rangle_{cyl}$
(5)

with $\Delta_{\sigma}$ the scaling dimension of scalar primary $\sigma$ [29].
This choice of coordinates is called radial quantization and is often useful
for scale-invariant theories as dilatations of $r$ on $\mathbb{R}^{3}$ are
mapped to shifts in cylinder time $t$. As will be explained in subsection
II.3, the determination of CFT data is also extremely facilitated by radial
quantization.

###  II.2 Conformal $n$-point functions

Due to the symmetries of CFTs, two-point functions of primary operators follow
power laws [28]

$\displaystyle\langle\mathcal{O}(\bm{x_{1}})\mathcal{O}(\bm{x_{2}})\rangle=\frac{C_{\mathcal{O}}}{|\bm{x_{1}}-\bm{x_{2}}|^{2\Delta_{\mathcal{O}}}}.$
(6)

Moreover, the operator product expansion (OPE) is exact, such that under the
path integral, a product of two primaries can be replaced by a sum over
“intermediate states” (schematically)

$\displaystyle\mathcal{O}_{i}(x_{1})\mathcal{O}_{j}(x_{2})=\sum_{k}f_{\mathcal{O}_{i}\mathcal{O}_{j}\mathcal{O}_{k}}C(x_{1}-x_{2})\mathcal{O}_{k}(x_{2}),$
(7)

with the OPE coefficients $f_{\mathcal{O}_{i}\mathcal{O}_{j}\mathcal{O}_{k}}$
as expansion coefficients. Thus, all $n$-point functions of a CFT can be
expressed as functions of the scaling dimensions $\Delta_{\mathcal{O}}$ and
the OPE coefficients $f_{\mathcal{O}_{i}\mathcal{O}_{j}\mathcal{O}_{k}}$ of
conformal primaries $\mathcal{O}$. Moreover, all those quantities can be
extracted by measuring different four-point functions.

In this work, we focus on the four-point function of four identical scalar
primaries $\sigma$, which - when divided by the corresponding two-point
functions - forms the amplitude $g(u,v)$. It is called the conformally
invariant amplitude as any factors Eq.(5) from Weyl rescaling the metric
Eq.(2) drop out. Thanks to the OPE, this amplitude can be expanded into a sum
over the primaries $\mathcal{O}$ of the CFT [30, 31]:

$\displaystyle g(u,v)$
$\displaystyle\equiv\frac{\langle\sigma(x_{1})\sigma(x_{2})\sigma(x_{3})\sigma(x_{4})\rangle}{\langle\sigma(x_{1})\sigma(x_{2})\rangle\langle\sigma(x_{3})\sigma(x_{4})\rangle}$
(8)
$\displaystyle=1+\sum_{\mathcal{O}}f_{\sigma\sigma\mathcal{O}}^{2}G_{\mathcal{O}}(\Delta_{\mathcal{O}};u,v).$
(9)

In this case, actually only the primaries $\mathcal{O}$ with even spin and
parity contribute because of the symmetries of the amplitude. The
$G_{\mathcal{O}}$ are the so-called conformal blocks, functions of the scaling
dimensions $\Delta_{\mathcal{O}}$ of primaries $\mathcal{O}$ as well as the
kinematic information encoded in the conformally invariant cross ratios $u$
and $v$

$\displaystyle\sqrt{u}=\frac{|\bm{x_{1}}-\bm{x_{2}}||\bm{x_{3}}-\bm{x_{4}}|}{|\bm{x_{1}}-\bm{x_{3}}||\bm{x_{2}}-\bm{x_{4}}|},\
\sqrt{v}=\frac{|\bm{x_{1}}-\bm{x_{4}}||\bm{x_{2}}-\bm{x_{3}}|}{|\bm{x_{1}}-\bm{x_{3}}||\bm{x_{2}}-\bm{x_{4}}|}.$
(10)

In odd dimensions, there is no known analytic form of the conformal blocks.
However, if we express the $G_{\mathcal{O}}$ as a function of the (also
conformally invariant) quantities

$\cosh(\tau)=\frac{1+\sqrt{v}}{\sqrt{u}},\ \ \ \
\cos(\alpha)=\frac{1-\sqrt{v}}{\sqrt{u}},$ (11)

they can be expanded in an absolutely convergent series in the Gegenbauer
polynomials [32, 33]:

$G_{\mathcal{O}}(\tau,\alpha)=\sum_{n\in
2\mathbb{N}_{0}}e^{-(\Delta_{\mathcal{O}}+n)\tau}\sum_{j}B_{n,j}(\Delta_{\mathcal{O}})C_{j}^{d/2-1}(\cos\alpha),$
(12)

where in $d=3$ dimensions the Gegenbauer polynomials correspond to the
Legendre polynomials $C_{j}^{d/2-1}(\cos\alpha)=P_{j}(\cos\alpha)$. The $j$ in
the second sum is determined by the usual rules for the addition of angular
momenta between the spin $l$ of the primary $\mathcal{O}$ and $n$,

$j\in\\{\max(0,l-n),\cdots,l+n\\},$ (13)

and the $B_{n,j}$ are normalization coefficients. In three dimensions, they
are normalized by

$B_{0,j}(\Delta_{\mathcal{O}})=4^{\Delta_{\mathcal{O}}}\frac{l!}{(1/2)_{l}}\delta_{j,l},$
(14)

where $()_{l}$ is the Pochhammer symbol. For their recursive calculation, we
used the Mathematica Notebook provided with Ref. [33].

###  II.3 The antipodal scalar four-point function on $\mathbb{R}\times S^{2}$

Figure 1: Schematic illustration of the antipodal frame on $\mathbb{R}\times
S^{2}$ in which we calculate the 4-point function

While in Euclidean space $\alpha$ and $\tau$ are not accessible, that changes
if we radially quantize the theory. If we measure the four-point function on
$\mathbb{R}\times S^{2}$ with $x_{i}=(t_{i},\bm{n_{i}}),\ i\in\\{1,2,3,4\\}$
placed such that they lie pairwise at identical times, $t_{1}=t_{2}$ and
$t_{3}=t_{4}$, and on antipodal points of $S^{2}$, $\bm{n}_{1}=-\bm{n}_{2}$
and $\bm{n}_{3}=-\bm{n}_{4}$ (see Fig. 1), we have [20]

$\bm{n_{1}}\cdot\bm{n_{4}}\equiv\cos(\theta)=\cos(\alpha),\ \ \
\cosh(t)=\cosh(\tau)$ (15)

such that we can directly measure

$\displaystyle
g(t,\cos\theta)=1+\sum_{\mathcal{O}}f_{\sigma\sigma\mathcal{O}}^{2}\sum_{n\in
2\mathbb{N}_{0}}e^{-(\Delta_{\mathcal{O}}+n)t}$ (16)
$\displaystyle\times\sum_{j}B_{n,j}(\Delta_{\mathcal{O}})P_{j}(\cos\theta).$

This is the big advantage of doing lattice calculations in radial quantization
compared to usual Euclidean lattices. Now $g(t,\cos\theta)$ has the form of a
partial wave expansion in the chosen antipodal frame

$g(t,\cos\theta)=1+\sum_{j}c_{j}(t)P_{j}(\cos\theta)$ (17)

with expansion coefficients

$\displaystyle
c_{j}(t)=\sum_{\mathcal{O}}f_{\sigma\sigma\mathcal{O}}^{2}\sum_{\begin{subarray}{c}n\in
2\mathbb{N}_{0}\\\
n\geq|j-l|\end{subarray}}^{\infty}B_{n,j}(\Delta_{\mathcal{O}})e^{-(\Delta_{\mathcal{O}}+n)t}$
(18)

so that each expansion coefficient consists of infinite towers of exponentials
for all even parity and even spin primaries. By fitting these expansion
coefficients, we can extract the quantities $\Delta_{\mathcal{O}}$ and
$f_{\sigma\sigma\mathcal{O}}$.

##  III Radially quantizing the 3d Ising CFT with Quantum Finite Elements

In Refs. [18, 20], Brower et al. succeeded at radial lattice quantization of
$\phi^{4}$ theory. For that, they approximated $S^{2}$ via a series of
simplicial lattices as shown in Fig. 2.

Figure 2: Construction of the simplicial lattice for refinement $L=3$

Starting from an icosahedron, its edges are first subdivided into $L$ equal
pieces, thereby introducing a finer triangulation. Subsequently, the
$N=10L^{2}+2$ vertices are projected onto the sphere. Using the finite element
method, the classical action can then be discretized on this simplicial
complex [20]:

$\displaystyle S=$ $\displaystyle\frac{1}{2}\sum_{t,x}\sum_{y\in\langle
x,y\rangle}\frac{l^{*}_{xy}}{l_{xy}}\left(\tilde{\phi}_{t,x}-\tilde{\phi}_{t,y}\right)^{2}+\frac{a^{2}}{4R^{2}}\sqrt{\tilde{g}_{x}}\tilde{\phi}_{t,x}^{2}$
(19)
$\displaystyle+\sqrt{\tilde{g}_{x}}\left[\frac{a^{2}}{a_{t}^{2}}\left(\tilde{\phi}_{t,x}-\tilde{\phi}_{t+1,x}\right)^{2}+m_{0}^{2}\tilde{\phi}_{t,x}^{2}+\lambda_{0}\tilde{\phi}_{t,x}^{4}\right],$

where we sum over $x\in\\{1,2,\cdots,N\\}$ on each sphere and
$t\in\\{1,2,\cdots,L_{t}\\}$ along the cylinder with periodic boundary
conditions. Here, $\tilde{\phi}_{t,x}$ is the dimensionless field, $m_{0}^{2}$
and $\lambda_{0}$ the dimensionless mass and coupling. $a$ and $a_{t}$ are the
average lattice spacings on the sphere and along the cylinder, respectively,
and $a$ is related to the radius $R$ of the sphere via

$a^{2}/R^{2}=\frac{8\pi}{\sqrt{3}N}=\frac{8\pi}{\sqrt{3}(10L^{2}+2)}.$ (20)

The bare speed of light $a/a_{t}$ is set to one, though it is renormalized by
interactions. For further details on the construction of this action as well
as on how the dimensions of the quantities are restored, see Refs. [18, 20].

While the position-dependent finite element weights $l^{*}_{xy}/l_{xy}$ and
$\sqrt{\tilde{g}_{x}}$ ensure the convergence of the action in Eq. (19) to the
spherically symmetric classical continuum theory [18], Brower et al. found
that for the quantum theory, they introduce fluctuations in the effective
cutoff that are amplified by UV divergences, preventing convergence to the
spherically symmetric quantum continuum theory [20]. To compensate for these
quantum UV defects, perturbative counterterms, the so-called Quantum Finite
Elements, were introduced to the lattice action:

$S_{QFE}=S-\sum_{t,x}\sqrt{\tilde{g}_{x}}\left(6\lambda_{0}\delta
G_{t,x}-24\lambda_{0}^{2}\delta G_{t,x}^{(3)}\right)\tilde{\phi}^{2}_{t,x},$
(21)

where $\delta G_{t,x}$ and $\delta G_{t,x}^{(3)}$ are calculated numerically
from the free lattice propagator $G_{t,x;t^{\prime},y}$ via

$\displaystyle\delta G_{t,x}$ $\displaystyle\equiv
G_{t,x;t,x}-\frac{1}{N}\sum_{x^{\prime}=1}^{N}\sqrt{\tilde{g}_{x^{\prime}}}G_{t,x^{\prime};t,x^{\prime}},$
(22) $\displaystyle\delta G_{t,x}^{(3)}$
$\displaystyle\equiv\sum_{t^{\prime},y}\sqrt{\tilde{g}_{y}}\left(G_{t,x;t^{\prime},y}^{3}-\frac{1}{N}\sum_{x^{\prime}=1}^{N}\sqrt{\tilde{g}_{x^{\prime}}}G_{t,x^{\prime};t^{\prime},y}^{3}\right).$
(23)

In order to tune this QFE action to the Ising critical surface, the Binder
cumulant was studied for fixed $\lambda_{0}=0.2$, finding
$m_{0}^{2}=-\mu_{0}^{2}=-0.27018(4)$ as the critical coupling. With the
parameters tuned to these values, the study of the scalar two-point function
carried out in [20] suggests that we reach the 3d Ising CFT in the continuum
limit $a/R\rightarrow 0$. In the following, we will fix the bare parameters to
these values for our analysis of the interacting theory.

##  IV  Calculating and fitting the antipodal lattice four-point function

To compute the antipodal scalar four-point function on the lattice, we used
our QFE action (21) and carried out Monte Carlo simulations with the Brower-
Tamayo cluster algorithm [34] combined with Metropolis [35] and overrelaxation
[36, 37, 38, 39]. After some initial equilibration sweeps,
$\mathcal{O}$(volume) four-point function configurations of our scalar lattice
fields $\tilde{\phi}$ were sampled in each sweep along with the corresponding
equal-time 2-point functions to obtain $g(u,v)$ as in Eq. (8). Subsequently,
we projected onto Legendre polynomials to obtain the expansion coefficients
$c_{j}(t)$ in Eq. (12). Note that for finite lattice spacing, the lattice
operator $\tilde{\phi}$ does not exactly correspond to the primary $\sigma$ as
in Eq. (12) but it is a mixture of all $\mathbb{Z}$-even scalar primaries of
the CFT. However, as $a/R\rightarrow 0$, the contribution of operators other
than $\sigma$ go to zero.

Such simulations were carried out for different lattice refinements
$L\in\\{24,28,32,36,40,44,48,56,64\\}$ in order to take the continuum limit
later by sending $L\rightarrow\infty$, corresponding to $a/R\rightarrow 0$. We
performed $N=8000$ independent simulations for $L=24$, $N=1600$ for $32\leq
L\leq 56$, and $N=800$ for $L=64$. The resulting data is freely available in
the Zenodo repository associated with this project [40]. For $L=64$, the
measured expansion coefficients up to $j=20$ are shown in Fig. 3 with the
error bars determined by the statistical error from averaging over the $N$
independent runs. This figure illustrates that even for the $L$-value with the
least amount of statistics, we have several timeslices of good data up until
high $j$, such that it should in principle be possible to extract information
about high-spin operators.

Figure 3: Monte Carlo results for the four-point function expansion
coefficients $c_{j}(t)$ from $N=800$ independent runs with statistical errors.
The coefficients are plotted for a $t$-range such that the relative error of
the effective mass does not exceed $\Delta m_{eff}/m_{eff}\leq 0.125$.

And in fact, by simply fitting unconstrained sets of exponentials to the
different $c_{j}$, we could extract five exponentials for $j=0$, and even for
$j=10$ we could fit up to 3 exponentials. However, the statistics did not
allow, especially for the higher $j$, to always resolve the contribution from
different primaries. For example, in the case of $j=2$, we expect from Eq.
(18) and the bootstrap results in Tab. 1 that the leading exponential $\propto
e^{-\Delta_{T}t}$ involves the scaling dimension $\Delta_{T}$ of the energy-
momentum tensor $T$ with $\ell=2$, while the next-lowest exponential $\propto
e^{-(\Delta_{\epsilon}+2)t}$ should stem from the $\ell=0$ operator
$\epsilon$. To be able to resolve these two contributions at $L=64$ without
having to introduce additional constraints in the fitting parameters such as
Eq. (18), we would likely have to increase the statistics by a factor 9-16,
which is computationally not feasible as already the $N=800$ runs took
$\mathcal{O}(6\times 10^{5})$ hours of computing time on single cores on the
BU shared computing cluster.

Operator $\mathcal{O}$ ($l^{P}$) | $\Delta_{\mathcal{O}}$ | $f_{\sigma\sigma\mathcal{O}}$
---|---|---
$\epsilon$ ($0^{+}$) | 1.41265(36) | 1.05185(12)
$\epsilon^{\prime}$ ($0^{+}$) | 3.82951(61) | 0.05304(16)
$T$ ($2^{+}$) | 3 | 0.32613776(45)
$T^{\prime}$ ($2^{+}$) | 5.499(17) | 0.01054(10)
Table 1: Conformal bootstrap values for scaling dimensions and OPE
coefficients of conformal primaries $\epsilon,T,\epsilon^{\prime}$, and $T$
[24, 23]

Thus, instead of fitting each $c_{j}$ independently with unconstrained
exponentials, we analyzed the Monte Carlo data by fitting the $c_{j}$
simultaneously for different $j$ with much more constrained functions based on
the expansion Eq. (18). This has the advantage that adding a new primary
operator $\mathcal{O}$ to the fit only introduces two new fit parameters in
total, namely $f_{\sigma\sigma\mathcal{O}}$ and $\Delta_{\mathcal{O}}$, while
describing (in principle infinite) towers of exponentials for each $c_{j}$.
The scaling dimension and OPE coefficients of a spin-$l$ operator are then
mainly determined by the operator’s contribution to $c_{\ell}$ while
appropriately contributing to higher-order terms in $c_{j}$ with $l\neq j$.
Overall, this procedure improves convergence and helps resolve the different
contributions.

Of course, strictly speaking, Eq. (18) is only true in the continuum limit as
our finite lattice spacing breaks conformal invariance. However, we still
performed fits based on the form of Eq. (18) because we expect that if $a/R$
is small enough, the corrections due to the breaking of lattice spacing can be
absorbed into corrections to our fit parameters that vanish as $a/R\rightarrow
0$.

In the present proof of principle, we focused on simultaneous fits to $c_{0}$
and $c_{2}$, for which the leading contributions come from the primaries
$\epsilon$ $(0^{+})$ and the energy-momentum tensor $T$ $(2^{+})$,
respectively. Furthermore, we included the first subleading operators for each
$\ell$, namely $\epsilon^{\prime}$ $(0^{+})$ and $T^{\prime}$ $(2^{+})$ in our
fits. We truncated the sums of exponentials at $n_{max}=20$ for $\epsilon$ and
correspondingly lower values for the other operators such that the largest
exponents are of the same order of magnitude. With these truncations, Eq. (18)
yields the following fit functions for $c_{0}(t)$ and $c_{2}(t)$ (in terms of
the lattice “time” $t$)

$\displaystyle c_{0}^{fit}(t)=$
$\displaystyle\sum_{n=0}^{n_{max}}f_{\sigma\sigma\epsilon}^{2}B_{n,0}(\Delta_{\epsilon})e^{-(\Delta_{\epsilon}+n)ta_{t}/R}$
(24)
$\displaystyle+\sum_{n=0}^{n_{max}-2}f_{\sigma\sigma\epsilon^{\prime}}^{2}B_{n,0}(\Delta_{\epsilon^{\prime}})e^{-(\Delta_{\epsilon^{\prime}}+n)ta_{t}/R}$
$\displaystyle+\sum_{n=2}^{n_{max}-4}f_{\sigma\sigma
T^{\prime}}^{2}B_{n,0}(\Delta_{T^{\prime}})e^{-(\Delta_{T^{\prime}}+n)ta_{t}/R}$
$\displaystyle+(t\rightarrow L_{t}-t),$ $\displaystyle c_{2}^{fit}(t)=$
$\displaystyle\sum_{n=0}^{n_{max}-2}f_{\sigma\sigma
T}^{2}B_{n,2}(\Delta_{T})e^{-(\Delta_{T}+n)ta_{t}/R}$ (25) $\displaystyle+$
$\displaystyle\sum_{n=0}^{n_{max}-4}f_{\sigma\sigma
T^{\prime}}^{2}B_{n,2}(\Delta_{T^{\prime}})e^{-(\Delta_{T^{\prime}}+n)ta_{t}/R}$
$\displaystyle+$
$\displaystyle\sum_{n=2}^{n_{max}}f_{\sigma\sigma\epsilon}^{2}B_{n,2}(\Delta_{\epsilon})e^{-(\Delta_{\epsilon}+n)ta_{t}/R}$
$\displaystyle+$
$\displaystyle\sum_{n=2}^{n_{max}-2}f_{\sigma\sigma\epsilon^{\prime}}^{2}B_{n,2}(\Delta_{\epsilon^{\prime}})e^{-(\Delta_{\epsilon^{\prime}}+n)ta_{t}/R}$
$\displaystyle+(t\rightarrow L_{t}-t),$

where the quantities $\Delta_{\mathcal{O}}$ and
$f^{2}_{\sigma\sigma\mathcal{O}}$ are taken as the fit parameters. Note that
we did not include contributions from the energy-momentum tensor $T$ to
$c_{0}(t)$ because for the theoretical continuum value $\Delta_{T}=3$, all
$B_{n,0}(3)$ vanish. In order to convert the lattice time of our data, which
is given in units of $a_{t}$, to the time on $\mathbb{R}\times S^{2}$ in Eq.
(18), which after the Weyl transform has units of $R$, we included factors
$a_{t}/R$ in the exponentials. This can be written as $a_{t}/R=a/(Rc_{R})$
where $a/R$ is known from Eq. (20) and $c_{R}$ is the renormalized speed of
light $c_{R}=a/a_{t}$. $c_{R}$ was determined non-perturbatively in Ref. [20],
yielding $c_{R}=0.996$. Furthermore, to take into account wraparound the
effects due to the periodicity along the cylinder to leading order, we added
for each exponential in our fit function the same term with $t$ substituted by
$L_{t}-t$.

Figure 4: Lattice data for $c_{0}$ and $c_{2}$ at a refinement of $L=64$
compared to theoretical curves based on truncations of Eq. (18). For each
curve, we include another primary operator contribution, corresponding to a
new line in Eqs. (24) and (24), or a higher $n_{max}$. The light blue curve
corresponds to the function we used for the fits in our analysis, including
operators $\epsilon,\epsilon^{\prime},T$, and $T^{\prime}$ with $n_{max}=20$.
Note that these are not fits to the data but theoretical curves with the
bootstrap values as scaling dimensions and OPE coefficients.

The light blue curves in Fig. 4 show the functions Eqs. (24) and (25) plotted
with the bootstrap values for $\Delta_{\mathcal{O}}$ and
$f_{\sigma\sigma\mathcal{O}}$, alongside our lattice data for $L=64$.
Comparison of these theoretical continuum curves with our lattice data shows
apparent consistency, justifying our use of fit functions of this form. The
other curves illustrate how the fit functions change depending on how many
operators and $n$ are included. Evidently, the operators $\epsilon$ and $T$
are most influential while including the subleading $\epsilon^{\prime}$ and
$T^{\prime}$ only slightly shift the theoretical curve, mostly at low
$t$-values. However, including these operators in fits decreased excited state
contamination of the fit parameters for the leading operators and
significantly increased the model probability of our fits, which we determined
using the Akaike Information Criterion as described in Ref. [41]. Including
even higher order operators like $\epsilon^{\prime\prime}$ led to either
unconstrained or unphysical fit parameters, however. Thus, for simultaneous
fits of $c_{0}$ and $c_{2}$, we determined the fits including
$\epsilon,\epsilon^{\prime},T$ and $T^{\prime}$ (Eqs. (24) and (25)) to be
optimal.

Because of the truncations and approximations made in our fit functions, we do
not expect them to describe our data within the entire $t$-range. Thus, we
performed a model averaging procedure as described in [41] in order to
determine a reasonable range of minimal timeslices
$(t_{0}^{min},t_{2}^{min})$. For each value of $L$, the fitting procedure can
then be summarized as follows:

  * •

First, we performed simultaneous fits to $c_{0}(t)$ and $c_{2}(t)$ with the
respective fit functions Eq. (24) and Eq. (25) with fixed maximal timeslices
$(t_{0}^{max},t_{2}^{max})$ and all possible combinations of starting times
$(t_{0}^{min},t_{2}^{min})$. Fig. 4 suggested the use of the bootstrap values
listed in Tab. 1 as initial guesses for our fit parameters to speed up the
convergence of the non-linear fitting procedure. However, we ensured that
there is no dependence of our results on the initial guesses by repeating the
fits with the first fit results as initial guesses and with checks perturbing
these initial guesses. $(t_{0}^{max},t_{2}^{max})$ were chosen such that the
relative error of the effective mass for $c_{0}$ and $c_{2}$ never exceeded
$12.5\%$. The simultaneous fits were carried out with the least squares method
using the L-BFGS-B optimization algorithm implemented in SciPy [42].

  * •

Then, we calculated the model probabilities for all fits with the Akaike
criterion according to [41].

  * •

Thereafter, we discarded fits with unphysical 111For $L=36$ some of the most
probable fits had unphysical values of
$\Delta_{\epsilon^{\prime}}<\Delta_{T}$, which we discarded. or unconstrained
parameters (i.e. parameters with errors higher than $50\%$) as well as fits
with $(t_{0}^{min},t_{2}^{min})$ far from the values with highest model
probability. After eliminating these fits, we renormalized the model
probability.

  * •

For each fit parameter, we lastly computed a weighted average over all fits
with renormalized model probability higher than $0.1\%$.

The results of the individual fits used in calculating the final model
averaged values for each $L$ are available in the Zenodo repository associated
with this project [40].

##  V Results

Figure 5: Model averaged fit results (green dots) for the scaling dimensions and OPE coefficients of the leading operators $\epsilon$ and $T$ as a function of the lattice spacing. We have performed extrapolations to $a/R\rightarrow 0$ of different functional forms (26)-(28) and including different ranges of $a/R$. The values shown as blue squares at $a/R=0$ correspond to the continuum results obtained after model averaging all different extrapolation fits. For each quantity, we also show the extrapolation fit with the highest model probability in orange - in (a) and (d), this is the finite-size scaling based fit, for (b) it is the linear extrapolation, and for (c) the quadratic fit, always including all data points. For comparison, the continuum bootstrap values are shown as dashed black lines. Quantity | Bootstrap | FSS Fit | Fit with highest $p_{model}$ | Fit Model Average | Deviation
---|---|---|---|---|---
$\Delta_{\epsilon}$ | 1.41265(36) | 1.3961(46) | 1.3961(46) | 1.3905(83) | 2.7 $\sigma$
$f_{\sigma\sigma\epsilon}^{2}$ | 1.10639(26) | 1.1625(31) | 1.1583(91) | 1.160(15) | 3.6 $\sigma$
$\Delta_{T}$ | 3 (exact) | 3.005(11) | 3.0021(86) | 3.002(12) | 0.2 $\sigma$
$f_{\sigma\sigma T}^{2}$ | 0.10636583(30) | 0.11261(58) | 0.11261(58) | 0.11280(76) | 8.5 $\sigma$
$\Delta_{\epsilon^{\prime}}$ | 3.82951(61) | 3.82(24) | 3.74(20) | 3.78(27) | 0.2 $\sigma$
$f_{\sigma\sigma\epsilon^{\prime}}^{2}$ | 0.0028102(59) | 0.0020(24) | 0.0027(20) | 0.0028(29) | 0.02 $\sigma$
$C_{T}/C_{T}^{(free)}$ | 0.946543(42) | 0.896(10) | 0.8961(81) | 0.896(11) | 3.9 $\sigma$
$f_{\sigma\sigma\epsilon}^{2}/f_{\sigma\sigma T}^{2}$ | 10.4017(24) | 10.312(79) | 10.19(12) | 10.30(16) | 0.6 $\sigma$
$C_{T}/C_{T}^{(free)}\cdot f_{\sigma\sigma\epsilon}^{2}/\Delta_{\epsilon}^{2}$ | 0.52478(30) | 0.5219(81) | 0.5219(81) | 0.528(16) | 0.2 $\sigma$
Table 2: Extrapolation results. We show the continuum values obtained with the
finite-size scaling based fits including all data points, the results from the
fits with the highest model probability, and the model averaged extrapolations
for all quantities that are also plotted in Figs. 5-8. These values are
compared to those obtained from the conformal bootstrap [24, 23], and in the
last column, we show the deviation of the model averaged extrapolation values
to the bootstrap results with respect to the statistical and fitting errors of
our lattice results.

The model averaged fit results for the scaling dimensions and OPE coefficients
of the leading operators are shown in Fig. 5 as functions of the dimensionless
lattice spacing $a/R$. The errors were calculated from the covariances of the
fits according to Eq. (17) in [41]. Thus, they do not capture possible
systematic errors in the setup of the QFE framework.

If we compare our lattice results to the values obtained by the conformal
bootstrap (shown in Tab. 1 and plotted as black dotted lines), it is apparent
that $\Delta_{\epsilon}$, $f_{\sigma\sigma\epsilon}^{2}$, and $\Delta_{T}$
trend towards the continuum bootstrap values as $a/R\rightarrow 0$. Only
$f_{\sigma\sigma T}^{2}$, the OPE coefficient associated with the energy-
momentum tensor, does not show such a clear trend.

For a quantitive comparison, we had to extrapolate our lattice results to the
continuum limit $a/R\rightarrow 0$. In Appendix B, we performed a finite-size
scaling analysis of $g(u,v)$, yielding that to leading order, the squared OPE
coefficients and scaling dimensions approach their continuum values with the
leading power of $(a/R)^{\Delta_{\epsilon^{\prime}}-3}\approx(a/R)^{0.83}$. To
extrapolate our data to the continuum, we therefore fit our lattice results
with functions of this form. Moreover, in the spirit of Taylor series
approximations to the true scaling form, we performed linear fits as well as
fits including a linear and a quadratic term. Thus, the different fit
functions we used for the continuum extrapolations are

$\displaystyle f_{\text{FSS}}$
$\displaystyle=c_{1}(\frac{a}{R})^{0.83}+c_{0},$ (26) $\displaystyle
f_{\text{linear}}$ $\displaystyle=c_{1}(\frac{a}{R})+c_{0},$ (27)
$\displaystyle f_{\text{quadratic}}$
$\displaystyle=c_{2}(\frac{a}{R})^{2}+c_{1}(\frac{a}{R})+c_{0}.$ (28)

Initially, we had also considered fits to general power laws,
$f_{\text{power}}=c_{1}(\frac{a}{R})^{c_{2}}+c_{0}$. However, our statistics
were not good enough to constrain all fit parameters. Similar problems arose
for the quadratic fit, as the only quantity for which such fits could be
constrained was $f_{\sigma\sigma\epsilon}^{2}$ due to its small error bars.
Thus, $f_{\sigma\sigma\epsilon}^{2}$ is the only quantity for which we
included such quadratic fits in determining our extrapolation results.

To account for the fact that for high $a/R$ (low $L$) we might not be close
enough to the continuum for our OPE-based fit functions to be applicable, we
also carried out all extrapolating fits excluding the first one or two data
points taken at the lowest lattice refinements ($L=24$ and $L=28$).
Subsequently, we used model averaging again to weigh the extrapolated values
of $\Delta_{\epsilon}$, $\Delta_{T}$, $f_{\sigma\sigma\epsilon}^{2}$ and
$f_{\sigma\sigma T}^{2}$ for the different fits by their model probability.
Performing this procedure, we saw that in general, the model probability
decreased the more low-$L$ data points we excluded. Moreover, we note that the
$f_{\text{FSS}}$ fit to all data points was among the two extrapolations with
the highest model probability for all quantities.

In Fig. 5, the most probable extrapolation fit (orange) as well as the model
averaged continuum extrapolation values (blue) are shown for the scaling
dimensions and OPE coefficients of the leading operators. It is apparent that
the relative errors of our lattice results are significantly higher for the
quantities associated with the $\ell=2$ operator $T$ compared to the leading
$\ell=0$ operator $\epsilon$ due to fact that we have better statistics for
$c_{0}$, in which $T$ does not appear, than for $c_{2}$. With these relatively
high error bars, all of our extrapolations for $\Delta_{T}$ (and, thus, in
particular the model averaged extrapolation) are in good agreement with the
exact value $\Delta_{T}=3$. For $\Delta_{\epsilon}$, the model averaged
extrapolation is also just in agreement with the bootstrap value with a
deviation of $2.7\sigma$, even with just taking into account statistical and
fitting errors.

Our continuum extrapolations of $f_{\sigma\sigma\epsilon}^{2}$ and
$f_{\sigma\sigma T}^{2}$, however, are significantly higher than the
respective bootstrap values. Given that the value of the OPE coefficients
could be sensitive to our normalization choice Eq. (14) for the conformal
blocks, we can explore the sensitivity to this potential source of systematic
error by studying the ratio of the OPE coefficients
$f_{\sigma\sigma\epsilon}^{2}$ and $f_{\sigma\sigma T}^{2}$ in Fig. 6, for
which we performed the same extrapolation procedure as for the other
quantities. If we take this ratio, all of our extrapolations are in
$1\,\sigma$ agreement with the bootstrap value. This suggests that the
normalization of conformal blocks at finite lattice discretization $L$
warrants further study.

Figure 6: Ratio of the model averaged fit results for the OPE coefficients of
the leading operators $\epsilon$ and $T$ (green dots) as a function of the
lattice spacing. We have performed extrapolations to $a/R\rightarrow 0$
according to (26) and (27) over different $a/R$-ranges. The extrapolation with
the highest model probability is the finite-size scaling based fit using the
lattice data for $L\geq 28$, shown in orange. The blue square at $a/R=0$
represents the continuum value obtained after model averaging all
extrapolations. For comparison, the continuum bootstrap value is shown as a
dashed black line. Figure 7: Model averaged fit results for the scaling
dimension and OPE coefficient of the subleading $\ell=0$ operator
$\epsilon^{\prime}$ (green dots) as a function of the lattice spacing. We have
performed extrapolations to $a/R\rightarrow 0$ according to (26) and (27) over
different $a/R$-ranges. The extrapolations with the highest model probability
are linear fits using the lattice data for all $a/R$, shown in orange. The
blue squares at $a/R=0$ represent the continuum values obtained after model
averaging all extrapolations. For comparison, the continuum bootstrap values
are shown as dashed black lines.

So far, we have only focused on the fit results for the quantities associated
with the leading operators $\epsilon$ and $T$ as these have the most signal
and the least excited state contamination. To demonstrate that our method also
works for the subleading operators, Fig. 7 shows our results for
$\Delta_{\epsilon^{\prime}}$ and $f_{\sigma\sigma\epsilon^{\prime}}^{2}$,
along with extrapolations to $a/R\rightarrow 0$. The error bars for those
quantities are even larger than for $\Delta_{T}$ but the extrapolations show
agreement with the bootstrap values. Here, we also do not have the issue of a
deviation of the OPE coefficient from the bootstrap value as the extrapolated
value for $f_{\sigma\sigma\epsilon^{\prime}}^{2}$ has an error of 100%. To get
more precise and reliable results for quantities associated with
$\epsilon^{\prime}$ and $T^{\prime}$, we would have to improve the statistics
and go to even higher $L$ in order to be able to include even more operators
in our fits to reduce the excited state contamination of these quantities.

Figure 8: (a) Central charge (green dots) in relation to $C_{T}^{(free)}=3/2$
as a function of lattice spacing calculated with our the model averaged fit
results for $f_{\sigma\sigma T}^{2}$ and $\Delta_{T}$, as well as
$\Delta_{\sigma}=0.518(2)$ as obtained in Ref. [20]. (b)
$C_{T}/C_{T}^{(free)}$ multiplied by our model averaged fit results results
for $f_{\sigma\sigma\epsilon}^{2}/\Delta_{\epsilon}^{2}$ as a function of the
lattice spacing. For both quantities, we have performed extrapolations to
$a/R\rightarrow 0$ according to (26) and (27) over different $a/R$-ranges. The
extrapolations with the highest model probability (orange) are a linear
extrapolation for $C_{T}/C_{T}^{(free)}$ and a finite-size scaling based
extrapolation for $C_{T}/C_{T}^{(free)}\cdot
f_{\sigma\sigma\epsilon}^{2}/\Delta_{\epsilon}^{2}$, in both cases using all
datapoints. The blue squares at $a/R=0$ represent the continuum values
obtained after model averaging all extrapolations. For comparison, the
continuum bootstrap values are shown as dashed black lines.

Lastly, we want to use our lattice results to determine the central charge of
the 3d Ising CFT, which can be calculated via

$C_{T}=\frac{\Delta_{\sigma}^{2}\Delta_{T}^{2}}{16f_{\sigma\sigma T}^{2}}.$
(29)

Our results for the central charge relative to the free value
$C_{T}^{(free)}=3/2$ can be seen in Fig. 8(a), where we used our lattice
values for $\Delta_{T}$ and $f_{\sigma\sigma T}^{2}$ to calculate $C_{T}$ at
different $L$. For $\Delta_{\sigma}$ we used our continuum value obtained in
Ref. [20]. Because our continuum extrapolation for $f_{\sigma\sigma T}^{2}$ is
significantly too high, it is not surprising that the extrapolation procedure
performed for the central charge yields a continuum value that is
significantly lower than the bootstrap value.

If we, however, take an appropriate ratio again, multiplying
$C_{T}/C_{T}^{(free)}\propto\Delta_{T}^{2}/f_{\sigma\sigma T}^{2}$ by
$f_{\sigma\sigma\epsilon}^{2}/\Delta_{\epsilon}^{2}$, the extrapolation yields
a continuum value in good agreement with the bootstrap (see Fig. 8(b)).

##  VI Free CFT analysis of lattice errors

Errors in the Monte Carlo simulations come from several sources: The
statistical Monte Carlo sampling, possible systematic errors with a non-
universal form of the action, and the finite lattice errors that need to be
removed by extrapolating to infinite lattice volume in the IR and zero lattice
spacing for the UV cut-off.

Clearly, the finite lattize errors can easily be evaluated for the exact free
scalar CFT, for which we can calculate the four-point amplitude on the same
lattices as in our current $\phi^{4}$ simulations for the interacting theory.
The comparison may be especially revealing here since the 3d Ising CFT is in
some ways a small perturbation of the free CFT exploited in the epsilon
expansion ($\epsilon=4-d$).

For the free theory, the exact continuum amplitude of the four-point function
as in Eq. 8

$\displaystyle g^{(free)}(u,v)=1$
$\displaystyle+\frac{2}{(2\cosh(t)-2\cos(\theta))^{\Delta_{\sigma}^{(free)}}}$
$\displaystyle+\frac{2}{(2\cosh(t)+2\cos(\theta))^{\Delta_{\sigma}^{(free)}}}$
(30)

is a sum over pairs of two-point functions expressed in our radial coordinates

$G^{(2)}_{\phi}(t,\cos\theta)=\dfrac{2}{\left(2\cosh
t-2\cos\theta\right)^{\Delta_{\sigma}^{(free)}}},$ (31)

where we now mean the primary operator $\sigma$ of the free CFT for which
$\Delta_{\sigma}^{(free)}=1/2$. The partial wave expansion Eq. (17) for the
free theory then has the expansion coefficients

$c^{(free)}_{j}(t)=8(2j+1)\sum_{n=0}^{\infty}e^{-(j+2n+1)t}\dfrac{2^{j}(j+n)!(2n-1)!!}{n!(2j+2n+1)!!}.$
(32)

for the even $j$ coefficients, while the odd $j$ coefficients are zero.
Moreover, as described in Appendix A, these equations may be modified to
include the finite length (or finite temperature) of the cylinder.

We computed the four-point amplitude Eq. (30) on our lattices, replacing the
free continuum propagator by the finite temperature lattice propagator,
calculated by inverting the quadratic action numerically. The goal was then to
analyze this synthetic free lattice data similarly to what was done for the
interacting case, in order to study how well our lattices enable
extrapolations to reproduce the exact free zero-temperature continuum OPE data
in Eq. (32).

Figure 9: Model averaged fit results for the scaling dimensions and operator
product expansion coefficients for the free operators $\epsilon$ ($0^{+}$) and
$T$ ($2^{+}$) as a function of the squared lattice spacing. The green results
are taken at fixed finite temperature, $L_{t}/L=16$, while the grey points in
(c) and (d) are extrapolations to zero temperature $L_{t}/L\rightarrow\infty$.
Extrapolations to the continuum (orange fit functions, blue squares at
$a/R=0$) are shown to compare with the exact continuum results at zero
temperature (black dashed line).

In particular, we calculated the finite temperature antipodal four-point
function for the free theory for the lattice refinements
$L\in\\{4,6,8,12,16,24,32\\}$ and different cylinder lengths $L_{t}$ such that
$L_{t}/L$ was between 4 and 16. Projecting this data on the Legendre
polynomials, we obtained the expansion coefficients $c^{(free)}_{j}(t)$ for
different values of $L$ and $L_{t}$. The data for the different
$c^{(free)}_{j}(t)$ is also available in the Zenodo repository associated with
this project [40]. For the two coefficients with $j=0$ and $j=2$ also
considered in the interacting theory, we then performed multi-exponential fits
of the form

$c^{fit}_{j}(t)=\sum_{k}a_{j}^{k}\left(e^{-E_{j}^{k}ta_{t}/R}+e^{-E_{j}^{k}(L_{t}-t)a_{t}/R}\right)$
(33)

in order to extract the leading coefficients and scaling exponents. In these
fits, the finite temperature is only taken into account via the addition of
the exponentials $e^{-E_{k}(L_{t}-t)a_{t}/R}$, just like it was done in the
interacting theory. We then used a model averaging procedure similar to what
was described for the interacting theory in section IV to average over fits to
different $t$-ranges and, in this free case, different numbers of included
exponentials.

Comparing Eq. (32) to (18), we could calculate the OPE coefficients
$f^{(free)}_{\sigma\sigma\epsilon}$ and $f^{(free)}_{\sigma\sigma T}$ for the
operators $\epsilon$ and $T$ of the free CFT from the leading fit parameters
$a_{j}^{0}$ and $E_{j}^{0}$ via

$\displaystyle f^{(free)2}_{\sigma\sigma\epsilon}$
$\displaystyle=\frac{a_{0}^{0}}{4^{E_{0}^{0}}}$ $\displaystyle
f^{(free)2}_{\sigma\sigma T}$ $\displaystyle=\frac{3}{8\cdot
4^{E_{2}^{0}}}\left(a_{2}^{0}-\frac{4E_{0}^{0}}{3\left(E_{0}^{0}+1\right)}a_{0}^{0}\right)$
(34)

Our lattice results for these OPE coefficients as well as the scaling
dimensions $\Delta_{\epsilon}^{(free)}$ and $\Delta_{T}^{(free)}$
corresponding to the leading exponents in $c_{0}$ and $c_{2}$, $E_{0}^{0}$ and
$E_{2}^{0}$, respectively, can be seen in Fig. 9. We plot them as a function
of $\left(\frac{a}{R}\right)^{2}$, which we know to be the leading scaling
behavior of lattice artifacts in the free case. The error bars stem from the
fitting errors. Note that we do not show the results for $L=4$ due to the fact
that for such small lattice refinements, the fitting errors are several orders
of magnitudes higher than for the other refinements.

Quantity | Exact value | $L_{t}/L=16$ | $L_{t}/L\rightarrow\infty$
---|---|---|---
$\Delta_{\epsilon}^{(free)}$ | 1 | 1.0000073(94) | -
$\Delta_{T}^{(free)}$ | 3 | 3.000030(35) | -
$f^{(free)2}_{\sigma\sigma\epsilon}$ | 2 | 1.998527(93) | 2.00014(15)
$f^{(free)2}_{\sigma\sigma T}$ | 0.09375 | 0.093681(16) | 0.093769(24)
Table 3: Extrapolation results for the free theory compared to the zero-
temperature exact values (compare Eq. (32)). For the OPE coefficients
$f^{2}_{\sigma\sigma\epsilon}$ and $f^{2}_{\sigma\sigma T}$ we list the
continuum extrapolations obtained from the data with $L_{t}/L=16$ fixed, as
well as for data which had been extrapolated to zero temperature,
$L_{t}/L\rightarrow\infty$, prior to the continuum extrapolation
$a/R\rightarrow 0$.

The results in green are taken at finite temperature with fixed $L_{t}/L=16$
just as in the interacting theory. Extrapolating these to the continuum by
fitting $f=c_{0}+c_{1}\left(\frac{a}{R}\right)^{2}$ leads to good agreement
with the analytic values in the case of the scaling dimensions, which we
expect to be $\Delta_{\epsilon}^{(free)}=1$ and $\Delta_{T}^{(free)}=3$ as
$a/R\rightarrow 0$.

The OPE coefficients, however, are affected by finite temperature effects.
Thus, for the fixed finite ratio $L_{t}/L=16$ corresponding to a finite
temperature, the extrapolations to $a/R\rightarrow 0$ with
$f^{(free)2}_{\sigma\sigma\epsilon}=1.998527(93)$ and
$f^{(free)2}_{\sigma\sigma T}=0.093681(16)$ slightly miss the zero temperature
continuum expectations from Eq. (32) corresponding to
$f^{(free)2}_{\sigma\sigma\epsilon}=2$ and $f^{(free)2}_{\sigma\sigma
T}=0.09375$. However, these deviations due to wraparound effects only amount
to 0.07%.

We can eliminate this systematic error by carrying out a double limiting
procedure, first extrapolating the lattice results to zero temperature
$L_{t}/L\rightarrow\infty$ via fit function $f(L_{t})=f+p_{1}\
e^{-p_{2}L_{t}}$ for each $L$, yielding the grey data in Fig. 9, and then
taking the limit $a\rightarrow 0$. Fig. 9 shows that if we perform this
procedure, our results perfectly agree with the continuum zero temperature
expectations.

From these results, we can draw two main conclusions. Firstly, as we know the
scaling with the lattice spacing exactly, we can perfectly extrapolate the
free theory data calculated on our simplicial lattices data to remove the IR
and UV cutoff, yielding the free continuum values. This might be an incentive
to study the cutoff dependence of the interacting theory to higher than just
leading order in order to improve extrapolations.

Secondly, there are wraparound effects influencing our data if we only carry
out calculations at finite T (finite $L_{t}/L$) and compare to $T\rightarrow
0$ ($L_{t}/L\rightarrow\infty$) continuum expectations. However, the errors
associated with these effects are around two orders of magnitude smaller than
the deviations of our OPE coefficient results from the bootstrap values that
we see in the interacting theory. Thus, we conclude that the improvement from
carrying out a limit $L_{t}/L\rightarrow\infty$ in the interacting theory
would be negligible compared to the effect of the systematic errors in the
interacting theory that we could not study with the free case.

##  VII Conclusion and Outlook

The Quantum Finite Element program seeks to extend Euclidean Monte Carlo
lattice methods to non-perturbative quantum field theory on curved manifolds.
In particular, the present goal is to demonstrate its advantage when applied
to radial quantization for $d\geq 3$ conformal field theory. For the prototype
lattice $\phi^{4}$ example, we have demonstrated that we can extract scaling
dimensions and OPE coefficients for leading as well as subleading operators of
the 3d Ising CFT from the four-point function on a
${\mathbb{R}}\times{\mathbb{S}}^{d}$ manifold. In particular, this is the
first time, to our knowledge, that the OPE coefficient $f_{\sigma\sigma T}$
has been determined from lattice Monte Carlo calculations.

Continuum extrapolations of our lattice simulations for the scaling dimensions
of operators $\epsilon$, $T$ as well as the subleading operator
$\epsilon^{\prime}$ agree with the values obtained with the conformal
bootstrap (and the exact value in the case of the energy-momentum tensor)
within the statistical and fitting errors. Extrapolations for the OPE
coefficients $f_{\sigma\sigma\epsilon}$ and $f_{\sigma\sigma T}$ deviate more
from the bootstrap values than is covered by these errors so that also our
value for the central charge is significantly too low in violation of the
bootstrap lower bound [22]. However, ratios of these quantities extrapolate to
values that are consistent with values obtained by the other two methods that
were able to extract both $f_{\sigma\sigma\epsilon}$ and $f_{\sigma\sigma T}$,
namely the conformal bootstrap and the fuzzy sphere, as is shown in the
comparison Tab. 4.

Quantity | This work | Bootstrap | Fuzzy sphere
---|---|---|---
$f_{\sigma\sigma\epsilon}^{2}/f_{\sigma\sigma T}^{2}$ | 10.30(16) | 10.402(24) | 10.53(23)
$C_{T}/C_{T}^{(free)}\cdot f_{\sigma\sigma\epsilon}^{2}/\Delta_{\epsilon}^{2}$ | 0.528(16) | 0.52478(13) | 0.5310(86)
Table 4: Comparison of the ratios obtained in this work with the corresponding
values obtained with the conformal bootstrap [23, 24] and the Fuzzy sphere
[27].

The method employed here follows a sequence of Quantum Finite Element
developments. First, there is the requirement to introduce UV counterterms in
2d on ${\mathbb{S}}^{2}$ and in 3d on ${\mathbb{R}}\times{\mathbb{S}}^{d-1}$
to modify the simplicial finite element method (FEM) to reach the Ising
critical surface. Next, the importance of the Ricci term in 3d was seen as it
dramatically improved the approach to continuum for the scalar two-point
function. While this simplicial QFE action has enabled us to study the OPE
expansion, there are clearly additional improvements to pursue, as is common
practice in lattice field theory.

Comparison of the free scalar theory on the simplicial lattice vs. the
continuum free CFT reveals two issues. First, the finite temperature (aka IR
finite length of the cylinder) is probably not a significant source of
systematic error. Second, the extrapolation to zero lattice spacing (i.e.
removing the UV cut-off) for the free CFT with $O(a^{2})$ scaling is under
good control. For the strongly-interacting CFT, however, the analytic
dependence on the cutoff is a more serious problem to address and finite-size
scaling analysis offers some guidance.

Within the current QFE framework presented here, there are still some
potential sources of systematic error that could be explored further. A slight
mistuning of the bare mass parameter $\mu_{0}^{2}$ could potentially lead to a
systematic shift in estimated OPE coefficients. This can only be checked by
new calculations at even larger values of $L$. Similarly, all calculations
were performed at only a single bare $\lambda_{0}$ with no attempt to
extrapolate to the bare lattice $\lambda_{0}\to 0$. Also, our determinations
of scaling dimensions and OPE coefficients were actually variational estimates
limited by our inability to model the contributions of higher primaries:
($\epsilon^{\prime\prime}$, $T^{\prime\prime}$, $\ell\geq 4$, etc). Any of
these studies would require a significant investment in software development
for increased parallelization of the code as well as at least an order of
magnitude increase in computational resources. Before making such an
investment, we would like to consider whether this framework is the best
option available.

Beyond cut-off effects, we continue to explore the validity of using
perturbative counterterms in super-renormalizable theories. Recall that our
counterterms only correct the mass parameter and do not modify the classical
FEM Laplace-Beltrami operator [44, 45] in Eq. (19). However, in a recent study
of the 2d Ising model on an affine lattice [46], equivalent to lattice
$\phi^{4}$ at bare lattice $\lambda_{0}\rightarrow\infty$, we demonstrated
that the critical theory required a modification of the kinetic term away from
the classical FEM operator. This provides yet further evidence that the way to
continue using our current simplicial action in Eq. 21 with perturbative
counterterms is to fix the dimensional renormalized coupling
($\lambda_{R}=\mathcal{O}(\lambda_{0}/a)$) instead of the bare $\lambda_{0}$
as we extrapolate to zero lattice spacing. This should reduce the quantum
correction of FEM Laplace-Beltrami operator to $O(a)$.

##  VIII Acknowledgements

The authors are pleased to acknowledge that the computational work reported on
in this paper was performed on the Shared Computing Cluster, which is
administered by Boston University’s Research Computing Services
(www.bu.edu/tech/support/research/). We thank Casey Berger and Andrew D.
Gasbarro, whose calculations this work builds upon. This work is supported by
the U.S. Department of Energy (DOE) under Award No. DE-SC0019061 for GTF
(Yale), Award No. DE-SC0015845 for RCB, and Award No. DE-SC0019139 for EKO.
AMEG acknowledges support from Cusanuswerk Bischöfliche Studienförderung and
the Baden-Württemberg Stiftung. This document was prepared by the Quantum
Finite Elements collaboration using the resources of the Fermi National
Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of
Science, Office of High Energy Physics HEP User Facility. Fermilab is managed
by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-
AC02-07CH11359.

##  Appendix A The free case

In the special case of the free theory, the propagator
$G^{(2)}_{\sigma}(t,\cos\theta)\equiv\langle\sigma(t_{1},\bm{n}_{1})\sigma(t_{2},\bm{n}_{2})\rangle$
for a free scalar field $\sigma(t,\bm{n})$ on $\mathbb{R}\times S^{2}$ is

$G^{(2)}_{\sigma}(t,\cos\theta)=\dfrac{2}{\left(2\cosh
t-2\cos\theta\right)^{\Delta_{\sigma}^{(free)}}}$ (35)

where $t=|t_{2}-t_{1}|$ and $\cos\theta=\bm{n}_{1}\cdot\bm{n}_{2}$. Expanding
as a series in Legendre polynomials this is

$G^{(2)}_{\sigma}(t,\cos\theta)=\sum_{j=0}^{\infty}2e^{-(j+1/2)t}P_{j}(\cos\theta)\;.$
(36)

where we have substituted the scaling exponent $\Delta_{\sigma}^{(free)}=1/2$.

The free antipodal four-point function is simply a sum of $s$, $t$, and
$u$-channels. Following the same notation used above this gives

$G^{(4)}_{\sigma}(t,\cos\theta)=G^{(2)}_{\sigma}(0,-1)^{2}+G^{(2)}_{\sigma}(t,\cos\theta)^{2}\\\
+G^{(2)}_{\sigma}(t,-\cos\theta)^{2}\;.$ (37)

Note that we have conveniently normalized the 2-point function (36) so that
the $s$-channel diagram is 1. Thus, the antipodal 4-point function in this
normalization already corresponds to the conformally invariant amplitude $g$.
Inserting Eq. (36) we obtain the somewhat complicated expression

$G^{(4)}_{\sigma}(t,\cos\theta)=1+4\left(\sum_{j=0}^{\infty}e^{-(j+1/2)t}P_{\ell}(\cos\theta)\right)^{2}\\\
+4\left(\sum_{j=0}^{\infty}e^{-(j+1/2)t}P_{j}(-\cos\theta)\right)^{2}$ (38)

which can again be expressed as a series in Legendre polynomials

$G^{(4)}_{\sigma}(t,\cos\theta)=1+\sum_{j=0}^{\infty}c^{(free)}_{j}(t)P_{j}(\cos\theta)$
(39)

where all of the odd $j$ coefficients are zero, and the even $j$ coefficients
are

$c^{(free)}_{j}(t)=8(2j+1)\sum_{n=0}^{\infty}e^{-(j+2n+1)t}\dfrac{2^{j}(j+n)!(2n-1)!!}{n!(2j+2n+1)!!}.$
(40)

As required, we recover the partial wave expansion for a free scalar four-
point function.

Finally, we consider corrections to (39) due to wraparound effects on a
periodic cylinder with length $L_{t}$ (i.e. finite temperature effects). In
the free case, the thermal 2-point function $\tilde{G}^{(2)}_{\sigma}$ in the
interval $0\leq t\leq L_{t}$ can be expressed as a sum over an infinite set of
image charges.

$\tilde{G}^{(2)}_{\sigma}(t,\cos\theta)=\sum_{n=-\infty}^{\infty}G^{(2)}_{\phi}(t+nL_{t},\cos\theta)\\\
\\\
=4\sum_{j=0}^{\infty}\dfrac{\cosh\left[(\Delta_{\sigma}^{(free)}+j)(t-L_{t}/2)\right]}{\sinh\left[(\Delta_{\sigma}^{(free)}+j)L_{t}/2\right]}P_{j}(\cos\theta)$
(41)

The four-point function can be evaluated as before, and to lowest order, the
finite temperature Legendre coefficients are related to the zero-temperature
coefficients by $\tilde{c}^{(free)}_{j}(t)\simeq
c^{(free)}_{j}(t)+c^{(free)}_{j}(L_{t}-t)$. The leading corrections occur near
$t=L_{t}/2$, giving

$\delta\tilde{c}^{(free)}_{0}(t\sim
L_{t}/2)\simeq\dfrac{8}{\sinh(\Delta_{\epsilon}^{(free)}L_{t}/2)}$ (42)

for $j=0$ and

$\delta\tilde{c}^{(free)}_{j}(t\sim L_{t}/2)\simeq
16\dfrac{\cosh\left[(\Delta_{\epsilon}^{(free)}+j-1)(t-L_{t}/2)\right]}{\sinh\left[(\Delta_{\epsilon}^{(free)}+j)L_{t}/2\right]}$
(43)

for all even $j\geq 2$.

##  Appendix B Finite-size scaling analysis for the interacting theory

To determine which fit function to use for extrapolating the lattice results
to the continuum $a/R\rightarrow 0$, we carried out a finite-size scaling
analysis based on Ref. [47].

###  B.1 Scaling of the free energy

The free energy is defined as

$\mathcal{F}(g_{\sigma},g_{\epsilon},\\{g_{\omega},\cdots\\},a/L)=\log Z$ (44)

where

$Z=\int D\phi e^{-S+h\cdot\phi}.$ (45)

Here, $g_{\sigma},g_{\epsilon}$ are relevant and $g_{\omega},\cdots$
irrelevant couplings of the action near the Wilson-Fisher fixed point, not the
bare couplings $h,\mu_{0},\lambda_{0}$ in the action. $a$ is the lattice
spacing, $L$ is the characteristic “finite size” of our system, and $a/L$ a
dimensionless ratio.

If we change the scale by a factor $\lambda$, the free energy renormalizes as
follows:

$\displaystyle\mathcal{F}(g_{\sigma},$ $\displaystyle
g_{\epsilon},\\{g_{\omega},\cdots\\},a/L)$ (46) $\displaystyle=$
$\displaystyle\lambda^{-d}F(\lambda^{y_{\sigma}}g_{\sigma},\lambda^{y_{\epsilon}}g_{\epsilon},\\{\lambda^{y_{\omega}}g_{\omega},\cdots\\},\lambda
a/L)+G(g_{\sigma},g_{\epsilon})$

Here, $F$ is the scaling part of $\mathcal{F}$, whereas $G$ is regular. The
exponents $y_{\mathcal{O}}$ are related to the scaling dimensions of the
respective operators via $y_{\mathcal{O}}=d-\Delta_{O}$

###  B.2 $n$-point functions

We calculate connected $n$-point functions of the scalar field $\phi$ from
$\mathcal{F}$ by taking derivatives with respect to the (bare) external field
$h$ and subsequently setting $h=0$:

$\displaystyle\left.\langle\phi(x_{1})\cdots\phi(x_{n})\rangle\right|_{conn.}=\left.\frac{\delta}{\delta
h(x_{1})}\cdots\frac{\delta}{\delta h(x_{n})}\mathcal{F}\right|_{h=0}$ (47)

Rewriting the derivatives with respect to $g_{\sigma}$, we get for the two-
point function

$\displaystyle\begin{split}\langle\phi(x_{1})&\phi(x_{2})\rangle|_{conn.}\\\
&\left.=\int_{x,y}\frac{\delta^{2}\mathcal{F}}{\delta g_{\sigma}(x)\delta
g_{\sigma}(y)}\frac{\delta g_{\sigma}(x)}{\delta h(x_{1})}\frac{\delta
g_{\sigma}(y)}{\delta h(x_{2})}\right|_{h=0}\end{split}$ (48)

where we have used that $g_{\sigma}$ is an odd function of $h$

$g_{\sigma}=h\alpha_{1}+h^{3}\frac{\alpha_{3}}{3!}+h^{5}\frac{\alpha_{5}}{5!}+\cdots$
(49)

such that terms like $\frac{\delta^{2}g_{\sigma}(x)}{\delta h(x_{1})\delta
h(x_{2})}$ vanish. For the four-point function, we get

$\displaystyle\langle\phi($ $\displaystyle
x_{1})\phi(x_{2})\phi(x_{3})\phi(x_{4})\rangle|_{conn.}$ $\displaystyle=$
$\displaystyle\int\limits_{x,y,u,w}\frac{\delta^{4}\mathcal{F}}{\delta
g_{\sigma}(x)\delta g_{\sigma}(y)g_{\sigma}(u)\delta g_{\sigma}(w)}$
$\displaystyle\qquad\qquad\times\frac{\delta g_{\sigma}(x)}{\delta
h(x_{1})}\frac{\delta g_{\sigma}(y)}{\delta h(x_{2})}\frac{\delta
g_{\sigma}(u)}{\delta h(x_{3})}\frac{\delta g_{\sigma}(w)}{\delta h(x_{4})}$
$\displaystyle+$
$\displaystyle\int\limits_{x,y}\frac{\delta^{2}\mathcal{F}}{\delta
g_{\sigma}(x)\delta g_{\sigma}(y)}\biggl{(}\frac{\delta g_{\sigma}(x)}{\delta
h(x_{1})}\frac{\delta^{3}g_{\sigma}(y)}{\delta h(x_{2})\delta h(x_{3})\delta
h(x_{4})}$ $\displaystyle\qquad\qquad+\frac{\delta g_{\sigma}(x)}{\delta
h(x_{2})}\frac{\delta^{3}g_{\sigma}(y)}{\delta h(x_{1})\delta h(x_{3})\delta
h(x_{4})}$ $\displaystyle\qquad\qquad+\frac{\delta g_{\sigma}(x)}{\delta
h(x_{3})}\frac{\delta^{3}g_{\sigma}(y)}{\delta h(x_{1})\delta h(x_{2})\delta
h(x_{4})}$ $\displaystyle\qquad\qquad+\frac{\delta g_{\sigma}(x)}{\delta
h(x_{4})}\frac{\delta^{3}g_{\sigma}(y)}{\delta h(x_{1})\delta h(x_{2})\delta
h(x_{3})}\biggl{)}$

Derivatives of the free energy with respect to $g_{\sigma}$ in the limit of
$h\rightarrow 0$ and thus $g_{\sigma}\rightarrow 0$ scale as

$\displaystyle\frac{\partial^{k}\mathcal{F}}{\partial
g_{\sigma}^{k}}=\mathcal{F}^{(k)}(g_{\epsilon},\\{g_{\omega},\cdots\\},a/L)$
(50)
$\displaystyle=\left(\frac{L}{a}\right)^{ky_{\sigma}-d}F^{(k)}\left(\left(\frac{L}{a}\right)^{y_{\epsilon}}g_{\epsilon},\left\\{\left(\frac{L}{a}\right)^{y_{\omega}}g_{\omega},\cdots\right\\},1\right)$
$\displaystyle\quad+G^{(k)}(g_{\epsilon})$

Near the critical point, we can expand these derivatives in Taylor series as
follows:

$\displaystyle F^{(k)}$ (51)
$\displaystyle=a_{k0}+a_{k1}(g_{\epsilon}-g_{\epsilon}^{*})\left(\frac{L}{a}\right)^{y_{\epsilon}}+a_{k2}(g_{\epsilon}-g_{\epsilon}^{*})^{2}\left(\frac{L}{a}\right)^{2y_{\epsilon}}+\cdots$
$\displaystyle\quad+b_{k1}(g_{\omega}-g_{\omega}^{*})\left(\frac{L}{a}\right)^{y_{\omega}}+b_{k2}(g_{\omega}-g_{\omega}^{*})^{2}\left(\frac{L}{a}\right)^{2y_{\omega}}+\cdots$

and

$\displaystyle
G^{(k)}=c_{k0}+c_{k1}(g_{\epsilon}-g_{\epsilon}^{*})+c_{k2}(g_{\epsilon}-g_{\epsilon}^{*})^{2}+\cdots$
(52)

where $g^{*}$ denotes the critical values of the couplings at the Wilson-
Fisher fixed point.

Putting everything together, we see that the two-point function scales as

$\displaystyle\left.\langle\phi\phi\rangle\right|_{conn.}\sim\left(\frac{L}{a}\right)^{2y_{\sigma}-2d}\alpha_{1}^{2}\left[F^{(2)}+\left(\frac{L}{a}\right)^{d-2y_{\sigma}}G^{(2)}\right]$
(53)

near the critical point, whereas the connected four-point function scales as

$\displaystyle\left.\langle\phi\phi\phi\phi\rangle\right|_{conn.}\sim\left(\frac{L}{a}\right)^{4y_{\sigma}-4d}\alpha_{1}^{4}\left[F^{(4)}+\left(\frac{L}{a}\right)^{d-4y_{\sigma}}G^{(4)}\right]$
$\displaystyle+4\left(\frac{L}{a}\right)^{2y_{\sigma}-4d}\alpha_{1}\alpha_{3}\left[F^{(2)}+\left(\frac{L}{a}\right)^{d-2y_{\sigma}}G^{(2)}\right]$
(54)

By tuning $\mu_{0}^{2}$ in our action to the critical surface, we eliminate
all terms $\propto(g_{\epsilon}-g_{\epsilon}^{*}=0)$, which diverge as
$L\rightarrow\infty$. The remaining terms then yield

$\displaystyle\left.\langle\phi\phi\rangle\right|_{conn.}$ (55)
$\displaystyle\sim\left(\frac{L}{a}\right)^{-2\Delta_{\sigma}}\alpha_{1}^{2}\left[a_{20}+b_{21}(g_{\omega}-g_{\omega}^{*})\left(\frac{L}{a}\right)^{d-\Delta_{\omega^{\prime}}}+\cdots\right.$
$\displaystyle\left.\quad+\left(\frac{L}{a}\right)^{2\Delta_{\sigma}-d}c_{20}\right]$

and

$\displaystyle\left.\langle\phi\phi\phi\phi\rangle\right|_{conn.}$ (56)
$\displaystyle\sim\left(\frac{L}{a}\right)^{-4\Delta_{\sigma}}\alpha_{1}^{4}\left[a_{40}+b_{41}(g_{\omega}-g_{\omega}^{*})\left(\frac{L}{a}\right)^{d-\Delta_{\omega^{\prime}}}+\cdots\right.$
$\displaystyle+\left.\left(\frac{L}{a}\right)^{-3d+4\Delta_{\sigma}}c_{40}\right]+4\left(\frac{L}{a}\right)^{-2\Delta_{\sigma}-2d}\alpha_{1}\alpha_{3}$
$\displaystyle\,\,\times\left[a_{20}+b_{21}(g_{\omega}-g_{\omega}^{*})\left(\frac{L}{a}\right)^{d-\Delta_{\omega^{\prime}}}+\cdots+\left(\frac{L}{a}\right)^{2\Delta_{\sigma}-d}c_{20}\right]$

We now have to take the ratio
$\frac{\left.\langle\phi\phi\phi\phi\rangle\right|_{conn.}}{\langle\phi\phi\rangle\langle\phi\phi\rangle}$
of Eq. (56) and Eq. (55) and expand in $a/L$. Using that for the $d=3$ Ising
model, the leading irrelevant operator $\omega$ corresponds to the conformal
primary $\epsilon^{\prime}$ and comparing with the values for
$\Delta_{\sigma}$ and $\Delta_{\epsilon^{\prime}}$ in [23], we find that to
leading order

$\displaystyle\frac{\left.\langle\phi\phi\phi\phi\rangle\right|_{conn.}}{\langle\phi\phi\rangle\langle\phi\phi\rangle}$
$\displaystyle=g(u,v)-1$ $\displaystyle\sim
C_{0}+C_{1}\left(\frac{a}{L}\right)^{\Delta_{\epsilon^{\prime}}-3}+\cdots$
(57)

such that the leading finite-size scaling exponent is
$\Delta_{\epsilon^{\prime}}-3\approx 0.83$. In our case, the characteristic
”finite size” of our system is the radius of the sphere $R$. Thus, we expect
the OPE coefficients from our fits at finite lattice spacing to receive
corrections

$f_{\sigma\sigma\mathcal{O}}^{2}(\frac{a}{R})=f_{\sigma\sigma\mathcal{O}}^{2}+c_{\mathcal{O}}^{\prime}\left(\frac{a}{R}\right)^{\Delta_{\epsilon^{\prime}}-3}+\cdots$
(58)

to leading order. Taking the logarithm, we find that also

$\Delta_{\mathcal{O}}(\frac{a}{R})=\Delta_{\mathcal{O}}+c_{\mathcal{O}}\left(\frac{a}{R}\right)^{\Delta_{\epsilon^{\prime}}-3}+\cdots$
(59)

Note, $C_{1}$ in Eq. (B.2) as the leading correction to scaling is a function
of the irrelevant couplings and its value can be set to zero by tuning the
bare couplings [48], thereby increasing the rate of convergence of the
discrete theory to the continuum limit. This is often referred to as non-
perturbative improvement of the discrete action.

## References

  * Andrews [1869] T. Andrews, The bakerian lecture: On the continuity of the gaseous and liquid states of matter, Philosophical Transactions of the Royal Society of London 159, 575 (1869).
  * Wilson [1971] K. G. Wilson, Renormalization group and critical phenomena. i. renormalization group and the kadanoff scaling picture, Phys. Rev. B 4, 3174 (1971).
  * Parola and Reatto [1995] A. Parola and L. Reatto, Liquid state theories and critical phenomena, Advances in Physics 44, 211 (1995), https://doi.org/10.1080/00018739500101536 .
  * Anisimov [1991] M. Anisimov,  _Critical Phenomena in Liquids and Liquid Crystals_ (Gordon and Breach Science Publishers, 1991).
  * Pelissetto and Vicari [2002] A. Pelissetto and E. Vicari, Critical phenomena and renormalization group theory, Phys. Rept. 368, 549 (2002), arXiv:cond-mat/0012164 .
  * Swendsen and Wang [1987] R. H. Swendsen and J.-S. Wang, Nonuniversal critical dynamics in Monte Carlo simulations, Phys. Rev. Lett. 58, 86 (1987).
  * Wolff [1989] U. Wolff, Collective Monte Carlo Updating for Spin Systems, Phys. Rev. Lett. 62, 361 (1989).
  * Hasenbusch [2010] M. Hasenbusch, Finite size scaling study of lattice models in the three-dimensional Ising universality class, Phys. Rev. B 82, 174433 (2010), arXiv:1004.4486 [cond-mat.stat-mech] .
  * Polyakov [1970] A. M. Polyakov, Conformal symmetry of critical fluctuations, JETP Lett. 12, 381 (1970).
  * Rose  _et al._ [2022] F. Rose, C. Pagani, and N. Dupuis, Operator product expansion coefficients from the nonperturbative functional renormalization group, Phys. Rev. D 105, 065020 (2022).
  * Hasenbusch [2018] M. Hasenbusch, Two- and three-point functions at criticality: Monte carlo simulations of the improved three-dimensional blume-capel model, Phys. Rev. E 97, 012119 (2018).
  * Herdeiro [2017] V. Herdeiro, Numerical estimation of structure constants in the three-dimensional ising conformal field theory through markov chain uv sampler, Phys. Rev. E 96, 033301 (2017).
  * Caselle  _et al._ [2015] M. Caselle, G. Costagliola, and N. Magnoli, Numerical determination of the operator-product-expansion coefficients in the 3d ising model from off-critical correlators, Phys. Rev. D 91, 061901 (2015).
  * Costagliola [2016] G. Costagliola, Operator product expansion coefficients of the 3d ising model with a trapping potential, Phys. Rev. D 93, 066008 (2016).
  * Cardy [1985] J. L. Cardy, Universal amplitudes in finite-size scaling: generalisation to arbitrary dimensionality, J. Phys. A 18, L757 (1985).
  * Fubini  _et al._ [1973] S. Fubini, A. J. Hanson, and R. Jackiw, New approach to field theory, Phys.Rev. D7, 1732 (1973).
  * Brower  _et al._ [2016] R. C. Brower, G. Fleming, A. Gasbarro, T. Raben, C.-I. Tan, and E. Weinberg, Quantum Finite Elements for Lattice Field Theory, PoS LATTICE2015, 296 (2016), arXiv:1601.01367 [hep-lat] .
  * Brower  _et al._ [2018] R. C. Brower, M. Cheng, E. S. Weinberg, G. T. Fleming, A. D. Gasbarro, T. G. Raben, and C.-I. Tan, Lattice $\phi^{4}$ field theory on Riemann manifolds: Numerical tests for the 2-d Ising CFT on $\mathbb{S}^{2}$, Phys. Rev. D 98, 014502 (2018), arXiv:1803.08512 [hep-lat] .
  * Brower  _et al._ [2017] R. C. Brower, E. S. Weinberg, G. T. Fleming, A. D. Gasbarro, T. G. Raben, and C.-I. Tan, Lattice dirac fermions on a simplicial riemannian manifold, Phys. Rev. D 95, 114510 (2017).
  * Brower  _et al._ [2021] R. C. Brower, G. T. Fleming, A. D. Gasbarro, D. Howarth, T. G. Raben, C.-I. Tan, and E. S. Weinberg, Radial lattice quantization of 3D $\phi$4 field theory, Phys. Rev. D 104, 094502 (2021), arXiv:2006.15636 [hep-lat] .
  * El-Showk  _et al._ [2012] S. El-Showk, M. F. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin, and A. Vichi, Solving the 3D Ising Model with the Conformal Bootstrap, Phys. Rev. D 86, 025022 (2012), arXiv:1203.6064 [hep-th] .
  * El-Showk  _et al._ [2014] S. El-Showk, M. F. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin, and A. Vichi, Solving the 3d Ising Model with the Conformal Bootstrap II. c-Minimization and Precise Critical Exponents, J. Stat. Phys. 157, 869 (2014), arXiv:1403.4545 [hep-th] .
  * Simmons-Duffin [2017] D. Simmons-Duffin, The Lightcone Bootstrap and the Spectrum of the 3d Ising CFT, JHEP 03, 086, arXiv:1612.08471 [hep-th] .
  * Reehorst [2022] M. Reehorst, Rigorous bounds on irrelevant operators in the 3d ising model CFT, Journal of High Energy Physics 2022, 10.1007/jhep09(2022)177 (2022).
  * Glück  _et al._ [2023] A.-M. E. Glück, G. T. Fleming, R. C. Brower, V. Ayyar, E. K. Owen, T. G. Raben, and C.-I. Tan, Computing the Central Charge of the 3D Ising CFT Using Quantum Finite Elements, PoS LATTICE2022, 370 (2023).
  * Han  _et al._ [2023] C. Han, L. Hu, W. Zhu, and Y.-C. He, Conformal four-point correlators of the 3D Ising transition via the quantum fuzzy sphere, (2023), arXiv:2306.04681 [cond-mat.stat-mech] .
  * Hu  _et al._ [2023] L. Hu, Y.-C. He, and W. Zhu, Operator product expansion coefficients of the 3d ising criticality via quantum fuzzy spheres, Phys. Rev. Lett. 131, 031601 (2023).
  * Di Francesco  _et al._ [1997] P. Di Francesco, P. Mathieu, and D. Senechal,  _Conformal Field Theory_, Graduate Texts in Contemporary Physics (Springer-Verlag, New York, 1997).
  * Rychkov [2016] S. Rychkov, EPFL Lectures on Conformal Field Theory in $D\geq 3$ Dimensions, (2016), arXiv:1601.05000 [hep-th] .
  * Dolan and Osborn [2001] F. A. Dolan and H. Osborn, Conformal four point functions and the operator product expansion, Nucl. Phys. B 599, 459 (2001), arXiv:hep-th/0011040 .
  * Dolan and Osborn [2011] F. A. Dolan and H. Osborn, Conformal Partial Waves: Further Mathematical Results, (2011), arXiv:1108.6194 [hep-th] .
  * Hogervorst and Rychkov [2013] M. Hogervorst and S. Rychkov, Radial Coordinates for Conformal Blocks, Phys. Rev. D 87, 106004 (2013), arXiv:1303.1111 [hep-th] .
  * Costa  _et al._ [2016] M. S. Costa, T. Hansen, J. a. Penedones, and E. Trevisani, Radial expansion for spinning conformal blocks, JHEP 07, 057, arXiv:1603.05552 [hep-th] .
  * Brower and Tamayo [1989] R. C. Brower and P. Tamayo, Embedded Dynamics for $\phi^{4}$ Theory, Phys. Rev. Lett. 62, 1087 (1989).
  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21, 1087 (1953).
  * Adler [1981] S. L. Adler, An Overrelaxation Method for the Monte Carlo Evaluation of the Partition Function for Multiquadratic Actions, Phys. Rev. D 23, 2901 (1981).
  * Whitmer [1984] C. Whitmer, Overrelaxation methods for Monte Carlo simulations of quadratic and multiquadratic actions, Phys. Rev. D 29, 306 (1984).
  * Brown and Woch [1987] F. R. Brown and T. J. Woch, Overrelaxed Heat Bath and Metropolis Algorithms for Accelerating Pure Gauge Monte Carlo Calculations, Phys. Rev. Lett. 58, 2394 (1987).
  * Creutz [1987] M. Creutz, Overrelaxation and Monte Carlo Simulation, Phys. Rev. D 36, 515 (1987).
  * Glück  _et al._ [2023] A.-M. E. Glück, V. Ayyar, R. C. Brower, G. T. Fleming, E. K. Owen, T. G. Raben, and C.-I. Tan, Dataset for ”The Operator Product Expansion for Radial Lattice Quantization of 3D $\phi^{4}$ Theory” (2023).
  * Jay and Neil [2021] W. I. Jay and E. T. Neil, Bayesian model averaging for analysis of lattice field theory results, Phys. Rev. D 103, 114502 (2021), arXiv:2008.01069 [stat.ME] .
  * Byrd  _et al._ [1995] R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, A limited memory algorithm for bound constrained optimization, SIAM Journal on Scientific Computing 16, 1190 (1995).
  * Note [1] For $L=36$ some of the most probable fits had unphysical values of $\Delta_{\epsilon^{\prime}}<\Delta_{T}$, which we discarded.
  * Strang and Fix [2008] G. Strang and G. Fix,  _An Analysis of the Finite Element Method 2nd Edition_, 2nd ed. (Wellesley-Cambridge, 2008).
  * Desbrun  _et al._ [2005] M. Desbrun, A. N. Hirani, M. Leok, and J. E. Marsden, Discrete Exterior Calculus, ArXiv Mathematics e-prints (2005), math/0508341 .
  * Brower and Owen [2023] R. C. Brower and E. K. Owen, Ising model on the affine plane, Physical Review D 108, 10.1103/physrevd.108.014511 (2023).
  * Blote  _et al._ [1995] H. W. J. Blote, E. Luijten, and J. R. Heringa, Ising universality in three dimensions: a Monte Carlo study, J. Phys. A 28, 6289 (1995), arXiv:cond-mat/9509016 .
  * Hasenbusch  _et al._ [1999] M. Hasenbusch, K. Pinn, and S. Vinti, Critical exponents of the three-dimensional Ising universality class from finite-size scaling with standard and improved actions, Phys. Rev. B 59, 11471 (1999), arXiv:hep-lat/9806012 .