A Monte Carlo method for critical systems in infinite volume: the planar Ising model
Victor Herdeiro and Benjamin Doyon
Department of Mathematics, King’s College London
Strand, London, U.K.
email: benjamin.doyon@kcl.ac.uk
In this paper we propose a Monte Carlo method for generating finitedomain marginals of critical distributions of statistical models in infinite volume. The algorithm corrects the problem of the longrange effects of boundaries associated to generating critical distributions on finite lattices. It uses the advantage of scale invariance combined with ideas of the renormalization group in order to construct a type of “holographic” boundary condition that encodes the presence of an infinite volume beyond it. We check the quality of the distribution obtained in the case of the planar Ising model by comparing various observables with their infiniteplane prediction. We accurately reproduce planar two, three and fourpoint functions of spin and energy operators. We also define a lattice stressenergy tensor, and numerically obtain the associated conformal Ward identities and the Ising central charge.
Contents
1 Introduction
Critical models of statistical physics, such as the Ising model at its Curie temperature, offer a unique opportunity to study the emergence of largescale behaviours. In this context, criticality occurs only on very small subsets of the parameter space of a given model. At such critical points, a second order phase transition takes place and the behaviour of the model changes drastically. Power laws appear in longrange correlations and in divergences of the free energy and susceptibilities, and with the universal amplitudes, they are described by the powerful theoretical framework of conformal field theory (CFT) [1, 2]. These effects are characteristic of the emerging universal collective behaviours, such as the large fluctuating loops separating ordered phases (a precise representation of critical bubbles of the nucleation theory), which have received attention recently especially in the context of conformal loop ensembles [3, 4] and its relation to CFT [5].
Monte Carlo simulations of statistical models are essential in order to explore with precision and from first principles these critical behaviours, and especially the emerging largescale objects in infinite volume. However, criticality, because of scale invariance, presents its unique numerical challenges. First, critical slowdown increases, in the most basic algorithms, the time necessary to numerically build the full longrange correlations. Second, at criticality, generic boundaries have effects that propagate well into the bulk in the form of conformal boundary conditions [6]. This effectively precludes the direct Monte Carlo numerical simulation of partial, finitedomain configurations in infinite systems. That is, it is hard to simulate the marginals, in the sense of probability theory, of the fluctuating degrees of freedom lying in finite domains with respect to the infinitevolume distribution. We will refer to such marginals as bulk marginals; they encode all local information (all correlation functions and loop distributions inside a given domain) of the infinitevolume model.
In this paper we propose a solution to the latter problem, illustrating it in the context of the planar Ising model. The Ising model is a milestone in the history of statistical physics [7]. While it was first introduced in one dimension by Lenz, Ising’s doctoral work [8] showed it could not order at large scales on such geometries, and Peierls arguments [9] made it the first model shown to exhibit a phase transition for dimensions greater than two. This breakthrough essentially opened the field of phase transitions. Onsager’s seminal contribution [10] rigorously proved important aspects of the phase transition, in particular the powerlaw behaviour of the magnetization as function of the external magnetic field, by giving a general solution in two dimensions. Since then, it has been the subject of very extensive research, a great part of the interest being due to its integrability.
In two dimensions and with zero external field, the model is known to go through a second order phase transition for some latticedependent critical temperature. At this point, the response functions and correlation length diverge and the system shows an infinite sensitivity to external perturbations. The fluctuations happen on every length scale, this being associated to the emergence of scale invariance. In general, within the framework of quantum field theory (QFT), emeging scale invariance implies invariance under conformal transformations, and combined with locality, leads to an infinite number of constraints on the correlation functions, algebraically encoded via the Virasoro algebra and its representations (the primary operators). This is CFT. These continuum, universal models are characterized, in particular, by their Virasoro central charge , and have been solved (e.g. analytical expressions or precise formulations for correlation functions) under the requirement of minimality [11, 1, 2]. The critical Ising lattice is the simplest example of such models, with central charge , and set of primary fields given by (the identity field), (the “spin” field – scaling limit of the binary lattice variable), and (the energy field – scaling limit of the local energy density).
On the numerical side, the Ising model has also been the subject of extensive research [12, 13, 14]. It was discovered that near the critical point the autocorrelations of local updates in Monte Carlo Markov chains diverge severely as , where is the linear size of the lattice, whereas the divergence is only far from the critical coupling. This phenomenon is associated to critical slowdown and has been explained as the difficulty for local updates to build the lowest energy excitations of the spectrum. The discovery or engineering of updates based on clusters or lattice flips [15, 16] has been a breakthrough in studying numerically the Ising model at its critical point, and largely solved the problem of critical slowdown.
However, at the critical point, because of the divergence of the correlation length and the response functions, the system also becomes extremely sensitive to its boundary conditions. Essentially, the universal effects of boundary conditions propagate well inside the bulk in the form of the allowed conformal boundary conditions of the CFT model [6]. A direct consequence is that it is presently impossible to sample the infinite volume Ising model even on a restricted finite subdomain: there is no known procedure to obtain, in an efficient fashion, samples representing bulk marginals.
The aim of this paper is twofold. First, we propose a Markov chain Monte Carlo (MCMC) sampler (which we refer to as a UV sampler) that produces, to good accuracy, samples representing bulk marginals: finitedomain marginals of critical bulk configurations. It does so by using scale invariance in order to generate, up to corrections which we attempt to characterize, a boundary condition that represents the presence of an infinite space beyond it. This may be seen as a “holographic” boundary condition, encoding the information of the infinite space beyond the finite domain. The UV sampler does not require the infinite conformal symmetry emerging in two dimensions, hence is generalizable to higher dimensions. It can also be adapted to study the full universal region near criticality (massive QFT).
The well known transfer matrix methods, not based on Monte Carlo sampling, also allow for very accurate studies of various bulk observables in twodimensional critical models (see for instance [17, 18]). However, these are fundamentally based on the infinite conformal symmetry present in two dimensions, are limited by the exponenrial growth of the time as function of the system’s size, and seem to be, at least up to now, more limited in the scope of the observables studied, as they do not actually produce planar bulk marginals.
Second, we present an extensive study of local observables that are natural in the CFT context: numerically evaluating some important bulk two, three and fourpoint functions, constructing the holomorphic stressenergy tensor, and verifying the conformal Ward identities. For completeness, we also present a study of the cluster boundaries, extracting some of their main properties. This in part serves to verify the accuracy of the proposed Markov chain, and in part to provide numerical results that may help in further characterizing the Ising critical point.
The paper is organized as follows. In section 2 we describe the theory behind the proposed method and its explicit implementation as a UV sampler in the Ising model, and we provide measures of the accuracy of the bulk marginal obtained. In section 3, we report on an extensive numerical study of Ising bulk observables, including multipoint functions and the stressenergy tensor, confirming the agreement with known analytical results. We conclude in section 4.
2 A Markov chain for generating bulk marginals
Consider the critical Ising model on a finite lattice. As the volume of the lattice is sent to infinity, marginals on fixed finite domains tend to limit distributions: these are the bulk marginals. However, for any large but finite lattice, generic lattice boundary conditions have effects that propagate well into the bulk, because of longrange powerlaw correlations. Thus the infinitevolume limit is very slowly approached. By renormalization group arguments, the powerlaw effects of generic boundary conditions are universal: they can be described by conformal boundary conditions of the associated CFT. Informally, this means that a generic lattice boundary condition is seen, from a large distance , as if it were a conformal boundary conditions of the CFT model at some other, effective distance . Two characteristics describe these universal effects: the conformal boundary condition itself, and the ratio determining the effective distance at which the conformal boundary condition appears to be. These are the emerging universal properties of the lattice boundary condition.
Clearly, a clever choice of the lattice boundary condition might correspond to a that is very large. With such a choice, bulk marginals are easily produced (as with large, it is not necessary to actually take the infinitevolume limit of the lattice model). A natural candidate for such a lattice boundary condition is the marginal on a codimensionone closed surface of lattice sites that separates the fixed domain from the rest of a very large lattice. Because of the locality of the Ising measure, such a boundary condition fully encodes the rest of the lattice: it is a “holographic” boundary condition.
The proposed MCMC sampler uses conformal invariance and the renormalization group in order to implement the above idea, sequentially adjusting the lattice boundary condition so as to make larger and larger. The method is an iterative procedure, where each iteration is composed of two steps, which are meant to effectively zoom in, in the sense of the renormalization group, onto a smaller and smaller region around a central point. The first step is a blowup procedure, which takes a small central region, of linear size times that of the original region, and blows it up to the linear size of the original region. This is the step that zooms in onto a small region, or, equivalently, that sends the original boundary further away. The distribution of the resulting boundary spins is the holographic boundary condition, encoding its exterior. The second step is a rethermalization, keeping fixed the new boundary spins. This effectively propagates the new boundary condition towards the inside of the region, recovering the information that was absent in the original small central region due to the finite lattice mesh. The implementation of this method using a Markov chain Monte Carlo (MCMC) sampler will be referred to as a UV sampler.
A UV sampler provides a way of performing Monte Carlo simulations of bulk configurations directly on small lattices, without the costly step of thermalizing on large lattices. The method appears to offer a significant advantage over the naïve way of generating thermalized configurations on a large lattice and extracting a bulk configuration by discarding all but a small finite sublattice. In this paper we concentrate on the twodimensional Ising model, but similar ideas can be used in any dimensionality.
In this section we give more details concerning the theory of this procedure, its MCMC implementation, and a basic analysis of the quality of bulk marginals obtained.
2.1 Renormalizationgroup steps towards the UV fixed point
In Euclidean CFT on the plane^{1}^{1}1Some of the arguments presented here generalize to QFT., one may represent, in an appropriate quantization scheme, any region outside a simply connected region , along with any set of local observables lying outside , by a state on the boundary . This state has enough information to reproduce, for instance, any correlation function with other observables lying in :
(1) 
for and ^{2}^{2}2To be precise, on the righthand side of (1) one uses the representation of the observable on the quantization space.. See Fig. 1.
The state corresponds to restricting on , or “integrating out” the exterior . If is a disk, one may use the radial quantization scheme, where the plane is foliated by concentric circles and time parametrizes their radii. In this scheme, states are seen as linear combinations of field configurations on circles, and is such a state on the circle , obtained from the original correlation function by a decomposition of the identity on . But the concept holds for any simply connected region .
Likewise, in a path integral formulation, with measure over field configurations on the plane, we have:
By ultra locality , and by locality (here is the closure of ), where, because of the microscopic connections (nonultra locality), both terms in the decomposition of depend on the field at the boundary . Thus the integral factorizes into:
(2) 
where
(3) 
carries all the information about operator insertions and the Gibbs measure beyond the boundary, encoding it into a function of the field configuration at the boundary . The function is the “wave function” associated to the state .
Naturally, in general the state is not a conformally invariant boundary state: it is not invariant under conformal transformations on . With trivial operators , it is, however, Möbius invariant.
The above can also be done on any domain of definition instead of the plane , and with any boundary state on . The state will then depend not only on but also on the domain and the boundary conditions on . Let us specialize to trivial exterior observables . We write for the state on characterized by the “exterior information” (the information outside ): this is the boundary state at , which we understand as implicitly containing the information of the domain of definition itself. On a simply connected domain with some boundary state , we have , and then the restriction to can be written as
(4) 
The following “projection” property is a simple consequence of these definitions:
(5) 
The projection property is very natural from the viewpoint of the path integral formulation, where can be evaluated by two embedded path integrals on complementary subdomains.
Let us now consider scale transformations, and restrict ourselves to convex domains. The scale transform of the state by a factor is a state on , of the form . The scale transform is defined by requiring scale invariance . The above then simply gives
(6) 
(where for instance a domain transforms as , and scaling fields as , with the scaling dimension of ). Rewriting, we have
(7) 
Let us define the operation
(8) 
This operation is a restriction to a smaller domain, , followed by a scaling by a factor . The map acts on, and generates, states on . Let us take . We now show that its power has the effect of scaling out the domain by a factor :
(9) 
Here, for lightness of notation, the symbol in the square brackets represents some (unspecified) conformal boundary condition on the boundary of the domain of definition . Clearly, (9) is true at . Assume it holds for some . Then
(by (7))  
(by (6))  
(by (5))  (10) 
which shows (9) by induction on . Thus,
(11) 
where on the righthand side, the symbol in the square brackets represents the exterior information of the domain of definition being the infinite plane. Therefore we see that the infinite power of on produces a boundary state on representing the integrating out of the exterior on the plane .
The above has a clear interpretation. The map represents, as a mapping of states on , the operation of scaling up the domain of definition by a factor , thus sending further away the domain boundary and increasing the space exterior to . It can be seen as a renormalizationgroup (RG) step on states on : tells us how the wave function of a state transforms when integrating out a shell of finite thickness near the boundary of the region, and scaling out the result back to the original domain (see Fig. 2).
This RG step flows towards the ultraviolet fixed point, as the boundary , which has the meaning of an infrared cutoff, is sent to infinity while the observables in are kept invariant; by scale invariance, this is equivalent to zooming onto a smaller and smaller region inside the original domain .
We therefore have an infinite chain of states on , starting with some conformal boundary condition , and whose step is an RG transformation towards the ultraviolet (see Fig. 3):
(12) 
This chain converges to the state . This is the state obtained by integrating out everything up to infinity, thus a bulk marginal. The chain allows us to reach this marginal from a finite domain by performing (an infinite number of) finitescale RG transformations.
2.2 Towards the UV fixed point on lattices: blowup and rethermalization
The above CFT description of boundary states representing exterior data has a counterpart in local statistical models. As an example, consider the twodimensional Ising model. Its configuration space is the set of all functions from a regular planar lattice embedded into the plane (say the triangular lattice), to the set . The (unnormalized) measure on a configuration is
(13) 
where is the set of edges of the lattice (and is the edge joining to ). For any domain (open set) (we assume that domain boundaries do not intersect any vertex of the lattice), we denote by the restriction of to the set of all vertices lying in , and by the restriction of to the set of “boundary vertices” (a boundary vertex is a vertex outside which is connected to at least one vertex in ). We will also use the notation . Because of the regular, planar structure of the lattice, the set of boundary vertices indeed stays uniformly near to the boundary for any domain (open set) . We also denote by the set of edges that have nonzero intersection with . We may then define the probability measure on any domain , with boundary conditions fixing , by , and the induced probability measures on subsets is .
Let be two domains, and consider the restricted probability : the probability of given . Since the measure factorizes on the set of edges (that is, the measure is a product of factors, one for each edge), this satisfies a “domain Markov property”: . That is, the information of , rather than the full , is sufficient in order to deduce the marginal on . Then, if is a random variable supported on and a random variable supported on , the expectation of on the domain with boundary condition can be written as
(14)  
where
That is, can be evaluated as an expectation on , with boundary condition on determined by the “wave function” . This is the lattice counterpart of (1), or of its pathintegral version (3)^{3}^{3}3Note that in this discrete version of the path integral formulation of the boundary state, the subtlety about and both depending on the boundary field becomes somewhat clearer..
This naturally suggests that we define the expectation with boundary condition (a function of ) as
(15) 
We then have the counterpart of (4), for any variable lying in ,
(16) 
where the new wave function is
(17) 
It is clear, from this, that the projection property (5) holds.
The operation of dilation in the discrete case is ambiguous. Let . Given a measure for the subset of the lattice lying in , a natural definition would be to construct a measure for the subset of the lattice lying in by attributing spins from as for all on the lattice, where represents the site on the lattice nearest to the coordinate . The main problem^{4}^{4}4Another potential problem is that may lie outside , but this small perturbation of the boundary is not important. is that this leaves “holes”: unattributed sites. This is a fundamental property of discrete systems: there is more information on larger region, hence the operation of scaling by a factor greater than one needs to be supplemented by additional information^{5}^{5}5The inverse definition does not leave holes, but it is effectively just a special prescription for filling in the holes, not expected to be particularly good as it tends to increase correlations (creating smallscale lumps)., which affects the definition of scaling. Certainly, all correlations of spins on the attributed sites are unambiguously the same, up to scaling, as those of their preimage, but other correlations depend on the choice of filling prescription. In order not to affect largedistance behaviours, the filling prescription should be local ( determined in terms of and few of its neighbours). Formally, we need such a prescription only for the boundary spins, as we need to describe dilations of states. Given a prescription that determines in terms of (or in terms of the spins on a fewsite neighbourhood of ), we define the scaled state as
(18) 
Suppose that the measure is chosen to be critical (for instance, the inverse temperature in the Ising model measure (13) is taken at its critical value, which depends on the lattice chosen). In this case, there is indeed scale invariance in the lattice model for largedistance observables: scale invariance is a property of emerging collective behaviors. Therefore, although it is not expected that the wave function , as a function of the spins , satisfy strict scale invariance properties, the largedistance features of the function are scale invariant. In the context of evaluating observables that do show scaleinvariant properties, we expect these largedistance features to dominate. Since the prescription only affects local features, lattice scale invariance implies
(19) 
In terms of Fourier modes along , we would expect the smallscale corrections only to have significant Fourier components at large frequencies.
The combination of (18) and (19) gives rise to (6), up to smallscale corrections. Thus, defining the RG operation in the lattice case as in (8), the derivation (10) can be reproduced, and we conclude that the chain (12) is a Markov chain for generating marginal configurations of a bulk critical lattice, with microscopic modifications due to the smallscale corrections.
The numerical implementation of each iteration of the Markov chain (12) involves two steps.

