A structural test for the conformal invariance of the critical 3d Ising model

How can a renormalization group fixed point be scale invariant without being conformal? Polchinski (1988) showed that this may happen if the theory contains a virial current -- a non-conserved vector operator of dimension exactly $(d-1)$, whose divergence expresses the trace of the stress tensor. We point out that this scenario can be probed via lattice Monte Carlo simulations, using the critical 3d Ising model as an example. Our results put a lower bound $\Delta_V>5.0$ on the scaling dimension of the lowest virial current candidate $V$, well above 2 expected for the true virial current. This implies that the critical 3d Ising model has no virial current, providing a structural explanation for the conformal invariance of the model.


Introduction
It is believed that the critical point of the 3d ferromagnetic Ising model is conformally invariant.
One strong piece of evidence is the excellent agreement between the critical exponents extracted from experiments and Monte Carlo simulations and from the conformal bootstrap [1,2]. Conformal invariance has been also checked directly on the lattice, by verifying functional constraints that it imposes on the shape of some correlation functions [3]. 1 In this paper we will provide another lattice test of this property, which is qualitatively different and in a sense more robust.
Any field theory coming from a local action, and in particular the 3d Ising model close to or 1 We would also like to point out a related lattice study of conformal invariance in 3d percolation [4].
at the critical temperature, has a local stress tensor operator T µν which is conserved: ∂ µ T µν = 0.
The structural property of conformally invariant local theories is that this local stress tensor operator is traceless: Our new test will probe this structural property, unlike previous lattice studies which tested its consequences.
The key question is: could the critical 3d Ising model be scale invariant (as befits any critical theory, being a fixed point of a renormalization group flow), but not fully conformally invariant? As was lucidly explained by Polchinski [5], 2 a theory will be scale invariant without being conformal if T µν is not traceless but its trace is a total divergence: where W µ is a vector operator, called the virial current, which is (a) not conserved and (b) not itself a total derivative. 3 Precisely this mechanism is responsible for scale without conformal invariance of the theory of elasticity, perhaps the simplest physically relevant example of this phenomenon [8]. 4 It's then natural to inquire if Eq. (1.2) can hold in the critical 3d Ising model, and we will show that it cannot. Our argument is based on the following simple observation: any operator W µ which is a candidate to appear in the r.h.s. of (1.2) must have two additional properties.
First of all, it should, just as T µν itself, be invariant under the internal symmetry of the model, Z 2 in the case of Ising. In addition, since T µν has canonical scaling dimension d, operator W µ should have dimension d − 1 = 2.
For the subsequent discussion, let us define V µ as the lowest Z 2 -even vector operator V µ , which is not a total derivative. If we manage to show that ∆ V > 2, this will imply that the model has no virial current candidates of appropriate dimension, and thus must be conformal.
Extending the discussion from d = 3 to the whole family of Z 2 -invariant Wilson-Fisher fixed points for 2 d 4, the dimension of V can be determined exactly in d = 2 and d = 4 (see appendix A). Namely, we have: 5 2 See also [6] for a review. Concerning the 3d Ising model, see especially section 4.2 of [7]. 3 If Wµ is a total derivative, the stress tensor can be "improved" to be traceless, so that Eq. (1.1) is satisfied for the improved Tµν . 4 It should be noted that this mechanism may be realized with a quirk in gauge theories. Namely it may happen that Eq. (1.2) holds but that the virial current is not a gauge invariant operator (and so is not a physical local operator). For example, this is how the 3d Maxwell theory avoids conformal invariance [9].
( 1.4) It also follows from the -expansion that the dimension of V in 4− dimensions will be 11±O( ). 6 Based on these results, one can expect that the dimension of V µ in critical 3d Ising model should be significantly larger than 2.
In this paper we will show, using lattice Monte Carlo simulations, that this expectation is correct. Namely, our analysis will imply a numerical lower bound on ∆ V : ∆ V > 5.0 (3d Ising) . (1.5) In particular, this proves that ∆ V > 2, and shows that the 3d Ising model has no candidates for W µ . This rules out the scale without conformal invariance scenario based on (1.2), and thus provides a new test of conformal invariance.
The paper is structured as follows. In section 2, we set up the lattice Monte Carlo simulation to measure a one-point function in a cubic lattice with peculiar boundary conditions (motivated in appendices C and D). Section 3 contains our numerical results that lead to (1.5). We conclude with a short discussion of the implications of our result. In appendix A, we compute ∆ V in the 2d Ising model and in the theory of a free massless scalar in d = 4. In appendix B, we summarize the general procedure for matching lattice operators with local operators of the critical field theory.
This is well known among the practitioners but we do not know any good pedagogical summary in the literature.

