Entanglement Entropy of Non-Unitary Integrable Quantum Field Theory

In this paper we study the simplest massive 1+1 dimensional integrable quantum field theory which can be described as a perturbation of a non-unitary minimal conformal field theory: the Lee-Yang model. We are particularly interested in the features of the bi-partite entanglement entropy for this model and on building blocks thereof, namely twist field form factors. Non-unitarity selects out a new type of twist field as the operator whose two-point function (appropriately normalized) yields the entanglement entropy. We compute this two-point function both from a form factor expansion and by means of perturbed conformal field theory. We find good agreement with CFT predictions put forward in a recent work involving the present authors. In particular, our results are consistent with a scaling of the entanglement entropy given by $\frac{c_{\text{eff}}}{3}\log \ell$ where $c_{\text{eff}}$ is the effective central charge of the theory (a positive number related to the central charge) and $\ell$ is the size of the region. Furthermore the form factor expansion of twist fields allows us to explore the large region limit of the entanglement entropy and find the next-to-leading order correction to saturation. We find that this correction is very different from its counterpart in unitary models. Whereas in the latter case, it had a form depending only on few parameters of the model (the particle spectrum), it appears to be much more model-dependent for non-unitary models.

At quantum critical points, the scaling limit of the EE has been widely studied within unitary models of CFT [39,40,7,8,41,42]. In particular, the combination of a geometric description, Riemann uniformization techniques and standard expressions for CFT partition functions is very fruitful. Recently [43], this was generalized to non-unitary CFT, where a general formula was obtained using such techniques. Near critical points, the scaling limit is instead described by massive quantum field theory (QFT), and geometric techniques relying on conformal mappings break down. As was found in [44,45,46], the most powerful way of studying the EE in unitary models of QFT is using an approach based on local branch-point twist fields. However, the question of the EE in non-unitary near-critical models is much more delicate, and standard arguments give little indications as to how to modify the field-theoretical approach. Importantly, the rigorous derivation presented in [43] provided a precise local-field description of the EE involving composite fields in the branch-point twist family, thus opening the door to its study in non-unitary QFT. In the present paper, using techniques of integrable QFT, we will study the scaling limit of the EE in the near-critical region of von Gehlen's model, described by the Lee-Yang QFT model. This paper is organized as follows. In section 2 we recall the main definitions and techniques, and provide a summary of our main results. In section 3 we introduce the Lee-Yang model and some general results on the form factor expansion of correlation functions, their logarithms and expectation values of local fields. In section 4 we review the twist field form factor equations and present solutions for the branch-point twist fields fields T and : T φ : in the Lee-Yang model. In section 5 we test our form factor solutions by performing a form factor expansion of the functions log : T φ : −2 : T φ : (r) :T φ : (0) and log T −2 T (r)T (0) and recovering the behaviours −4x :T φ: log(mr) and −4x T log(mr) for some constants x T , x :T φ: which we compare to CFT predictions. In section 6 we compare a form factor computation of the two-point functions above with a computation in zeroth order conformal perturbation theory. As a byproduct, we find general formulae for some of the CFT structure constants entering the OPEs of T andT and of : T φ : with :T φ :. In section 7 we present numerical results for the Rényi entropy near criticality and a detailed computation of the first three leading corrections to saturation of the EE. We find that the next-to-leading order correction to saturation is non-universal. In section 8 we present our conclusions and outlook. In appendix A we explain how the normalization and conformal dimension of the field : T φ : are fixed by CFT. In appendix B we present a detailed analysis of the one-particle form factor contribution to the two-point functions of T and : T φ :. A large n expansion of this function demonstrates that it provides a very substantial contribution to the power law behaviour of the two-point functions at short distances. In appendix C we present a computation of the three particle form factor of Lee-Yang twist fields. In appendix D we perform a computation of some of the structure constants entering the OPE of fields T andT and of fields : T φ : and :T φ : in CFT. In appendix E we present a computation of the numerical coefficient of the next-to-leading order correction to saturation of the EE in the Lee-Yang model.

