SciPost Phys. Core 7, 028 (2024)

Quantum Monte Carlo simulation of the
3D Ising transition on the fuzzy sphere
Johannes Stephan Hofmann1⋆ , Florian Goth2 , Wei Zhu3 ,Yin-Chen He4 and Emilie Huffman4†
1 Department of Condensed Matter Physics, Weizmann Institute of Science,
Rehovot, 76100, Israel
2 Institute for Theoretical Physics, University of Würzburg, Würzburg, 97074, Germany
3 School of Science, Westlake University, Hongzhou, 310024, P. R. China
4 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada⋆ johannes-stephan.hofmann@weizmann.ac.il ,† ehuffman@perimeterinstitute.ca

Abstract
We present a numerical quantum Monte Carlo (QMC) method for simulating the 3D
phase transition on the recently proposed fuzzy sphere [1]. By introducing an additional SU(2) layer degree of freedom, we reformulate the model into a form suitable
for sign-problem-free QMC simulation. From the finite-size-scaling, we show that thisQMC-friendly model undergoes a quantum phase transition belonging to the 3D Ising
universality class, and at the critical point we compute the scaling dimensions from thestate-operator correspondence, which largely agrees with the prediction from the conformal field theory. These results pave the way to construct sign-problem-free models
for QMC simulations on the fuzzy sphere, which could advance the future study on moresophisticated criticalities.
Copyright J. S. Hofmann et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

05-01-2024
Received
Accepted 25-04-2024
Check forupdatesPublished 09-05-2024
doi:10.21468/SciPostPhysCore.7.2.028

Contents
1 Introduction

2

2 Model and method
2.1 Review of fuzzy sphere regularization
2.2 The density form of interaction
2.3 Four component fuzzy sphere model

3
3
5
5

3 Results
3.1 QMC simulations
3.2 Finite-size-scaling
3.3 Dimensionless two-point correlator
3.4 Energy gaps and state-operator correspondence

6
6
8
8
10

4 Conclusions

10
1

SciPost Phys. Core 7, 028 (2024)

A Order parameter

12

B Extracting energy gaps

12

C Finite size scaling of energy gaps

15

References

161

Introduction

Critical phenomena emerging at classical and quantum phase transitions are of great interest
due to their experimental relevance and theoretical significance [2, 3]. Many critical phenomena are believed to be described by conformal field theories (CFTs), which are stronglyinteracting and pose challenges for studies in higher space-time dimensions beyond 2D (i.e.,
1+1D). A recent method known as fuzzy (non-commutative) sphere regularization [1] hasbeen invented to investigate 3D (i.e., 2+1D) critical phenomena governed by 3D CFTs on acylindrical geometry represented as S 2 × R. Compared to traditional lattice regularization,
the fuzzy sphere regularization offers numerous advantages in the study of 3D CFTs, primarily
due to the utilization of radial quantization in S 2 × R [4, 5] as well as the exact preservationof sphere SO(3) symmetry [6, 7], as convincingly demonstrated recently [1, 8–11].
Firstly, the fuzzy sphere enables direct access to information regarding the emergent conformal symmetry in the critical state [1, 10]. Secondly, it allows for the direct extraction ofvarious data of the CFTs, including numerous scaling dimensions of conformal primary operators [1, 10], operator product expansion coefficients [8], and four-point correlators [9].
For instance, scaling dimensions can be computed directly from excitation energies of the system, and their accuracy can be further improved using conformal perturbation [12]. Thirdly,
the fuzzy sphere scheme is applicable to a variety of 3D CFTs, including Ising [1], O(N)
Wilson-Fisher, SO(5) deconfined phase transition [10], critical gauge theories [10], and defect CFTs [11]. Lastly, the fuzzy sphere regularization exhibits an incredibly small finite-size
effect when the Hamiltonian is reasonably fine-tuned. These advantageous features of fuzzysphere regularization present an exciting opportunity to explore 3D CFTs with high efficiency,
accuracy, and comprehensiveness.
The fuzzy sphere regularization considers a microscopic quantum Hamiltonian modelingfermions (with multiple flavors) on continuous spherical space and projecting fermions into thelowest spherical Landau level [1,6,13]. In comparison with the regular lattice model, the fuzzysphere model preserves the continuous rotational symmetry exactly in the UV limit. Thanks tothe extremely small finite-size effect achieved through fine-tuning, numerical algorithms suchas exact diagonalization (ED) and density matrix renormalization group (DMRG) methods arehighly effective in studying the fuzzy sphere model of the 3D Ising CFT and SO(5) deconfinedphase transition. However, the computational cost of these two algorithms will eventuallygrow exponentially with the system size. More importantly, for cases involving a large numberof fermion flavors, the computational costs of ED and DMRG quickly surpass practical resource
and time limitations. In these cases, it would be helpful to be able to study models on the fuzzysphere using a method that scales polynomially in time, such as quantum Monte Carlo (QMC).
The aim of this paper is to utilize the 3D Ising CFT as an example to demonstrate theapplication of the QMC approach in studying 3D CFTs on the fuzzy sphere. A similar discussion
for the fuzzy torus model can be found in Ref. [13, 14]. In contrast to the fuzzy sphere Isingmodel introduced in Ref. [1], we introduce an additional flavor index to the fermions, which
results in no sign problem for the QMC simulations. As a benchmark, we provide numerical2

SciPost Phys. Core 7, 028 (2024)

results of finite-size scaling, indicating that this model also belongs to the 3D Ising universality
class. Furthermore, we introduce observables that enable the extraction of energy gaps in the
spectrum corresponding to specific symmetry quantum numbers. This allows us to investigate
the presence of conformal symmetry at criticality and extract scaling dimensions through thestate-operator correspondence. Our numerical results for energy gaps are consistent with theuniversality of the 3D CFT Ising model, albeit with a larger finite-size effect compared to thepreviously studied fuzzy sphere Ising model [1]. In summary, we believe the QMC enriches
the arsenal to study the fuzzy sphere model.
This paper is organized as follows: in Section II we introduce the model and its symmetries, and we discuss how it can be implemented in auxiliary-field QMC simulations. We
also argue for why the simulations are sign-problem-free. In Section III we discuss finite-sizescaling results and give evidence that the model is in the 3D Ising universality class, and we
discuss energy spectrum results and give evidence for emergent conformal symmetry. Section
IV contains our conclusions.

