Form Factors of the Tricritical Three-state Potts Model in its Scaling Limit

We compute the form factors of the order and disorder operators, together with those of the stress-energy tensor, of the two-dimensional three-state Potts model with vacancies along its thermal deformation of the critical point. At criticality the model is described by the non-diagonal partition function of the unitary minimal model $\mathcal{M}_{6,7}$ of conformal field theories and is accompanied by an internal $S_3$ symmetry. Its off-critical thermal deformation is an integrable massive theory which is still invariant under $S_3$. The presence of infinitely many conserved quantities, whose spin spectrum is related to the exceptional Lie algebra $E_6$, allows us to determine the analytic $S$-matrix, the exact mass spectrum and the matrix elements of local operators of this model in an exact non-perturbative way. We use the spectral representation series of the correlators and the fast convergence of these series to compute several universal ratios of the renormalization group.

We compute the form factors of the order and disorder operators, together with those of the stress-energy tensor, of the two-dimensional three-state Potts model with vacancies along its thermal deformation of the critical point.At criticality the model is described by the non-diagonal partition function of the unitary minimal model M6,7 of conformal field theories and is accompanied by an internal S3 symmetry.Its off-critical thermal deformation is an integrable massive theory which is still invariant under S3.The presence of infinitely many conserved quantities, whose spin spectrum is related to the exceptional Lie algebra E6, allows us to determine the analytic S-matrix, the exact mass spectrum and the matrix elements of local operators of this model in an exact non-perturbative way.We use the spectral representation series of the correlators and the fast convergence of these series to compute several universal ratios of the renormalization group.

I. INTRODUCTION
An interesting problem in theoretical physics consists of the study of two-dimensional statistical models which exhibit a second order phase transition for certain values of their coupling constants.As well known, the proper language to frame the scaling limit of these microscopical critical models is the one of conformal field theories (CFTs) [1] which, in two dimensions, have the peculiar feature of the underlying infinite-dimensional Virasoro algebra and the possibility to define consistent models employing only a finite number of families of irreducible representations of this algebra (the so-called degenerate fields).These 'minimal models' capture the universal scaling behavior of the two-dimensional critical models and are exactly solvable, in the sense that it is possible to obtain -at least in principle -closed expressions for the multi-point correlation functions of scaling operators, by exploiting properties of the Virasoro algebra and null-vector conditions of the corresponding degenerate fields [1,2].
However, no general framework has yet been formulated for generic off-critical models.Nonetheless powerful tools of analysis can be employed if the deformed theory turns out to be integrable, i.e. if it has infinitely many commuting conserved local quantities.When this happens, these conserved quantities provide severe constraints on the structure of scattering amplitudes, which are non zero only for the elastic processes and can be computed exactly [3][4][5][6][7][8][9].In practice it suffices to find all possible two-particle amplitudes, since all other scattering amplitudes involving n > 2 particles can be expressed in terms of their n(n − 1)/2 two-particle scatterings.
Once the S-matrix of an integrable theory is known, it is possible to study off-critical two-point correlation functions through the technique of form factors (FFs) [9,11,12].These quantities are defined as matrix elements of (semi-)local operators between the vacuum of the theory and multi-particle states.The novel contribution of this paper is the exact computation of FFs of the leading and sub-leading order and disorder operators in the thermal deformation of the TPM.We also compute some universal ratios of the renormalization group (see ref. [13] and references therein), which can be useful for numerical or experimental studies.
The q-state Potts model has been studied in great detail in the literature -see e.g.ref. [14] for a review.A generalization of this model, called 'dilute' due to the presence of vacancies, has been studied in [15].Its lattice realization has a Hamiltonian, invariant under the permutation group S 3 , of the form where each site i is associated with a Potts 'spin' s i = {0, 1, 2} and a vacancy variable t i = {0, 1}.In the latter equation K and F are interaction parameters associated to vacancies and Potts spins respectively, while D plays the role of chemical potential associated to vacancies.The usual Potts model is retrieved when the chemical potential D is negative and large in absolute value.Since the TPM possesses only two relevant S 3 invariant operators [2], the energy ϵ and the vacancy t, it is retrieved from eq. ( 1) in the limit of vanishing K.The TPM is also exactly solvable in a honeycomb lattice and proven to be self-dual, i.e. there is a mapping between its disordered (high-temperature) and ordered (low-temperature) phases [16].The content of this paper is organised as follows.In Section II we recall a few basic properties of minimal models of CFT, useful to set mostly the notation, followed by a thorough discussion of the TPM, in particular the operator content and the fusion rules of this model.Section III is dedicated to the thermal deformation of the TPM, regarded as an integrable field theory; the scattering theory and the mass spectrum are discussed there.The FFs are introduced in Section IV (their analytic properties are gathered in Appendix A).In Section V we review the FFs of the stress-energy tensor, originally obtained in [10].The FFs of the order and disorder operators of the thermal deformation of the TPM are presented in Section VI.The two-point correlation functions and the universal ratios of the renormalization group are discussed in Section VII.Some preliminary results coming from the Monte Carlo simulations of the model are presented in Section VIII.Our conclusions are finally summarised in Section IX.