General aspects and summary of main results
In order to provide a formal definition of the EE, the Hilbert space of an extended quantum system, such as a spin chain, is decomposed into a tensor product of local Hilbert spaces associated to its sites. Grouping together sites associated to the regions A andĀ, this gives: The EE in a state |ψ is the von Neumann entropy of the reduced density matrix ρ A associated to A: Another frequently used measure of entanglement is the Rényi entropy, which specializes to the von Neumann entropy at n = 1, We will study the ground state entanglement entropy in the scaling limit of infinite-length quantum chains. The scaling limit gives the universal part of the quantum chain behaviour near quantum critical points, described by 1+1-dimensional QFT. It is obtained by approaching the critical point while letting the length of the region A go to infinity in a fixed proportion with the correlation length ξ (measured in number of lattice sites). If ξ = ∞ from the start, the system is exactly at its critical point, and the scaling limit is described by CFT. In this case the entanglement entropy of unitary critical systems, as a function of , is divergent in a way which was first understood in [39,40], numerically confirmed in [7,8] and generalized and reinterpreted in [41,42]. The divergency is logarithmic with a proportionality constant depending on the central charge c of the CFT, and where ε is a non-universal ultraviolet cut-off (proportional to the lattice spacing) which is chosen so as to encode all o(1) corrections. The formulae above are easily adapted to the case of an infinite region = ∞ near criticality ξ < ∞, where is simply replaced by ξ in (5) [41,42].
In the full scaling limit, where and ξ are both large and in proportion to each other, there is a universal scaling function f ( /ξ) which interpolates between the two results, In this case the result is much less trivial and has been studied in unitary integrable [44,45] and non-integrable [46] models using massive QFT techniques.
We may ask how (if at all) the entanglement entropy is affected by non-unitarity. At criticality, it was shown in [43] that the entanglement entropy scales instead as S (n) where c eff := c − 24∆, and ∆ is the smallest (often negative in non-unitary models) scaling dimension of a primary field in the CFT. For the Lee-Yang model, for example, c eff = 4 5 as ∆ = − 1 5 . This result is not entirely surprising as the work of Itzykson, Saleur and Zuber [47] had previously shown that the effective central charge c eff also replaces c in the expression of the ground state free energy found by Affleck [48] and Blöte, Cardy and Nightingale [49]. However, the question of the entanglement entropy in non-unitary near-critical models is much more delicate. Importantly, the rigorous derivation of (7) presented in [43] has lead to new insights into the computation of entanglement entropy in non-unitary theories and its field theoretical interpretation, opening the door to its study away from criticality in QFT.
It is known since some time [39,40,41,42] that the bi-partite entanglement entropy in the scaling limit can be re-written in terms of more geometric quantities, using a method known as the "replica trick". The essence of the method is to "replace" the original QFT model by a new model consisting of n copies (replicas) of the original one. These are used to represent ρ n A when n is an integer, and then to evaluate Tr A ρ n A . The quantities S A and S (n) A for general n are then obtained by "analytic continuation" in n. The matrix multiplications in ρ n A and the trace operation give rise to the condition that the copies be connected cyclically through a finite cut on the region A. As a consequence, this trace is proportional to the partition function Z n (x 1 , x 2 ) of the original (euclidean) QFT model on a Riemann surface M n,x 1 ,x 2 with two branch points, at the points x 1 and x 2 in R 2 , and n sheets cyclically connected. The positions x 1 and x 2 of the branch points are dimensionful positions in the QFT model corresponding to the end-points of the region A in the scaling limit. This gives: Here, r := |x 1 − x 2 | is the euclidean distance between x 1 and x 2 . The above concepts hold, in principle, for any QFT model, unitary or not. In CFT, one may evaluate this by using the uniformization theorem: the Riemann surface M n,x 1 ,x 2 can be conformally mapped to the Riemann sphere with two punctures (or the cylinder) by using the map g reproduced in appendix A.
In the EE context, it was first noticed in [41,42] that the ratio of partition functions above can be reinterpreted as correlation functions of certain fields, which were not otherwise specified, in unitary CFT. This idea was then generalized to unitary massive theories in [44] and the fields where identified as branch-point twist fields T (x 1 ),T (x 2 ) characterized by their non-trivial exchange relations with other fields of the n-copy theory. These twist fields are defined only in the replica model (e.g. they become the identity field when n = 1), and are primary fields arising from the extra permutation symmetry present in the replica theory; they are associated to the Z n symmetry generators j → j + 1 mod n and j → j − 1 mod n respectively. In CFT, such twist fields and their relation to partition functions on Riemann surfaces were in fact studied much before their use in the computation of the EE was emphasized, see for instance [50]. In terms of these fields, the replica partition function is given by where T (x 1 )T (x 2 ) is a two-point function in the ground state of the replica theory. The branchpoint twist fields are chosen so as to have the CFT normalisation (e.g. the leading term in their OPE has coefficient 1). The constant Z n , with Z 1 = 1, is an n-dependent non-universal constant, ε is a short-distance cut-off which is scaled in such a way that dZ n /dn = 0 at n = 1, and, finally, ∆ T is the conformal dimension of the counter parts of the fields T ,T in the underlying n-copy conformal field theory, which can be obtained by CFT arguments [50,42,44]. It is easy to show that the formula (9) when inserted in (4) indeed reproduces (5) for CFT. The derivation above assumes unitarity of the theories under consideration. In such case ∆ T is by construction the lowest conformal dimension of any field in the replica theory which has the twist property. The CFT derivation of (5) has been generalized to the non-unitary case in [43] leading to the expressions (7). In this work it was also observed that the EE could be computed from a representation of the Z n -orbifold partition function of the theory via correlation functions involving certain Z n twist fields of the n-copy replica theory. This representation requires new twist fields : T φ : and :T φ :, obtained from the primary twist fields T andT as leading descendants in the product with the lowest-dimension field φ (of conformal dimension ∆). More precisely: and similarly for :T φ :. These composite fields were first introduced in [51] and further studied in [52]. The constant n 2∆−1 ensures conformal normalization, namely in CFT, where are the conformal dimensions of : T φ : and :T φ :. In the context of the study of the EE they were first obtained in [51]. However, as for many other quantities in this context, they had emerged previously in the study of orbifold CFT, see e.g. [53,54,55]. A detailed derivation of both the normalization constant and the power law in (11) is given in appendix A. The dimension ∆ :T φ: arises naturally in computations of the entanglement entropy in non-unitary CFT, and, as noticed in [43], suggests that for such theories, the partition function of the n-copy theory may be written instead as: where again Z n is such that it and its derivative at n = 1 are 1 and ε is a short distance cut-off. Compared to (9) the expression (14) involves not only a different twist field but also a normalization given by φ(x 1 )φ(x 2 ) n . For CFT it is easy to interpret this normalization as simply the norm of the ground state which in radial quantization is created by the action of the field φ on the conformal vacuum. As for the unitary case, it is easy to show that the formula (14) when inserted in (4) reproduces (7) for CFT. It is natural to assume that the same expression will hold beyond criticality. This paper is a first step towards putting this assumption to the test beyond criticality.

Summary of main results
From the formulae above it is clear that a study of the EE in massive QFT is in principle only possible by studying correlation functions of twist fields. This approach has been pursued successfully in several works [44,56,45] where the ratio of partition functions (9) at large distances r = |x 1 − x 2 | (the infrared (IR) region) has been studied for unitary 1+1-dimensional integrable QFTs. Integrability means that in these models there is no particle production in any scattering process and that the scattering (S) matrix factorizes into products of two-particle S-matrices which can be calculated exactly (for reviews see e.g. [57,58,59,60,61]). Although most of the integrable theories studied in this framework are unitary, well-known examples of non-unitary integrable QFTs exist. Best known among those examples is the Lee-Yang model whose exact S-matrix was first given in [62].
Taking the known S-matrix of an integrable model as input it is possible to compute the matrix elements of local operators (also called form factors). This is done by solving a set of consistency equations [63,64], also known as the form factor bootstrap program for integrable QFTs. It is this particular feature of integrable models which makes them interesting, as it means there is a systematic, non-perturbative way of computing multi-point functions of local fields. These computations are by no means easy, but often provide good numerical accuracy. In [44], the form factor program was generalised to branch point twist fields leading to the evaluation of (9) for various unitary models. In this paper we will pursue this program for the Lee-Yang model employing the formula (14). Our main results can be summarized as follows: 1) Form factor program for twist fields: We have found that the twist field form factor equations together with the requirement of form factor clustering are sufficient conditions to entirely fix all form factor solutions for any particle numbers. In particular, these constraints immediately give rise to two form factor families, naturally identifiable with the fields T and : T φ :. We have carried out a zeroth order perturbed CFT computation of the twist field two-point function for several values of n and compared this to a truncated form factor expansion of the same correlator. The former is expected to be accurate at short distances, the latter at large distances. Nevertheless, the agreement is relatively good, thus confirming the validity of the form factors found.
2) Saturation of the EE at large subsystem size: Let us absorb all non-universal o(1) constants of the short-distance behaviour of the Rényi EE into a short-distance cutoff n . Subtracting this non-universal contribution, the EE at large distances then saturates to a universal constant which can be calculated using QFT. More precisely, we find where the universal saturation U n is given by The constants K O are fundamental properties of QFT fields O, defined by where x O is the unique exponent making the limit finite and nonzero, and where O is the "conjugate" under internal symmetries ( φ = φ and : T φ : =:T φ :). In the unitary case, In the non-unitary case, both of these statements are modified. In particular, in the Lee-Yang model, these constants can be expressed in terms of massive QFT and of CFT data as where the vacuum expectation values are under CFT normalization, andC φ φφ andC φ 1 ···φn :T φ::T φ: are structure constants of conformal OPEs. Constants K O can also be expressed solely in terms of form factors of the massive model, as in (32).
3) Leading order correction to saturation: For unitary theories, one of the most interesting results [44,56] has been the identification of a universal leading order correction to the large-distance (large-r) saturation of the entropy of all unitary integrable theories. This exponentially decaying correction, of order o(e −2mr ), has a higher degree of universality than usual QFT quantities, as it only depends on the particle spectrum of the model. In [46] it was shown that, even more strikingly, this feature holds beyond integrability. This, however, seems to be broken in non-unitary models. For the Lee-Yang model, we found where a = −0.0769782... is a constant that is (a priori) model-dependent. This suggests that we may use this feature of the EE as a means to identify non-unitary critical points. Indeed, given a spin chain model whose critical point is not known, a study of entanglement at criticality will reveal the value of c eff . However, this does not say if c eff = c or not. Considering large size corrections away from criticality will reveal different types of exponential decay depending on whether or not the theory is unitary.
3 S-matrix and form factors in the Lee-Yang model