Rethermalization. Say the state is reached after iteration. For the next iteration, first the map is performed, and this involves a rethermalization step. Indeed, this requires extracting a spin configuration on , and thus requires propagating the boundary condition from to . This can be performed by thermalizing the model in with boundary : generating a typical configuration in with that boundary condition, from which the spin configuration on can be extracted. The thermalizing process may be done by a multitude of available numerical algorithms, and the details will be explained in subsection 2.3. This part of the numerical implementation is time consuming, but can be as precise as necessary, its accuracy being essentially limited only by the small fluctuations associated to microscopic thermalization.

Blowup. Second, a blowup step is performed. This step is the implementation of . This is very fast, but involves an adhoc prescription for fillingin the holes, as discussed above. Here we use the extra information gained by thermalizing on all of in the previous step, and implement this blowup not only for the boundary spins on but for all spins on . The result is a configuration on which has the correct boundary condition , and which is relatively near to the correct thermalization on the interior of , preserving the correct largedistance correlations. This gives an initial condition in that accelerates the subsequent rethermalization step in the next iteration of the chain. The subsequent rethermalization step is then mostly meant to reconstruct the short distances correlations damaged by the introduction of the holes.
It is a nontrivial matter to understand and assess the effects of the smallscale corrections and of the prescription for fillingin the holes in the above algorithm. We discuss these aspects in subsection 2.5.
2.3 Implementation
In this section, we present our proposal for a Markov chain Monte Carlo sampler implementing the function introduced above in the specific case of the lattice Ising model. For definiteness, we choose the triangular lattice, and the sublattice is rectangular shaped, of size . We note that we have also implementated the algorithm and performed verifications on the square lattice.
Our procedure is the following:

