1 Introduction

Topological phases are fascinating examples of quantum matter. In one spatial dimension, they can be stabilised if the Hamiltonian has a symmetry group. Gapped topological phases have been classified for both non-interacting fermionic systems (dubbed topological insulators or superconductors) [1,2,3,4,5,6] as well as general fermionic and bosonic systems (dubbed symmetry protected topological (SPT) phases) [7,8,9,10,11,12]. However, it has recently been realised that critical matter can also form distinct topological phases—even without gapped degrees of freedom in the bulk [13]. As in the gapped case, the topology manifests itself physically: for example, through exponentially localised zero-energy modes at the physical edges. As long as a symmetry is preserved, a topological invariant can prevent two critical systems from being smoothly connected. Relatedly, there is a lot of recent interest in topological critical phases which do have additional gapped degrees of freedom [14,15,16,17,18,19,20,21,22,23,24,25,26,27].

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 [4]) to gapless topological phases [13]. These are labelled by a topological invariant \(\omega \) (\(\in \mathbb {Z}\)) and the central charge c (\(\in \frac{1}{2} \mathbb {Z}_{\ge 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 \(\omega \) 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 \(\mathbb {C} \setminus \{0\}\). Remarkably, c and \(\omega \) can then be obtained by counting zeros of f(z) (see Fig. 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 remained an open question 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 \(\omega \) are easily derived from its zeros. Moreover, our earlier work left an uneasy tension: distinct critical phases could be distinguished by the topological invariant \(\omega \), 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. This is relevant also to the spin chains that are Jordan-Wigner dual to these fermionic chains: whilst the non-interacting classification is less natural there, the correlation functions we obtain contain useful physical information that can be related to an interacting classification.

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, Toeplitz determinant theory, has a long association with statistical mechanics (for a review, see [29]), and our analysis generalises the pioneering work of [30,31,32] to a wider class of physical models.

Fig. 1
figure 1

The Hamiltonians we consider can be expanded in a basis \(H_\alpha \) (defined below Eq. (2)). The physics is encoded in the meromorphic function f(z). The given definitions of c and \(\omega \) classify the phases of \(H_{\mathrm{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 in the bulk, and \(\omega \) the topological properties

Fig. 2
figure 2

Universal asymptotics of the ground-state correlation functions considered in this work. If \(c=0\) (i.e. the system is gapped), then \(\mathcal {O}_\alpha \) has exponentially decaying correlations with correlation length \(\xi _\alpha \). The ratios \(\xi _\alpha /\xi _\beta \) are a universal function of \(\omega \), with the global scale set by \(1/\xi :=\min _{\zeta \in Z[f]} |\log |\zeta ||\); see the discussion before Eq. (9). (There is long-range order, i.e. \(a_\alpha \ne 0\), if and only if \(\alpha = \omega \); 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 \(\varDelta _\alpha \); see Theorem 4. The dependence of \(\varDelta _\alpha \) on both c and \(\omega \) means that these parameters may be determined by measurements of scaling dimensions; see the discussion below Theorem 4. Note that there are exceptional cases that behave differently, as discussed in the text

Since topological phases cannot be distinguished locally, in this work we study the correlations of so-called string-like objects \(\mathcal {O}_\alpha \) (labelled by \(\alpha \in \mathbb {Z}\)), meaning that \(\langle \mathcal {O}_\alpha (1) \mathcal {O}_\alpha (N) \rangle \) involves an extensive (\(\sim N\)) number of operators. Using Wick’s theorem, these correlations reduce to \(N\times N\) determinants. We calculate their asymptotic behaviour using the theory of Toeplitz determinants [29], 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 [33,34,35,36,37], and we indeed prove that \(\mathcal {O}_\alpha \) has long-range order if and only if\(\alpha = \omega \). More surprising is that the ratios of the correlation lengths of these operators are universal, i.e. they depend on \(\omega \) 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 \ne 0\)), all correlations are algebraically decaying and we obtain the corresponding scaling dimensions of \(\mathcal {O}_\alpha \). It turns out that measuring these gives access to both c and \(\omega \). 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.

Additionally, we discuss these correlators in the spin chain picture that is obtained after a Jordan-Wigner transformation. For odd \(\alpha \), \(\mathcal {O}_\alpha \) becomes a local spin operator, and long-range order in \(\langle \mathcal {O}_\alpha (1) \mathcal {O}_\alpha (N) \rangle \) signals spontaneous symmetry breaking. However, for even \(\alpha \) these correlators are string order parameters for spin chain SPTs. Whilst there is no natural notion of ‘non-interacting spin chains’, our analysis may be helpful for determining the (interacting) classification of topologically non-trivial critical spin chains.

Since the physical consequences of our results can be understood without going into the mathematical details, we structure the paper as follows. First, in Sect. 2, we outline the model and state our main results. In Sect. 2.5 we discuss the dual spin chain; then in Sect. 2.6 we discuss connections to previous works. In Sect. 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 Sects. 4 and 5. The proofs of our results then follow in Sects. 6 and 7 for the gapped and critical cases respectively. Finally, in Sect. 8, we explain how our results may be extended in different directions.

2 Statement of Main Results

2.1 The Model

Consider a periodic chain where each site has a single spinless fermionic degree of freedomFootnote 1\(\{c^\dagger _n, c_n; n = 1\dots L\}\). For convenience define the Majorana modes on each site:

$$\begin{aligned} \gamma _n = c^\dagger _n + c_n, \quad \tilde{\gamma }_n = \mathrm {i}\big (c^\dagger _n - c_n\big ), \end{aligned}$$
(1)

where \(\{\gamma _n,\gamma _m\} = 2 \delta _{nm}\) and \(\{\gamma _n, \tilde{\gamma }_m\} = 0\). Our class of interest—time-reversal symmetric, translation-invariant free fermions with finite-range couplings—has Hamiltonian [13, 38,39,40]

$$\begin{aligned} H_{\mathrm{BDI}} = \frac{\mathrm {i}}{2}\sum _{\alpha =-\alpha _l}^{\alpha _r} \sum _{n \in \mathrm {sites}} t_\alpha \tilde{\gamma }_n \gamma _{n+\alpha }, \quad t_\alpha \in \mathbb {R}. \end{aligned}$$
(2)

This can be understood as an expansion in the basis . The coupling between sites has maximum range \(\alpha _{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 \(\gamma _n\) (\(\tilde{\gamma }_n\)) are called real (imaginary) since \(T\gamma _nT = \gamma _n\) and \(T\tilde{\gamma }_nT =-\tilde{\gamma }_n\). This class of models is also invariant under parity symmetry \(P = \mathrm {e}^{\mathrm {i}\pi \sum _j c^\dagger _j c_j}\). We study the thermodynamic limit \(L \rightarrow \infty \), with \(\alpha _{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 Sect. 8.1. This model was first analysed in its spin chain form in reference [41].

The coupling constants \(t_\alpha \) establish a one-to-one correspondence between \(H_{\mathrm{BDI}}\) and the complex functions

$$\begin{aligned} f(z) = \sum _\alpha t_\alpha z^{\alpha }. \end{aligned}$$
(3)

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(\mathrm {e}^{\mathrm {i}k})|\) gives the one-particle energy of a mode with momentum k. The phase \(\arg (f(\mathrm {e}^{\mathrm {i}k}))\) is the angle required in the Boguliobov rotation that defines these quasiparticle modes [13, 42]. 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=\mathrm {e}^{\mathrm {i}k})\) whenever we restrict z to the unit circle.

2.2 Phase Diagram

Our results characterise the correlations in the different phases of \(H_{\mathrm{BDI}}\), hence we give them context by describing the phase diagram. First, phases of matter are defined as equivalence classes of ground states under smooth changes of the Hamiltonian, where two states are equivalent if they can be connected without a phase transition (i.e. a sharp change in physical behaviourFootnote 2). In particular, we define two critical models to be in the same phase if physical quantities such as scaling dimensions vary smoothly. Smooth changes to the Hamiltonian (i.e. smooth changes to \(t_\alpha \), including increasing the finite range by tuning some \(t_\alpha \) off zero) are equivalent to smooth motions of the zeros of f(z), as we discuss in Appendix A.2.

\(H_{\mathrm{BDI}}\) has two invariants that label both gapped and gapless phases (see also Fig. 1):

$$\begin{aligned} c&= \frac{1}{2}\left( \# \text { zeros of}~ f(z)~\text {on the unit circle}\right) \end{aligned}$$
(4)
$$\begin{aligned} \omega&= N_z - N_p, \end{aligned}$$
(5)

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-degenerateFootnote 3. Note that \(\omega \) 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,\omega )\) specifies the phases of \(H_{\mathrm{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 \(\varSigma \), that are discussed in Appendix A.2.

Note that the equivalence between \(H_{\mathrm{BDI}}\) and f(z) allows us to easily find a Hamiltonian within each phase: \(H_\omega \) is a representative of the gapped phase with winding number \(\omega \) and \(H_\omega +H_{2c+\omega }\) is a representative of the gapless phase \((c,\omega )\).

2.3 String Operators

The above is already established in the literature. The results of the current work show that given f(z), one can ‘read off’ detailed information about two-point ground-state correlation functions of the operators \(\mathcal {O}_\alpha (n)\):

$$\begin{aligned} \mathcal {O}_\alpha (n)={\left\{ \begin{array}{ll} \mathrm {i}^x \left( \prod _{m=1}^{n-1} \mathrm {i}\tilde{\gamma }_m\gamma _m\right) \gamma _{n}\gamma _{n+1}\dots \gamma _{n+\alpha -1} \quad &{} \alpha >0 \\ (-\mathrm {i})^{x}\left( \prod _{m=1}^{n-1} \mathrm {i}\tilde{\gamma }_m\gamma _m\right) (-\mathrm {i}\tilde{\gamma }_{n})\dots (-\mathrm {i}\tilde{\gamma }_{n+|\alpha |-1}) \quad &{}\alpha <0\\ \prod _{m=1}^{n-1} \mathrm {i}\tilde{\gamma }_m\gamma _m \quad &{}\alpha =0\end{array}\right. } \end{aligned}$$
(6)

where \(x= |\alpha |/2\) for \(\alpha \) even and \(x=(|\alpha |-1)/2\) for \(\alpha \) odd (the phase factors make \(\mathcal {O}_\alpha \) hermitian). These operators are a cluster of \(|\alpha |\) 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 Sect. 3.1, see also [43]. There are two typical behaviours for these correlators. Let angle brackets denote ground-state expectation value, then in the gapped case we expect:

$$\begin{aligned} \langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \sim A_\alpha + B_\alpha N^{-\delta _\alpha }\mathrm {e}^{-N/\xi _\alpha }. \end{aligned}$$
(7)

The constants \(A_\alpha \) and \(\delta _\alpha \) as well as the correlation length \(\xi _\alpha \) do not depend on N. If \(A_\alpha \ne 0\) we have long-range order. \(B_\alpha \) is \(\varTheta (1)\) and may include an oscillation with N. Note that, by translation invariance, we could equally well have considered the correlation function \(\langle \mathcal {O}_\alpha (r)\mathcal {O}_\alpha (N+r)\rangle \) for any \(r \in \mathbb {Z}\)—we fix \(r=1\) throughout for notational convenience. Note that \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\beta (N+1)\rangle = 0\) if \(\alpha \ne \beta \) as a simple consequence of the Majorana two-point functions given in Sect. 4.

We will see below that the ground state expectation value \(\lim _{N\rightarrow \infty }\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) in the gapped phase with winding number \(\omega \) is non-zero only when \(\alpha = \omega \), and is hence an order parameter for that phase—this can be seen as an extension of the results of [43].

Because these correlators contain a string of fermionic operators of length of order N, these are called string order parameters with value \(A_\alpha \). Note that in the case that \(\mathcal {O}_\alpha \) is local, (as happens in the spin picture given in Sect. 2.5), it is usual to call the one point function \(\langle \mathcal {O}_\alpha (n)\rangle \) the order parameter. This is because in that case the ground state will spontaneous collapse such that \(\sqrt{A_\alpha } = \langle \mathcal {O}_\alpha (n)\rangle \). 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:

$$\begin{aligned} \langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \sim C_\alpha N^{-2\varDelta _\alpha } \end{aligned}$$
(8)

where \(\varDelta _\alpha \) is the smallest scaling dimension of a CFT operator that appears in the expansion of the continuum limit of \(\mathcal {O}_\alpha \). The prefactor \(C_\alpha \) may include spatial oscillations, and further details are given in Sect. 3.4. Surprisingly, the set \(\mathcal {O}_\alpha \) 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:

$$\begin{aligned} f(z) = \frac{\rho }{z^{N_p}} \prod _{j=1}^{N_z} \left( z-z_j\right) \prod _{j=1}^{2c} \left( z-\mathrm {e}^{\mathrm {i}k_j}\right) \prod _{j=1}^{N_Z} \left( z-Z_j\right) . \end{aligned}$$
(9)

\(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 numberFootnote 4 of zeros inside, on and outside the circle are denoted \(N_z\), 2c, \(N_Z\) respectively, and \(\rho \) is a real number. Since the \(t_\alpha \) 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 \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) form a complete set of order parameters for the gapped phases of \(H_{\mathrm{BDI}}\).

Theorem 1a

In the phase \((\omega , c=0, \varSigma )\) we have

$$\begin{aligned} \lim _{N\rightarrow \infty } |\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle | = \mathrm {const} \times \delta _{\omega \alpha }. \end{aligned}$$
(10)

The non-zero constant is given in Theorem 1b. The value of the sign \(\varSigma \) may be inferred by the presence or absence of a \((-1)^N\) oscillation in this correlator.

Theorem 1b

In the phase \((\omega ,c=0,\varSigma )\), the non-universal value of the order parameter is given by

$$\begin{aligned} \lim _{N\rightarrow \infty }|\langle \mathcal {O}_\omega (1)\mathcal {O}_\omega (N+1)\rangle | = \left( \frac{ \prod _{i_1,i_2=1}^{N_z} (1-z_{i_1} z_{i_2}) \prod _{j_1,j_2=1}^{N_Z} \left( 1-\frac{1}{Z_{j_1} Z_{j_2}}\right) }{\prod _{i=1}^{N_z}\prod _{j=1}^{N_Z} \left( 1-\frac{z_i}{ Z_j}\right) ^2}\right) ^{1/4}. \end{aligned}$$
(11)

Thus from the decomposition (9) we can read off \(\omega =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 \(\xi = 1/|\log |\zeta _\star ||\) where \(\zeta _\star \) is any zero that maximises the right hand side of that equation (see Fig. 2 for illustration). The set of \(\zeta \) that are optimal in this way we call closest to the unit circle; we will always mean this logarithmic scaleFootnote 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 \(\omega \) we then have that \(\xi _\alpha \), as defined in (7), is equal to \(\frac{\xi }{|\omega -\alpha |}\) (for \(\alpha \ne \omega \))—this is a consequence of:

Theorem 2

If the system is in the phase \((\omega ,c=0,\varSigma )\) then, in generic systems, we have the large N asymptotics

$$\begin{aligned} \langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle = \det ( M (N))\left( \lim _{M\rightarrow \infty }|\langle \mathcal {O}_\omega (1)\mathcal {O}_\omega (M)\rangle |\right) \mathrm {e}^{-N|\omega -\alpha |/\xi }\mathrm {e}^{\mathrm {i}\pi N m}\bigl (1 +o(1)\bigr ); \end{aligned}$$
(12)

\(m\in \mathbb {Z}\) is a known constant and M(N) is a known \(|\omega -\alpha | \times |\omega -\alpha |\) matrix. The elements of this matrix have magnitudes that depend algebraically on N—in particular, \(\det M(N) = \varTheta (N^{-\delta })\) for some \(\delta >0\). Generic systems are those where the nearest zero(s) 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 \(\xi _\alpha \), but if one has multiplicity, \(\xi _\alpha \)may be controlled by the next-closest zero. The \(\xi _\alpha \) are always upper bounded by \(\xi \), and in fact this bound is saturated in all exceptional cases except when there are mutually inverse closest zeros. See the discussion in Sect. 6.2 for full details.

The form of \(\det M(N)\) derived in Sect. 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 \(|\omega - \alpha |\), as demonstrated in Table 1 for the case that \(|\zeta _\star |<1\). If \(|\zeta _\star |>1\) then the second and fourth columns of Table 1 should be interchanged (and the definitions of \(\lambda \) and \(\kappa \) change in the obvious way based on the formulae in Propositions 1 and 3).

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 \(\kappa \) and \(\lambda \) are independent of N and defined 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 (\zeta _\star ))\) (these oscillations may, however, not appear in the leading order term of \(\det M(N)\)). Moreover, if \(|\zeta _\star |<1\) then \(\det M(N) = \varTheta (N^{-K|\omega -\alpha |})\), where \(K=1/2\) for \(\omega -\alpha > 0\) and \(K=3/2\) for \(\omega -\alpha <0\). The assignment of K is reversed when \(|\zeta _\star |>1\).

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 \(\xi _\omega = \xi /2\), following from:

Theorem 3

In generic systems that are in the phase \((\omega , c=0,\varSigma )\), we have for large N that

$$\begin{aligned} \langle \mathcal {O}_\omega (1)\mathcal {O}_\omega (N+1)\rangle =\left( \mathrm {e}^{\mathrm {i}\pi N m}\lim _{M\rightarrow \infty }|\langle \mathcal {O}_\omega (1)\mathcal {O}_\omega (M)\rangle |\right) \left( 1+B_N\frac{ \mathrm {e}^{-2N/\xi }}{N^2}\right) \left( 1+o(1)\right) . \end{aligned}$$
(13)

The factor \(B_N\) is given implicitly in the proof and satisfies \(|B_N| = O(1)\), \(m\in \mathbb {Z}\) is a known constant. Generic systems are defined as in Theorem 2.

The results of the discussion in Sect. 6.2 allow extension to non-generic systems. Given non-zero correlation lengths of \(\mathcal {O}_\alpha \) for \(\alpha \ne \omega \), the formula \(1/\xi _\omega = 1/\xi _{\omega -1}+1/\xi _{\omega +1}\) holds. This agrees with Theorem 3 in the generic case where \( \xi _{\omega +1} =\xi _{\omega -1}\).

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 \(\mathbb {Z}_{\ge 0}\times \mathbb {Z}\): they are labelled by the central charge \(c\in \frac{1}{2} \mathbb {Z}_{\ge 0}\) and topological invariant \(\omega \in \mathbb {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 \(\mathcal {O}_\alpha \). A graphical representation of this theorem is given in Fig. 2.

Theorem 4

Consider a critical chain in the phase \((\omega ,c>0,\varSigma )\) where the 2c zeros on the unit circle are non-degenerate. Let \(\tilde{\alpha } = \alpha - (\omega + c)\). Then the operator \(\mathcal {O}_\alpha \) has scaling dimension

$$\begin{aligned} \varDelta _\alpha (c,\omega ) = c \left( \frac{1}{4} + x^2 - (x - [x])^2 \right) \Big |_{x = \tilde{\alpha } / 2c} \end{aligned}$$
(14)

[x] denotes the nearest integer to x.

Note that \(\varDelta _\alpha \) explicitly depends on the topological invariant \(\omega \). 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 Sect. 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 \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) 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 Sect. 3.4. A similar result holds when we have degenerate zeros on the unit circle, as long as every degeneracy is odd. We do 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 \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) and learning \(\omega \)—one simply looks for the value of \(\alpha \) 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 at a transition between critical phases (discussed in detail below). In the critical case one can determine c and \(\omega \), but this requires information about more than one correlator. Inspecting the form of Eq. (14), displayed in Fig. 2, one concludes that it is not necessary to measure the scaling dimension of all \(\mathcal {O}_\alpha \) in order to determine the phase. One method would be to measure the scaling dimensions of \(\{\varDelta _\alpha ,\varDelta _{\alpha \pm 1},\varDelta _{\alpha +2}\}\) for some convenient \(\alpha \), and form the set of \(\delta _\alpha = \varDelta _{\alpha +1}-\varDelta _\alpha \). This difference is equal to \([(\alpha - \omega )/2c-1/2]\)—this means that \(\delta _\alpha \) is a constant integer on plateaus of widthFootnote 6 2c, and that neighbouring plateaus differ in value by one. If the \(\delta _\alpha \) are all differentFootnote 7 then we must have \(c=1/2\) and \(\omega \) can be determined easily using \(\omega = \alpha -\delta _\alpha \). Otherwise, one should then measure further scaling dimensions until the width of the constant plateau (equal to 2c) is found. Once c is known, \(\omega \) may be determined: on the edge of the plateau we have \(\omega = \alpha -2c \delta _\alpha \). 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 \(\mathcal {O}_\alpha \) are local for \(\alpha \) odd. In Appendix C we show that it is possible to recover both c and \(\omega \) 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 \(\omega \) from knowing only two correlation lengths; this is explained in Sect. 3.1.

2.5 The Dual Spin Chain

Our results apply not only to \(H_{\mathrm{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

$$\begin{aligned} X= \left( \begin{array}{cc}0~ &{}\quad 1 \\ 1~ &{}\quad 0\end{array}\right) ,\quad Y= \left( \begin{array}{cc}0 ~&{}\quad -\mathrm {i}\\ \mathrm {i}~ &{}\quad 0\end{array}\right) , \quad Z= \left( \begin{array}{cc}1 &{}\quad 0 \\ 0 &{}\quad -1\end{array}\right) . \end{aligned}$$
(15)

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

$$\begin{aligned} H_{\mathrm {spin}}&= \frac{t_0}{2} \sum _{n\in \mathrm {sites}} Z_n-\sum _{\alpha > 0} \frac{t_\alpha }{2} \left( \sum _{n\in \mathrm {sites}}X_n \left( \prod _{m=n+1}^{n+\alpha -1} Z_m\right) X_{n+\alpha }\right) -\nonumber \\ {}&\sum _{\alpha < 0} \frac{t_\alpha }{2}\left( \sum _{n\in \mathrm {sites}}Y_n \left( \prod _{m=n+1}^{n+|\alpha |-1} Z_m\right) Y_{n+|\alpha |}\right) . \end{aligned}$$
(16)

As before, we only allow a finite sum over \(\alpha \), have \(t_\alpha \in \mathbb {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_{\mathrm{spin}}\) to \(H_{\mathrm{BDI}}\). The main point is that our results for the behaviour of \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) apply equally well to the spin chain. The expressions for \(\mathcal {O}_\alpha \) in terms of spin operators are given in Table 2. Some of these operators appeared in the recent works [44, 45]. Note that, as displayed in Table 2, \(\mathcal {O}_\alpha \) is local for odd \(\alpha \) but remains a non-local string for even \(\alpha \). One can easily see that for odd \(\alpha \), \(\langle \mathcal {O}_\alpha \rangle \) is zero in any symmetric state; hence Theorem 1a implies that for odd \(\omega \) we have spontaneous symmetry breaking. For even \(\omega \), however, \(\mathcal {O}_\omega \) is a string order parameter for the spin chain. As we discuss further in Sect. 3.1, this is indicative of the spin chain forming an interacting SPT phase. Note that the two phases can coexist: for \(\alpha =3\), the order parameter is \(\mathcal {O}_3(n) = X_{n+1}Y_{n+2}X_{n+3}\); hence \(P = \prod _j Z_j\) and \(T=K\) have been broken. However, PT is preserved, and due to the similarity of \(\mathcal {O}_3\) with the cluster model Hamiltonian—\(\sum _j X_jZ_{j+1}X_{j+2}\)—the ground state is also an SPT phase protected by PT. Further details may be found in reference [40].

Table 2 Spin operators that are the Jordan-Wigner dual of the fermionic operators \(\mathcal {O}_\alpha \)

2.6 Relation to Previous Work

In reference [46], Lieb, Schultz and Mattis set the stage for the analysis of determinantal correlations in free fermion models and related spin chains. The key reference related to our results for gapped models is the classic paper of Barouch and McCoy [31]. There the authors study bulk correlations in 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 \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) for \(\alpha = 1,-1\) in the phases \(\omega =0,1\) that include what one would obtain from our theorems. Beyond that, the paper [47] includes a calculation of the value of the order parameter for \(\alpha = -1, 2\) in the special case that \(f(z) = z^3 - \lambda \). Some portion of the phase diagram for \(-2\le \omega \le 2\) is mapped out in reference [48] where order parameters are identified and calculated numerically. Several papers, for example [49, 50], study spin models with competing ‘large’ cluster term and Ising term (i.e. non-zero \(t_\alpha \), \(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 \(|\omega |\le 4\) in the spin language. Equivalent fermionic order parameters are easily found using the Jordan-Wigner transformation and the paper [43] includes the fermionic \(\mathcal {O}_{\alpha }\) for \(|\alpha |\le 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 [31] should be mentioned, along with [51], as seminal early works that derived critical behaviour for correlators \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) with \(|\alpha | \le 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 [52], and the same correlators were analysed in reference [53] using the mathematical methods found below. We also mention the quantum inverse scattering method as a tool for calculating scaling dimensions in certain cases [54] .

An isotropic spin chain is invariant under spin-rotation around the Z axis. In our fermionic model \(H_{\mathrm{BDI}}\), this manifests as invariance under spatial inversion \(H_\alpha \leftrightarrow H_{-\alpha }\), and hence is a model for which \(f(z) = f(1/z)\). This relation implies that \(\omega =-c\). The correlators \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) with \(|\alpha | \le 1\) in isotropic models with general c and \(\omega = -c\) were derived in references [55, 56] using the same methods as this paper. Our results go further by studying a wider class of models, including critical phases with general \((c,\omega )\), as well as a wider class of observables: \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) for all \(\alpha \). This allows us to observe that from knowledge of the scaling dimensions of these operators, one can identify the topological invariant \(\omega \).

3 Physical Context and Discussion

In this section we interpret our results and give them context. In Sect. 3.1 we discuss universality and its implications for spin chains, as well as the relation to symmetry fractionalisation. In Sect. 3.2 we relate correlation length in the bulk to localisation length at the edge and in Sect. 3.3 we derive critical exponents from our results. In Sect. 3.4 we connect our lattice results to continuum (CFT) models and finally in Sect. 3.5 we discuss entanglement scaling approaching a transition.

3.1 SPT Phases and String Orders: Universality and Symmetry Fractionalisation

In Sect. 2.4, we saw that for both gapped and critical systems in our class of models, one can measure the topological invariant \(\omega \) by looking at the string order parameters \(\mathcal {O}_\alpha \). For gapped phases, one needs to find the value of \(\alpha \) 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 non-interacting 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 Sect. 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, \(\mathcal {O}_\alpha \) has correlation length \(\xi _\alpha = \frac{\xi }{|\alpha -\omega |}\) (if \(\alpha \ne \omega \)), see Fig. 2 (we assume the generic case for this discussion). This means that although \(\xi \) depends on microscopic properties (like the position of the zeros of f(z)), the ratio \({\xi _\alpha }/{\xi _\beta }\) depends only on \(\omega \) 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 \(\alpha \) such that \( |\langle \mathcal {O}_\alpha (1) \mathcal {O}_\alpha (N)\rangle |\) tends to a non-zero limit as \(N \rightarrow \infty \). This requires going through an arbitrarily large set of observables. Surprisingly, it is sufficient to measure only, for example, two correlation lengths \(\xi _{\alpha _1}\) and \(\xi _{\alpha _2}\) (for the observables \(\mathcal {O}_{\alpha _1}\) and \(\mathcal {O}_{\alpha _2}\)) for any fixed choice of \(\alpha _1\) and \(\alpha _2\) satisfying \(|\alpha _1 - \alpha _2| \in \{1,2 \}\). To see this, note that there are three cases. Firstly, if one finds long-range order for either \(\alpha _1\) or \(\alpha _2\), then \(\omega \) is known. Secondly, if \(\xi _{\alpha _1} = \xi _{\alpha _2}\), one knows that \(\omega = \frac{\alpha _1 + \alpha _2}{2}\). In any other case, \(\alpha _1\) and \(\alpha _2\) will either both be larger or smaller than \(\omega \), such that \(\frac{\xi _{\alpha _1}}{\xi _{\alpha _2}} = \frac{\omega - \alpha _2}{\omega - \alpha _1}\). This can be uniquely solved, giving \(\omega = \frac{\alpha _1 \xi _{\alpha _1} - \alpha _2 \xi _{\alpha _2}}{\xi _{\alpha _1} - \xi _{\alpha _2}}\).

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 \(\alpha _1 = 1\) and \(\alpha _2 = -1\), then this corresponds to the correlation lengths \(\xi _X\) and \(\xi _Y\) of the local observables \(X_n\) and \(Y_n\). This fully determines the invariant

$$\begin{aligned} \omega = \left\{ \begin{array}{ccl} -1 &{} &{}\quad \text {if }\lim _{N \rightarrow \infty }\langle Y_1 Y_N \rangle \ne 0\\ 1 &{} &{}\quad \text {if }\lim _{N \rightarrow \infty }\langle X_1 X_N \rangle \ne 0\\ 0 &{} &{}\quad \text {otherwise and if } \xi _X = \xi _Y\\ \frac{\xi _X+\xi _Y}{\xi _X-\xi _Y} &{} &{}\quad \text {otherwise.} \end{array} \right. \end{aligned}$$
(17)

This means that one can distinguish, for example, the trivial paramagnetic phase from the topological cluster phaseFootnote 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.

3.1.2 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 [35, 57].

Fig. 3
figure 3

Consider a system which does not exhibit spontaneous symmetry breaking. a The ground state remains invariant when applying a global symmetry \(U = \prod _n U_n\). b When acting with \(U^X\) on a line segment X of length \(l \gg \xi \), then deep within that segment, the action is indistinguishable from the global symmetry operation U. Hence, the state can only be changed near the boundaries of X. In conclusion, effectively, \(U \approx U^L U^R\)

Let U be some on-site unitary symmetry, i.e. \([U,H] = 0\) with \(U = \prod _n U_n\). Consider the operator \(U^X = \large \prod _{n\in X} U_n\) where X is some large line segment of length l (see Fig. 3). If \(l \gg \xi \), 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 \(|\text {gs}\rangle \) is the ground state, then effectively \(U^X |\text {gs}\rangle = U^L U^R |\text {gs}\rangle \), 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 [35, 58]. This phenomenon is known as symmetry fractionalisation and is the essential insight that led to the classification of (interacting) SPT phases in 1D [9, 11, 12, 59]. To illustrate this, consider the case with a symmetry group \(\mathbb {Z}_2 \times \mathbb {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 \(U^L V = (-1)^{\frac{\omega }{2}} V U^L\) where \(\omega \in \{0,2\}\). 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 \(\mathbb {Z}_2 \times \mathbb {Z}_2\).) In fact, one can show that in the phase where \(\omega = 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, referring the interested reader to the review in reference [40]. Instead, we consider the effect on correlation functions.

One can consider the string correlation function \(\left\langle \text {gs} \right| U^X \left| \text {gs} \right\rangle = \langle U^X \rangle = \langle U^L U^R \rangle \). Due to locality, we have that \(\langle U^X \rangle \approx \langle U^L \rangle \langle U^R \rangle \) for \(l \gg \xi \). Since SPT phases do not spontaneously break symmetries, we have that \(\langle U^L \rangle = \langle U^L V \rangle = (-1)^\frac{\omega }{2} \langle V U^L \rangle = (-1)^\frac{\omega }{2} \langle U^L \rangle \). Hence, we conclude that the string correlation function \(\langle U^X \rangle \) has to be zero if \(\omega = 2\). Equivalently, measuring \(\langle U^X \rangle \ne 0\) implies that \(\omega = 0\), such that one calls this a string order parameter for the trivial phase (\(\omega = 0\)). Analogously, one can construct a string order parameter for the non-trivial phase: if \(\mathcal {V},\mathcal {W}\) are local operators anticommuting with V, then by repeating the above argument, we conclude that \(\langle \mathcal {V} U^X \mathcal {W} \rangle \ne 0\) implies that \(\omega = 2\) (with \(\mathcal {V}\) (\(\mathcal {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 = \large \prod _n Z_n\) and \(V = P_\text {odd} = \large \prod _{m\text { odd}} Z_m\). Two models with this symmetry are \(H_0 = \sum _n Z_n\) and \(H_2 = - \sum _n X_{n-1} Z_n X_{n+1}\). One can calculate that their symmetry fractionalisations are \(\omega =0\) and \(\omega = 2\), respectively. The above tells us that \(\prod _{n=1}^N Z_n\) is a string order parameter for the trivial phase. Similarly, taking \(\mathcal {V}_{n,n+1} = Y_n X_{n+1}\)—which indeed anticommutes with \(P_\text {odd}\)—then \(\mathcal {V}_{1,2} \left( \prod _{n=1}^N Z_n \right) \mathcal {V}_{N+1,N+2}\), or equivalently \(X_1 Y_2 \left( \prod _{n=3}^N Z_n \right) Y_{N+1} X_{N+2}\), is an order parameter for the topological phase \(\omega = 2\). In this case, the string order parameters we have derived—with respect to \(\mathbb {Z}_2\times \mathbb {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 [60]. 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. \(\gamma _n \gamma _{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 . 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 \(-\omega \) (and \(P_R\) charge \(\omega \)). 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 \(\langle \mathcal {V} \rangle = 0\) for any operator whose charge is non-zero. Hence, in this way we are naturally led to consider \(\gamma _1 \cdots \gamma _\alpha \left( \tilde{\gamma }_{\alpha +1} \gamma _{\alpha +1} \cdots \tilde{\gamma }_{N} \gamma _{N} \right) \tilde{\gamma }_{N+1} \cdots \tilde{\gamma }_{N+\alpha }\): this operator can only have long-range order if \(\alpha = \omega \).

The power of symmetry fractionalisation is that it is not specific to the non-interacting case. Does that mean that if one considers interacting Hamiltonians preserving the aforementioned structural symmetryFootnote 9, that we obtain an infinite number of distinct interacting phases, labelled by \(\mathbb {Z}\) and distinguished by the same string order parameters as the non-interacting case? This is not clear: the structural symmetry is not a conventional symmetry; in particular, the charge of this symmetry is not well-behaved under the multiplication of operators (which is usually essential to preserve charge under renormalization group flow). Hence, it is perhaps unlikely that this structural symmetry gives rise to non-trivial physics in the interacting case, but it might nevertheless be interesting to explore further.

3.2 Majorana Modes: Localisation Length Versus Correlation Length

In reference [13], we showed that if \(\omega >0\), then the system has \(\omega \) Majorana zero modes per edge. More precisely, to each of the \(\omega \) largest zeros \(\{z_i\}_{i=1,\cdots ,\omega }\) of f(z) within the unit disk, we associate a hermitian operator \(\gamma _L^i\), 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 \(\omega \) mutually anticommuting symmetries, from which one can show that to each edge we can associate a \(\sqrt{2}^\omega \)-fold degeneracy.Footnote 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\ne 0\)).

Reference [39] 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 \(\omega > 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 between 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 \(\omega >0\) holds also for \(\omega < -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 left-right inversion corresponds to \(f(z) \leftrightarrow f(1/z)\), which at the level the topological invariant corresponds to \(\omega \leftrightarrow -\omega -2c\). For any other value of \(\omega \), 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 \(\tau = \frac{T-T_c}{T_c}\). In the current quantum setting, there is an equally naturalFootnote 11 tuning parameter \(\varepsilon \in \mathbb {R}\): if f(z) represents a gapless Hamiltonian, then \(f_\varepsilon (z) := f(z(1-\varepsilon ))\) interpolates between two gapped phases. One can think of this as shrinking (\(\varepsilon <0\)) or expanding (\(\varepsilon >0\)) the radial component of the zeros of f(z). We will work in the case where the system is gapless at \(\varepsilon =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 / mdistinct 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, \(\eta \), defined by the scaling of the order parameter, \(\varPsi \), at the critical point (\(\langle \varPsi (1)\varPsi (N) \rangle \sim 1/N^\eta \)); \(\nu \), which encodes the divergence of the correlation length (\(\xi \sim |\varepsilon |^{-\nu }\)); the dynamical critical exponent, z, that relates how the correlation length \(\xi \) diverges relative to the characteristic time scale defined by \(\tau \), the inverse energy gap (\(\tau \sim \xi ^z\)); and \(\beta \) which relates to the decay of the order parameter (\(\varPsi \sim |\varepsilon |^\beta \)).

Exactly at the critical point, we have given explicit results above only for \(m=1\). In particular, Theorem  4 gives us that \(\eta = {c}/{2}\). In the special case of \(c={1}/{2}\), we recover the well-known result \(\eta = 1/{4}\). For \(m>1\) and odd, we can use Eq. (113) from Appendix F to easily derive that \(\eta = mc/2\). The other three exponents are defined away from criticality where we can allow for any \(m\ge 1\). The correlation length is determined by the nearest zero, such that \(\xi \sim 1/|\ln |1+\varepsilon || \sim 1/|\varepsilon |\). Hence, \(\nu = 1\); independent of c and m. The energy gap, \(\min _{z\in 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+\varepsilon ))^m\). The gap scales as \(\sim |\varepsilon |^m\), such that \(\tau \sim \xi ^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+\varepsilon )^{-2}|^{m^2/8} \sim |\varepsilon |^{m^2/8}\); all other terms do not go to zero with \(\varepsilon \). Combining 2c / m such factors, we have \(\varPsi \sim |\varepsilon |^{\frac{m^2}{8} \times \frac{2c}{m}}\). Hence, \(\beta = mc/4\). In the special case of the Ising transition, \(m=1\) and \(c={1}/{2}\), this reduces to the well-known result \(\beta = {1}/{8}\).

The above results are consistent with the well-known scaling relations [61]. In particular, we can straightforwardly confirm that \(2\beta = \nu (\eta +d-2)\), where our dimension is \(d=2\). In fact, the aforementioned relationship implies that \(\eta = {mc}/{2}\) for any multiplicity m. Moreover, such scaling relations can be used to derive other critical exponents: such as \(\gamma = \nu (2-\eta ) = 2-mc/{2}\), \(\alpha = 2-\nu d = 0\) and \(\delta = {\nu d}/{\beta } -1 = {8}/{mc}-1\). It is interesting to note that the critical exponent \(\gamma \) changes sign for \(mc = 4\). This data is summarised in Table 3.

Table 3 Summary of the critical exponents found in Sect. 3.3

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 \(\mathcal {O}_\alpha \).

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 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 [52, 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 \(\pi \)). The central charge is always equal to c.

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 \(\mathcal {A}_{\varvec{n}}(x)\) that create fixed numbers, \(n_j\), of quasiparticles of type j (i.e. excitations near momentum \(\pm k_j\)). For equal times, these take the form

$$\begin{aligned} \big \langle \mathcal {A}^\dagger _{\varvec{n}}(0)\mathcal {A}_{\varvec{n}}(x)\big \rangle = \sum _{\{m_j\}} C_{\mathcal {A}_{\varvec{n}}}x^{-2\sum _j \varDelta _{\varvec{n}}^{(j)}} \mathrm {e}^{ \mathrm {i}x \sum _j m_j k_j} \end{aligned}$$
(18)

for scaling dimensions \(\varDelta _{\varvec{n}}^{(j)}\), Fermi momenta \(k_j\) and the sum is over sets of \(m_j \in 2\mathbb {Z}\) . The amplitude \(C_{\mathcal {A}_n}\) depends on the appropriate form factor [65]. The oscillatory term comes from a multiplicative factor \(\mathrm {e}^{-\mathrm {i}x \delta p}\) when \(\mathcal {A}\) gives an intermediate excitation with momentum \(\delta 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 [54, 65, 66]. The \(\mathcal {O}_\alpha \) that we consider are ‘square-roots’ of local operators in the sense that the product \(\mathcal {O}_{\alpha }(n)\mathcal {O}_{\alpha }(n+m)\) (for small m) is a local fermionic operator on the lattice and has an expansion dominated by such \(\mathcal {A}\). One may then expect that correlations of \(\mathcal {O}_{\alpha }\) will have the same form as (18), but with \(m_j \in \mathbb {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 Sect. 2.5.

If our model has one zero on the unit circle (which must be at \(z=\pm 1\), corresponding to \(k=0, \pi \) respectively), then our continuum limit has \(c=1/2\), and is hence described by the Ising CFT [68]. The operator content of this theory is well understood and, amongst the primary operators, there are exactly two with scaling dimension 1 / 8. These are the ‘spin’ and ‘disorder’ operators, denoted by \(\sigma \) and \(\mu \) respectively. In the usual lattice Ising model, \(\sigma \) is the field corresponding to the local order parameter of the neighbouring ordered phase and \(\mu \) is a non-local string order parameter of the neighbouring paramagnetic phase. Moreover, \(\sigma \) is odd under parity symmetry (realised on the lattice as \(P = \prod _j Z_j\)), while \(\mu \) is even under P. Now, if we are in a model with \(\omega =0\) and \(c =1/2\), then we can continuously tune to a critical Ising chain \(H = \sum _j X_jX_{j+1} \pm Z_j\) for which the correspondence is standard and discussed in [68]. In particular, we will have \(X_n \rightarrow \sigma (n)\) and \(\prod _{j=-\infty }^n Z_j \rightarrow \mu (n)\); note that this is consistent with the aforementioned symmetry considerations. Moreover, models in our class with \(\omega =-1\) and \(c =1/2\) are in the same phase as a critical Ising chain \(H = \sum _j Y_jY_{j+1} \pm Z_j\) where the analogous standard correspondence is \(Y_n \rightarrow \sigma (n)\) and \(\prod _{j=-\infty }^n Z_j \rightarrow \mu (n)\), again consistent with locality and symmetry properties.

Our results on the lattice, in particular Theorem 4, allow us to extend this by showing that a system with \(c=1/2\) and winding number \(\omega \) has two operators with this dominant scaling dimension: \(\mathcal {O}_\omega \) and \(\mathcal {O}_{\omega +1}\). By considering locality and parity symmetry of the operators, as well as scaling dimension, we identify \(\mathcal {O}_\omega \) with \(\sigma \) when \(\omega \) is odd and \(\mu \) when \(\omega \) is even; \(\mathcal {O}_{\omega +1}\) is then the other field. Importantly, we conclude that the operators on the lattice that have overlap with the dominant primary scaling fields depend on \(\omega \). Indeed, along a path that connects a model with \(c=1/2\) and winding number \(\omega \) to a model with \(c=1/2\) and winding number \(\omega ' \ne \omega \) then we must encounter a multicritical point with \(c\ge 1\) where the pair of dominant operators in the \(c=1/2\) models will have degenerate scaling dimension and the dominant pair changes as we go through this point—we give an example of this in Sect. 3.4.4. Other operators \(\mathcal {O}_\alpha \) for \(\alpha \) neither \(\omega \) nor \(\omega +1\) will be dominated by descendants of \(\sigma \) for \(\alpha \) odd and of \(\mu \) for \(\alpha \) even (the lattice operators have dimension \(1/8 + j\) for some \(j\in \mathbb {Z}_+\), so we should take CFT descendants at level j). Note that in all correspondences we may have a \((-1)^n\) oscillatory factor.

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 \(\mathrm {e}^{\pm \mathrm {i}k_F}\), and a U(1) symmetry generated by (the standard model in this class is the XX spin chain [52], which in our language corresponds to \(H = H_{-1} + H_1\) with \(f(z) = z+1/z\)). These are isotropic models with \(f(z) = f(1/z)\), which implies that \(\omega =-c\); we discuss other values of \(\omega \) later.

The fermionic Hamiltonian may be bosonised as described in [63, 69, 70] and [52, Chap. 20], passing also to a continuum limit. We denote the resulting bosonic fields by \(\theta (x)\) and \(\varphi (x)\), such that

$$\begin{aligned}{}[\partial _x\varphi (x),\theta (y)]=[\partial _x\theta (x),\varphi (y)] = \mathrm {i}\pi \delta (x-y) \end{aligned}$$
(19)

where \(\delta (x)\) is the Dirac delta function. We now recall some standard results for the isotropic model  [52, 71]. Firstly, the vertex operator \(\mathrm {e}^{\mathrm {i}\theta (x)}\) creates a localised charge of the aforementioned U(1) symmetry, and \(\partial _x \varphi (x)\) is a density fluctuation: the total density is \(\rho (x) = (k_F + \partial _x \varphi (x))/\pi \). Vertex operators of the form \(\mathrm {e}^{\mathrm {i}(\lambda _1\theta (x) + \lambda _2 \varphi (x))}\) have scaling dimension \((\lambda _1^2+ \lambda _2^2)/4\). When \(\lambda _1 \in \mathbb {Z}\) and \(\lambda _2 = 2k\) for \(k \in \mathbb {Z}\), these operators are well-defined and mutually local. When \(\lambda _1 \in \mathbb {Z}\) and \(\lambda _2 = 2k +1\) for \(k\in \mathbb {Z}\), these operators have a branch cut extending to infinity. As discussed below, on the lattice this branch cut is related to the non-local Jordan Wigner strings. Notice that this is an asymmetry in the role of the fields \(\varphi \) and \(\theta \).

We now consider operator correspondences for \(\mathcal {O}_\alpha \). 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 in order to correctly correspond to the lattice model and hence to our results.

Firstly, following [63, 69], we note that \(\rho (x)\) generates the U(1) symmetry. Hence, we can make the formal identification

$$\begin{aligned} \mathrm {e}^{-\mathrm {i}\frac{\chi }{2} \sum _{j} Z_j} = \mathrm {e}^{\mathrm {i}{\chi } \sum _{j}(c^\dagger _jc_j-1/2)} \rightarrow \mathrm {e}^{\mathrm {i}{\chi } \int (\rho (x) -1/2) \mathrm {d}x}. \end{aligned}$$
(20)

The Jordan-Wigner strings follow by setting \(\chi = \pi \) and truncating the sum and integral respectively. We can then make the correspondence (see, for example, [69, 70] and [52, Chap. 20]):

$$\begin{aligned} \prod _{i=-\infty }^n Z_i \rightarrow a\mathrm {e}^{\mathrm {i}(k_F n + \varphi (n))}+ \overline{a}\mathrm {e}^{-\mathrm {i}(k_F n + \varphi (n))} \qquad (a \in \mathbb {C}). \end{aligned}$$
(21)

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

$$\begin{aligned} \psi _\pm (x) = \mathrm {e}^{\mathrm {i}\theta (x) \pm \mathrm {i}\varphi (x)} + \dots \quad \left( c_n \rightarrow \psi (n) \simeq \mathrm {e}^{\mathrm {i}k_F n}\psi _+(n) + \mathrm {e}^{-\mathrm {i}k_F n} \psi _-(n) \right) . \end{aligned}$$
(22)

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 \(\varphi \rightarrow \varphi \) and \(\theta \rightarrow -\theta \) under T. We can also confirm that the right-hand-side of (21) is hermitian and does not transformFootnote 12 under the U(1) or T, as required for the Z-string. The site at \(-\infty \) 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:

$$\begin{aligned} \bigg \langle \prod _{i=0}^n Z_i \bigg \rangle \rightarrow \bigg \langle \left( \mathrm {e}^{\mathrm {i}(k_F n + \varphi (n))}+ \mathrm {e}^{-\mathrm {i}(k_F n + \varphi (n))}\right) \left( \mathrm {e}^{\mathrm {i}\varphi (0)}+ \mathrm {e}^{-\mathrm {i}\varphi (0)}\right) \bigg \rangle , \end{aligned}$$
(23)

where we now suppress constant coefficients.

Now, by considering the U(1) action and time-reversal symmetry, we see that \(\sigma ^{\pm }_n \rightarrow \mathrm {e}^{\pm \mathrm {i}\theta (n)}\) (with real coefficients). We then have that

$$\begin{aligned} X_n\rightarrow & {} \cos (\theta (n)) \nonumber \\ Y_n\rightarrow & {} \sin (\theta (n)). \end{aligned}$$
(24)

Note that the relation between \(\sigma ^+\) and \(\psi _{\pm }\) in the CFT corresponds on the lattice to (103).

The above correspondences are well known [52, 69, 70]. In our notation, (21) corresponds to \(\mathcal {O}_0\) and (24) correspond to \(\mathcal {O}_{\pm 1}\). Operators \(\mathcal {O}_\alpha \) for \(|\alpha |>1\) are more involved. However, the correspondences that we have already established should allow us to find the continuum operator by taking operator products. We consider first the family of local spin operators, with \(\alpha \) odd. In particular, let us analyse \(\mathcal {O}_3(n) = X_nY_{n+1}X_{n+2}\). From Theorem 4, we know that this operator has scaling dimension 9 / 4. In the field theory the operators \(\mathrm {e}^{\pm 3\mathrm {i}\varphi }\) and \(\mathrm {e}^{\pm 3\mathrm {i}\theta }\) share this scaling dimension, as well as operators that include derivatives, for example \((\partial ^2 \theta ) \mathrm {e}^{\mathrm {i}\theta }\). We should exclude terms that include \(\mathrm {e}^{\mathrm {i}(2k+1) \varphi }\) on the grounds of locality. We now show that all operators depending solely on \(\theta \) and that have the correct scaling dimension may appear. Moreover, we give a method for determining the exact correspondence.

First, note that any product of three neighbouring X or Y operators can be expanded as a linear combination of terms \(\sigma _{n}^{\pm }\sigma ^{\pm }_{n+1}\sigma ^{\pm }_{n+2}\), with all possible sign combinations present. These products transform separately under U(1) with charge \(m\in \{-3,-1,1,3\}\) equal to the sum of the signs. The dominant terms would be \(\mathrm {e}^{\pm \mathrm {i}\theta (n)}\), with subdominant contributions from \(\mathrm {e}^{\pm 3 \mathrm {i}\theta (n)}\) and products of \(\mathrm {e}^{\pm \mathrm {i}\theta (n)}\) and \(\mathrm {e}^{\pm 3 \mathrm {i}\theta (n)}\) with derivatives of \(\theta (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 as stated above \(\mathcal {O}_{\pm 3}(n)\) has the same scaling behaviour as \(\mathrm {e}^{\pm 3 \mathrm {i}\theta (n)}\). Details of a formal calculation in terms of the CFT operator product expansion are givenFootnote 13 in Appendix B. Intuitively, \(\mathcal {O}_{\pm 3}\) is dominated by terms that create three charges, as well as the remnants of several terms that create one charge but (partially) destructively interfere with each other.

By generalising this idea to all odd \(\alpha \), we conjecture that \(\mathcal {O}_\alpha (n)\) has an expansion of the form

$$\begin{aligned} \mathcal {O}_\alpha (n) \rightarrow \sum _{m=0}^{|\alpha |} \mathrm {e}^{\mathrm {i}(|\alpha |- 2m)\theta (n)}\mathcal {D}_{(\alpha ,m)}(\theta (n))+\dots \quad \quad (\alpha ~ \mathrm {odd}). \end{aligned}$$
(25)

\(\mathcal {D}_{\alpha ,m}(\theta )\) is constant for \(m = \pm \alpha \), and for other values of m contains products of derivatives of \(\theta \) such that the scaling dimensions of each term matches the extremal terms \(\mathrm {e}^{\pm \mathrm {i}\alpha \theta (n)}\). This is consistent with both Theorem 10 and calculations of the type in Appendix B. Note that this conjecture includes all operators in the field theory that depend only on \(\theta \) and that have the correct scaling dimension; carefully taking operator products will give explicit expressions for the \(\mathcal {D}\).

The case of \(\alpha \) even follows the same pattern, except that we should always include a Jordan-Wigner string: the dominant term of which is given in (21). Hence, we arrive at

$$\begin{aligned}&\mathcal {O}_\alpha (n) \rightarrow \sum _{m=0}^{|\alpha |} \left( \mathrm {e}^{\mathrm {i}(|\alpha |- 2m)\theta (n)+\mathrm {i}(k_F n + \varphi (n))}+\mathrm {e}^{\mathrm {i}(|\alpha |- 2m)\theta (n)-\mathrm {i}(k_F n + \varphi (n))}\right) \mathcal {D}_{(\alpha ,m)}(\theta (n))\nonumber \\&\qquad \qquad \qquad \quad +\dots \quad (\alpha ~\mathrm {even}), \end{aligned}$$
(26)

where, as above, \(\mathcal {D}_{\alpha ,m}(\theta )\) is constant for \(m=\pm \alpha \) 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 cases 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 [40], generalised Kramers-Wannier dualities were discussed that map between our models. One class of transformations swap models such that \(f(z) \leftrightarrow z^m f(z)\) for some \(m \in \mathbb {Z}\). If f(z) is isotropic, then \(z^{m}f(z)\) is not—this allows us to extend the preceding correspondence to anisotropic models. Anisotropic models that are dual to isotropic models in this way have an appropriately transformed U(1) symmetry. First, let us separate two cases: m even and m odd. When m is even, the transformation will be local, and when m is odd it will be non-local (this is related to whether the neighbouring gapped phases have a local order parameter). Now, in both cases we should take the correspondences (25) and (26) and shift the label on the left-hand-side \(\mathcal {O}_\alpha \rightarrow \mathcal {O}_{\alpha +m}\) while not shifting the right-hand-sides. In the case that m is even, we take this as the correspondence. In the case that m is odd, we alter the right-hand-side by swapping \(\varphi \) and \(\theta \). For example, in the transition \(H=-\sum _n X_nZ_{n+1}X_{n+2} +\sum _n Z_n\) we identify \(X_n\) with \(\mathrm {e}^{\mathrm {i}(k_F n + \theta (n))}+ \mathrm {e}^{-\mathrm {i}(k_F n + \theta (n))}\). Notice that for odd m, the oscillatory factor \(\mathrm {e}^{\mathrm {i}k_F n}\) appears with the field \(\theta (x)\). This is because \(\partial _x \theta (x)\) is related to the density of the transformed U(1) symmetry. Moreover, this correspondence is consistent with the requirement that the non-local vertex operators with factors \(\mathrm {e}^{(2k+1)\varphi (x)}\) always correspond to non-local lattice operators.

The dualities discussed so far allow us to map an isotropic model to a representative of each phase \((c=1,\omega )\). As mentioned above, these representatives all have a U(1) symmetry and are thus not generic. 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 persistFootnote 14 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

$$\begin{aligned} f(z)&= (1-\lambda )z^2 + 2\lambda z - (1+\lambda )\nonumber \\&= (z-1)\bigl ((1-\lambda )z + (1+\lambda )\bigr ). \end{aligned}$$
(27)

Tuning \(-1\le \lambda \le 1\) we interpolate between two critical phases, with a transition at \(\lambda =0\). For \(\lambda =1\) the system is the standard critical Ising chain with \(H = -(\sum _j X_jX_{j+1} +Z_j)\). Models with \(\lambda >0\) are in the same phase: \((c=1/2,\omega =0)\). For \(\lambda =-1\) we have \(H = -(\sum _jX_jZ_{j+1}X_{j+2} - X_jX_{j+1})\). Models with \(\lambda <0\) can be smoothly tuned to this model and have \((c=1/2,\omega =1)\). For \(\lambda = 0\), we have the \(c=1\) transition between topologically distinct critical phases, with \(H = -\frac{1}{2}(\sum _jX_jZ_{j+1}X_{j+2}+Z_j)\). Table 4 gives the behaviour of the scaling dimensions of the most dominant \(\mathcal {O}_\alpha \) 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].

Table 4 Behaviour of scaling dimensions across a transition between critical phases—the model is \(f(z)= (1-\lambda )z^2 + 2\lambda z - (1+\lambda )\) and the dominant CFT fields associated to the dominant lattice operators are also given

3.4.5 CFT Operator Correspondence with c Complex Fermions

For higher values of c, we work formally with the linearised theory that consists of 2c real fermions. Note that it has been shown that spin models with \(f(z) = \pm z^\omega (z^{2c}\pm 1)\) are described at low energy by \(so(2c)_1\) WZW models  [48, 49], although a lattice-continuum operator correspondence has not been made. We can smoothly connectFootnote 15 any critical model in \(H_{\mathrm{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 bosonised as described above to give a set of fields \(\theta _j(x)\) and \(\varphi _j(x)\). The relevant vertex operators are of the form

$$\begin{aligned} \tau _{\varvec{\mu },\varvec{\nu }} (x)= \prod _{j=1}^c \mathrm {e}^{\mathrm {i}\left( (\mu _j + \nu _j) \theta _{j}(x) +(\mu _j - \nu _j) \varphi _{j}(x) \right) } \end{aligned}$$
(28)

where \(\mu _i\) and \(\nu _i\) are half-integer and have scaling dimension \(\varDelta _{\varvec{\mu },\varvec{\nu }} =\sum _j(\mu _j^2 +\nu _j^2)/2\) (we suppress Klein factors [64]). The half-integer condition makes them twist operators for the linearised fermion field [72]. They act nontrivially on all fermionic sectors, and by considering decoupled lattice models with Hamiltonian \(H_{-m}+H_{m}\) (where the bosonisation in each sector is clear), their locality properties correctly reflect those of the operators \(\mathcal {O}_\alpha \). Moreover, they have a minimal scaling dimension of c / 4 (notice that this coincides with the smallest scaling dimension of the \(\mathcal {O}_\alpha \), given in Theorem 4). As in the \(c=1\) case, when we have a U(1) symmetry, then \(\partial _x\varphi _j(x)\) or \(\partial _x \theta _j(x)\) will correspond to fluctuations in the charge density. Hence, the appropriate field will be accompanied by \(k_j x\)—this is \(\varphi (x)\) (\(\theta (x)\)) when the U(1) charge is local (nonlocal) on the lattice. Again, as above, these oscillatory factors persist away from these symmetric models.

Suppose now that we are in an isotropic model, with U(1) symmetry generated by \(S^z_{\mathrm{tot}}\). The conjectured expansion of lattice operators \(\mathcal {O}_\alpha \) goes through roughly as above. However, the identification of, say, \(\sigma ^+\) with a charge one operator is no longer so restrictive; this is because charges can have different signs in the different sectors and cancel. By considering the scaling dimensions, one concludes that the leading order terms have charge distributed evenly throughout the different sectors. For example, observe that the charge-two operator \(\mathrm {e}^{\mathrm {i}(\theta _1(x) + \theta _2(x))}\) with \(\varDelta = 1/2\) dominates the charge-two operator \(\mathrm {e}^{\mathrm {i}(3\theta _1(x) -\theta _2(x) )}\) with \(\varDelta =5/2\).

More generally, as argued in [56], in isotropic models \(\sigma ^+=\left( \mathcal {O}_1 + \mathrm {i}\mathcal {O}_{-1}\right) /2\) will be dominated by operators \(\tau _{\varvec{\mu },\varvec{\nu }}\) with \(\sum _j (\mu _j+\nu _j) = 1\) (charge condition) and \(|\mu _i -\nu _j |\le 1\) (dominance condition). These conditions give a sum of terms that are products of \(\mathrm {e}^{\pm \mathrm {i}\theta _j(x)}\) or \(\mathrm {e}^{\pm \mathrm {i}(\varphi _j(x)+k_j x)}\) in each sector, and hence that can be distinguished by the presence or absence of oscillatory factors \(\mathrm {e}^{\pm \mathrm {i}k_j x}\). This is analogous to the sum over Fisher-Hartwig representations (see Sect. 5) that we derive in Sect. 7.1, and the relevant oscillatory factors are confirmed in Theorem 10 (indeed, each such term is represented in the final result). Further, in [56] it was argued that \(\mathcal {O}_0\) will have \(\sum _j( \mu _j+\nu _j) = 0\) (charge condition) and \(|\mu _i -\nu _j |\le 1\). To extend this to operators with \(|\alpha |> 1\), consider first the operators \(\tau _{\varvec{\mu },\varvec{\nu }}\) with \(\sum _j(\mu _j +\nu _j) =|\alpha |\) (a maximal charge condition) and \(|\mu _i -\nu _j |\le 1\) (dominance condition). This gives a set of operators that we expect to dominate the continuum limit of \(\mathcal {O}_\alpha \). However, as in the \(c=1\) case, we expect that we should include terms where maximally charged operators \(\mathrm {e}^{\mathrm {i}R \theta (x)}\) are substituted with a series of terms \(\mathrm {e}^{\mathrm {i}(R-2m) \theta (x)}|_{m=1,\dots ,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.

To extend this to non-isotropic models, with \(\omega = -c+m\), we note that there will again be two cases—m even and m odd. For m even we simply shift the correspondences argued for the isotropic case. For m odd we shift and further swap \(\varphi _j(x)\) and \(\theta _j(x)\). These correspondences may also be derived from \(\sum _j(\mu _j +\nu _j) = |c+\omega -\alpha |\) (a maximal charge condition, when this is meaningful) and \(|\mu _i -\nu _j |\le 1\) (dominance condition)—although when m is odd we should swap \(\nu _j \rightarrow - \nu _j\) in these equations. Then, as mentioned above, we should include descendant operators with lower charge. These correspondences again agree with Theorem 10.

Example: Consider a model with \(c=\omega =2\) (for example, \(H= H_6+H_{2}\)). By solving \(\sum _{j=1}^2(\mu _j+\nu _j) = 3\) and then including descendants, we conjecture the correspondence:

$$\begin{aligned} \mathcal {O}_{1}(x)&\rightarrow \biggl (\mathrm {e}^{\mathrm {i}(2 \theta _1(x)+\varphi _1(x) + k_1 x)}+\mathrm {e}^{\mathrm {i}(2 \theta _1(x)-\varphi _1(x) -k_1 x)}+ \mathcal {D}(\theta _1)\left( \mathrm {e}^{\mathrm {i}(\varphi _1(x) + k_1 x)}+\mathrm {e}^{-\mathrm {i}(\varphi _1(x) +k_1 x)}\right) \nonumber \\&\quad \quad \,\, + \;\mathrm {e}^{-\mathrm {i}(2 \theta _1(x)+\varphi _1(x) + k_1 x)}+\mathrm {e}^{-\mathrm {i}(2 \theta _1(x)-\varphi _1(x) - k_1 x)}\biggr )\mathrm {e}^{\mathrm {i}\theta _2(x)}+ \left( 1\leftrightarrow 2\right) + \dots . \end{aligned}$$
(29)

Note that all terms have \(\varDelta = 3/2\), and all oscillate as either \(\mathrm {e}^{\pm \mathrm {i}k_1 x}\) or \(\mathrm {e}^{\pm \mathrm {i}k_2 x}\) as expected. Furthermore, our conjecture gives the same expansion for \(\mathcal {O}_{7}\), although the (suppressed) coefficients are not expected to be the same in general. Indeed, Theorem 10 indicates that the coefficients are different in the general case, since (95) is not symmetric under taking sign-reversed Fisher-Hartwig representations.

In the case that \(c\ge 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 [73]. 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 \(\rho _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 [42, 74]. 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 [75].

Theorem 5

(Its, Mezzadri, Mo 2008) Consider a sequence of gapped chains (as defined in Eq. (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 \(\zeta _j\), noting that \(\zeta _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 \rightarrow \infty \), has the following expansion as \(|\zeta _j| \rightarrow 1\):

$$\begin{aligned} S[\rho _A]= -\frac{1}{6} \sum _{j=1}^{2c} \log |\zeta _j - 1/\overline{\zeta }_j |+ O(1) . \end{aligned}$$
(30)

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: \(\{\mathrm {e}^{\pm \mathrm {i}\phi _1} \mathrm {e}^{1/\xi },\mathrm {e}^{\pm \mathrm {i}\phi _2} \mathrm {e}^{t_2/\xi ^{r_2}},\dots \mathrm {e}^{\pm \mathrm {i}\phi _c} \mathrm {e}^{t_c/\xi ^{r_c}} \}\), where \(\phi _i \ne \phi _j\) for \(i \ne j\), \(t_j>1\) and \(r_j <1\), and we approach the circle letting \(\xi \rightarrow \infty \). The conditions on \(t_j\) and \(r_j\) ensure that a closest zero is \(\mathrm {e}^{ \mathrm {i}\phi _1} \mathrm {e}^{1/\xi }\) for \(\xi \) large enough. Inserting into Theorem 5, we get that:

$$\begin{aligned} S[\rho _A] = \sum _{j=1}^c \frac{r_j}{3} \log (\xi ) +O(1). \end{aligned}$$
(31)

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 \(H= \sum _i X_iX_{i+1} - Y_iY_{i+1}\) which 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 [76] that applies far more generally and gives asymptotics as the lattice spacingFootnote 16\(a \rightarrow 0\)

$$\begin{aligned} S[\rho _A] = \frac{c}{3} \log (\xi '/a) +O(1), \end{aligned}$$
(32)

where c is the central charge of the underlying CFT and \(\xi '\) is the (fixed) correlation length of the system under consideration. This may also be interpreted as a scaling limit \(\xi ' \gg a\) [77], and the formula was confirmed in this sense for the XY model in [78]. Further relevant references are found in the review articles [77, 79]. We see that Eq. (30) is equivalent to formula (32) 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 renormalization group flow lines in parameter space.

4 String Correlators as Determinants

We now begin the analysis necessary to prove the results given in Sect. 2.4.

4.1 Fermionic Two Point Correlators

After defining f(z) as in Eq. (3), we have that:

$$\begin{aligned} H = \sum _k |f(k)|\eta ^\dagger _k \eta _k + \mathrm {const} \end{aligned}$$
(33)

where the Boguliobov quasiparticles are found by rotating the Bloch sphere vectorFootnote 17\((c_{-k},c^\dagger _k)\) through an angle f(k) / |f(k)| about the x-axis, giving

$$\begin{aligned} \eta _k = \frac{1}{2}\left( 1 +\frac{ {f(k)}}{|f(k)|} \right) c^\dagger _{k}+ \frac{1}{2}\left( 1 -\frac{ {f(k)}}{|f(k)|} \right) c_{-k}. \end{aligned}$$
(34)

The sum over k goes over momenta \(k_n = 2\pi n/L\), although we always work in the limit where this sum becomes an integral from 0 to \(2\pi \). Details of this diagonalisation may be found in, for example, [13, 41, 42]. The ground state, \(\left| \text {gs} \right\rangle \), is the vacuum for the quasiparticles \(\eta _k\), and from this we can easily calculate fermionic correlation functions—we refer the reader to [42] for details. We will use

$$\begin{aligned} \langle -\mathrm {i}\tilde{\gamma }_n \gamma _m \rangle&:= \left\langle \text {gs} \right| \left( -\mathrm {i}\tilde{\gamma }_n \gamma _m\right) \left| \text {gs} \right\rangle = \frac{1}{2\pi } \int _0^{2\pi } \frac{ {f(k)}}{|f(k)|} \mathrm {e}^{-\mathrm {i}(m-n) k}\mathrm {d}k \nonumber \\ \langle \tilde{\gamma }_n\tilde{\gamma }_m \rangle&= \langle \gamma _n\gamma _m \rangle = \delta _{nm} \end{aligned}$$
(35)

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 Sect. 6.2 allows us to find the large N asymptotics of \(\langle -\mathrm {i}\tilde{\gamma }_n \gamma _{n+N} \rangle \). As explained in the discussion around Sect. 6.2.4, in generic cases this correlator will be \(\varTheta (N^{-K} \mathrm {e}^{-N/\xi })\) where \(K\in \{1/2,3/2\}\) is easily determined. For critical chains f(k) has jump discontinuities. Decomposing as in Eq. (85) and integrating by parts, we have that the fermionic two point function is \(\varTheta (1/N)\)—this behaviour is as expected from the CFT description.

4.2 Wick’s Theorem

Because the Hamiltonian (2) is quadratic, ground state expectation values have a Pfaffian structure. More precisely, suppose that we have 2N distinct and mutually anticommuting operators, \(\mathcal {A}_n\), then:

$$\begin{aligned}&\langle \mathcal {A}_1\cdots \mathcal {A}_{2N}\rangle = \sum _{\mathrm{all ~pairings}} (-1)^\sigma \prod _{\mathrm{all ~pairs}~(m,n)} \langle \mathcal {A}_m\mathcal {A}_n \rangle \nonumber \\&\quad \bigg (=\langle \mathcal {A}_1\mathcal {A}_2 \rangle \langle \mathcal {A}_3\mathcal {A}_4 \rangle \cdots \langle \mathcal {A}_{2N-1}\mathcal {A}_{2N} \rangle -\langle \mathcal {A}_1\mathcal {A}_3 \rangle \langle \mathcal {A}_2\mathcal {A}_4 \rangle \cdots \langle \mathcal {A}_{2N-1}\mathcal {A}_{2N} \rangle +\dots \bigg ).\nonumber \\ \end{aligned}$$
(36)

\((-1)^\sigma \) 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 \(\langle \mathcal {A}_m\mathcal {A}_n \rangle \), and is a form of Wick’s theorem that is given in reference [46].

4.3 String Correlation Functions

Consider the two point correlation function of \(\mathcal {O}_\alpha \) for \(\alpha >0\):

$$\begin{aligned}&\displaystyle \langle \mathcal {O}_\alpha (1) \mathcal {O}_\alpha (N+1) \rangle&= (-1)^x \langle \gamma _1\dots \gamma _\alpha \left( \prod _{n=1}^{N} \mathrm {i}\tilde{\gamma }_n\gamma _n\right) \gamma _{N+1}\gamma _{N+2}\dots \gamma _{N+\alpha }\rangle \end{aligned}$$
(37)
$$\begin{aligned}&\displaystyle = (-1)^x\bigg \langle (-\mathrm {i}\tilde{\gamma }_1)(-\mathrm {i}\tilde{\gamma }_2)\dots (-\mathrm {i}\tilde{\gamma }_\alpha ) \left( \prod _{n=\alpha +1}^{N} \gamma _n(-\mathrm {i}\tilde{\gamma }_n)\right) \gamma _{N+1}\gamma _{N+2}\dots \gamma _{N+\alpha }\bigg \rangle \end{aligned}$$
(38)

We now transpose further terms to put unlike Majoranas as nearest neighbours and apply Wick’s theorem:

$$\begin{aligned} \langle \mathcal {O}_\alpha (1) \mathcal {O}_\alpha (N{+}1)\rangle= & {} (-1)^{N(\alpha -1)} \bigg \langle \prod _{n{=}1}^N(-\mathrm {i}\tilde{\gamma }_n \gamma _{n+\alpha })\bigg \rangle =(-1)^{N(\alpha -1)}\det (\langle -\mathrm {i}\tilde{\gamma }_n \gamma _{m+\alpha }\rangle )_{m,n=1}^{N} \end{aligned}$$
(39)
$$\begin{aligned}= & {} (-1)^{N(\alpha -1)}\det \left( \frac{1}{2\pi } \int _0^{2\pi } \frac{ {f(k)}}{|f(k)|}\mathrm {e}^{-\mathrm {i}\alpha k} \mathrm {e}^{-\mathrm {i}(m-n) k}\mathrm {d}k\right) _{m,n=1}^{N}. \end{aligned}$$
(40)

For \(\alpha <0\) an analogous calculation again leads to Eq. (40). Table 2 gives the spin operators Jordan-Wigner dual to the fermionic operators \(\mathcal {O}_\alpha (n)\) and Table 5 gives the equivalent spin correlators for all \(\alpha \). A derivation is given in Appendix A.4. Notice that for odd \(\alpha \) these operators and correlation functions are local in the spin variables, and for even \(\alpha \) they are nonlocal; they are always nonlocal for the fermions. Understanding the asymptotic behaviour of the determinant (40) is the key to the results given in Sect. 2.

5 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\times N\) Toeplitz matrix, T, takes the following ‘translation-invariant’ form:

$$\begin{aligned} (T)_{mn} = (t_{m-n}) = \left( \begin{array}{ccccc}t_0 &{}\quad t_{-1} &{}\quad t_{-2} &{}\quad \ldots &{}\quad t_{-(N-1)} \\ t_{1} &{}\quad t_0 &{}\quad t_{-1} &{}\quad \ldots &{}\quad t_{-(N-2)} \\ t_2 &{}\quad t_{1} &{}\quad t_0 &{}\quad \ldots &{}\quad t_{-(N-3)} \\ \vdots &{}\quad \vdots &{}\quad \vdots &{}\quad \ddots &{}\quad \vdots \\ t_{N-1} &{}\quad t_{N-2} &{}\quad t_{N-3} &{}\quad \ldots &{}\quad t_0 \end{array}\right) . \end{aligned}$$
(41)

This matrix can be thought of as the \(N \times 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\subseteq U \subseteq \mathbb {C}\). A function \(t: U \rightarrow \mathbb {C}\), integrable on the unit circle, generates a Toeplitz matrix through its Fourier coefficients:

$$\begin{aligned} t_{n} = \frac{1}{2\pi } \int ^{2 \pi }_{0} t\left( \mathrm {e}^{\mathrm {i}k}\right) \mathrm {e}^{- \mathrm {i}n k} \mathrm {d}k. \end{aligned}$$
(42)

We refer to such t(z) as the symbol of the Toeplitz matrixFootnote 18 and denote the Toeplitz determinant of order N that is generated by t as

$$\begin{aligned} D_{N}[t(z)] = \det (t_{m-n})_{m,n=1}^{N}, \end{aligned}$$
(43)

i.e. it is defined simply as the determinant of the \(N\times 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 \rightarrow \infty \). By inspecting Eq. (40), we see that \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) is, up to an oscillating sign, a Toeplitz determinant of order N generated by

$$\begin{aligned} t(z)= z^{-\alpha } \frac{f(z)}{|f(z)|} . \end{aligned}$$
(44)

To go further, we consider the symbol on the unit circle \(z = \mathrm {e}^{\mathrm {i}k}\) and attempt to factorise it as:

$$\begin{aligned} t(z) = \mathrm {e}^{V(z)}t_{\mathrm{singular}}(z). \end{aligned}$$
(45)

Here \(\mathrm {e}^{V(z)} \) is called the smooth part of the symbol, which we take to mean that V(z) is analytic on the unit circleFootnote 19. This implies that the winding number of \(\exp \left( V\left( \mathrm {e}^{\mathrm {i}k}\right) \right) \) is equal to zero. The Fourier coefficients of \(V\left( \mathrm {e}^{\mathrm {i}k}\right) \) are

$$\begin{aligned} V_{n} = \frac{1}{2\pi } \int ^{2 \pi }_{0} V\left( \mathrm {e}^{\mathrm {i}k}\right) \mathrm {e}^{- \mathrm {i}n k } \mathrm {d}k \end{aligned}$$
(46)

and we define the Wiener-Hopf factorisation of \(\mathrm {e}^{V(z)}\) as:

$$\begin{aligned} \mathrm {e}^{V(z)} =b_+(z){\mathrm {e}^{V_0}}b_-(z), \quad b_+(z)=\mathrm {e}^{\sum _{n=1}^\infty V_n z^n},\; b_-(z)=\mathrm {e}^{\sum _{n=1}^\infty V_{-n} z^{-n}}. \end{aligned}$$
(47)

In our work, we will have three families of symbol to consider. The first case is \(t_{\mathrm{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_{\mathrm{singular}}(z) = z^\omega \), this is needed to represent symbols t(z) that have an integral winding number \(\omega \). Finally, the third case represents symbols t(z) with sign-change jump discontinuities. Let \(\zeta = \mathrm {e}^{\mathrm {i}\theta }\) and consider the function on the unit circle:

$$\begin{aligned} g_{\zeta ,\beta } (z) = {\left\{ \begin{array}{ll} \mathrm {e}^{\mathrm {i}\pi \beta }, &{}\quad 0 \le \mathrm {arg} z< \theta , \quad \\ \mathrm {e}^{-\mathrm {i}\pi \beta }, &{}\quad \theta \le \mathrm {arg} z < 2 \pi . \end{array}\right. } \end{aligned}$$
(48)

For \(\beta \) half-integer this is piecewise proportional to \(\mathrm {i}\), with a sign change at \(z=\zeta \) and at \(z=1\). To represent a sign-change only at \(z=\zeta \), we put \(t_{\mathrm{singular}}^{\zeta ,\beta }(z) = z^\beta g_{\zeta ,\beta } (z) \), removing the jump at \(z=1\). Conversely a jump only at \(z=1\) would be represented simply by \(z^\beta \). Notice that any half-integer \(\beta \) represents the sign-change through \(g_{\zeta ,\beta }\), but the power of \(\beta \) that appears distinguishes the \(t_{\mathrm{singular}}^{\zeta ,\beta }(z)\). The singular part of a function with several sign-change jump discontinuities can be decomposed as a product

$$\begin{aligned} t_{\mathrm{singular}}(z)=\prod _jt_{\mathrm{singular}}^{\zeta _j,\beta _j}(z) = \prod _jz^{\beta _j} g_{\zeta _j,\beta _j} (z) = z^{\sum _j \beta _j} \prod _j g_{\zeta _j,\beta _j} (z) , \end{aligned}$$
(49)

where all \(\beta _j\) are half-integer, but note that now only the total \(\sum _j \beta _j\) is fixed by the symbol we wish to represent—this redundancy has important consequences. As an example, consider the symbol \(s(k) = \mathrm {sign}(\cos (k))\), this has jump discontinuities at \(z=\pm \mathrm {i}\). Hence we should represent it by two \(\beta \) half-integer singularities, and the fact that there is no overall winding implies that \(\beta _1 = -\beta _{2}\). This gives a family of representations

$$\begin{aligned} s(z) = \mathrm {const}\times g_{\mathrm {i},\frac{2n+1}{2}}(z)g_{-\mathrm {i},-\frac{2n+1}{2}}(z) \end{aligned}$$
(50)

where \(n\in \mathbb {Z}\) and the constant fixes the correct overall sign at \(z=1\).

With these ideas in place, notice that all three families of \(t_{\mathrm{singular}}(z)\) can be represented in the same way. If we use \(t_{\mathrm{singular}}(z) = z^\beta g_{\zeta ,\beta } (z) \) as a building block, then \(t_{\mathrm{singular}}(z) =1\) is the case \(\zeta =1, \beta =0\) and \(t_{\mathrm{singular}}(z) =z^\omega \) is the case \(\zeta =1, \beta =\omega \). 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:

$$\begin{aligned} t_{\mathrm{canon}}(z) = \mathrm {e}^{V(z)} z^{\sum _{j=1}^m \beta _j}\prod _{j=1}^m g_{z_j,\beta _j}(z) z_j^{-\beta _j}, \quad z=\mathrm {e}^{\mathrm {i}k}, \; k \in [0, 2\pi ); \end{aligned}$$
(51)

where for \(j=1, \ldots , m\) and \(0 \le k_{1}< \ldots< k_{m} < 2 \pi \), we have \(z_{j} = \mathrm {e}^{\mathrm {i}k_{j}}\), \(\beta _{j} \in \frac{1}{2}\mathbb {Z}\) and the function V(z) must be smooth as above. The factor \(\prod _j z_j^{-\beta _j}\) is just a multiplicative constant and is there to align notation with [80]. Any \(\beta _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,\beta _1} = \exp (-\mathrm {i}\pi \beta _1)\). Our notation deviates slightly from reference [80], where a \(\beta _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_{\mathrm{canon}}(z)\) to which it is equal. In fact, if we find a single representation with a set of \(\{\beta _j\}\), we can find another representation by shifting each \(\beta _j \rightarrow \beta _j +n_j\) such that \(\sum 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 \(\sum _j \beta _j^2\) 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 (51), replace all \(\beta _j\) by \(\tilde{\beta }_j=\beta _j+n_j\) such that \(\sum _j n_j = 0\). This new function is the FH-representation \(t(z; n_1, \dots ,n_m)\), defined relative to the representation \(t(z; 0, \dots ,0)=t(z)\). We then have the equality:

$$\begin{aligned} t(z;n_0,\dots ,n_m) = \prod _{j=1}^m z_j^{n_j} t(z), \end{aligned}$$
(52)

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 (43). 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 \in \{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]:

Theorem 6

(Szegő 1952) Let \(t(z) = \exp (V(z))\) be a symbol, with V(z) smooth as explained above and such that \(\sum _{n= -\infty }^\infty |n| |V_n|^2<\infty \). As the matrix dimension, N, goes to infinity:

$$\begin{aligned} D_{N}[t(z)] = \exp \left( N V_0 + \sum _{n=1}^{\infty } nV_nV_{-n}\right) (1+o(1)). \end{aligned}$$
(53)

If we have a symbol with an integral winding number, i.e. \(m=1\), \(k_1=0\), \(\beta _1 \in \mathbb {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)=\mathrm {e}^{V(z)}z^{-\nu }\), where V(z) satisfies the conditions for Theorem 6. Given \(b_\pm (z)\) as defined in (47), define the auxilliary functions:

$$\begin{aligned} l(z) = \frac{b_-(z)}{b_+(z)} \qquad m(z) = \frac{b_+(1/z)}{b_-(1/z)}, \end{aligned}$$
(54)

with associated Fourier coefficientsFootnote 20\(l_k,m_k\).

For \(\nu >0\) we have:

$$\begin{aligned} D_N[z^{-\nu }\mathrm {e}^{V(z)}] = (-1)^{N\nu } D_{N+|\nu |}[\mathrm {e}^{V(z)}]\times \det \left( \begin{array}{cccc}d_{N} &{} d_{N-1} &{} \dots &{} d_{N-\nu +1} \\ d_{N+1} &{} d_{N} &{} &{} d_{N-\nu +2} \\ \vdots &{} &{} &{} \vdots \\ d_{N+\nu -1} &{} d_{N+\nu -2} &{} \dots &{} d_{N}\end{array}\right) \end{aligned}$$
(55)

where \(d_k = l_{k} + \delta _k^+\). For \(\nu <0\) we instead have \(d_k = m_k + \delta _k^-\); (55) is otherwise unchanged.

General estimates for the error terms \(\delta _k^\pm \) are given in [82]—the only case we need is as follows. Suppose that the large Fourier coefficients of \(h(z)=\mathrm {e}^{V(z)}\) behave as \(|h_n |= O(\rho ^n)\) and \(|h_{-n} |= O(\sigma ^n)\) then for large k, \(\delta _k^+ = O(\rho ^{2k}\sigma ^k)\) and \(\delta _k^- = O(\rho ^{k}\sigma ^{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) = \mathrm {e}^{V(z)}\) satisfy the conditions for Theorem 6. Then we can write

$$\begin{aligned} \log D_{N}[t(z)] = N V_0 + \sum _{n=1}^{\infty } nV_nV_{-n} + E^{(1)}_N + E^{(2)}_N, \end{aligned}$$
(56)

where, for l(z) and m(z) defined above, we have \(E^{(1)}_N = - \sum _{n=1}^\infty n l_{N+n}m_{N+n}\). The error term \(E^{(2)}_N\) is subdominant—see [82] for general estimates. For the case relevant to us, with \(\rho \) and \(\sigma \) defined in Theorem 7, we have \(E^{(1)}_N = O(\rho ^N \sigma ^N)\) and \(E^{(2)}_N = O(\rho ^{2N} \sigma ^{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 [32]; 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.

Theorem 9

(Deift, Its, Krasovsky 2011) Consider a Toeplitz matrix generated by t(z) in the canonical form (51). Suppose \(\beta _j \not \in \mathbb {Z}\) for all j. Then, as the matrix dimension, N, goes to infinity:

$$\begin{aligned} D_{N}[t(z)]=\sum _{\begin{array}{c} \mathrm {Dominant}\\ \mathrm {FH-reps:}~\{n_j\} \end{array}}\left( \prod _{j=1}^m z_j^{n_jN}\right) \mathcal {R}(t(z;\{n_j\})(1+o(1)). \end{aligned}$$
(57)

Where:

$$\begin{aligned} \mathcal {R}(t(z;\{n_j\}))&=N^{-\sum _{j=1}^m \tilde{\beta }_j^2}\exp \left( NV_0 + \sum _{n=1}^\infty nV_nV_{-n}\right) \prod _{1\le i<j\le m} |z_i-z_j|^{2\tilde{\beta }_i\tilde{\beta }_j}\nonumber \\&\quad \times \prod _{j=1}^m b_+(z_j)^{\tilde{\beta }_j}b_-(z_j)^{-\tilde{\beta }_j}\prod _{j=1}^mG(1+\tilde{\beta }_j) G(1-\tilde{\beta }_j).\nonumber \\&\qquad \quad (\hbox {Recalling that}~~~\tilde{\beta }_j=\beta _j+n_j.) \end{aligned}$$
(58)

The \(V_n\) are unaltered when passing between FH-reps. Branches of \(b_{\pm }(z_j)^{\tilde{\beta }_j}\) are determined by \(b_\pm (z_j)^{\tilde{\beta }_j} = \mathrm {e}^{\tilde{\beta }_j\sum _{n=1}^\infty V_{\pm n} z_j^{\pm n}}\). G(z) is the Barnes G-function [84, Sect. 5.17]; given as a Weierstrass product by

$$\begin{aligned} G\left( z+1\right) =(2\pi )^{z/2}\mathrm {e}^{-\frac{1}{2}z(z+1)-\frac{1}{2}\gamma _E z^{2}} \prod _{j=1}^{\infty }\left( \left( 1+\frac{z}{j}\right) ^{j}\exp \,\left( -z+\frac{ z^{2}}{2j}\right) \right) , \end{aligned}$$
(59)

where \(\gamma _E\) is the Euler-Mascheroni constant. It is clear from Eq. (59) that G vanishes whenever the argument is a negative integer. Hence if \(\beta _j \in \mathbb {Z}\), the RHS of (58) 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

$$\begin{aligned} f(z)&= \rho \frac{1}{z^{N_p}} \prod _{i=1}^{N_z} (z-z_i) \prod _{j=1}^{N_Z} (z-Z_j) \end{aligned}$$
(60)
$$\begin{aligned}&= \left( \rho \prod _{j=1}^{N_Z} (-Z_j) \right) z^\omega \prod _{i=1}^{N_z} (1-z_i/z) \prod _{j=1}^{N_Z} (1-z/Z_j) =: \rho ' z^\omega f_0(z). \end{aligned}$$
(61)

\(\rho ' = \rho \prod _{j=1}^{N_Z} (-Z_j)\), 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

$$\begin{aligned} s =\mathrm {sign}(\rho )\times (-1)^{N_{Z}^+}. \end{aligned}$$
(62)

If we consider \((-1)^{N(\omega -1)}\langle \mathcal {O}_\omega (1)\mathcal {O}_\omega (N)\rangle \), then this is generated by \(t(z) = \mathrm {e}^{V(z)}\), where

$$\begin{aligned} V(z) - V_0 = \frac{1}{2} \left( \log f_0(z) - \log f_0(\overline{z})\right) \end{aligned}$$
(63)

for a continuous logarithm that could be found by integrating the logarithmic derivative of f. We instead jump to the following solution:

$$\begin{aligned} V(z) - V_0&= \frac{1}{2} \sum _{i,j} \mathrm {Log}(1-z_i/z)-\mathrm {Log}(1-z_i/\overline{z})+ \mathrm {Log}(1-z/Z_{j})-\mathrm {Log}(1-\overline{z}/{Z_j}) \end{aligned}$$
(64)
$$\begin{aligned}&= -\frac{1}{2} \sum _{n=1}^\infty \frac{1}{n}\left( \sum _{i=1}^{N_z} \left( \frac{z_i}{z}\right) ^n-\left( \frac{z_i}{\overline{z}}\right) ^n+ \sum _{j=1}^{N_Z}\left( \frac{z}{Z_j}\right) ^n-\left( \frac{\overline{z}}{Z_j}\right) ^n\right) , \end{aligned}$$
(65)

where the function \(\mathrm {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=\mathrm {e}^{\mathrm {i}k}\) we can put \(\overline{z}=1/z\) into (65).

This gives us an honest V(z) from which one can read off the Fourier coefficients:

$$\begin{aligned} V_0&= \log s = 0, \mathrm {i}\pi \end{aligned}$$
(66)
$$\begin{aligned} V_n&={\left\{ \begin{array}{ll} \frac{1}{2n}\left( \sum _iz_i^{n} - \sum _jZ_j^{-n}\right) &{}\quad n>0\\ -\frac{1}{2n}\left( \sum _i z_i^{n} - \sum _jZ_j^{-n}\right) &{}\quad n<0 . \end{array}\right. } \end{aligned}$$
(67)

Inserting into Theorem 6 we reach:

$$\begin{aligned} \det [t(z)] = s^N \exp \left( \sum _{n=1}^\infty - \frac{1}{4n}{\left( \sum _{i=1}^{N_z}z_i^{n} -\sum _{j=1}^{N_Z}Z_j^{-n}\right) }^{2}~ \right) . \end{aligned}$$
(68)

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 \(\mathrm {e}^{\mathrm {i}\pi N(\omega -1) + N\log (s)}\).

Note that with the Fourier coefficients of V in hand, we can find the Wiener-Hopf decomposition (47) of our symbol when z is on the unit circle.

$$\begin{aligned} b_+(z)&= \mathrm {e}^{\sum _{n=1}^\infty \frac{1}{2n}\left( \sum _iz_i^{n} - \sum _jZ_j^{-n}\right) z^n} =\prod _{i=1}^{N_z}\mathrm {e}^{-\frac{1}{2} \mathrm {Log}(1-zz_i)}\prod _{j=1}^{N_Z}\mathrm {e}^{\frac{1}{2} \mathrm {Log}(1-z/Z_j)} \nonumber \\ b_-(z)&= \mathrm {e}^{-\sum _{n=1}^\infty \frac{1}{2n}\left( \sum _i z_i^{n} - \sum _jZ_j^{-n}\right) z^{-n}} =\prod _{i=1}^{N_z}\mathrm {e}^{\frac{1}{2} \mathrm {Log}(1-z_i/z)}\prod _{j=1}^{N_Z}\mathrm {e}^{-\frac{1}{2} \mathrm {Log}(1-1/(zZ_j))}. \end{aligned}$$
(69)

Note that \(1/b_+(z)=b_{-}(1/z)\). Moreover,

$$\begin{aligned} b_+(z) = \sqrt{\frac{\prod _{j=1}^{N_Z}(1-z Z_j^{-1})}{\prod _{i=1}^{N_z}(1-zz_i)}} \end{aligned}$$
(70)

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\).

6.2 Correlation Lengths

We now use Theorem 7 to find the behaviour of the correlation function \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) in the gapped phase \(\omega \). For definiteness, let us label the zeros by proximity to the unit circle: \(|Z_i| \le |Z_j|\) and \(|z_i| \ge |z_j|\) for \(i<j\).

6.2.1 Asymptotics of l(z), m(z)

The key ingredient that we need are the asymptotically large Fourier coefficients of the auxilliary functions

$$\begin{aligned} l(z)&= \frac{b_-(z)}{b_+(z)} =\sqrt{\frac{\prod _{i=1}^{N_z}(1-zz_i)(1-z_i/z)}{\prod _{j=1}^{N_Z}(1-zZ_j^{-1}))(1-Z_j^{-1}/z)}} \end{aligned}$$
(71)
$$\begin{aligned} m(z)&= 1/l(z). \end{aligned}$$
(72)

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_i^{-1},Z_j,Z_j^{-1}\}\) 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:

$$\begin{aligned} l_n&= \frac{1}{2\pi } \int _0^{2\pi } l(k) \exp (-\mathrm {i}n k) \mathrm {d}k\nonumber \\&=\frac{1}{2\pi \mathrm {i}} \int _{S^1}\left( \frac{\prod _{i=1}^{N_z}(1-sz_i)(1-z_i/s)}{\prod _{j=1}^{N_Z}(1-sZ_j^{-1})(1-Z_j^{-1}/s)}\right) ^{1/2}~ s^{-(n+1)}\mathrm {d}s. \end{aligned}$$
(73)

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 Fig. 4—we call the nearest branch point(s) \(s_1\) (and \(\overline{s}_1\)), for \(\arg (s_1)\in [0,\pi ]\). We connect up the radial cuts outside a circle of large radius, the precise choice is unimportant.

Fig. 4
figure 4

Schematic for the two generic cases of the computation (73). 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. \(\times \) indicates the closest branch cuts to the unit circle that are outside the circle

In case A we consider the Hankel contour connecting infinity to the nearest real zero and back —this is exactly the relevant part of the snagged contour. After parameterising \(s=s_1\mathrm {e}^t\) for \(t \in \mathbb {R}_+\) and where \(\arg (t) =0\) below the axis and \(\arg (t) = -2\pi \) 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 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. Then,

$$\begin{aligned} l_n = {\left\{ \begin{array}{ll} - z_1^n\frac{1}{n^{3/2}}\underbrace{\frac{1}{2\sqrt{\pi }}\left( \frac{(1-z_1^2)\prod _{i=2}^{N_z}(1-z_i/z_1)(1-z_1z_i)}{\prod _{j=1}^{N_Z}(1-Z_j^{-1}/z_1)(1-z_1Z_j^{-1})}\right) ^{1/2}}_{\lambda }\left( 1+O(1/n)\right) \qquad &{}s_1=1/z_1\\ Z_1^{-n}\frac{1}{n^{1/2}}\frac{1}{\sqrt{\pi }}\left( \frac{\prod _{i=1}^{N_z}(1-Z_1z_i)(1-z_i/Z_1)}{(1-Z_1^{-2})\prod _{j=2}^{N_Z}(1-Z_1Z_j^{-1})(1-(Z_1Z_j)^{-1})}\right) ^{1/2}\left( 1+O(1/n)\right) \qquad &{}s_1=Z_1\end{array}\right. } \end{aligned}$$
(74)

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

$$\begin{aligned} l_n = {\left\{ \begin{array}{ll} \frac{1}{\sqrt{\pi }}{\text {Im}}(c_1 z_1^n )\frac{1}{n^{3/2}}\left( 1+O(1/n)\right) &{}\quad s_1=1/z_1\\ \frac{2}{\sqrt{\pi }}{\text {Im}}(c_2 Z_1^{-n})\frac{1}{n^{1/2}}\left( 1+O(1/n)\right) &{}\quad s_1=Z_1\end{array}\right. } \end{aligned}$$
(75)

for

$$\begin{aligned} c_1 =&-\left( -\frac{(1-z_1^2)(1-z_1\overline{z}_1)(1-\overline{z}_1/z_1)\prod _{i=2}^{N_z}(1-z_i/z_1)(1-z_1z_i)}{\prod _{j=1}^{N_Z}(1-Z_j^{-1}/z_1)(1-z_1Z_j^{-1})}\right) ^{1/2}\nonumber \\ c_2 =&-\left( -\frac{\prod _{i=1}^{N_z}(1-Z_1z_i)(1-z_i/Z_1)}{(1{-}Z_1^{-2})(1-Z_1/\overline{Z}_1)(1-(Z_1\overline{Z}_1)^{-1})\prod _{j=2}^{N_Z}(1{-}Z_1Z_j^{-1})(1{-}(Z_1Z_j)^{-1})}\right) ^{1/2}. \end{aligned}$$
(76)

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:

$$\begin{aligned} m_n = {\left\{ \begin{array}{ll} - Z_1^{-n}\frac{1}{n^{3/2}}\frac{1}{2\sqrt{\pi }}\left( \frac{(1-Z_1^{-2})\prod _{j=2}^{N_Z}(1-Z_1Z_j^{-1})(1-(Z_1Z_j)^{-1})}{\prod _{i=1}^{N_z}(1-Z_1z_i)(1-z_i/Z_1)}\right) ^{1/2}\left( 1+O(1/n)\right) \qquad &{}s_1=Z_1\\ z_1^{n}\frac{1}{n^{1/2}}\underbrace{\frac{1}{\sqrt{\pi }}\left( \frac{\prod _{j=1}^{N_Z}(1-Z_j^{-1}/z_1)(1-z_1Z_j^{-1})}{(1-z_1^2)\prod _{i=2}^{N_z}(1-z_i/z_1)(1-z_1z_i)}\right) ^{1/2}}_{\kappa }\left( 1+O(1/n)\right) \qquad &{}s_1=1/z_1.\end{array}\right. } \end{aligned}$$
(77)

For a complex conjugate pair of nearest singularities we have:

$$\begin{aligned} m_n = {\left\{ \begin{array}{ll} \frac{1}{\sqrt{\pi }}{\text {Im}}(c_2^{-1} Z_1^{-n} )\frac{1}{n^{3/2}}\left( 1+O(1/n)\right) &{}\quad s_1=Z_1\\ \frac{2}{\sqrt{\pi }}{\text {Im}}(c_1^{-1} z_1^{n})\frac{1}{n^{1/2}}\left( 1+O(1/n)\right) &{}\quad s_1=1/z_1,\end{array}\right. } \end{aligned}$$
(78)

where the \(c_i\) are defined in (76).

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 \(\mathrm {e}^{-n/\xi }\), using Cauchy’s theorem, with no algebraic prefactor (recall that \(\xi = 1/|\log |\zeta _\star ||\) where \(\zeta _\star \) 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 (73). Importantly, even in these exceptional cases, the nearest zero always sets the longest correlation length for the operators \(\mathcal {O}_\alpha \). 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 \(\xi \) 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 \(\mathrm {e}^{- n/\xi }\) 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 zeros which are mutually inverse. For example, if the closest zeros are at a and at \(1/a\in \mathbb {R}\). This is the only case where \(\xi \) defined in terms of one of these closest zeros 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_\pm (z)\) and so they do not contribute to the asymptotics of any\(\mathcal {O}_\alpha \). 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_\pm (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 \(\xi \). 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 \(\delta _N^\pm \). To do so, we need to find \(\rho \) and \(\sigma \) such that for \(h(z) = \mathrm {e}^{V(z)}\), \(|h_n |= O(\rho ^n)\) and \(|h_{-n}|= O(\sigma ^n)\) for large n. Recall that the relevant \(\mathrm {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, \(\rho = \sigma = 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 \(\langle \mathcal {O}_{\alpha }(1)\mathcal {O}_{\alpha }(N+1)\rangle \) —Proof of Theorem 2

Suppose that we are in the phase \(\omega \), then the generating function of the correlator is \(s\mathrm {e}^{V(z)}z^{\omega -\alpha }\). In the case \(\omega -\alpha >0\), using Theorem 7 we have that

$$\begin{aligned} \langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle&= (-1)^{N(\omega -1)} D_{N+\omega -\alpha }(s \mathrm {e}^{V(z)})\nonumber \\&\quad \times \det \left( \begin{array}{cccc}m_{N} &{} m_{N-1} &{} \dots &{} m_{N-(\omega -\alpha )+1} \\ m_{N+1} &{} m_{N} &{} &{} m_{N-(\omega -\alpha )+2} \\ \vdots &{} &{} &{} \vdots \\ m_{N+(\omega -\alpha )-1} &{} m_{N+(\omega -\alpha )-2} &{} \dots &{} m_{N}\end{array}\right) . \end{aligned}$$
(79)

The large determinant \(D_{N+\omega -\alpha }(s \mathrm {e}^{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 =\varTheta (\mathrm {e}^{-N/\xi })\). Then, in the position (ij) of the second matrix we have a factor of \(\mathrm {e}^{(-N+i -j)/\xi }\). The row and column index multiplicatively decouple, and so any individual term of the Laplace expansion of the determinant contains a factor of \(\mathrm {e}^{-N(\omega -\alpha )/\xi }\), hence we may factor this out and we have that:

$$\begin{aligned} \langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle&= \mathrm {e}^{\mathrm {i}\pi N(\omega -1)+N \log s}\left( \lim _{R\rightarrow \infty }|\langle \mathcal {O}_\omega (1)\mathcal {O}_\omega (R)\rangle |\right) \nonumber \\&\qquad \times \mathrm {e}^{-N(\omega -\alpha )/\xi }\underbrace{\det \left( (N+i-j)^{-K} \alpha _{N+i-j}\right) _{i,j=1}^{\omega -\alpha }}_{\det M(N)}. \end{aligned}$$
(80)

The matrix elements of M(N) are derived from the propositions above: i.e. \(K = 1/2\) or 3 / 2 and \(\alpha _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 \(\omega -\alpha <0\) the same calculation goes through with \(m_n\) replaced by \(l_n\) (and the second matrix has dimension \(|\omega -\alpha |\)). We have hence proved Theorem 2.

Now, putting together Theorems 1b and 2 prove Theorem 1a. In particular, we have shown that the correlators \(|\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle |\) do indeed form a set of order parameters that distinguish \(\omega \). 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 (80), this oscillation depends on both \(\omega \) and s, as defined in (62) (one must also take into account oscillations coming from \(\det M\)).

6.2.4 The Correlation Length of \(\langle \mathcal {O}_{\omega }(1)\mathcal {O}_{\omega }(N)\rangle \) in the phase \(\omega \)—Proof of Theorem 3

The proof follows from Theorem 8 and our calculations above. Firstly, using 6.2.2 we have that \(E^{(2)}_N\) is exponentially subdominant compared to \(E^{(1)}_N\). We do not evaluate \(E^{(1)}_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 [32]. Thus, in the generic case, we have \(E^{(1)}_N=O(|s_1|^{-2N}/N^2)\). 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 [31, 51] 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 (42) 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 [53] 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 \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) 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 \(\alpha \in \{-1,0,1\}\) in [56]. 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)^{N(\alpha -1)}\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \): \(z^{-\alpha } f(z)/|f(z)|\) and find the dominant Fisher-Hartwig representations. This goes as follows:

$$\begin{aligned} z^{-\alpha } f(z)/|f(z)|&= s z^{-\alpha } \frac{z^{N_z}}{z^{N_p}} \prod _{j=1}^{N_z} \frac{(1-z_j/z)}{|1-z_j/z|}\prod _{j=1}^{2c} \frac{(z-\mathrm {e}^{\mathrm {i}k_j})}{|z-\mathrm {e}^{\mathrm {i}k_j}|}\prod _{j=1}^{N_Z} \frac{(1-z/Z_j)}{|1-z/Z_j|} \end{aligned}$$
(81)
$$\begin{aligned}&= \underbrace{ sC\prod _{j=1}^{N_z} \frac{(1-z_j/z)}{|1-z_j/z|}\prod _{j=1}^{N_Z} \frac{(1-z/Z_j)}{|1-z/Z_j|}}_{\mathrm {e}^{V(z)}} \underbrace{C^{-1}z^{\omega -\alpha }\prod _{j=1}^{2c} \frac{(z-\mathrm {e}^{\mathrm {i}k_j})}{|z-\mathrm {e}^{\mathrm {i}k_j}|}}_{\mathrm{singular}}. \end{aligned}$$
(82)

We reemphasise that in this analysis we pick the phase of the complex zeros such that \(k_j \in [0,2\pi )\). The smooth part \(\mathrm {e}^{V(z)}\) can be analysed as in the gapped case—in particular, the Fourier coefficients for \(n\ne 0\) are given by (67) (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= \mathrm {e}^{\mathrm {i}k}\), the analysis of the singular part is split into three cases: a real zero at \(k = 0,\pi \), a pair of complex conjugate zeros at \(k = \phi \) and \(k=2\pi -\phi \), or a set of zeros of multiplicity greater than one. The third (fine-tuned) case is discussed in Sect. 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:

$$\begin{aligned} \frac{\exp (\mathrm {i}k)\pm 1}{|\exp (\mathrm {i}k)\pm 1|} = \exp (\mathrm {i}k/2) \times {\left\{ \begin{array}{ll} \cos (k/2)/|\cos (k/2)| \\ \mathrm {i}\sin (k/2)/|\sin (k/2)| = \mathrm {i}. \end{array}\right. } \end{aligned}$$
(83)

For a zero at \(-1\) we have a sign-change discontinuity at \(k = \pi \), and a zero at 1 has a sign-change-typeFootnote 21 singularity at \(k=0\). For a complex conjugate pair of zeros at \(\exp (\pm \mathrm {i}\phi )\) we have:

$$\begin{aligned} \frac{(\mathrm {e}^{\mathrm {i}k}-\mathrm {e}^{\mathrm {i}\phi })(\mathrm {e}^{\mathrm {i}k}-\mathrm {e}^{-\mathrm {i}\phi })}{|\mathrm {e}^{\mathrm {i}k}-\mathrm {e}^{\mathrm {i}\phi }||\mathrm {e}^{\mathrm {i}k}-\mathrm {e}^{-\mathrm {i}\phi }|} = \exp (\mathrm {i}k) \times \mathrm {sign}(\cos (k)-\cos (\phi )) . \end{aligned}$$
(84)

Since \(\phi \ne 0\) or \(\pi \), \(\mathrm {sign}(\cos (k)-\cos (\phi ))\) has sign-change discontinuities at \(k = \phi \) and \(k=2\pi -\phi \). We conclude that every zero contributes a factor \(\exp (\mathrm {i}k/2)\) as well as a sign-change at the location of the zero, which we can represent with a \(g_{k_j,\beta _j}(z)\) for \(\beta _j\) any half-integer.

Putting this information back into the symbol we reach

$$\begin{aligned} z^{-\alpha } f(z)/|f(z)|&= \underbrace{ s C\prod _{j=1}^{N_z} \frac{(1-z_j/z)}{|1-z_j/z|}\prod _{j=1}^{N_Z} \frac{(1-z/Z_j)}{|1-z/Z_j|}}_{\mathrm {e}^{V(z)}} \underbrace{z^{c+\omega -\alpha }\prod _{j=1}^{2c} g_{k_j,\beta _j}(z)\mathrm {e}^{-\mathrm {i}\beta _jk_j}}_{\mathrm{singular}},\nonumber \\ \end{aligned}$$
(85)

where the \(\beta _j\) are half-integer and

$$\begin{aligned} \sum _{j=1}^{2c} \beta _j = c+\omega -\alpha . \end{aligned}$$
(86)

We need to fix the multiplicative constant, C, by noting that the singular factors, isolated above, usually jump between \(\pm 1\), rather than \(\pm \mathrm {i}\), and we also need to include \(\prod \mathrm {e}^{-\mathrm {i}\beta _jk_j}\) in the singular part of (51). This leads to \(C=\prod _{j=1}^{2c}\mathrm {e}^{\mathrm {i}\beta _jk_j + \delta _{k_j,0}\mathrm {i}\pi /2}/g_{k_j,\beta _j}(1) \).

To find the asymptotics we need some minimal representation (a set of \(\beta _j\) minimising \(\sum _j \beta _j^2\)), 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+\omega -\alpha =(2m-1)c\), for \(m\in \mathbb {Z}\), where the minimal solution is unique: \(\beta _j = \frac{2m-1}{2}\) for all j. If we have \((2m-1)c<c+\omega -\alpha <(2m+1)c\) we form the set of minimal FH-reps by starting from \(\tilde{\beta }_j = \frac{2m-1}{2}\) and sending \(\tilde{\beta }_j \rightarrow \tilde{\beta }_j+1\) for

$$\begin{aligned} c+\omega -\alpha -(2m-1)c = (1-m)2c+\omega -\alpha \in \mathbb {Z} \end{aligned}$$
(87)

of the \(\tilde{\beta }_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+\omega -\alpha \) of the \(\beta _j\). There are

$$\begin{aligned} \left( {\begin{array}{c}2c\\ (1-m)2c+\omega -\alpha \end{array}}\right) = \frac{(2c)!}{((1-m)2c+\omega -\alpha )!(2mc-\omega +\alpha )!} \end{aligned}$$
(88)

minimal FH-reps in total. Given the parameters \(\omega \), \(\alpha \), c we get that \(m = 1 + \lfloor \frac{\omega - \alpha }{2c} \rfloor \), where \(\lfloor x \rfloor \) denotes the greatest integer less than or equal to x. Theorem 9 immediately gives the dominant scaling

$$\begin{aligned} |\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle | = \mathrm {const}\times N^{-2\varDelta }(1+o(1)), \end{aligned}$$
(89)

where

$$\begin{aligned} \varDelta _\alpha (c,\omega ) = \frac{1}{2}\left( (2cm-\omega +\alpha )\frac{(2m-1)^2}{4}+(2c(1-m)+\omega -\alpha )\frac{(2m+1)^2}{4}\right) . \end{aligned}$$
(90)

This formula for \(\varDelta _\alpha \) obscures some features of this function, to bring them out define \(\tilde{\alpha } = \alpha - (\omega +c)\), then one can show that

$$\begin{aligned} \varDelta _\alpha (c,\omega ) = c \left( \frac{1}{4} + x^2-(x - [x])^2 \right) \Big |_{x = \tilde{\alpha } / 2c} \end{aligned}$$
(91)

where [x] is the nearest integer to x. It is then clear that \(\varDelta _\alpha \) is symmetric under \(\tilde{\alpha } \leftrightarrow -\tilde{\alpha }\) 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 \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) at large N. Firstly, note that:

$$\begin{aligned} \mathrm {e}^{V(z)} = \underbrace{s \prod _{j=1}^{2c}\mathrm {e}^{\mathrm {i}\beta _jk_j}/g_{k_j,\beta _j}(1)}_{\mathrm {e}^{V_0}}\underbrace{\prod _{j=1}^{N_z} \frac{(1-z_j/z)}{|1-z_j/z|}\prod _{j=1}^{N_Z} \frac{(1-z/Z_j)}{|1-z/Z_j|}}_{\mathrm {e}^{V(z)-V_0}}. \end{aligned}$$
(92)

The first factor is a pure phase and contributes to the asymptotics as \(\mathrm {e}^{N V_0}\). By inspectionFootnote 22

$$\begin{aligned} V_0 = \mathrm {Log}(s) + \mathrm {i}\sum _{j=1}^{2c}(\beta _j(k_j-\pi ) + (2\beta _j\pi +\pi /2) \delta _{k_j,0}), \end{aligned}$$
(93)

it is important to emphasise that this is an imaginary number, and again recall that \(k_j \in [0,2\pi )\). The second factor contributes in two ways: firstly through \(\mathrm {e}^{\sum _n n V_nV_{-n}}\), exactly the quantity we calculated in Sect. 6.1. Secondly, we need powers of the Wiener-Hopf factors that were derived at the end of Sect. 6.1. Putting this all together we get:

Theorem 10

$$\begin{aligned} \langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle = N^{-2\varDelta _\alpha }\mathrm {e}^{\mathrm {i}N K}\sum _{\{n_j\}}\left( \prod _{j=1}^{2c} \mathrm {e}^{\mathrm {i}Nk_j n_j}\right) \mathcal {C}(\{\beta _j+n_j\}) (1+o(1)), \end{aligned}$$
(94)

where the sum is over all dominant FH-reps, these are parameterised by \(\{n_j\}\) and defined in Sect. 7.1. \(\varDelta _\alpha \) is \(\varDelta _\alpha (c,\omega )\) given in (91) and K is equal to \(-\mathrm {i}V_0 +\pi (\alpha -1)\). The representation dependent O(1) multiplier is given by

$$\begin{aligned} \mathcal {C}(\{\tilde{\beta }_j\})&= \left( \frac{ \prod _{i_1,i_2=1}^{N_z} (1-z_{i_1} z_{i_2}) \prod _{j_1,j_2=1}^{N_Z} \left( 1-\frac{1}{Z_{j_1} Z_{j_2}}\right) }{\prod _{i=1}^{N_z}\prod _{j=1}^{N_Z} \left( 1-\frac{z_i}{ Z_j}\right) ^2}\right) ^{1/4}\prod _{1\le i<j\le 2c} |\mathrm {e}^{\mathrm {i}k_i}-\mathrm {e}^{\mathrm {i}k_j}|^{2\tilde{\beta }_i\tilde{\beta }_j} \nonumber \\&\quad \times \prod _{j=1}^{2c} \left( \frac{\prod _{l=1}^{N_Z}(1-\mathrm {e}^{\mathrm {i}k_j} Z_l^{-1})(1-\mathrm {e}^{-\mathrm {i}k_j}Z_l^{-1})}{\prod _{i=1}^{N_z}(1-\mathrm {e}^{\mathrm {i}k_j}z_i)(1-\mathrm {e}^{-\mathrm {i}k_j}z_i)}\right) ^{\tilde{\beta }_j/2}\prod _{j=1}^{2c}G(1+\tilde{\beta }_j) G(1-\tilde{\beta }_j). \end{aligned}$$
(95)

In our construction of the dominant FH-reps we start from setting all \(\tilde{\beta }_j\) to be equal and then add 1 to a fixed number of them. This means that the difference \(\tilde{\beta }_i - \tilde{\beta }_j\) is either 0 , 1 or \(-1\) in all dominant FH-reps. Hence, pairs of complex conjugate zeros \(\mathrm {e}^{\mathrm {i}k_i }=\mathrm {e}^{-\mathrm {i}k_j}\) contribute \(\mathrm {e}^{\mathrm {i}n k_i N}\) to the oscillatory factor in (94), where \(n \in \{0,\pm 1\}\). We discuss the non-universal multiplier (95) 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

$$\begin{aligned} \left( \frac{(\mathrm {e}^{\mathrm {i}k}-\mathrm {e}^{\mathrm {i}\phi })(\mathrm {e}^{\mathrm {i}k}-\mathrm {e}^{-\mathrm {i}\phi })}{|\mathrm {e}^{\mathrm {i}k}-\mathrm {e}^{\mathrm {i}\phi }||\mathrm {e}^{\mathrm {i}k}-\mathrm {e}^{-\mathrm {i}\phi }|}\right) ^m = \exp (\mathrm {i}m k) \times \biggl (\mathrm {sign}(\cos (k)-\cos (\phi ))\biggr )^m . \end{aligned}$$
(96)

Equation (83) 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 Fisher-Hartwig canonical form with half-integer \(\beta \) singularities at \(\mathrm {e}^{\pm \mathrm {i}\phi }\). 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, there is no jump discontinuity and we do not analyse this here. The multicritical point in the XY model, with \(f(z) = (z+1)^2\) is an example of such a case.

8 Extensions of Our Results

8.1 Long-Range Chains

In this section we discuss the effect of allowing our model (2) to have non-zero 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^\infty \) smooth part, and a well defined winding number. We still have that \(\omega = N_z - N_p\), but note that poles are no longer restricted to the origin. The theory of Sect. 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

$$\begin{aligned} f(z) = \exp (V(z)) z^\beta , \end{aligned}$$

i.e. a symbol with a single Fisher-Hartwig jump singularity at \(z=1\), and demand that \(\exp (V(z))\) is a \(C^\infty \) function. Then:

$$\begin{aligned} D_{N}(f(z)) = \exp (N V_0) N^{-\beta ^2}(E+o(1)) \end{aligned}$$
(97)

where E is the constant defined implicitly in (58).

For nonvanishing E, or equivalently \(\beta \not \in \mathbb {Z}\), this is a special case of Theorem 9. However, for \(\beta \in \mathbb {Z}\) this gives us a concrete asymptotic bound on the Toeplitz determinant in the case of a symbol with a \(C^\infty \) 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^\infty \) symbol. In particular, in the phase \(\omega \) we have that \(\langle \mathcal {O}_\omega (1)\mathcal {O}_\omega (N+1)\rangle \) tends to a non-zero value that can be calculated using Szegő’s theorem. Moreover, Theorem 11 proves that \(\langle \mathcal {O}_\alpha (1)\mathcal {O}_\alpha (N+1)\rangle \) for all \(\alpha \ne \omega \) 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 \(\alpha \)th correlator in the phase \(\omega \) is \(O(\mathrm {e}^{-N/\xi _\alpha })\), where \(\xi _\alpha \) is defined as above and \(\xi \) 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 (51)). Note that a study of the critical scaling of entanglement entropy for the isotropic subclass of such chains is included in [42]. 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 \(\mathrm {e}^{\pm \mathrm {i}k}\) 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 expectFootnote 23 similar behaviour in more general transitions where, as well as the approaching zeros, there are additional ‘spectator’ zeros on the unit circle.

9 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 Fig. 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 \(\xi _\alpha \)—this allowed for the extraction of the topological invariant \(\omega \). This was more striking for the dual spin chains, where local observables can be used to measure \(\omega \). 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 \(\omega \) 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.Footnote 24 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.

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.