S-matrix
The Lee-Yang model is one of the simplest 1+1 dimensional integrable QFTs. From the CFT point of view, it may be regarded as a perturbation of the non-unitary minimal model associated with central charge c = − 22 5 . The primary operator content of the theory is very simple, consisting of the identity and a scalar field φ of conformal dimension ∆ = − 1 5 . Perturbing this CFT by the scalar field we obtain the massive Lee-Yang model. This theory has a single particle spectrum. The scattering amplitude corresponding to the scattering of two particles of the same type was found by Cardy and Mussardo [62] and can be written as It has a pole in the physical sheet at θ = 2πi 3 corresponding to the formation of a bound state, which in this case is the same fundamental particle of the theory. We note that the nonunitarity is manifested by the fact that the associated residue has the wrong sign. Nevertheless, the corresponding integrable massive model is well defined. The n-copy model, where the Z n twist fields live, possesses n particle species µ = 1, . . . , n, and a two-particle scattering matrix given by S µ 1 µ 2 (θ) = S(θ) δµ 1 ,µ 2 .

Form factor expansions of two-point functions
In this paper we study the correlators T (r)T (0) and, especially, : T φ : (r) :T φ : (0) and φ(r)φ(0) , as well as the associated entanglement entropy obtained via (14). It is well known that two-point functions of local operators in QFT can be expressed as infinite sums involving matrix elements of these operators. The matrix elements of relevance, also known as form factors, are defined as for a local field O. Here |0 represents the vacuum state and |θ 1 , . . . , θ k in µ 1 ,...,µ k are the physical "in" asymptotic states of massive QFT. They carry indices µ i , which are quantum numbers characterizing the various particle species, and depend on the real parameters θ i , which are called rapidities. The energy and momentum of a particle of mass m i are expressed in terms of its rapidity θ i as m i cosh θ i and m i sinh θ i , respectively. In terms of form factors, two-point correlation functions (in unitary models) may be expanded as As mentioned, the Lee-Yang model is non-unitary. As was noted in [65,66], non-unitarity affects the form factor expansion. A consequence of this is that many fields appear to be non-Hermitian under the Hilbert structure of asymptotic states. In the Lee-Yang model, an exact calculation of the form factors of the field φ shows that 0|φ(0)|θ 1 , . . . , θ k in * = in θ 1 , . . . , θ k |φ(0)|0 , where the right-hand side can be obtained by crossing symmetry. However, it turns out that the relation is surprisingly simple: As a consequence, the form factor expansion of the two-point function of the field φ, normalized by the square of the VEV, is a modification of (21) where sign factors (−1) k are included for the terms involving the k-particle form factors. This gives rise, in the single-copy Lee-Yang model, to: A natural way to understand this modification is through a discussion of the bound-state singularity occurring in the form factors. The additional (−1) k guarantees that the boundstate residue of the analytic continuation of the k-particle integrand, which, like that of the scattering matrix, has the wrong sign, is related to the k − 1-particle integrand in a way that would guarantee locality properties. As we will see below, form factors of twist fields : T φ :, :T φ :, T andT are subject to similar bound-state residue equations as those of φ. Hence, this interpretation suggests that a similar modification of (21) occurs for the form factor expansion of : T φ : (r) :T φ : (0) and of T (r)T (0) . That is, in the n-copy model, for both O = T ,Õ =T and O =: T φ :,Õ =:T φ :. Here we have used the fact that by symmetry under inversion of copies, T = T and :T φ : = : T φ : , and the form factor expansion includes sums over the copy numbers µ j . That this is the correct expansion follows from an equation similar to (22) for twist fields, see subsection 4.4 Finally, since the field φ is no longer Hermitian, its VEV is no longer expected to be real. As was shown by Zamolodchikov, φ is in fact purely imaginary -this can be explained by the fact that it occurs in the formal massive Lee-Yang action (written as a perturbation of the CFT action) with a purely imaginary coupling constant. A similar phenomenon makes the VEVs T and : T φ : not necessarily real. We will determine their phases (up to multiples of π) by evaluating analytically the normalization of their leading short-distance power-law, and by observing numerically that the right-hand side of (24) is positive for all mr.

Short-distance behaviour from form factors
Form factor expansions (21) and (23) are naturally large-distance expansions, in that they converge very rapidly for large values of rm. However, in many cases we want to explore small values of rm. In such cases two-point functions generally develop power-law behaviours in rm and it is very difficult to extract the precise power from a form factor expansion such as those above.
It was realized a long time ago [65] (see also [67] for a nice derivation and application to various models and [68] for a generalization to boundary theories) that if one is interested in the short-distance behaviour of correlators then an expansion of the logarithm of two-point function is more appropriate: The functions H O|µ 1 ,...,µ k k (θ 1 , · · · , θ n ) must of course be chosen so that the expansion (21) is recovered when exponentiating (25). This condition automatically implies for example that In general the H k functions can be interpreted as the "connected parts" of the F k functions (they are "cumulants" with respect to the rapidities). These are such that, if the clustering decomposition holds for the F k 's at large rapidities for all k, that is ∀k, ∈ N, then the H k 's vanish at large rapidities for all k. Thanks to this vanishing, for mr 1 we now expect each summand in the sum over k in the expression above to be dominated by a leading term proportional to log(mr). The constant coefficient of this term, summed over all particle contributions, will then give the power which governs the short-distance behaviour of the two point function. Let us call this power −4x O . Then, carrying out one integral in (25) and expanding the result for small mr we find [65,67,68] x This was used for the field φ in [66] and shown to agree well with conformal field theory results.
In addition, the proportionality constant (17) of the power law behaviour at short distances, can also be extracted from the form factor expansion. It was shown in [67] that by considering the leading correction to the log(mr) term in (25) one may also find a form factor expansion for the constant K O which is given by: 4 Twist field form factors

