Branch Point Twist Field Correlators in the Massive Free Boson Theory

Well-known measures of entanglement in one-dimensional many body quantum systems, such as the entanglement entropy and the logarithmic negativity, may be expressed in terms of the correlation functions of local fields known as branch point twist fields in a replica quantum field theory. In this"replica"approach the computation of measures of entanglement generally involves a mathematically non-trivial analytic continuation in the number of replicas. In this paper we consider two-point functions of twist fields and their analytic continuation in the 1+1 dimensional massive (non-compactified) free Boson theory. This is one of the few theories for which all matrix elements of twist fields are known so that we may hope to compute correlation functions very precisely. We study two particular two-point functions which are related to the logarithmic negativity of semi-infinite disjoint intervals and to the entanglement entropy of one interval. We show that our prescription for the analytic continuation yields results which are in full agreement with conformal field theory predictions in the short-distance limit. We provide numerical estimates of universal quantities and their ratios, both in the massless (twist field structure constants) and the massive (expectation values of twist fields) theory. We find that particular ratios are given by divergent form factor expansions. We propose such divergences stem from the presence of logarithmic factors in addition to the expected power-law behaviour of two-point functions at short-distances. Surprisingly, at criticality these corrections give rise to a log(logL) correction to the entanglement entropy of one interval of length L. This hitherto overlooked result is in agreement with results by Calabrese, Cardy and Tonni and has been independently derived by Blondeau-Fournier and Doyon (in preparation).