II. CRITICAL TPM: THE CFT APPROACH
As well known, minimal models M p,q of CFT are characterized by the following properties [1,2]: i) their central charge c is always less than one and can be parameterized by two coprime integers p, q (with q > p)  ii) there is a finite number of families of degenerate fields, both in the analytic and anti-analytic sector of the CFT.
In the analytic sector, a degenerate primary field ϕ r,s (z) has conformal dimension with 1 ≤ r ≤ q − 1 and 1 ≤ s ≤ p − 1, coincident to that of the anti-analytic counterpart ϕ r,s (z): ∆ r,s = ∆ r,s .Unitary minimal models are retrieved by setting in the previous formulas m = q = p + 1 with m ≥ 3.
iii) correlation functions of all conformal fields satisfy a set of linear differential equations that can be exactly solved, yielding in turn the structure constants of the operator algebra.
Another important feature of minimal models is that the operator product expansion (OPE) algebra of degenerate fields, in the analytic and anti-analytic sector separately, is 'closed': the OPE of two degenerate fields can be written in terms of degenerate fields and their descendants.In particular, the operators appearing in the OPE are obtained by iterative multiplication of the fields ϕ 2,1 and ϕ 1,2 .The general formula, which does not give information about the structure constants, is Finding the modular-invariant partition functions of a given universality class associated to a minimal model is equivalent to determining the operator content of the model.The problem of the modular-invariant partition functions for unitary theories was first addressed by Cardy [18,19], and the solution for minimal models with periodic boundary conditions was found soon after by Cappelli et al. [17].
The conformal dimensions of the fields in M 6,7 are summarized in Table I.The universality class M 6,7 contains two algebras of fields which close under the OPE: the pentacritical φ 10 Landau-Ginzburg field theory and the TPM.In correspondence to the minimal model M 6,7 there are then two modular-invariant partition functions: a diagonal one, which corresponds to the φ 10 Landau-Ginzburg theory, and an off-diagonal one relative to the TPM Here χ ∆ is the character of the irreducible representation of the Virasoro algebra of weight ∆, whose explicit formulas can be found in [20].
Using the notation of Table II for the various fields, it is possible to show that the fusion rules of Table III are compatible with a S 3 internal symmetry, so that operators can be organized according to irreducible representations of this discrete group, i.e. either in doublets (s, Z and A) or in the one-dimensional neutral representation (1, ϵ, t, X, Y and B).In the model there are no operators which transform according to the alternating representation.Table II: Operator content of the TPM at criticality.
Table III: Fusion rules of the TPM at criticality.Doublet operators like s and Z are charged under the generator Ω of the subgroup Z 3 .Let the components of the doublet corresponding to the leading operator s be respectively σ and σ.Similarly one defines the components of the subleading order operator Z to be Z and Z.As it happens in the case of the Ising model [21], the tricritical Ising model [22] and the three-state Potts model [23,24], in view of the self-duality of the lattice model we can introduce a leading disorder operator (µ, μ) associated to the leading order operator (σ, σ).In complete analogy to the leading counterpart, (ζ, ζ) is the subleading disorder operator associated to the subleading order operator (Z, Z).The conformal dimensions of these operators are: The transformation properties under S 3 of the order and disorder operators are as follows.Under the generator Ω of Z 3 the disorder operators are invariant and the order operators acquire a phase: where ω = exp(2iπ/3).On the other hand, the generator C of Z 2 acts like a charge conjugation operator: The assumption that only order operators are charged with respect to Z 3 is consistent with an alternative definition of TPM, in terms of a parafermionic algebra of central charge 6/7 -see eq.(A.8) of ref. [25].The parafermionic theory is defined to be invariant under the product group Z 3 × Z3 .The first Z 3 group acts nontrivially only on order operators, while the second only on disorder ones.However, the symmetry group of the TPM is S 3 , which contains Z 3 and Z 2 as subgroups.While Z 3 is the group of symmetry of the parafermionic theory, the second group Z 2 acts as a charge conjugation operator.
There are two important features of parafermionic theories.The first is that they are self-dual by definition: order and disorder operators in the ordered phase can be mapped into disorder and order operators of the disordered phase respectively.As a consequence, we will only focus on the FFs in the disordered phase.The second is that order and disorder operators are not mutually local, in general.In particular, for the TPM the non-locality factors are Identical results apply if we charge-conjugate both the first and second entry in the γ's.

