Entanglement entropy along a massless renormalisation flow: the tricritical to critical Ising crossover

We study the R\'enyi entanglement entropies along the massless renormalisation group flow that connects the tricritical and critical Ising field theories. Similarly to the massive integrable field theories, we derive a set of bootstrap equations, from which we can analytically calculate the twist field form factors in a recursive way. Additionally, we also obtain them as a non-trivial roaming limit of the sinh-Gordon theory. Then the R\'enyi entanglement entropies are obtained as expansions in terms of the form factors of these branch point twist fields. We find that the form factor expansion of the entanglement entropy along the flow organises in two different kind of terms. Those that couple particles with the same chirality, and reproduce the entropy of the infrared Ising theory, and those that couple particles with different chirality, which provide the ultraviolet contributions. The massless flow under study possesses a global $\mathbb{Z}_2$ spin-flip symmetry. We further consider the composite twist fields associated to this group, which enter in the study of the symmetry resolution of the entanglement. We derive analytical expressions for their form factors both from the bootstrap equations and from the roaming limit of the sinh-Gordon theory.


Introduction
Our understanding of many-body quantum systems both at and out equilibrium has been dramatically boosted in recent decades. An important part of this progress is due to the research carried out on entanglement in extended quantum systems. Being the fundamental feature of quantum mechanics, entanglement is responsible of a plethora of quantum phenomena, novel phases of matter, and collective effects [1][2][3][4]. Different quantities have been introduced to characterise and measure the amount of entanglement. The most prominent ones are the Rényi entanglement entropies. If we consider an extended quantum system in a pure state |Ψ⟩ and we take a spatial bipartition into subsystems A and B, then the Rényi entanglement entropies are defined as S n (ρ A ) = 1 1 − n log Tr(ρ n A ), (1) where ρ A is the reduced density matrix that describes the state of subsystem A, which can be obtained by taking the partial trace to the complementary subsystem B, ρ A = Tr B (|Ψ⟩ ⟨Ψ|). In the limit n → 1, Eq. (1) gives the von Neumann entanglement entropy S(ρ A ) = − Tr(ρ A log ρ A ).
The Rényi entanglement entropies are particularly interesting quantities to study in quantum field theories (QFT). In the ground state of two-dimensional conformal field theories (CFT), they grow logarithmically with the subsystem size, violating the area law, and they are proportional to the central charge of the theory [5,6]. In general, using the path integral approach, the ground state Rényi entanglement entropies can be cast as partition functions of the field theory on a Riemann surface. Alternatively, they can also be computed as correlation functions of branch points twist fields inserted at the end-points of subsystem A, which are spinless primaries in CFTs [6,7].
In the renormalisation group picture of QFTs as perturbed CFTs, an important result in two dimensions is the Zamolodchikov c-theorem [29,30], which describes the loss of information about the short-distance degrees of freedom along the flow. Employing the ∆-sum [31], it was found in Refs. [14,19] a function along the renormalisation flow associated to the branch point twist fields with the same qualitative behaviour as the Zamolodchikov c-function. This ∆-function monotonically decreases with the distance and it is equal to the scaling dimension of the twist fields at the IR and UV fixed points of the flow. A different c-function can also be directly constructed from the entanglement entropy [32].
In this paper, we investigate the ground state Rényi entanglement entropies in the massless QFT associated to the renormalisation group that connects the tricritical and critical Ising theories by perturbing the former with a relevant field. This theory is the simplest member of the well-known family of massless renormalisation group flows that have as UV and IR fixed points two consecutive A-series unitary conformal minimal models [33][34][35][36][37]. The form factor bootstrap program has been successfully applied in Ref. [38] to certain correlators along the tricritical-critical Ising flow. Here we extend it to the branch point twist fields and we obtain explicit expressions for the two and four-particle form factors. To this end, we follow the same strategy as in the massive case, we write the set of form factor bootstrap equations that take into account the particular exchange properties of the twist fields and we propose a general ansatz for their solution. Furthermore, we also derive the two and four-particle form factors along the massless flow from the roaming limit of the sinh-Gordon ones. By analytically continuing the scattering matrix of the sinh-Gordon model, one can find the Zamolodchikov's staircase model [39], a two-dimensional integrable scattering theory that describes a renormalisation group flow which interpolates between the successive A-series unitary conformal minimal models. It has been shown [40,41] that the form factors of different fields in the tricritical-critical Ising model flow can be obtained as roaming limits of certain form factors of the sinh-Gordon theory. We show here that a similar strategy holds for the twist field form factors. We also study the ∆-function associated to the twist fields along the flow, finding that it is monotonic and correctly reproduces their scaling dimension at the fixed points.
In the last years, one of the main research lines in the study of entanglement in extended systems has been its interplay with symmetries. If the system presents a global symmetry, the entanglement entropy can be further decomposed into the contribution of each symmetry sector [42][43][44]. The massless renormalisation flow between the tricritical and the critical Ising theories has a global Z 2 symmetry. Let us denote by Q the charge operator that generates the symmetry and assume that it is the sum of the charges in subsystems A and B, Q = Q A + Q B . In that case, the reduced density matrix admits the following decomposition in symmetry sectors where ρ A,q = Π q ρ A Π q /p(q), Π q is the projector onto the eigenspace of Q A with eigenvalue q, and p(q) = Tr(ρ A Π q ) guarantees the correct normalisation of ρ A,q , i.e., Tr(ρ A,q ) = 1.
The symmetry-resolved entanglement entropies S n (ρ A,q ) quantify the amount of entanglement in each charge sector. One usual way of calculating them is through the charged moments of ρ A , Z n (α) = Tr(ρ n A e iαQ A ), by applying the Fourier representation of the projector Π q . For a Z N symmetry group, q can only take values q = 0, . . . , N − 1 and Therefore, one can write The charged moments (3) were initially introduced in an independent way in Refs. [45,46] in the context of holography. After Ref. [43] pointed out its connection with the symmetry-resolved entropies (5), both quantities have been intensively studied in two-dimensional CFTs [43,44,[47][48][49][50][51][52][53][54][55][56][57][58] as well as in free and integrable QFTs [59][60][61][62][63][64][65][66]. Symmetry resolved entanglement has also been investigated in other systems such as lattice and spin models [67][68][69][70][71][72][73][74][75], as well as in ion trap and cold atom experiments [76][77][78][79]. Similarly to the neutral moments Z n (0) of ρ A , the charged ones can be expressed as partition functions of the field theory on a Riemann surface, but now in presence of an external magnetic flux. The charged moments Z n (α) can also be obtained from a correlator of composite branch point twist fields, which also take into account the non-trivial monodromy between the sheets of the Riemann surface due to the magnetic flux. The form factor bootstrap techniques developed for the standard branch point twist fields in integrable QFTs have been extended to the composite ones to study U (1) symmetries in free [60] and interacting QFTs such as sine-Gordon model [61], the Z 2 symmetry of massive Ising and sinh-Gordon models [62,63], and the Z 3 symmetry in the 3-state Potts model [64]. In this work, we consider the composite branch point twist fields associated to the Z 2 spin flip symmetry of the tricritical-critical massless flow, and we apply the bootstrap approach to obtain analytic expression for their two and four-particle form factors. Using them, we also investigate the corresponding ∆-function.
The paper is organized as follows. In Sec. 2 we introduce the system we are interested in, the massless renormalisation flow between tricritical and critical Ising theories. In Sec. 3, we review how entanglement entropies can be calculated in QFT as correlation functions of branch point twist fields and we discuss its spectral expansion in terms of form factors. In Sec. 4, we obtain the bootstrap equations for the standard branch point twist field form factors, we Table 1: Kac tables of the tricritical (Table 1a) and critical (Table 1b) Ising CFTs [80,81].
In each case, we report the conformal dimension of the primary fields ϕ r,s of the theory. The vertical and horizontal axes correspond to the s and r indices respectively.
propose a general ansatz for their solution, and we derive the explicit expressions for the two and four particle form factors. In Sec. 5, we repeat the same reasoning for the form factors of the composite twist fields associated to the Z 2 spin flip symmetry of the massless flow. In Sec. 6, we rederive the twist field form factors by taking the roaming limit of those of the sinh-Gordon model. In Sec. 7, we first study the ∆-function along the flow of the (composite) twist fields and we derive a cumulant expansion for the entanglement entropies. We end up in Sec. 8 with some conclusions and future prospects. We also include two appendices where we discuss in detail the derivation of some of the results presented in the main text.
2 The massless RG flow from the tricritical to the critical Ising theory In this paper, we investigate the ground state entanglement entropy along the renormalisation group flow that connects the tricritical and critical Ising CFTs, which are respectively the unitary minimal models M 4 and M 3 [35][36][37] with central charges [80,81] c UV = 7 10 , and c IR = 1 2 , and with Kac tables reported in Table 1. This integrable QFT, usually denoted as A 2 , is the simplest member of the infinite family of massless theories A p that interpolate between two consecutive A-series diagonal conformal minimal models M p+2 → M p+1 with central charges c UV = 1 − 6 (p + 2) (p + 3) , and c IR = 1 − 6 (p + 1) (p + 2) .
This family of integrable RG trajectories is obtained by deforming the UV CFT M p+2 with its relevant field ϕ 1,3 [33][34][35][36][37]. In particular, in the tricritical Ising CFT, ϕ 1,3 corresponds to the vacancy density field with conformal dimension h 1,3 = 3 5 (see Table 1a) [35][36][37]. In the Euclidean formalism, the action A flow of this flow takes the form where λ is a dimensionful coupling and, importantly, λ is positive, since for negative coupling a different massive integrable theory is obtained. Several other families of massless integrable flows have been identified as well [82][83][84][85][86][87][88][89]. The masslessness of the flow described by Eq. (8) can be understood by recalling that the tricritical Ising CFT M 4 is one of the simplest examples of superconformal theory [90][91][92]. The deformation with the vacancy field ϕ 1,3 leads to a spontaneous supersymmetry breaking [93] which gives rise to right-and left-moving massless Goldstone fermions ψ,ψ and ensures that the theory has vanishing mass gap.
In Ref. [93], it was shown that the low energy behaviour of the massless flow is described by the effective Lagrangian that is, the TT deformation of the critical Ising model. Notice that the Majorana fermions ψ,ψ of the Ising model are now identified with the Goldstone fermions of the massless flow A 2 , which are the only stable particles in this theory [94]. It is worth stressing that the massless flow at low-energies is described by a TT -deformed CFT. Such theories have been studied in great detail [95][96][97][98][99][100][101][102][103][104][105][106][107][108] and hence they provide non-trivial benchmarking for some of our results. The massless flow (8) as well as the effective Lagrangian (9) possess a mass scale M , which plays the role of the momentum scale at which non-trivial scattering happens between the fermions. We can parameterise the energy and momenta of the right-and left-moving Goldstone fermions in terms of a rapidity variable θ and of this mass scale M as [94] Since the massless fermions are the only stable particles, they form a complete basis of asymptotic states, which in the rapidity parameterisation read as which contains r right-moving and l left-moving fermions. If the rapidities are ordered as θ 1 > θ 2 > . . . > θ r and θ ′ l > . . . > θ ′ 2 > θ ′ 1 , then the set of states (11) corresponds to in-states, whereas the opposite ordering results in out-states. Different orderings are linked by scattering processes between the particles. Since the theory is integrable, the scattering of particles is completely elastic, preserves particle number and rapidity, and is fully characterised by the two-body S-matrices. Since the scattering of the particles is diagonal, the S-matrices are scalars and functions of the rapidity difference of the particles. In particular [94] S RR = S LL = −1 , that is, only the scattering between left-and right-movers is non-trivial. The massless flow (8) has been the subject of numerous studies. These involve its description in terms of the thermodynamic Bethe ansatz [35][36][37], or the determination of form factors, i.e., the matrix elements of the off-critical versions of the UV scaling fields, as well as certain correlation functions [38]. At the level of the free energy and form factors [40,41], the massless flow can also be recovered from the staircase model [39,82], which we shall introduce in Section 6. The model also shows interesting properties in inhomogeneous out-of-equilibrium situations as studied in [109] via generalised hydrodynamics.
To complete the brief review of this massless flow, we discuss its symmetry properties. Both the UV and IR limiting CFTs enjoy a spin-flip Z 2 symmetry under which the perturbing field also transforms trivally. Consequently, the massless flow inherits this symmetry as well. This fact can be made more transparent using the Landau-Ginzburg formalism, which allows the identification of the multicritical Ising CFTs with a Lagrangian [80,81]. In particular, the tricritical and critical Ising models can be described by the following actions in terms of the bosonic field φ where :: denotes normal ordering of the fields. The perturbing field of the UV theory ϕ 1,3 corresponds to : φ 4 : [80,81], which means that the action of the massless flow can be equivalently written as in which the invariance under the spin-flip Z 2 symmetry, which maps φ → −φ, is explicit.
Given the presence of a global Z 2 symmetry, a relevant question is whether the ground state entanglement entropy along the massless flow (8) can be resolved with respect to it. It is not immediately obvious if a reduced density matrix of the ground state of the theory commutes with the charge operator associated with the Z 2 symmetry. While for a continuous symmetry this is ensured by Noether theorem, in the case of discrete symmetries, the existence of a local charge density is not guaranteed. In order to justify the existence of such a local Z 2 charge, we can appeal to the defect line formalism. As understood in recent years, global symmetries in QFT are implemented by topological defects [110,111], which, in the case of CFT minimal models, correspond to the Verlinde lines operators. In particular, the spin-flip Z 2 symmetry is implemented by the Verlinde line associated with the primary operator ε. Such a defect line can be restricted to the subsystem A, with two disorder operators µ inserted at the end-points [110,111] (which we discuss in more detail in the next section). This formalism has been very recently used to study the symmetry resolution of entanglement in CFTs with respect to both continuous and discrete finite groups in Ref. [57]. Since the operators ε and µ exists along the entire massless flow, the previous construction may be extended outside the fixed points.

Entanglement entropy and Branch Point Twist Fields in QFT
In this section, we review the computation of the entanglement entropies in QFT as correlators of branch point twist fields. In QFT, the non-trivial task of computing entanglement entropies can be naturally formulated via the path integral approach. The main idea is that the moments of the reduced density matrix Tr(ρ n A ) and the charged moments Tr(ρ n A e iαQ A ) can be regarded as partition functions of the QFT on a Riemann surface consisting of n replicas of the space-time that are sewed along the subsystem A in a cyclical way [5,6].
Alternatively, one can take n copies of the QFT under analysis and quotient them by the Z n symmetry associated to the cyclic exchange of the copies. In (1 + 1)-dimensional relativistic QFTs, there exist local fields in the n-replica theory, called branch points twist fields (BPTF), that implement the boundary conditions imposed on the fields in the path integral on the n-sheeted Riemann surface. These twist fields can be generalised to cases in which the boundary conditions also involve additional phases, such as in the calculation of the charged moments Tr(ρ n A e iαQ A ), in which an Aharonov-Bohm flux is introduced between the sheets of the Riemann surface. In this setup, the corresponding twist fields are called composite BPTFs and they were originally introduced in other context [14,15]. Both types of fields are associated with particular symmetries of the replicated theory, which allows us to discuss them on the same ground. Therefore, for our purposes it is useful to distinguish the following twist fields: 1) the disorder field µ associated with the Z 2 spin-flip symmetry of the massless flow; 2) the standard BPTFs, T n and its conjugate T n , which are associated with the cyclic and the inverse cyclic permutation symmetry Z n among the copies in the n-replica massless flow. These fields play a central role in computation of the entanglement entropy; 3) the Z 2 -composite BPTFs, denoted as T µ n and T µ n , which are the result of fusing the former fields T µ n (x) = : T n µ : (x) , T µ n (x) = : T n µ : (x).
Therefore, they are associated both with the Z n symmetry under the cyclic permutation of the replicas and with the global Z 2 spin-flip symmetry present in the massless flow. These composite fields play the analogous role of the BPTF in the computation of the symmetry resolved entanglement entropies [43].
These twist fields are typically non-local or semi-local with respect to other quantum fields of the theory, in particular with respect to the fundamental field or to the interpolating field, which is associated with particle creation/annihilation. Non-locality can be formulated by non-trivial equal-time exchange relations between the two fields. Let us first consider the disorder operator µ and an operator O i living in the copy i of the replicated theory. Since the disorder field introduces an Aharonov-Bohm flux in the region y 1 > x 1 , the exchange relations of these two operators can be written as We refer to κ O as the charge of the operator O with respect to the Z 2 spin-flip symmetry. In particular, the Goldstone fermions ψ,ψ which generate the asymptotic states (11) have charge κ ψ = 1, i.e., they are odd under the spin-flip transformation. The action of the standard BPTFs when winding around a field is to cyclically map it from one replica to the next, as encoded in the equal time exchange relation In the case of the composite BPTFs, the winding around them further adds a phase e iκ O π . When considering discrete groups, as the Z 2 spin-flip symmetry of the tricritical-critical massless flow, we must be careful on how we include this phase. Unlike the continuous U (1) symmetry, discussed in Refs. [60,61], here we cannot distribute the flux uniformly in all the copies by inserting a phase e iκ O π/n when moving between replicas since this operation in not compatible with the properties of the Z 2 field µ. This issue can be addressed in two different ways. The first possibility is to insert a phase e iκ O π between all the copies; this corresponds to consider the exchange relation This approach was applied in Ref. [63], but it is only legitimate when we take an odd number of replicas n = 1, 3, 5, 7, . . . in which the identity e iπn = −1 clearly holds. The other approach consists of introducing the flux only between the last and the first replicas, in such a way that the phase e iπκ O only appears when a particles moves from the n-th copy to the 1-st one, that is for y 1 > x 1 and i ̸ = n , e iκ O π T µ n (x) O i+1 (y) , for y 1 > x 1 and i = n , This choice introduces a slight asymmetry between the replicas, but it is applicable to any number of replicas n. In Sec. 5, we discuss in more detail the effect of the two conventions (18) and (19), showing that they provide the same results for correlation functions under analysis. Analogous exchange relations can be formulated for the Hermitian conjugate fields T and T µ , with the difference that they move the field from the replica i to i − 1. In the following discussion, whenever we wish to treat both the standard and the composite twist fields at the same time, we use the notation T τ n , T τ n , where τ refers either to 'µ' for the composite or to the identity for the standard BPTF.
Using the (composite) BPTFs, one can switch from a path-integral to an operator formulation of both the neutral and charged moments of ρ A , which can be defined in terms of multi-point functions of the standard or the composite BPTFs in the replicated QFT inserted at the end-points of subsystem A. In particular, when A consists of a single interval, A = [0, ℓ], we have and The twist field formalism is especially useful at criticality, where conformal invariance fixes the properties of both the standard T n and the composite branch point twist field T µ n . In particular, in a unitary CFT with central charge c, the standard twist fields T n , T n are known to be primary operators with conformal dimension [6,7] We remind the reader that, in the A-diagonal unitary minimal models M p , the central charge is given by Eq. (7). In order to identify the conformal dimension of the composite twist fields T µ n , T µ n , one can use the fact that they are the fusion of T n , T n with the disorder field µ as shown in Eq. (15). In the tricritical and critical Ising models, the field µ is the Kramers-Wannier dual of the spin field σ = ϕ 2,2 and has the same conformal dimension reported in the Kac table in Table 1 [80,81] where we denote with UV the tricritical and with IR the critical Ising models respectively. Knowing the dimension of the disorder field, the one of the composite twist fields T µ n , T µ n is obtained as [43] In particular, for the tricritical and critical Ising models, Eq. (24) gives respectively The use of CFT techniques has provided exact results for Tr(ρ n A ) in many different situations [6,7]. On the other hand, away from criticality, the exact determination of the correlation functions of Eqs. (20) and (21) is known to be an extremely difficult task, except in the case of free theories [17,70]. In integrable QFTs, however, the form factor bootstrap approach provides a powerful tool to systematically investigate and construct (truncated) multi-point functions via form factors, namely matrix elements of generic local operators between the vacuum and the multi-particle states [12,13]. Although, in principle, all these matrix elements can be analytically computed, their resummation is an unsolved problem. Nevertheless, the multi-point correlation functions at large distances are generically dominated by the first few (lower-particle) form factors. For this reason this technique applies efficiently to the infrared properties of these theories as was first shown in [8] in the case of BPTFs and entanglement. As we shall see in this paper, the above considerations do not hold in massless theories, that is, when the IR limit of the QFT is described by a CFT as well. However, we show that it is possible to identify a subset of terms in the form factor expansion whose resummation reproduces the IR CFT results, while the remaining contributions yield non-trivial predictions for the behaviour of the entropies along the flow.