Form factor equations and minimal form factors
The form factor equations for Z n twist fields were derived in [44]. Details of the solutions procedure for two-particle form factors appeared there, and higher particle form factors of various models were computed in [45] and [69]. Interestingly, a very similar set of form factor equations had been derived much earlier [70] in a rather different context (e.g. the study of the response of an integrable QFT to a variation of the Unruh temperature). The details of the computation for the Lee-Yang model are very similar to those described in these works, with the only difference that the presence of the bound state pole in the S-matrix imposes further conditions on the form factors. In particular, bound state poles are present in addition to kinematic poles. The form factor equations only encode locality properties of fields, hence they are unchanged for form factors of any field in a the same Z n twist sector. In order to distinguish for form factors of T and : T φ :, we will impose additional conditions, and verify the correctness of the solutions by numerical comparisons with CFT predictions.
In what follows we will consider k-particle form factors F O|µ 1 ···µ k (θ 1 , θ 2 , . . . , θ k ) for a generic twist field O. We will later identify this field with T or : T φ : depending on various properties of the form factor solutions we obtain.
The two-particle form factor must satisfy (in the two-particle case we use the single argument and the kinematic residue equations Res Here and below we useμ = µ − 1 mod n. Higher particle versions of these equations read Res For this model there is the added difficulty of having to solve also the bound state residue equation associated to the scattering process a + a → a where a is the Lee-Yang particle on copy a. This takes the form where the so-called three-point coupling is fixed by and by choosing the negative imaginary direction: Γ = −i2 1/2 3 1/4 . For n = 1 this equation fixes the one particle form factor (which for spinless fields must be rapidity independent) through the equation These equations imply that the two-particle form factor solution given in [44] must be generalized to include the bound state pole. As discussed in [66] this may be done by defining a minimal form factor. A minimal form factor is a solution of (33) which has no poles in the (extended) physical sheet θ ∈ [0, 2πn) except possibly for bound state poles, and which tends to unity as |θ| → ∞. It turns out that this particular Riemann-Hilbert problem has a unique solution, and this solution possesses bound state poles with nonzero residues: where a(θ, n) encodes the bound state pole and f (θ) is given by the integral representation The latter function admits also a representation as an infinite product of gamma functions which was already given in [44] for the sinh-Gordon model (it suffices to take B = 2/3 and to invert the formula).
The expression (43) may be obtained as a solution to (33) using a similar integral representation of the two-particle scattering amplitude. In the absence of bound state poles, the resulting f (θ, n) would directly be the minimal two-particle form factor. In the present case, however, the function tends to 1 as |θ| → ∞ but has a simple pole at θ = 0. The factor a(θ, n) is the unique one that shifts this pole towards the position of the allowed bound-state singularity in the physical sheet, without affecting the large-|θ| behaviour.
Using the integral or Gamma-function representation, it may be shown that In order to compute higher particle form factors the following more general identities are important

Twist field one-and two-particle form factors
One expects that primary twist fields, with direct geometric meaning, would occur as solutions to (33) and to conditions of bound state and kinematical singularities with the additional requirement of convergence as |θ| → ∞. This additional requirement is expected to implement, in a path-integral picture, the least singular asymptotic condition possible at small distances near the position of the field. With these conditions, the most general form the two-particle form factor can take is where the first term is of the form required to solve the kinematic residue equation (and of the same form as for other theories previously studied [44]) and the second term is what is commonly termed a "kernel" solution of the kinematic residue equation (that is a solution without kinematic poles).
In general κ is an arbitrary constant, but it may be fixed by imposing the cluster decomposition property, namely where we have used the fact that lim θ→∞ F 11 min (θ) = 1. Then, the one-particle form factor on copy a may be fixed by combining this with equation (40), which translates into the following quadratic equation for F This leads to two possible solutions: where we have used the identity (44). The presence of two solutions immediately suggests the existence of two different leastsingular twist fields, by contrast to other models studied in the past. It is natural to conjecture that these are T and : T φ :, and given this, it is a simple matter to identify their respective form factor solutions. Indeed, the former specializes to the identity at n = 1, and the latter, to φ. We note that the solution with the negative sign specializes to 0 at n = 1, and that that with the positive sign specializes to the one-particle form factor of the field φ found in [66] (note that the constant v(0) in [66] is , 1) with our present notation, and that the coupling h was computed in [71]). These properties suggest the identifications The numerical results of sections 4, 5 and 6 provide further support for these identifications.

Higher particle form factors
Let us now consider only form factors of the form F , that is form factors involving only one particle type. This is sufficient as form factors involving other particles may be obtained from these by using the twist field form factor equations [44].
The higher particle form factors may be obtained by making the ansatz where x i = e θ i /n and α = e iπ/n . The functions Q k (x 1 , ..., x k ) are symmetric in all variables and have no poles on the physical sheet. This ansatz, as usual in the context of the computation of form factors of local fields (see e.g. [66,72]), expresses the form factors in such a way as to explicitly separate the part containing the poles from the part which has no singularities. In addition, the explicit presence of the minimal form factor and the symmetry in the variables x i automatically gives form factors which exhibit the correct monodromy properties in the rapidities. In the context of twist fields, this ansatz was used for the first time in [69].

Kinematic and bound state residue equations
Using (46), the kinematic residue equation with the ansatz (53) can be rewritten as (k ≥ 0): where P k is the polynomial where β = e − 2πi 3n and Denoting σ (k) i the i-th elementary symmetric polynomial on k variables x 1 , . . . , x k , which can be defined by means of the generating function, we can rewrite P k (x 0 , x 1 , . . . , x k ) as In the following we will omit the upper index (k) when there is no confusion possible. Besides (54), another equation that arises from the ansatz (53) is that using the bound state residue equation (38). The simplest case of this equation was given in (40) and this allowed us, in combination with the clustering property, to fix the one-particle form factor (49). For higher particles, using (45) we find (k ≥ 1) with From the ansatz (53) it follows that Q 1 = F O|1 1 .

Three-particle form factors
First let us analyze the two-particle case. We have by definition that F O 0 = Q 0 = O , and comparing to (47), we obtain the polynomial It is a simple matter to verify that this is indeed in agreement with the kinematic residue equation (54); given Q 0 this is the most general solution to (54) (k = 0), as was shown in [69] (in particular, the second term vanishes at x 1 = αx 2 ). Further, replacing (F O|1 1 ) 2 / O by the linear combination of the zero-and one-particle form factors, Q 0 and Q 1 (x 1 ), occurring via the quadratic equation (49), one can check that (62) is in agreement with (59). In fact, given arbitrary Q 0 and Q 1 (x 1 ), the resulting expression is the unique solution to (54) (k = 0) and (59) (k = 1).
As was shown above, the additional condition of clustering imposes the one-particle form factor to take only two possible values (proportional to the vacuum expectation value), according to (52). For n = 1 the solution (62) is either zero (if we take the first solution in (52)) or it reduces to Zamolodchikov's two particle solution for the Lee-Yang field [66] (if we take instead the second solution in (52)). This is in accordance with identifying the two-particle form factors with those of T and : T φ :, respectively.
Interestingly, it turns out that the above structure subsists to higher particles: given Q 2 (x 1 , x 2 ) and Q 1 (x 1 ), there is a unique solution to the kinematic and bound state residue equations (54) (k = 1) and (59) (k = 2) for the polynomial Q 3 (x 1 , x 2 , x 3 ). The solution has the following structure: where the parameters A i are complicated functions of n but rapidity-independent. The detailed computation of Q 3 (x 1 , x 2 , x 3 ) and the values of A i are reported in appendix C, and note in particular that the polynomials σ 6 1 and σ 4 1 σ 2 have vanishing coefficients. Again it is interesting to consider the limit n → 1 of the functions A i above. Using the two solutions (52), we now note that all constants vanish, A i = 0, when we consider that corresponding to the operator T (where F T |1 1 = 0 for n = 1), thus the three particle form factor also vanishes. On the other hand, if we consider the other solution in (52), which at n = 1 should correspond to the field φ, we find and a simple computation shows that our three-particle form factor, as expected, reduces to Zamolodchikov's solution [66]. It is tempting to use this benchmark (agreement with Zamolodchikov's solutions) to try and find the general solution for higher particle numbers. However, as the three-particle case shows, the reduction to n = 1 occurs thanks to great simplifications. At this stage, it is unfortunately not obvious at all how high-particle solutions may be constructed other than by brute force computation. The main reason for this is the presence of two (rather than one) kinematic pole in the form factor ansatz (53). This leads to polynomials Q O k (x 1 , . . . , x k ) of much higher degree than is the case in the standard form factor program.
Despite the complexity of the expression (162), there are certain simplifications that can be used to rewrite the three-particle form factor in a form which is more suitable for numerical computations. It turns out that where f 3 (θ 1 , θ 2 , θ 3 ) is the function that is obtained from Q 3 (θ 1 , θ 2 , θ 3 ) in (162) by setting all terms proportional to O −1 to zero. In other words, when divided by i<j (x i − αx j )(x j − αx i ), all those terms simplify giving just the second summand in the formula above. This summand represents a kernel solution to the form factor equations, in the sense already described in subsection 3.2. Finally, note that where we have used the property which can easily be derived from (49), (56) and (61). In other words, the three-particle solution automatically satisfies the clustering property. This is an extremely nontrivial check of the validity of the three-particle solution. This situation is in contrast to that of the sinh-Gordon model [69], where at each particle number, the clustering property has to be imposed in order to uniquely fix the solution. It also follows from the result above and the cluster property of the two-particle form factor that Properties (66) and (68) are very important as they insure the convergence of the integrals (25) for k = 3.