Lattice setup
We simulate the nearest-neighbor ferromagnetic 3d Ising model on the cubic lattice at the critical temperature. The Hamiltonian is We use the known critical temperature β = β c ≈ 0.2216546 [14,15].

Boundary conditions
Our lattice has spatial extent L × L × L sites. We set lattice spacing a = 1. Due to the difficulties of measuring a rather high scaling dimension ∆ V , we will only be able to go up to volumes L = 16.
We impose periodic boundary conditions in directions x 1 , x 2 , while at x 3 = 0 and x 3 = L − 1 we impose a mixture of fixed and free boundary conditions. Namely, for x 3 = 0 we impose the fixed s = +1 boundary condition for points with L/4 x 1 < 3L/4, while at x 3 = L − 1 we do the same 6 The coefficient of the O( ) correction term could be computed, but we don't need it.
for points with L/2 x 1 < L. The rest of the boundaries at x 3 = 0 and x 3 = L − 1 has free boundary conditions (see Fig. 1). The reasons for such a bizarre choice of boundary conditions will be explained shortly.
periodic periodic free s = +1 s = +1 free free x 1 x 2

Lattice operator
We will work with the lattice operator where x is a lattice point and is the symmetric lattice derivative in the ν direction.
Actually the precise form of the operator is unimportant, the only important thing is that O lat µ is not a total lattice derivative.