Introduction
The problem of quantifying the amount of entanglement which may be "stored" in the ground state of a many body quantum system has attracted the interest of the quantum information and theoretical physics communities for a long time. Measuring entanglement is of interest both if we are to employ entanglement as a quantum computing resource and if we want to learn more about the fundamental features of quantum states of highly complex quantum systems. Among such systems, 1+1-dimensional many body quantum systems have received considerable attention over the past decade. Much work in this area has been inspired by the results of Calabrese and Cardy [1] which used principles of Conformal Field Theory (CFT) to study a particular measure of entanglement, the entanglement entropy (EE) [2]. In this seminal work, they generalised previous results [3] and provided theoretical support for numerical observations in critical quantum spin chains [4]. Before we proceed any further a few definitions are in order: let | i be a pure state describing the ground state of quantum spin chain at zero temperature. Consider a bi-partition of the chain such as in Fig. 1(a) (suppose there are periodic boundary conditions). Then the entanglement entropy associated to region A may be expressed as S(`) = Tr(⇢ A log ⇢ A ) where ⇢ A = Tr B (| ih |) is the reduced density matrix associated to subsystem A and`is the subsystem's length.
One of the main results of [3,4,1] describes the entanglement entropy of 1+1 dimensional many body quantum systems (such as spin chains) in the continuous limit at criticality. Such systems are described by CFT and their EE displays universal features expressed by the now famous formula: S(`) = c 3 log✏ . That is, the EE of a subsystem of length`of an infinite critical system diverges logarithmically with the size of the subsystem, with a universal coe cient which is proportional to the central charge of the CFT, c. There are non-universal constant corrections to this leading behaviour which may be encoded by a short-distance cut-o↵ ✏. This behaviour has been numerically and analytically studied for a plethora of spin chain models in works such as [4,5,6,7,8,9,10,11,12,13]. Another popular measure of entanglement is the logarithmic negativity (LN) [14,15,16,17,18]. Consider again a quantum spin chain in a pure state | i and a partition such as depicted in Fig. 1(b). Then, the LN is a measure of the amount of entanglement between the two non-complementary sub-systems A and B. Its formal definition depends on the reduced density matrix ⇢ A[B as E(`1,`2,`3) = log(Tr|⇢ T B A[B |) where T B represents partial transposition with respect to subsystem B and |⇢| is the trace norm of ⇢, that is the sum of the absolute values of its eigenvalues. The LN of 1+1 dimensional critical systems has been studied numerically in [21,22,23] and more recently, both numerically and analytically exploiting fundamental conformal field theory principles, in [19,20]. Since then many particular models have been analysed (see e.g. [26,27,28,29]). However, for general configurations such as in Fig. 1(b) there is no known analytic formula for generic CFTs. There are however particular limits which are easier to treat such as the limit of adjoint intervals (`2 ! 0) and the limit of semi-infinite disjoint intervals (`1,`3 ! 1 keeping`2 finite). The former has been studied in [19,20] for generic CFT yielding the simple expression E(`1, 0,`3) = c 4 log`1`3 ✏(`1+`3) whereas the latter is harder to treat in critical systems but is of interest in the study of quantum systems near criticality. Such systems are described by 1+1 dimensional massive quantum field theories which, unlike CFT, allow for the existence of a finite correlation length. The negativity of such systems was first studied in [30] where new results for both of the limits above in near-critical systems were obtained.
In this paper we will be interested in a particular prescription for the calculation of both the EE of a single interval and the LN of semi-infinite disjoint regions. It turns out that both quantities may be expressed in terms of two-point functions of a particular class of fields known as branch point twist fields [1,32]. This relationship comes about through a technique commonly known as the "replica trick". The replica trick may be applied to both the computation of the EE and of the LN. It involves a rewriting of the definitions above as follows where the symbol n e in the second formula means n even, that is the limit n to 1 must be carried out by analytically continuing the function from even, positive values of n to n = 1. The representations above were used first in [3] for the EE and in [19,20] for the LN. The advantage of such representations is that both Tr(⇢ n A ) and Tr(⇢ T B A[B ) ne admit a natural physical interpretation as partition functions in "replica" theories. The replica theory is a new model consisting of n non-interacting copies of the original theory. In this context it is natural for n to take positive integer values. However, the definitions (1) require that such traces be analytically continued from n integer (and in the LN case, also even) to n real and positive. Hence, formulae (1) are advantageous in that partition functions in replica theories may be computed systematically by various approaches, but also disadvantageous because the analytic continuations involved are often very di cult to perform and there is no generic proof of existence and uniqueness.
It was first noted in [1] that the function Tr(⇢ n A ) may be expressed as a two-point function of fields with conformal dimension given by In fact such fields had been previously discussed in the context of the study of orbifold CFT where they emerge naturally as symmetry fields associated to the permutation symmetry of the theory [33,34]. In [32] such fields were named branch point twist fields and studied in the context of 1+1 dimensional massive QFT. Their connection to the cyclic permutation symmetry of the replica theory was made explicit by formulating their exchange relations with the fundamental fields of a generic replica QFT. For integrable QFT this allowed for the formulation of twist field form factor equations whose solutions are matrix elements of twist fields. Let T be a twist field associated to the cyclic permutation symmetry j 7 ! j + 1 andT its conjugate, associated with the permutation symmetry j 7 ! j 1 with j = 1, . . . , n. Then, we may write: Tr(⇢ n A ) = ✏ 4 n hT (0)T (`)i n and Tr(⇢ T B A[B ) n = ✏ 8 n hT ( `1)T (0)T (`2)T (`2 +`3)i n .
(3) At criticality, these formulae may be used directly to derive the expressions for S(`) and E(`1, 0,`2) given above. The same formulae may be used to study QFT beyond criticality as done in [32,30]. In this paper we will analyse the short-distance (e.g.`⌧ 1) behaviour of the correlators hT (0)T (`)i n and hT (0)T (`)i n = hT (0)T (`)i n in a massive free Boson theory. At short-distances we expect the massive QFT to be described by its corresponding ultraviolet limit (that is, the massless (non-compactified) free Boson CFT). Thus, we expect these two-point functions to exhibit power-law behaviours with powers related to the dimension of twist fields. Extracting these short-distance behaviours from a form factor expansion (which is eminently a large-distance expansion) is generally highly non-trivial and can seldom be done with precision for any fields. However, as we will see, this can be done with great precision for the massive free Boson, on account of the theory's simplicity and the special properties of the twist field form factors. For the massive free Boson all form factors of twist fields, that is objects such as are known explicitly. Here h0| represents the vacuum state and |✓ 1 , · · · , ✓ k i j 1 ...j k represents an in-state of k particles with rapidities ✓ 1 , . . . , ✓ k and quantum numbers j 1 . . . j k . In the free Boson case, these quantum numbers are just the copy number of the Boson in the replica theory. Here we have chosen to normalise all form factors by a constant (the vacuum expectation value of the twist field hT i n ). This will be convenient for later computations. By reconstructing the short-distance (power-law) behaviour of the correlators hT (0)T (`)i n and hT (0)T (`)i n for n 1, integer or not, we will provide strong evidence for our approach to performing the analytic continuation of the correlators in n. This will provide support for our methodology and will allow us to examine twist field two-point functions further and extract values of universal quantities such as expectation values and structure constants of twist fields.
The paper is organized as follows: In sections 2 and 3 we review basic CFT and QFT results, regarding the short distance behaviour of two-point functions of twist fields and how these twopoint functions may be expressed in terms of the form factors (4). In section 4 we show how the power-law decay of two-point functions of twist fields may be obtained exactly from form factors in the massive free Boson theory for n 1 real. In section 5 we provide form factor expansions for the constant (universal) coe cients that multiply the leading power-law in the two-point functions of twist fields. We employ these expansions to obtain numerical predictions for the ratio of the structure constant C T 2 T T and the expectation value hT i n , analytically continued from n odd and for the structure constant C T 2 T T analytically continued from n even. We compare our values of C T 2 T T for n even to analytical values obtained in [20] and find good agreement. We numerically examine the limit lim ne!1 + C T 2 T T and compare to an analytical prediction given in [20]. In section 6 we present an interpretation of the emergence of divergent sums in the representation of particular ratios of expectation values and structure constants of the massive free Boson theory. We propose that such divergences must be related to the presence of logarithmic corrections to the two-point functions at criticality. We conclude that such corrections will give rise to an additional log(log`) term in the EE and the Rényi entropy of one interval in the massless (noncompactified) free Boson theory. This is in full agreement with previous results for the LN [20] and the EE [24] of the compactified massless free Boson in the limit of infinite compactification ratio. For the EE the presence of such corrections has also been established analytically by a di↵erent method in [25] but had been overlooked in [31]. In section 7 we compare our numerical estimates of the value of lim ne!1 + C T 2 T T as well as the analytical value given in [20] to a value that can be read o↵ from numerical results in [35] for the LN of a harmonic chain out of equilibrium and their CFT interpretation [36]. We present our conclusions in section 8. Appendix A collects some useful summation formulae which feature in the form factor expansions of sections 4 and 5. Appendix B provides a discussion and assessment of the error of some of our numerical procedures.

Conformal Field Theory Recap
As described in the introduction, we wish to study the two-point functions hT (0)T (`)i n and hT (0)T (`)i n and examine their short-distance behaviour. This behaviour is entirely predicted by CFT and may be expressed as Similarly for n even (6) Note that by examining the next-to-leading order (`-independent) corrections above we may extract values for universal QFT quantities such as the twist field expectation value hT i n and the structure constants C T 2 T T and their ratios. These are di cult to compute by other methods, demonstrating once more that the form factor approach in particularly powerful in this context.
The di↵erence between the n odd and n even cases was first discussed in [19,20] and follows from the leading term in the conformal OPE of the field T with itself, which takes the form This leading term is characterized by a new twist field T 2 of conformal dimension T 2 which is associated with the permutation symmetry j 7 ! j + 2 for j = 1, . . . , n. As discussed in [19,20] the nature of this field is very di↵erent depending on whether n is odd or even. Whereas for n odd, the field T 2 is equivalent to the field T (the permutation j 7 ! j + 2 still allows for visiting all copies, albeit in a di↵erent order), for n even the permutation j 7 ! j + 2 divides even-and odd-labeled copies so that T 2 is equivalent to two copies of T acting on a n 2 -replica theory. Consequently the conformal dimension of T 2 is T 2 = n for n odd and T 2 = 2 n 2 for n even. For the same reasons hT 2 i n = hT i n for n odd and hT 2 i n = hT i 2 n 2 for n even. This simple interpretation also shows how the analytic continuations (1) from n even and n odd should be di↵erent. Note that, hT i 1 = 1 both for massive and massless theories as the twist field becomes the identity field at n = 1.
In massive theories, the correlator hT (0)T (`)i n encodes the`-dependent part of the negativity E(1,`, 1) of semi-infinite disjoint regions. This follows simply from the definition (3) and the factorization of correlation functions at large distances in massive QFT.
In this paper we will use a form factor expansion of these correlators to extract the leading term (the log`term). We will turn our attending to the next-to-leading order corrections in section 5.

Form Factor Expansion of two-Point Functions
In a massive integrable QFT such as the massive free Boson, the functions (5)-(6) admit a natural large m`expansion in terms of form factors. In general we have that the (normalized) logarithm of the two-point function of local fields O 1 , O 2 admits and expansion of the form log ✓ hO with where the functions h are given in terms of the form factors of the fields involved, N is the number of particles in the spectrum and p i represent the particle's quantum numbers. For example: and so on. Here we have used the generic property: The expansion (9) with (10) is usually referred to as the cumulant expansion of the two-point function (see e.g. [37,38,39]) and it is particularly well suited to extract the leading logb ehaviours seen earlier. If all form factors are know, this may be done by employing the fact that the operators O 1 , O 2 are spinless (this will be the case for twist fields) and thus relativistic invariance implies that all form factors depend only on rapidity di↵erences. In other words, one of the rapidities in the integrals (9) may be integrated over leading to where K 0 (x) is a Bessel function and Provided the functions h p 1 ...p j j (0, ✓ 2 , · · · , ✓ j ) vanish for large ✓ we may, for m`⌧ 1 expand the Bessel function as K 0 (m`d j ) = log` + log 2 log(md j ) + · · · where = 0.5772157... is the Euler-Mascheroni constant. For m`⌧ 1 we expect the behaviour log ✓ hO then, considering the leading term in the Bessel function expansion and summing the resulting series from (8) we have that In addition, the next-to-leading correction for small m`can also be obtained as shown in [39] and is given by

Form Factors in the Massive Free Boson Theory
It is now easy to adapt the definitions above to the two-point functions of interest. In our case we are considering a free Boson theory in a replica theory, so the particle number is N = n, where n is the number of replicas. The form factors of free Boson twist fields were first reported in [30] and they can be expressed in terms of the two-particle form factor For simplicity we will from now on call Form factors associated to other copy numbers can be simply obtained by employing the properties: A direct consequence of these properties is that for the free Boson FT ( ✓; n) as the scattering matrix is 1. A detailed derivation of (19) may be found in [32,40]. Similar properties can also be derived for higher particle form factors, so that every form factor ofT may be ultimately expressed in terms of form factors of T involving only particles in copy 1 of the theory [40]. In addition, due to the Z 2 symmetry of the free Boson Lagrangian, there are only non-vanishing even-particle form factors. Higher even-particle form factors may be simply obtained by employing Wick's theorem. In general they are given by [30] where S 2j represents the set of all permutations of {1, . . . , 2j} and ✓ ij := ✓ i ✓ j (a function with this combinatorial structure is know as a permanent in mathematics). For example: This formula can be easily generalised to generic particles (e.g. particles living in di↵erent replicas) by using the relations (19).

Form Factor Expansions in the Massive Free Boson Theory
Following the definitions above, let us write with and We now have all the formulae necessary to write down the functions c T T 2j (`; n) and c TT 2j (`; n) corresponding to (9) for the correlators (5)- (6). The simple structure of the form factors (20) coupled with the nature of the cumulant expansion (8)  (✓ 1 , . . . , ✓ 2j ) which are unique for free theories and have already been observed in previous work for the massive free Fermion [41]. (10), the first term contributing to each function h

As in the examples
where we uses the fact that T =T † , that is equation (11) and the definition (20). Employing the same equations, the first term contributing to the function h where the second equality follows from generalising equations (19) to higher particle form fac- j! terms. Therefore, their product will generate a sum of ( (2j)! 2 j j! ) 2 terms. However, many of these terms are identical up to integration in all rapidities. For example, the sum (25) for j = 2 is where the symbol = int means equality under integration in all rapidities. Employing the properties (19) this may be written as where ✓ p := ✓ + 2⇡ip. Finally, for the free Boson we also have that f (x p ; n) ⇤ = f (( x) p ; n). In general, it is easy to show that there are exactly (2j 1)! terms (identical under integration) which are so-called "fully-connected". In the j = 2 example above there are exactly 6 such terms, those in the first line of (27). Including the sum over all indices p i , these are terms of the form In (29) one sum has been carried out by simply setting p 2j = 1 multiplying by a factor n (since all copies are identical) and shifting all p 0 i s by 1. The crucial observation is that all terms which are not of this form (that is, they factorise into separate multiple sums such as the terms in the last line of (27)) are cancelled in the cumulant expansions (10). They generate precisely the products of h-functions on the r.h.s. of each definition. In summary, combining (25)-(26) with the properties (19) and employing the symmetry properties induced by the integrals in (23) and (24), we find that By entirely similar arguments it can be shown that The sums in p i that enter (30)-(31) can be computed exactly for the free Boson and they are given by the formula (113) in the appendix. This will allow us to easily analyse the short-distance behaviour of correlators, with the help of formulae (15)- (16). Let us consider each two-point function separately.  [41] where the free Fermion was considered. However, in [41] some of the computations were only presented in an appendix with limited detail thus we believe it instructive to revisit the steps involved.
Consider the expression (30) and employ the formula (113) to perform the sum over the p 0 i s. According to (113) the resulting function will depend on the sum of all rapidity dependencies of the functions involved, that is It is convenient to change variables as we have also that so that we may obtain the equivalent of (23) in terms on the new variables x i and obtain (15) and (16) by integrating over the variable x 2j . Interestingly, under this change of variables, the sum which involves only the odd-indexed variables and the di↵erence This means that the leading small`contribution to the function (30) after changing variables and integrating x 2j takes the expected form x TT log`with where The integral above may be factorised into two functions depending only on even-and odd-indexed variables, respectively. This may be achieved by introducing the new variable y = P j p=1 x 2p 1 (and eliminating the variable x 2j 1 ). In terms of this variable we may rewrite some of the cosh functions in the denominator as follows: With this change of variables we find that the integral (36) becomes It was shown in [41] that these integrals can be performed exactly giving Thus, the sum (40) may be written simply as Note that the integral representation (41) only strictly makes sense for j > 1, although the formulae (42) and (43) are valid for j 1 and indeed reproduce the original integral (36) for j = 1 and G 1 (y) = sech y 2 . Although (42) and (43) were already used in [41] it is worth briefly recalling how they follow from (41) and from each other. Equation (42) can be easily derived by computing the Fourier transform in the variable y of G j (y) from (41). Although (41) is a complicated expression, by Fourier transforming in y and then changing variables to ±y ! ±y P j 1 p=1 x p all j integrals readily factorize into Fourier transforms of the same function and one obtains from where (42) directly follows. This representation can then be employed recursively to obtain the closed formulae (43). Remarkably the computation of 2j 1 integrals in formula (36) is then reduced to computing a single integral, which may be easily done numerically.
Although each contribution to the sum (44) is just an integral of a simple function, it turns out that the sum itself is very slowly convergent for the massive free Boson. However, at least for small integer values of n it is possible to perform the sum very precisely. This is also helped by the fact that the function (37) takes particularly simple forms for n = 2, 3, 4 and 6 iF j (y, 2) sinh y = 2 2(j 1) , iF j (y, 6) sinh y = 1 6 Because of these simple, closed expressions we were able to evaluate the sum (44) Table 1: Numerical evaluation of the sum of (44) for n integer with truncation at j = 2000. The agreement with the predicted values 4 n (as given by (5)) is very good even though the sum (44) is very slowly convergent.

4
n for n integer with great precision (for the data in Table 1 the error remains below 2%). However, as discussed in [41], when n is non-integer, the integral (44) requires a non-trivial analytic continuation. In that case, additional terms need to be added to x TT which account for the residues of the poles of F j (y, n) that cross the real axis as n ! 1 + . The summand in the function F j (y, n) has kinematic poles at 2y ± (2p 1)i⇡ = (2kn + 1)i⇡ and 2y ± (2p 1)i⇡ = (2kn 1)i⇡ for k 2 Z.
This poles are due to the presence of kinematic poles of the two-particle form factor (18) at ✓ = i⇡ and ✓ = i⇡(2n 1), together with its periodicity property f (✓; n) = f ( ✓ + 2⇡in; n). This gives rise to four families of poles with corresponding residues of the function inside the sum (44) given by: Note that all these residues are zero for n integer (due to the presence of the sinh(i⇡kn) function) so that they only contribute for non-integer n. Once we have understood the pole structure of the integrand (44) we must then investigate which of these poles cross the real line in the limit n ! 1 + . This is relatively intricate as the position of each pole depends on n, k, j and p. To ease understanding let us consider a simple case as an example: n = 3 2 and j = 2 in the sum (44). We know that 4 3 2 = 0.14. If we simply evaluate (44) with as much precision as possible we obtain the value 0.0736 which strongly disagrees with the CFT formula. Moreover this disagreement cannot be entirely explained simply by the truncation of the sum (44). This disagreement is in fact due to the presence of poles of the function F 2 (y, 3/2) in (44) which cross the integration line (e.g. the real axis) as n approaches the value 3/2. If we now consider the generic poles (52) and the definition (37) we see that for j = 2 we can only have p = 1, 2. For p = 1 the four families of poles labeled by the integer k are: Note that the poles at ikn⇡ are not double, but arise as single poles of both summands in the function (37). It is clear that these poles are always above the real line (for k > 0) or below the real line (for k < 0), that is they never cross the real line, even if n is small. Similarly the poles at (kn ± 1)i⇡ remain above the real line whenever k > 0 or below the real line if k < 0 as n approaches 3 2 . Consider now the poles corresponding to p = 2. We now again have the following four families: We have already seen above that the poles y 1 and y 3 never cross the real line, so we may at most have some contributions from y 2 and y 4 . For k > 0 and n positive and large both families of poles are above the real line. However, for n = 3 2 we see that the pole (kn 2)i⇡ crosses the real line for k = 1. Similarly, for k < 0 and n positive and large all poles are in the lower half plane but the pole (kn + 2)i⇡ crosses the real line for 3 2 and k = 1. In summary, there are two poles for j = p = 2 located at ± i⇡ 2 that cross the real line as n ! 3 2 . The corresponding residue contributions are 2⇡i(R 2 (2, 2, 1, 3/2) R 4 (2, 2, 1, 3/2)) = 3i Therefore, the addition of the residua of these two poles improves the estimate of the conformal dimension from 0.0736 to the value 0.0885 (note that the formula (44) gives -4 n , hence the minus sign of (61)). Similarly, the addition of poles for higher values of j will bring this value ever closer to 4 3 2 = 0.14 as shown in Fig 2. In the general n case, in order to fully identify those poles that will cross the real line we find once more four cases: This gives the analytically continued values The shifts q 1 , q 2 take the value 1 when n[ p 1 n ] = p 1 and n[ p n ] = p, respectively and are zero otherwise. Here the symbol [.] represents the integer part. To conclude this section, we note once more that both the sequence (44) and (63) are very slowly convergent. Even after the inclusion of 2000 terms in Table 1 agreement with analytical results is not perfect. The values depicted in Fig. 2 show almost perfect agreement with the analytical result but only because we have managed to sum (44) and (63) almost exactly. We achieved this by first truncating each sum up to j = 150 and then carrying out a linear fit of the logarithm of individual terms from j = 20 to j = 150 against log j. Such fit is extremely precise and we could then use it to carry out the rest of the sum (from j = 151 to 1). This latter sum turns out to still give an important contribution to the final value (around 8%). This is rather surprising given that a previous investigation of the free Fermion, where very similar expressions emerge leads to rapidly convergent sequences and very accurate predictions, as shown in [41]. Despite this observation, the numerical results depicted in Fig. 2 provide strong evidence for (63) representing the correct analytic continuation to n non-integer. Despite the slight disagreement with the analytical formula, it is clear from Fig. 2 that (44) either under-or overstimates the value of 4 n if n is non-integer and that it has oscillations which are smoothed out by the addition of the residues associated with the poles (62) which cross the real line as n approaches 1.
As we will see, convergence issues appear to be a typical feature of the massive free Boson theory and will feature again when we compute other physical quantities. We will discuss their possible origin in sections 6, 8 and Appendix B.

The two-Point Function hT (0)T (`)i
n Once again we use the formula (113) to carry out the sum over the indices p i in (31). The result depends on the sum of all rapidity dependencies entering the two particle form factors f (✓; n) in the sums. In this case this leads to a remarkable simplification as by construction. This means that the value of the sum in (31) is given by the particular limiting case of (113), which after analytic continuation in n is given by (114). Thus we have that or, after introducing the variables x i defined earlier (33)- (34), integrating over the variable x 2j and expanding the Bessel function as in (12) we obtain The integrals above are special cases of formula (43) which allows for their direct evaluation. Note that they are entirely independent of the value of n which only enters through the function nh(j, n). It is easy to show that (this is just a special case of (43)) Substituting (67) into (66) we obtain the sum Employing the definition of h e (j, n) given in (116) we have that and similarly for n odd. We find .
Here the e and o superindices indicate analytic continuations of x T T from n even and odd, respectively. The values above are exactly those predicted by CFT as seen in the leading contributions to (6). The expected results are obtained for generic n, thus showing that the functions h e,o (j, n) indeed capture the correct analytic continuation from n integer and even or odd to n real and positive. In particular, by setting n = 1 in x e T T we recover the value 1 4 in line with CFT predictions for the logarithmic negativity [19,20,30]. It is worth emphasizing that the results (70) follow from exactly re-summing all the infinitely many terms resulting from a form factor expansion, something that can rarely be achieved for any QFT local fields.  (14) and (16). We may now evaluate K TT = 2 loghT i n by employing (16) and the results of the previous section. In particular, we will use the variables (33)- (34) in terms of which we obtain where u j (n) = and d j are the functions defined in (13) now expressed in terms of the variables x 1 , . . . , x 2j 1 as The integrals u j (n) can be computed by means of Monte Carlo integration procedures (except for n = 1 where u j (1) = 0 directly from the definition) leading for instance to the values depicted in Fig. 3 In principle we could use these fits to evaluate the sum (71) up to large values of j. However, it is clear from the fits (74) that which means that the sum (71) is divergent, even if each individual integral u j (n) takes a finite value. This is an a priori surprising result which needs to be physically understood. A discussion and interpretation of this result will be presented in section 6. We will show that despite the sum (71) being divergent, we may still extract useful information from it.

The two-Point Function hT (0)T (`)i
n Let us go back to formulae (6), (14) and (16) and let us examine the next to leading order correction to (6), that is the ratios of expectation values and three-point couplings of the twist field given in (6). According to (16) and employing once more the variables (33)-(34) we can write where and d j are the functions (73).   A crucial feature of the integrals v j is that they are n-independent. Besides the case j = 1 we have not found closed formulae for j > 1 but we have been able to compute the integrals very precisely through Monte Carlo integration procedures. Fig. 4 shows the numerically obtained values of v j for j  12. These values are very precisely fitted by the curve We may now use this fit to extrapolate to larger values of j (rather than carrying out the integrals). In this way, we will be able to perform the sum (76) up to very large values of j. For n odd, we obtain the values reported in Table 2.
hT in so that the formula (76) provides a prediction for a ratio of two universal QFT quantities. Further, because the full n-dependence is encapsulated by the function h o (j, n), we may also use the formula (76) for non-integer or even values of n. A graphical representation of the values obtained for n  7 is given in Fig. 5. As can be seen, we obtain a smooth function of n which displays linear behaviour for large n. In particular it is This result is exactly what we would expect since hT i 1 = 1 and 1 = 0 (for n = 1 we only have one replica so the twist field becomes the identity field). Also, for n odd the field T 2 = T and so C T 2 T T = C T T T . For n = 1 this is the structure constant associated with the identity field which should be also 1. We may attempt now to perform the same sum (76) employing the T T evaluated using formula (76) with h o (j, n) for various (non-integer and integer) values of n and summing as many terms as needed to ensure convergence. The blue line is the function 0.34 + 0.13 n + 0.215n. The fit is extremely good indicating that the ratio of C T 2 T T and hT i n decays exponentially for n odd. Here, as before, we have set the mass scale m = 1.
function h e (j, n) defined in (116). This should provide an analytic continuation from n even of the function Unfortunately, the sum (76) (similar to (71)) is divergent for n even. The di↵erence with respect to the n odd case is due to the asymptotic properties lim j!1 h e (j, n) = 1 n and lim It is however possible to evaluate C T 2 T T by subtracting the divergent sum (71) from (76) in such a way as to remove all dependence on the expectation values. In other words, we may compute In particular, for n = 2 we can employ the fact that u j (1) = 0 and u j (2) = 2 2(j 1) v j (this is due to the equality (46)) to find The result above follows trivially from the property h e (j, 2) = 1 2 8 j and gives a neat example of how the di↵erence of two divergent series may produce a convergent one. The sum above is identically zero (irrespective of the values of v j ), giving us the exact result lim This result is in agreement with what we expect from CFT considerations. It was first argued in [36] that the field T 2 is nothing but the identity field for n = 2 and so the result follows from the CFT normalization of correlators. For other values of n we rely on the numerical fits v fit j and u fit j (n) which are of course not exact. However, within the error of these fits we have been able to show that the sum (83) is indeed convergent. Fig. 6 shows our results for several even values of n. The solid (green) and dotted (black) lines presented in Fig. 6 are fits which provide a numerical analytic continuation from n even to n real and positive. In particular our numerical values for C T 2 T T are very well fitted by either C 1 (n) = e 1.074 1.064 n 0.274n or C 2 (n) = e 0.308 0.311n+0.456 log n and allow us to obtain the following values These results can be compared to an analytic prediction in [20] where the value of the structure constant was computed for the compactified free Boson in the double limit n ! 1 and ⌘ ! 1 where ⌘ is the compactification radius. The value predicted in [20] is where A = 1.2824... is Glaisher's constant and in [20] this number was called P 1 1 . This analytical value lies slightly above both values (86) and is closest to the value C 2 (1). This highlights the di culty of performing reliable analytic continuations based solely on a (small) number of numerically obtained values. The fact that the fit C 2 (n) seems to work best near n = 1 is natural once we notice that the analytical prediction (dashed line) also has an expansion of the form a + bn + c log n for large n (see discussion in section 6). T T for even values of n (dots). The solid and dotted lines provide two possible fits of the points obtained. The solid (green) curve is the function log(C 1 (n)) = 1.074 + 1.064 n + 0.274n and the dotted (red) curve is the function log(C 2 (n)) = 0.308 + 0.311n 0.456 log n. As we can see both are extremely good for the points we have and yet their curvatures around n = 1 are quite di↵erent. The dashed blue line represents an analytical prediction given in [20] (see section 6 for discussions on this comparison).

Interpretation of Divergent Series and log log-Corrections
In the previous sections we have shown that a form factor approach allows accurate access to leading and next to leading order short distance corrections to twist field two-point functions. Such corrections involve universal quantities which characterize both the massive theory and its conformal counterpart. They are given by expectation values hT i n and the three-point coupling C T 2 T T of twist fields, both of which are generally very hard to compute analytically, even for free theories. A feature of particular interest is that for the massive free Boson the form factor expansions of give rise to divergent sums. Our interpretation of such divergences is that they arise from the presence of (unaccounted for) logarithmic divergences of the corresponding correlators. In other words, the formulae (5) and (6) do not capture the true`-dependence of the correlators and this in turn means that the identification of (16) is not entirely justified. However, remarkably, our functions K T T and K TT still capture universal QFT information which is revealed when special divergence-canceling combinations of correlators are evaluated.
We propose that (5) for n even (90) where r 1 (n) and r 2 (n) are unknown functions and p is a constant. An analytic calculation for the massless free Boson showing the emergence of a log(log`) correction in (89) will be presented shortly in [25]. Obviously the presence of the constant p is equivalent to a redefinition of K T T and K TT and this means that there is naturally a certain ambiguity in the identification of the expectation values and three-point couplings through this approach.
From our form factor computation, the n-dependence of the functions r 1 (n) and r 2 (n) can be fixed by imposing the cancellation of divergences that we have numerically observed. There are three key observations that we may use: 1) The fact that the sum (83) is convergent implies that 2) Another numerical observation which was suggested by preliminary results of [25] is that the ratio 2 log hT i also admits a convergent form factor expansion representation even if the expansion of hT i n itself is divergent. The cancellation of divergences in (92) is equivalent to requiring r 1 (n) (n 1)r that is r 1 (n) = r(n 1) with r constant. 3) We may even determine the sign of this constant r even if form factors alone cannot fix its value. This is because the expansion (71) is not just divergent but tends to +1 (all functions involved in the sum are positive definite). This observation means that whatever the correction to the leading log`term is, its e↵ect should be to reduce its value (note that the factor K TT appears with a negative sign in the expansion (5). This means that r > 0.
The presence of logarithmic divergences in the correlators of the massive free Boson is not entirely surprising as we are dealing from the beginning with an underlying logarithmic CFT. It was shown in [43] in complete generality that an additional contribution to the Rényi entropy of the form r log(log`) will always emerge when dealing with logarithmic CFTs [42] (see eq. (13)). In this context, the coe cient r was shown to be a positive integer, which is related to the algebraic structure of the CFT.
In the specific case of the non-compactified massless free Boson, additional log log divergences of other twist field correlators at criticality are also found when studying the LN of adjacent regions in the compactified free Boson in the limit when the compactification radius ⌘ ! 1 [20] and also when studying the EE of two disconnected regions. The presence of a log(log`) correction in (89) can in fact be inferred directly from the results of [24] where the four-point function hT ( `1)T (0)T (`2)T (`2+`3)i n of the compactified massless free Boson was investigated. Combining Eq. (4) and Eq. (66) in [24] it was found that for large compactification radius ⌘ at criticality where g(`1,`2,`3) is a known ratio of lengths and x =`1`3 (`1+`2)(`2+`3) is the usual cross-ratio. It is easy to see that the leading term in the expansion of the functions above about`2 = 0 (or x = 1) is given by (95) Therefore we have that the von Neumann entropy diverges as Another correlator of interest was considered in [20]. It was shown that the LN of adjacent regions in the massless (non-compactified) free Boson behaves as where y =`1`2 1 +`2 and P 1 is the inverse three-point coupling (see [36] for a discussion of various equivalences between three-point couplings of twist fields such as the one used above). Note that once more a log log correction is present which appears with the same coe cient as in (90) when the same limit n e ! 1 + is taken. The identification (99) once more su↵ers from the ambiguity of whether or not the term 1/2 log 2 which is included in the double logarithm should be identified as part of the three-point coupling. The numerical comparison in Fig. 6 suggests that log P 1 should indeed be identified with log C T 2 T T in our setup. In fact, we can even compare our results to those given in [20] beyond the n = 1 analytic continuation. A full expression for the constant P n with n even was derived in [20] and is given by However, apart from the inherent ambiguity in the definition of the constants P n emerging from the presence of terms such as log 1 2 log` above, there is also another ambiguity emerging from the fact that the computations in [20] are done for a compactified free Boson and depend on the compactification radius through a factor ⌘ (n 1) 2 . In [20] it is argued that by taking first the limit n ! 1 and then ⌘ ! 1 results for the decompactified free Boson should be obtained, among them the constant P 1 . However, if we are to compare our numerical values in Fig. 6 to formula (100) then the presence of the factor ⌘ (n 1) 2 can play a role. It is of course rather hard to asses how this infinite constant (for n even and ⌘ ! 1) a↵ects the definition of C T 2 T T . Our benchmark has been to use the fact that lim ne!2 C T 2 T T = 1. It turns out that if we identify C T 2 T T with P 1 n as given above, then this condition is not satisfied. However, because of the intrinsic ambiguity on how ⌘ is defined we may argue that we could always scale P n by a factor of the form q (n 1) 2 where q is a constant (this would be equivalent to scaling ⌘ ! q⌘). We may then just pick q is such a way as to ensure that log P 2 = 0. From (100) we have that P 2 = p ⇡ 2 . Therefore, by choosing q = ⇡ 2 we may construct a scaled version of P n given bỹ which automatically has the desired properties It is this functionP n which is plotted in Fig. 6 (dashed blue line). The agreement with our data is reasonably good making this identification plausible. It would be nice to have an alternative analytical derivation of C T 2 T T directly for the non-compactified massless free Boson. For comparison, it is easy to carry out a large n expansion of the function logP n : the leading terms are 1.25 + 0.34n 0.5 log n. The linear term in n is in rather good agreement with the two fits presented in Fig. 6 as they both reproduce well the large n behaviour of our data. The log n term is well captured by the second fit.
In conclusion, our form factor-based numerical and analytical results lead us to conclude that the emergence of a log(log`)) correction in (90) is closely associated with a similar correction in (89). Combining our results with the expansion (96) that follows from [24] and the suggestion coming from [25] that the ratio (92) must be finite we propose that the Rényi entropy of the non-compatified massless free Boson has the behaviour Since the log(log`) correction is independent of n the same correction should also contribute to the von Neumann EE. It is worth noting that the presence of such subleading corrections was missed in the original treatment of the non-compactified massless free Boson [31]. When starting this investigation, the presence of log(log`) term in the Rényi entropy of the non-compactified massless free Boson was entirely unexpected and, as far as we know, had not been suggested by any previous studies. We have now found that a form factor computation combined with various other results strongly suggests its presence. It is worth considering whether or not such a term is universal in the same sense as the leading log`term is. In other words, is the coe cient +1 of log(log`) in (103) a universal number? Based on our understanding to date, there are strong hints that it is not, but that it may depend on the regularisation scheme used. In particular, we understand the the study [25] produces a di↵erent coe cient, both in sign and absolute value. On the other hand, unpublished numerical studies due to Andrea Coser and Cristiano de Nobili 2 employing an infinitely long harmonic chain and subsystems of sizes varying from few sites to thousands of sites, have found no evidence of such term. We do not yet understand how these various results can be reconciled but it is something we would like to investigate in the future.

