Massive Corrections to Entanglement in Minimal E8 Toda Field Theory

In this letter we study the exponentially decaying corrections to saturation of the second R\'enyi entropy of one interval of length L in minimal E8 Toda field theory. It has been known for some time that the entanglement entropy of a massive quantum field theory in 1+1 dimensions saturates to a constant value for m1 L<<1 where m1 is the mass of the lightest particle in the spectrum. Subsequently, results by Cardy, Castro-Alvaredo and Doyon have shown that there are exponentially decaying corrections to this behaviour which are characterised by Bessel functions with arguments proportional to m1 L. For the von Neumann entropy the leading correction to saturation takes the precise universal form -K0(2m1 L)/8 whereas for the R\'enyi entropies leading corrections which are proportional to K0(m1 L) are expected. Recent numerical work by P\'almai for the second R\'enyi entropy of minimal E8 Toda has found next-to-leading order corrections decaying as exp(-2m1 L) rather than the expected exp(-m1 L). In this paper we investigate the origin of this result and show that it is incorrect. An exact form factor computation of correlators of branch point twist fields reveals that the leading corrections are proportional to K0(m1 L) as expected.


Entanglement Entropy
Measures of entanglement, such as the entanglement entropy (EE), have attracted much attention in recent years, particularly in the context of one-dimensional many body quantum systems (see e.g. review articles in the special issue [1]). Among such systems, those enjoying conformal invariance in the scaling limit are of particular interest as they provide a theoretical and universal description of critical phenomena. In their seminal work Calabrese and Cardy [2] used principles of Conformal Field Theory (CFT) to study the (EE) [3] of quantum critical systems. Their results generalised previous work [4], provided theoretical support for numerical observations in critical quantum spin chains [5] and highlighted the fact that the EE encodes universal information about quantum critical points, such as the central change of the corresponding CFT. This information may be extracted numerically in a very efficient way, typically by employing Density Matrix Renormalization Group methods [6], and this has provided one of the main motivations to investigate measures of entanglement in critical and near-critical systems. From a mathematical physics viewpoint (the one taken in this paper) the investigation of the EE is driven by interest in developing a better (if possible, analytical) understanding of the universal properties of the ground state of extended many body quantum systems.
The EE is a measure of the amount of quantum entanglement, in a pure quantum state, between the degrees of freedom associated to two sets of independent observables whose union is complete on the Hilbert space. In the present paper, the two sets of observables correspond to the local observables in two complementary connected regions, A and B, of a 1+1-dimensional massive quantum field theory (QFT), and we will consider only the case where the quantum state is the ground state. Let |Ψ be such a ground state. Consider a space bi-partition of the theory as sketched in Figure 1. Then the EE associated to region A may be expressed as S( ) = −Tr(ρ A log ρ A ) where ρ A = Tr B (|Ψ Ψ|) is the reduced density matrix associated to subsystem A and is the subsystem's length.
The EE defined above is also known as von and have the property lim n→1 S n ( ) = S( ). In this paper we will consider the case n = 2 where S 2 ( ) = − log Trρ 2 A . We will refer to this quantity as the second Rényi entropy. We choose this particular value in order to compare with results obtained in [7] for the same quantity by a different method.
As mentioned earlier, at quantum critical points, the scaling limit of the EE has been widely studied in CFT [8,4,5,9,2,10] and in lattice realizations of critical systems such as quantum spin chains [11,12,13,14,15,16,17] and lattice models [18,19,20]. In particular, the combination of a geometric description, Riemann uniformization techniques and standard expressions for CFT partition functions is very fruitful. Recently [21], this was generalized to non-unitary CFT and to the EE of excited states [22,23], where general formulae were obtained using also 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. So far the most powerful way of studying the EE in massive QFT is by using an approach based on local branch-point twist fields [24,25,26]. This approach is very fruitful and complete as it allows both for numerical and analytical computations. In this context, the second Rényi entropy may be defined as: where T andT := T † are the branch point twist fields, is their conformal dimension (at criticality) [27,28,2] as a function of the central charge c, and is a non-universal short-distance cut-off. The expression T (0)T ( ) 2 above denotes the two-point function in the ground state for n = 2. An important subtlety is that branch point twist fields are local fields in a new QFT which is constructed as n non-interacting copies (in this case 2) of the original QFT. In this context, they are interpreted as symmetry fields associated to the cyclic permutation symmetry of the "replica" theory.
In this paper we aim to compare results based on a branch point twist field approach to recent numerical results by Tamás Pálmai [7] for the quantity (2). In [7] a new approach to the computation of the second Rényi entropy in massive integrable QFT was proposed which is based on the use of the Truncated Conformal Space Approach (TCSA) first proposed by Yurov and Zamolodchikov in [29]. The TCSA is based on Zamolodchikov's view of massive integrable models as massive perturbations of CFT [30]. It exploits the rich structure of the Hilbert space of CFT, perturbes and truncates the latter and then diagonalizes the "truncated" Hamiltonian. This provides a very successful way to access the low energy spectrum of massive integrable QFT with (a priori) any desired level of accuracy. The work [7] showed for the first time that TCSA can also be employed to access the quantity Trρ 2 A , and so may be applied to the study of measures of entanglement. This is a very interesting development which complements and enhances the existing twist field approach for massive QFTs.