Matching of the lattice operator with critical point operators
Close to the critical point, the lattice operator O lat µ can be expanded into a basis of local operators of the critical theory with well-defined scaling dimensions (see appendix B for a review): where O i is the critical theory operator which has a scaling dimension ∆ i , and c i are some latticedependent constants. Barring accidental cancellations, any lattice measurement related to O lat µ will be dominated by operators of lowest scaling dimensions appearing in the r.h.s. of (2.2). This is because the contribution of an operator of dimension ∆ i will be suppressed by 1/R ∆ i where R is a large distance scale (clearly we have to go to large distances to explore the critical point).
Notice that operators in the r.h.s. will have to be vectors, but they don't have to be primaries.
So, the total derivative terms involving derivatives of various Z 2 -even scalar operators which exist in the 3d Ising model (see Table 2 in [2]) are expected to appear in the r.h.s. of (2.2). The lowest of these are ∂ µ ε and ∂ µ ε , where ε, ε are the lowest-dimension Z 2 -even scalars, of dimension ∆ ε ≈ 1.41, ∆ ε ≈ 3.83. These derivative operators (especially ∂ µ ε) have rather low dimension.
Below we will introduce a trick which will allow us to project them out and focus on more interesting terms.
Crucially for us, since O lat µ is not a total derivative, the operator V µ we are interested in will appear in this expansion: is an unknown, non-universal, lattice quantity, and we will assume C = 0 since there is no reason to expect otherwise. The . . . include various terms which we are not interested in, and we should make sure that those terms do not mask the contribution of V µ .
Some of these terms involve operators of higher scaling dimension than V . The presence of those terms is harmless since their effect will be subleading in the large volume limit. More annoying are the total derivative terms involving derivatives of various Z 2 -even scalar operators which exist in the 3d Ising model (see Table 2 in [2]). Some of these have a rather low dimension and would mask V µ unless special care is taken. For example, we expect ∂ µ ε to appear in the r.h.s. of (2.2), where ε is the lowest-dimension Z 2 -even scalar, of dimension ∆ ε ≈ 1.41.
Another class of total derivative operators which we expect to appear are ∂ ν T µν , divergences of non-conserved spin-2 Z 2 -even operators. Assuming conformal invariance, the lowest such operator has dimension ∆ T ≈ 5.51 [2]. Divergences of higher spin operators are also expected in principle but will not play a role because of their even higher dimension.
In our study we will be able to filter out the contributions of derivatives of scalars (like ∂ µ ε) through the following trick, rendered possible by the periodic boundary conditions. We consider the average value of the x 1 -component of V lat µ integrated along a periodic circle in this direction: Integration kills off the derivatives taken in the direction of integration. As a result this integrated observable in the continuum limit does not couple to derivatives of scalars like ∂ µ ε. On the other hand divergences of spin-2 operators survive this projection, and their integral will contribute to I along with the integral of V µ . 7 We will measure the one-point (1pt) function of I. In infinite volume vector operators would have zero 1pt functions, but in finite volume with appropriate boundary conditions they can be nonzero. In our case we will have 4) with no dependence on x 2 due to the translation invariance in that direction. The scaling of this observable with L will be determined by the smaller of the two dimensions ∆ V , ∆ ∂T = ∆ T + 1: In this work we will only measure ∆ I , but we will not be able to determine which of the two operators V or ∂T dominates the scaling.
Another way to determine ∆ I would be to impose periodic boundary conditions also in the x 3 direction and to study finite size scaling for the 2pt function of I at separation L/2. This observable would scale as 1/L 2∆ I . We tried this strategy and found the signal completely swamped by noise, due to large ∆ I . Using the 1pt function improves the signal-to-noise ratio by a factor L ∆ I and will allow us to perform the measurement.
The . . . terms in (2.4) decay with a higher power of L. They originate from the higherdimension operators contributing to V lat µ as well as from corrections to scaling arising from the fact that in finite volume the theory is not exactly at the critical point but is still flowing to it in the renormalization group sense. Because of limited statistics, we will unfortunately be forced to simply neglect both of these corrections in our analysis.
This function will be measured in our simulation. To have nonzero f (t), the boundary conditions at x 3 = 0, L − 1 should break the flip symmetry in the x 1 direction: under which I changes sign. This is the case for our boundary conditions in Fig. 1. On the other hand, our boundary condition preserves the above x 1 flip accompanied by the x 3 flip: and a periodic shift of the x 1 direction by L/4. As a consequence, our function f (t) will be odd with respect to t = 1/2, and in particular f (1/2) = 0. 7 To kill all possible total derivatives, one could consider periodic conditions in all directions and to integrate over the whole volume. We do not currently have a concrete proposal implementing this idea. The main difficulty is that the one-point function of a vector operator vanishes on the 3-dimensional torus with periodic boundary conditions.
We have experimented with several other flip-breaking boundary conditions, and settled for the one in Fig. 1 because it gives rise to a particularly sizable f (t), thus further improving signalto-noise. See appendix C for a list of other possible boundary conditions, and appendix D for a heuristic procedure to quickly evaluate which boundary condition is expected to work best.
While it is not directly related to our computation, we would like to mention here one other instance where boundary conditions were used in lattice field theory to make a 1pt function of a tensor operator nonzero. Namely, in 4d lattice gauge theory, the 1pt function of the off-diagonal stress tensor component T 0x was measured imposing the "shifted" boundary conditions, when the fields are made periodic in the spatial directions, and periodic up to a coordinate shift in the Euclidean time direction [16]. This boundary condition is a particular case of the gluing boundary condition discussed in appendix C.

Choice of Monte Carlo algorithm
We perform Monte Carlo simulations using the single-spin-flip Metropolis algorithm. The choice of Monte Carlo algorithms plays a crucial role in the efficiency of the simulations. It is well known that the Wolff algorithm [17] is more efficient than the Metropolis algorithm at the critical temperature due to the scaling of the computational effort with the system size. However, even though the smaller critical slowdown exponent favors the Wolff algorithm for large systems, for small ones and for some statistical observables, the Metropolis algorithm may be more efficient.
This is what happened in our case.
To be more concrete, the standard measure of the simulation efficiency is based on the product of the algorithm execution time (τ CP U ) and the integrated autocorrelation time (τ c ). One reason to prefer the Metropolis algorithm is that in our case it led to very small integrated autocorrelation time of the vector operator sampling (this time scale depends on the statistical observable we are trying to measure).
Another important factor for this choice was the role of the boundary conditions. The use of