III. THERMAL DEFORMATION OF THE TPM
Consider the fixed point action A 0 of the TPM and define the action A of its thermal deformation as where ϵ is the energy operator and λ is its coupling constant.In the scaling limit of the lattice version of the model, i.e. the limit in which the temperature T approaches its critical value T c and the chemical potential is fixed at its critical value, we have λ ∝ (T − T c ).The action A corresponds to a massive theory and there is a non-vanishing trace of the stress-energy tensor Θ(z, z) given by In the following we will often denote the quantum field theory associated to the thermal deformation of the TPM as ϵ-TPM.In addition to the stress-energy tensor, in the ϵ-TPM there are infinitely many local conserved quantities which commute with each other and therefore the theory is integrable [3,7].The spins of conserved quantities for this model are given by the infinite sequence of residues {1, 4, 5, 7, 8, 11} mod h = 12 .
In addition, they coincide, modulo the Coxeter number h = 12, with the Coxeter exponents of the exceptional algebra E 6 .
The hidden E 6 structure of the model can be understood employing a different CFT description of the TPM, i.e. the so-called coset construction [26], which in this case consists in Let us discuss this point in more detail.The affine algebra (E 6 ) k (Fig. 1) has different integrable highest weights depending on the level k.Calling w i (i = 0, 1, . . ., 7) the fundamental weights of (E 6 ) k , which are the basis of the weight space in the natural basis w i = êi+1 (e.g.w 1 = (0, 1, 0, 0, 0, 0, 0)), the integrable highest weights at level k = 1 are w 0 , w 1 and w 5 , while at level k = 2 we have 2w 0 , 2w 1 , 2w 5 , (w 0 + w 1 ), (w 0 + w 5 ), (w 1 + w 5 ), w 2 and w 6 .Through the coset fields, we can recover the operator content of the TPM.As discussed for instance in [26], the procedure consists of finding the allowed branchings, the fractional dimensions of the integrable representations, conjugate representations (related to the symmetry group of the non-affine Dynkin diagram) and coset fields, i.e. classes of equivalent branchings related by automorphisms (Z 3 in the case of affine E 6 ).The results are summarized in Table IV.
The coset construction of the TPM only holds for the critical theory.Away from criticality there is however a correspondence between the thermal deformation of the TPM and the exceptional algebra E 6 and this involves the Toda field theories.These theories are built starting from a Lagrangian associated to a Lie algebra of finite rank (see, e.g.[9,27] and references therein).In [28] it was proven that the deformation of the Toda field theory of E 6 perturbed with the adjoint representation is integrable and has conserved charges whose spins coincide with those of the thermal deformation of the TPM.Thus the scattering matrix of this perturbed Toda field theory has the same minimal scattering matrix as our model: its mass spectrum contains six particles, two particle-antiparticle doublets, (l, l) and (h, h), and two self-conjugated particles, L and H, and the masses of these particles can be expressed as Cosets Operators where M (λ) is the mass of the lightest kink of the theory and depends on the coupling constant λ in eq. ( 11) as [29] M where In the latter equation γ(z) is a shorthand for Γ(z)/Γ(1 − z).
For the thermal deformation of the TPM, the excitations are ordinary particles only in the disordered phase λ > 0, while some of them become kinks in the ordered phase λ < 0. In a quantum field theory with localized interactions, particles are associated to asymptotic states which transform according to irreducible representations of the symmetry group of the theory.For the thermal deformation of the TPM, we have two S 3 -doublets and two singlets.Since the particles in the doublets (l, l) and (h, h) are generated by order (magnetization-like) operators, we define the nonlocality properties of the particles with respect to disorder operators to be Identical relations hold if l is replaced by h.This choice of non-local factors is consistent with the scattering theory of the ϵ-TPM, which will now be outlined.

A. Scattering Theory
In (1 + 1)-dimensional theories, all particles scatter along a straight line.It is useful to use the rapidity variable θ which parameterizes the on-shell momentum of the particle of type a and mass m a as p a (θ a ) = m a (cosh θ a , sinh θ a ) (18) and denote the one-particle state as |A a (θ)⟩.For multiparticle states we assume, as a convention, that all incoming particles can be written from left to right, with decreasing rapidity.On the other hand, outgoing states are written in the converse: particles with the highest rapidity on the right and particles with lowest rapidity on the left.Examples are provided in Table V.
Integrable models have commuting conserved quantities Q s that can be diagonalized simultaneously with the Hamiltonian.This condition imposes strong constraints on the scattering matrix.Apart from being unitary and satisfying crossing symmetry, the S-matrix is elastic and all n-particle processes can be factorized in terms of twoparticle scatterings [7,9].These two properties, along with the assumptions of bootstrap, i.e. treating bound states on the same footing as asymptotic particles, allow for the determination of the analytic two-body scattering amplitudes.In the case of the ϵ-TPM, the S-matrix [3,6] is diagonal, i.e. reflection amplitudes vanish, and the various amplitudes are reported in Table VI.They can be expressed in terms of the building block functions By looking at bound states with positive residue of two-particle scattering amplitudes it is possible to extract the mass spectrum of the model, which also coincides with the one of the E 6 Toda field theory eq. ( 14).
In the ϵ-TPM the degeneracy on the particle-antiparticle states (l, l) and (h, h) is lifted if we require that particle singlets and doublets transform according to the singlet and doublet irreducible representations of S 3 respectively.In particular, we suppose that particles l and h have Z 3 -charge ω, l and h have Z 3 -charge ω * , while L and H do not transform under the generator Ω of S 3 .On the other hand, Z 2 charge conjugation maps particles to their respective antiparticles.In the specific case of the two singlets L and H, particle and antiparticle coincide.
Table VI: S-matrix of the ϵ-TPM [3,6].Here [x] is a shorthand for s x/12 (θ) (see eq.( 19)).The superscript refers to the associated bound state.Notice that 12 is the Coxeter number of the exceptional algebra E 6 .