2

Model and method

2.1

Review of fuzzy sphere regularization

The fuzzy sphere regularization considers fermions moving on a sphere in the presence ofa magnetic monopole with 4πs flux sitting in the center of the sphere. In general, we canconsider multi-flavor fermions ψα with the flavor index α, described by a Hamiltonian,H = Hkin + Hint .

(1)

Here Hkin is the kinetic term of fermions, and Hint is an interaction which takes forms such asa density-density interaction,Zd 2 ⃗r1 d 2 ⃗r2 U(⃗r1 − ⃗r2 )na (⃗r1 )n b (⃗r2 ) ,(2)awhere na (⃗r ) = ψ† (⃗r )α ψ(⃗r )β Mαβand M a is a matrix defined in the fermion flavor space.
U(⃗r1 − ⃗r2 ) is a rotationally invariant interaction, and we take it to be short ranged such asδ(⃗r1 − ⃗r2 ) and ∇2 δ(⃗r1 − ⃗r2 ).
The energy levels of Hkin form quantized Landau levels, whose wave-functions are de(s)scribed by the monopole Harmonics (i.e. spin-weighted spherical Harmonics) Yn+s,m (θ , ϕ)[15], with n = 0, 1, · · · as the Landau level index and (θ , ϕ) as the spherical coordinates.
Each Landau level has an energy En = [(n + 1/2) + n(n + 1)/2s]ωc /2π, with the cyclotron
frequency ωc [6, 16]. The states of each Landau level have an angular momentum L = s + n,
hence they are (2s + 2n + 1)-fold degenerate, which can be labeled by the quantum numberof the z-component of the angular momentum L z , m = −s − n, −s − n + 1, · · · , s + n. Because
we may set the scale of Hkin as large as we like (ωc → ∞) relative to our scale of interest for
Hint , we may enforce that there are no fluctuations out of the lowest Landau level (LLL). The
physics of interest will come from the terms that make up Hint . Hence, we consider the limitthat Hkin ≫ Hint such that we can project the system into the LLL. The annihilation operator
ψ(θ , ϕ) on the LLL can be written ass

1 X (s)ψα (θ , ϕ) = p
Ȳs,m (θ , ϕ) cm,α ,N m=−s(3)where cm,α stands for the annihilation operator of Landau orbital m, and it is independent of
coordinates (θ , ϕ). N = 2s + 1 is the number of orbitals, playing the role of area of the 2D
3

SciPost Phys. Core 7, 028 (2024)p
space. The prefactor 1/ N ensures that the density operator,
1 X (s) (s) †ana (θ , ϕ) =Y Ȳ ccm ,β Mαβ,N m ,m s,m1 s,m2 m1 ,α 2
1

(4)

2

is an intensive quantity.
Under this LLL projection Eq. (3), any rotation invariant density-density interaction in theform of eq. (2) can be written as the Haldane pseudopotentials [6] in terms of second quantizedfermion operators,
X
†
a
†
b
Vm1 ,m2 ,m3 ,m4 (cm(5),α Mαβ cm4 ,β )(cm ,η Mηγ cm3 ,γ ) ,1

m1 ,m2 ,m3 ,m4

2

with
Vm1 ,m2 ,m3 ,m4 = δm1 +m2 ,m3 +m4

2sXVl (4s − 2l + 1)

l=0



s
×
m1

s
2s − lm2 −m1 − m2



s
m4


s
2s − l,m3 −m3 − m4

(6)




j1
j2
j3
where
is the Wigner 3 j-symbol. Vl are numbers whose values are specificallym1 m2 m3
depending on the form of the interaction U(⃗r1 − ⃗r1 ). For the remainder of this work, we focusgg
on U(⃗r1 −⃗r1 ) = U(Ω12 ) = N0 δ(Ω12 )+ N12 ∇2 δ(Ω12 ) with δ(Ω12 ) = δ(ϕ1 −ϕ2 )δ(cos θ1 −cos θ2 ).
The corresponding Haldane pseudo-potentials are
V0 =2s + 1s
g0 −g1 ,4s + 14s + 1

V1 =

s
g1 ,
4s − 1

Vl≥2 = 0 .

(7)

To realize the 2 + 1D Ising transition, Ref. [1] introduced a Hamiltonian that has two flavorsof fermions ψ† = (ψ†↑ , ψ†↓ ) with their interaction,Z
Hint = N 2 dΩ1 dΩ2 U(Ω12 ) n0 (θ1 , ϕ1 )n0 (θ2 , ϕ2 ) − nz (θ1 , ϕ1 )nz (θ2 , ϕ2 )Z
N dΩ n x (θ , ϕ) ,

−h

(8)

where Ω = (θ , ϕ) is a spherical coordinate and na (θ , ϕ) = ψ† (θ , ϕ)σ a ψ(θ , ϕ) is a local
density operator with σ x, y,z being Pauli matrices, σ0 = I2×2 . The first term behaves like an
Ising ferromagnetic interaction, while the second term is the transverse field.
It is straightforward to solve the second quantized Hamiltonian Eq. (5) using unbiased
numerical algorithm such as the ED and the DMRG, although their computational costs growexponentially with the system size N = 2s+1. So it is highly desirable to develop QMC method
for the simulation of a fuzzy sphere model, and it is the focus of this paper.
It is worth mentioning why the LLL projection leads to a fuzzy sphere. We can consider the
projection of the coordinates of a unit sphere, denoted as x = (sin θ cos ϕ, sin θ sin ϕ, cos θ ).
After the projection,
the coordinates are transformed into (2s + 1) × (2s + 1) matrices, where
R
(s)
(s)
(X )m1 ,m2 = dΩ x Ȳs,m (Ω)Ys,m(Ω). These matrices satisfy the following relations:
1

2

3
X

1
iεi jk X k ,[X i , X j ] =
s+1

i=1