At the starting point the system has all spins pointing up.

We treat the finite lattice as a torus and apply periodic boundary Wolff cluster flips [16]. With no external magnetic field the unique coupling is the parameter in (13), which is taken to be the infinitevolume critical value, here equal to
(20) Our simulation uses a value close enough to have a correlation length . This step is meant to bring the lattice close to the critical bulk with relatively little computational effort.

We then start a cycle by applying a dilation procedure. The “discrete” dilations will expand a central fraction of the lattice to the whole lattice. For a dilation by a factor , the spin assigned to the site of coordinates will be the value of the previous spin at location:
(21) with the floor function, and where is a uniformly distributed random parameter fixed for a given dilation (see the motivation in the caption of Fig. 5). Graphically and at the microscopic level, the operation looks like Fig. 4.
On much larger lattices, this discrete dilation is illustrated by Fig. 5.
Figure 5: A similar chronology as in Figure 4. From the first to the second graph we apply a dilation on a 512x512 lattice. The second graph shows the patterned distribution of the holes: they are distributed along a noisy square mesh, the noise coming from the random pick of a target spin when sharing its antecedent. The presence of these patterns motivated us to introduce the random parameter: for different values of this parameter the mesh of holes will be slightly displaced horizontally and vertically, this way uniformizing the distribution of holes on the boundary when averaging over and restoring an “average translation invariance”. In the continuous limit this (discrete) operation converges naturally to the dilation we are familiar with in CFT.
On a discrete domain there is always more information on than on and some spins may share the same preimage (21). We insist on not adding spurious information, and impose a onetoone rule: only one of them, picked at random after prioritizing boundary spins, will inherit the spin value localized at (21). Obviously this step will leave holes or unassigned spins. In order to fillin the holes, we apply a heatbath assignation: holes are filledin according to the Ising measure induced by their neighbouring spins. This prescription is meant to restore the first neighbour correlations close to the average bulk value, see Fig. 4 for a step by step graphical explanation and some additional details.