IV. GENERALITIES ON FORM FACTORS
The aim of this section is outlining the main properties of FFs, in particular the ones useful to address one and two-particle states.A further discussion of the analytic properties of such objects can be found in Appendix A.
For a relativistic theory, a generic matrix element of an Hermitian operator O between an incoming and outgoing state is given by FFs are nothing more than matrix elements between the vacuum of the theory and incoming multi-particle states, as pictorially described in Fig. 2. They are related to the previous equation by crossing, are analytic functions of the rapidity difference θ ij = θ i − θ j and are defined as For a scalar operator, knowledge of FFs allows us to write the spectral series for the correlation functions which, for the two-point correlators, reads [30] ⟨O An important feature of these spectral series is their rapid convergence [9,31].Practically, lowest mass particle FFs contribute more than higher mass particle ones, as it will be evidenced later by applications of the ∆-theorem [32] which is a sum rule for the conformal weight of an operator Φ  This sum rule is similar to Zamolodchikov's c-theorem sum rule for the UV central charge [33] For integrable theories, FFs of scalar operators can be parameterized as where D a1...an (θ 1 , . . ., θ n ) is a 2πi-periodic even function of the rapidity differences, Q O a1...an (θ 1 , . . ., θ n ) encodes the non-locality properties of O with respect to the operator generating particles, and F min aiaj (θ ij ) are the so-called 'minimal' FFs, which satisfy Being relevant to Section VI, the non-locality equations satisfied by FFs in the two-particle case read with θ = θ 1 − θ 2 .FFs can be found by solving a set of linear equations in correspondence of kinematic and dynamical poles, for a suitable parametrization of Q O a1...an (θ 1 , . . ., θ n ).Concerning bound state poles, e.g.
where Γ k ij is the real three-point vertex There are also residue equations coming from double poles of the S-matrix.The equation for a pole located at iϕ reads where γ = π − u ī cd − u j dē ( j denoting the antiparticle of j).The interested reader can find a thorugh discussion of these equations in Refs.[10,34].Kinematic poles are due to particle-antiparticle annihilation and in these cases two-particle FFs obey [35] −i lim where ⟨O⟩ is the vacuum expectation value (VEV) of the operator O. Table VII: Multi-particle states with lowest center of mass energy and nontrivial FFs with respect to the Z 3 -invariant operators.
In the two-particle case, Q O aiaj (θ) is a polynomial in e θ , which satisfies an upper bound on its maximal degree [34].Define y O to be the coefficient multiplying θ i as it approaches infinity, i.e.
The upper bound on the degree of the FF reads where ∆ O is the conformal weight of the operator O in the UV CFT.The latter bound and the aforementioned linear equations are in most cases sufficient to fix the FFs as functions of the VEV of some operator.However, it may happen that they are not enough.In computing the FFs of order and disorder, leading and subleading, operators of the thermally deformed TPM, we also made use of second order equations coming from the asymptotic factorization of FFs [32].