XiXi =s
12s+1 .
s+1

(9)

The fact that the three coordinates satisfy the SO(3) algebra formally defines a fuzzy sphere[7]. It is interesting to note that in the limit as s → ∞, the fuzziness disappears and a unit
sphere is recovered.
4

SciPost Phys. Core 7, 028 (2024)

2.2

The density form of interaction

To facilitate QMC simulation, we would like to write the Hamiltonian in terms of the density
operator in the angular momentum space nal,m , defined as,na (θ , ϕ) =1X a mn Y (θ , ϕ) .
N l,m l,m l(10)Here Ylm (θ , ϕ) is the spherical harmonics, with m = −l, −l + 1, · · · , l and l ∈ Z. nal,m can beobtained using the spherical harmonic transformation,

nal,m = N=NZdΩ Ȳlm (θ , ϕ)na (θ , ϕ)vst 2l + 1 X
4π

3s+m1

(−1)

m1 =−ss
To have the term
−m1



s
−m1


l
s
†
a
t jsls−s0scm
cm1 −m,β Mαβ.
1 ,αm m1 − m

(11)


l
s
non-vanishing, we should have l ≤ 2s. One can show
m m1 − mn†l,m = (−1)m nl,−m .
P
In this context, it is convenient to decompose the potential U(θ12 ) = l 2l+1
4π Ul Pl (cos θ12 )
Pl
4π
m
m
using the Legendre polynomials, Pl (cos θ12 ) = 2l+1 m=−l Ȳl (Ω1 )Yl (Ω2 ), such that the interaction terms take the form
Z
2s
l
X
X
N 2 dΩ1 dΩ2 U(θ12 )na (θ1 , ϕ1 )n b (θ2 , ϕ2 ) =b(nal,m )† nl,m
,

Ul

l=0

(12)

m=−l

with the coefficients Ul = g0 /N − l(l + 1)g1 /N 2 .

2.3

Four component fuzzy sphere model

In comparison to Ref. [1], we consider four flavors of fermions, ψ† = (ψ†↑,+ , ψ†↑,− , ψ†↓,+ , ψ†↓,− ),i.e., we introduce an additional “layer” degree of freedom (+, −). The Pauli-matrices σ i andτi act on the spin (↑, ↓) and layer indices, respectively. Let us define the operators n0l,m andnzl,m according to Eq. (11) with M 0 = σ0 τ0 and M z = σz τ0 , respectively. The Hamiltonian
reads
2s
l 
X
X
X

† x 0
(13)
Hint =
Ul(n0l,m )† n0l,m − (nzl,m )† nzl,m + hcm
σ τ cm ,l=0

m=−l

and the interaction favors a ferromagnetic state for g0 , g1 > 0.
There are four symmetries of this model, which are
1. Ising Z2 symmetry: cm → σ x τ0 cm .
2. SO(3) sphere rotation symmetry: cm=−s,...,s form the spin-s representation of SO(3).
†3. Particle-hole symmetry: cm → iσ y τ0 cm, and i → −i.

4. Layer SU(2): generated by cm → σ0 τ x, y,z cm .

5

SciPost Phys. Core 7, 028 (2024)

4.0

N=8
N=10
N=12
N=14

3.5

3.0

0.4

0.2

2.5

0.17

0.18

0.19

0.20

0.21

0.22

0.23

0.24

1.5

4

3

2

2.0

0.3

N=8
N=10
N=12
N=14

5

(

0.5

(

U4

)

M2 / N 4 2

0.6

c.

N=8
N=10
N=12
N=14
)

0.7

b.

M2 / N 4 2

a.
0.8

0.19

0.20

0.21

h

0.22

0.23

0.24

h

1

0.3

0.2

0.1

0.0

(h hc) N 3
(

0.1

0.2

)

Figure 1: Order parameter data for V0 = 0.5564 and V1 = 0.1, which shows evidence
for a continuous phase transition consistent with that of the 3D Ising Model. (a)Binder ratio data showing a crossing that drifts somewhat in system size: N = 12 andN = 14 cross around h = 0.21, whereas N = 8 and N = 10 cross closer to h = 0.20.
(b) Magnetization data showing consistent crossing between h = 0.21 and h = 0.22
for N = 10, 12, 14. Ising ∆σ = 0.518 is assumed. (c) Magnetization data plottedalong with a universal scaling function fit. Fixing ∆σ = 0.518 and ∆ε = 1.41 (asshown in this figure) yields a good fit (χ 2 = 1.305) with hc = 0.2129(8), consistent
with 3D Ising universality. Fitting using an h = 0.21 estimate consistent with the
Binder ratio and magnetization crossing data gives ∆σ = 0.49(2) and ∆ε = 1.28(6).
The first three of these symmetries are the same as those of the two-flavor model studiedin [1]. The layer symmetry is an additional symmetry for the four flavors, which allows for
sign-problem-free QMC simulations of this model. At the Ising transition, the layer SU(2)degrees of freedom need to be gapped. We have verified this in Appendix C.
Before moving on, we remark that the four component fuzzy sphere model Eq. 13 is not a
simple product of the two-component model Eq. 8, so the phase transition point hc of the four
component model is distinct from that of the two-component model. While the interaction ischosen to be of density form in the layer degree of freedom, i.e., using τ0 , to disfavor layermagnetism, the two layers are coupled and spontaneous layer-symmetry breaking cannot ruledout in general. Nevertheless, we ensure the universality of the four component model falls in
the 3D Ising class, as shown below.

3

Results

3.1

QMC simulations

We simulate the model, eq. (13), using projector auxiliary Quantum Monte Carlo (AFQMC).
To fit this goal, we rewrite the sums of quartic terms in the following way
Ul

l 
X
m=−l

nal,m

†

nal,m = Ul

l
X



(−1)m nal,−m nal,m

m=−l
l 
X

(14)

2
Ul
=
(1 + i)nal,m + (1 − i)(−1)m nal,−m ,4 m=−lwhere in the first equality we make use of the density operator identity
nl,m † = (−1)m nl,−m .
6(15)

SciPost Phys. Core 7, 028 (2024)