Form factors of the fieldsT and :T φ :
In the previous subsections we have concentrated our analysis on computing the form factors of the fields T and : T φ :. However, the correlators we are interested in also involve the fieldsT and :T φ : thus their form factors are also required. In fact the form factors of all these fields are not independent from each other. We may think of T and : T φ : and ofT and :T φ : as twist fields associated to the two opposite cyclic permutation symmetries i → i + 1 and i + 1 → i (i = 1, . . . , n, n + 1 ≡ 1). From the additional symmetry under the inversion of copy numbers it follows that F and similarly for : T φ : and :T φ :. At the same time, as already explained in subsection 3.2, from the non-unitarity of the theory we would expect that (note that T = T ). These equations both define the form factors ofT and impose the condition expressed by the last equality above on the form factors of T . We have verified that this is satisfied for all our solutions, and that similar equations hold for : T φ :. These equations are the counter-part of (22) for twist fields, and show that the form factor expansion (24) is correct. Finally, another important relation which we have used in subsequent computations is the following identity which allows us to express any form factor in terms of form factors involving only the particle living in copy 1. The same equation holds for the field : T φ :.

Identification of twist field operators: numerical results
In previous sections we have provided compelling evidence for the identification of the two families of form factor solutions that we have obtained with the twist fields T and : T φ :. This evidence is based on the (highly non-trivial) fact that the one-particle and higher form factors of the field we identified as T vanish at n = 1 whereas those of : T φ : reduce to the form factors of φ obtained in [66]. Further evidence may be gathered by, for example, examining the short distance behaviour of the correlators T (r)T (0) and : T φ : (r) :T φ : (0) . We must therefore first understand what the expected behaviour of such correlators should be for the theory at hand.
Let us first consider the conformal field theory. In CFT such correlators are expected to converge at small distances as and Indeed note that the powers above are positive for the Lee-Yang model as both c and ∆ are negative (see section 3). This is of course a consequence of non-unitarity.
In the massive theory, however, we expect that the leading short distance behaviours of these correlators should be described by a different power law: and The reason for this is entirely analogous to the observation made in [66] regarding the correlator φ(r)φ(0) . It was found that for short distance in the massive theory the leading behaviour of this correlator was r −2∆ rather than the conformal behaviour r −4∆ . Zamolodchikov argued that this was due to the fact that the leading behaviour of the conformal OPE comes from the field φ rather than the identity. In the massive theory the expectation value φ = 0 and therefore the contribution to the OPE from the field φ itself becomes the dominating term in the short distance expansion of the two-point function.
Similarly, it is possible to argue that the leading contribution to the OPEs of T andT and of : T φ : and :T φ : corresponds to the field φ 1 φ 2 . . . φ n where φ i represents the field φ on copy i. This field has dimension n∆, and it is the field of smallest (most negative) conformal dimension that can be constructed in the n-copy Lee-Yang model. Since its expectation value is nonzero, it thus gives the leading contribution at short distances. Massive OPEs of twist fields will be discussed in more detail in section 6.
Thus, by employing a form factor expansion we may check whether the expected behaviours are indeed recovered from our form factor solutions. We will include up to three particle form factors as done in [66]. We have performed a numerical evaluation of the formula (30) including up to three particle form factors for the twist fields T and : T φ :. We confirm with good accuracy that the twist fields exhibit the behaviours (74) and (75) Table 1: Study of the two-point function T (r)T (0) of the n-copy Lee-Yang theory at short distances. Near the critical point we expect this correlator to exhibit a power-law behaviour of the form r −4x T where x T = ∆ T − n∆ 2 = − n 12 + 11 60n . This value should be best reproduced in the massive theory the more form factor contributions are added. The data above show that this expectation is indeed met by considering up to three-particle form factors.  Table 2: Study of the two-point function : T φ : (r) :T φ : (0) of the n-copy Lee-Yang theory at short distances. Near the critical point we expect this correlator to exhibit a power-law behaviour of the form r −4x :T φ: where x :T φ: = ∆ :T φ: − n∆ 2 = − n 12 − 1 60n . The data above show good agreement with CFT by considering up to three-particle form factors.

Comparison with perturbed conformal field theory results
A further consistency check of our form factor solutions may be carried out by comparing a form factor expansion of the correlators T (x 1 )T (x 2 ) and : T φ : (x 1 ) :T φ : (x 2 ) to its counterpart in perturbed conformal field theory.