V. FORM FACTORS OF THE STRESS-ENERGY TENSOR
This section outlines the computation of the FFs of the stress-energy tensor Θ(x), as originally obtained in [10].Since the FFs of the stress-energy tensor corresponding to the same particle states have all the same pole structure, the two-particle FFs of the trace Θ can be parametrized by the following polynomial: where the degree is determined by the asymptotic behavior eq.(33).Notice that when particles have the same mass, i.e. in the case of a particle-antiparticle doublet or a singlet, the prefactor in the latter equation is completely absorbed in the unknown coefficients c (i) ab .Being conserved quantities of the theory, the FFs of Θ do not obey the usual residue equations for particle-antiparticle annihilation, rather: Apart from these properties, Θ is a Z 3 -invariant operator, like µ and ζ.The problem of finding the one and two-particle FFs is therefore similar for those three operators.In the following we'll present the parametrization and the equations satisfied by the one-particle and first two-particle FFs (Table VII).The first two-particle FFs of Θ are parametrized as follows Solutions of this system of equations are provided in Table VIII, and reproduce the ones of [10].As a further check of consistency, we computed the contributions to the c-theorem eq. ( 24) (Table IX).Since the fundamental mass m l obeys eq. ( 35), it is possible to compute ⟨Θ⟩ from the cluster equation [36] lim θ→∞ where w LL is a phase, which can be either 1 or −1.In order to make the ∆-theorem contributions for the energy operator ϵ positive, we choose the phase w LL such that ⟨Θ⟩ is negative, yielding Table X: ∆-theorem contributions for the energy operator ϵ of the first five particle states (∆ ϵ = 1/7 ≈ 0.14286).
Combining Eqs. ( 12) and ( 15) we get which coincides with the VEV of the energy operator in the thermal perturbation of the whole minimal model M 6,7 in the conformal normalization [37].
The FFs of ϵ retain the same pole structure and polynomial dependence of the ones of Θ, with modified coefficients; in particular, the FFs are parametrized as follows with the unknown coefficients satisfying For the sake of completeness, their values are summarized in Table VIII and the contributions of the respective FFs in the ∆ sum rule are shown in Table X.

VI. FORM FACTORS OF ORDER AND DISORDER LEADING AND SUBLEADING OPERATORS
This section is dedicated to the derivation of the first FFs of magnetic operators in the thermally deformed TPM.Let us start from the leading disorder operator.

A. Leading Disorder Operator µ
Because of C-invariance, ⟨μ⟩ and ⟨µ⟩ coincide.In Table VII are presented the first states with lowest mass and nontrivial FFs.Since the computation of these FFs involves non-locality factors, the derivation will be outlined step-by-step.
Let us start from F µ l l (θ).Since in the end we are concerned with computing the off-critical two-point correlation functions of the operators, it is important to notice that the only contribution that needs to be evaluated numerically in eq. ( 22) is that of F µ l l (θ).Indeed, because of C-invariance, the FFs of the operator μ are read from that of µ: where the last equality is a consequence of eq. ( 26).
Since the disorder operator has nonzero non-locality factor with respect to particles l and l, the FF will also possess an additional term which takes into account particle-antiparticle annihilation.Because of the asymptotic behaviors of F min l l (θ) and D l l(θ), Q µ l l(θ) satisfies the following inequality on its degree: Notice the absence of the absolute value on the first term in the l.h.s..The polynomial Q µ l l(θ) is not symmetric in the exchange θ → −θ, since the non-locality factors of µ with respect to particles l and l are nonzero.The non-locality equation that must be satisfied is with γl ,µ = −1/3.Solutions are sought of the form exp(αθ), provided that they satisfy e αθ = −e −i2π/3 e α(θ+2πi) .
Therefore the parametrization for this FFs is Pole equations allow us fixing all three unknowns -a 0 l l, a 1 l l and F µ L -as functions of the VEV of the leading disorder operator µ (Table XI).
As far as F µ LL (θ) is concerned, the absence of non-locality phases simplifies the parametrization of the FF: However, the corresponding S-matrix has a second-order pole; therefore the set of pole equations is Again, all three unknowns (F µ H , a 0 LL and a 1 LL ) can be fixed by simply solving these linear equations (Table XI).The first nontrivial check of the consistency of these solutions is given by the second order cluster equation satisfied when the phase w LL = 1.
The FF F µ l h(θ) is obtained along the lines of F µ l l (θ).In particular, its parametrization reads and satisfies the pole equations   Notice that the latter two pole equations have a peculiar dependence on the two-particle FFs F µ ll (θ) and F µ l l (θ).Solutions are reported in Table XI.We also computed the first contributions to the ∆-theorem (Table XII).