Results
We eraging over x 2 ). Since our lattice operator (2.1) has range 3, we only did the measurement for The total number of such decorrelated spin configurations that we generated was 2.4 × 10 12 (resp. 3.5 × 10 13 ) for L = 12 (resp. L = 16). A much smaller number sufficed for L = 8. For N = L 3 /4 spin flips between the two measurements, the integrated autocorrelation time between the subsequent measurements of Obs(x 3 ) was close to 1 for every x 3 .
Our simulations were parallelized on a cluster and took a total of about 300 CPU-years.
The numerical results of these measurements are given in table 1, and are shown in plots below as a function of t = x 3 /(L − 1). 8 In these plots we show the data multiplied by (L/12) ∆ for various values of ∆. According to (2.4), the curves for different L are supposed to collapse if ∆ = ∆ I . At least this is supposed to happen for sufficiently large L, when contributions from the subleading terms . . . in (2.4) become unimportant.
In Fig. 2 we take ∆ = 2, the value needed for a virial current candidate. Clearly the curves show no collapse, ruling out the existence of the virial current.
A side remark: as mentioned in the previous section, the function f (t) should be odd with respect to t = 1/2 for our choice of the boundary conditions. This antisymmetry is indeed satisfied within error bars, as can be seen in the figures. 9 To assign an error to our determination of ∆ I , we propose the following heuristic procedure.
We vary ∆ around 6 and see when the curves clearly deviate from the collapsing behavior in the interval 0.2 t 0.8, judging by the eye. One way to quickly perform this analysis is to use the Manipulate function of Mathematica. This way we arrive at our confidence interval:

Fig. 3:
In this plot ∆ = 6, which is our central value for ∆ V .
See Fig. 4 for what the collapse plots look like at the extreme ends of the confidence interval. 10 While the "judging by the eye" procedure may seem subjective and ad hoc, we don't believe a much better statistical procedure can be advocated given our limited amount of data.
We have cross-checked our determination of ∆ I by focussing on the three points we get the same answer ∆ I = 6 ± 1. 10 If we omit the L = 8 datapoints from our analysis (e.g. if one is worried that these points are still significantly affected by the subleading . . . corrections in (2.4)), then we get ∆I = 5.5 ± 1.5 using the same procedure. We quote this number only for comparison, as we do not feel that completely discarding the L = 8 points is justified.