Conformal perturbation theory and twist fields structure constants
We may regard the action of the integrable quantum field theory as a perturbation of Lee-Yang CFT action by a term proportional to a coupling constant λ and the CFT field φ(x,x) of conformal dimension ∆, and compute correlators by performing perturbation theory about the conformal critical point on the coupling λ [73]. As is well-known, in the massive theory this coupling is related to the mass scale m as λ ∝ m 2−2∆ (with a known proportionality factor [71]). The massive correlators can then be obtained by using OPEs where conformal structure constants are modified to structure functions of mr that can be evaluated perturbatively in λ (they are convergent series in integer powers of λ), and where vacuum expectation values, which are non-perturbative in for n ≤ 11. The squares, circles and triangles, represent the up to one-, two-and three-particle form factor contributions. The black solid line represents the exact values at criticality. All curves clearly show strong linearity in n which is consistent with the CFT behaviour, where the coefficient of n (e.g. slope of the curves) approaches the CFT value as more form factor contributions are added. The agreement with CFT gets worse as n increases. This is also to be expected as the larger n is, the larger the contribution of higher particle form factors becomes (all form factor contributions are in fact proportional to n).
λ, are nonzero. The same type of comparison between a form factor and a perturbed CFT computation was carried out in [66] for the two-point function of the field φ in Lee-Yang. Let us now consider the OPEs of T withT and of : T φ : with :T φ :. They involve only fields in the non-twisted sector (we mean by this all fields constructed by considering n non-interacting copies of the fields of the original theory) and by construction they must be invariant under cyclic permutation of the copies. Let us consider the following primary, cyclically invariant, homogeneous fields, composed of multilinears in the fields φ i on the various copies: we label them by sets {k 1 , . . . , k J } of J different integers in [1, n] for J = 1, 2, . . . , n (we may take k 1 < · · · < k J ), and take them to be Φ k 1 ,...,k J := φ k 1 · · · φ k J + cyclic permutations S k 1 ,...,k J .
The symmetry factor S k 1 ,...,k J is equal to the order of the subgroup of the cyclic replica permutations which preserve the sequence k 1 , k 2 , . . . , k J of replica indices. That is, Φ k 1 ,...,k J is the sum, over all elements σ ∈ Z n in the cyclic replica permutation group Z n , of σ(φ k 1 · · · φ k J ), divided by the order of the stabilizer, in Z n , of φ k 1 · · · φ k J . This definition guarantees that in Φ k 1 ,...,k J , every independent multilinear term, including the initial term φ k 1 · · · φ k J itself, appears with coefficient 1. The number of independent multilinears in Φ k 1 ,...,k J is n/S k 1 ,...,k J . The symmetry factors for low values of J can be written explicitly: S 1,k = 2 (n even, k = n/2 + 1) 1 (otherwise) The fields Φ k 1 ,...,k J have conformal dimensions J∆. In order to have a basis of primary, cyclically invariant homogeneous fields of dimension J∆, we need to further restrict the indices k 1 , . . . , k J . We may certainly fix k 1 = 1, and further restrictions hold due to the residual equivalence relation generated by {1, . . . , k J } ∼ {1, n+2−k J , n+1+k 2 −k J , . . . , n+1+k J−1 −k J }. More generally, in the set of all replica-index sets {k 1 , . . . , k J }, there is a foliation by Z n orbits, and a basis of fields Φ k 1 ,...,k J can be taken as fields parametrised by single representatives of each Z n orbit.
The OPEs in the massive theory can be regarded as "deformations" of the conformal OPEs such that the structure constants are replaced by functions of mr. Denoting by O andÕ any given pair of conjugate (i.e. whose twist actions cancel out) twist fields, it takes the form where r := |x 1 − x 2 |, m is the physical mass of the Lee-Yang model. The functions admit an expansion in integer powers of the coupling λ, hence in powers of (mr) 2−2∆ , and the constantsC Φ k 1 ,...,kp OÕ are the structure constants of the CFT. In our analysis we will in fact only consider the leading term (the CFT contribution) to these structure functions, that is, we will only carry out zeroth order perturbation theory whereby the mass dependence is introduced through the non-vanishing expectation values of OPE fields. The analysis is still non-trivial because of the presence of nonzero expectation values. Note that with the definition (11) and the standard definition of T andT we have the conformal normalizatioñ The conformal OPEs (and structure constants) of the branch point twist field T have been studied in several places in the literature. The most general study can be found in Appendix A of [74] where general formulae for the structure constants associated to the OPE of T withT in general (unitary) CFT are given. Structure constants have also played an important role within the study of the entanglement of disconnected regions [75,76]. More recently the structure constants of other types of twist fields which arise naturally within the study of the negativity have been studied in [77,78]. However we do not know of any studies of the OPE and structure constants of composite fields such as : T φ :. Here we provide explicit step-by-step computations of the conformal structure constantsC Φ 1 TT ,C Φ 1 :T φ::T φ: ,C (see appendix D for details) which are proportional to one-, two-, three-and four-point functions of the field φ (other structure constants would involve higher-point functions, which are harder to access). Our computations focus on the Lee-Yang model but could be easily generalized to other minimal models (for the field T the ingredients needed for such generalization are already provided in [74]). Other structure constants and massive corrections thereof will involve higher point functions. The results are: where κ is a model-dependent function which characterizes the four-point function of fields φ Other structure constants may be computed in terms of higher-point functions so that in general we expectC But for the last line in (82) and for (85), all formulae above are particular cases of those given in [74].
The difficulty of calculating such terms is then reduced to the difficulty of obtaining higherpoint functions in CFT. Such higher-point functions will also be required in order to obtain most massive corrections to the CFT structure constants, a problem which we will not be addressing in this work.

The case n = 2
As explained earlier, obtaining the CFT structure constants becomes a difficult problem for the field : T φ : as soon as we consider OPE terms involving products of more than two fields and for the field T when we consider products involving more than four fields. For this reason, the case n = 2 is particularly interesting as in this case the leading contribution to the OPE is given by the bilinear fields Φ 1,2 = 2φ 1 φ 2 defined earlier. The leading expansions in the massive theory are Here the numerical coefficients arise from the total numbers of independent multilinears in Φ 1 and in Φ 1,2 in the case with n = 2, which are n/S 1 = 2 and n/S 1,2 = 1 respectively. All subleading terms correspond to Virasoro descendants and massive corrections to the structure constants, hence are suppressed by positive powers. The structure constants arẽ In the Lee-Yang model all the constants (89) can be computed. The CFT structure constant C φ φφ can be found for instance in [66], The four point function of the Lee-Yang model has been studied in [79,80,22]. Following [22] we can write the four point function as in (83) with where A simple calculation then gives which gives an approximation of the two-point function at zeroth order in perturbed CFT. In order to compare this with the form factor expansion, we need to fix the vacuum expectation values of T and : T φ :. Recall that we used the CFT normalization to set the coefficients of r 11 10 and r 3 2 , respectively equal 1. This in principle uniquely fixes the expectation value. Although the resulting expectation value is not known explicitly, we may use the form factor expansion (32) to estimate it. From CFT the leading behaviours should be : We observe from Figure 3 that the truncated form factor expansions of the two-point functions potentially change sign (pass by the value 0) only at short distances, at positions that become smaller as more particles are added. Since they approach the CFT form at short distances, this implies that the full two-point function never becomes zero. Since the ratios T (r)T (0) tend to unity at large distances, they are then positive for all values of r. Hence, we find T 2 < 0 and : T φ : 2 > 0.
In fact, from (96) and (97), we have that The constants K T and K :T φ: as expressed in (32) are necessarily positive, and the fact that twopoint functions never become zero is related to the convergence of the series (32). A numerical evaluation of (32) including up to three-particle form factors yields

The cases n = 3 and n = 4
For the fields T andT it is also possible to compute the two-point function in the zeroth order approximation for n = 3, 4. It is given by for n = 3 (where we have used the numerical coefficients n/S 1,2 = 3 and n/S 1,2,3 = 1), with whereas for n = 4 we have where we can again easily compute This gives T (r)T (0) = r +· · · , (108) for n = 4.
Like for the n = 2 case, it is possible to compare these results to a form factor expansion once the expectations values of T and : T φ : have been obtained by using (32). In the three-particle approximation we find K T = 2.02966 and K :T φ: = 2.60713 for n = 3, K T = 2.89127 and K :T φ: = 3.48758 for n = 4, giving and T 2 = −41.5266 for n = 4. (112)

Summary and discussion
In this section we have studied the two-point functions : T φ : (r) :T φ : (0) and T (r)T (0) by using two well-known approaches: a form factor expansion (up to 3 particles) and perturbed CFT (at zeroth order) calculation. Examining Figures 3 and 4 we can say that agreement between both approaches is good in terms of the range of values that the correlators take but not particularly good if we compare the slope and precise values the functions take at particular points. This level of agreement (and disagreement) is not entirely surprising given the expected range of validity of each approach: the form factors approach is eminently a large mr expansion and although considering contributions up to three particles should provide a relatively good description for small values of mr we do not expect it to be very precise for very short distances. Conformal perturbation theory works best near criticality, that is for very small values of mr, exactly where form factors should be less accurate. Besides, we have carried out perturbed CFT at zeroth order so the expectation is that this should really only be accurate for very small values of mr. Finally, a numerical comparison between CFT and form factors is only possible if the form factor normalization constant (that is the vacuum expectation value of the field) is known. In our case we can only access these expectation values approximately through yet again a form factor expansion. This introduces a further error (the vacuum expectation values obtained this way are smaller in absolute value than their exact values) which results in an overall shift of the form factor points.
Overall the results we obtain are not dissimilar to Zamolodchikov's results [66] for the twopoint function φ(r)φ(0) in Lee-Yang. Agreement with CFT was slightly better in [66] as first order corrections in perturbed CFT were also included and the exact value of φ 2 was known from an independent thermodynamic Bethe ansatz computation.
Despite the many limitations described above, it is still the case that agreement between form factor numerics and zeroth order perturbed CFT is better for some particular correlators than for others. We do not have a good physical explanation as to why this should be the case but it appears to depend on the particular functional form of the perturbed CFT curve obtained for each case, that is the relative weight of the various contributing terms and the region of values of r where the term with the lowest power of r is leading. In the previous section we established that although the form factor expansion is only rapidly convergent for mr 1, it does still provide a good estimate of the short distance behaviour of correlators. An alternative way of testing this results is by performing a computation of the Rényi entropy as defined in (14). This involves also the computation of φ(r)φ(0) , which was first obtained in [66] (and which can be obtained from : T φ : (r) :T φ : (0) by setting n = 1). Figure 5 shows the results of such a computation for n = 2, 4, 6 and 8.

