Entanglement entropies of minimal models from null-vectors

We present a new method to compute R\'enyi entropies in one-dimensional critical systems. The null-vector conditions on the twist fields in the cyclic orbifold allow us to derive a differential equation for their correlation functions. The latter are then determined by standard bootstrap techniques. We apply this method to the calculation of various R\'enyi entropies in the non-unitary Yang-Lee model.


Introduction
Ideas coming from quantum information theory have provided invaluable insights and powerful tools for quantum many-body systems. One of the most basic tools in the arsenal of quantum information theory is (entanglement) entropy [1]. Upon partitioning a system into two subsystems, and , the entanglement entropy is defined as the von Neumann entropy S( ) = −Tr ρ log ρ , with ρ being the reduced density matrix of subsystem .
Entanglement entropy (EE) is a versatile tool. For a gapped system in any dimension, the entanglement entropy behaves similarly to the black hole entropy : its leading term grows like the area of the boundary between two subsystems instead of their volume, in a behavior known as the area law [2][3][4][5][6][7][8] : where α is a non-universal quantity. Quantum entanglement -and in particular the area law -has led in recent years to a major breakthrough in our understanding of quantum systems, and to the development of remarkably efficient analytical and numerical tools. These methods, dubbed tensor network methods, have just begun to be applied to strongly correlated systems with unprecedented success [9][10][11].
For critical systems, a striking result is the universal scaling of the EE in one-dimension [12][13][14]. For an infinite system, with the subsystem being a single interval of length , one has where c is the central charge of the underlying Conformal Field Theory (CFT). This result is based on a CFT approach to entanglement entropy combined with the replica trick, which maps the (Rényi) EE to the partition function on an N -sheeted Riemann surface with conical singularities. In some particular cases -essentially for free theories -it is possible to directly calculate this partition function [15][16][17][18][19][20] using the general results from the 1980's for free bosonic partition functions on Riemann surfaces [21][22][23][24]. In most cases however this is very difficult. An alternative approach is to replicate the CFT rather than the underlying Riemann surface. Within this scheme one ends up with the tensor product of N copies of the original CFT modded out by cyclic permutations, and the conical singularities are mapped to twist fields, denoted as τ. These theories are known as cyclic orbifolds [25][26][27][28]. Within this framework the Rényi EE boils down to a correlation function of twist fields in the cyclic orbifold [29,30]. The case of a single interval is particularly simple as it maps to a two-twist correlation function. When the subsystem is the union of m > 1 disjoint intervals most results are restricted to free theories [16][17][18][19][20]31], and much less is known in general [32]. In the orbifold framework, this maps to a 2m-twist correlation function, which is of course much more involved to compute than a simple two-point function.
In this article we report on a new method to compute twist fields correlation functions.
Our key ingredients are (i) the null-vector conditions obeyed by the twist fields under the extended algebra of the cyclic orbifold and (ii) the Ward identities obeyed by the currents in this extended algebra. Note that the null-vector conditions for twist fields were already detected in [26], but until now then they have only been exploited to determine their conformal dimension. Our method is quite generic, the only requirement being that the underlying CFT be rational (which in turn ensures that the induced cyclic orbifold is rational). This approach provides a rather versatile and powerful tool to compute the EE that is applicable to a variety of situations, such as non-unitary CFT, EE of multiple intervals, EE at finite temperature and finite size, and/or EE in an excited state.
We illustrate this method with the most basic minimal model of CFT: the Yang-Lee model. This model has only two primary fields: the identity 1 and the field φ. However, the simplicity of this situation -in particular, the nice form of the null-vector conditions obeyed by the identity operator -comes with a slight complication: the model is not unitary, and φ has a negative dimension. Hence, the vacuum |0〉 and the ground state |φ〉 are distinct (i.e. the vacuum is not the state with lowest energy), which implies that the ground state breaks conformal invariance. This leads to an important modification in the path integral description used in the replica trick : the boundary conditions at every puncture must reflect the insertion of the field φ (and not the identity operator). In practice, this means that the twist field τ must be replaced by τ φ ∝ :τφ: as noted in [33], but also that the correlation functions of these twist fields must be evaluated in the ground state |φ〉 rather than in the vacuum |0〉.
Hence we see that, for the Yang-Lee model at finite size, even the single-interval entropy requires the computation of a four-point function, and this is where the full power of nullvector equations can be brought to bear.
The plan of this article is as follows. In Sec. 2 we review the cyclic orbifold construction of [26][27][28], and its relation to Rényi entropies. In Sec. 4 we describe a basic example where the null-vector conditions on the twist field only involve the usual Virasoro modes, and thus yield straightforwardly a differential equation for the twist correlation function. In Sec. 5, we turn to more generic situations, where the null-vector conditions involve fractional modes of the orbifold Virasoro algebra: we first introduce the Ward identities for the conserved currents T (r) (z) in the cyclic orbifold, and use them to derive differential equations for a number of new twist correlation functions. Finally, in Sec. 6 we describe a lattice implementation of the twist fields in the lattice discretisation of the minimal models, namely the critical Restricted Solid-On-Solid (RSOS) models. We conclude with a numerical check of our analytical results for various EEs in the Yang-Lee model.

Entanglement entropy and conformal mappings
Consider a critical one-dimensional quantum system (a spin chain for example), described by a conformal field theory (CFT). Suppose that the system is separated into two parts : and its complement . The amount of entanglement between and is usually measured through the Von Neumann entropy. If the system is in a normalised pure state |ψ〉, with density matrix ρ = |ψ〉〈ψ|, its Von Neumann entropy is defined as: The Rényi entropy is a slight generalisation, which depends on a real parameter N : In the limit N → 1, one recovers the von Neumann entropy: For integer N , a replica method to compute this entropy was developed in [14] (see [29] for a recent review). The main idea consists in re-expressing geometrically the problem. The partial trace ρ acts on states living in and propagates them, while tracing over the states in . It can be seen as the density matrix of a "sewn" system kept open along but closed on itself elsewhere.
When is a single interval ( = [u, v]), the resulting Riemann surface is conformally equivalent to the sphere. It can be unfolded (mapped to the sphere) using a change of variable of the form: When |ψ〉 = |vac〉 is the vacuum state of the CFT, this change of coordinates allows [29] to compute the entropy of a single interval in an infinite system, with the well-known result: Throughout this paper, we shall rather consider the case of an interval = [0, ] in a finite system of length L with periodic boundary conditions. In this case, one has [29]: where we have slightly changed the notation to indicate that the total system is of finite size L.
This type of calculations becomes more complicated for the entropy of other states than the vacuum: two operators then need to be added on each of the sheets of Σ N . This is one of the main limitations of the method based on conformal mapping : a lot of the structure of the initial problem disappears after the conformal map. In this case, a one-variable problem (the size of the interval) becomes a 2N -variable problem. These complicated correlation functions have only been computed for free theories [34,35], and have been used in various contexts since then [36][37][38]. Moreover, if Σ is the initial surface where the system lives, then the genus of Σ N is g(Σ N ) = N g(Σ) + (N − 1)(p − 1), where p is the number of connected components of . Hence if is not connected or if the initial surface is not the Riemann sphere, one has to deal with CFT on higher-genus surfaces.