Discussion and conclusions
One goal of this paper was to emphasize that there is a simple and robust way to check the conformal invariance of any critical lattice model, which requires the measurement of the lowest non-derivative vector operator V which is a singlet under all global symmetries. This operator can play the role of the virial current, and potentially cause scale without conformal invariance, but only if its dimension is exactly d − 1.
In this paper we considered this strategy in the critical 3d Ising model. Since the dimension of V appears to be large, to carry out our measurement we had to introduce several tricks increasing the efficiency of Monte Carlo simulations. In particular, we had to consider an integrated lattice operator to decouple some uninteresting total derivative terms, and to optimize boundary conditions to maximize the (integrated) 1pt function of V , which was our Monte Carlo target.
Further boundary condition optimization is likely possible (see appendix D) and might allow to reduce the error bars in future studies.
The main limitation of our approach to measuring ∆ V is that while it decouples total derivatives of scalars, it does not do so for divergences of spin-2 operators. As a result we measure not where T is the lowest non-conserved Z 2 even spin-2. So, our result ∆ I = 6 ± 1 only implies a lower bound ∆ V 5.0 on the dimension of V . Still, the virial current value ∆ V = 2 is soundly ruled out by this lower bound. This confirms that the 3d Ising model is conformally invariant.
Now assuming conformal invariance, we know from the conformal bootstrap that ∆ T ≈ 5.51 [2]. This suggests that our measurement of ∆ I was dominated by ∆ T + 1, while V itself may be much higher. This scenario appears likely also in light of extremely high values of ∆ V in d = 2, 4 reported in the Introduction.
In this paper we have not carried out any correction-to-scaling analysis. It would be interesting to repeat the simulation in the Blume-Capel model which is in the same universality class as the Ising model but has a free parameter allowing to drastically reduce corrections to scaling [15].
It would be also interesting to determine or bound the dimension of V for the O(N ) and other models.
Finally, we would like to comment on the determination of ∆ V using the conformal bootstrap.
The numerical conformal bootstrap has determined scaling dimensions of about 100 operators of the critical 3d Ising model [2]. The operators which have been determined appear in the operator product expansions (OPEs) of σ × σ, ε × ε and σ × ε, where σ and ε are the lowest dimension Z 2 -odd and Z 2 -even scalars. The OPEs σ × σ and ε × ε, being OPEs of identical scalars, contain only operators of even spin. The OPE σ × ε contain only Z 2 -odd operators. The operator V , being a Z 2 -even vector, does not appear in these OPEs, and therefore it has not been so far probed by the conformal bootstrap. In the future, the OPEs σ × σ and ε × ε , where σ and ε are the subleading Z 2 -odd and Z 2 -even scalars, will hopefully be included in the bootstrap analysis.
These OPEs contain V and can be used to determine its dimension.
Of course, determination of ∆ V using the conformal bootstrap already presupposes that the model is conformally invariant. This has to be distinguished from the lower bound on V obtained in our paper, which is valid independently of conformal invariance, and so allowed us to test this property.
Note added. In the first arXiv version of this paper [10] the reader will find an appendix criticizing the argument in [11] for conformal invariance of the critical 3d Ising model. We consider the objections raised there still valid, and the rebuttal [12] unsatisfactory. However, we removed the appendix to keep the focus on the positive results obtained in our own work.
As a consistency check, we have determined ∆ V = 11 using an alternative method. We performed the conformal block decomposition of the four-point function 11 where G ∆,s stands for the conformal block of dimension ∆ and spin s (corresponding to the SO (4) irreducible representation ( s 2 , s 2 )). Again we find the first vector primary at dimension 11. One can also see that the vector primary operator we identified is parity-even. This follows immediately because parity odd vector primary operators cannot appear in the OPE of two scalars (like φ 2 and φ 4 ) in a parity symmetric theory. In addition, it is easy see that the vector operator contains 6 fields φ and 5 derivatives. 13 We also studied the conformal character decomposition 11 We normalized the operators φ 2 and φ 4 to have unit two-point function. 12 We use the standard conformal block as defined in [20,21]. 13 The φ content of each primary can be obtained by studying the partition function where r is a fugacity for the number of φ's in each local operator.
of the free massless scalar in d = 3. The lightest vector primary still contains 6 fields φ and 5 derivatives, which leads to ∆ V = 8 in d = 3.
The conclusion that the lowest Z 2 even vector primary has dimension 11 was reached independently by Marco Meineri [22]. He used a different approach, which also provides an explicit expression for this primary in terms of φ and its derivatives. In d = 4 − , this vector primary operator will get an O( ) anomalous dimension, computable starting from an explicit expression in [22]; this will not be done here.
One potential worry could be the recombination of this multiplet with a short multiplet when > 0. However, it is well known (see e.g. [23] for a discussion) that the only multiplets that recombine are the multiplet of φ with the one of φ 3 and the multiplets χ short 2+2n,n,n (conserved currents of spin 2n) with χ 3+2n,n− Notice that in all the above discussion we set ∂ 2 φ = 0 in 4d, eliminating operators involving the letter "∂ 2 φ" from consideration. When we go to (4 − ) dimensions, we will have the equation of dimension 14 + O( ). This operator is not a primary [26], so the lowest evanescent vector primary is still somewhere higher. We conclude that the evanescent operators cannot compete with the 11 + O( ) primary that we found above.
The vector quasiprimaries can also be found by studying the (global) conformal block decomposition of a four-point function involving two different scalar operators. In 2d Ising, the simplest choice is (with ∆ = 1) and TT (with ∆ = 4). Such correlation functions can be easily computed using the conformal Ward identities. In particular, we obtained The conformal block expansion in the z,z → 0 channel is given by (shifts in the familiar exponents w.r.t. β/2 due to unequal dimensions of external scalars). This confirms that ∆ V = 14 in the 2d Ising CFT.

B. Comments on operator matching
Here we collect some well known facts about operator matching between UV theory and its IR fixed point. UV theory may be a lattice spin model, a field theory with cutoff, or a continuum limit field theory.