Bi-partite entanglement entropy of large subsystems
In this section we will used the form factors previously obtained to study the bi-partite entanglement entropy of the Lee-Yang model paying special attention to the region mr >> 1.

Saturation
Using (14), the entanglement entropy of non-unitary theories is where we have used the short-hand notation Observe that lim n→1 A(r, n) = B(r) and that lim mr→∞ A(r, n) = lim mr→∞ B(r) = 1. 1 + 1 n log(mε), and evaluated in logarithmic scale using form factors. The form factor contributions up to one-, two-and three-particles are considered both for the correlators : T φ : (r) :T φ : (0) and φ(r)φ(0) . The solid line represents the CFT prediction c eff (n+1) 6n log(mr) (note that for mr < 1, this is negative). All graphs show a clear logarithmic divergence at mr = 0 (as expected). Additional form factor contributions (also as expected) improve agreement with CFT.
The expression above can be written as The constant is a convenient short-distance cutoff related to ε by c eff 3 log(m ) = c eff 3 log(mε) + lim n→1 d dn Z n − (C φ φφ ) n /C Φ 1,...,n :T φ::T φ: . The dimensionless, universal saturation constant is ..,n :T φ::T φ: where K O was defined in (17) (equivalently (31)). The meaning of U is clear from a comparison of the small-and large-distance behaviour of the entanglement entropy: These easily generalize to the Rényi entanglement entropy at arbitrary n, giving the saturation behaviour (15) with the universal saturation constant U n expressed in (16), and in particular U = U 1 .

Leading order correction to saturation
It is easy to see that where A := dA/dn. We now need to compute the objects above, all of which are given in terms of various two-point functions and their limits at n = 1. As we have seen in the introduction, both the two-point functions A(r, n), B(r) and the logarithm log B(r) admit expressions in terms of form factors. Here we only want to investigate the first and second order corrections to saturation of the entanglement so we will consider only up to the two-particle contribution to the form factor expansion. By doing so we have that A(r, n) = 1 + A 1 (r, n) + A 2 (r, n) + · · · (120) where is the one-particle contribution, and is the two-particle contribution. Similarly, where the ratio F φ 1 φ was given in (51), and The form factor F φ 2 (θ) was given by Zamolodchikov in [66] and can be written as Thus, we find that where the first two terms will give the next-to-leading order contribution to the entanglement entropy (i.e. the leading correction to its saturation value) and the remaining terms give the next-to-next-to leading order contibution. We will now analyse this expression in more detail.
In appendix E we show that A 1 (r, 1) is given by We also need A 2 (r, 1) which is given by This simple result was established in [44] for all integrable quantum field theories and even beyond integrability [46]. Then, the expression above simplifies to lim n→1 d dn Thus, the von Neumann entropy of the Lee-Yang model takes the form where U is the model-dependent constant (117) and and In contrast to results found for unitary theories [44,46], the results above suggest that the leading and next-to-leading order correction to saturation of the entropy of large blocks are strongly model-dependent. In particular, the leading correction is proportional to the constant a which clearly depends on specific features of the model under consideration (that is, the one-particle form factor). This term is directly related to the one-particle form factor and in particular to its value and the value of its derivative at n = 1. The fact that both these quantities are non-zero for : T φ : is special for this field -they would have been zero if we had used Tand we are tempted to conclude that this phenomenon is related to the non-unitary nature of the model. It would be interesting to test these results numerically for example by studying the spin chain model considered in [19,20,29,43].