Correlation functions of twisted operators
The Rényi entropies can alternatively be interpreted as correlation functions of twist operators. We consider a system of finite length L with periodic boundary conditions, in the quantum state |ψ〉. In the scaling limit, this corresponds to a CFT on the infinite cylinder of circumference L, with boundary conditions specified by the state ψ on both ends of the cylinder. After the conformal mapping z → exp 2πz L , one recovers the plane geometry. The Rényi entropy of a single interval = [0, ] in the pure state |ψ〉 is given as a correlation function in the N orbifold CFT (see below): where Ψ = ψ ⊗N corresponds to N replicas of the operator ψ at a given point. The twist operators τ and τ implement the branch points a the ends of the interval. Since these branch points introduce singularities in the metric, one has to choose a particular regularisation of the theory at each branch point: each choice of regularisation corresponds to a choice of primary twist operator τ. The classification of primary twist operators is obtained by the induction procedure (see Sec. 3.3), which uniquely associated any primary operator φ of the mother CFT to a twist operator τ φ , with dimension h φ = (N − 1/N )c/24 + h φ /N . In a unitary CFT, the most relevant operator is the identity (i.e. the conformally invariant operator), and the correct choice for the twist operator in (6) is τ = τ 1 . In Sec. 6 we introduce the construction of a lattice regularisation scaling to any given primary twist operator τ φ in a minimal model of CFT. More generally, if is a union of p ≥ 1 disjoint intervals: then one may define the p-interval correlation function: with and any choice of twist operators (τ 1 , . . . τ p ) and ( τ 1 , . . . τ p ).