The Model
In this paper we consider an integrable massive QFT sometimes referred to as the critical Ising model in a magnetic field (IMMF) and also known as the minimal E 8 Toda field theory. In the spirit of Zamolodchikov's work [30], the theory can be described as a massive perturbation of the conformal Ising model whose operator content consists of simply three fields: the identity, the energy field ε and the spin field σ. It is well-known that a massive perturbation by the energy field gives rise to an integrable QFT known as the massive Ising model. This theory has a single particle and the two-body scattering matrix is simply S(θ) = −1 as a function of the rapidity θ. Surprisingly, perturbing with the spin field σ instead gives rise to a much more complex but still integrable interacting QFT, the IMMF [30,31]. The theory consists of 8 self-conjugate particles of different masses. All of the particles can also be formed as bound states of two other particles in the spectrum, that is, the corresponding two-body scattering matrices have a rich pole structure in the physical sheet with poles of up to order 12. Following Zamolodchikov's work, a plethora of papers by many authors led to the realization that the IMMF is but a particular case of a much wider family of integrable QFTs known as minimal Toda field theories (a detailed historical account of these findings can be found in [32,35] and references therein). These in turn are "simplified" versions of another class of models, the Affine Toda field theories (ATFTs), in the sense that the S-matrices of ATFTs are equal those of minimal Toda theories, up to coupling-dependent multiplicative factors which have no poles in the physical sheet. ATFTs have been studied since a long time and have played a prominent role in the development of the field of integrable field theories [36,37]. Based on the IMMF example and on extensive work on classical Toda theory it came to be expected that a different theory should exist for each simple Lie algebra and that this Lie algebraic structure would be crucial in the understanding of these models. Subsequently, a lot of work was carried out in order to compute the S-matrices of ATFTs related to each simple Lie algebra. Universal formulae for all simply-laced ATFTs were given in [32]. A universal description of the S-matrices of non-simply-laced ATFT based on q-deformed Coxeter elements was first proposed in [34] and further studied in [35] where a universal integral representation for all ATFT S-matrices was given. In this context, the eight particles in the spectrum of the IMMF are in one-to-one correspondence with the simple roots of E 8 (see Fig. 2). Indeed, their values are the Perron-Frobenius eigenvector of the corresponding Cartan matrix. A detailed account of the masses and scattering matrices of the theory can be found for instance in the review [38] and also in [39]. Here we will only report the data that we need for the present paper. We will only require the values of the masses of the four lightest particles in the spectrum where m 1 is the mass of the lightest particle. We note that the masses of the IMMF satisfy the inequality m i > 2m 1 for i > 3 whereas m 1,2,3 < 2m 1 . In the following we will also require some of the two-particle scattering amplitudes as functions of the rapidity variable θ. They generally have the structure: (5) has four simple poles at 4πi 5 , 3πi 5 , 7iπ 15 and 4iπ 15 corresponding to the formation of bound states 1, 2, 3 and 4. Associated to these simple poles are the three-point couplings Γ c ab defined as i(Γ c ab ) 2 = Res θ=iπα S ab (θ).
In particular will enter in some of the equations we will see later.