B.1. Matching in the lattice spin model
We consider first the lattice spin model case, and will explain the necessary modifications to UV representations. The lattice theory has lattice operators which form multiplets under the lattice symmetry group (cubic group). The critical theory is sometimes called CFT, but here we will avoid using this terminology since we don't want to assume conformal symmetry from the start.
The important point is that critical theory correlators are defined at all distances 0 < r < ∞, while correlators of the lattice theory are defined at discrete distances r a.
How to recover parameters of the critical theory in a lattice simulation? Two issues complicate this extraction. The first issue is that operators of the lattice theory, naturally given in terms of lattice variables, do not have well-defined scaling dimension, but should be thought of as linear combinations of such operators. The second issue is that the lattice theory, even with couplings finetuned to the second-order phase transition, does not sit precisely at the fixed point, but only flows to it at large distances. Let us consider in turn how these issues manifest themselves.
Consider the simplest lattice operator, spin S lat (x). We should expand it in critical theory operators. The appearing terms will have to be, as S lat (x), Z 2 -odd cubic group singlets. The expansion (sometimes referred to as matching) will have the form: There are infinitely many terms but we only wrote the first few representative ones. σ and σ are the first two Z 2 -odd scalars of the critical theory (of dimension ∆ σ ≈ 0.518, ∆ σ ≈ 5.29).
Derivatives of these operators with indices contracted so that they are scalars can also appear (∂ 2 σ being shown as a representative case). In addition scalar derivatives of tensor Z 2 operators are also expected to appear, the representative case being the divergence of some Z 2 odd vector On a lattice with spacing a, all coefficients A i in this expansion will be given by withÃ i a dimensionless number and ∆ i the dimension of the critical operator multiplied by the corresponding coefficient. On a lattice of unit spacing they will be simply O(1) numbers.
With the expansion (B.1), correlators of S lat (x) in the lattice theory, can be matched with sums of correlators of operators in the critical theory. For example, for the 2pt function we have: Here the correlator in the l.h.s. can be measured in a lattice simulation, and by this equation it should be equal to a sum of critical correlators in the r.h.s. Consider for example correlators in infinite volume. The critical theory correlators are expressed in terms of scaling dimensions of the fields. For scalars: where 1 is just a normalization. For a vector operator we would have Here the constant α equals −2 in a CFT with R µ a vector primary, but in a scale invariant theory but non-conformal theory it could be different. Also in a non-conformal theory there could be nonzero 2pt functions between operators of unequal scaling dimension which then have to be added to the r.h.s. of (B.2). In any case, according to this discussion, and taking into account the expected size of coefficients A i , the r.h.s. of (B.2) contains a series of terms decaying with the distance as const.(a/r) p i where the powers p i are simply related to scaling dimensions of operators appearing in the r.h.s. of (B.1). We see that only dimensionless ratios of distances enter into this expression. If we go to distances r a, then the lowest power p 1 = 2∆ σ will dominate and the first correction will be suppressed by two more powers of the distance. The terms involving d µνλσ tensor will have nontrivial angular dependence, a sign of rotational symmetry breaking.
The leading such term will appear from the crossterm σ∂ µ ∂ ν ∂ λ ∂ σ σ and will be tiny, suppressed by 4 powers of the distance.
To complete the just given discussion, we need to address the above-mentioned second issue, taken into account by perturbing the action of the critical theory by irrelevant operators. More precisely, we can describe the system by the action This indicates that the r.h.s. can contain vector critical operators (V µ ), derivatives of scalars (∂ µ ) and divergences of tensors (∂ µ T µν , excluding the stress tensor T µν as it is conserved), as well as 15 The idea of improved lattice actions is to use models that allow to tune to zero the couplings of the first few leading irrelevant operators in (B.5). For example, the Blume-Capel model used in [15] allows to set g1 = 0 thus removing the leading corrections to scaling due to ε . 16 The coefficients Ai are of course not the same as for S lat (x).
rotation-invariance breaking terms involving higher-rank tensors contracted with special tensors like d µνλσ , to get objects which transform correctly under the cubic group.
In the generic case we expect all A i = O(a ∆ i ) as for S lat (x). In the special case of O lat µ (x) being a total lattice derivative, 17 we will have A 1 = A 4 = 0 and only the terms like A 2 , A 3 could contribute.
We emphasize that all we know of the coefficients A i on general grounds is that they are O(a ∆ i ) numbers. There is no simple theoretical way to determine these numbers apart from a lattice simulation. All operators which are allowed by lattice and internal symmetries (and total lattice derivative constraints) will appear in the r.h.s. The problem of determining these coefficients is a "long distance" problem: it has to do with how the microscopic theory approaches the IR fixed point at long distances.