The squared operators in the second line of (14) are Hermitian, and thus AFQMC as implemented in [17] is applicable. The projector we use is the half-filled solution to thewhen
P model† x 0
g0 = g1 = 0, where the Ising spins are polarized by the transverse field term h cm
σ τ cm .
Now we show, that the QMC simulation of this model is sign-problem-free.
Afp
ter the Hubbard-Stratonovich transformation, we have a prefactor of−∆τUl /4. If
g0 − g1 l(l + 1)/(2s + 1) of the expression in (13) is always positive, then we get an extra
factor of i for the n0 terms, which picks up a sign under antiunitary transformations. The
antiunitary particle-hole transformation P ,†cm→ iσ y τ0 cm ,i → −i ,nzl,m → (−1)m nzl,−m ,n0l,m → −(−1)m n0l,−m ,(16)

is a symmetry with P 2 . Combined with the SU(2) layer symmetry, it guarantees the absence ofthe sign-problem in this model [18]. Hence, the computational complexity scales polynomiallywith system size N . However, in this basis, the auxiliary fields couple to the operators nal,m oftypical rank N , compared to the usual rank of O(1) in lattice models. Therefore, the compute
time of the algorithm scales as N 4 instead of the conventional N 3 .
In this particular model, we focus on the critical point that occurs in a regime where both
g0 , g1 > 0. We have not proven in the discussion above that the absence of a sign problem
occurs when instances of g0 − g1 l(l + 1)/(2s + 1) is positive for small l but is negative for largel, yet we encounteredproblem in our simulations. One explanation may be that the no sign s l s
decays exponentially in l, causing a suppression of terms that
Wigner-3j prefactor−s 0 schange the overall prefactor signs in (13), and so the smallness of the couplings of these terms
may be important.
Due to the nonlocal nature of the operators in (13), controlling Trotter discretization errors
becomes a more demanding task, as observed in [14]. We alleviate some of these effects
by adopting a stabilized second-order Trotter decomposition developed by Blanes et al., as
discussed in [19]. The effectiveness of these alternate splitting schemes in the realm of AFQMC
was shown in [20]. Furthermore, to implement the Wigner-3j prefactors efficiently, we utilize
the software package detailed in [21].
Here we utilize QMC to compute the evolution of the order parameter and CFT dimensionless two-point correlators across the transition point, as well as extract energy gaps forthe lowest lying states using time-displaced correlation functions (see Appendix B). While theLowest Landau level basis has already been used for a QMC study in [22], the CFT-inspired
use of time-displaced correlation functions and two-point correlators on the fuzzy sphere is
new for QMC studies.
Below we will set V1 = 0.1, V0 = 0.5564, and tune h to realize a 2+1D Ising transition. We
find this choice of (V0 , V1 ) has a smaller finite size effect. By inverting the equations in (7), wehave
4s + 14s − 1
g0 =
V0 +
V1 ,
2s + 12s + 1(17)
4s − 1g1 =V1 ,sso this corresponds to the region g0 , g1 > 0.

7

SciPost Phys. Core 7, 028 (2024)

3.2

Finite-size-scaling

To look for a phase transition, we begin with the order parameter for the Ising phase transition,
which in the Landau Level basis (see Appendix A) is given by
M=

s
X

† z
cm
σ cm .

(18)

m=−s

In Figure 1(a), we have plotted the Binder cumulant, given by
U4 = 1 −M4
3 〈M 2 〉2,(19)

and we see a crossing in the vicinity of 0.20 − 0.21, that drifts to larger couplings with larger
N = 2s + 1. With this evidence of there being a quantum phase transition, we can find further
evidence that the phase transition is in the Ising universality class by assuming that ∆σ is equal
to 0.518, as is known for the universality class [23], and checking the the magnetization data,as seen in Figure 1(b). Here we see a good crossing for N = 10, 12, 14, which is consistent
with the choice of ∆σ . The crossing is at a larger h than the Binder cumulant crossing, thatis because the Binder cumulant suffers from larger finite size effects. Similar finite size effects
have also been observed in the two-layer model [1].
To see that the data is consistent with both η(= 2∆σ −1) and ν(= 1/(3−∆ε )), critical exponents in the 3D Ising universality class, we have performed a data collapse to a universal scalingp 4−2∆σfunction, assuming that M 2 / N
has a functional form of f0 + f1 x + f2 x 2 + f3 x 3 ,p 3−∆ε
where x = (h − hc ) N. In addition to the fit parameters, we pay attention to the quantity χ 2 /dof, with numerator χ 2 given byP
2n (On − Fn )
2
χ =
,
(20)
σ2n
where On is a measurement, Fn is the expected value from the fit, and σn is the variance for themeasurement. The quantity χ 2 /dof is reduced by the number of degrees of freedom (dof) for
the fitted data, which is the number of measurements minus the number of free parameters.
A good fit that captures the features of the data but does not overfit is close to 1. When we fix∆σ = 0.518 and ∆ε = 1.41 and leave the other five parameters free, we get a fit for the data
N = 10, 12, 14, with χ 2 /dof = 1.305 and estimate for the critical coupling of hc = 0.2129(8),as seen in Figure 1(c). If instead we fix h = 0.21, as suggested by the Binder ratio, andleave six parameters including the critical exponents free, a fitting (χ 2 /dof = 1.491) gives us∆σ = 0.49(2) and ∆ε = 1.28(6), values consistent with Ising universality for these relatively
small system sizes.

3.3

Dimensionless two-point correlator

To take the advantage of fuzzy sphere regularization, below we compute CFT dimensionless
two-point correlators on a sphere [9] at equal time,
Gφφ (θ ) = 〈φ(θ = ϕ = 0)φ(θ , ϕ = 0)〉1=,(2 sin (θ /2))2∆φ(21)

where φ is a CFT primary operator, and (θ , ϕ) are the spherical coordinates specifying thepositions of the two operators. We mainly focus on the lowest Z2 -odd primary σ, which can
8SciPost Phys. Core 7, 028 (2024)

a.
N=14
N=12
N=10
N=8
N=6

0.65

b.

0.4