Non-unitary models
Although the goal of the present paper is not to study specifically entanglement in non-unitary models, some emphasis is put on the Yang-Lee singularity model. The reason for this is that the corresponding minimal model has the simplest operator algebra (it has only two primary fields), which makes calculations more tractable and easy to present. However it should be stressed that what we are computing are partition functions on N -sheeted surfaces. For a unitary model this corresponds to Rényi entropies, and for that reason we chose to refer to these partition functions as "entropies" even in the non-unitary case. This is just a matter of terminology, and we do not claim that they provide a good measure of the amount of entanglement. The problem of entanglement entropy in non-unitary models has already been addressed in various contexts [33,[39][40][41]. For comparison with the existing literature on the subject, we clarify in this section the specific choices and observations that we made for non-unitary models. We refrain from using the bra/ket notations to avoid any possible source of confusion.
Consider a Hamiltonian H acting on a vector space E. The transpose operator t H acts in the dual space (consisting of all linear forms) E * as for any linear form w. We assume that H is diagonalizable with a discrete spectrum and eigenbasis {r j } H r j = E j r j .
The dual basis {w j }, which is defined by w i (r j ) = δ i j , is an eigenbasis of t H and the Hamiltonian can be written as A possible definition for the density matrix of the system at inverse temperature β is In particular at zero temperature this yields where r 0 denotes the ground state of H. Assuming a decomposition E = E ⊗ E one can then trace over to define ρ . Let { f j } be a basis of E and { f * j } the dual basis, the trace over is defined as Note that tracing over is independent of the basis { f j } chosen, and does not require any inner product. With ρ at hand one then defines the Von Neumann and Rényi entropies in the usual way.
The main advantage of this construction is that the corresponding (Rényi) entropy Tr(ρ N ) maps within the path-integral approach to an Euclidean partition function on an N -sheeted Riemann surface. Underlying this result is the identification of the reduced matrix ρ with the partition function on a surface leaving open a slit along . Note that such a partition function can be computed purely in terms of matrix elements of the transfer matrix, and therefore it does not involve any inner product structure.
The disadvantages of this construction are twofold. The main one is that the reduced density matrix (and hence the entanglement entropy) may not be positive. While this may seem like a pathological property, loss of positivity in a non-unitary system might be acceptable depending on the context and motivations. The other one is that this definition only applies to eigenstates of H (and statistical superposition thereof). This stems for the fact that there is no canonical (i.e. basis independent) isometry between E and E * . In practice this means that knowing the ground state r 0 is not enough to compute the entanglement entropy, one also needs to know H to characterize w 0 .
Consider now an inner-product structure on E, i.e. a non-degenerate hermitian form 1 〈·, ·〉. By virtue of being non-degenerate, this inner product induces a canonical isometry between linear forms and vectors. For every vector v ∈ E, denote by v † the linear form defined by v † (x) = 〈v, x〉, x ∈ E .
Every element in E * can be written in this form, and the map I : v → v † is an antilinear isometry from E to E * . In particular one can associate a vector l j to every linear form w j such that w j = l † j . The vectors l j are what is commonly referred to as left eigenvectors of H. These are nothing but the eigenvectors of H † , the hermitian adjoint of H, which is characterized by the following property The transpose t H and the Hermitian adjoint H † are closely related: In particular the relation t H(w j ) = E j w j becomes While the previous prescription for the density matrix amounts in this context to an alternative prescription isρ For many non-unitary models there exist a natural notion of inner product that makes the Hamiltonian self-adjoint and is compatible with locality, 2 at the cost of not being definitepositive. For such an inner-product left and right eigenvectors coincide and it follows that This is typically the case within the CFT framework : the standard CFT inner product is such that L † n = L −n , and in particular the Hamiltonian is self-adjoint: L 0 = L † 0 . This is also the case for the loop model based on the Temperley-Lieb algebra [41] or for the Yang-Lee spin chain (see appendix E).
Let us now assume that the inner product is definite-positive. For a unitary system H = H † : left and right eigenvectors coincide, and both prescriptions yield the usual notion of density matrix and entanglement. For a non-Hermitian Hamiltonian operator H in general l j and r j are different (even if E j is real). If H is symmetric but not real (in some orthonormal basis), then the eigenvectors have non-real components, and are related through complex conjugation: On a more fundamental level this illustrates the fact that the canonical map between linear forms and vectors is antilinear.
The prescriptionρ, together with a positive-definite inner product, seems to be physically more natural than ρ as it yields a positive entanglement entropy (as can be seen from the Schmidt decomposition). Moreover it does not depend on H, only on the state considered and on the inner product. However this quantity is very much sensitive to the inner product chosen, and for a non-Hermitian Hamiltonian there is no canonical choice of a positive-definite inner product.
On a more technical side, when computing any quantity involvingρ in the path-integral formalism one needs to implement explicitly the (inner-product dependent) time-reversal operation r 0 → l 0 (e.g. r 0 → r * 0 for symmetric H) in order to get a consistent Euclidean description. Such a time-reversal defect can be thought of as a specific boundary condition in the tensor product C F T ⊗2 , which is typically a difficult problem.
It has been argued in [33] that for P T -symmetric Hamiltonian, the left and right ground states r 0 and l 0 coincide (while working with a definite positive inner product). This would circumvent this difficulty. However we found that this is not the case for the Yang-Lee model in finite size. Moreover assuming r 0 = l 0 immediately yields positive entanglement entropies, which again we found is not the case (both within our numerical and analytical calculations, see Figure 3).
In the following, when the model considered is non-unitary, we will choose (13) as the density matrix so that the Euclidean path-integral formalism described in Sec. 2.1, i.e. the interpretation of Rényi entropies as partition functions on a replicated surface with branch points, can be used straightforwardly. This is also the choice made in [39,41].
Within the Euclidean path-integral formalism an additional fact to take into account when studying non-unitary models is the existence of a primary state φ with a conformal dimension lower than the CFT vacuum h = 0 (i.e. the conformally invariant state). As was first pointed out in [33], this has a dramatic effect on the twist field : the most relevant twist operator is no longer τ 1 , but rather τ φ . Repeating the steps of section 2.2, the one-interval Rényi entropy in the ground state |φ〉 is mapped within the orbifold approach to where Φ = φ ⊗N . In [33] it was further claimed that the entanglement entropy in a non-unitary model behaves as where c eff = c − 24h φ is the effective central charge. However this result was based on an incorrect mapping to an Euclidean partition function, namely instead of (18). We claim that the behavior (19) is incorrect, and the Cardy-Calabrese formulas (4-5) for the entanglement entropy 3 cannot be applied, even with the substitution c → c eff .

The cyclic orbifold
The expression (4) suggests that the partition function on Σ N can be considered as the twopoint function of a "twist operator" of dimension Indeed, this point of view corresponds to the construction of the cyclic orbifold. Mathematically, one starts with N copies of the original CFT model (called the mother CFT, with central charge c, living on the original surface Σ) then mod out the N symmetry 4 (the cyclic permutations of the copies). This cyclic orbifold theory was studied extensively in [26][27][28]. We give an overview of the relevant concepts that we shall use.

The orbifold Virasoro algebra
All the copies of the mother CFT have their own energy-momentum tensor T j (z) [andT j (z) for the anti-holomorphic part]. Their discrete Fourier transforms (in replica space) are called T (r) (z), r ∈ {0, · · · , N − 1} and are defined by: The currents T j are all energy-momentum tensors of a conformal field theory, so their Operator Product Expansion (OPE) with themselves is: For two distinct copies, T j (z 1 ) T k (z 2 ) is regular ; on the unfolded surface, even when z 1 → z 2 the two currents are at different points. With that in mind, the OPE between the Fourier transforms of these currents can be written: where the indices r and s are considered modulo N . The modes of the currents are defined as: In the untwisted sector of the theory the mode indices m have to be integers since the operators T (r) (z) are single valued when winding around the origin. In the twisted sector however the operators T (r) (z) are no longer single-valued, and the mode indices m can be fractional. Generically in the cyclic N orbifold we have 3 When discussing the result of [33] the distinction between ρ = r 0 l † 0 and ρ = r 0 r † 0 is irrelevant since they argue that l 0 = r 0 at criticality. 4 For simplicity, in the following, we consider only the case when N is a prime integer.
The actual values of m appearing in the mode decomposition are detailed below: see (35) and (38). From the OPE (24) one obtains the commutation relations: where (m, n) ∈ ( /N ) 2 . The actual energy-momentum tensor in the orbifold theory is T orb (z) = T (0) (z). It generates transformations affecting all the sheets in the same way, so in the orbifold it has the usual interpretation (derivative of the action with respect to the metric). Correspondingly, the integer modes L (0) m∈ form a Virasoro subalgebra. The T (r) (z) for r = 0 also have conformal dimension 2, and play the role of additional currents of an extended CFT with internal N symmetry.

Operator content of the N orbifold
The untwisted sector Let z be a regular point of the surface Σ N . A generic primary operator at such a regular point, which we shall call an untwisted primary operator, is simply given by the tensor product of N primary operators φ 1 , . . . , φ N of dimensions h 1 , . . . , h N in the mother CFT, each sitting on a different copy of the model: In the case when Φ includes at least one pair of distinct operators φ i = φ j , it is also convenient to define the discrete Fourier modes where the indices are understood modulo N . The normalisation of Φ (r) (z) is chosen to ensure a correct normalisation of the two-point function: where h Φ = N j=1 h j . In particular, for a primary operator φ in the mother CFT, if one sets φ 1 = φ h and φ 2 = · · · = φ N = 1, one obtains the principal primary fields of dimension h: The OPE of the currents with generic primary operators are: where we have introduced the notations This expression reduces to a simple form in the case of a principal primary operator: From the expression (32), the product of T (r) (z) with an untwisted primary operator is singlevalued, and hence, only integer modes appear in the OPE: The twisted sectors The conical singularities of the surface Σ N are represented by twist operators in the orbifold theory. A twist operator of charge k = 0 is generically denoted as τ [k] , and corresponds to the end-point of a branch cut connecting the copies j and j + k. If A j denotes the j-th copy of a given operator A of dimension h A , one has: This relation can be considered as a characterisation of an operator τ [k] of the k-twisted sector. As a consequence, the Fourier components T (r) have a simple monodromy around τ [k] : and similarly for the primary operators Φ (r) and φ (r) h . Hence, the OPE of T (r) (z) with a twist operator can only include the modes consistent with this monodromy: If one supposes that there exists a "vacuum" operator τ [k] 1 in the k-twisted sector, one can construct the other primary operators in this sector through the OPE: For convenience, in the following, we shall use the short-hand notations: In a sector of given twist, most fractional descendant act trivially: Hence the short-hand notation:

Induction procedure
Suppose one quantises the theory around a branch point of charge k = 0 at z = 0. After applying the conformal map z → w = z 1/N from Σ N to a surface where w = 0 is a regular point, the currents T (r) (z) transform as: Accordingly, one gets for the generators: where the L n 's are the ordinary Virasoro generators of the mother CFT. It is straightforward to check that these operators indeed obey the commutation relations (27). Similarly, the twisted operator τ φ (39) maps to the primary operator φ h of the mother CFT: The relations (42)(43) are called the "induction procedure" in [28]. Using (42)(43) for This formula first appeared in [42,43], and, in the context of entanglement entropy, in [44].

Null-vector equations for untwisted and twisted operators
We intend to fully use the algebraic structure of the orbifold. If the mother theory is rational (i.e. it has a finite number of primary operators), then so is the orbifold theory. Also, from the induction procedure, we shall find null states for the twisted operators in the orbifold. In our approach, these null states are important, as they are the starting point of a conformal bootstrap approach.
Non-twisted operators. In the non-twisted sector, the null states are easy to compute. A state Φ in the non-twisted sector is a product of states of the mother theory. If one of these states (say, on the j th copy), has a null vector descendant in the mother theory, then the modes of T j (z) generate a null descendant. After an inverse discrete Fourier transform, these modes are easily expressed in terms of the orbifold Virasoro generators L (r) n . For instance, take the mother CFT of central charge c, and consider the degenerate operator φ 12 . We can parameterise the central charge and degenerate conformal dimension as and the null vector condition then reads: For a generic number of copies N we have and hence for N = 2, and Φ = φ 12 ⊗ φ 12 , we obtain More generally, any product of degenerate operators from the mother CFT is itself degenerate under the orbifold Virasoro algebra.
Twisted operators. The twisted sectors also contain degenerate states, which are of great interest for the following. For example, let us take k = 1, and characterise the degenerate states at level 1/N . A primary state τ obeys L (r) m τ = 0 for any r and positive m ∈ + r/N . For τ to be degenerate at level 1/N , one needs to impose the additional constraint: Using the commutation relations (27), we get: Thus, the operator τ is degenerate at level 1/N if and only if it has conformal dimension h τ = c 24 N − 1 N . This is nothing but the conformal dimension (21) of the vacuum twist operator τ 1 . Hence, one always has A more generic method consists in using the induction procedure. First, (51) can be recovered by applying (42)(43): One can obtain the other twisted null-vector equations by the same induction principle. Let us give one more example in the k = 1 sector: the case of τ φ 12 . The relations (42)(43) give: and hence τ φ 12 is degenerate at level 2/N . Note the insertion of some factors N in the nullvector of τ φ 12 as compared to the null-vector equation (46) of the mother CFT.

Yang-Lee two-interval correlation function
We consider the CFT of the Yang-Lee (YL) singularity of central charge c = −22/5, where the primary operators are the identity 1 = φ 11 = φ 14 with conformal dimension h 1 = 0, and φ = φ 12 = φ 13 with h φ = −1/5. We shall compute the following four-point function in the 2 orbifold of the YL model: In the N = 2 cyclic orbifold of the YL model, the untwisted primary operators are 1, Φ = φ ⊗φ and the principal primary fields φ (r) with r = 0, 1. They have conformal dimensions, respectively, h 1 = 0, h Φ = −2/5 and h φ = −1/5. Note that for N = 2 the only twisted sector has = 1, and hence τ ≡ τ. For the same reason, we shall sometimes omit the superscripts on the generators L (r) n , as r = 0, 1 are the only possible values. The twisted primary operators are τ 1 and τ φ , with conformal dimensions h 1 = −11/40 and h φ = −3/8. Here we have used the standard convention 〈ψ(∞) . . .〉 := lim R→∞ R 4h ψ 〈ψ(R) . . .〉 . Geometrically, this correlation function correspond to the partition function of the Yang-Lee model on a twice branched sphere, which can be mapped to the torus.
The identity operator of the YL model satisfies two null-vector equations: Through the induction procedure, this yields null-vector equations for the twist operator τ 1 : where the Fourier modes are r = 1 and r = 0 respectively, as required from (38). The first equation of (56) is generic for all N = 2 orbifolds, and determines the conformal dimension of the τ 1 operator. In contrast, the second equation is specific to the YL model. It only involves the integer modes, which all have the usual differential action when inserted into a correlation function. Hence, due to the second equation of (56), the derivation of G(x,x) is very similar to the standard case of a four-point function involving the degenerate operator φ 12 (see appendix B). The conformal block in z → 0 have the expression: with I 1 (x) = 2 F 1 (7/10, 11/10; 7/5|x) , and the total correlation function can be written:

OPE coefficients
The coefficients X j and Y j give access to the OPE coefficients in the 2 orbifold of the YL model: Recalling Mapping to the torus The mapping from the torus (with coordinates t) to the branched sphere (with coordinates z) is: This maps 0, x, 1, ∞ ← 1 2 , 1 2 (1 + τ), τ 2 , 0. The relation between x and the nome q = e 2iπτ is given by: Mapping the torus to the branched sphere, the partition function transforms as: The torus partition function of the Yang-Lee model involves two characters, χ 1,1 (τ) (1), and χ 1,2 (τ) (φ).
The characters of minimal models have the well-known expression χ r,s (τ) = K r,s (τ)− K r,−s (τ), with: where η is the Dedekind eta function. We expect the following relation between the conformal blocks of the correlation function and the characters of the theory: Those two relations are not trivial, and the simplest way to prove them seems to be by showing that the right-hand terms are vector modular form, with the same modular transformations as χ (see by example [45] for details). An easy check consist in expanding the right-hand side in power of q, confirming the equality for first orders: The relation 62 between the partition on the torus and G(x,x) is verified:

Yang-Lee one-interval correlation function
With minimal modifications to the previous argument, we can also compute the following four-point function : Which, physically, is related to the generalised Rényi entropy S N =2 (x, φ, τ 1 ).
Technically, it is more convenient to work with twist operators located at 0 and ∞, so we introduce Using the suitable projective mapping, one has the relation: Then, through the null-vector of τ 1 , we can obtain a differential equation of order two for this correlation function: The Riemann scheme of this equation is: which is consistent with the OPEs: By appropriately shifting the function F , F (x,x) = |x| 4/5 |1 − x| 8/5 f (x,x), we can turn 68 in 174. At that point we can simply re-use the results of the last sections with parameters: The final result for the four-point function G(x,x) (66) is G(x,x) = |x| 4/5 |1 − x| 11/5 × X 1 2 F 1 (4/5, 7/10; 11/10|x) 2 +X 2 x −1/10 2 F 1 (3/5, 7/10; 9/10|x) 2 , where X 1 , X 2 are given in (188), and the parameters a, b, c, d are given in (73). Using the identity (166) on hypergeometric functions, we see that the solution (74) for G(x,x) agrees with the direct computation given in Appendix D.

Orbifold Ward identities
Generically the null vectors for a twist operator can involve some generators L (r) m , with r = 1 and fractional indices m ∈ + r/N [see (38)], which do not have a differential action on the correlation function. In this situation, we shall use the extended Ward identities to turn the null-vector conditions into a differential equation for the correlation function.
Let us consider the correlation function: where ( 2 , 3 ) are any two operators and (| 1 〉, | 4 〉) are any two states of the cyclic orbifold. Each operator j or state | j 〉 can be in a twisted sector [k j ] with k j = 0 mod N , or in the untwisted sector (k j ≡ 0 mod N ). The overall N symmetry imposes a neutrality condition: k 1 + k 2 + k 3 + k 4 ≡ 0 mod N . Let C be a contour enclosing the points {0, x, 1}. Then this is a closed contour for the following integral: where m j ∈ + r k j /N for j = 2, 3, 4. Then, by deforming the integration contour to infinity, we obtain the following identity: where and the coefficients a p , b p , c p , d p are defined by the Taylor expansions: In (77) we have used the notation: where C x is a contour enclosing the point x. If all the j 's are chosen among the primary operators or their descendants under the orbifold Virasoro algebra, then the sums in (77)

Ising two-interval ground state entropy
The Ising model is the smallest non-trivial unitary model, with three fields {1, σ, ε}, of conformal charges 0, 1 /16, 1 /2. Its central charge is 1 2 . By using the correspondance between the Ising model and free fermions, the two-interval case was computed in [16], for any value of N . Here, as an example, we compute the N = 2 two-interval entropy, with our method, and using only the null vector of the model. The N = 2 orbifold of the Ising model will have central charge 1. The first field in the twisted sector is τ 1 , the identity twist, of charge h τ 1 = 1 /32.
The two null vector of the identity in the mother theory are: And, through the induction process, we obtain null vector equations for τ 1 .

Yang-Lee one-interval ground state entropy
Definitions. As argued above, the N = 2 ground state entropy of the YL model is related to the correlation function: where Φ = φ ⊗ φ. It will be convenient to work rather with the related function: Null vectors and independent descendants. In the mother theory, the primary field φ has two null-descendants: Through the induction procedure, this implies the following null vectors for the twist field τ φ : The Ward identities (77) will give relations between descendants at different levels. In order to get an idea of the descendants we need to compute, we need to know the number of independent descendants at each level. The number of Virasoro descendants at a given level k is equal to the number of integer partitions of k. In a minimal model not all those descendants are independent. The number of (linearly) independent descendants at level k of a primary field φ is given by the coefficients of the series expansion of the character associated with φ.
Explicitly, if τ is the modular parameter on the torus, the character of a field φ rs in the minimal model p,p is given by: The coefficient of order k in the series expansion with parameter q of the character gives the number of independent fields at level k. Through the induction procedure (42)(43), the module of τ φ rs under the orbifold algebra has the same structure as the module of φ rs under the Virasoro algebra. Hence this coefficient also gives the number of descendants at level k /N in the module of τ φ rs . For example, for N = 2 and (p, p ) = (5, 2) (Yang-Lee), the numbers of independent descendants for the field τ φ = τ φ 12 are given in the following So up to level 3, all descendants at an integer level (even formed of non-integer descendants) can be re-expressed in terms of integer modes. One gets explicitly: Ward identity. We now use (77) with the indices (m 1 , . . . , m 4 ) = (1/2, 0, 0, −5/2) and: Recalling that L (1) 0 Φ = 0, the only surviving terms are: Inserting the explicit expressions for the coefficients d p , we get: Differential equation. Combining (90) and (91), one gets a linear relation involving only the integer modes L (0) m . Using (169), this leads to the third-order differential equation: The Riemann scheme of this equation is: m k with r 1 + · · · + r k ≡ 0 mod N . With this choice, the toroidal partition function of the orbifold has a diagonal form Z = j |χ j | 2 in terms of the characters (see [27,28]). The OPEs under the invariant subalgebra 0 are then: Note that L −1/2 τ φ is a primary operator with respect to 0 . The local exponents (93) are consistent with these OPEs.

Holomorphic solutions.
To express the solutions in terms of power series around x = 0, it is convenient to rewrite (92) using the differential operator θ = x∂ x . On has, for any k ∈ : This yields the new form for (92): where: The key identity satisfied by the operator θ is, for any polynomial P and any real α: Hence, if we choose α to be a root of P 0 (i.e. one of the local exponents at x = 0), there exists exactly one solution of the form: a n x n , with a n = 1 .
For example, the conformal block corresponding to α = 1/2: is given by The series converges for |x| < 1, and may be evaluated numerically to arbitrary precision using (101).

Numerical solution of the monodromy problem.
To determine the physical correlation function F (x,x) (87), we need to solve the monodromy problem, i.e. find the coefficients for the decompositions: where (I 1 , I 2 , I 3 ) are the holomorphic solutions of (92) with exponents (1/2, 2/5, 9/10) around x = 0, and (J 1 , J 2 , J 3 ) are the holomorphic solutions with exponents (4/5, 2/5, 3/5) around x = 1. In the present case, we do not know the analytic form of the 3 × 3 matrix A for the change of basis: However, this matrix can be evaluated numerically with arbitrary precision, by replacing the The coefficients giving a monodromy-invariant F (x,x) consistent with the change of basis (105) are: X 1 = 30.6594 , We have set Y 1 to one because it corresponds to the conformal block with the identity as the internal operator.

OPE coefficients.
The coefficients X i and Y j are related to the OPE coefficients in the N = 2 orbifold of the YL model as follows: Our numerical procedure can be checked by comparing some of the numerical values (107) to a direct calculation of the OPE coefficients done in Appendix C. They match up to machine precision.

Excited state entropy for minimal models at N = 2
The same type of strategy can be used to compute excited state entropies in other minimal models, for some specific degenerate states, and small values of N . In this section we solve explicitly the case of N = 2 for an operator degenerate at level 2.

Definitions.
We consider the minimal model , where p, p are coprime integers. Using the Coulomb-gas notation, φ 21 is one of the states of the mother CFT which possess a null vector at level 2. In terms of the Coulomb-gas parameter g = p/p , the central charge and conformal dimensions of the Kac table read: In particular, we have h 21 = (3g − 2)/4. The null-vector condition for φ 21 reads: As shown before, in a unitary model (when q = p−1), the entanglement entropy of an interval in the state |φ 21 〉 is expressed in terms of the correlation function Null vectors. The null vectors of τ 1 and Φ can be obtained through the induction procedure: Ward identity. Using the Ward identity (77) with the choice of indices (m 1 , · · · , m 4 ) = (0, −1 /2, −1 /2, −1) for the function −1 |Φ〉, and the null-vector condition (118), we obtain the relation: Using the orbifold Virasoro algebra (27) and the explicit expression for the d p 's, we get: Differential equation. By inserting the first null-vector condition of (119) into G(x,x), we get: Then, the substitution of the ( L (1) −1 ) 2 term by (121) gives a linear relation involving only the modes L (0) m , which have a differential action on G(x,x). This yields the differential equation for G(x,x): Its Riemann scheme is: These exponents correspond to the fusion rules in the channels x → 1 and x → 0, respectively.

Solution.
The problem can be solved as in Sec. B.1. If we multiply the correlation function by the appropriate factor, equation (123) becomes the hypergeometric differential equation: Following exactly the reasoning in Sec. B.2, we find the correlation function: where d = c − a − b and Of course, we may compare this solution to the one obtained by a conformal mapping of the four-point function 〈φ 21 φ 21 φ 21 φ 21 〉. The latter also satisfies a hypergeometric equation, and using the transformation 166, the two solutions can be shown to match.

Excited state entropy for minimal models at N = 3
We now consider the correlation function in the N = 3 orbifold of the minimal model (p, p ): The conformal mapping method would result in a much more complicated 6-point function, which is not the solution of an ordinary differential equation. But through the orbifold Virasoro structure, we can obtain such an equation. The method is similar in spirit to what was done for N = 2, finding the null vector conditions on the field Φ, then using the Ward identities.
Null vectors. We will need the null vectors of Φ 3 up to level 4: A similar relation can be found for the correlation functions of the form by using other Ward identities.

Differential equation.
Putting everything together we find the following differential equation: where θ = x∂ x and: The Riemann scheme is given by: Interpretation in terms of OPEs. The conformal dimensions of the internal field in the channels x → 1 and x → 0 are respectively: Since 〈τ 1 . τ 1 .(φ 31 ⊗ 1⊗1)〉 ∝ 〈φ 31 〉 = 0, the conformal block with internal dimension h 31 is in fact not present in the physical correlation function. In the channel x → 0, this is mirrored by the presence of two fields separated by an integer value : there is a degeneracy for the field τ φ 41 .
This can be verified by applying the bootstrap method. The equation has 4 linearly independent solutions, which can be computed by series expansion around the three singularities. Like in 5.3, the monodromy problem can be solved by comparing the expansions in their domain of convergence. Up to machine precision the coefficient corresponding to h φ 3,1 vanishes for all central charges. Nevertheless, the other structure constants converge, and we can still compute the full correlation function through bootstrap. For the non-zero structure constants, we checked that they were matching their theoretical expressions for simple minimal models (Yang-Lee and Ising).
The effective presence of only three conformal blocks also seem to imply that we should have been able to find a degree three differential equation, instead of four, for this correlation function. However we have not managed to derive such a differential equation.
In this section we describe a lattice implementation of the twist fields in the lattice discretisation of the minimal models, namely the critical Restricted Solid-On-Solid (RSOS) models. Entanglement entropy in RSOS models has already been considered in [47] for unitary models and in [39] for non-unitary models, but for a semi-infinite interval and away from criticality.

The critical RSOS model
Let us define the critical RSOS model with parameters (m, k), where m and k are coprime integers and k < m. Each site r of the square lattice carries a height variable a r ∈ {1, 2, . . . , m}, and two variables a and b sitting on neighbouring sites should differ by one : |a − b| = 1. The Boltzmann weight of a height configuration is given by the product of face weights: where the crossing parameter λ is The quantum model associated to the critical RSOS model is obtained by taking the very anisotropic limit u → 0 of the transfer matrix. For periodic boundary conditions, one obtains a spin chain with basis states |a 1 , a 2 , . . . , a L 〉, where a i ∈ {1, . . . , m}, and |a i − a i+1 | = 1, and the Hamiltonian is where e i only acts non-trivially on the heights a i−1 , a i , a i+1 : and the indices i ± 1 are considered modulo L. For simplicity, we now consider the RSOS model on a planar domain. The lattice partition function Z RSOS and the correlation functions admit a graphical expansion [48] in terms of nonintersecting, space-filling, closed loops on the dual lattice. The expansion of Z RSOS is obtained by associating a loop plaquette to each term in the face weight (137) as follows: Then, after summing on the height variables, each closed loop gets a weight β = 2 cos λ. Furthermore, following [49], correlation functions of the local variables also fit well in this graphical expansion. Let us recall, for example, the expansion of the onepoint function 〈ϕ q (a r )〉. For any loop which does not enclose r, the height-dependent factors from (137) end up to sin λb/ sin λa, where a (resp. b) is the outer (resp. inner) height adjacent to the loop. Thus, summing on the inner height b gives the loop weight: where we have introduced the adjacency matrix A ab = 1 if |a − b| = 1, and A ab = 0 otherwise. For the loop enclosing r and adjacent to it, the factor ϕ q (b) should be inserted into the above sum, which gives: where β q = 2 cos πq Repeating this argument recursively, in the graphical expansion of 〈ϕ q (a r )〉, one gets a loop weight β q for each loop enclosing the point r. The N -point functions of the ϕ q 's are treated similarly, through the use of a lattice Operator Product Expansion (OPE) [49]. This critical RSOS model provides a discretisation of the minimal model (p, p ), with central charge and conformal dimensions: with the identification of parameters: The operator ϕ q changes the loop weight to β q : thus, in this sector in the scaling limit, the dominant primary operator, which we denote φ q , has conformal dimension It is then easy to show 5 that h φ q is one of the dimensions of the Kac table (147). Note that h φ q=k = 0 corresponds to the identity operator.

Partition function in the presence of branch points
We consider the reduced density matrix ρ for the interval = [i, j]. A generic matrix element 〈a i , a i+1 , . . . , a j |ρ |a i , a i+1 , . . . , a j 〉 corresponds to the partition function of the lattice shown in Fig. 1, with the heights a i , . . . a j and a i , . . . a j fixed, and the other heights summed over. In this convention, a branch point (or twist operator) sits on a site r of the square lattice, and is denoted t( r). Computing the n-th Rényi entropy (2) amounts to determining the partition function Z (n) RSOS on the surface obtained by "sewing" cyclically n copies of the diagram in Fig. 1 along the cut going from a i to a j . The graphical expansion of the partition function on this surface with two branch points is very similar to the case of a planar domain. The only difference concerns the loops which surround one branch point. Since such a loop has a total winding ±2πn instead of ±2π, the height-dependent factors from (137) end up to (sin λb/ sin λa) n , where a (resp. b) is the external (resp. internal) height adjacent to the loop. Since (sin λb) n is not an eigenvector of the adjacency matrix A, the sum over b does not give a well-defined loop weight. For this reason, we introduce a family of modified lattice twist operators: With this insertion of ϕ q at the position of the twist, the sum over the internal height gives: and hence any loop surrounding t q gets a weight β q (145). Thus, the scaling limit of t q is the primary operator τ φ q , which belongs to the twisted sector, and has conformal dimension In particular, since β k = β = 2 cos λ, one has φ k = 1, and the lattice operator t k corresponds to τ 1 in the scaling limit. Note that the "bare" twist operator t is itself a linar combination of the t q 's: The scaling limit of t is thus always determined by the term t 1 , since it has the lowest conformal dimension. In the case of a unitary minimal model (k = 1), this corresponds to τ 1 . In contrast, for a non-unitary minimal model (p, p ), t scales to the twist operator τ φ 1 , where φ 1 is the primary operator with the lowest (negative) conformal dimension in the Kac table:

Rényi entropies of the RSOS model
When defining a zero-temperature Rényi entropy, two distinct parameters must be specified: 1. The state |ψ〉 R in which the entropy is measured (or equivalently the density matrix ρ = |ψ〉 R 〈ψ| L ).
2. The local state φ q of the system in the vicinity of branch points. This determines which twist operator t q should be inserted. In the case of the physical Rényi entropy defined as (2), one inserts the linear combination t( r) = m q=1 x q t q ( r). In the following, we will be interested in the Rényi entropy of an interval of length in the spin chain H RSOS (139) of length L with periodic boundary conditions. This corresponds to the lattice average value: where |Ψ〉 = |ψ〉 ⊗N . The "physical" Rényi entropy (2) is related to: The average values (154-155) scale to correlation functions on the cylinder {z| 0 ≤ Im z < L}: where u = i , and similarly for t q replaced by t. Using the conformal map x = exp(2πu/L), these are related to the correlation functions on the complex plane: In the case of the (generalised) Rényi entropy in the vacuum, we have ψ = 1, and this becomes a two-point function, which is easily evaluated: In particular, when φ q = 1, one recovers the result from [14]: For a generic state |ψ〉 however, the entropy S N (x, ψ, τ φ q ) remains a non-trivial function of , and does not reduce to the simple form (158).

Numerical setup
We have computed some Rényi entropies (154)  A lattice eigenvector (scaling either to |1〉, or to |φ〉) of H RSOS with periodic boundary conditions is obtained by exact diagonalisation with the QR or Arnoldi method, and then used to construct the reduced density matrix ρ , where = [0, ]. For the computation of S N (x, 1, t q ) and S N (x, ψ, t q ), the factor ϕ q (a 0 ) ϕ q (a ) [see (150)] is inserted into the trace (2) of ρ N . For the computation of S N (x, φ, t), no additional factor is inserted. From the above discussion, we expect S N (x, φ, t) to be described by the insertion of the dominant twist operators τ φ .

Results for entropies at N = 2 in the Yang-Lee model
Here we present our numerical results obtained with the procedure described above. In all the cases considered the cylinder correlation functions have been rescaled with the factor L 4h , where h is the appropriate twist field conformal dimension. The collapse of various finite size data further confirms the correct identification of the twist field (τ φ or τ 1 ).
The results obtained are in excellent agreement with our CFT interpretation (18) and with our analytical results. Moreover they are clearly not compatible (see Fig. 3) with the claim which can be found in the literature [33,41] (the effective central charge of the Yang-Lee model is c eff = 2/5).

Conclusion
In this article we have studied the Rényi entropies of one-dimensional critical systems, using the mapping of the N th Rényi entropy to a correlation function involving twist fields in a N cyclic orbifold. When the CFT describing the universality class of the critical system is rational, so is the corresponding cyclic orbifold. It follows that the twist fields are degenerate : they have null vectors. From these null vectors a Fuchsian differential equation is derived, although this step can be rather involved since the null-vector conditions generically involve fractional modes of the orbifold algebra. The last step is to solve this differential equation and build a monodromy invariant correlation function, which is done using standard bootstrap methods. We have exemplified this method with the calculation of various one-interval Rényi entropies in the Yang-Lee model, a two-interval entropy in the Ising model and some one-interval entropies computed in specific excited states for all minimal models.
We have also described a lattice implementation of the twist fields in the lattice discretisation of the minimal models, namely the critical Restricted Solid-On-Solid (RSOS) models. This allows us to check numerically our analytical results obtained in the Yang-Lee model. Excellent agreement is found.
The main limitation of our method is that its gets more involved as N increases, and as the minimal model (p, p ) under consideration gets more complicated (i.e. as p and/or p increases). For this reason, we have limited our study to N = 2 and N = 3 in the Yang-Lee model (5,2). However, this method is applicable in a variety of situations where no other method is available, for instance when the subsystem is not connected (e.g. two-intervals EE). We will address this case in a following paper.
Another interesting research direction would be to develop a Coulomb Gas formalism for the cyclic orbifold, as it would provide an efficient tool to solve the twist-field differential equations à la Dotsenko-Fateev. Indeed, the Coulomb Gas yields a very natural way to write down conformal blocks (in the form of closed contour integrals of screening operators), to compute their monodromies, and from there to solve the bootstrap. is given by: • Under the transformation x → 4 x/(1 + x) 2 , we have: (166)

B Four-point function satisfying a second-order differential equation
In this appendix we compute the correlation function in the 2 orbifold of the YL model. It follows from the null-vector This is the standard form of a null-vector at level 2, which yields in the usual way a second order differential equation whose solutions are hypergeometric functions. For completeness we recall the key steps in computing G(x,x).

B.1 Differential equation
A standard CFT argument yields, for any n ∈ , any primary operators ( 2 , 3 ), and any states (| 1 〉, | 4 〉): Then (56) translates into the ordinary differential equation for G ( x,x): This equation has the following Riemann scheme, giving the local exponents, i.e. the allowed power-law behaviours at the three singular points 0, 1 and ∞: In the limits x → 0, 1, ∞, we have the OPE: So the local exponents (171) are consistent with the internal states {1, Φ} of the conformal blocks for the channels. If we perform the appropriate change of function to shift two of these local exponents to zero: then (170) turns into the hypergeometric differential equation: with parameters: It is also convenient to introduce the parameter d = c − a − b. If one repeats the argument with the anti-holomorphic generatorsL n , one obtains the same equation as (174), with (x, ∂ x ) replaced by (x,∂ x ).

B.2 Determination of the four-point function
A basis of holomorphic solutions to (174) is given by: where 2 F 1 is Gauss's hypergeometric function (161). The basis I j has a diagonal monodromy around x = 0: Similarly, by the change of variable x → 1 − x, one obtains a basis of solutions with a diagonal monodromy around x = 1: We shall construct a solution of the form From the properties of the operators τ 1 and Φ (see Sec. 2), G(x,x) should be single-valued, which imposes the form X i j = δ i j X i for the coefficients in (180). The solution should also admit a decomposition of the form: G(x,x) = |x| 11/10 |1 − x| 11/10 Again, single-valuedness of G(x,x) imposes the form Y k = δ k Y k . The key ingredient to determine the coefficients X j and Y j is the matrix for the change of basis between {I 1 (x), I 2 (x)} and {J 1 (x), J 2 (x)}. Using the properties (162-163) of hypergeometric functions, one obtains: where A is given in (165). Comparing (180) and (181), we get the matrix relations: Imposing a diagonal form for X and Y yields two linear equations on (X 1 , X 2 ): Since the entries of A are real, these two relations are equivalent. Similarly, one gets a linear relation between Y 1 and Y 2 . Finally, one gets the ratios: The symbol Γ denotes Euler's Gamma function, and we also introduced the short-hand notation: Moreover, the term |J 1 (x)| 2 in (181) corresponds to the OPE τ 1 × τ 1 → 1, which fixes Y 1 = 1.

C Direct computation of OPE coefficients of twist operators
In this appendix, we perform the computation of the structure constants appearing in the Yang-Lee model on the N = 2 orbifold. They provide a non-trivial check of the validity of our method based on solving the differential equation for conformal blocks. In the following, 〈. . .〉 denotes the average in the orbifold theory, whereas 〈. . .〉 Σ 2 (resp. 〈. . .〉 ) denotes the average in the mother theory on the two-sheeted Riemann surface (resp. on the Riemann sphere).
Some of those results were already obtained in [33], a generic way of computing those threepoint functions can be found in [50]. In the specific case of the Ising model similar three-point functions were found in [19] and [32].
For three-point functions involving only untwisted operators, the correlation function decouples betwenn the N copies. For instance: • C τ φ , Φ, L −1 /2 L −1 /2 τ φ : we also need this structure constant which involves a descendant state. The behaviour of the descendants states during the unfolding is given by the induction procedure 42: Hence, for the three-point function: This four-point function is computed in D. To compute the structure constant, we also need to normalize the descendant state:

E Quantum Ising chain in an imaginary magnetic field
We consider the Hamiltonian: with periodic boundary conditions, with h and λ real, in the regime 0 < λ < 1.
Within the usual inner product all operators σ a j are self-adjoint, and it is clear that H is not (the matrix representing H in the usual basis is symmetric but not real). Alternatively one can work with a different hermitian form, namely According to this hermitian form -which is not definite positive -the Hamiltonian density σ x j σ x j+1 + σ z j + ihσ x j (and therefore H itself) is self-adjoint. With the usual inner product, note that P H P = H † , so P maps right eigenvectors to left eigenvectors. In particular in the -unbroken phase, we have where |r 0 〉 is the ground state. Then |r 0 〉 = |l 0 〉 iff |r 0 〉 is an eigenstate of P. For small systems (L = 1, 2) one can check analytically that this is not the case. We have observed numerically that this trend persists for larger systems. A curious observation is that for a single site (L = 1), the Hamiltonian is not diagonalizable at the transition : the two lowest eigenvalues E 0 and E 1 merge into a non-trivial Jordan block. It would be interesting to study whether this is also the case for larger systems, as it would suggest some logarithmic behavior in the continuum.
Despite being non-hermitian, the eigenvalues of H are either real, or they appear in pairs of complex conjugates (E, E * ). This can be seen using symmetry [51], or simply by noting that after the unitary similarity transformation H = U H U † with U = L j=1 exp(iπσ z j /4), one gets a real, nonsymmetric operator H.
For h = 0, the Hamiltonian is Hermitian, and thus its spectrum is real. The regime 0 < λ < 1 and h = 0 corresponds to the (anisotropic limit of) the 2d Ising model in the high-temperature phase, where the correlation length ξ is finite. When h is increased while λ is kept constant, the ground state and first excited energies remains finite, up to a threshold value h c (λ, L), where they "merge" into a complex conjugate pair: see Fig. 6. The point h c corresponds to the vanishing of the partition function in the 2d Ising model. In the scaling limit, it converges to a critical point h c (λ, L) → h c (λ), called the Yang-Lee edge singularity. The finite-size study of the Yang-Lee edge singularity through the model (204) is rather subtle: for a given system size of L sites, one should first determine the threshold value h c (λ, L), and then approach this value from below. We have computed numerically the oneinterval ground state N = 2 Rényi entropy S 2 for the model (204) at λ = 0.8 and system sizes L = 12, 14, 16, 18 sites. The density matrix is defined as in the rest of the paper as ρ = r 0 l † 0 , so the quantity we compute corresponds at criticality to (18). These numerical calculations lead to the following observations, depicted in Fig. 7. In the off-critical regime h h c (λ, L), the entropy S 2 has a concave form. Then, when increasing the value of h and approaching h c (λ, L) from below, the function undergoes a crossover to the convex form predicted by CFT (92). While not being positive, the entanglement entropy defined using ρ = r 0 l † 0 is surprisingly effective at detecting the phase transition.