B.2. Matching in the lattice field theory
It is instructive to consider what changes when we replace the spin model by the latticized φ 4 field theory, defined by the lattice action For each value of the quartic coupling λ > 0 we can find a value of the mass parameter corresponding to a second-order phase transition. For this value m 2 * (λ) the theory flows at large distances to the critical theory, which does not depend on λ and is actually the same as for the Ising spin model. The operators of the UV theory can be then expanded in critical theory operators. For example, we can write an expansion for φ(x) of the same form (B.1) as for the spin operator S lat (x). The symmetry reasoning which led to this expansion remains the same, and the same operators will appear in the r.h.s. However, the discussion of the size of coefficients A i has to be slightly modified.
We say that the φ 4 theory is strongly coupled at the lattice scale if the quartic coupling λ is not small. The appropriate dimensionless condition in 3d is λa 1. 18 The effects of such largish quartic coupling are strongly felt already at the lattice scale (and a fortiori at all longer distance scales). Because of this, the RG flow will converge to the IR fixed point at distances r not much higher than a. The matching coefficients in the strongly coupled latticized φ 4 theory will thus be of the same generic size as for the spin Ising model, i.e. A i = O(a ∆ i ). 17 Total lattice derivatives are local operators that when summed over a region of the lattice, reduce to operators at the boundary of that region. For example, ∇µφ(x) = [φ(x + aeµ) − φ(x − aeµ)]/(2a) is a lattice derivative. 18 Notice that the lattice field φ(x) has dimension 1/2 like a free scalar field in 3d. This implies that the quartic coupling λ has mass dimension 1.
If on the other hand the quartic satisfies λa 1, the starting point of RG flow finds itself not far from the gaussian UV fixed point (UVFP). The RG trajectory can then be divided into two parts (see Fig. 6). In this case we say that the UV lattice theory is 'weakly coupled'. The first part of the RG flow happens in the neighbourhood of the UVFP. It corresponds to distances 0 , where 0 = 1/λ a. The second part starts at distances ∼ 0 where the flow transitions from the neighbourhood of free UVFP to the strongly interacting IRFP.

Rotational invariant theories
for weakly coupled flows we haveũ 2 ∼ λa ∼ a/ 0 1. Furthermore, because we tuned the mass term to its critical value we also haveũ 1 1. The first term breaking rotational invariance has The second part of the flow starts at the scale 0 = 1/λ a. Therefore, the scale 0 plays the role of UV cutoff for the second part of the flow. It is then useful to write u i =ū i to define dimensionless couplingsū i with respect to the UV cutoff for the second part of the flow. This givesū 2 = O(1) for the quartic coupling andū 3 ∼ (a/ 0 ) 2 1 for the leading irrelevant 19 The couplings of irrelevant operators that involve more than two powers of φ are also suppressed by the small parameter λa because at λ = 0 the lattice path integral is exactly gaussian.
coupling that breaks rotational symmetry. The second part of the flow can then be described using the action (B.5) with dimensionless couplingsg i defined by g i =g i . We expect We thus see that the second part of RG flow starts with some irrelevant operators in the action having dimensionless couplings much smaller than the other ones. This effect was absent in the spin lattice model case, where all irrelevant operators were expected to be present at the cutoff scale with O(1) coefficients in lattice units. As a consequence, rotation breaking in the IR, already small in the spin model case, will be even further suppressed in the weakly coupled lattice field theory case. Now let us discuss matching of operators, which also happens in two stages. First we expand lattice field theory operators into operators of the UVFP. E.g. we will have Coefficients of this expansion have a power series expansion in λ. For example we expect . Then we have to expand UVFP operators in IRFP operators.
This matching is done at the scale 0 . E.g. we have: Since this matching is done at the scale where the flow is strongly coupled, the coefficients B i cannot be easily predicted and are expected to be O(1). Combining the two matchings, we will get expressions for lattice field theory operators in terms of IRFP operators.

C. Possible boundary conditions
One can imagine modifying our setup described in the main text, by changing the boundary conditions at x 3 = 0, L − 1. The purpose would be to find boundary conditions which lead to an even larger f (t) and thus improve the signal-to-noise ratio. It makes sense to keep translation invariance in the x 2 direction, so that I(x 2 , x 3 ) is x 2 independent and can be averaged in this direction.
As discussed in the main text, we have to break the x 1 flip symmetry. One way to do this is to choose different boundary conditions for different parts of the x 3 = 0, L − 1 boundaries, depending on x 1 .
In addition to the free and fixed boundary conditions (b.c.) described in the main text, there are two other imaginable types of b.c. worth discussing.
The gluing b.c. changes topology of our manifold, by gluing one part of the boundary to another.
For example, one can imagine gluing the gray parts of the x 3 = 0, L − 1 boundaries in Fig. 1, instead of imposing the fixed b.c. there. In practice, gluing is achieved by identifying points pairwise or, equivalently in the large L limit, by creating links joining the points being glued. In the just mentioned example, we would be identifying points Gluing does not have to preserve order, for example we could have instead chosen to glue the gray parts of the boundaries while simultaneously flipping the x 1 coordinate. Such a reversed gluing would be a different boundary condition.
One can even glue parts of the same boundary, e.g. the lower and upper white parts of the Fig. 1 (again, in the direct or the reversed x 1 order).