h = 0.215
h = 0.212
h = 0.209
h = 0.206
h = 0.203

0.3

0.2

0.40
0.195

0.200

0.205

0.210

0.215

0.0

h

exact

3

0.50
0.45

N = 14
N = 12
N = 10
N=8
N=6

4

0.5

0.55

c.

5

0.6

proxy

0.60

proxy

0.7

G

0.70

0.1

0.2

1/ N

0.3

0.4

0.5

2

1

0

0

1

2

3

4

5

6

proxy

Figure 2: Data from CFT dimensionless two-point correlators. (a) Values of ∆σ
proxyfrom Gσσ using the equal-time correlators. The crossing of the ∆σvalues throughthe ∆σ = 0.518 critical exponent value occurs at larger h-values as N increases, consistent with the finite-size-scaling drift that was observed earlier. (b) Extrapolation
to infinite lattice size ∆σ from Gσσ correlation functions for different values of h and
N = 6, 8, 10, 12, 14. Linear fits suggest that, from this dataset, h = 0.212 is closest to
criticality. (c) Plotting of the angle-dependence of the Gσσ correlator with data fromN = 6, 8, 10, 12, 14 at h = 0.212.
z
be well approximated by the UV operatorpp n [8, 9], up to a non-universal normalization (sayA) and higher order corrections O(1/ N ) from operators with higher scaling dimensions.
So we can first compute the equal-time two-point correlator,f (θ ) = 〈nz (θ = ϕ = 0)nz (θ , ϕ = 0)〉
=

2s
X

¶
¬
Ȳl,m=0 (θ , 0)Yl,m=0 (0, 0) nzl,0 nzl,0 ,

(22)

l=0

and thenp Gσσ (θ ) = Af (θ ) + O 1/ N ,(23)

where A is a nonuniversal number. Because we have an explicit expression for f , we know
the exact values for f (θ = π) and f ′′ (θ )|θ =π , where the derivatives are taken in θ . Then by
assuming that Gσσ has the critical scaling form of (21), we can solve the following system,proxy

Af (π) = 1/(2 sin(π/2))2∆σ ,proxy ∂2 
2∆σAf ′′ (θ ) θ =π =
1/(2
sin(θ
/2))
.
θ =π
∂ θ2
proxy

(24)

∆σ
is a number that will extrapolate to the universal ∆σ at the critical point as N → ∞.
Second derivatives are used for the second equation in (24) because the first derivatives in θ
are zero for both f and critical Gφφ at θ = π.
proxyFigure 2(a) shows the extracted ∆σvalues for different values of h in the vicinity ofproxy
the hc determined by finite-size-scaling. Here we see that the ∆σindeed crosses throughthe 3D Ising ∆σ = 0.518 value in this region, and furthermore we see that the h at which thiscrossing occurs increases with system size N , which is consistent with the drift that we saw inthe finite-size-scaling. Moreover, the drift appears to be slowing with increasing N , anotherconsistency with finite-size scaling.
We can see more consistencies with finite-size-scaling from the results in Figure 2(b), which
proxylinearly extrapolate the values of ∆σas N → ∞ for different values of h. Here we see that9

SciPost Phys. Core 7, 028 (2024)

for the infinite N extrapolation, the h = 0.212 data is closest to the critical ∆σ , which isconsistent with the earlier universal scaling fit of hc = 0.2129(8). We use this h = 0.212 data
to show the calculation of Gσσ (θ ) as a function of θ for QMC data at N = 6, 8, 10, 12, 14 andhow it approaches the exact expression as N increases. The QMC data-derived expressions forGσσ (θ ) are given in Figure 2(c) and are rescaled so that they are equal to the exact expressionof Gσσ (θ = π).

3.4

Energy gaps and state-operator correspondence

Next we turn to the state-operator correspondence [4, 5] on the sphere, namely, the scaling
dimensions ∆n are related to energy gaps byδEn = En − E0 =

v
∆n ,
R

(25)

where R is the radius of the sphere and v is the model-dependent velocity of light.
While we are unable to get the full low lying energy spectrum directly using QMC, we
are able to obtain energy gaps for the lowest lying states in each symmetry quantum number
sector by using time-displaced correlation functions (see Appendix B). For an operator OS with
the quantum number S, we have:X(26)〈OS (τ) OS (0)〉 =an2 e−τ( ES,n −E0 ) ,nwhere E0 is the ground state energy, ES,n represents the energies of eigenstates |ψn 〉 in thequantum number sector S, an is an operator OS and state |ψn 〉 dependent non-universal factor.
At long time τ ≫ 1, the lowest energy will dominate and can be extracted by fitting the
exponential decay.
In the data that follows, we will use density operators nil,m to measure energy gaps indifferent quantum number sectors. Specifically,1) nzl,m can measure gaps in the Z2 -odd, parity-even, and angular momentum (i.e. Lorentzspin) l sector,x2) nl,mcan measure gaps in the Z2 -even, parity-even, and angular momentum l sector,

3) n0l,m can measure gaps in the Z2 -even, parity-odd, and angular momentum l sector.
Figure 3 shows QMC data at the critical point h = 0.212. The energy gaps are scaled
such that the gap measured from n0l=2,m=0 is rescaled to 4, the scaling dimension of the lowest
parity-odd descendent of the energy-momentum tensor, ∆ϵνρη ∂ρ Tµν . In doing so, we find thegaps measured from other operators to be consistent with the scaling dimensions of primaryxand descendant operators of the 3D Ising CFT. The density operator (nl=2,m=0) we measured
does not seem to have an overlap with the state of stress tensor (with ∆ Tµν = 3). Instead it
gives the level-2 descendant of ε primary, i.e., ∂µ ∂ν ε.

4Conclusions

We have introduced a model that is amenable to using sign-problem free quantum Monte
Carlo to simulate the (2+1)−D transverse Ising model on a fuzzy sphere. Through finite-size
scaling we have found data consistent with the model’s phase transition being in the 3D Ising
universality class, and we also have shown that we can recover the same critical exponents
from the model’s energy spectra, which is evidence of emergent conformal symmetry.
10

SciPost Phys. Core 7, 028 (2024)