The second step of the cycle is a rethermalization. As the system is not precisely at criticality (which is impossible on a discrete lattice), the above dilation, which can be seen as an RG transformation towards the UV, took it further away. In the rethermalization, we thus apply Monte Carlo evolution steps with the aim to bring the sample closer to criticality. This rethermalization propagates the information contained in the new, dilated boundary condition towards the inside of the lattice. Our choice of evolution algorithm here is that of fixedboundary SwendsenWang (SW) flips, see Appendix A for implementation details and motivations.

Finally, we repeat the last two steps for as long as desired. As we show in the next subsection, three cycles seem to be enough to generate the correct bulk spin correlations.
We note that criticality of the initial configuration, obtained at the second point in the above procedure, is not essential for the procedure to work. Indeed, if, before the dilation, the sublattice exhibited longdistance correlations up to distances of order , these correlations are stretched by the blowup steps up to orders (and, as mentioned, the rethermalization steps are meant to recover the correct shortdistance correlations lost by the introduction of holes). From a massive QFT point of view, after a large number of iterations, the mass gap has been closed. Indeed, the philosophy of the proposed Markov chain is to aim at scaling down or reducing any energy scale, such as the mass gap or the inverse effective length scale of the whole system, in order to reach a region where the theory exhibits full scale invariance. It has been proven that in two dimensions scale invariance gets promoted to conformal invariance [19, 20]. In the next section, after introducing local field observables, numerical checks of the conformal invariance of the sample will be presented.
2.4 Quality of bulk correlations
As a measure of the quality of the sublattice marginal, we propose a quantity which addresses the question as to if bulk correlations have been established. Consider the statistics of the product for positions and in some central subdomain of the sublattice. In any sublattice marginal of an infinite system, this statistics will generate, at large enough distances, an exact powerlaw correlation. The spinspin correlator in the planar Ising model at criticality is well known to display the power law profile
(22) 
at large enough . We may then do a leastsquare linear fit of the set of points obtained by plotting with respect to . Two important data emerge: the fitted power , and the leastsquare uncertainty . The quality of the power law is measured by the smallness of . Contributions to include both statistical noise, as well as spurious correlations leading to systematic departures from a power law. Suprious correlations may be either of universal nature, such as those induced by a conformal boundary or a periodic boundary condition, or nonuniversal, such as microscopic effects. For the spin variable , microscopic effects are minimal, and the large number of samples mean that statistical noise is also small. Hence we propose the quantity
as a simple measure of the quality of the bulkmarginal universal correlations. A high number represents a good quality. See Fig. 6 for the numerical results on the value . The increase in quality shows a decrease of powerlaw fitting uncertainty by a factor of about 10. After three cycles, an optimal quality seems to have been reached. Comparison with the periodicboundary case indicates that indeed is dominated by universal correlations, and not by noise or microscopic effects.
More precise stepbystep numerical results can be visualized in Fig. 7, where the power of the power law and fitting uncertainty are displayed along the timeline. The power measured for the torus indicates a slower decay of correlations on the torus (it is clear that the torus has an excess of correlations over larger distances, coming from information going around it), but after three cycles it agrees very well with the expected bulk value. We observe that each dilation shows an abrupt jump in the direction of a lack of correlations, as the holes introduced are at best correlated to only first neighbours and thus noisy from a longer distance point of view.
Two remarks are in order. First, although the power law behaviour (22) of the correlations is expected to be true for large separations , while the short distance behaviour would be predicted to be dominated by microscopic effects, we will see later that the spin operator  the most local operator on the lattice  shows essentially no microscopic effects in its correlation. This justifies why we measure its correlations starting at distance 1. Second, we are restricting ourselves to observations on the center of the sublattice. As we will see later and as we have argued above, the borders show some departure from the bulk expectations, which we will attempt to characterize.
As a second check of the quality of the bulk configurations obtained, we consider the average of the energy density, defined as:
where is the number of lattice sites. To the best of our knowledge, the expected value of is not known in the bulk planar triangular lattice. We may obtain an approximate value by measuring it on the torus and extrapolating (), giving us:
(23) 
The numerical results are displayed in Fig. 8. We clearly see, again, that the dilation induces a jump in the value of the energy density: even though we are using a heatbath assignation method, there are holes that happen to be next to each other (recall that here so that at least of the sites post dilations are unassigned). The jump happens in the upward direction, because of the decrease of (here represents any vector pointing to a nearest neighbour), which is a loss in nextneighbour correlations. Despite that jump, the evolution very quickly stabilizes, and after three dilations we can already see a convergence to a thermalized value of sitting on the critical value (23).
It is interesting to remark here that in both spinspin correlations and energy density average, the thermalization seems to happen passed the third dilationrethermalization step: as if first neighbour () and further neighbours () thermalize almost simultaneously. Both sets of data show how dilations and lattice flips work hand in hand in the algorithm in order to wash away the initial boundary information and to bring the sublattice distribution closer to a marginal of the planar Ising model.
We have not observed any significant autocorrelations between realizations separated by a single dilationrethermalization cycle. This indicates some effectiveness of the Markov chain at sampling the indenpendent bulk marginals; although autocorrelations are likely to appear as .
Different runs were also made for other values of . The observed convergence is qualitatively unchanged, up to the general rule that the number of dilations needed to achieve mixing to a bulk marginal increases as lambda decreases (no trivial formula seems to describe this relation).
The same Markov chain has been implemented on a square Ising lattice at its critical coupling. The results on the monitoring of the convergence to a bulk marginal are very similar. This indicates that, correctly implemented, the map is a universal method to sample finite domains of a critical lattice system.
Finally, we remark that many “arbitrary” choices were made in our implementation of the Markov chain, including that of the starting point, of the holefilling prescription and of the evolution algorithm. These choices are the result of many trials and errors and were motivated by the criteria of computational speed and quality of sampling.
2.5 Remnant boundary effects
Residual effects from holes being assigned on the boundary cannot be fully corrected by any prescription or rethermalization. There are thus remnant boundary effects, which it is important to attempt to characterize.
Conformal boundaries (infinite magnetic field or infinite temperature applied on the boundary spins) or periodic boundaries are known to induce a departure of the average energy density, decaying as a power law (with decay exponent 1) in the distance to the nearest boundary or as a power law in the lattice size respectively. Our data shows a behaviour similar to these “canonical setups” with the boundary condition (after the step) inducing a power law on the expectation value of the lattice operator
It is detailed in (24) as the fluctutation of the lattice energy operator, which renormalizes into the energy operator from the Ising Kac table. On a planar lattice we would expect .
Numerically, see Fig. 10, we observe that this induces an expectation value on the lattice energy fluctuation field . This is expected to be a consequence of the unthermalized spins on the boundary (the number of which is proportional to the dilation parameter ). Indeed, recall that the blowup procedure included a prescription for holes left within the region (after blowup from ). The prescription used is to locally thermalize using the information of the assigned neighbours. Links to unassigned neighbours (holes) are effectively set to 0. Holes are thus assigned random locallythermalized values, but with a modified Hamiltonian, corresponding to the introduction of link defects (see Fig. 9). In the interior of , this is of little consequence because the rethermalization step guarantees that shortdistance correlations are correctly reconstructed and these defects are washed away. However, at the boundary, these defects remain. Since interior rethermalization is performed with fixed boundaries, these defects are expected to have an effect on interior distribution.
The profile we measure for this expectation value is a decaying power law in the distance to the closest border. Interestingly, the decay exponent seems to be dependent on the latest dilation parameter, and to take a continuum of values. For the values of taken, the exponent always lies between 1 and 2. Since it is larger than 1, the boundary effects are smaller than those coming from a conformal boundary condition. The decay exponent increases, and the amplitude decreases, as is brought nearer to 1, and these effects can be made to be very small at distances greater than about 20 sites from the boundary.
It has been observed in different setups [21, 22, 23] that the critical Ising model with a line of defects shows a continuous spectrum of perturbations depending on the defect strength. This defect strength can be a departure from the critical coupling such as in [21]. In this light, it is tempting to identify our observations with this phenomenon: the line of defects can be interpreted as the sublattice boundary and the density of cut bond as a defect strength, an effective departure from the critical coupling.
We also provide in Appendix B a basic CFT analysis of the effects of the lattice on the Markov chain. This analysis is expected to give only a maximal value for the decay exponent of the boundary effect, predicted there to be 2 (in agreement with the observed values, all lesser than 2).
3 Numerical verifications
3.1 Fractal dimensions of the loops
At criticality, the system is scale invariant and it is known that the clusters of equal spins and their boundaries form fractal sets [24]. The values of the fractal dimensions are well known from Coulomb Gas techniques or more recently derived from SLE arguments [25]. The possibility of measuring the effective fractal dimension of the loops’ boundaries and of the domains lying between the loops offer a good way of verifying that the samples produced from the proposed Markov Chain are indeed critical Ising samples. We emphasize that these are effective fractal dimensions, as objects defined on a finite lattice cannot be rigorously fractal; nevertheless they still very nearly exhibit, for instance, fractional scaling behaviours.
The sample used for the estimations below were taken in a pool of size sublattices of size 512x512. Each was generated by the Markov chain detailed in section 2, stored after each completion of the rethermalization steps. The dilation parameter used along the chain was . The two fractal dimensions – of the cluster boundaries and of the clusters themselves – were checked independently employing two different methods.
Box counting
Box counting is an effective method for computing the fractal dimensions of any object by measuring the scaling exponent of the number of its intersections with a mesh of varying size. For sufficiently small meshes, it is expected the the number of intersections scale as , where is the mesh size and the fractal dimension of the geometric object (see Fig. 11). In the sampled sublattices, looking only at the loops with diameter (defined as the maximum distance between two polygon edges) greater or equal to 40 lattice units, we find the average value of the fitted exponents to be . It is known from that the fractal dimension of Ising is given by with . Our estimation is well centered on the theoretically known value although with a quite large uncertainty.
By the same algorithm, we estimate the fractal dimension of the cluster domains to be , to be compared with the theoretically predicted value: . Here we selected clusters whose exterior bounding loop is of diameter at least 30. Again, our estimation suffers from a relatively large fit uncertainty. The culprit seems to be in the fact that box counting method offers a rather poor powerlaw aspect, with important statistical deviations due to the finiteness of the mesh size.
Bulk finitesize scaling
The same exponents have long been evaluated on the torus by the use of finitesize scaling. This is based on the assumption that the mean looplength and clustermass of the longest loop or heaviest cluster, respectively, scales as a power law in the periodic lattice size , with exponent the associated fractal dimension, at least at leading order [26].
Focussing on the fractal dimension of the length of loops, on a torus of size x, we have the following:
where runs over the set of closed loops in each toroidal lattice.
To our knowledge no similar estimator is known inside the bulk. By scale invariance arguments, we should expect that an estimator looking for the longest loops on a finite x sublattice of the infinite planar lattice should display the same power law at leading order, as the fractal dimension is not dependent on the boundary conditions. However, as we measure only the loops entirely inside the subsection, for small values of L, the impact of excluding the loops touching the borders becomes significant. Our estimation for was unsatisfying. To counter this issue, we decided to refine the estimator by introducing an arbitrary parameter . The recipe is the following: in order to avoid the small boundary effects of the bulk marginal, we cut a subsection inside a larger bulk marginal of size . The subsection is of size for and , and in it, we run over all boxes. The estimator will look at the longest loop fitting in the running box inside . For a fixed value of , we vary and expect to have an average maximal length scaling as .
This new estimator, named , should factorize into functions of and separately. Indeed, under a scale transformation , we should have:
constraining:
The evaluation of this new estimator for the length and mass exponents – on the same samples on which we ran the box counting estimation – using a value gives Fig. 12.
Just as for estimations using box counting, we find agreement with the theoretically known value for the fractal dimension of the boundaries, the progress is in the reduction of its uncertainty by one order of magnitude. Regarding the mass exponent, we find here better agreement with a fitted value of . This is a significant improvement from the box counting estimation, and, to our appreciation, sufficient on its own to convince of the relevance of adapting finitesize scaling estimators to bulk subsections. Here as well, we gained an order of magnitude on the precision. It is interesting to stress that this implementation of finitesize scaling in the bulk seems to be matching the precision magnitude of finitesize scaling on the torus [26].
These results on their own they are not enough proof that our samples have features close to those expected in the bulk, as fractal dimensions are not expected to be affected by conformal boundaries. They are only verifications that the samples are very close to the critical point. The next subsections provide indepth numerical analyses of proper bulk properties: bulk correlation functions.
3.2 Two and threepoint functions of spin and energy fields
As emphasized in subsection 2.4, twopoint functions are good indicators of how close to critical Ising bulk marginals the MCMC sampler gets. Here we look closer at the correlations over larger distances.
Recall that CFT scaling fields have twopoint functions of the general form
where is the spin of the field and its scaling dimension (here and ).
First consider the spin variable . This scales to a CFT primary field as
The results for the absolute value of its twopoint function are displayed in Fig. 13. This was obtained with up to 200 lattice spacings in every possible direction, and averaging over approximately 20 000 samples of size 512x512, taken after the rethermalization steps (after having achieved complete mixing) and separated by at least one dilation. Excellent agreement is found with a power law behaviour, and one can see that the results are much nearer, especially at large distances, to those expected for bulk marginals than equivalent results numerically evaluated in a torus geometry.
We also verified bulk rotation invariance. The spin variable is spinless, , by definition (more precisely, it is invariant under the lattice rotation symmetry group; of course this is expected to be promoted to full rotation invariance in CFT). However, discrepancies from bulk marginals may result in residual boundary effects from the squareshaped boundary. We estimated such effects by fitting the angular dependence of the correlator onto a cosinus function in order to estimate an effective . We looked at the correlator in three different directions (which are related by lattice rotations), followed by a reweighting of the initial data as to suppress the dependence on the distance . The fit gave the estimate , in good agreement with a theoretical value of 0.
The energy density operator may be defined by:
(24) 
where , that is, the sum runs over the positions that are neighbours of the site . The second term aims at suppressing the nonuniversal expectation value of the first term^{6}^{6}6The authors have not found a reference for this value in the literature on the Ising model, but are suspecting an exact value equal to 4., so it behaves as a CFT primary operator at large distances,
The energy field has spin 0 and scaling dimension 1. We analyzed the twopoint function with data collected over the same samples as those discussed above. We restricted to separations smaller or equal to 50 lattice units as beyond these, the noise was becoming dominant (it seems that the variance on the energy correlations is stronger than on the lattice spin correlations). We also excluded from the fitted data all distances smaller than 6, as they seemed to suffer from some microscopic effects. Our fit gives a scaling weight of , falling satisfactory close to the theoretical value of . See Fig. 14.
We also repeated the threedirectional fit we had done for the lattice spin variable, and could fit an effective spin for the energy density to a value of ; once again in good agreement with a theoretical value of 0.
From the powerlaw fits, the nonuniversal normalization factors and are evaluated to:
(25)  
(26) 
This allows us to estimate a first “dynamical” element of the CFT data (quantities that are not directly constrained by behaviours under global conformal symmetry): the structure constant .
We evaluated the correlator in the limit where the position of the first spin variable is close to the position of the energy variable. A result of an operator product expansion in CFT is:
with and fixed and where is the lattice spacing.
Looking at this observable on horizontal directions in 512x512 lattices for , excluding distances under 30 lattice units from the boundaries in every direction, and running over 5 000 lattice samples, gives the output displayed in Fig. 15. There, the ratio with the twopoint spin correlator has been taken, thus taking away the factor . The power law behaviour is obvious in this loglog plot, from the shortest distance (that is, it seems the requirement is satisfied even at distances ) and even up to large distances. Fitting on a power law gives a scaling weight of , for the remaining power of , once again in satisfactory agreement with the theoretical value.
The fit of the offset , with the previous estimations of , allows us a numerical estimation of the structure constant , in good agreement with the theoretically known value of (see e.g. [2]).
This coefficient had already been derived numerically in a Monte Carlo estimation [27], though using a different method with free boundaries and a sizedependent coupling tuned to exhibit scale invariant correlations. To our knowledge, our estimation is the first to be directly in bulk marginals, and brings one additional order of magnitude in precision to the numerical check.
Another interesting element of the CFT data is the threepoint energy field coefficient . This vanishes because is odd under the duality operation of the Ising model swapping the spin and disorder fields. This duality maps the hightemperature expansion of the partition function to the lowtemperature expansion. The critical point is selfdual, thus giving rise to an extra symmetry that is not apparent in the Ising measure or the MCMC sampler. This duality is broken by the initial state of the Markov chain, and thus needs to be built by its mixing. We also note that the variable transforms multiplicatively under duality only after appropriate shifting by the nonuniversal expectation value (24), again pointing to the nontriviality of this extra symmetry. The check consists of a fit of the lattice energy density three point function on:
(again with , and ), and using insertions with distances in a range from 7 to 40. This gives us the following numerical estimation:
Fitting a vanishing coefficient with more precision is a real challenge as it is dominated by noise. Here, our best effort gives us an estimation quite supportive of the duality invariance constraint.
3.3 The holomorphic stressenergy tensor and the conformal Ward identities
One of the most fundamental fields in CFT is the holomorphic stressenergy tensor (see e.g. [2]). Its existence indicates the presence of local conformal invariance: the fact that Ising measures defined on conformally equivalent domains (not necessarily simply connected), when conformally transported, give rise to the same scaling limits (see [1, 2] for a fieldtheoretic explanation of this relation, and [3, 4, 28] for a conformalloopensemble demonstration).
The main property of the stressenergy tensor is that it has dimension and spin . It is in fact the leading coefficient with spin 2 in shortdistance expansions of generic observables in CFT. For instance, the spinspin operator product expansion, in the identity operator channel, gives [2]:
which can be inverted by means of a Fourier transform:
(27) 
One can also argue that the holomorphic stressenergy tensor should be expected to be reproduced by the most local spin2 lattice variable.
In order to define the stressenergy tensor on a triangular lattice, we may thus extend the definition used in [29] from square to triangular lattices:
(28) 
Note that here (and contrary to the square lattice case) the variable takes complex values.
Up to a normalization, this definition generalizes (27) to the lattice, as the integral of becomes a sum over the symmetry directions. Since, as we saw, the spinspin correlator exhibits a CFT behaviour already from a distance of 1, it indeed appears to be sufficient to take a sum over first neighbours in order to implement (27) on the lattice.
A more refined definition might be obtained by using instead the nexttonearest, or farther, neighbours. The price would be to losing locality, and thus increasing microscopic effects and requiring larger distances in correlation functions in order to reproduce CFT behaviours. However, clearly definition (28) has spin 2 under the lattice rotation symmetry group, and it is the most local spin2 variable that can be constructed. Independently from its relation with the CFT formula, this definition may therefore be expected to reproduce the holomorphic stressenergy tensor. As we see below, this is indeed the case.
The relation with the CFT holomorphic stressenergy tensor is then expected to be
(29) 
We repeated the estimation for the power law of the absolute value of the twopoint correlator . This numerical estimation was significantly more challenging than in the cases of the spin or energy variables, since its correlations decay much faster with distance (as ), the signal being rapidly dominated by noise. We had to extend to approximately 670 000 uncorrelated lattice samples. On these, we looked in a restricted subdomain excluding 30 lattice sites from the boundaries, and used operator insertions along three different directions . This choice was motivated by computational efficiency. See Fig. 16.
We clearly see a ‘bump’ above the fitted power law behaviour, very strong at distance 2 (note that and overlap by sharing a lattice spin) and disappearing for distances larger than 7. Because of the difficulty we had in fitting this data on a power law, we decided to add a prescription in the fit and only include the set of points inducing the smallest uncertainty on the fitted exponent value, or in other words the distances for which the power law behaviour was the most manifest. This gave us points at distances 8, 11 and 12; for a fitted exponent value of . This is in satisfying agreement to the CFT value.
We also repeated the estimation of the spin . Instead of looking at directions , we decided to perform this consistency check by looking at its correlation to its neighbours in the directions  including neighbours with distance up to 11. Splitting the real and imaginary part and fitting them on a cosinus and sinus respectively gave us the numerical spin value of , in good agreement with the CFT expectation.
3.4 “Numerical” central charge
Although not identifying it uniquely, the central charge is a key characteristic of a CFT model: for instance, its monotonic behaviour under RG flows has clear physical interpretations, and its involvement in the Virasoro Algebra encodes how descendant operators transform under conformal mappings. However, numerically, it is harder to unravel as it is not a direct observable.
Different methods for computing it have been used in the past: for instance, putting the theory on a cylinder and looking at the scaling behaviour of the free energy, or, more recently, using the twist operator in order to construct conical singularities (related to the entanglement Rényi entropy).
Here, we derive it in bulk marginals simply by continuing our study of correlation functions. It is known from CFT that the stressenergy tensor has the following two point correlator:
In order to use this formula, we must numerically estimate the normalization constant in (29). We therefore first look at the stressenergy tensor insertion within a twospin correlator, which by Ward identities must obey:
where is the scaling weight of the holomorphic part of the spin operator.
In order to evaluate , we assumed invariance under translations and restricted ourselves to horizontal insertions (for simplicity, and thus only keeping a dependence on two distances). The lattice stressenergy tensor is inserted inbetween the spins with a maximal separation of 45 lattice units to each of them. The insertions were also constrained to a minimal distance to the borders of 40 so as to avoid boundary induced departures from the bulk marginal. The estimation ran on 3 000 samples of size 512x512. See Fig. 17.
The proportionality between the two graphs is apparent (as well as the overall lack of noise in retrieving this pure CFT result). We decided to fit the coefficient by optimization, looking at the minimum of:
with, again, , and . The sum is restricted to points with , thus excluding points more sensitive to noise, and also leaving out distances between insertions closer than 3 lattice units, as they are too sensitive to the microscopy. The above quantity has a paraboloid aspect (see Fig. 18) and shows a clear minimum at
(30) 
Using (25), this gives .
We may now evaluate the central charge by estimating the amplitude of the power law . In order to obtain the most precise fit of the corresponding offset in a loglog fit, we decided to run the calculation of the correlator from scratch on a set of approximatively 500 000 samples, excluding as far as 50 lattice units from the boundaries and looking at distances between insertions of 5, 6, 7 and 8.
The fit of the offset gives us
(31) 
(fitting between distances 4.72 and 8, for which the goodness of the fit is optimal). Using (25) and (30), we get the following numerical estimation for the central charge:
This estimate is not directly centered on the theoretically expected value , but is within one standard deviation. To our knowledge, at least one comparable numerical study shows a slightly more precise estimation of the quantity [29], by using a finite size scaling method on the torus.
3.5 Fourpoint functions
The quantities that encode most of the information of a CFT model are the fourpoint correlation functions. The lattice spin fourpoint function is known from CFT[1] to have the following analytical form, solution of the corresponding second order differential ‘null state’ equation specific to the minimal model describing the critical Ising model:
where , , and the “dynamical factor” is:
This “signal” does not decrease as quickly as the energy or stressenergy correlators but the challenge appearing here is the large freedom in which we can make our insertions. Our first numerical estimation was carried in the limit where two pairs of close insertions are well separated, here . The four point function simplifies to:
On a sample of 512x512 sublattices we estimated the quantity using insertions at least 60 lattice units from the boundaries and pairs within which the separation did not exceed 60 lattice units and respectively separated by at least 160 lattice sites. We also restricted ourselves to horizontal directions for simplicity. Fitting our data on a double power law gave us the following results for different cut off distances inside the two pairs:
CutOff  

60  
40  
20  
5 