C.2. Changing the strength of boundary interactions
We may change the strength of interaction among spins belonging to some part of the boundary to β bdry = β c . Two particularly interesting values of β bdry are as follows.
• β bdry = β sp ≈ 0.33302. This fixes β bdry to the value corresponding to the "special" boundary phase transition. Recall that the special transition separates the "ordinary" boundary behavior for which the boundary remains disordered at the critical temperature, from the "extraordinary" one when the boundary is ordered at the critical temperature. The ordinary (extraordinary) behavior is realized at β bdry < β sp (β bdry > β sp ). The β sp for the 3d Ising model given above was determined in [28]. Since the boundary points have fewer neighbors than the bulk points, β bdry = β c belongs to the "ordinary" phase, and this explains why β sp > β c .
• β bdry = ∞. This enforces that all spins are equal along a part of the boundary, which is the maximally efficient way to enforce the "extraordinary" boundary behavior. Notice that unlike the fixed boundary condition, the spins can still fluctuate between ±1, but only all at once. This difference may seem minor, but it has the following practical consequence.
The fixed b.c. can be used if the simulations are performed using the Metropolis algorithm, as in the main text. On the other hand, if the simulations are performed using cluster algorithms, it leads to lowering the acceptance rate since clusters which touch the boundary cannot be flipped. The β bdry = ∞ boundary condition does not have this difficulty.
There are many imaginable combinations of the four boundary condition types which break symmetries of the lattice in a way which makes f (t) nonzero. It is tedious to simulate one by one all possible combinations for the Ising model and see which one gives the largest f (t). It would be nice to have a way to guess a good boundary condition. A heuristic method is described in the next appendix.

D. Heuristic optimization of boundary conditions
Consider the free massless scalar theory on the cubic lattice, described by the action: We consider in this theory a lattice operator V lat µ given by the same equation (2.1) with φ(x) instead of s(x). We make a heuristic hypothesis that one can get an idea about the size of I in the critical Ising model by measuring the same quantity in the free scalar theory on the same cubic lattice. One motivation for this hypothesis is that in d = 4 the two theories are actually identical. We won't attempt to justify this hypothesis any further. It's amusing that empirically it seems to work. Once the b.c. is so heuristically guessed, the actual hard computation will be an honest Monte Carlo simulation in the 3d Ising.
To use the heuristic, we have to establish a correspondence between boundary conditions for the two models. This correspondence is as follows: 1. The free b.c. in the Ising corresponds to the Dirichlet b.c. for the free scalar. Indeed, the free b.c. in Ising leads to the "ordinary" boundary behavior, where the order parameter is effectively zero on the boundary [29].
2. The gluing b.c. in Ising clearly corresponds to the same gluing for the free scalar.
3. β bdry = ∞ for the Ising corresponds to imposing that φ(x) remains constant on this part of the boundary for the scalar.
5. The fixed 3d Ising boundary condition can be modeled by adding a constant magnetic field (linear in φ(x) term) on the boundary, pushing the free scalar in the needed direction.
We won't give full details on how one actually performs the calculation for the free scalar.
This calculation is inexpensive since one is computing a gaussian path integral. One constructs the lattice action, evaluates the Green's function, and finally evaluates the observable. The computation is done numerically and takes only a few seconds for a given boundary condition.
The most expensive step is the Green's function evaluation which requires to invert an L 3 × L 3 matrix.
After playing with the free scalar, we concluded that the boundary condition in Fig. 1 is particularly promising. Notice that since we have the same fixed b.c. on two parts of the boundary, and since we measure a Z 2 -even observable, for the purpose of the heuristic computation we could replace the fixed boundary condition with β bdry = ∞.
Before we discovered the heuristic optimization trick, we tried other boundary conditions in the 3d Ising, but they led to a smaller f (t).
We could have just postulated the boundary condition in Fig. 1, but we prefer to play in the open. This is because we have not performed exhaustive optimization. Even better b.c. likely exist, and our heuristic may be helpful to search for them.