Structure of the Paper
This paper is organised as follows: in section 2 we review the form factor approach for branch point twist fields. We propose expressions for some of the two-particle form factors in the IMMF as well as a set of consistency conditions that allow us to determine also the one-particle form factors of the four lightest particles in the spectrum. In section 3 we explain how the second Rényi entropy may be expressed in terms of twist field form factors and write down an expression including the six leading form factor corrections to its saturation value. We determine the precise coefficients of these corrections by solving the form factor equations proposed in section 2. We confirm the presence of exponentially decaying corrections, led by e −m 1 . In section 4 we compare our results to those obtained in [7] by employing the TCSA approach and discuss their level of agreement. We present our conclusions and outlook in section 5.

Generalities
In 1+1 dimensional integrable QFT the most successful approach to computing multi-point functions of local operators is by expressing them in terms of form factors of individual fields.
Form factors of local fields in integrable QFT can generally be computed exactly by pursuing the so-called form factor programme [40,41] which was extended to the treatment of branch point twist fields in [24]. The programme has been carried out for countless models and fields and provides extremely accurate results for correlators, particularly two-point functions. Here we are interested in the correlator (2) in the IMMF. Let (a, j) represent the quantum numbers of a particle of type a = 1, . . . , 8 living in copy j = 1, . . . , n (we will specialize to n = 2 later on). We may employ the so-called cumulant expansion [42,43,44]: with where the functions h a 1 ...a k j 1 ...j k (θ 1 , · · · , θ k ; n) are given in terms of the form factors of the fields involved. In particular, are the normalized one-and two-particle form factors. Here 0| represents the vacuum state and |θ (a,j) , |θ 1 θ 2 (a 1 ,j 1 )(a 2 ,j 2 ) represent in-states of 1 and 2 particles, respectively. The states are characterized by the rapidities θ i and particle quantum numbers (a i , j i ). In this paper we will only be interested in one and two particle form factor contributions which provide the most important contributions to (11) for large . For spinless operators in relativistic theories we have that the one-particle form factors are rapidity independent. Therefore, from now on we will simply write F T |(a,j) 1

One-and Two-Particle form Factors of the IMMF
Form factors of the IMMF where constructed in [39] for the stress-energy tensor (e.g. the spin field) and in [45] for the energy field. This construction can be easily adapted to the twist field by employing the programme proposed in [24]. In addition, here we only want to study a subset of the one-and two-particle form factors of the twist field which means that we do not need to engage into solving the complicated recursive equations that arise for higher particle numbers. We will also avoid the consideration of higher order poles which has been extensively discussed in [38,39]. A particular feature of the computation of form factors of twist fields is that many of the basic formulae are very similar to those used in the construction of form factors of standard local fields, particularly when it comes to constructing the two-particle minimal form factor and therefore, what follows is strongly guided by the analysis of [39]. A minimal solution to the two-particle form factor equation will be denoted by This minimal form factor satisfies a twist field version of Watson's equations which may be summarised as [24] f a 1 a 2 (θ, Let S a 1 a 2 (θ) have an integral representation of the form where s a 1 a 2 (θ) is a function which depends of the theory and the particles a 1 , a 2 . Employing this integral representation it is easy to show that the solution to (15) may be written as It is also well-known that if S a 1 a 2 (θ) = −1 in (15) then the free Fermion solution f a 1 a 2 (θ) = −i sinh θ 2n is obtained. Combining these two results and comparing with the formulae provided in [39] for the minimal form factors of the IMMF we find where the exponential may be also expressed as the infinite product of gamma functions: and the set S a 1 a 2 has been defined after (5). Following [24] and [39], the most generic two-particle form factor takes the form with θ 12 := θ 1 − θ 2 , and u(2k The function K j 1 j 2 (θ; n) encodes the full kinematic pole structure of the form factor, having kinematic poles at θ = iπ and θ = iπ(2n − 1) in the extended physical strip Im(θ) ∈ [0, 2πn], whereas B α (θ; n) encodes the bound state pole structure, as characterised in [39]. The minimal form factors in (20) can be easily obtained from (18) by employing standard relations which can be found for instance in [24]. Finally, the functions Q j 1 j 2 a 1 a 2 (θ; n) are solutions to the equations: namely, they are linear combinations of functions of the type cosh k θ n for k = 1, 2, . . ., similar to the ansatz already employed in [39]. Which values of k are involved will be determined by constraints on how Q 11 a 1 a 2 (θ; n) behaves as θ → ∞ which we will discuss in the next subsection. In this paper we will only require the two-particle form factors F (1,1)(1,1) 2 (θ 1 , θ 2 ; n) and F (1,1)(2,1) 2 (θ 1 , θ 2 ; n). These are special cases of (20) explicitly given by and F T |(1,1)(2,1) 2 (θ 1 , θ 2 ; n) = Q 11 12 (θ 12 ; n) 2n with

