Energy trapped Ising model

In this paper we have considered the 3D Ising model perturbed with the energy operator coupled with a non uniform harmonic potential acting as a trap, showing that this system satisfies the trap-size scaling behavior. Eventually, we have computed the correlators $\langle \sigma (z) \sigma (0)\rangle$, $ \langle \epsilon (z) \epsilon (0)\rangle$ and $\langle \sigma (z) \epsilon (0)\rangle$ near the critical point by means of conformal perturbation theory. Combining this result with Monte Carlo simulations, we have been able to estimate the OPE coefficients $C^{\sigma}_{\sigma\epsilon}$, $C^{\epsilon}_{\sigma\sigma}$ and $C^{\epsilon}_{\epsilon\epsilon}$, finding a good agreement with the values obtained in [1,2].


I. INTRODUCTION
In recent years, conformal data for several CFTs have been determined thanks to the conformal bootstrap program [2][3][4][5][6]. In addition to that, combining this numerical high-precision technique to analytical methods developed in the framework of Conformal Perturbation Theory (CPT) [7][8][9], it is possible to determine the behavior of the off-critical correlators of many different systems. This approach has been applied successfully to the well known 3D Ising model, by adding perturbations proportional to the spin and the energy operator [10,11].
Starting from a slightly different perspective, CPT can be also combined to Monte Carlo simulations to get insight both on the behavior of the correlators outside the critical point and on the CFTs data at the critical point. In [1], the author followed this approach to study the Ising model perturbed by a confining potential coupled to the spin operator. This model is particularly interesting because the behavior of the 1-point expectation values can be argued just by applying simple renormalization group arguments [14][15][16], and depends only on the trap potential parameters (trap size scaling (TSS) behavior). Moreover, many experiments involving Bose-Einstein condensates and cold atoms show a critical behavior even in the presence of a trapping potential [12,13].
In this paper we pursue further this program, studying the Ising model perturbed by a trapping potential coupled to the energy operator in 3 dimensions. There are many reasons to investigate this system. From a purely theoretical point of view, one can wonder if the TSS argument still holds if the trapping potential is coupled to the energy operator instead of the spin operator. Moreover, studying the effects of the energy-trapping potential on the 2-point functions out of criticality provides an alternative method to estimate the CFT data at the critical point.
Finally, the study of the correlation functions out of the critical point is relevant also from the experimental point of view. Indeed, a trapping potential coupled to the energy operator can be effectively seen as a perturbation of the system by a non-uniform thermal gradient, a thermal trap. This setup might be implemented in real system experiments and the knowledge of the correlators is fundamental in order to understand the behavior of the observables of this system out of the critical point.

II. THE MODEL AND THE TRAP SIZE SCALING
We consider the Ising model perturbed by a confining potential coupled to the energy operator: where S cf t is the d-dimensional Ising model action, z is the radial coordinate and U (r) = ρ|z| p is the trap potential. In this paper we will consider p ≥ 2, focusing mostly on the harmonic potential case p = 2. The parameter ρ is the trap parameter, which is related to the characteristic trap length l −p ≡ ρ defined in [14], and determines the shape of the trap.
Here we will study the large-trap case, namely the small ρ regime, where the CPT approach can be safely applied.
A. The one point functions As shown in [14], where the Trap Size Scaling (TSS) ansatz has been introduced, the behavior of the 1-point function can be extracted using renormalization group arguments.
In fact, it can be shown that near the critical point the one-point functions for spin and energy, defined in the center of the trap (z = 0), are where ∆ σ , ∆ are the scaling dimensions of the operators σ and respectively, A σ and A are non universal constants and the exponent θ is the characteristic trap exponent. This exponent can be determined using scaling arguments if one notices that the perturbation has to be scale invariant. Rescaling the radial coordinate z by a factor b, z → z b , the perturbation transforms as where ∆ ρ = 1 θ . Eventually the scale invariance condition yields

B. Two-point functions
Regarding the two-point functions, we can make use of the Operator Product Expansion (OPE) to express them as series involving 1-point expectation values: Each of the Wilson coefficient C k ij (ρ, z), evaluated outside the critical point, can be expanded in series of the trap characteristic parameter ρ, namely: As shown in [8], the series expansion asymptotically converges and all the coefficients are  (19)). Out of the critical point, the correlators can be expanded as a series of the parameter ρ in the following way: As said before, the derivatives of Wilson coefficient can be written in terms of known quantities [7][8][9]. For instance, ∂ ρ C I σσ (z 1 ) reads: This integral can be evaluated using a Mellin transform technique (see appendix A for details). In particular, the second term is just a regulator needed to cancel the IR-divergent part, meaning that only the first term contributes to the final result. Expanding the first term in (2.11) in terms of the known correlation function at the critical point we find: where ∆ t = 3−∆ and y = z 2 /z 1 . We refer to Appendix A for the details of the computation.
The final result is: where I(p) is numerical factor that can be expressed in terms of Gamma functions and the relevant parameters of the model, as shown in (A8). In what follow we are going to consider mostly the harmonic potential case, p = 2, for which I(2) −8.4448.
Following the same procedure we can also evaluate the derivative of C I : Putting all together, the expressions (2.8)-(2.10) can be expressed as: As usual in this approach, the asymptotic convergence of the series expansion is guaranteed for distances (measured from the center of the trap) less than about one correlation length. In the last equation the symbol # stands for the numerical value of the coefficient . The computation of this coefficient within the CPT framework involves the use of a 4-point correlation function at the critical point: is not known analytically at the critical point, (2.18) can not be evaluated exactly. However, as we will show later, this term can not be neglected and it will be determined a posteriori using Monte Carlo simulations.