B. Subleading Disorder Operator ζ
The computation of FFs for the subleading disorder operator ζ presents little differences if compared to the case of µ.First of all, since the ζ has the same transformation properties as µ, the non-vanishing FFs of ζ are the same as the ones of µ (Table VII).The only difference in computing the FFs between the leading and subleading disorder operator lies in their different anomalous dimensions, which in turn affects the upper bound on the FFs degree.On the other hand non-locality coefficients remain unchanged in passing from the leading to the subleading case.These considerations are accounted for in the parametrization of some of the Q polynomials:   While the second solution is associated to the disorder operator µ, the first one is acceptable: its contribution to the ∆-theorem is given in Table XII and it is compatible with the conformal weight of the operator.Notice that convergence is slower than the leading case, as it generally happens when one considers less relevant operators -see also Ref. [22].Solving the system of equations, together with eq. ( 63), yields the results of Table XI.
C. Leading Order Operator σ and Subleading Order Operator Z The FFs of the leading and subleading order operators are simpler to compute if compared to their disorder counterparts, because there are no non-locality factors involved.The first nontrivial FFs are summarized in Table XIII.It turns out that the solution of their pole equations require the knowledge of the FFs of disorder operators.
The parametrizations of the first FFs of σ are Table XIV: Parameters of the first FFs of the leading order operator σ and subleading order operator Z.
and the set of pole equations is These equations are linearly dependent and F σ l cannot be fixed in this way.This value, however, can be obtained employing the cluster equation This equation allows us to fix F σ l up to a sign, which however is not relevant to the computation of two-point correlation functions.The parameters are summarized in Table XIV.
For the subleading order operator, the parametrization of the FFs is and the set of pole equations is the same as the one of σ, provided that we substitute σ → Z and charge-conjugate all particles, since σ and Z have complex conjugate Z 3 phases.It turns out that, in solving the set of linear equations, F Z l and F Z h remain unknown.Conjointly solving allows us finding F Z l and F Z h without any phase ambiguity between them -F Z l is fixed up to an overall sign, still.The parameters of the FFs of Z are summarized in Table XIV.