Conclusions and outlook
In this paper we have provided an in-depth study of the two-point functions of twist fields in the massive Lee-Yang model and their application to the computation of the bi-partite entanglement entropy. The main tools used for our study are branch-point twist fields and the relationship between their correlation functions in replica theories and the bi-partite entanglement. For massive unitary theories this connection was established and explored in [44], and the present work addresses the problem for a massive non-unitary model for the first time. Representing the EE using correlation functions of twist fields indeed provides the only known method so far for performing computations of the bi-partite entanglement in massive QFT models. The non-unitarity of the theory has important consequences for the computation of entanglement and several stark differences are found with respect to the unitary case.
The twist field T and its conjugateT considered in [44] are not the right operators to consider in non-unitary models when performing entropy computations. Instead, as first proposed in [43] for CFT, one must consider the two-point function of suitably normalized composite fields : T φ : and :T φ : introduced in [51], defined as leading contributors to the OPE T (x)φ(0) andT (x)φ(0) where φ is the lowest-dimension primary field of the model.
In the present work we find the exact form factors of T and : T φ : up to three particles. Remarkably we find that the form factor equations together with the requirement of clustering are sufficient to entirely fix all form factors and to provide in a natural way two families of solutions corresponding to the two twist fields T and : T φ :. We also give numerical evidence that the resulting correlation functions agree, at short distances, with CFT results. This is done by numerically evaluating truncated form factor expansions, or the logarithm thereof, of correlation functions at short distances, and comparing with a zeroth order perturbed CFT computation of the twist field two-point function, in the spirit of Zamolodchikov's work [66]. The CFT computation also provides some of the first general results regarding OPEs of the composite twist fields : T φ : and :T φ :. Finally these results are used to compute the Rényi and von Neumann entropy, in particular in the limit of large blocks, which is well described by the form factor expansion. For large blocks we find that the corrections to saturation are strongly model-dependent for non-unitary theories. This is in contrast to the very universal form of the leading correction found for unitary models [44,46], only depending on the particle spectrum.
It would be interesting to compare the present results about entanglement entropy with a numerical evaluation in the quantum Ising model with imaginary transverse magnetic field spin first considered by von Gehlen [19,20], whose near-critical universal region is expected to be described by the Lee-Yang QFT [21,22]. This would provide the first strong evidence beyond criticality supporting the conjecture that the composite fields : T φ : and :T φ : are the correct fields for representing branch points, or conical singularities, in non-unitary models.
The ideas in the present work and in [43], in particular the form (14) of the twisted replica partition function in non-unitary QFT, lead us to speculate a relation between correlators of composite fields in non-unitary QFT and correlators of physical fields in its "unitary counterpart". Following ideas from the field of PT-symmetric quantum mechanics [23,24], given that the non-unitary theory (i.e. with a non-Hermitian hamiltonian) considered has a real energy spectrum, we may infer that there must be a similarity transformation which maps the Hamiltonian and correlators of the Lee-Yang model to the Hamiltonian and correlators of some unknown unitary theory. It is tempting to propose that the operation of taking composites with the lowest-dimension field φ implements, up to normalization, such a similarity transformation.
It would be very interesting to test this and related ideas further. Several future directions of research follow naturally from this work: a more detailed study of the OPEs of twist fields for arbitrary n is desirable, not only to better understand the properties of replica CFTs but also as building blocks for perturbed CFT computations for larger values of n. A systematic understanding of higher particle form factor solutions is still missing as is the study of the entanglement entropy of excited states and multipartite regions in massive (unitary or not) QFT.
A Definition of the field : T φ : The power a is fixed by requiring that the limit exist, and the normalization A is determined by requiring conformal normalization of the resulting field. These, as well as structure constants studied in Appendix D, may be evaluated by using standard methods of CFT [72]. Correlation functions with twist field insertions at y 1 and y 2 in the n-copy model are interpreted as correlation functions on a n-sheeted Riemann surface M n,y 1 ,y 2 with branch points in place of the twist fields, and conformal uniformization to the sphere is used. For the uniformization step, one makes use of the conformal map g : M n,y 1 ,y 2 →Ĉ \ {0, ∞}, g(z) = z − y 1 z − y 2 with ∂g ∂z := ∂g = 1 n In order to compute A and a, we compute the following ratio of correlation functions: The new ratio of correlators involved is interpreted as a correlator of φ j 1 (x 1 )φ j 2 (x 2 ) on the Riemann surface M n,y 1 ,y 2 , and can be computed by using the conformal map above to relate them to correlators in the complex plane. Thus Hence we must set We find the known result [51] that the dimension of : T φ : is and we have the correct normalization B Large n expansion of the one-particle form factor contribution The one particle contribution to the powers x T and x :T φ: as defined in section 5 listed in the tables 1 and 2 may be computed exactly as is simply given by the function (see eq. (30)) for O = T or O =: T φ :. It is easy to show that this function admits a large n expansion in powers of 1/n starting with a term linear in n. Recall the expressions (52). Combining those with (44) we can rewrite the expectation values as 2π sin π 6n sin π 2n f (iπ, n) For large n we find that sin π 3n 2π sin π 6n sin π 2n cos π 3n + 2 sin 2 π 6n 2 = 2n and sin π 3n 2π sin π 6n sin π 2n cos π 3n − 2 sin 2 π 6n 2 = 2n We now study the expansion of f (iπ, n) −1 for n large. It is possible to show that f (iπ, n) −1 = ∞ k=0 f k (n) with f k (n) given by the following expression Γ kn + 1 2 Γ kn + 7 6 Γ kn + 4 3 Γ (k + 1)n − 1 2 Γ (k + 1)n + 1 6 Γ (k + 1)n + 1 It is easy to see that the leading contribution for n-large comes from the n-independent part of the k = 0 term in the product. We have that and so on. This gives where the coefficient of 1/n 2 is not exact but has been obtained by considering contributions up to f 5 (the next term would have given a correction of 2 × 10 −3 to the coefficient). The expansion above indeed shows that the structure of the one-particle form factor contribution closely matches what is expected from CFT since Indeed, comparing coefficients we find that the one-particle contribution provides 56% of the coefficient of n for both T and : T φ :, 65% of the 1/n coefficient for the field : T φ :, and 49% of the same coefficient for the field T . In short, the one-particle form factor provides a very substantial contribution to the two-point function of twist fields, both for short and long distances.

D Conformal structure constants of twist fields
In this appendix we present detailed computations of the conformal structure constantsC O TT andC O :T φ::T φ: for different choices of the local field O. These structure constants are used in section 6 in zeroth order perturbed CFT computations. The general strategy relies upon the fact that correlation functions of twist fields in CFT may be computed in two different ways: on the one hand we may treat the twist fields as standard local fields in the n-copy model on the manifold M n,x 1 ,x 2 (as defined in 138); on the other hand we may conformally map the correlation function to the complex plane by using the map (138)

TT
The CFT structure constantC Φ 1 TT may be computed as follows. We may select out the term proportional to Φ 1 = n j=1 φ j in the OPE of T (x 1 )T (x 2 ) by evaluating the three point function below, where a single field φ 1 is inserted: On the other hand, we identify the ratio of correlators on the left-hand side as a correlator of φ 1 (x 3 ) on the manifold M n,x 1 ,x 2 , and use the conformal map g to relate this to the one-point function φ(g(x 3 ) on R 2 . Since this one-point function is zero in the complex plane, we have that, in generalC The third term in the OPE T (x 1 )T (x 2 ) contains bilinears of the fields φ j , with the constraint that they be cyclically symmetric. The only possibility are the fields Φ 1,k defined earlier, with k = 2, . . . , [n/2] + 1; the restriction on k is to avoid over-counting, as Φ 1,k = Φ 1,n−k+2 . The couplingC Φ 1,k :T φ::T φ: may be computed exactly as in the previous subsection. We consider, for some k ∈ {2, . . . , [n/2] + 1}, the following ratio of correlators, which we evaluate using the OPEs in order to extract the structure constant: T (x 1 )T (x 2 ) In the last step, we have used the fact that, by definition, every independent bilinear in Φ 1,k occurs with coefficient 1. We can then evaluate this explicitly by conformally mapping to the complex plane: (where the power functions are on their principal branch), thus, comparing both formulae we findC We consider now the next correction to the OPE of T andT , involving the fields Φ 1,k,j with k > j > 1. Again, the ranges of k and j must be further restricted in the OPE in order not to overcount the fields. We do not need to discuss this in general; we just note that in both cases n = 3 and n = 4 there is a single field to count, Φ 1,2,3 = φ 1 φ 2 φ 3 (for n = 3) and Φ 1,2,3 = φ 1 φ 2 φ 3 + φ 2 φ 3 φ 4 + φ 3 φ 4 φ 1 + φ 4 φ 1 φ 2 (for n = 4). As usual we first consider the consequence of the OPE, We then perform the calculation of the correlation function by mapping to the sphere, n g(x 3 ))φ(e 2πik n g(x 4 ))φ(e 2πij n g(x 5 )) = n −6∆ |x 2 − x 1 | 6∆ φ(e 2πi n g(x 3 ))φ(e 2πik n g(x 4 ))φ(e 2πij n g(x 5 )) (|x 3  As before, we compute first and we may compute the same three point function by using the conformal map (138) together with the definition (11) : T φ : (x 1 ) :T φ : ( φ(e 2πij 1 n g(y 1 ))φ(e =C φ φφ lim Hence we concludeC T (x 1 )T (x 2 ) We then calculate: φ(e 2πij 1 n g(y 1 ))φ(e 2πij 2 n g(y 2 ))φ(e 2πi n g(x 3 ))φ(e 2πik n g(x 4 )) = n 4∆−2 n 1−4∆ lim whence we concludeC where κ is a model-dependent function which characterizes the four-point function E Computation of A 1 (r, 1) We have seen that A 1 (r, n) = − n π F :T φ:|1 1 : T φ : The n-dependence of this expression is contained on the one-particle form factor, so we need to compute lim n→1 d dn   n F :T φ:|1 1 : T φ : where we have used the fact that The derivative above can be computed from the integral representation (43)  Simplifying we obtain, with B 1 (r) defined in (124).