III. CONVERSION TO THE LATTICE AND NUMERICAL RESULTS
The model previously described can be solved on a lattice in order to verify the validity of the CPT expansion and to get some insights in the numerical factors which we have not been able to determine analytically. The Hamiltonian of the system on a cubic lattice can be expressed in the following form: where σ i is the spin field, r i is its distance from the center of the confining potential and h is a possible magnetic field perturbation whose importance will be shortly explained. The conformal point is recovered for h = 0. To get a more precise physical intuition about the trapping effect, it is convenient to perform the transformation σ i = 1 − 2ρ i . Then, the new variable ρ i can only assume two values (0 and 1) and it can be thought as a density of particles in a d-dimensional gas. Eventually, the Hamiltonian reads: where µ = 2h − 4qJ is the chemical potential and q is the coordination number (q = 6 in three dimensions). The main advantage of this transformation is that, since the potential U (r i ) diverges at large r i , it makes it apparent that the only way to prevent the last term in (3.2) to diverge is to set either ρ i = 1 or ρ i = 0 for all i far from the center of the trap.
The first condition is not physically acceptable (all the particles running away to infinity) and it can be avoided by inserting a small and positive magnetic field h in eq. 3.1, namely: This leaves us only with the second possibility, which is equivalent to require a null density of particles ( ρ i = 0) far from the center of the lattice, which means that the system is trapped.
As the system breaks translational invariance, we may wonder G σ to be different from

B. One-point functions
Since the potential is coupled to the temperature, which in the lattice is non-zero at the critical point, the effective scaling parameter on the lattice to be compared with analytical prediction is ρ lat ≡ β c ρ. Thus, the one-point functions on the lattice are:  [14] also in the present case. Eventually, we fix the exponents to the value 3.7-3.8 and we repeat the fit with only two free parameters to find the remaining constants, obtaining A lat σ = 1.6390(13) and A lat = 2.226 (11).

C. Two-point functions
With the definitions 3.4-3.6 at hand, and taking into account the bulk contribution to the energy operator on the lattice E cr , the two-point functions on the lattice (denoted with with the average ... lat ) assume the following form: In order to make contact with the CPT theoretical results (2.8)-(2.10), we must consider the lattice conversion factors R σ and R . Regarding the first, R σ = 0.55245 (13) according to [21]. Estimates of R vary from 0.2306(38) [21] to 0.2377(9) [1]. This is the largest source of systematic uncertainty in our simulations. For this reason we will adopt the average R = 0.2341 with a variation ±0.0030 to evaluate the final systematic error. Finally, the structure constant on the lattice (C σσ ) lat must be converted taking into account the conversion rules for and ρ, namely lat = R and ρ lat = R −1 ρ. Eventually, combining (2.8)-(2.10) with (3.9)-(3.11) we obtain: (3.14) The parameter b in the second term of (3.14) is related to the coefficient (2.18), which, as already mentioned, we have not been able to compute analytical using CPT. This parameter will be evaluated a posteriori by fitting the numerical results.
We can now insert the lattice quantities A lat and A lat σ calculated in section III B, and directly fit the continuum structure constants C σσ and C . Fit results, reported in the table II, are in good agreement with the known values: C σσ = 1.0518537, C = 1.532435 [4].  one can see from the tables, once C σ σ is left as a free parameter its value agrees better with the known result if we take into account the parameter b. This is confirmed in Figure   3, where it is evident that the presence of b significantly improves the agreement between the theoretical prediction and the numerics. This proves that the second term in (3.14) is definitely important and must be taken into account.  14 and fixing C σ σ to the known value (third column), and fit of the structure constant C σ σ setting b = 0 (fifth column). It is evident from the data that b contributes non-trivially to the correlator, as our numerical results do not match the known value for C σ σ if we set b = 0.

IV. DISCUSSION
In this paper we have further developed the program of studying systems in their offcritical scaling regime, using the consolidated approach based on the OPE and the possibility   Fig. 1 confirm once again the validity of the TSS ansatz [14] as a powerful tool to determine the behavior of the expectation values of the model outside The integral 2.11 can be evaluated using a Mellin transform technique, following what has been done in [7][8][9]. In particular, it is convenient to introduce the quantity where Θ(|m|z) = e −m|z| is an IR-regulator needed to guarantee the convergence of the is essentially the Mellin transform of g up to angular coefficients. This means that in order to find an expression for the derivatives of the Wilson coefficients, one just needs to evaluate the Mellin transform of the function g(z).
In our case, as it can be seen from Eq. 2.12, The Mellin transform can be evaluated by performing the angular integral and rewriting the result in terms of beta-functions as follows: Then, we are ready to extract the m ∼ 0 behavior from A3. The only contribution comes from the residue at s = 0, so that which for p = 2 gives I(2) −8.4448.