CFT, primary and descendents
CFT, primary and descendents
CFT, T primary and a descendent
QMC, 2 odd and parity even
QMC, 2 even and parity even
QMC, 2 even and parity odd

4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0

1

l

2

Figure 3: Comparison between 3D Ising CFT data and rescaled energy gaps for N = 8
and h = 0.212 measured by density operators. The energy gaps are rescaled by thefactor that makes the lowest Z2 -even, parity-odd gap at l = 2 equal to the lowest
parity-odd descendent of the energy-momentum tensor, 4.0.
While these calculations are not competitive with ED and DMRG for small lattices, thiswork opens the door to larger scale calculations for models where there are too many degrees
of freedom to make ED/DMRG calculations infeasible, or for when large sizes are desiredfor more accurate determination of critical exponents. Here, we were able to introduce anadditional layer degree of freedom to avoid the sign-problem, and chose the interaction to
disfavor spontaneous layer-symmetry breaking energetically. However, such instabilities cannot be ruled out in general. However, there are many interesting critical phenomena which
naturally involve multiple flavors; one interesting target is the critical gauge theories proposed
in Ref. [10].

Acknowledgments
We thank Chao Han and Zheng Zhou for stimulating discussions.
Funding information Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Industry Canada andby the Province of Ontario through the Ministry of Colleges and Universities. WZ is supported by National Natural Science Foundation of China (No. 92165102) and foundation of
the Westlake University. JH was supported by the European Research Council (ERC) under
grant HQMAT (Grant Agreement No. 817799), the Israel-US Binational Science Foundation(BSF), and by a Research grant from Irving and Cherna Moskowitz. FG acknowledges financial support through the German Research Foundation, project-id 258499086 - SFB 1170
‘ToCoTronics’. This work used Bridges 2 at the Pittsburgh Supercomputing Center throughallocation PHY170036 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants#2138259, #2138286, #2138307, #2137603, and #2138296.

11

SciPost Phys. Core 7, 028 (2024)

A

Order parameter

The Ising order parameter is nz , and for QMC simulations we measure the two-point correlation function,X′′′〈nz (θ , ϕ)nz (θ ′ , ϕ ′ )〉 =〈nzl,m nzl′ ,m′ 〉Ylm (θ , ϕ)Ylm′ (θ , ϕ )l,m,l ′ ,m′X=(−1)m 〈nzl,0 nzl,0 〉Ylm (θ , ϕ)Yl−m (θ ′ , ϕ ′ ) .

(A.1)

l,m

The last equation comes from the conservation of angular momentum,


l
l′ 0
z
z
〈nl,m nl ′ ,m′ 〉 =Om m′ 0 l




l l 0
l
l
0
z
z
〈nl,0 nl,0 〉/= δl,l ′ δm,−m′
0 0 0m −m 0= (−1)m δl,l ′ δm,−m′ 〈nzl,0 nzl,0 〉 .

(A.2)

Therefore, we need to evaluate 〈nzl,0 nzl,0 〉 for each l. To do the finite-size-scaling, we calR
culate the order parameter 〈M 2 〉, with M = dΩ nz (θ , ϕ),2〈M 〉 ==

ZdΩdΩ′ 〈nz (θ , ϕ)nz (θ ′ , ϕ ′ )〉Z
X(−1)m 〈nzl,0 nzl,0 〉dΩdΩ′ Ylm (θ , ϕ)Yl−m (θ ′ , ϕ ′ )l,m= 4π〈nz0,0 nz0,0 〉 .

(A.3)

R
p
The last equation comes from dΩYlm (θ , ϕ) = 4πδl,0 δm,0 . Using Wigner-3j identities, wefind that,
v



s
t 1 X
s
0 s
s 0 s † zz
3s+m1
n0,0 = (2s + 1)(−1)
c σ c m1
−m1 0 m1
−s 0 s m1
4π m =−s
1
v
s
t 1 X
=
c † σ z c m1 .
(A.4)
4π m =−s m1
1

Therefore, the order parameter M 2 is〈M 2 〉 =s
X
m1 ,m2 =−s

B

†
†
〈(cm
σz cm1 )(cm
σz cm2 )〉 .
1

2

(A.5)

Extracting energy gaps

In projector QMC, we are able to get the energy gap between the first excited state in symmetry
sector S and the ground state in the following way. Assuming a trial wavefunction, |ψ0 〉, an
operator that creates overlap between
the states in symmetry sector S and the ground state,P
OS , and a complete set of states n |n〉 〈n|, where |n〉 is an eigenstate with energy eigenvalue

12

SciPost Phys. Core 7, 028 (2024)

a.

100

fit
fit-comp
data
obs

10 2
10 3

10 1
0

5

10

15

t

10 4

20

c.

10 3

fit
fit-comp
data

10 1

0

5

10

15

20

t

d.

25

30

fit
fit-comp
data

10 4

obs

10 2

obs

fit
fit-comp
data

10 1

obs

100

b.

10 3

10 5

10 4

10 6
0.0

2.5

5.0

7.5

10.0

t

12.5

15.0

17.5

0

20.0

1

e.

3

t

4

5

6

7

fit
fit-comp
data

10 2

obs

2

10 3
10 4
10 5
0

2

4

6

t

8

10

12

Figure 4: Fits to one exponential using semilog plots. The vertical lines show thelocations of the endpoint and startpoints for the data. Plots (a), (b), and (c) show
the data from nz0,0 , nz1,0 , and nz2,0 , respectively. Plot (d) shows the data from n02,0 andxplot (e) shows the data from n2,0.
En , we have that
〈OS (τ) OS (0)〉 ==〈ψ0 | e−(β−τ)H OS e−τH OS |ψ0 〉〈ψ0 | ψ0 〉X 〈ψ | e−(β−τ)H O |n〉 〈n| e−τH O |ψ 〉
0

S

n

=

S

X 〈ψ | O |n〉 〈n| O |ψ 〉0n0

〈ψ0 | ψ0 〉S

〈ψ0 | ψ0 〉

S

0

(B.1)

e−β E0 e−τ(En −E0 ) .