Form factors and spectral representations of BPTF correlation functions
From the knowledge of the exchange relations (17) satisfied by the BPTFs, one can formulate bootstrap equations for their FF in integrable QFTs [8][9][10], generalising the standard form factor program for local fields [12,13], which for the tricritical-critical Ising flow (8) has been investigated in Ref. [38].
Let us consider the two-point correlation function of the (composite) BPTFs in the ground state of the theory and insert the set of asymptotic states (11), which form a complete basis, where τ = 0, µ corresponds to the standard or the Z 2 -composite BPTF respectively. In the multi-particle states of the n-replica theory, the subindex γ i = R, L specifies if the particle with rapidity θ i is a right-(R) or left-mover (L). Moreover, each particle is labelled by an extra index ν i which indicates the copy where the particle lives; therefore, it takes values from 1 to n and it is identified up to ν i ∼ ν i + n.
In the n-replica theory, the S-matrix connects non-trivially only particles living on the same replica, while particles in different copies do not interact and no scattering events occur between them. In light of this, the S-matrix of the replicated model takes the form [8] where S γ i ,γ j (θ) is the S matrix of the original theory, which for the massless flow (8) is reported in Eq. (12).
Since the vacuum of the theory is invariant under space and time translations, we can rewrite the spectral expansion in Eq. (26) as where E i and p i are the single particle energies and momenta reported in Eq. (10). The elementary FFs of a generic (semi-)local operator O(x, t) are their matrix elements between the vacuum and the asymptotic multi-particle states (27), i.e.
Therefore, as clear from Eq (29), the spectral sum representation of the twist field correlation function can be rewritten in terms of FFs in the following way where we switched to the Euclidean formalism for simplicity and ℓ denotes the Euclidean distance. We can see from the above formula that the computation of the twist field correlation functions can be naturally formulated by means of FFs via the insertion of a complete set of asymptotic multi-particle states. Crucially, the form factors in integrable QFTs can often be determined exactly, giving access to the corresponding correlation functions. In the following, we review some basic properties of the twist field FFs and present the bootstrap equations from which their analytic expressions can be obtained.