Fixing One-and Two-Particle Form Factors
The equations (25)-(26) give the two-particle form factors of interest up to the functions Q 11 11 (θ; n) and Q 11 12 (θ; n). As anticipated earlier, these functions may be determined by employing additional constraints. In particular, the kinematic and bound state residue equations for the two-particle form factors require that: and lim In the IMMF all particles are self conjugate so that the form factor (25) must satisfy the condition (29), giving the constraint: 9,14,16,21,25,30 sin πα 30n .
The same form factor possesses three bound state poles related to the formation of bound states, 1, 2 and 3. This means that it satisfies three versions of equation (30), giving three additional constraints: 5,9,11,16,20,25 sin απ 30n , The form factor (26) has no kinematic poles but has four bound state poles associated to the formation of bound states 1, 2, 3 and 4. They give the additional constraints:    (n) with a = 1, 2, 3, 4 as well as the polynomials Q 11 11 (θ; n) and Q 11 12 (θ; n). We know from the two particle form factor equations that they must be even functions of θ and we would also like to require the cluster decomposition property which has been discussed in detail in [49] and observed for numerous models (a particularly rich example can be found in [50]), namely that lim and lim These properties provide very strong constraints for the functions Q 11 11 (θ; n) and Q 11 12 (θ; n), as they allow us to determine the highest power of e θ/n that can be involved. We have that It is also easy to determine the leading behaviours of f 11 (θ; n) and f 12 (θ; n) as θ → ∞. In this limit, integrals of the type may be approximated by changing variables to x = tθ and then expanding for small values of x θ . At leading and next-to-leading order, this yields the simple integral: Numerical evaluation of the minimal form factors for θ large confirms the behaviour above, up to n-dependent proportionality constants which we, unfortunately, have been unable to find a convergent analytic expression for Nonetheless, the numbers v(n) and u(n) can be numerically estimated for every n. This means that, in order to satisfy (39)-(40) we need thus, in general Q 11 11 (θ; n) = A 11 11 (n) + B 11 11 (n) cosh θ n + C 11 11 (n) cosh 2 θ n , and Q 11 12 (θ; n) = A 11 12 (n) + B 11 12 (n) cosh θ n + C 11 12 (n) cosh 2 θ n .
And the cluster decomposition equations (39)- (40) can be written as It is worth noting that the same conclusions regarding the form of the functions Q 11 11 (θ; n) and Q 11 12 (θ; n) can be reached by appealing to a well-known argument presented for instance in [39] according to which the form factors of unitary operators (in particular, the twist field) can diverge at most as e ∆nθ i when one of the rapidities θ i → ∞ and where ∆ n is the conformal dimension (3) of the twist field. In fact, for c = 1 2 it turns out that ∆ n < 1 for n < 49 and therefore, at least for a wide range of values of n the form factors must tend to a constant as any of the rapidities they depend upon tends to infinity. The property of cluster decomposition, additionally establishes what this constant must be.
Putting together equations (31)- (38) and (48) we end up with 10 equations for 10 unknowns: A 11 11 (n), B 11 11 (n), C 11 11 (n), A 11 12 (n), B 11 12 (n), C 11 12 (n) and the four one particle form factors F (n) and F T |(4,1) 1 (n). This means we are now in a position to determine all these functions and to investigate how our results apply to the study of the second Rényi entropy.

Second Rényi Entropy of the IMMF
From the definition (2) and the expansion (10) we know that the Rényi entropy may be expressed as an infinite sum involving integrals over the form factors of the twist field. In particular the first few contributions can be written as and S 2 (∞) = −4∆ 2 log(m ) − 2 log T 2 is the non-universal saturation constant. The expansion above describes corrections to saturation which are exponentially decaying for large m 1 . If we consider all terms above, it becomes quickly apparent that some of the one-particle form factor contributions are subleading compared to some of the two-particle form factor contributions, for large. This is because, as mentioned earlier in the paper, the masses of particles 4, . . . , 8 are all larger than twice the mass of particle 1. In summary, this means that, the first six leading form-factor corrections to saturation in order of importance are and we have carried out one of the integrals in (49). For m 1 1 the leading contribution to the first integral in (51) can be written as: and similarly for the last integral where the remaining integrals are convergent and can be easily evaluated numerically.