This term has contributions from all eigenstates that have the symmetry S. The higher energy
states will have gaps that will be suppressed relative to the smallest energy gap, and so we canapproximate
01〈OS (τ) OS (0)〉 = C1 e−τ( ES −E0 ) + C2 e−τ( ES −E0 ) ,(B.2)
where ES0 and ES1 are the lowest energy and second lowest energy corresponding to states insymmetry sector S, respectively.
In practice, we found that in order to extract ES0 , sometimes a fit to the two exponentials
with prefactors C1 and C2 in (5) is necessary, but sometimes a fit to a single exponential (whichassumes C2 = 0) is more appropriate. The procedure for fitting to one versus two exponentials
involves the following steps:
13

SciPost Phys. Core 7, 028 (2024)

fit
fit-unrestrictedfit-part
fit-part
fit-part_ur
fit-part_ur
data

obs

10 1

10 2

b.

fit
fit-unrestrictedfit-part
fit-part
fit-part_ur
fit-part_ur
data

10 1
10 2

obs

a.

10 3
10 4

10 3
0

2

4

6

t

8

10

12

14

0

2

4

6

8

t

10

12

14

16

Figure 5: Fits to two exponentials using semilog plots. The vertical solid lines show
the locations of the endpoint and midpoint for a restricted fit–where the C1 exponential and the C2 exponential are fitted separately but iteratively using informationfrom previous results to arrive at a final answer. The dashed lines give an interval of
midpoints that were tested in order to estimate the error bar. The solid diagonal lines
give the two-exponential restricted and unrestricted fits (“fit” and “fit-unrestricted”).
The dashed diagonal lines are fits for each of the two exponentials, both for restrictedx(“fit-part”) and unrestricted (“fit-part_ur”) fits. Plot (a) shows the data from n0,0
, andx
plot (b) shows the data from n1,0 .
1. Find the τ interval where the data is distinguishable from zero according to error bars.
The largest time in this interval is the initial “endpoint” guess.
2. Initially guess that the “midpoint” in time–where one exponential versus the other exponential dominates–is 30% of the full time interval.
3. Test a fit to a single exponential–if the initial data point is smaller than the t = 0 value
for the fitted function, gradually adjust the midpoint and endpoint guesses down until
this is not the case.
4. If the initial data point is still within errors of the t = 0 value for the single-exponential
fitted function, fit the data to a single exponential. If not move on to a two-exponential
fit and then revise the midpoint guess such that the value of the larger exponential in
the fit is negligible compared to that of the smaller exponential at the midpoint.
x
For the operators nzl,0 , n02,0 , and n2,0
, the fit to a single exponential ends up being more apx
propriate. Figure 4 shows the fits to a single exponential (in the cases of n22,0 and n2,0
, wechose a single exponential because there was very little small τ data to fit to a higher expox
x
nential). However, the n0,0
, n1,0observables require two exponentials. Figure 5 shows thetwo exponential fits for these operators at coupling h = 0.212 and s = 3.5.
For the two-exponential fits, we first use the midpoint data to iteratively fit one exponential
at a time: midpoint to endpoint is the fit for the lower energy and then the startpoint to the
midpoint is the fit for the higher energy. We alternate fitting one exponential versus the otherwhile fixing the parameters of the nonfitted exponential according to the previous fit. This
gives us the “restricted” fit listed as “fit” in Figure 5. We then use these fitted energy values as
initial guesses for an “unrestricted” fit that fits both exponentials at once. This gives the “fitunrestricted” in Figure 5. The value of the energy estimate is the mean of these restricted andunrestricted fit energies. Finally, we obtain error bars by performing fits to a single exponential

14

SciPost Phys. Core 7, 028 (2024)

Layer SU(2) Singlet

0.9

2even, PH even, l=2
2even, PH even, l=1
2 even, PH even, l=02odd, PH even, l=22odd, PH even, l=12odd, PH even, l=0

0.8
0.7

ES E0

0.6
0.5
0.4
0.3
0.2
0.1
0.0
0.00

0.05

0.10

0.15

0.20

0.25

1/ N

0.30

0.35

0.40

Figure 6: Data showing
the energy gaps obtained from the SU(2) singlet operators aspa function of 1/ N for N = 8, 10, 14. This data is in the vicinity of the critical point
at V1 = 0.1, V0 = 0.5564, h = 0.2. The gaps decrease with system size, as expected.
Extrapolations are shown for the Z2 -odd sector, but not the even sector since theerror bars are so large.
for the smaller energy from a midpoint to the endpoint, where we calculate the midpoint asτ=−1(ES1 − ES0 ) ln



εC1
C2

,

(B.3)

1
where ε is a small number representing the time when the value of C2 e−τ( ES −E0 ) is ε times0
C1 e−τ( ES −E0 ) . We take a range of ε ∈ {0.01, 0.1} to fit ES0 and use this range of ES0 values toestimate the error for the energy. The boundaries of this range of midpoints are given by thedashed vertical lines in Figure 5.

C

Finite size scaling of energy gaps

Because the QMC model studied has an additional SU(2)-layer symmetry, one check to makeis whether the degrees of freedom in the layer SU(2) non-singlet representations are gappedat the phase transition. All the operators in the Hamiltonian are of the form σ i τ0 , which are
the singlets of the layer SU(2)psymmetry. Figure 6 shows the energy gaps obtained from these
operators as a function of 1/ N at h = 0.2, which is in the vicinity of the critical point. Allgaps are decreasing with system size and the gaps for the Z2 -odd sector seem to be trending
linearly towards the origin, as required by the state-operator correspondence. The Z2 -evensector gaps are also decreasing with increasing system size, but they have larger error bars
due to interference with other higher energy descendants in their spectrum, and the details of
their fits are given in Appendix B of the Supplementary Material. On the other hand, the layer
SU(2) non-singlet, e.g., layer triplet gaps should be finite in the thermodynamic
limit. Figurep
i x, y,z
7 shows layer triplet energy gaps measured by σ τ, as a function of 1/ N , and from here
we see that these excitations appear to be gapped.
15

SciPost Phys. Core 7, 028 (2024)