VII. UNIVERSAL RATIOS OF THE RENORMALIZATION GROUP
In the following we are going to consider some universal ratios of the renormalization group associated to the thermal deformation of TPM, i.e. to the action eq.(11).Consider an operator ϕ j of the theory which has at criticality the anomalous dimension X j = 2∆ j and therefore renormalization group eigenvalue y j = 2 − 2∆ j .The coupling λ of the theory is a dimensional quantity related to the lightest excitation of the theory m l as in eq.(15).Denoting by K j the metric factor associated to the operator ϕ j , on general ground the scaling form of the VEV ⟨ϕ j ⟩ can be written as (see [13] for a more detailed discussion) with , are related to two-point connected off-critical correlators through the fluctuationdissipation theorem Universal ratios are defined as ratios involving the correlation length ξ, the expectation values ⟨ϕ j ⟩ and generalized susceptibilities Γϵ jk such that the dependence on the metric factors K j cancels out [13].These quantities are associated to the class of universality and therefore are the same for different microscopic realizations of the model.The first trivial universal ratio is the ratio between the correlation lengths in the disordered and ordered phase: In order to find nontrivial universal ratios starting from the FFs computed in previous sections, eq. ( 81) must be rewritten as to make the dependence on λ manifest.Notice however that some susceptibilities are identically zero, because of Z In general, FFs depend naturally on the VEV of some operator F ϕj a1,...,an (θ 1 , . . ., θ n ) = f ϕj a1,...,an (θ 1 , . . ., θ n )⟨ϕ where ϕ ′ j is often coincident with ϕ j , like in the case of µ, ϵ and ζ, but not always, as it happens for σ and Z. Manipulations of eq. ( 22) yield where Here β aj = m aj /m l .By approximating the two-point correlation function with the sum of the contributions associated to the one and two-particle states with lowest center-of-mass energy, we present hereafter some universal ratios of the Ref. [40] where they were estimated by means of a different method, i.e. the transfer-matrix technique.We performed simulations at 81 equally spaced values of K in the [K tc − 0.04, K tc + 0.04] interval; note that values of K larger than K tc correspond to thermal deformations driving the system into a broken-symmetry phase, in which the magnetization is expected to be non-zero in the thermodynamic limit.The square lattices that we simulated had sizes (L/a) 2 = 64 2 , 128 2 , 256 2 , and 512 2 .For each value of K and each lattice size, after at least 10 4 complete updates for reaching Monte Carlo thermalization, we collected 8 • 10 5 configurations for the lattices of sizes (L/a) 2 = 64 2 and 128 2 , 4 • 10 5 configurations for the lattices of size (L/a) 2 = 256 2 , and a number of configurations between 387610 and 751375 for the lattices of size (L/a) 2 = 512 2 .We used these ensembles of configurations to compute the expectation values of all observables described below, with the exception of the spin-spin correlation function, which was estimated from ensembles consisting of 4 • 10 5 configurations for each value of K on the lattices of sizes (L/a) 2 = 64 2 , 128 2 , and (L/a) 2 = 256 2 , and a number of configurations between 201835 and 386155 for each K value on the lattices of size (L/a) 2 = 512 2 .
The observables that we evaluated in our Monte Carlo simulations include: • The density of vacancies, which is the fraction of sites in the s i = 0 state • The contribution to the energy (in units of k B T ) coming from spin-spin interactions, divided by the number of lattice sites, which is denoted as In the following (with a slight abuse of terminology), we will also refer to this quantity simply as 'energy'; in addition, we considered the associated susceptibility, too.
• The magnetization in each configuration, which is defined as where f max denotes the fraction of spins in the majority state in the given configuration, • Finally, the two-point spin-spin correlation function, which is defined as where û0 and û1 respectively denote the positively oriented unit vectors in the t and x directions of the lattice.
For all of these variables, we computed the integrated autocorrelation time (for a discussion about the numerical estimate of this quantity in Monte Carlo simulations, see, e.g., Ref. [41]), and grouped the whole set of 'measurements' of each observable obtained in every simulation into bins of length equal to the integrated autocorrelation time.We verified that, as expected, the average over each bin produced data with negligible residual autocorrelation.Figure 3 shows how the distribution of vacancies extracted from our simulations varies as a function of K tc − K.The first observation is that the results from simulations on lattices of the different sizes that we considered exhibit a remarkable collapse of data, indicating that finite-volume corrections are under control.Furthermore, the data are consistent with an inflection point at K = K tc , and the value of ρ at that point is in full agreement with the result (denoted by the black cross) obtained through the transfer-matrix technique in Ref. [40].
Next, we consider the average value of ϵ, which is plotted in Fig. 4, together with the results for the associated susceptibility (displayed in the inset plot, with a logarithmic scale on the vertical axis).As for the vacancy density, also for ⟨ϵ⟩ we observe a clear collapse of the data obtained from simulations on lattices of different area, and an  inflection point at a location consistent with K = K tc .For this observable, the Monte Carlo results are expected to be modelled by a function of the form: Fitting the data in the range 0.005 ≤ K tc − K ≤ 0.04 to eq. ( 101), one obtains e 0 = 2.036( 4) and e 1 = −1.6124(23).
The curve obtained from the fit is shown by the black line in Fig. 4. The good quality of this fit, in which the exponent of (K tc − K) is fixed to the value predicted analytically, confirms the validity of the conformal field theory description of the results obtained from Monte Carlo simulations, at least for the range of K values and for the lattice sizes under consideration in this study.101) is also shown.The inset plot shows the results for the susceptibility associated with ϵ (rescaled by the number of lattice sites), which are displayed with a logarithmic scale for the vertical axis.
Our Monte Carlo results for the average magnetization are displayed in Fig. 5; note that, in the broken-symmetry phase, the results obtained from simulations on lattices of different area are consistent with each other.On the other hand, in the symmetric phase, the results obtained from smaller lattices are systematically larger than those from larger lattices: this is an artefact of the definition of the magnetization in eq. ( 98), which is based on the fraction of spins in the majority state in each configuration: by definition, for a system with a finite number of degrees of freedom, in the symmetric phase f max is always larger than 1/3 -except in the exceedingly rare configurations in which the total number of variables in the spin states 1, 2, and 3 are exactly equal to each other.Nevertheless, our Monte Carlo results show that the magnetization in the symmetric phase vanishes when it is extrapolated to the thermodynamic limit.
In the broken-symmetry phase, fitting the results for the magnetization in the range −0.04 ≤ K tc − K ≤ −0.005 to one obtains d 0 = 0.9526 (15) and d 1 = 0.1185 (13).The result of this fit is shown by the black line in Fig. 5.As for ⟨ϵ⟩, the fact that one obtains a good modelling of Monte Carlo results using this simple functional form, with the exponent of |K tc − K| fixed to the value that is predicted by conformal field theory, confirms the validity of this description of Monte Carlo results.
-0.04 -0.02 0 0.02 0.04 Figure 5: Average magnetization, plotted against K tc − K, from our Monte Carlo simulations on lattices with the same linear sizes as in Fig. 3 and in Fig. 4. The curve obtained fitting our data in the broken-symmetry phase to eq. ( 102) is also shown (black line).
Next, we show examples of our results for the two-point spin-spin correlation function in Fig. 6.In particular, it is interesting to study the universal ratio between the correlation lengths in the symmetric and in the symmetrybreaking phases.Sufficiently close to the tricritical point, it is expected that the ratio of correlation lengths obtained at values of K that are symmetric with respect to K tc should become consistent with the theoretical prediction ξ + /ξ − = √ 2 = 1.41421....As an example to illustrate this type of study, the data displayed in Fig. 6 are obtained from two test ensembles of Monte Carlo simulations on a lattice of size (L/a)2 = 256 2 at K = K tc ± 0.01, i.e., close to the tricritical point. 1xpressing all distances in units of the lattice spacing, the correlators can be fitted to It should be noted that the extraction of the correlation length from the Monte Carlo results for the correlators is non-trivial, since it is based on a three-parameter, non-linear fit of cross-correlated data, and the fitted parameters are quite strongly correlated with each other, too.This, in turn, leads to a likely overestimate of the uncertainties affecting the fitted parameters, including, in particular, the correlation length.However, we followed a very conservative approach to the error budget, without attempting to reduce these uncertainties.Figure 7 shows the correlation lengths computed from Monte Carlo simulations on lattices of different sizes and at different K values.In particular, the main plot (in which a logarithmic scale is used for both axes) displays the data in the symmetric phase and their fit to the power-law dependence on |K tc − K| predicted from the theory: ξ = a 0 |K tc − K| −7/12 . (105) The plot makes manifest that the simulation results precisely follow the power-law behavior predicted by conformal field theory.The only free parameter of the fit to eq. ( 105) is the overall amplitude, for which we obtain a 0 = 0.356 (1).Analogous results in the broken-symmetry phase are shown in the inset plot (where a linear scale is used for both axes).