Form factors of the branch point twist field in the massless flow
Given the exchange properties of the standard BPTFs (16), it is possible to write down the bootstrap equations for the form factors (30) associated with these fields in integrable QFTs. Relying on earlier works [8][9][10], we can immediately specify these equations for our massless theory (8). If we denote the FFs of T n by F T |ν γ (θ; n), then their bootstrap equations can be written as [8][9][10] − i Res − i Res where we introducedγ andγ i denotes the anti-particle of γ i (which coincides with the particle in the theory under consideration). Here θ and γ, ν are shorthands for (θ 1 , θ 2 , . . . , θ k ) and (γ 1 , γ 2 , . . . , γ k ), (ν 1 , ν 2 , . . . , ν k ) respectively, with γ = R, L andR = R,L = L. In the argument of the S-matrices, θ ij = θ i − θ i .
In the massless flow (8), two particles of any type cannot form a bound state. It is also easy to see that the one-particle FFs of BPTF are vanishing. The reason for this is that these fields are neutral w.r.t. Z 2 charge. This implies that only FFs with an even number of R and an even number of L particles are non-vanishing and, consequently, the odd-particle FFs are zero.
Moreover, relativistic invariance imposes that where Σ is the Lorentz spin, which is Σ = 0 for the twist field. Another important property of form factors which will be useful in our analysis is the cluster property, studied in detail in Ref. [31] and recognised in different models, see e.g. [19,[112][113][114][115]. In the limit in which the difference between the particle rapidities diverges, the form factors factorise in the product of form factors with a lower number of particles. In our model, the clusterisation of the different particle species can be phrased as where θ and ν stand for the rapidities and replica indices of the 'R' particles, and θ ′ and ν ′ for the 'L' particles. The cluster property for particles of the same species is instead written as with an analogous expression for the clustering of 'L' particles.
Let us now use the previous axioms to construct a set of solutions of the bootstrap equations (32)-(35) for the BPTF form factors. To fix the ideas, we first place every particle on the first replica ν i = 1. A convenient ansatz for the form factors is [38] where we have r right-moving and l left-moving particles and we have defined x i = e θ i /n , y i = e −θ ′ i /n and ω = e iπ/n . Notice that we simplified our notation by omitting the reference to the replica indices. In the ansatz (40), Q T r,l are polynomials of their variables and f RR = f LL and f RL are the minimal form factors. In Eq. (40), the kinematical singularity of the FF (see Eq. (34)) comes entirely from the denominators and therefore the cyclic permutation and the exchange axioms, Eqs. (32) and (33), are automatically satisfied requiring the following identities for the minimal form factors: By prescribing that the minimal form factor f RR has no poles and has the mildest asymptotic behaviour, we end up with the unique solution and f RR = f LL , which is identical to the minimal form factor of the massive Ising theory [8].
For f RL , the defining equations are whose solution can be explicitly given based on the knowledge of the Fourier representation of the non-trivial S-matrix S RL in Eq. (12). In particular, we can write the solution as using an integral representation, or, alternatively, in terms of a mixed product integral representation which is more convenient for numerical evaluation. Notice that, for n = 1, the minimal form factor (44) reduces to the known result for a single replica obtained in Ref. [38]. Moreover, In order to fix this normalisation along the massless flow, we compute the value N −1 n of f RL in the limit θ → ∞, which reads where we have defined the sequence which is equal to Catalan's constant G for n = 1, recovering the normalisation for f RL found in the non-replicated theory in Ref. [38]. With the choicẽ we fix all the constants in the form factors for the massless flow. An important property that the f RL minimal form factor satisfies is which shall be very useful in the rest of the section. The ansatz (40), with the definitions for the minimal FFs f RR (42) and f RL (45), satisfies all the axioms for the BPTF FFs. The eventual determination of F T R,L (θ, θ ′ ; n) can be done recursively. In fact, by applying the residue axiom in Eq. (34) to the ansatz (40), one can derive recursive equations for the unknown Q T r,l (x, y; n) functions that relate Q T r+2,l (x, y; n) or Q T r,l+2 (x, y; n) to Q T r,l (x, y; n), that is, to Q T r,l functions with fewer particles. In the next subsections and in App. A.1, we explicitly demonstrate how the determination of higher-particle FFs is carried out by solving the recursive equations for the polynomials Q T r,l .

Two-particle form factors and form factors with only one species
Since the Lorentz spin of the BPTFs is zero, their two-particle FFs only depend on one rapidity variable (37), that is, the rapidity difference θ 1 − θ 2 . Recall that, because of the spin-flip symmetry, we can only have 'RR' and 'LL' form factors, which means that these quantities coincide with those of the massive Ising QFT (c.f. Eqs. (40) and (42)) up to the vacuum expectation value ⟨T n ⟩. These quantities, nevertheless, can also be easily obtained from the bootstrap equations (32), (33). For the two-particle form factors, they imply that In this case, the kinematic residue equation (34), connects the two-particle FFs and the vacuum expectation value of the twist field. We can therefore write If this formula is recast in the form of the ansatz (40), then we have the equivalent expression in which we identify where σ j is the fully symmetric polynomial of degree j in the variables x 1 and x 2 . Since in the formula above both particles live in the first replica, we slightly changed the notation, namely we denote the form factor corresponding to two right-moving particles living in the first replica F T |11 RR as F T 2,0 . In the following, we shall use this convention whenever all the particles are on the first replica. The 'LL' form factor can be obtained by replacing x 1 and x 2 by y 1 and y 2 in Eq. (54).
From F T |11 γγ (θ; n), we can obtain the form factors F T |jk γγ (θ; n) corresponding to particles in different replicas from [8] The form factorsF of the antitwist field T n can be simply obtained from those of T n through the relation As we already said, the only non-vanishing FFs with higher-particle number are those containing an even number of 'R' and 'L' particles. It is easy to see that, in the particular case of form factors only containing an even number of particles of the same type, that is, the 'RR...RR' and 'LL...LL' form factors, they exactly coincide with the standard BPTF FFs of the massive Ising theory up to the vacuum expectation value ⟨T n ⟩, similarly to the two-particle case discussed above. These form factors can be easily obtained from the two-particles ones. In particular, the form factor with 2k particles of the same type is given by for If the ordering of the indices ν i is not the canonical one, using the exchange axiom (32) one can reshuffle the particles and their rapidities to satisfy ν 1 ≥ ν 2 ≥ . . . ≥ ν 2k and apply (59). In particular, for the 'RRRR' or 'LLLL' FFs with all the particles in the same replica, we have the simple formula

Solution for the four particle 'RRLL' form factor
The first non-vanishing form factors that contain both 'R' and 'L' particles appear at the four-particle level: with any permutation of 'R' and 'L'. Similarly to the other FFs previously discussed, it is sufficient to determine only the 'RRLL' form factor with all the particles on the first replica. Using then the exchange relation (32) we can readily obtain any other sequence of the particle species, and, applying the cyclic permutation axiom (33), we can obtain FFs for particles living on different replicas. Following the notation introduced for the form factors with all the particles on the first replica, we will denote F T |1111 RRLL as F T 2,2 . In this case, the ansatz (40) takes the form Applying now the residue axiom (34) to Eq. (62), we can derive recursive equations for the H T 2,2 normalisation factor and the Q T 2,2 function. The detailed solution of this equation for the case of four-particles (RRLL) is presented in App. A.1 and here we report the results of the calculations.
The normalisation factor reads where N n is given by Eq. (47), while for the polynomial Q T 2,2 we obtain where σ i , i = 1, 2 denotes the completely symmetrical polynomial of degree i in two variables. Using these results, the final solution for the full FF is which we can also rewrite as We remark that the form factor in Eq. (65) is one of the main results of this paper. As we will show in Sec. 7, it will provide the leading correction to the IR expressions for the entanglement entropy.

Form factors of the Z 2 -composite branch point twist field in the massless flow
In this section, we derive the bootstrap equations for the form factors of the Z 2 -composite BPTFs associated to the disorder field µ along the massless flow (8) and we obtain their explicit solution for the two and four-particle cases. Similarly to the standard BPTFs discussed in Sec. 4, from the exchange properties of the Z 2 -composite twist fields (16), we can easily write down their form factor bootstrap equations. Importantly, these equations include the non-trivial phase e iπκ O in the monodromy properties due to the insertion of the disorder field µ. The asymptotic states (27) that enter in the definition of the twist field FFs are constructed from the fields ψ,ψ, which are odd under the Z 2 transformation, i.e. κ ψ = 1, and therefore we must take into account a phase e iπ when moving between replicas. However, as we discussed around Eqs. (18) and (19), we have two different ways to introduce it, either as a whole phase e iπ in each replica, which is valid only for odd n, or inserting it only in the last one. These two approaches lead to slightly different form factor bootstrap equations. In this section, we comment both choices. In particular, we will show that the two conventions give the same result for the form factors up to some (−1) factors which do not influence the final physical result.
Let us denote as F T µ |ν γ (θ, n) the form factors of the composite twist fields T µ n . If we introduce the phase e iπ on the last replica only, that is taking the exchange relations (19), the bootstrap equations take the form − i Res − i Res On the other hand, if we introduce the same flux between all the copies, we have − i Res − i Res where notations are the same as for the standard BPTF discussed in Sec. 4; in particular, we recall that γ i = R, L. Both the Lorentz spin and the Z 2 charge of the composite BPTFs are zero. Observe that the phase (−1) in Eqs. (72) and (74) as well as in Eqs. (68) and (70) is due to the non-trivial monodromy of the fields ψ,ψ with T µ n (compare with the analogous axioms for the standard BPTF in Eqs. (33) and (35)).
Similarly to the standard BPTF, only FFs with an even number of 'R' and 'L' particles are non-vanishing and, consequently, the odd-particle FFs are zero. Additionally, the FFs of the composite BPTF satisfy the momentum space clustering property in the same form as the FFs of the standard BPTF in Eqs. (38) and (39).
Analogously to what we have done in Sec. 4, let us assume the following ansatz for the composite twist field FFs in which, for simplicity, we place every particle in the first replica where we have r right-mover and l left-mover particles, and x i = e θ i /n , y i = e −θ ′ i /n and ω = e iπ/n as previously. The cyclic permutation and the exchange axioms can automatically be satisfied if the equalities are imposed, that is, the minimal form factors satisfy the non-trivial monodromy due to the insertion of the external flux. The solution of Eq. (76) can be easily obtained from the standard minimal form factor in Eq. (42) by simply introducing a factor 2 cosh(θ/2n) which changes the monodromy properties [63] f µ γγ (θ; n) = 2 cosh For f µ RL (θ; n) instead we have two possible choices. We might either choose the unaltered equation without the '−1' monodromy with the solution f µ RL = f RL , or we can also introduce the monodromy such that the solution becomes As we will later see in Sec. 6, the exponential factor in Eq. (80) also appears in the roaming limit approach. Importantly, the two choices for the minimal form factor f µ RL in Eqs. (78) and (80) are completely equivalent because for the composite BPTFs the number of 'R' and 'L' particles is always even. This implies that, in a FF, we always have the product of an even number f µ RL terms, which implies that the (−1) phases always mutually cancel. In order to connect in a clearer way with the roaming limit that we later discuss in Sec. 6, we choose Eq. (80) as the minimal form factors in the ansatz (75) for the composite BPTF. If we had taken (78), we would have got different expressions for the functions Q T µ r,l , which would differ only by products of x i and y j with the same integer powers.
We remark that, in contrast to what happened in Sec. 4, the ansatz (75) does not guarantee that Q T µ r,l is actually a polynomial. In fact, as we will explicitly show, this function is in general a rational function. The reason for this is the monodromy changing factor introduced in the minimal form factors in Eqs. (77) and (80). These terms possess additional zeros that cancel out with the denominator of the function Q T µ r,l , guaranteeing that the pole structure remains compatible with the bootstrap axioms.

Two-particle form factors and form factors with only one species
Similarly to the standard BPTFs, for the composite BPTFs the only non vanishing form factors at the two-particle level are those containing a pair of 'R' or 'L' particles, which coincide with the analogous expressions of the massive Ising QFT [8]. Alternatively, they can easily be obtained from the bootstrap equations, either from Eqs. (67), (68) or from Eqs. (71), (72). For the two-particle form factors, the bootstrap equations imply that The kinematic residue equations (69) or (73) relate the FFs to the vacuum expectation value of T µ n as −i Res The solution for the equations above can be immediately written by plugging in the two-particle FF of the standard twist field (53) the minimal form factor of Eq. (77) that takes into account the non-trivial monodromy of T µ n , obtaining where, for simplicity, we have placed every particle on the first replica. Notice that Eq. (83) is not in the form of our ansatz (75), but it can be recast accordingly as where an analogous expression for the 'LL' form factor holds upon replacing x i with y i . Since in the above formula each particle lives on the first replica, we again used the simplified notation to denote the FF, namely we write F T µ 2,0 which indicates that we have two right-moving particles on the first replica. In the following, we shall use this convention whenever all the particles are on the first replica.
The two-particle FFs with arbitrary replica indices can be straightforwardly obtained from the result (84) with all the particles on the first replica. Importantly, the different flux convention in Eqs. (67)- (70) or in Eqs. (71)- (74) only differ in some (−1) factors. In particular, if the flux is only inserted on one replica, we have The FFs of the anti-twist field T µ n denoted by F T µ a (θ, n) can be simply written as [60] If the flux is instead introduced on each replica, we have while the FFs of the anti-twist field T µ n satisfy Eq. (86). In the computation of the symmetry resolved entropy, the additional factor (−1) k−j always cancels out, leading to the same value for both choices.
Similarly to the treatment of the standard BPTFs, it is easy to see that, in the particular case of form factors that only contain an even number of particles of the same type -that is the 'RR...RR' and 'LL...LL' form factors-, they exactly coincide with the Z 2 -composite BPTF FFs of the massive Ising theory [8], as occurs in the two-particle case discussed above. Assuming that ν 1 ≥ ν 2 ≥ . . . ≥ ν 2k , they can be written in terms of a Pfaffian involving the two-particle FFs as where W µ is an anti-symmetric matrix with entries For a different ordering of the replica indices ν i , we can apply the exchange axiom (67) to reorder them in the form ν 1 ≥ ν 2 ≥ . . . ≥ ν 2k and then use Eq. (88). In particular, for the four-particle 'RRRR' or 'LLLL' FF with all particles on the same replica, Eq. (88) takes the form

Solution for the four particle 'RRLL' form factor
To obtain the first non-zero form factors that couple right-and left-moving particles, we have to move to the four-particle level, in which we find F T µ |ν 1 ν 2 ν 3 ν 4 RRLL and all the possible permutations of 'R' and 'L'. As for the standard BPTFs, it is sufficient to determine only the 'RRLL' form factor with all the particles on the first replica. In fact, using the exchange relation (67) we can directly get any other sequence of the particle species and, applying the cyclic permutation axiom (68), we can find the FFs for particles living on different replicas. If we denote the form factor F T µ |1111 RRLL as F T µ 2,2 , then it reads according to the ansatz (75).
Applying now the residue axiom to Eq. (91) we can derive recursive equations for the normalisation factors H T µ 2,2 and the Q T µ 2,2 functions in a similar way as we did for the standard BPTFs in Sec. 4.2. In App. A.2, we find the solution for the functions H T µ 2,2 , and Q T µ 2,2 , As we explain in Appendix A.2, the set of equations that allows to obtain Q T µ 2,2 recursively is under-determined. This ambiguity in the solution can be fixed by requiring that the form factor F T µ 2,2 reduces to the one of the disorder field µ in the single replica limit n → 1. One can further check that the normalisation term H T µ 2,2 also matches the one of µ in that limit. In Sec. 7, we use the ∆-sum rule to provide an additional test of the validity of our solution.

Roaming limit of twist field form factors
In the previous sections, we computed the form factors of the twist fields along the tricritical Ising massless flow directly from the solution of their bootstrap equations. In this section, in order to provide a non-trivial check of our expressions, we present an alternative derivation based on the roaming limit of the sinh-Gordon model. After reviewing the general notions of this approach, we will then use it to recover the form factors in the massless flow as the limit of those in the sinh-Gordon theory.
Let us first briefly introduce the sinh-Gordon (ShG) model. This theory is defined via the Euclidean action This is the simplest interacting integrable relativistic QFT and has been the subject of an intense research activity since many decades, see, e.g., [81,[116][117][118][119][120][121][122][123]. The spectrum of the model consists of multi-particle states of a massive bosonic particle with the dispersion relation E = m cosh θ, p = m sinh θ, where m is the particle mass. The two-particle S-matrix is given by [116] S ShG (θ) = tanh 1 2 where B is defined in terms of the coupling g appearing in the action in Eq. (94) as For the ShG model, the form factors of various operators are known [112,117,118], including the standard and the Z 2 -composite BPTFs in the n-replica theory [8,19,63]. It was observed in Ref. [39] that the S-matrix of the sinh-Gordon model can be analytically continued from the self-dual point B = 1 to complex values and resulting S-matrix is a new perfectly valid scattering theory, which has been called the staircase or roaming trajectories model. Using Bethe ansatz, it was found that, as the real parameter θ 0 increases, the c-function shows a 'staircase' of defined plateaux with values equal to the central charges of the M p unitary diagonal minimal models and, in the intervals between the plateaux, the flow was found to approximate the A p massless crossovers M p+2 → M p+1 generated by the perturbing field ϕ 1,3 discussed in Sec. 2. Therefore, in the roaming limit θ 0 → ∞, the staircase model describes a renormalisation group flow that passes by the successive minimal models M p . The final point of the flow is a massive Ising theory. In another work [40], it was shown that the c-function defined by the c-theorem [29,30] using a spectral series in terms of the form factors of the trace of the stress-energy tensor Θ [112,117] presents the same behavior. In addition, it was explicitly demonstrated that the FFs of the of the stress-energy tensor for the A p massless flows can be reconstructed from that of the ShG model. Importantly, for this construction to work, the rapidities in the FFs have to be also shifted by ±kθ 0 /2 with specific integers k. A follow-up publication targeted specifically the A 2 tricritical-critical Ising flow (8), and showed that the form factors of the order and disorder operators along the flow can also be obtained via the roaming limit of the appropriate ShG FFs and, although not published, the correspondence holds for the ε field of the flow as well. As we have said, the staircase model also incorporates the massive Ising field theory, which is regarded in this context as a flow from the critical Ising fixed point to a massive one, and where the consecutive RG flows between the multicritical Ising CFTs terminate. Accordingly, it was demonstrated in [119] that the FFs of the massive Ising theory can be obtained from the ShG FFs by merely taking the limit of Eq. (97), i.e., scaling the rapidity variables withing the FFs. In contrast, for other than the A 2 flow and massive Ising QFT, only the FFs of the field Θ were found to be reproduced by the roaming limiting procedure, and hence the validity of this approach is not a priori obvious and well understood.
Regarding the replicated staircase model, in Ref. [19] the form factors of the standard BPTFs in the sinh-Gordon have been computed up to the four-particle order. While the explicit roaming limit of these form factors was not carried out, they were used in the computation of the conformal dimension of the BPTFs applying the ∆-sum rule [31], which we discuss in more detail in Sec. 7. In particular, it was found that the two-particle contribution correctly reproduces the first 'step' of the staircase, from the critical Ising CFT to the massive theory, while the four-particle one gives the result for the massless flow A 2 from tricritical to critical Ising [19]. This result reveals that the roaming limit also holds for the branch point twist fields of the replicated theory. In the following, we make a step further, by explicitly performing the roaming limit of the form factors of both the standard and the composite BPTFs up to the four-particle order, showing that they reduce to the exact expressions in the A 2 flow (8) obtained via the bootstrap program in Secs. 4 and 5. Proving the correspondence in the first few non-trivial particle levels provides strong evidence that the roaming limit for standard and composite BPTFs is valid for any (composite) BPTF form factors in the A 2 massless flow. In the ShG model, the k-particle form factors of the BPTFs can be parameterised in the usual fashion, that is [19], where each particle is put on the first replica and the superscript τ = 0, µ denotes the standard or the composite BPTF respectively. The minimal form factor for the standard twist field f ShG (θ; n) is given by while the one for the composite field is obtained by including an appropriate monodromy changing factor analogously to what we have done in Eq. (77) for the massless flow. The minimal FF f ShG in Eq. (99) is normalised in such a way that f ShG (±∞, B; n) = 1. The roaming limit construction of the twist field FFs in the massless flow can then be formulated as where we split the rapidities in the sinh-Gordon FF into r right-moving (θ) and l left-moving (θ ′ ) ones, which we shift by θ 0 and −θ 0 respectively. The function B(θ 0 ) is defined in Eq. (97).
In the rest of this section, we explicitly demonstrate the validity of the limit in Eq. (101) up to the four-particle level. Let us first focus on the ShG minimal FFs. Based on Ref. [40], it can be shown that, in the roaming limit (101), the minimal form factor in Eq. (99) reduces to where f RL is the minimal form factor (45) in the massless flow and the normalisation constant N n was found in Eq. (47). Similarly, for the composite twist field one has with f µ RL given by Eq. (80). Note that, according to the definition (101), only the above cases are the relevant limits for the minimal form factor. Some of them involve an exponential factor e ±θ 0 /(2n) but we anticipate that similar factors originate from other terms of the entire FF and they eventually cancel. From Eqs. (102) and (103), it is easy to see that the roaming limit in the two-particle case correctly reproduces the form factors of the massless flow. This is clearer when the two-particle ShG FF is rewritten as a function of the rapidity difference as In the limit (101) this expression reproduces either Eq. (53) (for the standard BPTF) or Eq. (83) (for the composite BPTF) for the 'RR' and 'LL' cases, while it vanishes in the 'RL' case because of the diverging denominator.

Roaming limit of the four-particle FFs of the standard BPTF
It is also not difficult to show that the four-particle 'RRRR' and 'LLLL' FFs are provided by the roaming limit (101). If we consider the standard BPTFs, using as Q T 4 the polynomial determined in [19] and reviewed in App. B.1, we can proceed in the following way. According to Eq. (101), the 'RRRR' or 'LLLL' form factors that only contain right or left movers are given by In this limit, the denominator of the ShG FF (98) does not change but acquires the diverging factor e 6θ 0 /n , whereas for the polynomial Q T 4 we obtain the following lengthy expression which diverges exponentially as e 8θ 0 /n when θ 0 → ∞. In addition, taking into account the limit of the minimal form factor reported in Eq. (102), we have and for the normalisation factor H T n , we find H T n = 2 sin(π/n)ω 2 nf ShG (iπ, B; n) 2 ω 2 −→ e θ 0 /n sin(π/n)ω 2 n sin(π/2n) 2 ω 2 = e θ 0 /n 4 e i 6π n n 2 cos 2 π 2n .
Counting the divergent factors e θ 0 /n in final expressions of Eqs. (106), (107), (108) and (109), we can conlude that the 'RRRR' ('LLLL') roaming limit form factor of the ShG twist field is finite. In fact, putting all the above results together it is straightforward to check that the limit (105) works and Eq. (61) is exactly reproduced.
Turning to the case of the 'RRLL' form factor, we have to consider For the denominator, the limit gives whereas for the polynomial Q T 4 we obtain the following expression which we can rewrite as (1 + ω)y 2 1 y 2 For product of the minimal FFs, we find The limit of the normalisation factor H T n gives the same result as in Eq. (109). Combining (111), (113), (114), and the normalisation (109), it is immediate to see that the divergent exponential factors e θ 0 mutually cancel and that the roaming limit yields Eq. (65), confirming the validity of Eq. (110).

Roaming limit of the four-particle FFs of the composite BPTF
Unlike the four-particle form factor of the standard twist field, the one of the composite twist field was not previously known in the sinh-Gordon theory. In App. B.2, we compute this form factor by constructing and solving the bootstrap equations, starting from the usual ansatz in Eq. (98). Notice that, as we discuss in App. B.2, now the function Q T µ k is not a polynomial but a rational function. At the four-particle level, the explicit expressions of the normalisation H Let us first consider the form factors 'RRRR' and 'LLLL', containing only either rightor left-moving particles. Following Eq. (101), we see that we need to compute the limit of F T µ 4 (θ + θ 0 /2, B(θ 0 ); n). As in Sec. 6.1 for the standard twist field, we study separately this limit for the different terms that constitute the composite ShG form factor in Eq. (98). Applying the limit of the minimal composite form factor reported in Eq. (103), we have where we have rewritten it in terms of x i = e θ i /n . It is convenient to take the limit of the function Q T µ 4 , reported in Eqs. (213), (214), and of the denominator of the ansatz (98) together with the one of the minimal form factor. We find that this limit reproduces the form factor in Eq. (90) up to a normalisation with an exponential e −θ 0 /n Finally, we see that the normalisation term H T µ n , whose explicit expression is given in Eq. (205), becomes cancelling precisely the multiplicative factor in Eq. (116), such that the roaming limit correctly reproduces the 'RRRR' (or 'LLLL') form factor in Eq. (90), as expected. Considering now the 'RRLL' form factor, we can see that it can be obtained from the limit of Eq. (101) in the particular case The joint limit of the denominator of the ansatz (98) and of the polynomial Q T µ 4 reported in Eqs. (213), (214) gives where Q T µ 2,2 flow is the polynomial in the massless flow of Eq. (93). The normalisation H T µ Eq. (103), we have Putting all together, we find that limit of the 'RRLL' form factor is again finite as expected, confirming the validity of the roaming limit also for the composite twist field T µ n .

Standard and symmetry resolved entropies for the massless flow
In this section, we use the form factors computed in the previous sections to study the behaviour of the correlation functions of the standard and composite twist fields. After calculating the running dimension of the field along the renormalisation flow, we investigate the entanglement entropy, comparing it with expected results.

Running dimension from the ∆-sum rule
As we discussed in Sec. 2, the model under examination interpolates between the tricritical Ising CFT M 4 in the UV and the Ising CFT M 3 in the IR, providing the simplest example of a massless renormalisation flow between two A-series diagonal minimal models [35][36][37]. In both fixed points, the properties of the standard twist field T n and the Z 2 composite one T µ n are known from conformal invariance [6], as we reviewed in Sec. 3. In particular, the conformal dimension of the standard twist field is given by Eq. (22) while the dimension of the composite one is in Eqs. (24), (25) for the fixed points of interest. The knowledge of the exact conformal dimensions of the fields in the IR and the UV fixed points of the massless flow provides a non-trivial check of the correctness of the form factors via the ∆-sum rule [31]. Let us start by considering the twist field T n . Along a renormalisation group flow, the difference of conformal dimensions of the field T n in the IR and in the UV is given by an integral of the two-point function between T n and the trace of the stress-energy tensor Θ [31] h UV In order to compute the ∆-sum rule (122) Table 2: Comparison of the difference of the conformal dimensions in the UV and IR fixed points h UV − h IR with the results of the ∆-sum rule, for both the standard twist field T n and the composite one T µ n . The 'CFT' columns collect the exact result fixed by conformal invariance in Eqs. (22), (25), while '∆-sum rule' is the result of the ∆-sum rule truncated at four-particle order, reported in Eq. (125). The column 'n' indicates the number of replicas. We can see that at the four-particle order we already find good agreement for all the number of replicas considered.
which in the case of the (non-replicated) massless tricritical flow have been obtained in [38]. In particular, in a massless model, all the form factors of Θ containing either only left-('L') or right-movers ('R') identically vanish. When considering the replicated theory, we have to take the sum of Θ in each the copy. Therefore, the only non-vanishing form factors are the ones with identical replica indices F Θ|11...1 r,l = F Θ r,l . After integrating out the distance t in the spectral expansion of the ∆-sum rule (122), we finally find [8,31] where E is the energy (reported in Eq. (10) for a massless model).
The leading non-trivial form factor of Θ is the four-particle 'RRLL' one, coupling two rightand two left-movers [38] where γ is Euler-Mascheroni's constant and f RL (θ) = f RL (θ; n = 1) is the minimal form factor in Eq. (45) for a single replica n = 1. Since all form factors have an even number of left-and right-moving particles, Eq. (124) is the only contribution at the four-particle level [38]. We can then consider the approximation where F Θ 2,2 is given in Eq. (124) and F T |1111 2,2 is the twist field FF that we obtained in Eqs. (65), (66). Analogous expressions hold for the ∆-sum rule of the composite twist field T µ n , replacing F T |1111 2,2 with the form factor F T µ |1111 2,2 of the composite field reported in Eqs. (91), (92), (93). In Table 2, we compare the exact difference of conformal dimensions of both the standard and the composite twist fields with the result of the ∆-sum rule at the four-particle order (125) for n = 2, 3, 4 replicas. The integral in Eq. (125) has been computed numerically using the Divonne routine of the library Cuba [124] for the software Mathematica, using a cut-off θ j ∈ [−60, 60] for the rapidities. Already at the four-particle order we find a good agreement between the exact CFT result and the ∆-sum rule, confirming the correctness of the form factors computed in Sec. 4 and the relatively small weight carried by the higher order FFs. This is consistent with Ref. [19], where it was found that, for the staircase model (reviewed in Sec. 6), the four-particle contribution obtained in the roaming limit reproduces the difference in conformal dimensions of the standard twist field along the massless flow (8).
The ∆-sum rule (122) can be modified to give a running dimension of the (composite) twist fields along the flow [19,89] where now the integral over the distance t starts from a finite length ℓ. As we did before in Eq. (123), we expand the two-point function in form factors and we integrate over the distance t, obtaining [19,89] where again, in the massless flow, the leading contribution is given by the 'RRLL' form factors, A running ∆-theorem (126) can also be formulated for the composite twist field by considering the appropriate form factor F T µ |1111 2,2 of that operator obtained in Eqs. (91)- (93). In Ref. [14], it was argued that the running dimension h(ℓ) of the branch point twist field T n provides an entropic c-function which is monotonically decreasing along the flow. In Fig. 1a we report the result of the numerical integration of the running Delta theorem in Eq. (128) for the standard branch point twist field T n taking n = 2, 3, 4 replicas. We observe that, already at the four-particle order, the running dimension monotonically decreases with ℓ for all the number of replicas considered. In Fig. 1b, we plot the running dimension of the composite twist field T µ n at four-particle order. In particular for n = 2 replicas, the dimensions of the twist fields and of the charge operator conspire to give the same ultraviolet and infrared conformal dimensions for the composite twist field. Remarkably, we see that along the flow the running dimension varies and is not monotonic in ℓ, differently from the standard twist field. In the inset, we zoom in the region of small running dimension which shows that also for larger number of replicas n the behaviour is non-monotonic.

Cumulant expansion of the entanglement entropy
As a main result of this paper, in this section we discuss the form factor expansion of the entanglement entropy along the massless renormalisation group flow. As we will show, the formal expressions require a suitable regularisation, after which the form factors containing particles with the same chirality reproduce the logarithmic entanglement entropy of the infrared Ising CFT, while those that include particles of different chirality provide the corrections along the flow.  Figure 1: Semi-logarithmic plot of the running conformal dimension h T τ obtained using the ∆-sum rule at four-particle order reported in Eq. (128). In the plot on the left, we show the result for the standard twist field T n , while on the right for the composite one T µ n . The dotted gray lines indicate the exact difference between the UV and IR conformal dimension obtained from Eq. (22) in the plot on the left and from Eqs. (25) in the plot on the right. As expected, for small distances ℓ, the running conformal dimension approaches the exact UV results, as also reported in Table 2. On the left, we see that the one of the standard twist field decreases monotonically, consistently with its behaviour as an entropic c-function. On the right, in the inset we zoom on the running dimension of the composite twist field (for n = 2 we have . The running dimension of the composite field is not monotonic along the flow. Instead of studying the correlator of the twist field, we find more convenient to directly apply its form factor expansion to the Rényi entanglement entropies defined in Eq. (1). Plugging in Eq. (20) the spectral series (31) of the twist field correlator and expanding the logarithm for the Rényi entropy order by order in the number of particles, we obtain the cumulant expansion [18,125] where, in analogy with Ref. [125], we have introduced the cumulants  , obtained by subtracting all the possible clusterisations for large rapidities [18,125]. Recall from the discussion around Eqs. (38), (39) that, due to the clustering property, at large rapidity differences between the particles, the form factor factorises in the product of form factors with less particles. For example, up to the four-particle level, the connected components take the form f (131) By definition, the connected form factors f T |j 1 ... r,l vanish for large rapidities. As we will see, this improves the convergence of the integral in Eq. (130).
In the expansion (129), we recognise two different kinds of cumulants. Those containing only form factors diagonal in the chiralities, c T r,0 , c T 0,l , which we will call non-interacting cumulants, and the ones that couple left-and right-movers, which we will call interacting. In the rest of the section, we treat the two kinds of terms separately since, as we will see, they give different contributions to the entanglement entropy (129).

Non-interacting cumulants
Let us first focus on the non-interacting cumulants. As we saw in Sec. 4, in the massless flow, the form factors containing either only right-or only left-movers are identical to those of the massive Ising theory except for the vacuum expectation value ⟨T n ⟩, implying that their connected components are identical f Given this identity, we can analyse them by applying the same strategy as in Ref. [18] for the massive Ising theory, which we also report in App. C.
Using the Pfaffian structure of the form factors in Eq. (59), it was shown that the noninteracting cumulants c T r,0 have the general expression [10,17,18,62] c T k,0 (M ℓ; n) = n k (2π) k n j 2 ,...,j k =1 where θ ij = θ i − θ j , we have summed over j 1 , and we have introduced the notation From the form of Eq. (134), with all terms cyclically connected [18,62], it is clear why they are known as fully connected. Importantly, Eq. (134) holds for both the massless tricritical-critical flow and for the massive Ising theory. The only difference between the cumulants in these two models is the form of the energy E appearing in the exponential factor. This difference has however a major effect in the integral in Eq. (134). For simplicity, we can start by analysing the two-right-mover cumulant c T 2,0 . The generalisation to higher particles will be straightforward. In our massless flow, as already recognised in Ref. [38] for a different correlation function, the two-fold integral in Eq. (134) is IR divergent due to the absence of a mass gap. In fact, in the relative and center-of-mass coordinates, θ 12 = θ 1 − θ 2 and A = (θ 1 + θ 2 )/2, the energy (10) of two right-moving particles takes the form For right-movers, the IR region E → 0 corresponds to large and negative center-of-mass rapidity A → −∞. Since the form factors do not depend on A, we see that the integrand of Eq. (134) tends to a non-zero constant for A → −∞, leading to a divergence when the integral in A is performed.
In order to cure this IR divergence, we introduce a cut-off Λ in the center-of-mass rapidity A. Since the form factors do not depend on A, the resulting integral can be cast in terms of the exponential integral function Ei(x), We see that Λ plays the role of a cut-off at large distances with M ℓ ≪ Λ. In this limit, using the expansion we obtain a logarithmic dependence in the interval length ℓ, where z 2 is the function Remarkably, up to an additive constant and the large distance cut-off Λ, the sum of the left-and right-moving two-particle cumulants in our massless flow, c T 2,0 + c T 0,2 , is equal to the UV limit of the two-particle cumulant of the massive Ising model (cf. Eq. (221) in App. C). This is consistent with the expectation that, in the IR, the contributions of the interacting cumulants vanish because the flow leads to the critical Ising fixed point. As such, we expect that for large distances the non-interacting cumulants completely reproduce the logarithmic entanglement entropy of the Ising CFT.
Moving to higher particle cumulants c T r,0 , we expect a similar structure. In the presence of more than two particles, a convenient set of coordinates is again provided by the center-of-mass rapidity A = 1 k k j=1 θ j and the difference between the rapidities θ j,j+1 = θ j − θ j+1 , with Jacobian equal to one. For convenience, we further define the rapidities in the center-of-mass frame of reference which can be shown to depend only on the rapidity differences. In the massless flow, the r-fold integral in Eq. (134) is divergent in the large negative center-of-mass rapidity region A → −∞.
It is important to stress a subtle point. Due to the clustering property (see Eqs. (38) and (39)), the integral of the form factor is divergent in the direction of the sum of any two rapidities θ j . However, in the cumulants, the non-connected factorised component is subtracted as in, e.g., Eqs. (132) and (133), guaranteeing that the integral of the connected part converges in those directions. The only remaining divergence is the one in the direction of large negative center-of-mass A, as it happens for the two-particle cumulant (137).
In the center-of-mass coordinates defined before Eq. (141), the energy of r right-moving particles in the massless flow takes the form As already done in Eq. (137) for the two-particle case, we again introduce a cut-off Λ on the large negative center-of-mass rapidity A and we write the integral over it in terms of the exponential integral function Ei(x), In the large cut-off limit Λ ≫ M ℓ, we can approximate the cumulant using the expansion of the exponential integral in Eq. (138), obtaining the expected logarithmic behaviour with z k (n) equal to w(−θ 2l,2l+1 + 2πi (j 2l − j 2l+1 )) w(θ 2l+1,2l+2 + 2πi (j 2l+1 − j 2l+2 )) .
As happens with the two-particle cumulants, the sum c T r,0 +c T 0,r in the massless flow coincides with the UV limit of the r-particle cumulant of the massive Ising theory up to additive constants (see Eq. (225) in App. C). In Ref. [18], the resummation of the z r (n) terms is carried out. Taking Eq. (144) and applying their result, we find that r even c T r,0 (M ℓ, Λ; n) + l even This shows that, up to additive constants, the sum of the non-interacting left-and right-movers contribution to the entanglement entropy (129) in the massless flow gives the entropy of the Ising CFT in the IR fixed point.

Interacting cumulants
As shown in the previous discussion, the non-interacting cumulants contribute to the entropy of the IR Ising CFT; hence the corrections for smaller distances are provided by the interacting cumulants c T r,l , which couple left-and right-movers. In this section, we study the only interacting cumulant at four-particle level, namely c T 2,2 c T 2,2 (M ℓ; n) = n 2 × 2 (2π) 4 where the 'RRLL' form factor F T 2,2 is given in Eq. (65). As for the non-interacting cumulants, the integral of F T 2,2 is divergent in the IR limit. However, unlike the previous discussion, now the subtraction of the clusterisation in the connected component f T |1j 2 j 3 j 4 2,2 ensures that the cumulant is convergent and a regularisation is not needed.
In Fig. 2, we report the result of the numerical integration of the cumulant c T 2,2 for different values of M ℓ and for n = 2, 3 replicas, performed using the Divonne routine of the library Cuba [124]. In the UV region M ℓ ≪ 1, we expect a leading logarithmic behaviour, since the sum of the interacting and non-interacting form factors should reproduce the logarithmic entanglement entropy of the tricritical Ising UV fixed point. In Fig. 2a, we plot the interacting cumulant c T 2,2 in Eq. (148) for n = 2, 3 replicas and we compare it with the fit of the numerical points to a logarithmic function −α n log M ℓ + C n . We perform the fit for M ℓ ⩽ 2 × 10 −4 when n = 2 and for M ℓ ⩽ 6 × 10 −4 when n = 3, obtaining the parameters From Fig. 2a, we see that for M ℓ ≪ 1 the cumulant is in good agreement with the expected logarithmic behaviour.
To understand the behaviour in the IR (M ℓ ≫ 1), recall from the general introduction of Sec. 2 that, near the IR fixed point, the effective theory describing the massless flow is the TT deformation of the critical Ising CFT, as shown in Eq. (9). The entanglement entropies of generic TT -deformed CFTs have been heavily studied in recent years, see e.g. [95][96][97][98][99][100][101][102][103][104][105][106][107][108]. In particular, in Ref. [97], the entropy of an interval of length ℓ in a system of finite size L has been computed perturbatively for a generic TT -deformed CFT. Let the deformed action be where the TT coupling g has dimensions of an inverse mass squared, g ∝ M −2 . The first perturbative correction to the Rényi entanglement entropy of the IR CFT was found to be [97] (1 − n) δS (1) n (ℓ, L, g) = − where ϵ is a non-universal UV cut-off. Here we are interested in the entanglement entropy in the thermodynamic limit L → ∞ of Eq. (152), Comparing the effective Lagrangian (9) with the generic one in Eq. (151) and taking into account that in our case T = − 1 2 ψ∂ψ andT = − 1 2ψ∂ψ , we can conclude that g = − 4 π 2 M 2 . Therefore, since the central charge of our IR point is c = 1 2 , Eq. (153) specialised to our massless flow gives Observe that in the prediction of Eq. (154) the leading correction is of the form A n ℓ −2 + B n ℓ −2 log ℓ. The coefficient A n is not universal due to the presence of the U V cutoff ϵ, while the factor B n is. In particular, for n = 2, 3 replicas, its numerical value is = 0.02096 . . . , for n = 3.
Note also that the leading correction in Eqs. (152), (153), (154) is non-zero only for for n ⩾ 2 Rényi entropies, while it vanishes in the replica limit n → 1 [97]. It is worthwhile to compare the first-order perturbative prediction in Eq. (154) with the leading correction that we obtain here from the form factor cumulant expansion (129), which is given by the interacting cumulant c T 2,2 . In Fig. 2b we study this cumulant for n = 2, 3 replicas as a function of M ℓ and we perform a best fit of the numerical points to a function A n ℓ −2 + B n ℓ −2 log ℓ, for M ℓ ⩾ 50 when n = 2 and for M ℓ ⩾ 20 when n = 3. We obtain For large distances, we find a good qualitative agreement with the functional form predicted in Eq. (154). However, while for n = 2 replicas the numerical result of the fit for B 2 is close to the predicted value in Eq. (155), this is not the case for n = 3. A possible explanation of this discrepancy is that the higher-particle interacting cumulants c T r,l , which we are neglecting, also contribute to the term log ℓ/ℓ 2 and their contribution depends on the number of replicas n.

Entanglement entropy
Finally, we can put together the results obtained in Secs. 7.2.1 and 7.2.2 to get the total entanglement. In Fig. 3, we consider the Rényi entanglement entropy in the massless flow as a function of M ℓ for n = 2 and 3. The results in this figure (represented as symbols) have been obtained with the cumulant expansion (129), including the first 30 non-interacting cumulants c T r,0 and c T 0,l and the leading interacting cumulant c T 2,2 . In the plot, we also report the expected behaviour of the entanglement entropy when approaching the UV (dashed curves) and IR (dotted curves) fixed points. When approaching the IR, M ℓ ≫ 1, we find a very good agreement between the truncated cumulant expansion and the expected IR asymptotics for both values of n. In the UV, M ℓ ≪ 1, while the truncated expansion presents a behaviour compatible with a logarithmic divergence in M ℓ, it does not quantitatively agree with the UV entropy. This is expected, since the UV limit is the regime where the higher-particle interacting cumulants c T r,l contribute more significantly and here we are only considering the four-particle one, c T 2,2 . It would be interesting to include higher order terms, but this is a challenging task due to the difficulty of computing higher-particle cumulants, which involve an increasing number of multidimensional integrals. Nevertheless, in the light of Fig. 3, we can conclude that only including the leading interacting cumulant is enough to qualitatively observe the crossover between the IR and UV regimes.
Before concluding this section, let us comment on the symmetry resolved entanglement entropy. Also in this case, a cumulant expansion analogous to Eq. (129) holds true, by replacing the form factors with appropriate ones for the composite twist field that we have determined in Sec. 5. In particular, the expansion of the symmetry resolved entanglement entropy contains both the non-interacting cumulants reproducing the Ising CFT and the interacting ones providing the corrections, analogously to what happen for the standard entropy. The symmetry resolved entropy in the massive Ising model has been recently studied in Ref. [62] where it was found that (differently from what happens for the standard BPTF) the cumulants of the composite twist fields are divergent although the theory is massive; consequently a regularisation was required. This fact suggests that also for the massless flow, the regularisation employed for the total entropy is not sufficient to find a finite result for the composite twist field. Resolving such a regularisation is a problem that goes beyond the scope of this paper and we hope to return on the issue in the future.

Conclusions
In this paper, we investigated the ground state Rényi entanglement entropies of a single interval in the massless QFT associated to the renormalisation group flow connecting the tricritical and critical Ising CFTs. The corresponding two-point correlation function of branch points twist fields admits a form factor expansion along the flow. We showed that these form factors can be calculated in two different and independent ways. On the one hand, we have directly applied the bootstrap approach of Ref. [8] for massive integrable QFTs: based on the symmetries of the theory and the exchange properties of the twist fields, we obtained a set of equations for the form factors and we found a general ansatz that solves them. Alternatively, we obtained the form factors using the Zamolodchikov's staircase model, an extension of the sinh-Gordon theory with complex couplings that includes the tricritical-critical renormalisation flow. In this framework, the form factors of several fields in this massless flow have been obtained as the roaming limit of those in the sinh-Gordon theory. We showed that the same strategy works for the branch points twist fields; we derived explicit expressions for the two and four-particle form factors, from which the higher-particle ones can be recursively derived. Obviously the two approaches gave identical results.
The form factor expansion of the entanglement entropy can be rearranged order by order in the number of particles in terms of cumulants, which are given by the connected part of the form factors. In this cumulant expansion, we distinguished free and interacting cumulants. The former only contain particles of the same chirality and give the entropy of the IR Ising CFT. In fact, we found that, after a proper regularisation, they are equal to those that appear in the massive Ising theory. On the other hand, the interacting cumulants, which contain particles with different chiralities, describe the behaviour of the entanglement entropy along the flow. In particular, we checked that the lowest-particle interacting cumulant yields in the UV limit the expected logarithmic behaviour in the subsystem size. The IR limit can be described by a TT deformation of the Ising CFT. We showed that the lowest-particle interacting cumulant expansion approaching the IR point qualitatively reproduces the result for a generic TT perturbation at first order [97].
The massless flow (8) that we studied here is also connected with the SU (3) 2 -homogeneous sine-Gordon (HSG) model [87]. As shown in Refs. [19,88,89], along the renormalisation group flow, the central charge and the twist field dimension of this theory present two plateaux, analogously to the behaviour of the staircase model considered here. For certain values of the parameters, it was shown that one of these plateaux corresponds to the massless flow from tricritical to critical Ising [19,88]. Since the form factors of the standard twist field in the SU (3) 2 -HSG model have been obtained in Ref. [19] up to the four-particle order, it would be interesting to recover our results from an appropriate limit of the HSG expressions.
The massless flow connecting the tricritical and critical Ising CFTs enjoys a global Z 2 symmetry. In this work, we also considered the composite branch point twist fields associated to this symmetry. Their two-point functions give the charged moments of the reduced density matrix from which one can determine the symmetry-resolved entanglement entropy. Similarly to the standard twist fields, we obtained their bootstrap equations, which now include the non-trivial monodromy due to the insertion of the charge, and we found a general ansatz for their solution, which allows to obtain the higher-particle form factors recursively. We further derived them as the roaming limit of the composite twist field form factors of the sinh-Gordon theory. Remarkably, the latter were neither known in the literature, a gap which we also filled here.
Our goal now is to identify the four-particle 'RRLL' form factors using the well-known two-particle quantities. For simplicity we place every particle on the first replica and specify Eq. (159) to the case of interest Applying now the residue axiom in Eq. (34) to the ansatz (160) we can derive recursive equations for the H T 2,2 normalisation factors and the Q T 2,2 functions. Let us first also recall that the minimal form factor f RL satisfies the identity The residue of the denominator of the ansatz (160) takes the form − i Res from which we can obtain the residue of the entire expression (160) as − i Res where we used Eqs. (161) and (162). Via algebraic manipulations we can simplify the above formula to − i Res Following the residue axiom in Eq. (34), the residue of the kinematical pole in Eq. (164) has to reproduce the two-particle form factor in Eq. (53). We first recast it in the shape of our ansatz as where we have defined Q T 0,2 (y 1 , y 2 ) = σ 2 (y 1 , y 2 ) = y 1 y 2 .
Comparing the residue in Eq. (164) with the two-particle form factor in Eq. (165) leads to the recursion equations for H T 2,2 and for the polynomial Q T 2,2 which we separate as Q T 2,2 (ωx, x, y 1 , y 2 ; n) = (ω 1/2 xy 1 − 1)(ω 1/2 xy 2 − 1) = 1 − ω 1/2 x (y 1 + y 2 ) + ωx 2 y 1 y 2 . (170) We postulate a solution to Eq. (170) completely symmetrical in x 1 , x 2 , y 1 , y 2 and hence writing which is the most general expression compatible with the fact that the form factor has zero Lorentz spin, and that sending each rapidity to ±∞ the entire FF converges to zero. Posing x 1 = ωx, x 2 = x, we get a unique solution for the unknown constant A, namely (172) This means that the entire solution can be written as which we can also rewrite as

A.2 Form factors of the Z 2 -composite BPTF
Once again, we start our derivation by recalling and repeating the ansatz for Z 2 -composite BPTF where we have r right-mover and l left-mover particles and again x i = e θ i /n and y i = e −θ ′ i /n . The cyclic permutation and the exchange axioms are already satisfied since is fulfilled via f µ γγ (θ; n) = 2 cosh(θ/(2n))f γγ (θ; n) = sinh(θ/n) . (177) for γ ′ different from γ, we similarly satisfy In full analogy to what we have done in the previous section for the standard twist field, we apply the residue axiom Eq. (69) to the ansatz (175) in order to derive recursive equations for the normalisation factors H T µ r,l and the Q T µ r,l functions. Since the denominator of the ansatz (175) is the same as the one for the standard twist fields in Eq. (175), we can reuse the same residue we have computed in Eq. (162). Using again the property (161) of the minimal form factor, the residue of the ansatz with 4 particles yields − i Res Again, from the residue axiom (69), this expression must be compared with the two-particle FF, which we can rewrite as where y 2 2 − y 2 1 /(2y 1 y 2 ) = sinh((θ ′ 1 − θ ′ 2 )/n). We then end up with the equation for Q T µ 2,2 as well as the normalisation which we can separate as Notice that, differently from what happened in Eq. (170) for the standard twist field, now the function Q T µ 2,2 is not a polynomial but a rational function. We write the solution to Eq. (184) as which is the most general expression compatible with (i) the form factor has zero Lorentz spin, and (ii) sending each rapidity to ±∞ the entire FF converges to a constant. When setting x 1 = ωx and x 2 = x we can obtain the same solution for A as for the case of the standard BPTF, namely A = − 1 2 cos π 2n and hence The ansatz with the above fraction of polynomial Q T µ 2,2 (ωx 1 , x 2 , y 1 , y 2 ; n) (0) satisfies all the FF axioms. Notice that while the n → 1 limit of the standard BPTF is not well defined, the FFs of the composite twist field reduce to those of the disorder field 1 + x 1 x 2 y 1 y 2 As pointed out in the main text, the solution of the bootstrap equation is in general not unique, since we can often add to our polynomial Q also a (non-trivial) kernel solution, that is, another polynomial (or fraction of polynomials) Q (k) which satisfies the homogeneous equation Polynomial kernel solutions at the two-and four-particle level have been identified in [19]. In particular, the two-particle kernel solution reads as from which the required four-particle kernel solution for the flow can be constructed by squaring the expression due to the the anticipated symmetry between the variables of the RR and LL particles. Based on the above consideration, we can write the eventual kernel as Q T µ 2,2 (x 1 , x 2 , y 1 , y 2 ; n) (k) = ω 2 n 2 8 cos 3 π 2n x 1 x 2 − 1 4 sec 2 π 2n (x 1 + x 2 ) 2 (x 1 x 2 )(x 1 + x 2 )(y 1 + y 2 ) × × y 1 y 2 − 1 4 sec 2 π 2n (y 1 + y 2 ) 2 y 1 y 2 (y 1 + y 2 ) (190) that is, taking the product of (189) and additionally, by also renormalising the expression with (x 1 x 2 y 1 y 2 )(x 1 + x 2 )(y 1 + y 2 ) which does not spoil the kernel property. We chose the pre-factor in a way that the entire expression for the polynomial Q T µ 2,2 (x 1 , x 2 , y 1 , y 2 ; n) = Q T µ 2,2 (x 1 , x 2 , y 1 , y 2 ; n) (0) + Q T µ 2,2 (x 1 , x 2 , y 1 , y 2 ; n) (k) gives (1 + x 1 x 2 y 1 y 2 )/(x 1 x 2 y 1 y 2 ) in the n → 1 limit, which reproduces Q µ 2,2 . The normalisation factors match as well, since H µ 2,2 = −4N 4 1 = 2 e −4G/π .

B Form factor bootstrap for branch point twist fields in the sinh-Gordon model
In this appendix, we first report the known results for the four-particle form factor of the standard twist field in the sinh-Gordon model and we then derive the previously unknown form factor of the composite one.

B.2 Form factors of the Z 2 -composite BPTF
As we mentioned in the main text, differently from the the four-particle form factor of the standard twist field T n , in the sinh-Gordon model the one of the composite field T µ n was not previously known in the literature. In this appendix we compute this form factor by constructing and solving the bootstrap equations, in full analogy to what we have done in Sec. 5 and in App. A.2 in the case of the massless flow.
Plugging the function Q T µ 4 in Eqs. (212)-(214) and the normalisation H T µ 4 from Eq. (205) in the ansatz (196) we finally obtain the four-particle form factor for the composite twist field that we used in Sec. 6.2.

C Cumulant expansion of the entanglement entropy in the massive Ising theory
In this appendix, we review the known results for the form factor expansion of the entanglement entropy in the massive Ising model, obtained in Refs. [18]. In particular, we find a direct relation between UV limit of the cumulant expansion of the entropy in the massive Ising and the non-interacting part of the expansion in the massless flow, studied in Sec. 7.2.1.
In the massive Ising theory, if we denote by m the mass gap, the ground state Rényi entanglement entropy admits the following cumulant expansion [18], These cumulants can be reexpressed as in Eq. (134) and, therefore, the k-particle cumulant c T k, Ising is similar to the k-rightor k-left-mover non-interacting cumulants c T k,0 , c T 0,k in the massless flow (130), differing only in the energy E in the exponential factor. In the massive Ising theory, the energy of k particles is E(θ 1 , . . . , θ k ) = k i=1 m cosh(θ i ) .

(219)
As shown in Ref. [8], in the UV limit mℓ ≪ 1, the expansion of the Bessel function K 0 reproduces the expected UV logarithmic behaviour of the entanglement entropy up to an additive constant [8] c T 2, Ising (mℓ; n) ≈ mℓ≪1 −z 2 (n) log mℓ + const, where the function z 2 (n) was introduced in Eq. (140). We can now investigate the higher-particle cumulants c T k, Ising . If we write them in terms of the center-of-mass coordinates, we can apply the integral identity and the fact that the form factors only depend on the relative rapidities to integrate out the center-of-mass rapidity A. We then obtain and ξ j are defined in Eq. (141). In the UV limit mℓ ≪ 1, by expanding the Bessel function using Eq. (220), we get at leading order the logarithmic behaviour of the entropy expected in the Ising CFT up to an additive constant c T k, Ising (mℓ; n) ≈ mℓ≪1 −z k (n) log(mℓ) + const, where the coefficient z k is the same as in Eq. (145). Comparing Eq. (225) with the analogous formula in Eq. (144), we can immediately see that the UV limit of the k-particle massive Ising cumulants is twice the k-right-mover cumulants of our massless flow. Notice that the factor 2 comes from the expansion of the Bessel function and ultimately its origin is the difference in the energy of the two models. Before concluding this appendix, let us make a remark on the computation of the coefficients z k (n). The expression in Eq. (145) contains k − 1 integrals and, therefore, it is not practical for numerical calculations. In Ref. [18], the analytic continuation of Eq. (145) was carried out for n ⩾ 1 replicas, writing z k (n) as a single integral for any k (see also [17,62]) where, for k = 2p, for p even, and w(θ; n) is given in Eq. (135). Eq. (226) is efficient for numerical calculations. We employed it to compute the first 30 non-interacting cumulants in the truncated expansion of the entropies of the massless flow plotted in Fig. 3.