Asymptotic Correlations in Gapped and Critical Topological Phases of 1D Quantum Systems

Topological phases protected by symmetry can occur in gapped and---surprisingly---in critical systems. We consider the class of non-interacting fermions in one dimension with spinless time-reversal symmetry. It is known that the phases in this class are classified by a topological invariant $\omega$ and a central charge $c$. Here we investigate the correlations of string operators in order to gain insight into the interplay between topology and criticality. In the gapped phases, these non-local operators are the string order parameters that allow us to extract $\omega$. More remarkable is that the correlation lengths of these operators show universal features, depending only on $\omega$. In the critical phases, the scaling dimensions of these operators serve as an order parameter, encoding both $\omega$ and $c$. More generally, we derive the exact long-distance asymptotics of these correlation functions using the theory of Toeplitz determinants. We include physical discussion in light of the mathematical results. This includes an expansion of the lattice operators in terms of the operator content of the relevant conformal field theory. Moreover, we discuss the spin chains which are dual to these fermionic systems.

In a previous work, we extended the well-known classification of the gapped topological phases of quadratic fermionic Hamiltonians with spinless time-reversal symmetry [28] ('BDI class' of Altland and Zirnbauer's tenfold way [29]) to gapless topological phases [13]. These are labelled by a topological invariant ω (∈ Z) and the central charge c (∈ 1 2 Z ≥0 ) of the conformal field theory (CFT) that describes the continuum limit if the model is critical. If the system is gapped, we say that c = 0 and ω reduces to the well-known winding number of the BDI class [28]. What allowed for a complete analysis was the fact that each Hamiltonian in this class can be efficiently encoded into a holomorphic function f (z) on the punctured complex plane C \ {0}. Remarkably, c and ω can then be obtained by counting zeros of f (z) (see Figure 1). This rephrasing allowed us to argue that two critical models in this class can be smoothly connected if and only if they have the same topological invariants and central charges.
What remains an open question, however, is the extent to which the topological nature of these gapped and gapless phases is reflected in their correlation functions. Relatedly, it is natural to ask how the correlations are encoded in f (z)-especially since c and ω are easily derived from its zeros. Moreover, our earlier work left an uneasy tension: distinct critical phases could be distinguished by the topological invariant ω, yet it was not clear to what extent this lattice quantity is related to the CFT in the continuum. Hence, bridging this gap in terms of a lattice-continuum correspondence is desirable. More generally, since these models are exactly solvable we can hope to obtain a lot of information, and perhaps uncover unexpected features.
The aim of this work is twofold: on the one hand, we focus on answering the aforementioned questions conceptually, linking universal properties of correlations to the function f (z) and shedding light on the interplay of criticality and topology. On the other hand, since our models allow for a rigorous analysis, we give derivations of exact asymptotic expressions for important correlators. The method we use, Figure 1. The Hamiltonians we consider can be expanded in a basis H α (defined below equation (2)). The physics is encoded in the meromorphic function f (z). The given definitions of c and ω classify the phases of H BDI , where Z[g] denotes the (multi)set of zeros of g (with multiplicity) and N p the order of the pole at the origin. Physically, c encodes the low-energy behaviour, and ω the topological properties. Universal asymptotics of the ground-state correlation functions considered in this work. If c = 0 (i.e. the system is gapped), then O α has exponentially decaying correlations with correlation length ξ α . The ratios ξ α /ξ β are a universal function of ω, with the global scale set by 1/ξ := min ζ∈Z[f ] |log|ζ||; see the discussion before equation (9). (There is long-range order, i.e. a α = 0, if and only if α = ω; see Theorem 1b.) If c = 0 and the zeros on the unit circle have multiplicity one (i.e. the bulk is described by a CFT with central charge c), then the correlation functions obey a power law with universal scaling dimension ∆ α ; see Theorem 4. Note that there are exceptional cases that behave differently, as discussed in the text.
Toeplitz determinant theory, has a long association with statistical mechanics (for a review, see [30]), and our analysis generalises the pioneering work of [31][32][33] to a wider class of physical models. Since topological phases cannot be distinguished locally, in this work we study the correlations of so-called string-like objects O α (labelled by α ∈ Z), meaning that O α (1)O α (N ) involves an extensive (∼ N ) number of operators. Using Wick's theorem, these correlations reduce to N × N determinants. We calculate their asymptotic behaviour using the theory of Toeplitz determinants [30], phrasing the answers in terms of the zeros of f (z). Figure 2 summarises some of the main results. In the gapped case (c = 0), it is well-known that SPT phases can be distinguished by string order parameters [34][35][36][37][38], and we indeed prove that O α has long-range order if and only if α = ω. More surprising is that the ratios of the correlation lengths of these operators are universal, i.e. they depend on ω only. Moreover, the largest correlation length has a universal relationship to the zero of f (z) which is nearest to the unit circle. In the critical case (c = 0), all correlations are algebraically decaying and we obtain the corresponding scaling dimensions of O α . It turns out that measuring these gives access to both c and ω. Moreover, we propose a continuum-lattice correspondence for these operators. We expect that this correspondence will prove useful in exploring the effect of interactions on the phase diagram.
Since the physical consequences of our results can be understood without going into the mathematical details, we structure the paper as follows. First, in Section 2, we outline the model, state our main results, and discuss connections to previous works. Then in Section 3 we give further details of how our results fit into the broader physical context. In particular, we discuss general approaches to string order parameters and the consequences of universality in the gapped phases, give a CFT analysis of long-distance correlations and also show how our results allow us to deduce critical exponents. Only after this do we give the mathematical preliminaries in Sections 4 and 5. The proofs of our results then follow in Sections 6 and 7 for the gapped and critical cases respectively. Finally, in Section 8, we explain how our results may be extended in different directions.

Statement of main results
2.1. The model. Consider a periodic chain where each site has a single spinless fermionic degree of freedom 1 {c † n , c n ; n = 1 . . . L}. For convenience define the Majorana modes on each site: where {γ n , γ m } = 2δ nm and {γ n ,γ m } = 0. Our class of interest -time-reversal symmetric, translationinvariant free fermions with finite-range couplings -has Hamiltonian This can be understood as an expansion in the basis H α = i /2 n∈sitesγ n γ n+α . The coupling between sites has maximum range α l/r to the left and right. This model has an antiunitary symmetry, T , that acts as complex conjugation in the occupation number basis associated to the fermions c n and satisfies T 2 = 1. The Majorana operators γ n (γ n ) are called real (imaginary) since T γ n T = γ n and Tγ n T = −γ n . This class of models is also invariant under parity symmetry P = e iπ j c † j c j . We study the thermodynamic limit L → ∞, with α l/r finite but not fixed -i.e. we will consider models with differing maximum range. The results given in this section are all for such finite-range chains, but we discuss the extension to long-range chains in Section 8.1. This model was first analysed in its spin chain form (see Section 2.5) in reference [39].
The coupling constants t α establish a one-to-one correspondence between H BDI and the complex functions This is a holomorphic function away from a possible pole at the origin. By the fundamental theorem of algebra, f (z) is specified by the degree of this pole and a multiset of zeros (up to an overall multiplicative constant) . The basic relevance of f (z) is that |f (e ik )| gives the one-particle energy of a mode with momentum k. The phase arg(f (e ik )) is the angle required in the Boguliobov rotation that defines these quasiparticle modes [13,40]. Remarkably, many other physical questions can be answered through simple properties of this function. Note that we will consistently abuse notation f (k) := f (z = e ik ) whenever we restrict z to the unit circle. H BDI has two invariants that label both gapped and gapless phases (see also Figure 1): where N z is the number of zeros of f (z) inside the unit disk and N p is the degree of the pole at the origin. If c = 0, the model is gapped. For gapless models, c is the central charge of the low energy CFT when the zeros on the circle are non-degenerate 3 . Note that ω is an invariant since it cannot change under smooth motion of the zeros without changing c. It is moreover topological: it cannot be probed locally, but distinguishes phases and manifests itself through protected edge modes [13]. That the pair (c, ω) specifies the phases of H BDI was shown in reference [13]. If in addition to the symmetries P and T that stabilise the aforementioned phases, one also enforces translation symmetry, then there are additional invariant signs, denoted Σ, that are discussed in Appendix A.2. Note that the equivalence between H BDI and f (z) allows us to easily find a Hamiltonian within each phase: H ω is a representative of the gapped phase with winding number ω and H ω + H 2c+ω is a representative of the gapless phase (c, ω).
where x = |α|/2 for α even and x = (|α| − 1)/2 for α odd (the phase factors make O α hermitian). These operators are a cluster of |α| Majorana operators to the right of site n multiplied by an operator giving the parity of the number of fermions to the left of n. Such operators appear naturally as we discuss in Section 3.1, see also [41]. There are two typical behaviours for these correlators. Let angle brackets denote ground-state expectation value, then in the gapped case we expect: A α and δ α as well as the correlation length ξ α do not depend on N . If A α = 0 we have long-range order. B α is Θ(1) and may include an oscillation with N . Note that, by translation invariance, we could equally well have considered the correlation function O α (r)O α (N + r) for any r ∈ Z -we fix r = 1 throughout for notational convenience. Note that O α (1)O β (N + 1) = 0 if α = β as a simple consequence of the Majorana two-point functions given in Section 4.
We will see below that the ground state expectation value lim N →∞ O α (1)O α (N + 1) in the gapped phase with winding number ω is non-zero only when α = ω, and is hence an order parameter for that phase.
Because these correlators contain a string of fermionic operators of length of order N , these are called string order parameters with value A α . Note that in the case that O α is local, (as happens in the spin picture given in Section 2.5), it is usual to call the one point function O α (n) the order parameter. This is because in that case the ground state will spontaneous collapse such that In this work we prefer to use a single convention and always refer to the two point function as 'the' order parameter.
At critical points with a low-energy CFT description we expect: where ∆ α is the smallest scaling dimension of a CFT operator that appears in the expansion of the continuum limit of O α . The prefactor C α may include spatial oscillations, and further details are given in Section 3.4. Surprisingly, the set O α also act as order parameters for critical phases in a sense that we explain following Theorem 4.

2.4.
Main results. To fix notation, let us write: N p is the order of the pole at the origin, which is also the range of the longest non-zero coupling to the left. The number 4 of zeros inside, on and outside the circle are denoted N z , 2c, N Z respectively, and ρ is a real number. Since the t α are real, all zeros are either real or come in complex conjugate pairs. We first state results for the gapped case. Firstly we have that the correlators O α (1)O α (N + 1) form a complete set of order parameters for the gapped phases of H BDI .
Theorem 1a. In the phase (ω, c = 0, Σ) we have The non-zero constant is given in Theorem 1b. The value of the sign Σ may be inferred by the presence or absence of a (−1) N oscillation in this correlator.
Theorem 1b. In the phase (ω, c = 0, Σ), the non-universal value of the order parameter is given by Thus from the decomposition (9) we can read off ω = N z − N p and calculate the order parameter through the detailed values of the zeros. We discuss the mathematical form of the order parameter in Appendix E.
The next results show that the length scale in gapped phases is set by ξ = 1/|log|ζ || where ζ is any zero that maximises the right hand side of that equation (see Figure 2 for illustration). The set of ζ that are optimal in this way we call closest to the unit circle; we will always mean this logarithmic scale 5 when we talk about distance from the unit circle. The following results will be stated for 'generic cases' -we argue that these cases are typical in Appendix A.2. Now, in the phase ω we then have that ξ α , as defined in (7), is equal to ξ |ω−α| (for α = ω) -this is a consequence of: Theorem 2. If the system is in the phase (ω, c = 0, Σ) then, in generic systems, we have the large N asymptotics m ∈ Z is a known constant and M (N ) is a known |ω − α| × |ω − α| matrix. The elements of this matrix have magnitudes that depend algebraically on N -in particular, det M (N ) = Θ(N −δ ) for some δ > 0. Generic systems are those where the nearest zero(es) to the unit circle is either a single real zero, or are a complex conjugate pair of zeros.
The analysis we give extends to exceptional cases -more than two closest zeros will almost always give the same ξ α , but if one has multiplicity, ξ α may be controlled by the next-closest zero. The ξ α are always upper bounded by ξ, and in fact this bound is saturated in all exceptional cases except when there are mutually inverse closest zeros. See the discussion in Section 6.2 for full details.
The form of det M (N ) derived in Section 6.2.3 allows for some further general statements. Firstly, if there is one real zero nearest to the circle, then det M (N ) is real and does not oscillate with N . The algebraic factor depends non-trivially on |ω − α|, as demonstrated in Table 1 for the case that |ζ | < 1. If |ζ | > 1 then the second and fourth columns of Table 1 should be interchanged (and the definitions of λ and κ change in the obvious way based on the formulae in Propositions 1 and 3).
If there are two complex zeros nearest to the unit circle then det M (N ) is real but can contain O(1) oscillatory terms such as sin(N arg(ζ )) (these oscillations may, however, not appear in the leading 4 We consider a multiset of zeros {ζj} and allow ζj with different index to coincide. This makes the counting unambiguous. 5 That is, 1/|log |ζi||. This gives the natural length scale set by each zero, since the set of these lengths is invariant Table 1. The value of det(M (N )) in the case that there is one zero closest to the unit circle, and that zero is inside the circle. The constants κ and λ are independent of N and defined in Propositions 1 and 3.
We complete our analysis of gapped models with a result for the asymptotic approach to the value of the order parameter. In particular, we prove that ξ ω = ξ/2, following from: Theorem 3. In the phase (ω, c = 0, Σ) and in generic systems, we have for large N The factor B N is given implicitly in the proof and satisfies |B N | = O(1), m ∈ Z is a known constant.
We now discuss results for the gapless phases. In critical chains the phases in the BDI class described in the bulk by a CFT and connected to a stack of translation invariant chains with arbitrary unit cell are classified by the semigroup Z ≥0 × Z: they are labelled by the central charge c ∈ 1 2 Z ≥0 and topological invariant ω ∈ Z. The proof, using the f (z) picture, is given in [13]. Our present interest is confined to translation-invariant Hamiltonians that lie in one of these phases, and our next result gives the scaling dimension of the infinite class of operators O α . A graphical representation of this theorem is given in Figure 2.
[x] denotes the nearest integer to x.
Equation (14) is independent of the choice in rounding half-integers, although for later notational convenience we define it to round upwards in that case. In Section 7 we prove Theorem 4 on the way to the more detailed Theorem 10. That theorem gives the full leading order term in the asymptotic expansion of O α (1)O α (N + 1) at criticality, including nontrivial oscillatory factors that are helpful in identifying lattice operators with fields in the CFT description. We give a discussion of this CFT description in Section 3.4. A similar result holds when we have degenerate zeros on the unit circle, as long as every degeneracy is odd. The theory does not give results for the case that we have any zero of even degeneracy.
In the gapped case, Theorem 1a makes a simple link between measuring O α (1)O α (N + 1) and learning ω -one simply looks for the value of α with long-range order. It is not immediately obvious how to generalise this to critical models. Theorem 4 shows that, as in the gapped phases, the behaviour of correlation functions allows us to see marked differences between different critical phases. In particular, the link between lattice operator and the operators that dominate its CFT description changes discontinuously at a transition between critical phases (discussed in detail below). In the critical case one can determine c and ω, but this requires information about more than one correlator. Inspecting the form of equation (14), displayed in Figure 2, one concludes that it is not necessary to have access to the scaling dimensions of all O α in order to determine the phase. One method would be to measure the scaling dimensions of {∆ α , ∆ α±1 , ∆ α+2 } for some convenient α, and form the set of δ α = ∆ α+1 − ∆ α . This difference is equal to [(α − ω)/2c − 1/2] -this means that δ α is a constant integer on plateaus of width 6 2c, and that neighbouring plateaus differ in value by one. If the δ α are all different 7 then we must have c = 1/2 and ω can be determined easily using ω = α − δ α . Otherwise, one should then measure further scaling dimensions until the width of the constant plateau (equal to 2c) is found. Once c is known, ω may be determined: on the edge of the plateau we have ω = α − 2cδ α . Inferring the critical phase through these scaling dimensions is analogous to distinguishing the gapped phases through the string order parameter. If our model is taken to represent a spin chain then the O α are local for α odd. In Appendix C we show that it is possible to recover both c and ω using scaling dimensions of local operators on the spin chain. Moreover, in gapped chains one can use the universality of the gapped correlations to similarly infer ω from knowing only two correlation lengths; this is explained in Section 3.1.
2.5. The dual spin chain. Our results apply not only to H BDI but also to certain spin-1/2 chains. We briefly review this correspondence so that the reader can have both pictures in mind, and to help us make links to the literature in the next section. We write the Pauli operators as Define X n , Y n , Z n as the operators X, Y , Z acting on the nth site (and tensored with identity on all other sites). A class of translation-invariant spin chains is given by Hamiltonians of the form As before, we only allow a finite sum over α, have t α ∈ R and take periodic boundary conditions. This is the class of generalised cluster models. Note that this includes the quantum Ising, XY and cluster models as special cases. In Appendix A.3 we give details of the Jordan-Wigner transformation that relates H spin to H BDI . The main point is that our results for the behaviour of O α (1)O α (N + 1) apply equally well to the spin chain. The expressions for O α in terms of spin operators are given in Table 2. Some of these operators appeared in the recent works [42,43].
2.6. Relation to previous work. The key reference related to our results for gapped models is the classic paper of Barouch and McCoy [32]. There the authors study the XY model which is the spin model equivalent to (2) with non-zero t 0 , t 1 , t −1 only (and hence f (z) depends on two zeros). The section of that paper on zero temperature correlations contains results for O α (1)O α (N + 1) for α = 1, −1 in the phases ω = 0, 1 that include what one would obtain from our theorems. Beyond that, the paper [44] includes a calculation of the value of the order parameter for α = −1, 2 in the special case that f (z) = z 3 − λ. Some portion of the phase diagram for −2 ≤ ω ≤ 2 is mapped out in reference [45] where order parameters are identified and calculated numerically. Several papers, for example [46,47], study spin models with competing 'large' cluster term and Ising term (i.e. non-zero t α , t −1 and t 0 ). In these cases winding numbers are identified, but not order parameters or their values. Our computation giving Theorem 1b is novel, extending previous calculations by addressing the full set of translation invariant models in the BDI class which require f (z) with an arbitrary (finite) set of zeros. Moreover, this generality shows the robustness of these order parameters throughout the phase diagram.
As mentioned, several papers have identified the form of the order parameters for |ω| ≤ 4 in the spin language. Equivalent fermionic order parameters are easily found using the Jordan-Wigner transformation and the paper [41] includes the fermionic O α for |α| ≤ 2 as well as discussing the general case. In our work we prove that the intuitive general case holds by linking these order parameters to the generating function f (z) and matching the winding number of f (z) to the 'unwinding number' of each correlator.
There are many works that study correlations in particular quantum phase transitions in our model. Again reference [32] should be mentioned, along with [48], as seminal early works that derived critical behaviour for correlators O α (1)O α (N + 1) with |α| ≤ 1 at the c = 1/2 Ising transition. Transitions with higher c include the c = 1 XX model that is a standard model in physics [49], and the same correlators were analysed in reference [50] using the mathematical methods found below. We also mention the quantum inverse scattering method as a tool for calculating scaling dimensions in certain cases [51] .
An isotropic spin chain is invariant under spin-rotation around the Z axis. In our fermionic model H BDI , this manifests as invariance under spatial inversion H α ↔ H −α , and hence is a model for which f (z) = f (1/z). This relation implies that ω = −c. The correlators O α (1)O α (N + 1) with |α| ≤ 1 in isotropic models with general c and ω = −c were derived in references [52,53] using the same methods as this paper. Our results go further by studying a wider class of models, including critical phases with general (c, ω), as well as a wider class of observables: O α (1)O α (N + 1) for all α. This allows us to observe that from knowledge of the scaling dimensions of these operators, one can identify the topological invariant ω.

Physical context and discussion
3.1. SPT phases and string orders: universality and symmetry fractionalisation. In Section 2.4, we saw that for both gapped and critical systems in our class of models, one can measure the topological invariant ω by looking at the string order parameters O α . For gapped phases, one needs to find the value of α for which there is long-range order (Theorem 1a), whereas in the critical case, one uses the scaling dimensions (see Theorem 4 and the discussion following it). The existence of topological string order parameters for critical phases is novel. However, even for the gapped phases that we consider, the string order parameters are unusual. This is for two reasons. Firstly, the usual justification for string orders relies on the concept of symmetry fractionalisation, which arises in the classification of interacting SPT phases and is usually not employed in the classification of noninteracting topological insulators and superconductors. Secondly, even in the interacting case, phases which are protected by anti-unitary symmetries do not give rise to the kind of string order parameters we discuss in this work. Bridging this gap is the purpose of Section 3.1.2. However, first we discuss the remarkable result that the correlations in the gapped phases exhibit universal properties.
3.1.1. Universality. In the gapped phases, O α has correlation length ξ α = ξ |α−ω| (if α = ω), see Figure  2 (we assume the generic case for this discussion). This means that although ξ depends on microscopic properties (like the position of the zeros of f (z)), the ratio ξ α /ξ β depends only on ω and is hence constant in each phase. This has interesting consequences. In principle, to determine the topological invariant of a gapped phase, one has to find an α such that | O α (1)O α (N ) | tends to a non-zero limit as N → ∞. This requires going through an arbitrarily large set of observables. Surprisingly, it is sufficient to measure only, for example, two correlation lengths ξ α 1 and ξ α 2 (for the observables O α 1 and O α 2 ) for any fixed choice of α 1 and α 2 satisfying |α 1 − α 2 | ∈ {1, 2}. To see this, note that there are three cases. Firstly, if one finds long-range order for either α 1 or α 2 , then ω is known. Secondly, if ξ α 1 = ξ α 2 , one knows that ω = α 1 +α 2 2 . In any other case, α 1 and α 2 will either both be larger or smaller than ω, such that This can be uniquely solved, giving ω = . The above shows that, using universality, one can replace an infinite number of observables by just two. However, it has an even more surprising consequence in the spin language. If we choose α 1 = 1 and α 2 = −1, then this corresponds to the correlation lengths ξ X and ξ Y of the local observables X n and Y n . This fully determines the invariant This means that one can distinguish, for example, the trivial paramagnetic phase from the topological cluster phase 8 by measuring the decay of correlation functions of local observables. This is truly unusual and presumably an artifact of looking at spin models that are dual to non-interacting fermions. It would be interesting to investigate such ratios between correlation lengths in interacting models and determine whether this is a measure of the interaction strength between quasi-particles.

Symmetry fractionalisation.
To contrast our analysis to the standard justification for string orders, we briefly repeat how string order parameters arise within the context of symmetry fractionalisation. It is worth emphasising that the known constructions for string order parameters of the type that we discuss are only for SPT phases which are protected by unitary symmetries [36,54]. Let U be some on-site unitary symmetry, i.e. [U, H] = 0 with U = n U n . Consider the operator U X = n∈X U n where X is some large line segment of length l (see Figure 3). If l ξ, then deep within X, U X looks like a bona fide symmetry operator. Hence, it is only near the edge of X that U X can have a non-trivial effect. In other words, if |gs is the ground state, then effectively U X |gs = U L U R |gs , where U L and U R are operators that are exponentially localised near the boundary of the region X. This can be made rigorous using matrix product states [55,56]. This phenomena is known as symmetry fractionalisation and is the essential insight that led to the classification of (interacting) SPT phases in 1D [9,11,12,57]. To illustrate this, consider the case with a symmetry group Z 2 × Z 2 , generated by global on-site unitary symmetries U and V . Since U n and V n commute on every site, we have that U X V = V U X . Moreover, U X = U L U R (when acting on the ground state subspace), implying that U L and V have to commute up to a complex phase. Using that V 2 = 1, we arrive at This defines a discrete invariant which allows to distinguish two symmetry-preserving phases. (One says that the phases are labelled by the inequivalent classes of projective representations of Z 2 ×Z 2 .) In fact, one can show that in the phase where ω = 2, the negative sign implies degeneracies both in the entanglement spectrum and, for open boundary conditions, in the energy spectrum [8].
Here we will not go into such details, but refer the interested reader to the review in reference [58]. Instead, we consider the effect on correlation functions.
One can consider the string correlation function gs| U X |gs = U X = U L U R . Due to locality, we have that U X ≈ U L U R for l ξ. Since SPT phases do not spontaneously break symmetries, we Hence, we conclude that the string correlation function U X has to be zero if ω = 2. Equivalently, measuring U X = 0 implies that ω = 0, such that one calls this a string order parameter for the trivial phase (ω = 0). Analogously, one can construct a string order parameter for the non-trivial phase: if V, W are local operators anticommuting with V , then by repeating the above argument, we conclude that VU X W = 0 implies that ω = 2 (with V (W) localised near the left (right) of region X). Note that these string order parameters always work only in one direction: there is no information if one measures them to be zero. This is in striking contrast with the string order parameters we found for our non-interacting class of models.
Let us make this discussion more concrete with an example, where U = P = n Z n and V = P odd = m odd Z m . Two models with this symmetry are H 0 = n Z n and H 2 = − n X n−1 Z n X n+1 . One can calculate that their symmetry fractionalisations are ω = 0 and ω = 2, respectively. The above tells us that N n=1 Z n is a string order parameter for the trivial phase. Similarly, taking V n,n+1 = Y n X n+1 -which indeed anticommutes with P odd -then , is an order parameter for the topological phase ω = 2. In this case, the string order parameters we have derived -with respect to Z 2 × Z 2 -for the trivial phase (connected to H 0 ) and the topological cluster phase (connected to H 2 ) happen to be the same as we encountered in the non-interacting case -with respect to the P and T symmetries -see Table 2.
Can we make a connection with our non-interacting classification and the concept of symmetry fractionalisation? For this it is easiest to work in the fermionic language. It is known that if one studies the fractionalisation of only the P and T symmetries, then there are only eight distinct phases [59]. However, since our model is non-interacting, the P and T symmetries imply an additional structural symmetry: the Hamiltonian can only contain terms which have an equal number of real and imaginary Majorana modes. This implies that if we have any operator which has a well-defined number of real minus imaginary Majorana operators (e.g. γ n γ n+1 would have 'charge' two), then the Hamiltonian time evolution would conserve this. To see how this is useful, consider a fixed-point model H α = i /2 n∈sitesγ n γ n+α . It is a simple exercise to check that for the symmetry fractionalisation of P = P L P R , we have that P L has charge −ω (and P R charge ω). By the aforementioned argument and the concept of adiabatic connectivity, these charges of P L and P R should be stable throughout each gapped phase. It is easy to see that V = 0 for any operator whose charge is non-zero. Hence, in this way we are naturally led to consider γ 1 · · · γ α (γ α+1 γ α+1 · · ·γ N γ N )γ N +1 · · ·γ N +α : this operator can only have long-range order if α = ω.
We thus see that the notion of symmetry fractionalisation does give a good plausibility argument for Theorem 1a. Do note that the latter result is stronger, in particular it shows the converse (i.e. that the long-range order can never become accidentally zero). However, the power of symmetry fractionalisation is that it implies that these string orders, edge degeneracies and other topological properties should be stable against any interactions which respect the above symmetry 9 . This would be interesting to explore further.

3.2.
Majorana modes: localisation length versus correlation length. In reference [13], we showed that if ω > 0, then the system has ω Majorana zero modes per edge. More precisely, to each of the ω largest zeros {z i } i=1,··· ,ω of f (z) within the unit disk, we associate an hermitian operator γ i L , all of them mutually commuting and squaring to the identity. Moreover, these operators commute with T and are exponentially localised near the left edge, with respective localisation lengths −1/ ln |z i |. (The same is true for the right edge, where they anticommute with T .) The crucial property which makes them so-called zero-energy modes, is that they commute with the Hamiltonian (up to a finite-size error which is exponentially small in system size). Hence we have ω mutually anticommuting symmetries, from which one can show that to each edge we can associate a √ 2 ω -fold degeneracy 10 . This is to be contrasted with the fact that the ground state is unique for periodic boundary conditions. This is a characteristic property of topological insulators (and more generally, symmetry-protected topological phases) in one spatial dimension. This is well-known for the gapped phases of the BDI class, but the proof in reference [13] shows that this analysis goes through when the bulk is critical (i.e. c = 0). Reference [60] notes the link between the localisation length of the Majorana edge mode and the behaviour of bulk spin correlations in a model equivalent to the XY model, and conjectured that this is a general phenomenon. Here we simply point out that the largest localisation length of a Majorana mode (if present) need not coincide with the bulk correlation length. Indeed, the latter is determined by the zero of f (z) closest to the unit circle. In particular, if ω > 0, the localisation lengths of the Majorana modes are determined by zeros within the unit disk, whereas it could certainly be that a zero outside the unit circle dominates the bulk correlation length. This disagreement with the two quantities is consistent with the observation in [13] that one can tune a gapped phase to a critical point whilst (some of) the edge modes remain exponentially localised. This discussion for ω > 0 holds also for ω < −2c (if we interchange the words 'left' and 'right'), with the edge modes now being associated to zeros outside the unit circle. As mentioned above, spatial leftright inversion corresponds to f (z) ↔ f (1/z), which at the level the topological invariant corresponds to ω ↔ −ω − 2c. For any other value of ω, there are no edge modes [13].
3.3. Critical exponents. Critical exponents encode how different physical quantities diverge upon approaching a phase transition. In the classical case, the tuning parameter is usually the renormalised temperature τ = T −Tc Tc . In the current quantum setting, there is an equally natural 11 tuning parameter ε ∈ R: if f (z) represents a gapless Hamiltonian, then f ε (z) := f (z(1 − ε)) interpolates between two gapped phases. One can think of this as shrinking (ε < 0) or expanding (ε > 0) the radial component of the zeros of f (z). We will work in the case where the system is gapless at ε = 0 with 2c zeros on the unit circle, allowing for a multiplicity m (which we take to be the same for every distinct zero). This means we have 2c/m distinct zeros on the unit circle. We emphasise that the expressions in the remainder of this section are derived only in the case of uniform multiplicity. Our results allow for an analysis of the general case, but we do not wish to pursue this here.
We derive four critical exponents: the anomalous scaling dimension, η, defined by the scaling of the order parameter, Ψ, at the critical point ( Ψ(1)Ψ(N ) ∼ 1/N η ); ν, which encodes the divergence of the correlation length (ξ ∼ |ε| −ν ); the dynamical critical exponent, z, that relates how the correlation length ξ diverges relative to the characteristic time scale defined by τ , the inverse energy gap (τ ∼ ξ z ); and β which relates to the decay of the order parameter (Ψ ∼ |ε| β ).
Exactly at the critical point, we have given explicit results above only for m = 1. In particular, Theorem 4 gives us that η = c/2. In the special case of c = 1/2, we recover the well-known result η = 1/4. For m > 1 and odd, we can use equation (113) from Appendix F to easily derive that η = mc/2. The other three exponents are defined away from criticality where we can allow for any m ≥ 1. The correlation length is determined by the nearest zero, such that ξ ∼ 1/| ln |1 + ε|| ∼ 1/|ε|. Hence, ν = 1; independent of c and m. The energy gap, min z∈S 1 |f (z)|, depends on the location of all zeros. However, close to criticality, we need to care only about the zeros describing the transition. Moreover, each of the distinct zeros has a local minimum, which all scale the same way. So without loss of generality, we can consider f (z) = (z − (1 + ε)) m . The gap scales as ∼ |ε| m , such that τ ∼ ξ m with dynamical critical exponent z = m. This is consistent with m = 1 being described by a CFT with central charge c, since that implies z = 1. Lastly, we consider the order parameter given in Theorem 1b. For each distinct zero, the order parameter has a factor |1−(1+ε) −2 | m 2 /8 ∼ |ε| m 2 /8 ; all other terms do not go to zero with ε. Combining 2c/m such factors, we have Ψ ∼ |ε| In the special case of the Ising transition, m = 1 and c = 1/2, this reduces to the well-known result β = 1/8.
3.4. CFT and continuum-lattice correspondence. We now explain how certain features of Theorem 4 (and Theorem 10) fit in to a CFT analysis of the critical point. This section is not rigorous, the aim is to complement the mathematical proofs with a perhaps more intuitive physical picture. We also use Theorem 10 to make claims about passing from the lattice to the continuum description of the operators O α .
If our system has 2c non-degenerate zeros on the unit circle, then one can argue that the appropriate low-energy theory is a CFT built from 2c real, massless, relativistic free fermions [62]. Briefly, one 11 Moreover, one can check that this agrees with τ under the quantum-classical correspondence.
can linearise the dispersion relation |f (k)| about all of its zeros on the circle (Fermi momenta), and each local linearised mode is such a fermion. Moreover, since |f (k)| = |f (−k)| we can combine the real fermions from a pair of complex conjugate zeros to form a complex Fermion with central charge one. This is helpful as complex fermions can be bosonised using the methods given in [49,63,64]. In general, then, we have a low energy theory of an even number of complex fermions and either 0, 1 or 2 real fermions (located at k = 0 or π). The central charge is always equal to 2c.
3.4.1. High level analysis. In reference [65], general CFT considerations are applied to integrable models to find asymptotic correlation functions of local operators A n (x) that create fixed numbers, n j , of quasiparticles of type j (i.e. excitations near momentum ±k j ). For equal times, these take the form for scaling dimensions ∆ (j) n , Fermi momenta k j and the sum is over sets of m j ∈ 2Z . The amplitude C An depends on the appropriate form factor [65]. The oscillatory term comes from a multiplicative factor e −ixδp when A gives an intermediate excitation with momentum δp -these are non-zero when we have particle-hole excitations where the particles and holes are at different Fermi points. Relevant discussion is found in references [51,65,66]. The O α that we consider are 'square-roots' of local operators in the sense that the product O α (n)O α (n+m) (for small m) is local on the lattice and has an expansion dominated by such A. One may then expect that correlations of O α will have the same form as (18), but with m j ∈ Z. This is verified in Theorem 10, apart from the possibility of an additional (−1) x oscillation. This is needed when the low energy degrees of freedom are modulated by this oscillation. The constant in Theorem 10 implicitly gives form factors of the relevant fields [65,67].

3.4.2.
CFT operator correspondence for one real fermion. We now consider more details of the CFT description and the correspondence with canonical continuum operators. The following discussion will be in terms of the spin chain dual to the fermionic system, as discussed in Section 2.5.
If our model has one zero on the unit circle (which must be at z = ±1), then our continuum limit has c = 1/2, and is hence described by the Ising CFT [68]. This contains two primary operators with scaling dimension 1/8, the 'spin' and 'disorder' operators, denoted by σ and µ respectively. In the usual lattice Ising model, σ is the field corresponding to the local order parameter of the neighbouring ordered phase and µ is a non-local string order parameter of the neighbouring paramagnetic phase. Our results on the lattice, in particular Theorem 4, show that a system with winding number ω has two operators with this dominant scaling dimension: O ω and O ω+1 . We identify O ω with σ when ω is odd and µ when ω is even; O ω+1 is the other field. This attribution is consistent with locality and parity symmetry of the operators. Importantly, we conclude that the operators on the lattice that have overlap with the dominant primary scaling fields depend on ω. Indeed they change discontinuously at a transition between two critical phases described by c = 1/2 CFTs with distinct values for the topological invariant ω -we give an example of this in Section 3.4.4. Other operators O α for α neither ω nor ω + 1 will be dominated by descendants of σ for α odd and of µ for α even (the lattice operators have dimension 1/8 + j for some j ∈ Z + , so we should take CFT descendants at level j).

3.4.3.
CFT operator correspondence for one complex fermion. Let us now consider the case c = 1 with a complex conjugate pair of zeros, at e ±ik F , and a U(1) symmetry generated by S z tot = 1 /2 i Z i (the standard model in this class is the XX spin chain [49]). These are isotropic models with f (z) = f (1/z), which implies that ω = −c; we discuss other values of ω later.
The fermionic Hamiltonian may be bosonised as described in [63,69,70] and [49,Chap. 20], passing also to a continuum limit. We denote the resulting bosonic fields by θ(x) and ϕ(x), such that where δ(x) is the Dirac delta function. The vertex operator e iθ(x) creates a localised charge of the aforementioned U(1) symmetry, and ∂ x ϕ(x) is a density fluctuation: the total density is ρ(x) = (k F + ∂ x ϕ(x))/π. Vertex operators of the form e i(λ 1 θ(x)+λ 2 ϕ(x)) (well defined for λ i ∈ Z) have scaling dimension (λ 2 1 + λ 2 2 )/4 [49,71]. We now consider operator correspondences. Bosonisation will not fix constant coefficients of the operators, so we will usually suppress them in the following. Note, however, that hermiticity and symmetries constrain certain coefficients. We may also need to include an additional antiferromagnetic oscillation.
Firstly, following [63,69], we note that ρ(x) generates the U(1) symmetry. Hence, we can make the formal identification The Jordan-Wigner strings follow by setting χ = π and truncating the sum and integral respectively. We can then make the correspondence (see, for example, [69,70] and [49,Chap. 20]): The fermionic creation operator creates a U(1) charge and is multiplied by a Jordan-Wigner string, this leads to standard expressions for the dominant contributions to the right (+) and left (−) moving continuum fermion fields Time-reversal symmetry swaps these right and left moving fields (one may check the lattice expansion given in, for example, [70]) and so we have that ϕ → ϕ and θ → −θ under T . We can also confirm that the right-hand-side of (21) is hermitian and does not transform 12 under the U(1) or T , as required for the Z-string. The site at −∞ may appear problematic in isolation, but we always consider two-point correlators and thus the infinite string to the left will drop out (similarly in the continuum one may take correlation functions to avoid considering the boundary). The meaningful correspondence here is: where we now suppress constant coefficients. Now, by considering the U(1) action and time-reversal symmetry, we see that σ ± n → e ±iθ(n) (with real coefficients). We then have that X n → cos(θ(n)) Y n → sin(θ(n)).
Operators O α for |α| > 1 are more involved. We consider first the family of local spin operators, with α odd. In particular, let us analyse O 3 (n) = X n Y n+1 X n+2 . Note that any product of three neighbouring X or Y operators can be expanded as a linear combination of terms σ ± n σ ± n+1 σ ± n+2 , with all possible sign combinations present. These products transform separately under U(1) with charge m ∈ {−3, −1, 1, 3} equal to the sum of the signs. The dominant terms would be e ±iθ(n) , with subdominant contributions from e ±3iθ(n) and products of e ±iθ(n) and e ±3iθ(n) with derivatives of θ(n). It is possible that the coefficient of the dominant term and the first few subdominant terms vanish due to destructive interference -Theorem 4 verifies that this occurs and in fact that O ±3 (n) has the same scaling behaviour as e ±3iθ(n) . Details of a formal calculation in terms of the CFT operator product expansion are given 13 in Appendix B. Intuitively, O ±3 is dominated by terms that create three charges, as well as the remnants of several terms that create one charge but destructively interfere with each other. We can generalise to all odd α the conjecture that O α (n) has an expansion of the form D α,m (θ) is constant for m = ±α, and for other values of m contains products of derivatives of θ such that the scaling dimensions of all terms match the extremal terms e ±iαϕ(n) . This is consistent with both Theorem 10 and calculations of the type in Appendix B.
12 Under a U(1) rotation through χ ∈ R, lattice operators are conjugated by exp(iχ n Zn/2) and fields transform as . 13 Parenthetically, this calculation indicates that other triples of neighbouring X and Y operators will scale as e iθ ; we have confirmed this in numerical simulations.
Even α follow the same pattern, except that we should always include a Jordan-Wigner string: the dominant term of which is given in (21). Hence we take O α (n) → |α| m=0 e i(|α|−2m)θ(n)+i(k F n+ϕ(n)) + e i(|α|−2m)θ(n)−i(k F n+ϕ(n)) D (α,m) (θ(n)) + . . . α even, (26) where, as above, D α,m (θ) is constant for m = ±α and contains appropriate numbers of derivatives such that all terms have the same scaling dimension. This is consistent with Theorem 10 and CFT calculations. In all case we should allow multiplication by a global antiferromagnetic oscillation as discussed in the previous section.
The preceding analysis required the U(1) symmetry of isotropic models as a starting point. In reference [58], generalised Kramers-Wannier dualities were discussed that map between our models. One class of transformations swap models such that f (z) ↔ z n f (z) for some n ∈ Z. If f (z) is isotropic, then z n f (z) is not -this allows us to extend the preceding correspondence to anisotropic models. We should shift the labels α → α + n while not shifting 14 the right hand sides of the correspondences (25) and (26). For example, in the cluster model H = − n X n Z n+1 X n+2 + n Z n we identify X n with e i(k F n+ϕ(n)) + e −i(k F n+ϕ(n)) . Moreover, this duality allows us to map an isotropic model to a representative of each phase (c = 1, ω). Anisotropic models that are dual to isotropic models in this way have an appropriately transformed U(1) symmetry. It is then nontrivial to extend this analysis to general anisotropic models. Theorem 10 indicates, however, that this correspondence should continue to hold. Since the above argument does not make use of the fact that our underlying lattice model is non-interacting, we expect the correspondence to persist 15 in interacting models if the U(1) symmetry is preserved. However, if the U(1) symmetry is broken, we see no reason to expect it to continue to hold (in contrast to the non-interacting case).

3.4.4.
Example: transition between topologically distinct critical phases. To complement the discussion so far, we consider the example Tuning −1 ≤ λ ≤ 1 we interpolate between two critical phases, with a transition at λ = 0. For λ = 1 the system is the standard critical Ising chain with H = −( j X j X j+1 + Z j ). Models with λ > 0 are in the same phase: (c = 1/2, ω = 0). For λ = −1 we have H = −( j X j Z j+1 X j+2 − X j X j+1 ). Models with λ < 0 can be smoothly tuned to this model and have (c = 1/2, ω = 1). For λ = 0, we have the c = 1 transition between topologically distinct critical phases, with H = − 1 2 ( j X j Z j+1 X j+2 + Z j ). Table 3 gives the behaviour of the scaling dimensions of the most dominant O α as the system crosses the transition. We emphasise that while both sides of the transition are described by an Ising CFT, the scaling dimensions of the lattice operators change discontinuously. Note that the c = 1 model has two real zeros so the CFT discussion above does not quite apply -a similar bosonisation scheme does work, as applied to a doubled Ising model in reference [68].
3.4.5. CFT operator correspondence with c complex fermions. For higher values of c, we work formally with the operator content of the product of 2c copies of the Ising CFT. Note that it has been shown that spin models with f (z) = ±z ω (z 2c ±1) are described at low energy by so(2c) 1 WZW models [45,46], although a lattice-continuum operator correspondence has not been made. We can smoothly connect 16 any critical model in H BDI to that subset of models [13].
Let us suppose for now that we have no real zeros and c complex conjugate pairs of zeros (we order the zeros so that k i = −k 2c−i ). Then we have c canonical complex fermions which can each be 14 Attributing the mutually conjugate fields ϕ(n) and θ(n) requires some additional convention. Note that their scaling dimensions do not distinguish them since our models have Luttinger parameter equal to one [49]. We fix that ϕ(n) is the field that is attached to oscillations e ±ik F n on the lattice. 15 Note that the scaling dimensions of the vertex operators will depend on the interactions, so the continuum operators in our expansions (including derivative terms) will not necessarily all have the same dimension. The dominant operators should be identified as a subset of these. 16 That is, along a path where the CFT data varies smoothly.
bosonised as described above to give a set of fields θ i (x) and ϕ i (x). The relevant vertex operators are of the form where µ i and ν i are half-integer and ∆ µ,ν = j (µ 2 j + ν 2 j )/2 (we suppress Klein factors [64]). These operators act nontrivially on all fermionic sectors, and have a minimal scaling dimension of c/4 (notice that this coincides with the smallest scaling dimension of the O α , given in Theorem 4). The operator ϕ i (x) appears physically as an integrated charge density, so it is accompanied by k i x. The conjectured expansion of lattice operators O α goes through roughly as above, however, the leading order terms have charge distributed evenly throughout the different sectors. In particular, observe that the charge-two operator e i(θ 1 (x)+θ 2 (x)) with ∆ = 1/2 dominates the charge-two operator e i(2θ 1 (x)+ϕ 2 (x)) with ∆ = 5/4 (and indeed the charge-two operator e i2θ 1 (x) with ∆ = 1). Then, as argued in [53], in isotropic models σ + = O 1 + iO −1 will be dominated by operators τ µ,ν with j µ j + ν j = 1 (charge condition) and |µ i − ν j | ≤ 1 (dominance condition). These conditions give a sum of terms that are products of e ±iθ or e ±i(ϕ+k j x) in each sector, and hence that can be distinguished by the presence or absence of oscillatory factors e ±ik j x . This is analogous to the sum over Fisher-Hartwig representations (see Section 5) that we derive in Section 7.1, and the relevant oscillatory factors are confirmed in Theorem 10 (indeed, each such term is represented in the final result).
To extend this to models with non-zero ω and to operators with |α − ω| > 1, we consider first the operators τ µ,ν with j µ j + ν j = c + ω − α (maximal charge condition) and |µ i − ν j | ≤ 1 (dominance condition). This gives a set of operators that we expect to dominate the continuum limit of O α . However, as in the c = 1 case, we expect that we should include terms where maximally charged operators e iRθ are substituted with a series of terms e i(R−2m)θ | m=1,...,R multiplied by derivatives, such that each term has the same scaling dimension. The relevant scaling dimensions and oscillatory factors are confirmed in Theorem 10.
Example: Consider the case c = 2 with ω − α = 1. We conjecture the correspondence: Note that all terms have ∆ = 3/2, and all oscillate as either e ±ik 1 x or e ±ik 2 x as expected. Although O ±α contain the same dominant operators, the coefficients are not expected to be the same in general. Theorem 10 indicates that the coefficients are different in the general case, since (96) is not symmetric under taking sign-reversed Fisher-Hartwig representations.
In the case that c ≥ 1 and any zero is real, we do not conjecture the operator correspondence. We expect that similar arguments could work after bosonising a doubled model -this is performed for the c = 1/2 case in [72]. Note that Theorem 10 does not distinguish the case of two real zeros from two complex conjugate zeros at the level of scaling dimension.
3.5. Entanglement scaling. The entanglement entropy of a subsystem is another physically important quantity. Let ρ A be the ground state reduced density matrix on sites 1 up to N and consider asymptotics in large N after taking the length of the (periodic) chain to infinity. The most general results for isotropic critical chains in our class are given in [40,73]. Having identified the correlation length in gapped chains, derived from the nearest zero to the unit circle, it is interesting to consider the following theorem adapted from [74]. (2)) such that 2c of the zeros approach the unit circle, and that the limiting chain has no degenerate zeros on the circle. We label these approaching zeros by ζ j , noting that ζ j can be either inside or outside the circle, and is either real or a member of a complex conjugate pair. Then the entanglement entropy of a subsystem of size N , in the limit N → ∞, has the following expansion as |ζ j | → 1:

Theorem 5 (Its, Mezzadri, Mo 2008). Consider a sequence of gapped chains (as defined in equation
Note that the O(1) term is constant with respect to all the zeros that approach the unit circle (which are allowed to approach independently). Now, let us consider a sequence of models with a set of 2c zeros that approach the unit circle; for notational convenience let us fix them to be complex zeros outside the unit circle, other cases lead to the same result. Let this set of approaching zeros be specified by: {e ±iφ 1 e 1/ξ , e ±iφ 2 e t 2 /ξ r 2 , . . . e ±iφc e tc/ξ rc }, where φ i = φ j for i = j, t j > 1 and r j < 1, and we approach the circle letting ξ → ∞. The conditions on t j and r j ensure that a closest zero is e iφ 1 e 1/ξ for ξ large enough. Inserting into Theorem 5, we get that: limit. Having different rates of approach to the circle means that we are necessarily approaching a multicritical point and we see crossover behaviour in the entanglement scaling. A simple example that allows this behaviour is the approach to the c = 1 critical point with is infinitesimally close to a c = 1/2 line of transitions in the phase diagram of the XY model. This is reminiscent of the Calabrese-Cardy formula [75] that applies far more generally and gives asymptotics as the lattice spacing 17 a → 0 where c is the central charge of the underlying CFT and ξ is the (fixed) correlation length of the system under consideration. This may also be interpreted as a scaling limit ξ a [76], and the formula was confirmed in this sense for the XY model in [77]. Further relevant references are found in the review articles [76,78]. We see that equation (31) is equivalent to formula (33) in the vicinity of a regular critical point. At multicritical points the path approaching the transition is important, and the Calabrese-Cardy formula is expected to hold along RG flow lines in parameter space.

String correlators as determinants
We now begin the analysis necessary to prove the results given in Section 2.4. (3), we have that:

Fermionic two point correlators. After defining f (z) as in equation
where the Boguliobov quasiparticles are found by rotating the Bloch sphere vector 18 (c −k , c † k ) through an angle f (k)/|f (k)| about the x-axis, giving The sum over k goes over momenta k n = 2πn/L, although we always work in the limit where this sum becomes an integral from 0 to 2π. Details of this diagonalisation may be found in, for example, 17 Elsewhere in our work, units are fixed such that a = 1. 18 The c k are the Fourier transform of the lattice fermions from which we built the γn in equation (1). [13,39,40]. The ground state, |gs , is the vacuum for the quasiparticles η k , and from this we can easily calculate fermionic correlation functions -we refer the reader to [40] for details. We will use −iγ n γ m := gs| (−iγ n γ m ) |gs = 1 2π γ nγm = γ n γ m = δ nm as elementary correlation functions in the rest of the paper, noting that it is arg(f (z)), on the unit circle, that controls these correlations.
As an aside, note that for gapped chains the analysis of Section 6.2 allows us to find the large N asymptotics of −iγ n γ n+N . As explained in the discussion around Section 6.2.4, in generic cases this correlator will be Θ(N −K e −N/ξ ) where K ∈ {1/2, 3/2} is easily determined. For critical chains f (k) has jump discontinuities. Decomposing as in equation (86) and integrating by parts, we have that the fermionic two point function is Θ(1/N ) -this behaviour is as expected from the CFT description. (2) is quadratic, ground state expectation values have a Pfaffian structure. More precisely, suppose that we have 2N distinct and mutually anticommuting operators, A n , then:

Wick's theorem. Because the Hamiltonian
(−1) σ is the sign of the permutation that reorders the operators into each particular pairing. This expression is proportional to the Pfaffian of the antisymmetric matrix A m A n , and is a form of Wick's theorem that is given in reference [79].

String correlation functions.
Consider the two point correlation function of O α for α > 0: We now transpose further terms to put unlike Majoranas as nearest neighbours and apply Wick's theorem: For α < 0 an analogous calculation again leads to equation (41). Table 2 gives the spin operators Jordan-Wigner dual to the fermionic operators O α (n) and Table 4 gives the equivalent spin correlators for all α. A derivation is given in Appendix A.4. Notice that for odd α these operators and correlation functions are local in the spin variables, and for even α they are nonlocal; they are always nonlocal for the fermions. Understanding the asymptotic behaviour of the determinant (41) is the key to the results given in Section 2.

Toeplitz determinants
Several theorems for the asymptotic behaviour of large Toeplitz determinants are required to prove our results, hence we use this section to review them in detail. This section is intended to not only state the results but to give an exposition of how to use them in practice. The reader already familiar with these ideas can hence skip this section and refer back to it where necessary. Note that we reformulate and simplify the statement of some theorems appropriately for our application, the most general statements are available in the given references.
First, recall that an N × N Toeplitz matrix, T , takes the following 'translation-invariant' form: This matrix can be thought of as the N × N truncation of an infinite matrix, with element t −n on the nth descending diagonal. Consider a region of the complex plane, U , such that S 1 ⊆ U ⊆ C. A function t : U → C, integrable on the unit circle, generates a Toeplitz matrix through its Fourier coefficients: We refer to such t(z) as the symbol of the Toeplitz matrix 19 and denote the Toeplitz determinant of order N that is generated by t as i.e. it is defined simply as the determinant of the N × N truncated matrix generated by t. It is the analytic properties of t that govern the form of the asymptotics of this determinant as N → ∞. By inspecting equation (41), we see that O α (1)O α (N + 1) is, up to an oscillating sign, a Toeplitz determinant of order N generated by To go further, we consider the symbol on the unit circle z = e ik and attempt to factorise it as: Here e V (z) is called the smooth part of the symbol, which we take to mean that V (z) is analytic on the unit circle 20 . This implies that the winding number of exp V e ik is equal to zero. The Fourier coefficients of V e ik are and we define the Wiener-Hopf factorisation of e V (z) as: In our work, we will have three families of symbol to consider. The first case is t singular (z) = 1, which works when our symbol t(z) is smooth enough that its logarithm gives us an appropriate V (z). The second case is t singular (z) = z ω , this is needed to represent symbols t(z) that have an integral winding number ω. Finally, the third case represents symbols t(z) with sign-change jump discontinuities. Let ζ = e iθ and consider the function on the unit circle: For β half-integer this is piecewise proportional to i, with a sign change at z = ζ and at z = 1. To represent a sign-change only at z = ζ, we put t ζ,β singular (z) = z β g ζ,β (z), removing the jump at z = 1. Conversely a jump only at z = 1 would be represented simply by z β . Notice that any half-integer β represents the sign-change through g ζ,β , but the power of β that appears distinguishes the t ζ,β singular (z). 19 We will always go in this direction: from symbol to matrix. The reverse is possible providing the ti decay fast enough. 20 This smoothness requirement has been weakened by many authors; a strong result is given in [80], to which we refer the interested reader. For our purposes the strong condition of analyticity is acceptable.
The singular part of a function with several sign-change jump discontinuities can be decomposed as a product where all β j are half-integer, but note that now only the total j β j is fixed by the symbol we wish to represent -this redundancy has important consequences. As an example, consider the symbol s(k) = sign(cos(k)), this has jump discontinuities at z = ±i. Hence we should represent it by two β half-integer singularities, and the fact that there is no overall winding implies that β 1 = −β 2 . This gives a family of representations where n ∈ Z and the constant fixes the correct overall sign at z = 1.
With these ideas in place, notice that all three families of t singular (z) can be represented in the same way. If we use t singular (z) = z β g ζ,β (z) as a building block, then t singular (z) = 1 is the case ζ = 1, β = 0 and t singular (z) = z ω is the case ζ = 1, β = ω. Motivated by this discussion, we write down the canonical form of reference [80] for a symbol that is non-vanishing on the unit circle and has sign-change jump discontinuities: where for j = 1, . . . , m and 0 ≤ k 1 < . . . < k m < 2π, we have z j = e ik j , β j ∈ 1 2 Z and the function V (z) must be smooth as above. The factor j z −β j j is just a multiplicative constant and is there to align notation with [80]. Any β j in this expression must be nontrivial, hence the symbol has m jump discontinuities. Note that we allow m = 0 when the symbol is simply exp(V (z)), and the edge case z 1 = 1 has g 1,β 1 = exp(−iπβ 1 ). Our notation deviates slightly from reference [80], where a β 0 is associated to z = 1 even if there is no singularity there -this does not affect the adapted theorems we quote below. Now, as explained above, for a symbol t(z) with multiple jump discontinuities, there is an infinite class of different t canon (z) to which it is equal. In fact, if we find a single representation with a set of {β j }, we can find another representation by shifting each β j → β j + n j such that n j = 0; however, we may have to amend our choice of V (z) to include an additional multiplicative constant. We are interested in representations where j β 2 j is minimal -these will contribute to the leading-order asymptotics and so we refer to them as dominant.
Following [80], in order to write down the dominant asymptotics, it is helpful to introduce the notion of FH-representations. Given a symbol t(z) written in canonical form (52), replace all β j bỹ β j = β j + n j such that j n j = 0. This new function is the FH-representation t(z; n 1 , . . . , n m ), defined relative to the representation t(z; 0, . . . , 0) = t(z). We then have the equality: this means that, in general, the FH-representation differs from a canonical form for the symbol by a multiplicative constant. We illustrate this by example in Appendix D. An algorithm is given in [80] to find the finite number of dominant FH-representations, where it is shown that all of these contribute to the leading asymptotics of the determinant (44). For our purposes, finding a dominant representation will be simple; and given one dominant FH-representation of f (z) for which we define n j = 0, all other dominant FH-representations have n j ∈ {1, −1, 0}.
We now recall theorems relevant to the three cases introduced above. Szegő's strong limit theorem [81] gives the dominant asymptotics for matrices generated by smooth symbols with no winding, i.e. the case m = 0. We use a form adapted from reference [80]: (Szegő 1952). Let t(z) = exp(V (z)) be a symbol, with V (z) smooth as explained above and such that ∞ n=−∞ |n||V n | 2 < ∞. As the matrix dimension, N , goes to infinity: If we have a symbol with an integral winding number, i.e. m = 1, k 1 = 0, β 1 ∈ Z, the next theorem, adapted from a result of Fisher and Hartwig [82], allows us to reduce it to the product of a determinant that can be evaluated by Szegő's theorem and another small determinant.
Theorem 7 (Fisher, Hartwig 1969). Let t(z) = e V (z) z −ν , where V (z) satisfies the conditions for Theorem 6. Given b ± (z) as defined in (48), define the auxilliary functions: with associated Fourier coefficients 21 l k , m k . For ν > 0 we have: is otherwise unchanged. General estimates for the error terms δ ± k are given in [82] -the only case we need is as follows. Suppose that the large Fourier coefficients of h(z) = e V (z) behave as |h n | = O(ρ n ) and |h −n | = O(σ n ) then for large k, δ + k = O(ρ 2k σ k ) and δ − k = O(ρ k σ 2k ). Given the definitions in the above theorem, we can also state a formula from [82] for the leading order correction to Theorem 6.
Theorem 8 (Fisher, Hartwig 1969). Let t(z) = e V (z) satisfy the conditions for Theorem 6. Then we can write where, for l(z) and m(z) defined above, we have E is subdominant -see [82] for general estimates. For the case relevant to us, with ρ and σ defined in Theorem 7, we have E 2N σ 2N ). The final theorem we need is the generalised Fisher-Hartwig conjecture. The asymptotics for symbols with fractional jump discontinuities was initially conjectured by Fisher and Hartwig in [33]; this conjecture was then generalised to the class of symbols that we need by Basor and Tracy [83]. This generalised case was proved by Deift, Its and Krasovsky in [80], and we give a simplified form of their result relevant to our work.
The V n are unaltered when passing between FH-reps. Branches of b ± (z j )β j are determined by b ± (z j )β j = eβ j ∞ n=1 V ±n z ±n j . G(z) is the Barnes G-function [84, §5.17]; given as a Weierstrass product by where γ E is the Euler-Mascheroni constant. It is clear from equation (60) that G vanishes whenever the argument is a negative integer. Hence if β j ∈ Z, the RHS of (59) vanishes and this is not the first term of the asymptotic expansion (instead we should use Theorem 7).
6. Gapped chains -analysis 6.1. Closed form for the order parameter -Proof of Theorem 1b. Since c = 0, the complex function is given by , and it is only the sign of this real number that is important; moreover, since the Z j come in complex conjugate pairs, the sign only depends on N + Z , the number of zeros on the positive real axis and outside the unit circle. For bookkeeping purposes, define for a continuous logarithm that could be found by integrating the logarithmic derivative of f . We instead jump to the following solution: where the function Log(z) is the principal branch of the complex logarithm -this is clearly smooth and recovers f (z) when we take the exponential. Note that we used that the zeros are either real, or occur in complex conjugate pairs. On the unit circle z = e ik we can put z = 1/z into (66). This gives us an honest V (z) from which one can read off the Fourier coefficients: Inserting into Theorem 6 we reach: On expanding the square and interchanging the finite sums with the sum over n in the exponent, we can then perform the sum over n leading to Theorem 1b. The term under the fourth root is always a positive real, and the principal logarithm implies that we take the positive root. For completeness, note that the oscillatory factor multiplying the order parameter is given by e iπN (ω−1)+N log(s) .
Note that with the Fourier coefficients of V in hand, we can find the Wiener-Hopf decomposition (48) of our symbol when z is on the unit circle. .
where the square-root is continuous on the unit circle and the branch is fixed as the positive root of a positive real at z = 1.

Correlation lengths.
We now use Theorem 7 to find the behaviour of the correlation function O α (1)O α (N + 1) in the gapped phase ω. For definiteness, let us label the zeros by proximity to the unit circle: |Z i | ≤ |Z j | and |z i | ≥ |z j | for i < j.

Asymptotics of l(z), m(z).
The key ingredient that we need are the asymptotically large Fourier coefficients of the auxilliary functions Note that so far l and m are defined only on the unit circle and with the principal branch of the square-root (in fact, due to the complex-conjugate pairs of roots, the arguments of the square-root are strictly positive). For the purposes of this calculation, we assume the generic problem where the branch points in R = {z i , z −1 i , Z j , Z −1 j } are all distinct, we will comment later on the effect of multiplicity. We need the dominant asymptotic term of the nth Fourier coefficient of l(k) for large n: We analytically continue l(k) off the unit circle into the complex s-plane. The idea is to move the contour of integration out to infinity, where the s −n term in the integrand will cause the integral to vanish there. The integrand has branch cuts on which the contour gets snagged, and the dominant contribution will come from the nearest branch points outside the unit circle -this is the Darboux principle [85]. By inspection we have either a square-root or inverse square-root branch point at every element of R. If there are an odd number of such points inside (and therefore, by symmetry, outside) the unit circle, then zero and infinity are also branch points -hence there are always an even number of branch points both inside and outside the unit circle. We choose any branch cut pattern inside the unit circle (where no cut crosses the circle). Outside the unit circle we order the branch points by radial distance from the origin. In generic circumstances there will be either one real branch point (case A), or a complex-conjugate pair of branch points (case B), closest to the origin. Choose the cuts to be leaving all branch points radially. An example for each of the two cases is depicted in Figure 4 -we call the nearest branch point(s) s 1 (and s 1 ), for arg(s 1 ) ∈ [0, π]. We connect up the radial cuts outside a circle of large radius, the precise choice is unimportant.
In case A we consider the Hankel contour connecting infinity to the nearest real zero and backthis is exactly the relevant part of the snagged contour. After parameterising s = s 1 e t for t ∈ R + and where arg(t) = 0 below the axis and arg(t) = −2π above the axis, this integral obeys the conditions for Watson's lemma for loop integrals -see, for example, [86, §15.6.1] and [87]. This gives us an Figure 4. Schematic for the two generic cases of the computation (74). Blue (wavy) lines indicate branch cuts in the integrand. The black curve is the initial integration contour S 1 , and the red (lighter) curve is the deformed contour. × indicates the closest branch cuts to the unit circle that are outside the circle.
asymptotic series of which we need only the first term. Recall that we have ordered our zeros so that s 1 is either 1/z 1 or Z 1 -then we have Proposition 1. Suppose there is a single real root closest to the unit circle.
where the square-root is principal (with positive real argument).
This follows from the above discussion after using the same method to estimate the contribution of all other snagged contours -these are bounded above by |z * | −n where |z * | > s 1 , and are thus exponentially subdominant.
In case B we use the same method but now sum over the dominant contributions coming from the two branch points. This leads to Proposition 2. for This constant is the first term of the Taylor series of the regular part of the integrand at the branch point, and the square root is continuously connected to the principal branch on the real axis.
Note that in the case where we have only two roots, and they form a conjugate pair (as happens in the XY model), the constants are evaluated with the principal square root.
The exceptional cases where Propositions 1 and 2 do not apply are when f (z) has zeros with multiplicity, more than a pair of zeros closest to the unit circle, or both. We discuss these cases below.
We also need the asymptotic behaviour of m n . Fortunately no further analysis is needed: m(z) and l(z) share the same structure but are mutually inverse. Hence we have: Proposition 3. In the case of a nearest singularity s 1 on the real axis we have: For a complex conjugate pair of nearest singularities we have: where the c i are defined in (77). Now, if a zero has multiplicity two then we get either a simple pole of l(z) (and hence a zero of m(z)) or a zero of l(z) (and hence a simple pole of m(z)). A simple pole will give an exponential decay e −n/ξ , using Cauchy's theorem, with no algebraic prefactor (recall that ξ = 1/| log |ζ || where ζ is (any) one of the zeros of f (z) closest to the unit circle). A zero of l(z) is not a singularity so our contour will not be snagged there -we must hence look at the next-nearest singularity to the unit circle. Higher order multiplicities will give branch points, higher-order poles or higher-order zeros, and the calculations similarly go through. Higher-order poles will never have a vanishing residue for all n, and in fact for large n the dominant term in the residue will come from derivatives of s −(n+1) in (74). Importantly, even in these exceptional cases, the nearest zero always sets the longest correlation length for the operators O α . This is because, from the discussion above, either l n or m n has asymptotic decay controlled by the nearest zero (and hence there is an observable with correlation length ξ which follows from the calculation below).
Having more than two equidistant singularities requires summing over the contributions from each of them; this will give an e −n/ξ decay for zeros of multiplicity one (the coefficient must be calculated in each case, and for higher multiplicity one sums the contributions outlined above) -there may be destructive interference for certain values of n. This can include equidistant singularities coming from zeros both inside and outside the unit circle. Another exceptional case of this type is two closest zeros both on the real axis (i.e. at a and −a). Again we sum over the contributions which are given explicitly by the formulae in Propositions 1 and 3.
The final exceptional case is where we have degenerate closest zeroes which are mutually inverse. For example, if the closest zeros are at a and at 1/a ∈ R. This is the only case where ξ defined in terms of one of these closest zeroes is not realised as the longest correlation length (although it is still an upper bound) -the contribution of the mutually inverse zeros cancels in the definition of b ± (z) and so they do not contribute to the asymptotics of any O α . In such a case, the longest correlation length is set by the closest zero of f (z) whose inverse is not a zero. The starkest examples of this behaviour are in isotropic models, where b ± (z) = 1 and the correlation length is zero for all observables! This also follows from the observation that the ground state of a gapped isotropic model in our class is always a product state.
In summary, we have, in generic cases, that l n and m n decay exponentially with correlation length ξ. In exceptional cases their decay is at least this fast. Generically, if the nearest zero is inside the circle, we have an algebraically decaying prefactor n −3/2 for l n and n −1/2 for m n , this assignment is reversed if the nearest zero is outside. Moreover, if the nearest zero is complex then l n and m n have an oscillatory prefactor.

6.2.2.
Error terms in Theorem 7. In order to use Theorem 7 we need to estimate the errors δ ± N . To do so, we need to find ρ and σ such that for h(z) = e V (z) , |h n | = O(ρ n ) and |h −n | = O(σ n ) for large n.
Recall that the relevant e V (z) = e V 0 b + (z)b − (z) -this has exactly the same singularities as l(z) and m(z) up to exchanging square-roots with inverse square-roots. The analysis above goes through and we see that, in all non-degenerate cases, ρ = σ = 1/|s 1 |. We thus have that either d n = l n + O(l 2n m n ) or d n = m n + O(l n m 2n ). For n large this means that we can replace the matrix elements d n of the small determinant in Theorem 7 with either l n or m n without affecting the leading order behaviour. 6.2.3. The asymptotics of the correlator ( O α (1)O α (N + 1) -Proof of Theorem 2. Suppose that we are in the phase ω, then the generating function of the correlator is se V (z) z ω−α . In the case ω − α > 0, using Theorem 7 we have that The large determinant D N +ω−α (se V (z) ) is of Szegő form, and is, to leading order, equal to the result of Theorem 1b -i.e. the value of the order parameter. Inserting the dominant term of m N as found in the previous section, the second determinant may be evaluated directly to find the leading order term of the correlator. We have almost proved Theorem 2, but need to do some further analysis to isolate the exponential decay. This is the point where we specialise to generic situations, so that we are guaranteed that m N = Θ(e −N/ξ ). Then, in the position (i, j) of the second matrix we have a factor of e (−N +i−j)/ξ . The row and column index multiplicatively decouple, and so any individual term of the Laplace expansion of the determinant contains a factor of e −N (ω−α)/ξ , hence we may factor this out and we have that: .
The matrix elements of M (N ) are derived from the propositions above: i.e. K = 1/2 or 3/2 and α n are the coefficients that can oscillate with n. Hence, det M (N ) will contribute an algebraic dependence on N (and not affect the exponential scaling). For ω − α < 0 the same calculation goes through with m n replaced by l n (and the second matrix has dimension |ω − α|). We have hence proved Theorem 2. Now, putting together Theorems 1b and 2 prove Theorem 1a. In particular, we have shown that the correlators | O α (1)O α (N + 1) | do indeed form a set of order parameters that distinguish ω. The sign of f (z) (an invariant of our model) may be inferred by the presence or absence of (−1) N oscillation. As can be seen in (81), this oscillation depends on both ω and s, as defined in (63) (one must also take into account oscillations coming from det M ). N in closed form, but need that the first term in the sum (−l n+1 m n+1 ) gives the dominant scaling, as claimed in [33]. Thus, in the generic case, we have E (1) To see this, one needs to consider the different orders in the full asymptotic expansion of l n and m n , as given by Watson's lemma [86]. In particular, one can factor out the dominant term from |l n+1 m n+1 | and E (1) N then becomes a sum of many convergent series multiplied by non-positive powers of N (along with exponentially subdominant contributions coming from other singularities further from the circle than s 1 ). One of these convergent series is O(1) and we denote it by B N -this will oscillate with N if we have oscillation in l N and m N . Putting this together one reaches Theorem 3. The constant B N , along with further corrections, are evaluated in [32,48] for correlators that are equivalent to X and Y correlators in the XY model. 6.2.5. Possible alternative proof. An alternative approach to proving Theorem 2 would be to use the Fisher-Hartwig conjecture. The idea of such a proof is given in [88] -one should expand the Fourier contour defining the Toeplitz matrix (43) out to the nearest singularity in the generating function, and then rescale back to the unit circle. The deformed symbol is then singular on the unit circle (by construction), and, if it can be written in Fisher-Hartwig form, then Theorem 9 can be used to derive the leading order asymptotics. This method is applied in [50] to X and Y correlation functions in the XY model. 7. Gapless chains -analysis 7.1. Scaling dimensions. In this section we calculate the large N asymptotics of O α (1)O α (N + 1) for a system described by (9) with non-zero c. This was solved for isotropic models (i.e. models where f (z) = f (1/z)) and for α ∈ {−1, 0, 1} in [53]. We now explain how to use Theorem 9 to find the answer in the general case. The idea is simple; take the symbol corresponding to (−1) z −α f (z)/|f (z)| and find the dominant Fisher-Hartwig representations. This goes as follows: We reemphasise that in this analysis we pick the phase of the complex zeros such that k j ∈ [0, 2π). The smooth part e V (z) can be analysed as in the gapped case -in particular, the Fourier coefficients for n = 0 are given by (68) (the phase factor C is needed to put the singular part in canonical form, and this shifts V 0 ). Turning to the unit circle, z = e ik , the analysis of the singular part is split into three cases: a real zero at k = 0, π, a pair of complex conjugate zeros at k = φ and k = 2π − φ, or a set of zeros of multiplicity greater than one. The third (fine-tuned) case is discussed in Section 7.3, we ignore it for now. Note that we explicitly exclude such cases in the statement of Theorem 4 where we limit ourselves to chains described at low energy by a CFT. Now, for the real zero we have: = exp(ik/2) × cos(k/2)/| cos(k/2)| i sin(k/2)/| sin(k/2)| = i.
For a zero at −1 we have a sign-change discontinuity at k = π, and a zero at 1 has a sign-change-type 22 singularity at k = 0. For a complex conjugate pair of zeros at exp(±iφ) we have: Since φ = 0 or π, sign(cos(k) − cos(φ)) has sign-change discontinuities at k = φ and k = 2π − φ. We conclude that every zero contributes a factor exp(ik/2) as well as a sign-change at the location of the zero, which we can represent with a g k j ,β j (z) for β j any half-integer. Putting this information back into the symbol we reach where the β j are half-integer and 22 I.e. we are in the limiting case where g 1,β is actually constant, but the contribution to the Toeplitz determinant is as if it were a sign-change singularity.
We need to fix the multiplicative constant, C, by noting that the singular factors, isolated above, jumped between ±1, rather than ±i, and we also need to include e −iβ j k j in the singular part of (52). This leads to C = 2c j=1 e iβ j k j /g k j ,β j (1). To find the asymptotics we need some minimal representation (a set of β j minimising j β 2 j ), from which we can generate a set of minimal FH-reps to insert into Theorem 9. We find the solutions by first considering the cases c + ω − α = (2m − 1)c, for m ∈ Z, where the minimal solution is unique: for all j. If we have (2m − 1)c < c + ω − α < (2m + 1)c we form the set of minimal FH-reps by starting fromβ j = 2m−1 2 and sendingβ j →β j + 1 for of theβ j . We will consider our starting FH-rep (with all n i = 0) to be the one where we shift the first (1 − m)2c + ω − α of the β j . There are minimal FH-reps in total. Given the parameters ω, α, c we get that m = 1 + ω−α 2c , where x denotes the greatest integer less than or equal to x. Theorem 9 immediately gives the dominant scaling where This formula for ∆ α obscures some features of this function, to bring them out defineα = α − (ω + c), then one can show that where [x] is the nearest integer to x. It is then clear that ∆ α is symmetric underα ↔ −α and that the minimal scaling dimension of our operators is c/4.

7.2.
The dominant asymptotic term. We can go further with Theorem 9 to get the first term in the asymptotic expansion of O α (1)O α (N + 1) at large N . Firstly, note that: The first factor is a pure phase and contributes to the asymptotics as e N V 0 . By inspection 23 it is important to emphasise that this is an imaginary number, and again recall that k j ∈ [0, 2π). The second factor contributes in two ways: firstly through e n nVnV −n , exactly the quantity we calculated in Section 6.1. Secondly, we need powers of the Wiener-Hopf factors that were derived at the end of Section 6.1. Putting this all together we get: 23 The asymmetric second term in the sum is necessary because a singularity at kj = 0 is an edge case where where the sum is over all dominant FH-reps, these are parameterised by {n j } and defined in Section 7.1. ∆ α is ∆ α (c, ω) given in (92) and K is equal to −iV 0 + π(α − 1). The representation dependent O(1) multiplier is given by In our construction of the dominant FH-reps we start from setting allβ j to be equal and then add 1 to a fixed number of them. This means that the differenceβ i −β j is either 0 , 1 or −1 in all dominant FH-reps. Hence, pairs of complex conjugate zeros e ik i = e −ik j contribute e ink i N to the oscillatory factor in (95), where n ∈ {0, ±1}. We discuss the non-universal multiplier (96) in Appendix E. 7.3. Degenerate zeros on the unit circle. In the case that some of the zeros on the unit circle are degenerate, the analysis of the singular part given above follows through by raising to the power of the relevant multiplicity. Conjugate pairs of zeros must have the same multiplicity so contribute to the singular part as Equation (84) is similarly raised to the power m. We see an important difference between odd and even multiplicity. For odd m the degenerate zeros behave as above and we have a valid Fisher-Hartwig canonical form with β singularities at e ±iφ . In the case that m is odd for all zeros on the unit circle we can derive an analogue of Theorem 10, the steps are given in Appendix F. For even m at any zero, we do not have a Fisher-Hartwig canonical form and the analysis fails. The multicritical point in the XY model, with f (z) = (z + 1) 2 is an example of such a problematic case.

Extensions of our results
8.1. Long-range chains. In this section we discuss the effect of allowing our model (2) to have nonzero coupling constants between sites at 'long-range' -i.e. that there is no finite constant beyond which all couplings vanish. First, consider the case that the couplings decay with an exponential tail at large distances. This means that f (z) has a C ∞ smooth part, and a well defined winding number. We still have that ω = N z − N p , but note that poles are no longer restricted to the origin. The theory of Section 5, with care, may still be used to reach the same broad conclusions as in the finite-range case. In addition, we need the main result of [89]: Theorem 11 (Erhardt, Silbermann 1996). Take a symbol of the form i.e. a symbol with a single Fisher-Hartwig jump singularity at z = 1, and demand that exp(V (z)) is a C ∞ function. Then: where E is the constant defined implicitly in (59).
For nonvanishing E, or equivalently β ∈ Z, this is a special case of Theorem 9. However, for β ∈ Z this gives us a concrete asymptotic bound on the Toeplitz determinant in the case of a symbol with a C ∞ smooth part.
Szegő's theorem along with Theorem 11 allows us to extend the classification of gapped phases via string order parameters to long-range chains with C ∞ symbol. In particular, in the phase ω we have that O ω (1)O ω (N + 1) tends to a non-zero value that can be calculated using Szegő's theorem. Moreover, Theorem 11 proves that O α (1)O α (N + 1) for all α = ω tends to zero at large N . This proves that Theorem 1a remains valid for long-range chains with exponentially decaying couplings. In fact, one can go further than Theorem 11 and use the methods of [82] to give an analogue of Theorem 2. The αth correlator in the phase ω is O(e −N/ξα ), where ξ α is defined as above and ξ is derived from the singularity of the symbol closest to the circle (this singularity will come from either a zero or a pole of f (z)).
In critical chains, we may use the Fisher-Hartwig conjecture to derive the scaling dimensions exactly as in the finite-range case, on condition that there are finitely many zeros of f (z) that are on the unit circle, and that they remain well separated (this means that we may write our symbol in the canonical form (52)). Note that a study of the critical scaling of entanglement entropy for the isotropic subclass of such chains is included in [40]. While the winding number remains well defined, further analysis must be given to extend the results of [13] to long-range chains.
Models where couplings have algebraic tails are also physically relevant, and of topical interest [90,91]. In this case, f (z) will no longer be analytic and so singularities occur in the symbol distinct from Fermi points (zeros on the circle) and winding number (discontinuities in the logarithm). As f (z) is continuous, the winding number remains geometrically well-defined for gapped models. The theory of Toeplitz determinants may still be used in this case, and is deserving of a detailed analysis.
8.2. Uniform asymptotics approaching transitions. Our results give asymptotic correlations at particular points in the phase diagram. One may also be interested in how these correlations change along a path in parameter space, particularly where this path crosses a transition. This problem was studied analytically in reference [92] for the 2D classical Ising model (and hence the 1D quantum XY model). There are two cases where relatively recent 'black-box' results in the literature can be applied to a broader class of models. Firstly, consider a generalised Ising transition where we begin in a general gapped phase and a single zero approaches the unit circle. The relevant Toeplitz determinant asymptotics are given in reference [93]. Secondly, consider the case of two zeros e ±ik that come together. This is a generalisation of the approach to the multicritical point in the XY model along the isotropic critical line. The relevant Toeplitz determinant asymptotics are given in reference [94]. In both cases the crossover is controlled by a solution to the Painlevé V equation (althought a different one in each case). Due to the multiplicative nature of contributions to Toeplitz asymptotics, one would expect 24 similar behaviour in more general transitions where, as well as the approaching zeros, there are additional 'spectator' zeros on the unit circle.

Conclusion
Using Toeplitz determinant theory, we have investigated string-like correlation functions in a wide class of gapped and critical topological models. The salient features of their asymptotics can be deduced from the zeros of the associated complex function f (z). For example, the location of the zeros in the complex plane allow us to deduce whether the system is gapped or critical, furthermore giving access to correlation lengths and universal scaling dimensions (as summarised in Figure 2). Even detailed information, like the exact value of the order parameter, is a simple function of the zeros of f (z). The generality of these results allowed us to derive lattice-continuum correspondences, critical exponents and order parameters for the topologically distinct gapless phases. We now mention a few interesting paths to explore.
One surprising result was the universality of the ratios between the correlation lengths ξ α -this allowed for the extraction of the topological invariant ω. This was more striking for the dual spin chains, where local observables can be used to measure ω. It would be interesting to explore what happens upon introducing interactions. One possible scenario is that ratios of distinct correlation lengths give a measure of the interaction strength between the quasi-particles created by the corresponding operators.
One of the motivations of this work was to study how the invariants c and ω are reflected in physical correlations. The full classification of topological gapless phases within this symmetry class was obtained in the non-interacting case in reference [13]. Since this relied on concepts that are well-defined only in the absence of interactions, it does not directly generalise 25 . However, correlation functions and their symmetries are much more general concepts, and having now characterised the topology in terms of them, a natural next step is to use this to extend the classification to the interacting case.
Another interesting extension to the interacting case was touched upon in Section 3.1.2. This suggested that the Z classification of gapped phases should be stable under interactions if we allow 24 We are grateful to the participants of the AIM workshop 'Fisher-Harwig asymptotics, Szegő expansions and statistical physics' for discussions on this point. 25 However, numerical simulations indicated the stability away from the non-interacting limit.
only terms that contain the same number of real and imaginary Majorana modes. Such a structural property is not a conventional type of symmetry, and it would be interesting to investigate it, and its consequences, further. Lastly, as discussed in the previous section, the exact solvability and Toeplitz theory extend to cases with long-range couplings. This would certainly be interesting to explore, as removing constraints on f (z) leads to new asymptotic behaviours of the correlation functions beyond those that we have analysed in this paper. Table 4. Spin correlation functions that are the Jordan-Wigner dual of the fermionic string correlators O α (1)O α (N + 1) .
for some cases with degenerate zeros on the unit circle in Section 7.3. Then the dispersion relation |f (k)| cannot be linearised and we do not have a CFT description. It is clear that even conditioning on having many zeros on the unit circle, these are rare points in parameter space. Finally we mention the extra signs Σ. In the gapped case, the sign of f (1) is invariant -it must be real so can only change by passing through zero and hence closing the gap. A gapped model in the phase ω can be smoothly connected to f (z) = ±z ω . In reference [13] we showed that there are two invariant signs when c > 0 and the model is described by a CFT -in that case we can continuously connect any model to one with f (z) = ±z ω (z 2c+ω ± 1), the two signs cannot be removed without a phase transition. We hence have a description of the phase diagram that labels both gapped and critical phases by the triple (ω, c; Σ) where Σ ∈ Z 2 for c = 0 and Σ ∈ Z 2 × Z 2 for c > 0 gives the relevant signs. This sign information is easy to keep track of, so we classify phases including this sign -one is free to discard the extra information this gives.
A.3. The spin model. We now go into more detail related to Section 2.5. First a note on the Hilbert space of the spin chain. It is formally similar to that of the fermionic chains, as both are built from a set of two-dimensional Hilbert spaces. They differ in that the mathematical structure is simpler: M n=1 H n , where the local Hilbert space H n C 2 -in contrast to the fermions, operators localised on distinct sites commute.
We now define a Jordan-Wigner transformation that allows us to (almost) map H BDI into H spin and back. Let (iγ m γ m )γ n .
transform fermions into spins. The inverse transformation is given by Applying this transformation to H spin gives us (2), except that for all couplings extending over the final bond between sites L and L + 1 ≡ 1, we have a multiplicative factor of (−1) F -the total fermionic parity. Since the Hamiltonian (2) is quadratic, it preserves the parity, and so we can solve (2) in two total parity sectors. Details may be found in [39], where it is shown that we get two copies of (2), one with periodic and one with antiperiodic boundary conditions. Since we will be interested in bulk correlation functions, which will be independent of boundary conditions in the L → ∞ limit, we claim that simply using our results for the periodic fermion chain will be enough to understand these correlations in the periodic spin chain.
A.4. O α as a spin operator. In this section we explain how to derive the contents of Table 2, from which Table 4 follows. The quickest way to proceed in all cases is to use the nearest-neighbour substitutions : X n Y n+1 = γ n iγ n γ nγn+1 = −iγ nγn+1 (104) Y n X n+1 =γ n iγ n γ n γ n+1 = +iγ n γ n+1 .

Appendix E. Discussion of nonuniversal factors
It is interesting to note that the order parameter given in Theorem 1b is a symmetric function in the variables {z i } and, separately, {1/Z i } (listing zeros with multiplicity as distinct symbols). Another way to see why this occurs is through noting that a Toeplitz determinant generated by t(e ik ) is the same as the average of τ (e ik j ) := j t(e ik j ) over the group U(N ) (with eigenangles labelled by k j ) [40,96]. From the analysis of Section 6.1, for t(z) = e V (z) we have that τ is a symmetric function separately in the arguments {z j }, {1/Z j } and e ik j ; so can be expanded in a basis of symmetric functions. Let us write τ ({z j , Z j , e ik j }) = k a k s (1) k ({z j })s (2) k ({1/Z j })s k ({z j })s (2) k ({1/Z j }) can be factored out for each k, leaving us with a result that is a sum over products of symmetric functions and so the determinant is a symmetric function in the appropriate variables.
In the critical case we can rewrite the Θ(1) multiplier (96) in a way that gives a structure similar to the order parameter. In particular, when we have that all |β j | = 1/2, then Notice that, up to the normalising G-functions, this constant is built from terms of the form (1−ζ ±1 i ζ ±1 j ) where the ζ i are zeros of f (z). The sign ofβ j somehow tells us whether the jth zero on the unit circle acts as if it is inside the unit circle, or acts as if it is outside. Indeed, we see that if the ith and jth zero are both inside or both outside, we get a positive power of the term (1 − e ik j e −ik i ) and if one is in and one is out it is a negative power -this mirrors the behaviour of the factors coming from the z i and Z j . In the second line we have factors mixing zeros on the circle with zeros inside and outside the circle. Since all zeros on the circle come in complex conjugate pairs, terms of the form (1 − e ik l z j ) appear twice -however depending on the relative sign of β j and β j for this pair these factors can either cancel, or give a square. A similar squared term appears in the factor matching the z i to the Z j on the first line. For |β j | = n/2 for n > 1 we have a multiplicative effect where the contribution from the jth zero on the circle is counted n times. This is reminiscent of the CFT description where operators that involve many excitations give multiplicative contributions from the same Fermi point (which is located at some momentum k j ).
Appendix F. Critical models with degenerate roots on the unit circle As explained in Section 7.3, if we have any zeros of even multiplicity on the unit circle then we cannot proceed. Hence, consider f (z) with zeros of odd multiplicity m j at e ik j . The index runs over i = 1 . . . N 0 and so the total number of zeros on the circle is given by 2c = N 0 j=1 m i . Note that by symmetry we must have equal multiplicities at complex conjugate zeros.
The main difference in the analysis is that there is only one β for each unique zero (i.e. the degenerate zeros at that point correspond to only one FH singularity, but contribute to the winding multiple times) -this alters the sum rule (87), which becomes A method for solving (113) is to first solve (87) as in Section 7.1: assigning a half-integerβĵ wherê j = 1 . . . 2c. We then group these as: As all multiplicities are odd, this will lead to a canonical form for the symbol, but not necessarily a dominant FH-rep -this is because we minimised 2ĉ j=1β 2 j whereas we need to minimise N 0 j=1 β 2 j . We proceed by adding one to the smallest β j and subtracting one from the largest β j until the distance between smallest and largest is equal to zero or one. With this set of β we can construct a dominant canonical form as in Section 7.1, noting that the sum in the definition of V 0 (94) should now range over unique zeros only (i.e. goes from 1 to N 0 ). Moreover, if the β j are not all equal, we construct the other dominant FH-repsβ j = β j + n j as described in Section 5. We then have thatβ j ∈ { c+ω−α N 0 − 1/2, c+ω−α N 0 +1/2}, and the scaling dimension follows. Theorem 9 leads again to a variant of Theorem 10, where we sum over the dominant FH-reps just described, and where all products over the k j are over unique zeros only.