IX. DISCUSSION AND CONCLUSIONS
In this work we have addressed the computation of the form factors of the relevant operators of the TPM once this model is perturbed away from its critical point by varying the temperature.The corresponding massive theory is integrable and has a distinguished fingerprint associated to the E 6 exceptional Lie algebra.We have determined the FFs associated to the leading and subleading order and disorder operators and we have checked their correct identification using the ∆ sum rule.We have shown how to theoretically compute several universal ratios of the renormalization group associated to the class of universality of the model in any of its microscopic realizations.We also presented some results coming from a Monte Carlo study of the lattice model, which provided clear evidence that, in the range of parameters that we considered, the results are in agreement with the expectations from CFT, and that, in addition, the universal ratio of the masses in the high-and low-temperature phases of the model is fully consistent with our theoretical prediction.The main plot refers to simulations at K < K tc , whereas the inset plot shows analogous results in the broken-symmetry phase.

Figure 1 :
Figure 1: Dynkin diagram of the affine E 6 algebra.The first term in parenthesis labels the root, while the second its comark.

Figure 2 :
Figure 2: Form factor of the operator O.

Figure 3 :
Figure3: Average density of vacancies, displayed as a function of K tc − K, obtained from our simulations on square lattices of different sizes (denoted by symbols of different shapes and colors).The plot also shows the average density of vacancies at K tc (shown by the black cross) that was computed in Ref.[40] using a different numerical technique.

Figure 4 :
Figure4: Average contribution to the energy density from spin-spin interactions, plotted against K tc − K, from Monte Carlo calculations on square lattices of linear size ranging from L/a = 64 to L/a = 512; the curve obtained by fitting the data for 0.005 ≤ K tc − K ≤ 0.04 to eq. (101) is also shown.The inset plot shows the results for the susceptibility associated with ϵ (rescaled by the number of lattice sites), which are displayed with a logarithmic scale for the vertical axis.

GFigure 6 :
Figure 6: Two-point spin correlation functions computed from Monte Carlo simulations on a lattice of size (L/a) 2 = 256 2 at two different K values, shown as a function of the spatial separation r in units of the lattice spacing a.The dashed lines are the curves obtained from the fits to eq. (103).

Figure 7 :
Figure7: Correlation lengths computed from Monte Carlo simulations at different K values and different lattice sizes, and their fit to the power-law behavior predicted from conformal field theory (black line), given by eq.(105).The main plot refers to simulations at K < K tc , whereas the inset plot shows analogous results in the broken-symmetry phase.

Table I :
Kač table of the universality class M 6,7 .

Table VIII :
Parameters of the first FFs of the trace of the stress-energy tensor Θ and the energy operator ϵ.

Table XI :
Parameters of the first FFs of the leading disorder operator µ and subleading disorder operator ζ.

Table XIII :
[36]i-particle states with lowest center of mass energy and nontrivial FF with respect to the order operator σ.The ones of the subleading order operator Z are retrieved by charge-conjugating all particles.The set of linear equations for ζ remains the same as the ones for µ, except for a change in label µ → ζ in all of them.However, it turns out that there is one linearly dependent equation, which doesn't allow fixing one of the unknowns, e.g.F ζ L .Studying the cluster equation[36]

Table XV :
3-symmetry as, for instance, Γϵ σµ , Γϵ σZ and Γϵ σ σ .Charge conjugation also reduces the number of quantities First contributions for the generalized susceptibilities of Z 3 -invariant operators.

Table XVI :
First contributions for the generalized susceptibilities of magnetization operators.to compute, e.g.Γϵ μμ = Γϵ µµ .In Table XV are presented the generalized susceptibilities involving Z 3 -invariant operators (µ, ϵ and ζ), while Table XVI is dedicated to magnetization ones (σ and Z).
on High-Performance Computing, Big Data and Quantum Computing (ICSC) funded by MUR (M4C2-19) -Next Generation EU (NGEU), by the Italian PRIN "Progetti di Ricerca di Rilevante Interesse Nazionale -Bando 2022", prot.2022TJFCYB, and by the "Simons Collaboration on Confinement and QCD Strings" funded by the Simons Foundation.The simulations were run on CINECA computers.We acknowledge support from the SFT Scientific Initiative of the Italian Nuclear Physics Institute (INFN).