Three Point Couplings and out of Equilibrium Negativity
Another way of obtaining the value lim ne!1 + C T 2 T T is to compare to other existing numerical results. In particular, in [35] the negativity of the harmonic chain, a discrete system whose continuous limit is described by a massless (non-compactified) free Boson, was numerically investigated. In particular, the LN of the harmonic chain out of equilibrium was numerically evaluated. The set up is one in which two harmonic chains are independently thermalized at temperatures 1 r and 1 l and then connected and let to evolve unitarily (e.g. quenched) at time t = 0. The time evolution of the LN is then investigated. In this context, the authors obtained very nice numerical results which were later shown to be in full agreement with predictions based on CFT [36]. Fig. 6 in [35] is of particular interest as the "staircase" pattern of the LN, as well as the height of the steps, have a CFT interpretation [36]. In this figure the LN of adjacent regions at the same temperature 1/ is presented as a function of time t for di↵erent choices of . Consider equation (90) in [36]. This equation essentially says the following: if we take the data in Fig. 6 in [35] we will see that there is an initial region around t = 0 where the negativity grows logarithmically as and c 1 is a non-universal constant. Then, the negativity reaches a plateau which is temperaturedependent and is described by formula (88) in [36]. The value of the negativity at the plateau is given by where c 2 is another non-universal constant. The key observation made in [36] is that although (104) and (105) involve non-universal constants, the di↵erence between these constants is a universal CFT quantity related to the three point coupling lim ne!1 C T 2 T T . More precisely 2 log( lim From Fig. 6 in [35] it is easy to obtain an approximate value of c 2 as it is determined by the heights of the first plateau in each curve. Unfortunately we have not had access to the raw data so we could only determine the heights approximately. Considering the four curves in the figure we find that for = 5 the plateau is located around E 2 = 0.936, for = 10 we have a plateau at or lim ne!1 + C T 2 T T ⇡ 1.38. Given the approximate values we have used and the fact that the results of [35] are numerical (not exact), this estimate is in very good agreement with (87). This agreement is also remarkable because the CFT interpretation given in [36] did not consider the possibility of log log corrections to the LN yet, like the form factor approach, it seems to still capture universal information about the CFT. We speculate that the subtraction of the constants c 1 and c 2 has a similar divergence-canceling e↵ect as described in section 6.

Conclusions and Outlook
In this paper we have studied the short-distance behaviour of the normalised two-point functions of branch point twist fields hT (0)T (`)in hT i 2 n and hT (0)T (`)in hT i 2 n in the replica massive 1+1 dimensional free Boson theory. Our work is based on the use of the form factor approach which allows us to obtain the exact matrix elements of the fields T andT up to their expectation values. For this reason, it is natural to consider the ratios above, where the dependence on the expectation values is e↵ectively canceled out. From our numerical and analytical results we conclude that In this paper we have shown that a form factor approach can provide extremely accurate predictions for the powers b(n) and d e,o (n). Whereas b(n) may be expressed as an infinite sum of simple terms (44), d e,o (n) may be computed from and exact resummation of a form factor expansion (70). In addition, we have performed the analytic continuation from n integer to n 1 and real, so that b(n) and d e,o (n) may be obtained from form factor expansions also for noninteger values of n. This provides a powerful test of our approach to the analytic continuation of correlators of twist fields, a problem which is of key importance in their applications to measures of entanglement.
Remarkably, the form factor approach also allows us to obtain infinite-sum representations for some of the ratios (109). Interestingly the sums associated to a(n) and c e (n e ) turn out to be divergent, whereas the sum representation of c o (n o ) is not only rapidly convergent but may also be analytically continued to real n 1. In this paper we argue that the divergences we have found may be explained by the presence of the log`powers in the denominators of (108). The presence of similar corrections was noted in studies of the LN of the massless free Boson [20]. For the von Neumann and Rényi entropies of one interval, they were implicit in the results of [24] and they have now been independently derived in [25]. However, they were missed in the original treatment of the non-compactified massless free Boson [31]. The presence of log`terms in (108) directly leads to the log(log`) term in (103), that is, it leads to the prediction that all Rényi entropies as well as the von Neumann entropy of the non-compactified free Boson should be corrected by a log(log`) term. This result is surprising and it is worth considering whether or not such extra terms are universal. The evidence we have so far suggest that they are not. It appears that the computations to appear in [25] also find a log(log`) term at criticality, albeit with a di↵erent coe cient. At the same time, unpublished numerical studies due to Andrea Coser and Cristiano de Nobili have found no evidence of such terms in an infinite harmonic chain for a wide range of sub-system sizes. We must therefore conclude that log(log`) corrections to the entanglement entropy are probably non-universal and we still need to understand how and why this is the case. Interestingly, the occurrence of log(log`) corrections to the Rényi and von Neumann EE of a single interval in logarithmic CFTs [42] was shown in [43] and it would be nice to understand better how the non-compactified massless free Boson falls (or not) within this class of theories.
Despite the fact that a(n) and c e (n e ) are given by divergent sums in our form factor approach, we have numerically observed that certain combinations of these sums are convergent. In particular, it is possible to obtain convergent expressions for and a(n e )c e (n e ) a(n e /2) = C T 2 T T for n even.
These numerical observations (in some cases backed by complementary results [25] from other approaches) have allowed us to actually fix the powers of log`in (108). Employing these convergent series we have obtained numerical values of log C T 2 T T for n even from n = 2 to n = 14. We have then attempted to find a good interpolation of the points obtained that would allow us to find the value of lim ne!1 + C T 2 T T . We observed that di↵erent fits can give very di↵erent predictions for this constant (perhaps not-surprisingly as we only had few points). More generally, this provides an instructive example of the di culty of performing the analytic continuation numerically. Luckily we were able to compare these fits to an analytic prediction from [20]. We have also been able to provide a further numerical estimate of the value lim ne!1 + C T 2 T T by combining numerical results for the LN out of equilibrium in a harmonic chain [35] and their CFT interpretation [36]. Remarkably, the value obtained also agrees rather well with the analytic prediction [20].
In [30] an analytic formula for hT i n was proposed and therefore it is natural to ask whether or not our results for the ratio (92) are matched by this formula. It turns out that the agreement is poor. There are now indications [25] that the formula given in [30] was not correct mainly because the presence of log`corrections to the two-point function hT (0)T (`)i n had not acknowledged in the original computation. We hope that a comparison between our results and those of [25] will be possible in the near future.
The current work has demonstrated that the massive free Boson theory allows for a form factor treatment which we can hardly hope to emulate to interacting theories. This is on account of the simplicity of twist field form factors (and the unusual fact that they are all known). For this reason this is an ideal model for which detailed three-and four-point function form factor computations may be feasible leading to new insights into the properties of the LN and the EE of disconnected regions in gapped systems. These are interesting problems which have not yet been addressed for massive models and we hope to return to them in the future. At the same time, despite it being a non-interacting theory, the massless limit of the massive free Boson is a logarithmic CFT and as a consequence, "unusual" logarithmic divergences are present in the correlators of branch point twist fields. This gives rise to very interesting structures and in particular to log(log`) corrections to the Rényi and EE of one interval, in agreement with various predictions [43,25,20,24]. It would be interesting to investigate these unusual corrections further and to establish more rigorously whether or not they are universal and/or numerically observable in some discrete realization of the non-compactified massless free Boson theory. most likely non-universal. Finally, O.A. Castro-Alvaredo is also grateful to the Isaac Newton Institute (Cambridge) for hospitality during the scientific programme Mathematical Aspects of Quantum Integrable Models in and out of Equilibrium, January 2016, where some preliminary results were presented and useful discussions took place.
In [30] several summation formulae involving two-particle form factors were obtained for generic 1+1-dimensional QFTs. Let f (✓; n) be the two-particle form factor as defined in (18), then these formulae specialize as follows to the massive free Boson theory The identity above is obtained by analytic continuation in n using a "cotangent trick" to turn the sum into a contour integral and then use the kinematic singularities of the function f (✓; n) to explicitly compute this integral. The formula (112) can be easily generalised (by induction) to include additional sums. A similar procedure was employing in [41] for the free Fermion theory. The formulae for the free Boson are almost identical, up to a few sign changes. The case of interest here corresponds to performing and odd number of sums. Employing once more the notation x j := x + 2⇡ij, we find An important property of (113) is its behaviour in the limit P 2j i=1 x i = 0. Although the presence of the sinh-function in the numerator suggests the sum should be zero, this is not the case as the sum in j develops kinematic poles in the same limit. Those contributions arise in three particular instances of the sum: First, when p = 1 second, when p = kn, and thirdly when p = kn + 1 with k 2 Z. In each case a kinematic pole is captured. How many of these poles contribute to the sum depends therefore on the relative values of p and n. In general we may write that: with h(j, n) : thus, for n > j only the first term contributes. In many computations it will be important to analytically continue h(j, n) from n even or n odd. It is natural to define and to be the analytic continuations of h(j, n) from n even and from n odd, to n real and positive, respectively.

B Numerical Precision
Many of the results of this project are, at least in part, based on the use of numerical algorithms and di↵erent types of approximations. For this reason we think it is important to give a brief discussion of how we think these procedures a↵ect the precision of our results. From the scattering theory point of view, we have dealt with the simplest theory we could possibly imagine: non-interacting particles with two-particle scattering matrix S(✓) = 1. In that respect, the number and nature of the challenges we have faced when attempting to evaluate physical quantities numerically has been somewhat surprising. There were two main challenges: 1) Physical quantities are given in terms of (infinite) slowly convergent sums: this is particularly striking for the series (44) and for its analytically continued version (63). However, it also emerges less visibly in the series representations (69). In this case, the slow convergence is less apparent because we are able to perform all sums analytically, however it is easy to see that convergence of this sequence is also slow. In all cases, it can be shown that the terms in the series involved decay roughly as 1/j s for s > 1.
2) Physical quantities are given in terms of divergent sums: in some cases, the situation is even worse than 1) and we actually have divergent representations for quantities of physical interest such as the expectation value (71) and the ratio (81). In fact, we are not aware of any other form factor calculation where such sequences emerge. As discussed in the paper, we have discovered that convergent sequences may still be constructed by combining several divergent ones in a physically meaningful way.
Interestingly, although the free Fermion theory which was studied in [41] is in many ways very similar in structure and in the sort of computations involved, it turns out that none of these issues arise. We believe that divergence and poor convergence of for factor expansions are both related to the presence of log`divergences in the correlators of the massless free Boson, a feature which is of course absent for the free Fermion (and most other theories of interest). The cancellation of divergences by combining several divergent sum is then equivalent to the statement that although some correlators diverge as r x (log`) y for some x, y > 0 the log`divergences can be canceled by computing instead ratios of several correlators. For those sequences which are convergent and which we could evaluate, we used a range a methods that allowed us to sum if not the full sequence at least a large number of terms. The only sequence which we could sum exactly was (69). A particularly neat example is the sequence (76) for n odd. For all the values of n we considered, the sum has fully converged after 200 terms (which is fast, compared for example to (44) and (63)). However, even to evaluate this sequence we have employed the fit (79), meaning that the values we summed are the result of approximating the integrals (77) by a fit (rather than evaluating (77) for all j up to 200). Instead the integrals (77) were evaluated using a Monte Carlo algorithm for j = 1, . . . , 12 (these are the dots in Fig. 4) and then a numerical fit of v j was performed based on these first 12 values. The error on the values depicted in Fig. 5 which is induced by this procedure is hard to estimate, although based on the errors derived from the Monte Carlo (which are very small) and the errors of the fitting v fit j , which are also small, we expect the results in Fig. 5 to be very accurate.
This same technique of using fits to be able to sum further terms has been employed in our evaluation of the structure constant C T 2 T T for n even in (83). Our objective was to find an appropriate fit of those values that would allow us to carry out the numerical extrapolation to n = 1. As shown in Fig. 6 we were only partly successful in this. Indeed, various sensible interpolating functions turned out to exhibit very di↵erent behaviours precisely around n = 1. This is due to a large extent to the fact that, unfortunately we only had very few values to fit. This in turn is due to the di culty of finding accurate values for the integrals u j (n) through Monte Carlo for large j. The numerical values depicted in Fig. 6 were obtained by summing (83) and employing interpolating functions for the integrals v j and u j (n). Again, it is hard to estimate the percentage of error induced by employing these fits, mainly because the functions (74) where obtained by interpolating with only few points (see Fig. 3).
There is another approximation involved in evaluating (83) and that is the fact that the divergent parts of the three sums involved do not (numerically speaking) cancel exactly. For all values of n we find that the divergent terms in the sum, which diverge as 1/j, have coe cients of the order of 10 2 (which are neglected to achieve convergence).
Let  n be the percentage error associated with the truncation and cancellation of divergences in the sum (83) for each value of n. In order to extrapolate to n = 1 we employed a fit of the form: log C T 2 T T ⌘ y fit (n) = a + bn + c n Let y(n) be the numerical values to be fitted. We define the error (n) =  n y(n) where  n is the percentage error on the value y(n) stemming from the various approximations discussed above (e.g if the error on the y(n) were of 10% then  n = 0.1). In general we would expect that error to be approximately the same for every n. However, in the case of the data in Fig. 6 we have seen that the value y(2) = 0 is actually exact (see eq. (84)). We may assume that all  n are the same for n > 2 and  2 = 0. In order to compute the values (and their relative errors) of the fitting constants a, b and c we perform a least-squared fitting [44]. In particular, the constants a, b, and c are such that minimise the quantity: 2 = X n even (y fit (n) y(n)) 2 (n) 2 , Taking derivatives with respect to the coe cients a, b, and c we can minimise 2 : 8 > > > > < > > > > : Since only the y(n) values are a↵ected by error, the associated error of the a, b, and c constants is given by: and their errors. Assuming for instance that  n = 0.1 for all n > 2 (that is a 10% error on each of the numerical values y(n)) we obtain errors: a = 0.157838... b = 0.0201026... and c = 0.240441...
And the value of C T 2 T T in n = 1 is then given by: where the error at n e = 1 is given by: C T 2 T T = e (a+b+c) (| a| + | b| + | c|) = 0.506659...
A similar study could be performed for the other fit presented in Fig. 6. For a fit of the form y fit (n) = a + bn + c log n and assuming again  n = 0.1 for n > 2 we obtain a = 0.308 ± 0.028, b = 0.312 ± 0.031, and c = 0.456 ± 0.119.
As we can see not only the second fit was in better agreement with the analytical prediction around n = 1 but it entails in general a much smaller error.