Numerical Computation of the One-Particle Form Factors
Although we are now in a position to solve equations (31)- (38) and (39)- (40) and therefore obtain explicit formulae for the one-and two-particle form factors of interest, in practice these equations are rather cumbersome and finding exact formulae is extremely difficult (even for n = 2). They can however be very easily solved numerically for any given value of n. An interesting observation from solving the equations numerically is that the solution is not unique. This is mainly due to the equations (39)- (40) which are quadratic in the one-particle form factors. There are in fact three solutions for each value of n. This obviously poses the question as to which of these solutions is the correct one. It also indicates that the form factor  (51)).
equations allow for several twist field solutions. This is not surprising as all operators enjoying the twist-property should have form factors satisfying the same set of equations. Among those operators there is the branch point twist field we are interested in, but also other twist fields, such as the composite twist fields introduced in [46] and studied in [47,21,48]. These are fields which are defined (at criticality) as the leading field in the OPE of the standard branch point twist field and any other local field of the replica theory.
Fortunately, there is one simple way of telling such composite twist fields and the twist fields we are interested in apart: composite twist fields have form factors which do not vanish at n = 1 whereas the twist field T must reduce to the identity field at n = 1, hence all its form factors vanish at n = 1. Imposing this condition we find a single solution with the desired property of having vanishing form factors at n = 1. In addition, there are certain consistency checks that we may further apply to test this solution. A common test is the ∆-sum rule proposed in [49] which may be written as where Θ(r) is the energy-momentum tensor. Employing a standard form factor expansion, it is then possible to recover the conformal dimension of the twist field (3) from its two-point function with the energy-momentum tensor. This computation is possible thanks once more to the results of [39] where the form factors of the energy-momentum tensor were obtained, in particular all one-and two-particle form factors. Fig. 3 shows the numerically obtained values of ∆ n employing (55). As can be seen it is possible to compute these values also for non-integer n as the form factors as well-defined for all values of n. It is also noticeable that the saturation of the ∆-sum rule becomes worse the larger n is. This is a rather common feature of the branch point twist field which is due to the fact that all form factor contributions to the expansion are proportional to n and so the weight of further form factor corrections is increased as n increases. We have also observed that the numerical determination of the constants u(n) and  Figure 4: One particle form factors associated to particle 1 (circles), particle 2 (squares), particle 3 (triangles) and particle 4 (rombi).
v(n) in the asymptotics (44) becomes more difficult for larger n. Both constants are the result of the numerical integration of a decaying but wildly oscillating function which is delicate. A less precise knowledge of the values u(n) and v(n) may well also contribute to the worsening agreement observed in Fig. 3. Fortunately though we have very good agreement for n = 2 which both confirms the one-particle form factors are correct and shows that including the six first contributions to the form factor expansion leads to nearly exact results. The first four one-particle form factors F T |(a,1) 1 (n) for a = 1, 2, 3, 4 are all real numbers and their values are presented in Fig. 4. All values are rather small and, in absolute value, are smaller the higher the particle label. In particular: (2) = 0.0230515.
We can also numerically determine the values of the coefficients A 1a 11 (n), B 11 1a (n) and C 11 1a (n) for a = 1, 2. For n = 2 they are:

Exact Corrections to Saturation
Putting together all the numerical values found in the previous section, we see that the formula (51) may be expressed as: √ m 1 , respectively. From (59) it is clear that the one-particle form factor contributions are led by small coefficients (as the one-particle form factors are all small numbers). However, the two-particle contributions are characterized by even smaller coefficients and are in fact strongly suppressed (they are essentially negligible for m 1 above 1), as shown in Fig. 5. Therefore the suggestion The dashed curve is the contribution −0.0182441K 0 (m 1 ), that is the leading contribution to (51). Clearly the contribution proportional to K 0 (m 1 ) is leading for the full range of values of considered here. put forward in [7] that the one-particle form factors are zero and therefore not detectable in TCSA is mistaken. In the next section we will review how this suggestion was arrived at in [7] and clarify the origin of this apparent mismatch between TCSA and form factor results.