Layer SU(2) Triplet

0.9

2 even, PH even, l=2
2 even, PH even, l=1
2 even, PH even, l=02odd, PH even, l=22odd, PH even, l=12odd, PH even, l=0

0.8
0.7

ES E0

0.6
0.5
0.4
0.3
0.2
0.1
0.0
0.00

0.05

0.10

0.15

0.20

0.25

1/ N

0.30

0.35

0.40

Figure 7: Data showing
the energy gaps obtained from operators of the form σ i τzp
as a function of 1/ N for N = 8, 10, 14 (operators of the form σ i τ x, y would givethe same states, making this the SU(2) triplet symmetry class). This data is in thevicinity of the critical point at V1 = 0.1, V0 = 0.5564, h = 0.2. These energies appear
to be gapped out.

References[1] W. Zhu, C. Han, E. Huffman, J. S. Hofmann and Y.-C. He, Uncovering conformal symmetry
in the 3D Ising transition: State-operator correspondence from a quantum fuzzy sphereregularization, Phys. Rev. X 13, 021009 (2023), doi:10.1103/PhysRevX.13.021009.
[2] S. Sachdev, Quantum phase transitions, Cambridge University Press, Cambridge, UK,
ISBN 9780511973765 (2011), doi:10.1017/CBO9780511973765.
[3] J. Cardy, Scaling and renormalization in statistical physics, Cambridge University Press,
Cambridge, UK, ISBN 9781316036440 (1996), doi:10.1017/CBO9781316036440.
[4] J. L. Cardy, Conformal invariance and universality in finite-size scaling, J. Phys. A: Math.
Gen. 17, L385 (1984), doi:10.1088/0305-4470/17/7/003.
[5] J. L. Cardy, Universal amplitudes in finite-size scaling: Generalisation to arbitrary dimensionality, J. Phys. A: Math. Gen. 18, L757 (1985), doi:10.1088/0305-4470/18/13/005.
[6] F. D. M. Haldane, Fractional quantization of the Hall effect: A hierarchy of incompressible
quantum fluid states, Phys. Rev. Lett. 51, 605 (1983), doi:10.1103/PhysRevLett.51.605.
[7] J. Madore, The fuzzy sphere, Class. Quantum Gravity 9, 69 (1992), doi:10.1088/02649381/9/1/008.
[8] 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),
doi:10.1103/PhysRevLett.131.031601.

16

SciPost Phys. Core 7, 028 (2024)[9] C. Han, L. Hu, W. Zhu and Y.-C. He, Conformal four-point correlators of the 3D
Ising transition via the quantum fuzzy sphere, Phys. Rev. B 108, 235123 (2023),
doi:10.1103/PhysRevB.108.235123.
[10] Z. Zhou, L. Hu, W. Zhu and Y.-C. He, The SO(5) deconfined phase transition under the
fuzzy sphere microscope: Approximate conformal symmetry, pseudo-criticality, and operator
spectrum, (arXiv preprint) doi:10.48550/arXiv.2306.16435.
[11] L. Hu, Y.-C. He and W. Zhu, Solving conformal defects in 3D conformal field theory using
fuzzy sphere regularization, (arXiv preprint) doi:10.48550/arXiv.2308.01903.
[12] B.-X. Lao and S. Rychkov, 3D Ising CFT and exact diagonalization on icosahedron, SciPost
Phys. 15, 243 (2023), doi:10.21468/SciPostPhys.15.6.243.
[13] M. Ippoliti, R. S. K. Mong, F. F. Assaad and M. P. Zaletel, Half-filled Landau levels: A
continuum and sign-free regularization for three-dimensional quantum critical points, Phys.
Rev. B 98, 235108 (2018), doi:10.1103/PhysRevB.98.235108.
[14] Z. Wang, M. P. Zaletel, R. S. K. Mong and F. F. Assaad, Phases of the (2 + 1) dimensional
SO(5) nonlinear sigma model with topological term, Phys. Rev. Lett. 126, 045701 (2021),
doi:10.1103/physrevlett.126.045701.
[15] T. T. Wu and C. N. Yang, Dirac monopole without strings: Monopole harmonics, Nucl. Phys.
B 107, 365 (1976), doi:10.1016/0550-3213(76)90143-7.
[16] M. Greiter, Landau level quantization on the sphere, Phys. Rev. B 83, 115129 (2011),
doi:10.1103/PhysRevB.83.115129.
[17] F. Assaad et al., The ALF (algorithms for lattice fermions) project release 2.0. Documentation for the auxiliary-field quantum Monte Carlo code, SciPost Phys. Codebases 1 (2022),
doi:10.21468/SciPostPhysCodeb.1.
[18] C. Wu and S.-C. Zhang, Sufficient condition for absence of the sign problem inthe fermionic quantum Monte Carlo algorithm, Phys. Rev. B 71, 155115 (2005),
doi:10.1103/PhysRevB.71.155115.
[19] S. Blanes, F. Casas and J. M. Sanz-Serna, Numerical integrators for the hybrid Monte Carlo
method, SIAM J. Sci. Comput. 36, A1556 (2014), doi:10.1137/130932740.
[20] F. Goth, Higher order auxiliary field quantum Monte Carlo methods, J. Phys.: Conf. Ser.
2207, 012029 (2022), doi:10.1088/1742-6596/2207/1/012029.
[21] H. T. Johansson and C. Forssén, Fast and accurate evaluation of Wigner 3j, 6j, and 9j
symbols using prime factorization and multiword integer arithmetic, SIAM J. Sci. Comput.
38, A376 (2016), doi:10.1137/15m1021908.
[22] B.-B. Chen, X. Zhang, Y. Wang, K. Sun and Z. Y. Meng, Phases of (2 + 1)D SO(5) non-linear
sigma model with a topological term on a sphere: Multicritical point and disorder phase,(arXiv preprint) doi:10.48550/arXiv.2307.05307.
[23] M. Hasenbusch, Finite size scaling study of lattice models in the three-dimensional Ising
universality class, Phys. Rev. B 82, 174433 (2010), doi:10.1103/PhysRevB.82.174433.

17