Comparison to TCSA Results
A good way to gain at least a qualitative understanding of how the TCSA results and the form factor results compare is to display them in the same graph. This is what is shown in Fig. 6. More precisely, the next-to-leading order corrections to the second Rényi entropy obtained from TCSA and from form factors are presented. Let us describe the data in some detail: • The black solid line represents such corrections as are obtained from an exact form factor calculation which includes the first six terms of the form factor expansion as shown in equation (51). As shown in Fig. 5 out of these six contributions, the first one is very dominant, with subsequent contributions being strongly suppressed. This suggests that the solid curve in Fig. 6 provides a nearly exact description of the corrections to saturation of the entropy of a subsystem of size m 1 in an infinite system.
• The numerical data in the Fig. 6 represented by squares, circles and rombi are the TCSA results for finite systems of lengths m 1 = 6, 7 and 8, respectively. These are the same data as employed in [7] but where the saturation value has been subtracted 1 . Since the TCSA approach (by construction) can only deal with finite systems, the expectation is that a match with a form factor computation will only be achieved for large system sizes (how large will typically depend on the quantity that is being computed). Observing Fig. 6 we can conclude that agreement with the form factor data is indeed better as volume increases (as expected) and is already very good for the data corresponding to m 1 = 8.
• The final element of Fig. 6 is the dashed curve which depicts the function This function is the result of fitting the numerical TCSA data corresponding to m 1 = 6 with a decaying exponential 2 . A similar fit of the m 1 = 5 data was used in [7] to deduce the rate of decay of the exponential corrections to saturation. This fit plays an important role as it is from this single evidence that the author of [7] concludes there are no corrections to entanglement of the form e −m 1 and therefore the one-particle form factors must be zero.
From Fig. 6 and the points above it would seem that two contradictory conclusions follow: the TCSA data and the form factor data agree rather well, yet the exponential decay identified in [7] disagrees with the form factor calculation. Where does this mismatch come from?
This question is easy to answer and the answer lies in the manner in which the exponential decay of the TCSA data was (wrongly) identified. This wrong identification is due to the fact that the fit (60) is a good fit of the m 1 = 6 data but clearly not a good fit of the larger system size data. The correct manner of identifying the exponential decay from TCSA would have been to first find a reasonable infinite volume extrapolation of the data and then find a fit of that extrapolation. We would expect such a procedure to give a result much more in tune with the form factor analysis.

Conclusions
In this paper we have carried out an in depth study of the second Rényi entropy in the transverse field Ising model or E 8 minimal Toda field theory employing the branch point twist field approach developed in [24]. This work, follows on from a stream of works where the entanglement of various particular integrable models has been studied [24,51,25,52,53,54] by the same method. Our results provide the fist detailed form factor study of the entanglement of an interacting theory with a complicated particle content and no internal symmetries so that all twist field form factors are non-vanishing. We have compared our results to those obtained by a new method for the evaluation of measures of entanglement developed by T. Pálmai [7]. This new method is based on the use of the truncated conformal space approach [29] to massive integrable quantum field theories.
One of the motivations for this work was to check the veracity of the claim that the leading corrections to saturation of the second Rényi entropy are not those predicted by the twist field approach. This claim was put forward in [7] based on the analysis of some TCSA data. In this paper we have shown that the leading corrections to saturation are indeed those expected from a twist field form factor analysis and that there was instead a fault in the interpretation of the TCSA data. The main conclusion about the exponential decay of entanglement had been reached based on data for relatively small system sizes and, not surprisingly, these data did not reproduce correctly the infinite size behaviour described by the form factor approach. Once this point is clarified there is in fact no contraction between the form factor and TCSA approaches and, based on the data at our disposal, we can say that they agree reasonably well already for systems sizes of the order of m 1 = 8.
We conclude from this analysis that it is actually rather difficult to correctly identify subleading exponential corrections to entanglement purely from a TCSA analysis, as this requires the ability to extrapolate numerical data to the infinite size limit with great precision. In cases where this is possible, TCSA is a good numerical alternative to the form factor approach which, although leading to exact results, can be hard to use in complicated models. The lengthy computations needed to arrive at (51) illustrate this point rather well.