]Exact correlations in topological quantum chains

N. G. Jones
Mathematical Institute, University of Oxford, Oxford, OX2 6GG, UK

and the Heilbronn Institute for Mathematical Research, Bristol, UK
R. Verresen
Department of Physics, Harvard University, Cambridge, MA 02138, USA

Although free-fermion systems are considered exactly solvable, they generically do not admit closed expressions for nonlocal quantities such as topological string correlations or entanglement measures. We derive closed expressions for such quantities for a dense subclass of certain classes of topological fermionic wires (classes BDI and AIII). Our results also apply to spin chains called generalised cluster models.  While there is a bijection between general models in these classes and Laurent polynomials, restricting to polynomials with degenerate zeros leads to a plethora of exact results: (1) we derive closed expressions for the string correlation functions—the order parameters for the topological phases in these classes; (2) we obtain an exact formula for the characteristic polynomial of the correlation matrix, giving insight into ground state entanglement; (3) the latter implies that the ground state can be described by a matrix product state (MPS) with a finite bond dimension in the thermodynamic limit—an independent and explicit construction for the BDI class is given in a concurrent work (Jones, Bibo, Jobst, Pollmann, Smith, Verresen, Phys. Rev. Res. 3 033265 (2021)); (4) for BDI models with even integer topological invariant, all non-zero eigenvalues of the transfer matrix are identified as products of zeros and inverse zeros of the aforementioned polynomial. General models in these classes can be obtained by taking limits of the models we analyse, giving a further application of our results. To the best of our knowledge, these results constitute the first application of Day's formula and Gorodetsky's formula for Toeplitz determinants to many-body quantum physics.

§ INTRODUCTION

Modelling certain many-body quantum systems in terms of non-interacting fermions has proved to be unreasonably effective. Historically forming the basis of Fermi liquid theory <cit.>, it has in recent times led to the discovery of topological insulators and superconductors <cit.>. Free-fermion chains—dual to famous spin models such as the two-dimensional classical Ising model and one-dimensional quantum XY model—have also been fertile ground for mathematical innovations. Whilst correlation functions of local fermionic operators are simple, non-local string correlation functions—essential for calculating certain local spin correlations as well as understanding the topological nature of fermionic systems—required and indeed stimulated significant developments in the theory of the asymptotics of Toeplitz determinants  (determinants of matrices that are constant along each diagonal) <cit.>.

Considering only the two-dimensional classical Ising model as an example, there are a vast number of remarkable exact results for spin correlations <cit.>. However, it is unusual to find simple closed formulae for these correlators. Typically, in analysing string correlations or the entanglement properties of fermionic systems, one finds closed expressions for dominant terms in the asymptotic expansion for, respectively, large string length <cit.> or large subsystem size (for the entanglement entropy) <cit.>. A rare and famous exception[Another exceptional case, less relevant to this work, is the spin-spin correlator `along the diagonal' at $T=T_c$ in the 2D classical Ising model (equivalently, at the phase transition in the 1D quantum transverse-field Ising model). This can be expressed in an exact closed form as a ratio of Barnes G-functions <cit.>.] that does admit simple closed formulae for arbitrary distances is the `disorder line' in the XY spin chain—this is dual to a nearest-neighbour free-fermion model and there are two symmetry-broken product state ground states <cit.>.

In the present work, we consider Hamiltonians for one-dimensional classes of spinless fermions with phases classified by a topological winding number (commonly referred to as the BDI and AIII classes <cit.>). The BDI class consists of superconducting time-reversal invariant spinless fermions, while the AIII class consists of charge-conserving spinless fermions with a sublattice symmetry; these classes are non-interacting instances of symmetry-protected topological (SPT) phases <cit.>. We show that for a simply characterised subclass of both the BDI and AIII classes, we can find closed expressions for certain physically relevant string correlation functions for models within this subclass. Moreover, we give an exact formula for the characteristic polynomial of the correlation matrix—intimately related to entanglement properties of the system—where this formula depends only parametrically on the size of the matrix. We further show that our subclass of interest is dense in the whole class, meaning that any Hamiltonian in the BDI or AIII class can be obtained as the limit of a sequence of Hamiltonians studied in this work.

The above results allow us to prove the existence of an exact matrix product state (MPS) representation <cit.> of the ground state for this dense subclass.
An MPS is one of the simplest tensor networks for describing many-body wavefunctions and has proven to be a valuable concept both analytically (e.g. in the discovery of fixed-point SPT states <cit.>) and as a numerical ansatz <cit.>. In a concurrent work <cit.>, we give an explicit construction for the MPS ground state in this subclass of the BDI class (up to a measure zero set of exceptional cases), where this subclass is referred to as the MPS skeleton. Particular examples of models in this subclass were previously studied in <cit.>. While the explicit MPS representation is a valuable tool, it is not straightforward to extract closed analytic expressions for the string correlations and entanglement spectrum, motivating the analysis in this work. Moreover, we prove the existence of an MPS for a subclass of class AIII, where the construction in <cit.> does not directly apply. It is interesting that the techniques used here, from Toeplitz determinant theory, are rather different to those used in the explicit construction, which is based on Witten's conjugation argument for frustration-free models <cit.>. Both are in turn different to approaches to free-fermion MPS based on Gaussian states <cit.>. One can think of the MPS ground state as two $\chi \times \chi$ matrices, where $\chi$ is the bond dimension of the MPS. As well as proving that $\chi < \infty$ (even for infinite system sizes) in both classes, for the BDI class we find a rigorous lower bound on $\chi$.  The explicit construction in <cit.> gives an upper bound on $\chi$ and in certain cases this upper bound coincides with the lower bound given here, proving that this is the optimal bond dimension in these cases.

We now expand on the above claims. As we review in detail below, any (finite-range) Hamiltonian in theses classes is equivalent to a Laurent polynomial $f(z) = \sum_\alpha t_\alpha z^\alpha$, (where the coupling constants $t_\alpha \in \mathbb{R}$ for BDI and $t_\alpha \in \mathbb{C}$ for AIII). This Laurent polynomial is characterised by the degree of the pole (or zero) at $z=0$ and a finite number of zeros away from $z=0$ . If no zero is on the unit circle, then the corresponding Hamiltonian is gapped and the ground state phase diagram is labelled by a winding number <cit.>. Physically, a non-zero winding number is topologically non-trivial, as evidenced—for instance—by topologically-protected zero-energy edge modes <cit.>.
Each gapped phase has a string order parameter <cit.>. These are non-local string correlation functions defined as a ground state expectation value of a product of an extensive number ($\Theta(N)$)
of fermionic operators. In the phase with winding number $\omega$, the corresponding string order parameter is non-zero in the limit of large $N$, while the string correlation functions that are order parameters for other phases decay with $N$. The subclass of models that we consider in this work is the case where every zero of $f(z)$ away from $z=0$ has even multiplicity. An exact formula for all string correlation functions for each model in this subclass constitutes our first main result.

For free-fermion systems, the entanglement spectrum of a subsystem is simply calculated from the eigenvalues of the correlation matrix—a matrix with elements consisting of two-point fermionic correlators <cit.>. From these eigenvalues one can easily calculate the entanglement entropy. Our second main result gives a method for finding a closed formula for the characteristic polynomial of the correlation matrix for a subsystem of size $N$, with explicit $N$ dependence. As mentioned, this result allows us to prove the existence of an exact MPS representation of the ground state. We also use this to give the characteristic polynomial in a series of examples.

Given an MPS, a useful construction in analysing correlations is the transfer matrix. Our results on string correlation functions allow us to derive properties of this transfer matrix, without needing the MPS itself. We focus on the BDI class; then, given a property that is generically satisfied in our class of models, we use the transfer matrix to give a lower bound on the bond dimension of any MPS representation of the ground state. The proof of this lower bound is through identifying eigenvalues of the transfer matrix that appear in correlation functions. As mentioned, in certain cases this lower bound coincides with an upper bound derived in <cit.>, giving the optimal bond dimension of the exact MPS ground state. In these cases, given the upper bound, we can find the entire spectrum of the transfer matrix—the eigenvalues are products of zeros and inverse zeros of $f(z)$ (multiplied by a sign that we determine). We also show how our results constrain the eigenvectors of the transfer matrix.

The class of models that we consider is a class of exceptional cases in the full class of BDI or AIII Hamiltonians. However, in this work we show that any (finite-range) model in these classes can arise as a limit of Hamiltonians considered in this paper. We remark that as we take the limit, we allow longer and longer range couplings. We illustrate how taking the limit of the string order parameter recovers the result for general models obtained in <cit.>. The concurrent work <cit.> discusses a concrete example where the transverse-field Ising chain is approximated by a sequence of models in this subclass.

The key object determining ground state correlations for these models is the symbol, or generating function, $\sqrt{f(z)/\overline{f}(1/z)}$ (note: $\overline{f}(z) = \sum_\alpha \overline{t_\alpha} z^\alpha$). In particular, Fourier coefficients of this function give the fermionic two-point function, and other correlations follow from this through Wick's theorem. In the general case, this is a multivalued function that we can analyse with branch cuts in the complex plane. For example, we can calculate asymptotic Fourier coefficients using the Darboux principle <cit.>. In our case, since every zero has even multiplicity, we have a rational function. This allows us to find a closed form for the Fourier coefficients (this observation has been used for identifying the correlation length of Gaussian MPS in <cit.>). Moreover, our result on string correlation functions is an application of Day's formula for Toeplitz determinants with rational symbol[This means that the matrix elements are Fourier coefficients of a rational function, in this case related to $\sqrt{f(z)/\overline{f}(1/z)}$. We review this terminology and the relevant theory in Section <ref>.] <cit.> and our result on the correlation matrix is an application of Gorodetsky's formula for block Toeplitz determinants with rational symbol <cit.>. Toeplitz determinants have long been associated with problems in statistical physics (see, for example, the review <cit.>), and Toeplitz determinants with rational symbols have previously appeared in problems in statistical mechanics and quantum many body physics <cit.>. However, to the best of our knowledge, this is the first application of Day's formula and Gorodetsky's formula in this context.

We note that similar conditions on particular generating functions being rational have previously appeared in works on Gaussian MPS—Gaussian states with finite bond dimension <cit.>. One characterisation of these Gaussian MPS is that the correlation matrix in Fourier space is a rational function of $\rme^{\rmi k_j}$, for momenta $k_1,\dots,k_d$ (this result applies in any spatial dimension $d$). Hence, in the models we consider, the correlation matrix having a rational symbol is a necessary condition that we have an MPS ground state. Our results show that in the one-dimensional BDI and AIII classes this is also sufficient.

The paper is organised as follows. In Section <ref> we introduce the BDI class in detail. For each gapped phase of the model, we have a string order parameter that can be viewed as a correlation function of two string operators. Section <ref> gives our results for the value of the correlation function of all such string operators in any given model in the given subclass. Section <ref> then gives our results for the characteristic polynomial of the correlation matrix of our model. We give examples of how to use this to find the eigenvalues of the correlation matrix. In Section <ref> we introduce the AIII class and explain how the main results of Section <ref> and Section <ref> apply to a subclass of these models. In Section <ref> we discuss how our results allow us to understand properties of the ground state transfer matrix in class BDI. Then
Section <ref> shows how to obtain general models in both classes as limits of models studied in this work; we moreover demonstrate how this can be used to obtain the formula for the order parameter in general BDI models. The remainder of the paper is dedicated to proving these results. We introduce the needed theory of Toeplitz determinants in Section <ref>. In Section <ref> we prove the results of Section <ref> and in Section <ref> we prove the results of Section <ref>.

§ THE BDI CLASS

Consider a periodic one-dimensional chain with $L$-sites, where for each site we have a spinless fermionic degree of freedom $c^\dagger_n$ (for $n=1\dots L$) that satisfies $\{ c_n, c_m \}= 0$ and $\{c_n^\dagger ,c_m \} = \delta_{nm}$. We will analyse this model in the thermodynamic limit $L\rightarrow \infty$. Define the Majorana operators
\begin{align}
\gamma_n = c^\dagger_n + c_n \qquad \tilde \gamma_n = \rmi (c^\dagger_n -c_n),
\end{align}

which satisfy $\{ \gamma_n,\gamma_m \} = \{ \tilde \gamma_n, \tilde \gamma_m \} = 2 \delta_{nm}$ and $\{ \gamma_n, \tilde \gamma_m \} = 0$.
The class of BDI Hamiltonians is defined as the vector space of all Hamiltonians that are quadratic in fermionic operators and that are symmetric under the anti-unitary involution $T \gamma_n T = \gamma_n$ and $T \tilde \gamma_n T = -\tilde \gamma_n$ <cit.>.
The most general translation-invariant one-dimensional BDI Hamiltonian is of the form:

\begin{align}
H = \frac{\rmi}{2}\sum_\alpha \sum_{n \in \mathrm{sites}}t_\alpha \tilde \gamma_n \gamma_{n+\alpha}, \qquad t_\alpha \in \mathbb{R}. \label{HBDI}
\end{align}

This model was first analysed in <cit.>. We will work in settings where $t_\alpha$ is only non-zero for a finite number of $\alpha \in \mathbb Z$ (physically, this corresponds to local Hamiltonians).
Using the Jordan-Wigner transformation, the above Hamiltonian is equivalent to a spin-12 chain—details are given in Appendix <ref>.
As discussed above, it is helpful <cit.> to introduce the complex function

\begin{align}
f(z) = \sum_\alpha t_\alpha z^\alpha. \label{fz}
\end{align}

The Hamiltonian (<ref>) is diagonalised as follows. For each momentum $k \in [0,2\pi)$, $\lvert f(\rme^{\rmi k})\rvert$ is the one-particle energy of the mode labelled by $k$, while the mode itself is given by

\begin{align}
\eta_k = \frac{1}{2}\left(1 +\frac{ {f(\rme^{\rmi k})}}{|f(\rme^{\rmi k})|} \right)c^\dagger_{k}+ \frac{1}{2}\left(1 -\frac{ {f(\rme^{\rmi k})}}{|f(\rme^{\rmi k})|} \right)c_{-k},\label{modes}
\end{align}

where $c_{k}$ is the Fourier transform of $c_n$.
The eigenstates of the many-body Hamiltonian then corresponding to filling these modes. Other properties of the model follow from knowing the zeros of this function. For a finite-range Hamiltonian, $f(z)$ is a Laurent polynomial: asymptotic correlations for that case were analysed in <cit.> and asymptotics of the characteristic polynomial of the correlation matrix and the entanglement entropy were analysed for gapped Hamiltonians in <cit.> and for certain[Subject to the restriction that $f(z) = f(1/z)$—physically this means an additional symmetry under spatial inversion.] gapless Hamiltonians in <cit.>.

If the Hamiltonian corresponding to (<ref>) is gapped, then the ratio appearing in (<ref>) satisfies:

\begin{align}
\frac{ {f(z)}}{|f(z)|}=\sqrt{\frac{f(z)}{f(1/z)}}. \label{sqrtfz}
\end{align}

In this paper we are interested in continuous families of `exceptional cases' for this model. In particular, we restrict to the case where (<ref>) is a rational function. This means that we can take[For completeness, let $h(z)$ denote any Laurent polynomial that satisfies $h(z)=h(1/z)$, has no zeros on the unit circle and has positive constant term. Then we can multiply (<ref>) by $h(z)$ without changing the ground state—see the discussion in <cit.>, as well as comments in Remark <ref>.] every zero to have even multiplicity, i.e., $f(z)$ has the following form:

\begin{align}
f(z) = \sigma \frac{1}{z^{n_P}} \prod_{j=1}^{n_z} \left(z-z_j\right)^2\prod_{k=1}^{n_Z} \left(z-Z_k\right)^2, \qquad\qquad \sigma \in \{\pm 1\} \label{fzcanon}
\end{align}

where $n_P\in \mathbb{Z}$ and we have that $0<\lvert z_j \rvert <1 <\lvert Z_k \rvert$. (In Section <ref> we will see that general models in the BDI class can be obtained from limits of these exceptional cases.) We have implicitly normalised $f(z)$, or equivalently the Hamiltonian, since this overall normalisation does not affect the ground state. Note that since $t_\alpha\in\mathbb{R}$, zeros of $f(z)$ are either real, or come in complex-conjugate pairs. Given this form, we will assume the generic case, that  $\{z_{j_1}^{\vphantom {-1}},z_{j_2}^{-1},Z_{k_1}^{\vphantom {-1}},Z_{k_2}^{-1}\}_{j_1,j_2=1\dots n_z; k_1,k_2=1\dots n_Z}$ are pairwise distinct.
Assuming this generic case is mainly for ease of presentation. In fact, our results on the correlation matrix (implying the existence of an exact MPS ground state) can be applied directly to the non-generic case, while other results can be generalised to non-generic cases by taking an appropriate limit, this is discussed in Remark <ref>. It is also helpful to index each set of zeros by proximity to the unit circle, i.e. $\lvert Z_1\rvert \leq\lvert Z_2\rvert \leq \dots \leq \lvert Z_{n_Z}\rvert$ and $\lvert z_1\rvert \geq \lvert z_2 \rvert \dots \geq \lvert z_{n_z}\rvert$.
Unless explicitly mentioned, we fix $\sigma=1$—see Remark <ref> for the case $\sigma=-1$. Since this model has no zeros on the unit circle, it is gapped[Our analysis can be applied to gapless models where the multiplicity of zeros on the unit circle is even. As discussed in Remark <ref> and <cit.> these models have the same ground state as a gapped model in our class.] and the gapped phase diagram is labelled by a winding number, given by
\begin{align}
\omega =2n_z-n_P.
\end{align}

Ground state correlations are determined, using Wick's theorem, by the fermionic two-point correlation function <cit.>. This is calculated in Appendix <ref> for the class of models given by (<ref>).

§ RESULTS: CORRELATIONS OF STRING OPERATORS IN BDI

The fermionic string operators that form the string order parameter for the phase $\omega = \alpha$ are given by[The phase factor ensures that $\mathcal{O}_\alpha=\mathcal{O}_\alpha^\dagger$. Note that there is an erroneous absolute value in the corresponding term in <cit.>. Other formulae in that paper are unchanged given the correct definition (<ref>).  ] <cit.>:

\begin{align}\label{Odef}
\mathcal{O}_\alpha(n)=\begin{cases} \rmi^{\lfloor \alpha/2 \rfloor} \left(\prod_{m=1}^{n-1} \rmi\tilde\gamma_m\gamma_m\right)\gamma_{n}\gamma_{n+1}\dots\gamma_{n+\alpha-1} \quad& \alpha>0 \\
(-\rmi)^{\lfloor \alpha/2 \rfloor}\left(\prod_{m=1}^{n-1} \rmi\tilde\gamma_m\gamma_m\right)  (-\rmi\tilde\gamma_{n})\dots(-\rmi\tilde\gamma_{n+\lvert\alpha\rvert-1}) \quad &\alpha<0\\
\prod_{m=1}^{n-1} \rmi\tilde\gamma_m\gamma_m \quad &\alpha=0.\end{cases}
\end{align}

The order parameter in the phase $\omega$ is $\lim_{N\rightarrow\infty}\lvert\langle\mathcal{O}_\omega(1)\mathcal{O}_\omega(N+1)\rangle \rvert ,$ where angle brackets denote the ground state expectation value.
Our first result concerns the behaviour of these correlation functions for finite $N$. The case with $n_z=0$ and $n_Z=0$ is trivial: in that case $f(z) = 1/z^{n_P}$ and the string correlators satisfy
\begin{align}
\langle\mathcal{O}_\alpha(1) \mathcal{O}_\alpha(N+1)\rangle = \delta_{\alpha,-n_P}(-1)^{N(n_P-1)}.\end{align}

For all other cases we have an exact formula in terms of the zeros of $f(z)$. In fact, the correlators are zero for $N$ sufficiently large, apart from the cases $\omega - n_z\leq \alpha \leq \omega +n_Z$. For these cases, the formula is a linear combination of terms $r_M^N$, where each $r_M$ is a product of zeros and inverse zeros. The precise statement is quite involved, so let us first give as an example the case $f(z) = (z-a)^2(z-b)^2$, where either $a,b \in \mathbb{R}$ or $a =\overline{b}$.

Let $f(z) = (z-a)^2(z-b)^2$ where $\lvert a\rvert<1$ and $\lvert b\rvert>1$. Since $n_z =n_Z= 1$, we are in the phase $\omega =2$, and the only correlators that are not exactly zero (for sufficiently large $N$) are those with $1\leq \alpha \leq 3$. We have that:

\begin{align}
\langle\mathcal{O}_1(1) \mathcal{O}_1(N+1)\rangle  &= (-1)^N \frac{\left(b^2-1\right) (a b-1) }{a b^2 (b-a)}a^{N} \qquad &N\geq 1\\
\langle\mathcal{O}_2(1) \mathcal{O}_2(N+1)\rangle  &= (-1)^N \bigg(\underbrace{\frac{\left(1-a^2\right) (1-1/b^{2}) }{(1-a/b)^2}}_{\textrm{Order parameter}} + \frac{(1-ab)^2}{(b-a)^2}\left(\frac{a}{b}\right)^N \bigg)\quad &N\geq 1 \label{exampleO2}\\
\langle\mathcal{O}_3(1) \mathcal{O}_3(N+1)\rangle  &= (-1)^N \frac{\left(1-a^2\right) (1-a b) b}{(b-a)}b^{-N} \qquad &N\geq 1.
\end{align}

Let $f(z) = (z-a)^2(z-b)^2$ where $\lvert a\rvert<1$ and $\lvert b\rvert<1$. Since $n_z = 2$, $n_Z= 0$, we are in the phase $\omega =4$, and the only correlators that are not exactly zero (for sufficiently large $N$) are those with $2\leq \alpha \leq 4$. We have that:

\begin{align}
\langle\mathcal{O}_2(1) \mathcal{O}_2(N+1)\rangle  &= (-1)^N \left(ab\right)^{N} \qquad &N\geq 1 \label{O2inside}\\
\langle\mathcal{O}_3(1) \mathcal{O}_3(N+1)\rangle  &= (-1)^N \frac{1}{b-a}\bigg(b (1-a^2)(1-ab)b^N -a (1-b^2)(1-ab)a^N \bigg) \qquad &N\geq 1 \label{O3inside}\\
\langle\mathcal{O}_4(1) \mathcal{O}_4(N+1)\rangle  &= (-1)^N \underbrace{{\left(1-a^2\right) (1-b^2) (1-ab)^2}}_{\textrm{Order parameter}} \qquad &N\geq 2 \label{O4inside} .
\end{align}

Let $f(z) = (z-a)^2(z-b)^2$ where $\lvert a\rvert>1$ and $\lvert b\rvert>1$. Since $n_z = 0$, $n_Z= 2$, we are in the phase $\omega =0$, and the only non-zero correlators (for sufficiently large $N$) are those with $0\leq \alpha \leq 2$. We have that $\langle\mathcal{O}_0(1) \mathcal{O}_0(N+1)\rangle $ is given by (<ref>), $\langle\mathcal{O}_1(1) \mathcal{O}_1(N+1)\rangle$ is given by (<ref>) and $\langle\mathcal{O}_2(1) \mathcal{O}_2(N+1)\rangle$ is given by (<ref>) under the replacements $a\rightarrow1/a$, $b\rightarrow 1/b$.

The correlations take the form of a sum of terms where each term is a coefficient that is a rational function of zeros and inverse zeros of $f(z)$, multiplied by $r_M^N$, where $r_M$ is a product of zeros from inside the unit circle and inverse zeros from outside unit the circle. This term $r_M$ is labelled by a set $M$ that tells you which zeros and inverse zeros contribute.
The order parameter term appears for the particular set where $r_M=1$. Informally, for values of $\alpha$ outside the range of $\alpha$ given, we `run out' of zeros to include in the product, and so the correlator is zero.
Notice that the formulae given are symmetric functions of zeros inside the unit circle and, separately, the zeros outside the unit circle. Furthermore, we are able to define the non-trivial correlators for $\lvert a\rvert,\lvert b\rvert>1$ in terms of the non-trivial correlators for $\lvert a\rvert,\lvert b\rvert<1$, under the transformation $a\rightarrow1/a$, $b\rightarrow 1/b$. Moreover, for $\lvert a\rvert <1$ and $\lvert b\rvert>1$, the formula for $\langle\mathcal{O}_1(1) \mathcal{O}_1(N+1)\rangle$ is identical to the formula $\langle\mathcal{O}_3(1) \mathcal{O}_3(N+1)\rangle$ under the replacements $a\rightarrow 1/b$ and $b\rightarrow 1/a$. These relationships apply more generally, due to an identity between the $ \langle\mathcal{O}_\alpha(1) \mathcal{O}_\alpha(N+1)\rangle$ in the chain corresponding to $f(z)$ and  $ \langle\mathcal{O}_{-\alpha}(1) \mathcal{O}_{-\alpha}(N+1)\rangle$ in the chain corresponding to $f(1/z)$. We will see in the proof that depending on the sign of $n_z+n_Z-n_P-\alpha$ the calculation uses  either $f(z)$ or $f(1/z)$ (this inverts the zeros), giving the two cases in the following general result:

Consider a chain with a non-trivial $f(z)$ as defined in (<ref>), assuming the generic case, and for each $\alpha \in \mathbb{Z}$ define $N_{\alpha} = \max\{\lvert n_z+n_Z-n_P-\alpha \rvert,1\}$. Then:

* for $\alpha < \omega - n_z$ and $\alpha > \omega +n_Z$

\begin{align}
\langle\mathcal{O}_\alpha(1) \mathcal{O}_\alpha(N+1)\rangle  &= 0 \qquad N\geq N_{\alpha};
\end{align}

* for $\omega - n_z\leq \alpha \leq \omega +n_Z$
\begin{align}
\langle\mathcal{O}_\alpha(1) \mathcal{O}_\alpha(N+1)\rangle  &= (-1)^{N(n_P+1)}\sum_{M} C_M^{\vphantom N} r_M^N \qquad N\geq N_{\alpha},\label{eq:thm1}
\end{align}

where for $n_z>0$ and $n_Z>0$:
* for $\omega - n_z\leq \alpha\leq \omega +n_Z-n_z$ the sum is over sets $M$ that label sets of zeros, with the non-zero constants $C_M$ and $r_M$ defined in Case 1 below.

* for $\omega +n_Z-n_z \leq \alpha\leq \omega +n_Z$ the sum is over sets $M$ that label sets of inverse zeros, with the non-zero constants $C_M$ and $r_M$ defined in Case 2 below.
If either $n_z=0$ or $n_Z=0$:

* if $n_z =0$, and $\alpha<\omega+n_Z$, then put $n_z=1$, $\alpha \rightarrow \alpha+2$ and take the limit $z_1 \rightarrow 0$ with Case 1. If $n_z=0$ and $\alpha = \omega+n_Z$ then the correlator is given by $(-1)^{N(\omega-1)} \prod_{k=1}^{n_Z} Z_k^{-N}$.

* if $n_Z =0$ and $\alpha>\omega-n_z$, then put $n_Z=1$, and take the limit $Z_1 \rightarrow \infty$ with Case 2. If $n_Z =0$ and $\alpha=\omega-n_z$, then the correlator is given by: $(-1)^{N(\omega-1)} \prod_{j=1}^{n_z} z_j^N$.

In both cases, in the formula (<ref>), the sum is over all subsets $M \subseteq \{1,\dots,n_z+n_Z\}$ of some fixed size $\lvert M \rvert$. For each such $M$, we also define:

\begin{align}
M^c &=  \{1,\dots,n_z+n_Z\}\setminus M \label{case1eq1} \end{align}
and further define $M_z = \{1,\dots, n_z\}$ and $M_Z= \{1,\dots, n_Z\}$. Then we have:

Case 1:
Define $\tau_j = z_j$ for $1\leq j \leq n_z$ and $\tau_{n_z+j}=Z_j$ for  $1\leq j \leq n_Z$.
Then the sum in (<ref>) is over all subsets $M$ of size $\lvert M \rvert =\alpha+n_P-n_z$, and $C_M$ and $r_M$ are defined by:

\begin{align}
r_M&=  \frac{\prod_{k \in M^c} \tau_k} {\prod_{j=1}^{n_Z}Z_j},\nonumber\\\nonumber
C_M &= \prod_{k \in M^c,m\in M_Z} (\tau_k - Z^{-1}_m) \prod_{l \in M_z, j\in M} (z_l^{-1} - \tau_j) \prod_{l\in M_z}z_l^{-(n_z+n_Z -n_P-\alpha)} \\&\times\prod_{l \in M_z, m\in M_Z} (z_l^{-1} - Z^{-1}_m)^{-1} \prod_{k \in M^c,j\in M} (\tau_k - \tau_j)^{-1}  \prod_{k\in M^c}\tau_k^{-(n_z+n_Z -n_P-\alpha)}\label{case1eq2} .
\end{align}

Case 2:
Define $\tau_j = z_j^{-1}$ for $1\leq j \leq n_z$ and $\tau_{n_z+j}=Z_j^{-1}$ for  $1\leq j \leq n_Z$.
Then the sum in (<ref>) is over all subsets $M$ of size $\lvert M \rvert =2n_z+n_Z-n_P-\alpha$, and $C_M$ and $r_M$ are defined by:
\begin{align}
r_M&=  \prod_{j=1}^{n_z}z_j{\prod_{k \in M^c} \tau_k} ,\nonumber\\
C_M &= \prod_{k \in M^c,m\in M_z} (\tau_k - z_m) \prod_{l \in M_Z, j\in M} (Z_l - \tau_j)\prod_{l \in M_Z} Z_l^{n_P+\alpha-n_z-n_Z}\nonumber\\&\times \prod_{l \in M_Z, m\in M_z} (Z_l - z_m)^{-1} \prod_{k \in M^c,j\in M} (\tau_k - \tau_j)^{-1} \prod_{k \in M^c}\tau_k^{-(n_P+\alpha-n_z-n_Z)}. \label{case2eq2}
\end{align}

Note that in both Case 1 and Case 2, we can have $\lvert M \rvert =0$ or $\lvert M^c \rvert =0$. For $\lvert M \rvert =0$, we define double products with one variable in $M$ that appear in the above formulae, such as $\prod_{l \in M_z, j\in M} (z_l^{-1} - \tau_j)$, to be equal to one. Similarly double products with one variable in $\lvert M^c \rvert =0$ are equal to one.  Note also that $\alpha=\omega +n_Z-n_z$ can be evaluated using both Case 1 and Case 2.

Theorem <ref> is proved by considering a particular limiting case of Day's formula for Toeplitz determinants (Theorem <ref>). The limit in question arises by writing $z^k = \lim_{\eps\rightarrow 0} \prod_{j=1}^k (z-\eps x_j)$, for $x_j$ pairwise distinct, so that we can apply Day's formula. The limit of the final formula simplifies for $N\geq N_{\alpha}$ leading to the result above. However, for any $N$ such that $1\leq N<N_{\alpha}$ the steps outlined in the proof can be used to evaluate the correlator.
One evaluation of this limit for $N<N_{\alpha}$ is the following: for $n_z\neq0$ and $n_Z=0$ we have that for any $\nu>0$:

\begin{align}
\langle \mathcal{O}_{\omega+\nu}(1) \mathcal{O}_{\omega+\nu}(N+1) \rangle &=\begin{cases} (-1)^{n_z(\omega-1)}\left(\prod_{j=1}^{n_z} z_j^\nu\right)  \prod_{j_1=1}^{n_z} \prod_{j_2=1}^{n_z} \left(1-z_{j_1}z_{j_2}\right) \qquad& N=  n_z \\
0 \qquad & N>n_z\end{cases}.
\end{align}

This follows from the results of Section 7 of <cit.>. Note that in this case, Theorem <ref> gives that this correlator is zero for $N\geq n_z + \nu$. There is an analogous formula for the case $n_Z\neq 0$ and $n_z = 0$.

For $\alpha = \omega$, there is always one choice of $M$ where $r_M=1$. For this choice of $M$, $C_M$ is the order parameter. We can evaluate this $C_M$ using Theorem <ref> to reach:

\begin{align}
\lim_{N\rightarrow\infty}\lvert\langle\mathcal{O}_\omega(1)\mathcal{O}_\omega(N+1)\rangle \rvert = \frac{ \prod_{j_1,j_2=1}^{n_z} (1-z_{j_1} z_{j_2})
\prod_{k_1,k_2=1}^{n_Z} \left(1-\frac{1}{Z_{k_1} Z_{k_2}}\right)}{\prod_{j=1}^{n_z}\prod_{k=1}^{n_Z} \left (1-\frac{z_j}{ Z_k}\right)^2}. \label{eq:orderparameter}
\end{align}

This formula also holds in the non-generic case—see Remark <ref>. Note that this result, following from Day's formula, agrees with the general formula for the order parameter given in <cit.> as applied to our case. The general formula is proved using Szegő's theorem <cit.>. In fact, by taking a limit of models, one can use (<ref>) to recover the result of <cit.> for the general case (see Section <ref>).

For $\alpha = \omega +n_Z$ and $\alpha = \omega-n_z$ we see that there is only one choice of $M$, and hence the correlator has a particularly simple form.
Indeed, one has:

\begin{align}
\langle\mathcal{O}_{\omega +n_Z}(1) \mathcal{O}_{\omega +n_Z}(N+1)\rangle  &= (-1)^{N(n_P+1)} C_M \left(\prod_{k=1}^{n_Z} Z_k \right)^{-N}  &N \geq \max\{n_z,1\}~ \label{eq:exampleforasymptotics}\\
\langle\mathcal{O}_{\omega-n_z}(1) \mathcal{O}_{ \omega-n_z}(N+1)\rangle  &= (-1)^{N(n_P+1)} C_M \left(\prod_{j=1}^{n_z} z_j \right)^N  &N \geq \max\{n_Z,1\}.
\end{align}
For $n_z=0$ or $n_Z=0$ these formulae apply and this means that the order parameter reaches its limiting value for $N\geq N_{\alpha}$ with no correction term. This property was emphasised for the case $f(z) = (z-a)^2$ in <cit.>. While the result for $n_z=0$ or $n_Z=0$ follows from our analysis, it may also be derived as an application of the results of Section 7 of <cit.>.

The large $N$ asymptotics for each correlator can be deduced easily from Theorem <ref>. To find the dominant term, one chooses the subset $M$ such that the corresponding $r_M$ is as large as possible. Note that choosing $M$ of size $m$ is equivalent to choosing $M^c$ of size $n_z+n_Z-m$. We maximise $r_M$ by choosing the $n_z+n_Z-m$ largest zeros (or inverse zeros) to make up $M^c$. Concretely, for a correlator in Case 1 we choose $M^c = \{Z_{n_Z},Z_{n_Z-1}, \dots \}$ of the appropriate size. Since all zeros are either real or come in complex-conjugate pairs, this is not necessarily a unique maximum. If we find a dominant $r_M$ as described and it contains a complex zero but not its conjugate, then there will be a corresponding $M'$ with the zero replaced by its conjugate such that $\lvert r_M \rvert = \lvert r_M' \rvert$. In general, we could have more than two zeros of the same absolute value. The dominant asymptotic term will then be a sum over contributions from all $M$ with the maximal $r_M$. The result is in agreement with Theorem 3 of <cit.> (where it is assumed that there is no degeneracy of the smallest zero that appears in $r_M$ and so we have a single dominant contribution).

The second most dominant term will be where we take the dominant $M^c$, remove the smallest zero that appears, and replace it by the largest zero that did not feature in the dominant term (up to accounting for zeros of the same absolute value as above). In general, by ordering the size of $\prod_{k \in M^c}\tau_k$, one can find as many terms in the large $N$ asymptotics as desired.

We can define the correlation length, $\xi_\alpha$ for the operators $\mathcal{O}_\alpha$ with $\alpha \neq \omega$ by:

\begin{align}
\langle\mathcal{O}_{\alpha}(1) \mathcal{O}_{\alpha}(N+1)\rangle = O^*(\rme^{-N/\xi_\alpha})
\end{align}
as $N\rightarrow\infty$. Following <cit.>, we write $f(x)=O^*(g(x))$ if $g(x)$ is a tight upper bound for $f(x)$; i.e., $f(x)=O(g(x))$ but $f(x)\neq o(g(x))$. This allows for oscillatory cases where $\xi$ gives the exponentially decaying envelope, but $f(x)$ may vanish for some values of $x$.
Theorem <ref> and Remark <ref> lead to:

Consider a chain with a non-trivial $f(z)$ as defined in (<ref>). Then the correlation length is given by:

\begin{align}
\xi_\alpha^{-1} = \begin{cases}
\sum_{j=1}^{\omega-\alpha} \log\lvert z_j^{-1}\rvert \qquad & \omega-n_z\leq\alpha<\omega\\
\sum_{j=1}^{\alpha-\omega} \log\vert Z_j\rvert \qquad & \omega<\alpha\leq\omega+n_Z\\
\infty \qquad & \alpha<\omega-n_z\quad \mathrm{or} \quad \alpha>\omega+n_Z.
\end{cases}
\end{align}

In <cit.> an expression for the large $N$ asymptotics of each correlator is given. The models we consider here arise as limiting cases of models in that paper, so it is interesting to compare our results. We firstly point out that the method used in <cit.> for evaluating the asymptotics requires the evaluation of the determinant of an $\lvert\omega-\alpha\rvert\times \lvert\omega-\alpha\rvert$ matrix, where the matrix elements are, up to some error terms, the asymptotically large Fourier coefficients of functions $l(z)$ and $m(z)$ that are defined in Appendix <ref> <cit.>. Analysing this determinant for the cases studied here is non-trivial, except when $\lvert\omega-\alpha\rvert=1$. Firstly, let us suppose that $n_Z \geq 1$ and $n_z \geq 1$. In Appendix <ref> we show that $l_N =O^*(\lvert Z_1\rvert ^{-N})$ and $m_N =O^*(\lvert z_1\rvert ^{N})$. Then we have that for $\omega-\alpha=1$, the correlator behaves like $m_N$ and for $\omega-\alpha=-1$ the correlator behaves like $l_N$. This agrees with our results here. For example, suppose the order parameter can be calculated using Case 1. Then for $\alpha = \omega$ the dominant term comes from  $M^c = \{Z_{n_Z},Z_{n_Z-1}, \dots Z_1\}$. For $\alpha = \omega-1$, $\lvert M \rvert$ is one smaller, so  $\lvert M^c \rvert$ is one larger and so dominant terms come from $ M^c = \{Z_{n_Z},Z_{n_Z-1}, \dots Z_1,z_j\}$ where $\lvert z_j\rvert =z_1$—each of these sets has $\lvert r_M\rvert =\lvert z_1\rvert$, and hence the asymptotic behaviour $O^*(\lvert r_M\rvert^N)$ agrees with $m_N$. We give an explicit example in Appendix <ref>, showing that the terms match exactly.

An important difference with the results of <cit.> is that, in the more general class considered there, generically in that class the correlation length of each string operator is given by $\xi_\alpha^{-1} = \xi^{-1}\lvert \omega-\alpha \rvert$; where the basic correlation length $\xi$ is given by
\begin{align} \xi^{-1} = \min\{\log\lvert z_1^{-1}\rvert,\log\lvert Z_1\rvert). \end{align}

This means that the decay of these correlators depends only on the zero closest to the unit circle. By comparing this result with Corollary <ref>, for $\lvert \omega-\alpha \rvert>1$ we see that generically in the class of models given by (<ref>) we have faster decay than in the generic class of models considered in <cit.>.

Recall that the one-particle energies of our model are given by $\lvert f(\rme^{\rmi k})\rvert$. Allowing for a non-zero temperature, $T$, the following equation appears in studying zeros of the partition function <cit.> (often one considers a large but finite system, so that $k$ is quantised):

\begin{align}
\lvert f(\rme^{\rmi k})\rvert = \rmi(2n+1)\pi T \qquad n\in\mathbb{Z}. \label{eq:partitionfunction}
\end{align}

For example, the Fisher zeros correspond to solutions of this equation[More generally, for a system of size $L$, Fisher zeros are solutions of $Z(L,\beta)=0$. In our case, we have that $Z(L,\beta)=\prod_{k_m}(1+\rme^{-\beta \lvert f(k_m)\rvert})$ (where we have fixed the ground state to have zero energy), leading to (<ref>).] in the complex $\beta$ plane for $\beta = 1/T$ <cit.>. If we instead take this equation and complexify $z=\rme^{\rmi k}$, we see that at zero temperature the zeros of $f(z)$ and zeros of $f(1/z)$ coincide with the set of zeros of (<ref>); this method is used to study disorder lines in <cit.>. Note that to analytically continue the left-hand-side of (<ref>) we can write $\lvert f(z) \rvert = \sqrt{f(z)f(1/z)}$ on the unit circle.
The subclass of models considered in this paper then corresponds to degeneracy of solutions to (<ref>)—
it would be of interest to understand how our results generalise to the finite-temperature case.

Another important correlator is the emptiness formation probability <cit.>, given by:

\begin{align}
P(N)=\left\langle \prod_{j=1}^N c^{\vphantom \dagger}_j c^\dagger_j \right \rangle .
\end{align}

This is a string correlator for complex fermions, rather than the Majorana fermions considered above, and is the probability, in the ground state, that $N$ consecutive fermionic sites are unoccupied (equivalently, in the spin chain picture, that there are $N$ consecutive up-spins). A closed formula for $P(N)$ was found for the disorder line in the XY model in <cit.>, using the simple formula for the ground state in that case. We remark that closed formulae for $P(N)$ can be found for the class of models considered in this paper, since this correlator may also be evaluated using Day's formula. We give an example below, and in Section <ref> outline the general approach.

Let $f(z) = z^{-2} (z-a)^2(z-b)^2$, with $\lvert a\rvert<1$ and $\lvert b\rvert>1$. Then:

\begin{align}
P(N) = \left(\frac{b^{-1}-a}{2}\right)^N\left(\frac{(1-a)(1+b)}{2(b-a)}+ (-1)^N\frac{(1+a)(1-b)}{2(a-b)} \right) \qquad N\geq 1.
\end{align}

This is derived using Day's formula in Section <ref>.

§ RESULTS: CORRELATION MATRIX IN BDI

Let $a_{2n-1} = \gamma_n$ and $a_{2n} =\tilde\gamma_n$. The correlation matrix for a subsystem of size $N$ is defined by:

\begin{align}
A_N= \left(  \rmi \langle a_j a_k \rangle - \rmi\delta_{jk}\right)_{j,k=1\dots 2N} .
\end{align}

The matrix $A_N$ has $2N$ imaginary eigenvalues, $\{\pm\rmi\nu_1\dots \pm\rmi\nu_N\}$. $A_N$ always has an even number of eigenvalues that are equal to zero. If there are $2m$ such eigenvalues then $\nu_1,\dots, \nu_m=0$ and the remaining $\nu_j$ are non-zero. With a slight abuse of terminology to simplify our discussion, we will call the $N$ non-negative imaginary parts $\{\nu_1\dots \nu_N\}$ `the eigenvalues' of the correlation matrix (these eigenvalues are zeros of the characteristic polynomial $\det(\rmi\lambda-A_N)$).  Let $\rho_N$ be the density matrix of a subsystem of our chain of size $N$, the eigenvalues of $\rho_N$ can be computed from the eigenvalues of the correlation matrix. Indeed, the following formula (given in, for example, <cit.>) allows us to write the eigenvalues of the reduced density matrix in terms of the eigenvalues of the correlation matrix:

\begin{align}
\lambda_{x_1\dots x_N} =\prod_{j=1}^N \frac{1+(-1)^{x_i}\nu_j}{2}, \qquad x_j=0,1.
\end{align}

We go over all values $x_i$, so have $2^N$ eigenvalues of $\rho_N$ from the $N$ eigenvalues of the correlation matrix $\{\nu_1\dots \nu_N\}$. That the $\lambda$ are eigenvalues of a density matrix means that each $\nu_j \leq 1$.  An eigenvalue $\nu_j=1$ we call trivial. Moreover, the von Neumann entanglement entropy ($S(N)=-\tr(\rho_N\log{\rho_N})$), or more generally the Rényi entropy ($S_\alpha(N)=\frac{1}{1-\alpha}\log\tr(\rho_N^\alpha)$), has a simple expression in terms of the non-trivial eigenvalues $\nu_j$ <cit.>.
In our class of interest, there are only a finite number of non-trivial eigenvalues as $N\rightarrow \infty$. For $f(z)=z^k$, then $\det(\rmi\lambda-A_N)=(-\lambda^{2})^{\lvert k\rvert}(1-\lambda^2)^{N-\lvert k\rvert}$. In all other cases, the claim follows from:

Consider a chain with a non-trivial $f(z)$ as defined in (<ref>) such that $n_P=n_z+n_Z$. We allow the non-generic case where zeros may coincide, but assume that $z_j\neq Z_k^{-1}$ for any $j,k$. Denote the correlation matrix on a subsystem of size $N$ by $A_N$. Then:

\begin{align}
\det(\rmi \lambda-A_N)= (1- \lambda^2)^N  \left(\frac{\prod_{j=1}^{n_z} z_j^2}{\prod_{k=1}^{n_Z} Z_k^2}\right)^N \frac{\det{M(N,\lambda)}}{\det{M(0,\lambda)}}, \label{thm2formula}\end{align}

where $M(n,\lambda)$ is a $4(n_z+n_Z)\times 4(n_z+n_Z)$ matrix that can be determined from the zeros of $f(z)$ and the Smith canonical form of a known $2\times 2$ matrix polynomial (defined in (<ref>)-(<ref>)). Moreover, there exist $d', d'' \in \mathbb{N}$ (independent of $N$) such that:

\begin{align}
\det{M(n,\lambda)} = \mu \left(\frac{\prod_{j=1}^{n_z} z_j^2}{\prod_{k=1}^{n_Z} Z_k^2}\right)^{-n} \prod_{j_1=1}^{d'}(\lambda^2-\tilde x_{j_1}(n))\prod_{j_2=1}^{d''}(\lambda-c_{j_2}) , \label{thm2formula2}\end{align}

where the constants $\mu$ and $c_j$ are independent of $n$ and $\tilde x_j(0)=1$, $\tilde x_j(n)\leq 1$ for $n\geq 1$ and $\tilde x_j(n)<1$ for some $n>0$.

For any $f(z)$ as defined in (<ref>) such that $n_P \neq n_z+n_Z$, the characteristic polynomial of the correlation matrix can be found by taking a limit of a related case where $n_P = n_z+n_Z$ that can be evaluated using the above. Details of this limit are discussed in Remark <ref>.

As we explain in Remark <ref>, the assumption that $f(z)$ cannot contain mutually inverse zeros can be trivially accounted for. Then we have:

Given any $f(z)$ as defined in (<ref>), with no assumptions of a generic case, Theorem <ref> leads to the following expression for the characteristic polynomial of the correlation matrix:

\begin{align}
\det(\rmi\lambda-A_N)= (1- \lambda^2)^{N-d}  \times \prod_{j=1}^{d}(\tilde x_j(N)-\lambda^2),\label{characteristicpolynomial}
\end{align}

where $d\in\mathbb{N}$ is independent of $N$.

By computing $\det{M(N,\lambda)}$ we can determine $d$, the number of non-trivial eigenvalues of the correlation matrix. This is illustrated in Example <ref> and Example <ref>.

Let $f(z) = \frac{1}{z}(z-b)^2$ with $\lvert b\rvert <1$. This corresponds to the aforementioned disorder line of the XY model <cit.>. We can evaluate the terms appearing in Theorem <ref> as follows. Using the definition of $M(N,\lambda)$ given in Section <ref> and the Smith canonical form for this example given in Appendix <ref> we have:

\begin{align}
\begin{array}{cccc}
\lambda(1 -b^2)   & 0 & 0 & 0 \\
\lambda  \left(\lambda ^2-3\right)b & \rmi \lambda ^2 & \lambda  b^N & 0 \\
0 & 0 & 0 & -\rmi  \lambda ^2\left(b^2-1\right)b^{-(N+2)} \\
0 & -\rmi \lambda ^2 & -\lambda ^3 b^{-N} & -\rmi \lambda ^2 b^{-(N+1)} \left(\left(b^2-1\right) N-2\right) \\
\end{array}
\right)\\
\det{M(N,\lambda)}&=\left(-b^{-2}\left(b^2-1\right)^{6} \lambda^{6}\right) b^{-2N}\left(\lambda^2 -b^{2N}\right)\\
\frac{\det{M(N,\lambda)}}{\det{M(0,\lambda)}}&=\frac{b^{-2N} \left(\lambda^2 -b^{2N}\right)}{\lambda^2-1}
\\
\det(\rmi\lambda-A_N)&= (1- \lambda^2)^{N-1}(b^{2N}-\lambda^2). \end{align}

This means that there is one non-trivial eigenvalue of the correlation matrix for a subsystem of size $N$: $\nu_1=\lvert b\rvert^{N}$. As $N\rightarrow \infty$ the limit gives one zero eigenvalue. This agrees with taking the limit as we approach the disorder line in the XY model—there the infinitely many correlation matrix eigenvalues in the limit $N\rightarrow \infty$ are known for all values of the couplings (away from the disorder line these models are outside the class analysed here) <cit.>.

For $f(z) = \frac{1}{z}(z-b)^2$ with $\lvert b\rvert >1$, there is one non-trivial eigenvalue for a subsystem of size $N$: $\nu_1=\lvert b\rvert^{-N}$. This follows from the result for $\lvert b\rvert <1$ and Remark <ref> below.

Let $f(z) = \frac{1}{z^2}(z-a)^2(z-b)^2$. First consider the case $\lvert a\rvert <1$ and $\lvert b\rvert <1$.  We can evaluate the terms appearing in Theorem <ref>, see also Appendix <ref>, to give:

\begin{align}
\frac{\det{M(N,\lambda)}}{\det{M(0,\lambda)}}&=\frac{a^{-2N} b^{-2N}}{(\lambda^2-1)^2}\nonumber\\
&\times\left(\lambda^4 + \lambda^2 \frac{a^N b^N}{(a-b)^2}\left( 2(1-a^2)(1-b^2)-(1-a b)^2\left(\left(\frac{a}{b}\right)^N +\left(\frac{b}{a}\right)^N\right)\right) +a^{2N}b^{2N}\right) \label{example3eq}.
\end{align}

Thus we have two non-trivial eigenvalues of the correlation matrix for a subsystem of size $N$, given by the zeros of (<ref>). As $N\rightarrow \infty$ the coefficients of $\lambda^2$ and $\lambda^0$ go to zero, so in this limit we have two non-trivial eigenvalues equal to zero.

Now consider the case $\lvert a\rvert <1$ and $\lvert b\rvert >1$.  We can evaluate the terms appearing in Theorem <ref> to give:

\begin{align}
\frac{\det{M(N,\lambda)}}{\det{M(0,\lambda)}}=\frac{a^{-2N} b^{2N}}{(\lambda^2-1)^2}\Bigg(\lambda^4 &+\frac{ \lambda^2 }{(b-a)^2}\left( 2(1-a^2)(1-b^2)-(1-a b)^2\left( a^{2N}+b^{-2N}\right)\right)\nonumber\\
& +\left( \frac{a^Nb^{-N}(1-ab)^2-(1-a^2)(1-b^2)}{(b-a)^2}\right)^2\Bigg).\label{example3eq2}
\end{align}

Taking the limit $N\rightarrow \infty$ gives us that the two non-trivial eigenvalues of the correlation matrix are degenerate and given by:

\begin{align}
\nu_1^2=\nu_2^2 = {\frac{(1-a^2)(1-1/b^2)}{(1-a/b)^2}}.
\end{align}

Note that these limiting eigenvalues are equal to each other and to the order parameter, given in (<ref>)—these points are discussed in Remarks <ref> and <ref>.

The final case $\lvert a\rvert >1$ and $\lvert b\rvert >1$ is equivalent to the first case, where we make the replacement $a\rightarrow 1/a$ and $b \rightarrow 1/b$ in all formulae. We explain this in general terms in Remark <ref>.

Corollary <ref> tells us that as we take a subsystem of size $N\rightarrow \infty$, there are a finite number of non-zero eigenvalues of the reduced density matrix. This in turn directly implies that the ground state for $f(z)$ of the form (<ref>) can be written as an MPS of fixed bond dimension. Except for a measure zero set of cases, an MPS is explicitly constructed in <cit.> using different methods. The results of that paper give us an upper bound on $d$, the number of eigenvalues that depend on $n$ in (<ref>). The bond dimension, $\chi$, of the optimal MPS representation of the ground state is related to $d$ by $2^d \leq \chi^2$. The upper bound is determined by the range of the Hamiltonian: $d \leq  2\log_2(\chi) \leq 2\lceil \mathrm{range}(H)/2\rceil$, where $\mathrm{range}(H)$ is the largest value of $\lvert \alpha \rvert$ such that $t_\alpha\neq0$. This holds in Examples <ref> and <ref>. It would be interesting to derive this bound in general using the methods of this work.

The correlation matrix is of block Toeplitz form, defined in Section <ref>.
The Szegő-Widom theorem for block Toeplitz determinants <cit.> tell us that as $N \rightarrow \infty$
\begin{align}
\det(\rmi\lambda-A_N) = (1- \lambda^2)^N E_W(1+o(1)),\end{align}

where $E_W$, Widom's constant, may be zero (see, for example, <cit.>).
This theorem is discussed in the context of the correlation matrix of a quantum chain in <cit.>.

Applying Theorem <ref> allows us to derive formulae for $E_W$ in our case. For example, for $f(z)=\frac{1}{z^2}(z-a)^2(z-b)^2$ with $\lvert a \rvert <1$ and $\lvert b\rvert >1$, by taking the limit of (<ref>) we have that $E_W$ is given by:

\begin{align}
\end{align}

We prove Theorem <ref> using Gorodetsky's formula (Theorem <ref>). Note that <cit.> gives a detailed discussion of asymptotics of determinants of block Toeplitz matrices that can be analysed using Gorodetsky's formula. As with our Theorem <ref>, to use these results to derive the non-trivial part of the asymptotics, $E_W$, we require an expression for the matrices in the Smith canonical form of the relevant matrix polynomial.

The models we consider are limiting cases of those studied in <cit.>. Another approach to finding the asymptotics would be to take a limit of the results of that paper.

Given $f(z)$ defined in (<ref>), and suppose that $n_P=n_z+n_Z-k$ for $k> 0$. We can write this as:

\begin{align}
f(z) &= \frac{z^k}{z^{n_z+n_Z}}\prod_{j=1}^{n_z} \left(z-z_j\right)^2\prod_{k'=1}^{n_Z} \left(z-Z_{k'}\right)^2\\
&=\lim_{\eps \rightarrow 0}\frac{1}{z^{n_z+n_Z+k}} \prod_{j=1}^{n_z+k} \left(z-z_j\right)^2\prod_{k'=1}^{n_Z} \left(z-Z_{k'}\right)^2
\end{align}

where $z_j = \eps x_j$ for $n_z+1 \leq j \leq n_z+k$. Before taking the limit $\eps \rightarrow 0$, we can define $n'_z=n_z+k$ and then apply Theorem <ref> with $n_P=n'_z+n_Z$. Since we are evaluating the finite determinant $\det(\rmi\lambda-A_N)$, taking the limit $\eps \rightarrow 0$ of this result gives us the characteristic polynomial of the correlation matrix for cases with $n_P=n_z+n_Z-k$.

We also have that the eigenvalues of the correlation matrix for the system corresponding to $f(z)$ are identical to the eigenvalues of the correlation matrix for the system corresponding to $f(1/z)$. This follows from the block Toeplitz form defined in Section <ref>. Moreover, for the purposes of calculating the eigenvalues of the correlation matrix, this transformation amounts to replacing all zeros by their inverse, and replacing
\begin{align}n_P \rightarrow 2(n_z+n_Z)-n_P.\end{align}

This allows us to analyse $n_P=n_z+n_z+k$ for $k>0$ as follows:

\begin{align}
f(1/z) &= \frac{z^k}{z^{n_z+n_Z}}\prod_{j=1}^{n_z} \left(z-1/z_j\right)^2\prod_{k'=1}^{n_Z} \left(z-1/Z_{k'}\right)^2.
\end{align}

Using the same reasoning as above, we then apply Theorem <ref> with $n_P= n'_z+n'_Z$, where $n'_z=n_Z+k$ and $n'_Z=n_z$, and then take the limit where $k$ zeros inside the unit circle go to zero.

We can use Remark <ref> and Example <ref> to find the correlation matrix eigenvalues for $f(z) = (z-b)^2$. For $\lvert b \rvert <1$ we take the limit $a\rightarrow 0$ in (<ref>). For $\lvert b \rvert >1$ we take the limit $a\rightarrow 0$ in (<ref>). Moreover, the correlation matrix eigenvalues for $f(z) =\frac{1}{z^2} (z-b)^2$ can be found by replacing $b \rightarrow 1/b$ appropriately in the above results.

In Example <ref>, on taking the subsystem size $N \rightarrow \infty$, we had two non-trivial eigenvalues such that $\nu_1=\nu_2$. Such a symmetry, where we have a twofold degeneracy of non-trivial correlation matrix eigenvalues, was observed for the infinitely many non-trivial $\nu_m$ in the XY model in <cit.>. In fact in that case $\nu_1$ is not degenerate, but for all $m>1$ we have $\nu_{2m}=\nu_{2m+1}$ in the limit. This has a straightforward physical origin, as pointed out in <cit.>.
Let us consider the spin chain picture, with the model given in Appendix <ref>. If the winding number is even, there is no symmetry breaking, and the correlation length is thus finite. Hence, for a large enough block, the two edges will decouple. Moreover, there is a symmetry between the two edges[E.g., note that bond-centered inversion swaps the two edges and maps $f(z) \to f(1/z)$; an on-site rotation over an angle $\pi$ also maps $f(z) \to f(1/z)$. In combination, we thus have a symmetry of the model which exchanges the two edges of the block.]. All correlation eigenvalues (in the limit $N \rightarrow \infty$) must thus come in degenerate pairs, corresponding to the correlations associated to each edge and the exterior of the block. For odd winding numbers, the ground state spontaneously breaks spin-flip symmetry. Hence, a symmetry-preserving ground state can be seen as a macroscopic `cat state' superposition of two ground states which obey the aforementioned doubled correlation matrix eigenvalues, resulting in one additional, non-degenerate correlation eigenvalue. This explains the observations in Example <ref> and <cit.>. On these physical grounds, we conjecture this property to hold much more generally, although we do not have a proof of this even for the class considered in this paper.

We have the following relationship between the correlation matrix and a string correlation function that may be evaluated with Theorem <ref>:

\begin{align}
\det(A_N) = \prod_{j=1}^N \nu_j^2=  \Big(\langle \mathcal{O}_0(1) \mathcal{O}_0(N+1) \rangle\Big)^2;\label{eq:detAN}
\end{align}

this follows from the block Toeplitz structure, see also (<ref>). This relationship holds for all $N$, and one can check, for example,
that the product of the zeros of (<ref>) is equal to the square of the right-hand-side[Example <ref> has $n_P=0$, while Example <ref> has $n_P=2$, so this is the relevant correlator.] of (<ref>). Moreover, suppose that we have two non-trivial correlation eigenvalues, and assume the symmetry argued in Remark <ref>, so that $\nu_1=\nu_2$ in the limit $N\rightarrow \infty$. Then
it follows immediately from (<ref>) that the limiting correlation eigenvalues are equal to the order parameter, as observed in Example <ref>. Equation (<ref>) also shows that, in the limit $N\rightarrow \infty$, the correlation matrix has at least one zero eigenvalue if and only if the winding number is non-zero. Equivalently, we have a degeneracy of the entanglement spectrum if and only if the trivial string order parameter $\mathcal O_0(N)$ does not have long-range order, consistent with the physical fingerprints associated to non-trivial SPT phases <cit.>.

Throughout we have assumed that $f(z)$ defined in (<ref>) has a positive overall sign. If the overall sign is negative, the correlators $\langle\mathcal{O}_\alpha(1) \mathcal{O}_\alpha(N+1)\rangle $ are multiplied by $(-1)^N$—this is shown in Section <ref>. Moreover, the correlation matrix eigenvalues are unchanged—see Section <ref>.

We define the generic case for (<ref>) to be where $\{z_{j_1}^{\vphantom {-1}},z_{j_2}^{-1},Z_{k_1}^{\vphantom {-1}},Z_{k_2}^{-1}\}_{j_1,j_2=1\dots n_z; k_1,k_2=1\dots n_Z}$ are pairwise distinct. If we allow mutually inverse zeros, i.e. some $z_j = Z_k^{-1}$, then it turns out that the ground state is equivalent to a related generic model: we simply remove all of these mutually inverse pairs of zeros and shift $n_P$. This is also the case when we allow zeros on the unit circle of even multiplicity. A full analysis of these cases can be found in the paper <cit.>. We can further relax the condition that the zeros are pairwise distinct by writing, say, $z_j=z_{j'}+\eps$ and taking the limit $\eps\rightarrow 0$, as in Remark <ref> and Remark <ref>.

Let us consider more carefully the case of degenerate zeros with $z_j=z_{j'}+\eps$, or $Z_j = Z_j+\eps $ in Theorem <ref>. An example would be given by the correlator $\langle\mathcal{O}_3(1)\mathcal{O}_3(N+1)\rangle$ given in (<ref>), where we set $b=a+\eps$. Then for $\eps\rightarrow 0$:

\begin{align}
\langle\mathcal{O}_3(1)\mathcal{O}_3(N+1)\rangle = (-1)^N (1-a^2)\left((1-a^2)N+1+a^2)\right)a^N (1+o(1)). \label{eq:nongenericcorrelator}
\end{align}

This limit is non-trivial in that the denominator of (<ref>) is order $\eps$—the constant term in the numerator cancels to give a finite limit. This illustrates the general case. Suppose $M$ contains $z_j$ but not $z_j+\eps$, then the term $\prod_{k \in M^c,j\in M} (\tau_k - \tau_j)^{-1}$ that appears in $C_M$ diverges as $\eps \rightarrow 0$. However, there will be a corresponding set $M'$ that is the same as $M$ except that it contains $z_j+\eps$ in place of $z_j$. Then $C_{M'}$ also diverges, but we have that  $C_M^{\vphantom N}r_M^N+C_{M'}^{\vphantom N}r_{M'}^N$ is finite. This discussion further generalises to multiple nearly degenerate zeros. This cancellation must occur in general since Day's formula is an exact formula for the determinant, and by considering the matrix elements this limit is well behaved (see the discussion in Sections <ref> and <ref>).

This analysis allows us to deduce that the formula for the order parameter given in Remark <ref> applies also for the case of degenerate zeros. In fact, this formula holds in all non-generic cases. Recall that if we have mutually inverse zeros, this model has the same ground state as the model with those zeros removed (and with a shift of $n_P$ that fixes the winding number) <cit.>. By removing the zeros, we are in a case that we have already analysed and so can evaluate the order parameter with (<ref>). However, this result agrees with keeping the zeros and evaluating (<ref>), any terms involving zeros that are mutually inverse will cancel.

§ RESULTS: CORRELATIONS IN AIII

§.§ The AIII class

We now consider a different class of models on a one-dimensional fermionic chain. Let us take $2L$ spinless fermionic degrees of freedom $c^\dagger_{A,n}$ and $c^\dagger_{B,n}$ for $n=1\dots L$. The class of AIII Hamiltonians <cit.> are charge-conserving and are symmetric under the anti-unitary involution $T$ such that $Tc^\dagger_{A,n} T=c_{A,n}$ and $T c^\dagger_{B,n}T=-c_{B,n}$ (physically $T$ is the sublattice symmetry). The class of models which are translation-invariant with respect to the two-site unit cell is given by:

\begin{align}
H = \sum_{\alpha,n} \tau_\alpha c^\dagger_{B,n}c^{\vphantom \dagger}_{A,n+\alpha} + \overline{\tau}_\alpha c^\dagger_{A,n+\alpha}c^{\vphantom\dagger}_{B,n}. \label{H_AIII}
\end{align}

We can diagonalise $H=\sum_k \lvert f(\rme^{\rmi k})\rvert\left(\eta^\dagger_{+,k}\eta^{\vphantom\dagger}_{+,k} -  \eta^\dagger_{-,k}\eta^{\vphantom\dagger}_{-,k}\right)$ with bands given[The branch of the square-root is chosen so that this is equal to $f(z)/\lvert f(z)\rvert$ for $z$ on the unit circle.] by:

\begin{align}
\eta_{\pm,k}^\dagger = \frac{1}{\sqrt{2}}\left(c^\dagger_{A,k} \pm \sqrt{\frac{f(z)}{\overline{f}(1/z)}}c^\dagger_{B,k}\right) \label{AIIImode}
\end{align}

where $f(z) =\sum_{\alpha} \tau_\alpha z^\alpha$ and $\overline{f}(z) =\sum_{\alpha} \overline{\tau}_\alpha z^\alpha$. As in BDI, the gapped phase diagram is characterised by $\omega$, the winding number of $f(\rme^{\rmi k})$. Since $f(z)$ is a Laurent polynomial, we have that $\omega = N_z-N_P$ where $N_z$ is the number of zeros inside the unit circle, and $N_P$ is the degree of the pole at zero. Examples in this class include the Su-Schrieffer-Heeger model <cit.>. Non-trivial ground state correlations may be easily derived from (<ref>) as:

\begin{align}
\overline{\langle c^\dagger_{B,n-\alpha}c^{\vphantom \dagger}_{A,n} \rangle}= \langle c^{ \dagger}_{A,n}c^{\vphantom\dagger}_{B,n-\alpha} \rangle = -\frac{1}{4\pi\rmi}\int\sqrt{\frac{f(z)}{\overline{f}(1/z)}}z^{-\alpha-1} \rmd z \label{AIIItwopoint}
\end{align}

and $\langle c^\dagger_{A,n}c^{\vphantom \dagger}_{A,n-\alpha}\rangle=\langle c^\dagger_{B,n}c^{\vphantom \dagger}_{B,n-\alpha}\rangle=\frac{1}{2}\delta_{\alpha,0}$. For the subclass of models we focus on, we give these two-point functions in Appendix <ref>.

A natural change of variables is defined by:

\begin{align}
c_{A,n}=\frac{1}{2}\left(-\gamma_{2n} + \rmi\gamma_{2n-1} \right) \qquad
c_{B,n}=\frac{1}{2}\left(\tilde\gamma_{2n-1} + \rmi\tilde\gamma_{2n} \right)
\end{align}

where the $\gamma_n,\tilde\gamma_n$ are Majorana operators as introduced in Section <ref>. The Hamiltonian (<ref>) becomes:

\begin{align}
H = \frac{\rmi}{2} \sum_{\alpha,n}\Big( \Re(\tau_\alpha) \tilde{\gamma}_n \gamma_{n+2\alpha} +  \Im(\tau_\alpha) (-\tilde{\gamma}_{2n-1} \gamma_{2n+2\alpha}+\tilde{\gamma}_{2n}\gamma_{2n+2\alpha-1})\Big)\label{H_AIII_majorana}.
\end{align}

The sublattice symmetry $T$ acts as $T \gamma_n T =\gamma_n$, $T\tilde\gamma_nT=-\tilde \gamma_n$, and we see that (<ref>) is a BDI Majorana chain, with two-site unit cell. If $\tau_\alpha \in \mathbb{R}$ then the model is actually translation-invariant, and is a special case of (<ref>), while for complex $\tau_\alpha$ this class of models is distinct.

§.§ Correlations of string operators

The string order parameters in this class will, as before, be given by string correlators of fermionic parity strings with certain end-point operators:

\begin{align}
\tilde{\mathcal{O}}_0(n) &= \exp\left(\sum_{m=1}^{n-1}\rmi\pi\left(c^\dagger_{A,m}c^{\vphantom \dagger}_{A,m}+c^\dagger_{B,m}c^{\vphantom \dagger}_{B,m}\right)\right) =(-1)^{n-1}\prod_{m=1}^{n-1} \rmi\tilde{\gamma}_{2m-1}\gamma_{2m-1}\rmi\tilde{\gamma}_{2m}\gamma_{2m}\nonumber\\
\tilde{\mathcal{O}}_\alpha(n) &=\exp\left(\sum_{m=1}^{n-1}\rmi\pi\left(c^\dagger_{A,m}c^{\vphantom \dagger}_{A,m}+c^\dagger_{B,m}c^{\vphantom \dagger}_{B,m}\right)\right)\prod_{j=n}^{n+\alpha-1}(1-2c^\dagger_{A,j}c^{\vphantom \dagger}_{A,j})\nonumber \qquad \textrm{for~}\alpha>0\\& =(-1)^{n-1}\Bigg(\prod_{m=1}^{n-1} \rmi\tilde{\gamma}_{2m-1}\gamma_{2m-1}\rmi\tilde{\gamma}_{2m}\gamma_{2m} \Bigg)\prod_{j=n}^{n+\alpha-1}(-\rmi\gamma_{2j-1}\gamma_{2j})\nonumber\\
\tilde{\mathcal{O}}_{\alpha}(n) &= \exp\left(\sum_{m=1}^{n-1}\rmi\pi\left(c^\dagger_{A,m}c^{\vphantom \dagger}_{A,m}+c^\dagger_{B,m}c^{\vphantom \dagger}_{B,m}\right)\right)\prod_{j=n}^{n+\lvert\alpha\rvert-1}(1-2c^\dagger_{B,j}c^{\vphantom \dagger}_{B,j})\nonumber\qquad \textrm{for~} \alpha<0\\&=(-1)^{n-1}\Bigg(\prod_{m=1}^{n-1} \rmi\tilde{\gamma}_{2m-1}\gamma_{2m-1}\rmi\tilde{\gamma}_{2m}\gamma_{2m} \Bigg)\prod_{j=n}^{n+\lvert\alpha\rvert-1}(-\rmi\tilde\gamma_{2j-1}\tilde\gamma_{2j}). \label{stringcorrelatorsAIII}
\end{align}

String order parameters are usually defined for interacting SPT phases <cit.>, our case is different since we have a free-fermion system—see the discussion in <cit.>. In Theorem <ref> below, we prove that these string correlators are order parameters for the different phases; i.e., $\tilde{\mathcal O}_\alpha$ has long-range order if and only if $\alpha = \omega$, the latter being the topological invariant defined above. Special cases of these fermionic string order parameters[For example, the correlator of $\left(\tilde{\mathcal{O}}_{-1}(n)-\tilde{\mathcal{O}}_{1}(n)\right)/2$ is used to detect a phase with two Majorana edge modes. This would correspond to winding number $\pm1$ in our model, consistent with $\tilde{\mathcal{O}}_{\pm1}(n)$ being the order parameter for those phases.] appear in <cit.>.  We can evaluate the ground-state correlation $\big\langle \tilde{\mathcal{O}}_\alpha(1)\tilde{\mathcal{O}}_\alpha(N+1)\big\rangle$ using Wick's theorem <cit.> and the two-point correlators given above. We show in Appendix <ref> that these correlators are equal to the square of the absolute value of a Toeplitz determinant (see also Section <ref>). This determinant can be evaluated using the same techniques as in the BDI case. We now justify that these are indeed the AIII analogues of the string order parameters in the BDI class—this follows by generalising the analysis of <cit.> to this class of models. We will not pursue this fully, but give the following result:

Consider a general model in the AIII class. This corresponds to
\begin{align}
f(z) =  \rme^{\rmi \theta} \frac{1}{z^{N_P}} \prod_{j=1}^{N_z} \left(z-z_j\right)\prod_{k=1}^{N_Z} \left(z-Z_k\right)\label{fz_genericAIII}
\end{align}

where $\lvert z_j\rvert<1$, $\lvert Z_k\rvert>1$ and $\theta\in [0,2\pi)$. The winding number is given by $\omega=N_z-N_P$ and the string correlators satisfy

\begin{align}
\lim_{N\rightarrow\infty}\big\langle\tilde{\mathcal{O}}_\alpha(1)\tilde{\mathcal{O}}_\alpha(N+1)\big\rangle = \mathrm{const} \times \delta_{\alpha,\omega}
\end{align}

The non-zero constant is the value of the order parameter, given by:

\begin{align}
\lim_{N\rightarrow\infty}\big\langle\tilde{\mathcal{O}}_\omega(1)\tilde{\mathcal{O}}_\omega(N+1)\big\rangle = \Bigg \vert\frac{ \prod_{j_1,j_2=1}^{N_z} (1-z_{j_1} \overline{z}_{j_2})
\prod_{k_1,k_2=1}^{N_Z} \left(1-Z_{k_1}^{-1} \overline{Z}_{k_2}^{-1}\right)}{\prod_{j=1}^{N_z}\prod_{k=1}^{N_Z}  \lvert 1-\frac{z_j}{ Z_k}\rvert^2}\Bigg \vert^{1/2}.
\end{align}

We outline the proof, using Szegő's theorem <cit.> and results of Hartwig and Fisher <cit.>, in Appendix <ref>.

For models with $f(z)$ of the form:

\begin{align}
f(z) = \rme^{\rmi\theta} \frac{1}{z^{n_P}} \prod_{j=1}^{n_z} \left(z-z_j\right)^2\prod_{k=1}^{n_Z} \left(z-Z_k\right)^2,  \label{fzcanonAIII}
\end{align}
we can proceed using the methods of this paper, leading to an analogue of Theorem <ref>. The generic case here is that  $\{z_{j_1}^{\vphantom {-1}},\overline{z}_{j_2}^{-1},Z_{k_1}^{\vphantom {-1}},\overline{Z}_{k_2}^{-1}\}_{j_1,j_2=1\dots n_z; k_1,k_2=1\dots n_Z}$ are pairwise distinct.

Consider a chain with a non-trivial $f(z)$ as defined in (<ref>), assuming the generic case, and for each $\alpha \in \mathbb{Z}$ define $N_{\alpha} = \max\{\lvert n_z+n_Z-n_P-\alpha \rvert,1\}$. Then:

* for $\alpha < \omega - n_z$ and $\alpha > \omega +n_Z$

\begin{align}
\big\langle\tilde{\mathcal{O}}_\alpha(1) \tilde{\mathcal{O}}_\alpha(N+1)\big\rangle  &= 0 \qquad N\geq N_{\alpha};
\end{align}

* for $\omega - n_z\leq \alpha \leq \omega +n_Z$
\begin{align}
\big\langle\tilde{\mathcal{O}}_\alpha(1) \tilde{\mathcal{O}}_\alpha(N+1)\big\rangle  &= \Big\lvert \sum_{M} C_M^{\vphantom N} r_M^N \Big\rvert^2 \qquad N\geq N_{\alpha},\label{eq:thm1b}
\end{align}

where for $n_z>0$ and $n_Z>0$:
* for $\omega - n_z\leq \alpha\leq \omega +n_Z-n_z$ the sum is over sets $M$ that label sets of zeros, with the non-zero constants $C_M$ and $r_M$ defined in Case 1 below.

* for $\omega +n_Z-n_z \leq \alpha\leq \omega +n_Z$ the sum is over sets $M$ that label sets of inverse zeros, with the non-zero constants $C_M$ and $r_M$ defined in Case 2 below.
If either $n_z=0$ or $n_Z=0$:

* if $n_z =0$, and $\alpha<\omega+n_Z$, then put $n_z=1$, $\alpha \rightarrow \alpha+2$ and take the limit $z_1 \rightarrow 0$ with Case 1. If $n_z=0$ and $\alpha = \omega+n_Z$ then the correlator is given by $ \prod_{k=1}^{n_Z} \lvert Z_k\rvert^{-2N}$.

* if $n_Z =0$ and $\alpha>\omega-n_z$, then put $n_Z=1$, and take the limit $Z_1 \rightarrow \infty$ with Case 2. If $n_Z =0$ and $\alpha=\omega-n_z$, then the correlator is given by: $\prod_{j=1}^{n_z} \lvert z_j\rvert^{2N}$.

In both cases, in the formula (<ref>), the sum is over all subsets $M \subseteq \{1,\dots,n_z+n_Z\}$ of some fixed size $\lvert M \rvert$. For each such $M$, we also define:

\begin{align}
M^c &=  \{1,\dots,n_z+n_Z\}\setminus M  \end{align}
and further define $M_z = \{1,\dots, n_z\}$ and $M_Z= \{1,\dots, n_Z\}$. Then we have:

Case 1:
Define $\tau_j = z_j$ for $1\leq j \leq n_z$ and $\tau_{n_z+j}=Z_j$ for  $1\leq j \leq n_Z$.
Then the sum in (<ref>) is over all subsets $M$ of size $\lvert M \rvert =\alpha+n_P-n_z$, and $C_M$ and $r_M$ are defined by:

\begin{align}
r_M&=  \frac{\prod_{k \in M^c} \tau_k} {\prod_{j=1}^{n_Z}\overline{Z}_j},\nonumber\\\nonumber
C_M &= \prod_{k \in M^c,m\in M_Z} (\tau_k - \overline{Z}^{-1}_m) \prod_{l \in M_z, j\in M} (\overline{z}_l^{-1} - \tau_j) \prod_{l\in M_z}z_l^{-(n_z+n_Z -n_P-\alpha)} \\&\times\prod_{l \in M_z, m\in M_Z} (\overline{z}_l^{-1} - \overline{Z}^{-1}_m)^{-1} \prod_{k \in M^c,j\in M} (\tau_k - \tau_j)^{-1}  \prod_{k\in M^c}\tau_k^{-(n_z+n_Z -n_P-\alpha)}.
\end{align}

Case 2:
Define $\tau_j = z_j^{-1}$ for $1\leq j \leq n_z$ and $\tau_{n_z+j}=Z_j^{-1}$ for  $1\leq j \leq n_Z$.
Then the sum in (<ref>) is over all subsets $M$ of size $\lvert M \rvert =2n_z+n_Z-n_P-\alpha$, and $C_M$ and $r_M$ are defined by:
\begin{align}
r_M&=  \prod_{j=1}^{n_z}z_j\prod_{k=1}^{n_Z}\frac{Z_k}{\overline{Z}_k}{\prod_{k \in M^c} \tau_k} ,\nonumber\\
C_M &= \prod_{k \in M^c,m\in M_z} (\tau_k - \overline{z}_m) \prod_{l \in M_Z, j\in M} (\overline{Z}_l - \tau_j)\prod_{l \in M_Z} \overline{Z}_l^{n_P+\alpha-n_z-n_Z}\nonumber\\&\times \prod_{l \in M_Z, m\in M_z} (\overline{Z}_l - \overline{z}_m)^{-1} \prod_{k \in M^c,j\in M} (\tau_k - \tau_j)^{-1} \prod_{k \in M^c}\tau_k^{-(n_P+\alpha-n_z-n_Z)}.
\end{align}

Much of the discussion in Section <ref> applies analogously here. We will simply give one further example.

Let $f(z) = (z-a)^2(z-b)^2$ where $\lvert a\rvert<1$ and $\lvert b\rvert>1$. Since $n_z =n_Z= 1$, we are in the phase $\omega =2$, and the correlator $\langle\tilde{\mathcal{O}}_2(1) \tilde{\mathcal{O}}_2(N+1)\rangle$ is given by:

\begin{align}
\big\langle\tilde{\mathcal{O}}_2(1) \tilde{\mathcal{O}}_2(N+1)\big\rangle &=\Bigg\vert \frac{\left(1-\lvert a\rvert^2\right) (\lvert b\rvert^{2}-1) }{\lvert b-a\rvert^2}\left(\frac{b}{\overline{b}}\right)^N + \frac{\lvert 1-b\overline{a}\rvert^2}{\lvert b-a\rvert ^2}\left(\frac{a}{\overline{b}}\right)^N \Bigg\vert^2 \qquad\qquad N\geq 1\nonumber\\
&=\underbrace{\left(\frac{\left(1-\lvert a\rvert^2\right) (\lvert b\rvert^{2}-1) }{\lvert b-a\rvert^2}\right)^2}_{\textrm{Order parameter}}+\frac{\left(1-\lvert a\rvert^2\right) (\lvert b\rvert^{2}-1)\lvert 1-b\overline{a}\rvert^2 }{\lvert b-a\rvert^4}\left(\left(\frac{a}{b}\right)^N +\left(\frac{\overline{a}}{\overline{b}}\right)^N \right)\nonumber\\
& \hspace{5cm}  +\frac{\lvert 1-b\overline{a}\rvert^4}{\lvert b-a\rvert ^4}\left(\frac{\lvert a\rvert^2}{\lvert b \rvert^2}\right)^N \qquad\qquad N\geq 1.
\end{align}

Note that for $a,b\in\mathbb{R}$ this example reduces to a translation-invariant BDI model. Given the form of (<ref>), we see that we have two decoupled BDI chains, each with $f(z) =(z-a)^2(z-b)^2 $. The correlator $\langle{\mathcal{O}}_2(1) {\mathcal{O}}_2(N+1)\rangle$ can be calculated for each of these decoupled chains using Example <ref> and we see that in that case $\langle{\mathcal{O}}_2(1) {\mathcal{O}}_2(N+1)\rangle^2=\big\langle\tilde{\mathcal{O}}_2(1) \tilde{\mathcal{O}}_2(N+1)\big\rangle $.

§.§ Correlation matrix

For the AIII class, the correlation matrix for a subsystem containing $N$ two-site unit cells is defined by:

\begin{align}
C_N= \left(\langle b^\dagger_j b^{\vphantom{\dagger}}_k \rangle\right)_{j,k=1\dots 2N},
\end{align}

where we define $b_{2j-1}=c_{A,j}$ and $b_{2j}=c_{B,j}$. Following the discussion in <cit.>, to find the eigenvalues of the reduced density matrix or, equivalently, the entanglement spectrum, we want the eigenvalues of the matrix $\hat{A}_N= 1-2C_N$. The eigenvalues, $\lambda$, of $\hat{A}$ are real and satisfy $-1\leq \lambda \leq 1$. The values $\lambda =\pm 1$ are trivial eigenvalues that do not correspond to entanglement. Our main result for class AIII is that the proof of Theorem <ref> carries over, leading to:

Consider a Hamiltonian of the form (<ref>) with a non-trivial $f(z)$ of the form (<ref>) such that $n_P=n_z+n_Z$. Recall that in this class zeros need not come in conjugate pairs, and we allow the non-generic case where zeros may coincide, but assume without loss of generality[As in Remark <ref>, if there were such zeros they can be removed without affecting the ground state.] that $z_j\neq \overline{Z}_k^{-1}$ for any $j,k$. Denote the correlation matrix on a subsystem of size $N$ by $C_N$ and define $\hat{A}_N=1-2C_N$. Then:

\begin{align}
\det(\lambda-\hat{A}_N)=   (\lambda^2-1)^N \left(\frac{\prod_{j=1}^{n_z} \overline{z}_j^2\prod_{j=1}^{n_Z} \overline{Z}_j^2}{\left(\prod_{j=1}^{n_Z} \lvert Z_j\rvert^2\right)^2}\right)^N \frac{\det{M(N,\lambda)}}{\det{M(0,\lambda)}}, \end{align}

where $M(n,\lambda)$ is a $4(n_z+n_Z)\times 4(n_z+n_Z)$ matrix that can be determined from the zeros of $f(z)$ and the Smith canonical form of a known $2\times 2$ matrix polynomial (defined in (<ref>)-(<ref>)). Moreover, there exist $d', d'' \in \mathbb{N}$ (independent of $N$) such that:

\begin{align}
\det{M(n,\lambda)} = \mu \left(\frac{\prod_{j=1}^{n_z} \overline{z}_j^2\prod_{k=1}^{n_Z} \overline{Z}_k^2}{\left(\prod_{k'=1}^{n_Z} \lvert Z_{k'}\rvert^2\right)^2}\right)^{-n} \prod_{j_1=1}^{d'}(\lambda^2-\tilde x_{j_1}(n))\prod_{j_2=1}^{d''}(\lambda-c_{j_2}) , \end{align}

where the constants $\mu$ and $c_j$ are independent of $n$ and  $\tilde x_j(n)$ satisfies $\tilde x_j(0)=1$, $\tilde x_j(n)\leq 1$ for $n\geq 1$ and $\tilde x_j(n)<1$ for some $n>0$.

For any $f(z)$ as defined in (<ref>) such that $n_P \neq n_z+n_Z$, then the characteristic polynomial of the correlation matrix can be found by taking a limit of a related case where $n_P = n_z+n_Z$ that can be evaluated using the above. Details of this limit are discussed in Remark <ref>—in particular, this discussion does not rely on zeros of $f(z)$ appearing in complex-conjugate pairs.

For further details see Section <ref> and Appendix <ref>. As in the BDI case, this results allows us to conclude that there exists an exact MPS representation of the ground state.

To find the entanglement spectrum requires using the definition of $M(n,\lambda)$ given in Section <ref> (with the changes for class AIII given in Appendix <ref>), as well as finding the Smith canonical form of the relevant matrix. By analogous calculations[The matrix from the Smith canonical form needed for this calculation is given in Appendix <ref>.] to those in Example <ref>,
we find that for $f(z) = \frac{1}{z}(z-b)^2$ with $\lvert b\rvert <1, b\in \mathbb{C}$:

\begin{align}
\det(\lambda-\hat{A}_N)&= ( \lambda^2-1)^{N-1}(\lambda^2-\lvert b\rvert^{2N}). \end{align}

§ THE TRANSFER MATRIX

As discussed in Remark <ref>, Theorems <ref> and <ref> allow us to deduce the existence of an exact MPS representation for the ground state of BDI and AIII models with $f(z)$ of the form given in (<ref>) and (<ref>) respectively. Moreover, an exact MPS representation is constructed for the BDI class in <cit.>. Given an MPS, it can be put into a canonical form <cit.>—we will henceforth suppose that the MPS under consideration is in this form. From the MPS, we can define a transfer matrix and in this section we show how Theorem <ref> allows us to deduce properties of this transfer matrix for the BDI case. As an application of this, we can find a lower bound on the bond dimension of an MPS representation of the ground state—this complements the upper bound found in <cit.>. In certain cases, the upper and lower bounds coincide and allow us to give the optimal bond dimension of such an MPS. Throughout this section we will work with the spin chain model (see Appendix <ref>), and use standard graphical notation to give an intuitive illustration of some formulae <cit.>. While we focus on the BDI case, an analogous discussion based on Theorem <ref> would allow us to deduce properties of the transfer matrix in the AIII class.

at (-4.5,0) a);

(nodeA) at (-3,0.5) [thick,minimum size=.15cm,circle,draw] $A$;
at (-3,-.5) $j$;
at (-4,0.5) $\alpha$;
at (-2,0.5) $\beta$;
[rounded corners] (node cs:name=nodeA, angle = 180)–+(-.3,0.);
[rounded corners] (node cs:name=nodeA, angle = 0) –+(.3,0.);
[rounded corners] (node cs:name=nodeA, angle = -90) –+(0,-0.3);

at (1,0) b);

(nodeE) at (2.5,0) [thick,minimum size=.15cm,circle,draw] $E_\mathbb{X}$;
[rounded corners] (node cs:name=nodeE, angle = 70) – +(0,.5)–+(.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -70) – +(0,-.5)–+(.75,-0.5);
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.5)–+(-.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-.75,-0.5);

at (4,0) $=$;

(nodeA1) at (5.5,1) [thick,minimum size=.15cm,circle,draw] $A$;
[rounded corners] (node cs:name=nodeA1, angle = 180)–+(-.5,0.);
[rounded corners] (node cs:name=nodeA1, angle = 0) –+(.5,0.);
(nodeX) at (5.5,0)  [thick,minimum size=.15cm,circle,draw] $\mathbb{X}$;
(nodeA2) at (5.5,-1) [thick,minimum size=.15cm,circle,draw] $\overline{A}$;
[rounded corners] (node cs:name=nodeA2, angle = 180)–+(-.5,0.);
[rounded corners] (node cs:name=nodeA2, angle = 0) –+(.5,0.);
[rounded corners] (node cs:name=nodeA1, angle = -90) –+(0,-0.25);
[rounded corners] (node cs:name=nodeA2, angle = 90) –+(0,0.25);

at (1,0) c);
[scale=2] at (1.5,0) $\Bigg($;

(nodeE) at (3,0) [thick,minimum size=.15cm,circle,draw] $E_\mathbb{I}$;
[rounded corners] (node cs:name=nodeE, angle = 70) – +(0,.5)–+(.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -70) – +(0,-.5)–+(.75,-0.5);
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.5)–+(-.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-.75,-0.5);
[scale=2] at (4.5,0) $\Bigg)^n$;

at (5,0) $=$;
[scale=1] at (6.5,0.) $\sum_{k=1}^{\chi^2} \mu_k^n$;

(noder) at (8,0) [thick,minimum size=.15cm,circle,draw] $r_k$;
[rounded corners] (node cs:name=noder, angle = 90) – +(0,.5)–+(-.75,0.5);
[rounded corners] (node cs:name=noder, angle = -90) – +(0,-.5)–+(-.75,-0.5);
(nodel) at (9,0) [thick,minimum size=.15cm,circle,draw] $l_k$;
[rounded corners] (node cs:name=nodel, angle = 90) – +(0,.5)–+(.75,0.5);
[rounded corners] (node cs:name=nodel, angle = -90) – +(0,-.5)–+(.75,-0.5);

Graphical notation: a) is the MPS tensor with bond indices $\alpha, \beta$; b) is the generalised transfer matrix; c) is a formula for $E_\mathbb{I}^n$ when the transfer matrix is diagonalisable.

§.§ Overview

For a periodic spin-12 chain with $L$ sites, a translation-invariant MPS representation of a state is of the form

\begin{align}
\ket{\psi} = \sum_{j_1, \dots j_n = \pm 1/2} \tr(A_{j_1}\dots A_{j_L})\ket{j_1,\dots,j_L} .\label{MPSdef}
\end{align}

We will take $\ket{\psi}$ to be the ground state of a model defined in (<ref>).
For a fixed value of $j$, $A_{j}$ is a $\chi \times \chi$ matrix, and one can think of $A_{j}$ as a $\chi\times\chi\times2$ tensor. Given $A_{j}$, one can compute various quantities for a system of size $L$ and then take the limit $L\rightarrow \infty$. A useful operator is the (generalised) transfer matrix:

\begin{align}
E_\mathbb{X} = \sum_{j,k=\pm 1/2} \mathbb{X}_{kj} A_{j}  \otimes \overline{A}_{k} \quad\qquad \textrm{for} \quad \mathbb{X}_{kj} = \bra{k}\mathbb{X}\ket{j}.
\end{align}

We can view this as a $\chi^2 \times \chi^2$ matrix, where we group indices connecting sites to the left and sites to the right[If we put the indices in we have: $E^{(\alpha\alpha';\beta\beta')}_\mathbb{X} = \sum_{j,k=\pm 1/2} \mathbb{X}_{kj} A^{(\alpha;\beta)}_{j}  \overline{A}^{(\alpha';\beta')}_{k}$.]
(see Figure <ref>).  The transfer matrix, $E_\mathbb{I}$, can be used to calculate the norm of the state, through

\begin{align}
\braket{\psi}{\psi}&=\lim_{L\rightarrow \infty} \tr\left(E_\mathbb{I}^L\right). \label{MPSnorm}\end{align}

Note that $A_j$ can be rescaled so that the state is normalised in the limit, the canonical form fixes a rescaling <cit.>. Moreover, correlations of single-site operators $\mathbb{X}_1$ and $\mathbb{X}_{N+1}$ take the form:

\begin{align}
\bra{\psi}\mathbb{X}_1\mathbb{X}_{N+1}\ket{\psi}&=\lim_{L\rightarrow \infty} \tr\left(E_\mathbb{X}^{\vphantom{N-1}}E_\mathbb{I}^{N-1}E_\mathbb{X}^{\vphantom{N-1}}E_\mathbb{I}^{L-(N+1)}\right).\label{MPScorr}\end{align}

It is simple to extend this to more general correlators, for example
\begin{align}
\bra{\psi}\mathbb{X}_1\mathbb{X}'_2\mathbb{X}_{N+1}\ket{\psi}&=\lim_{L\rightarrow \infty} \tr\left(E_\mathbb{X}^{\vphantom{N-1}}E_\mathbb{X'}^{\vphantom{N-1}}E_\mathbb{I}^{N-2}E_\mathbb{X}^{\vphantom{N-1}}E_\mathbb{I}^{L-(N+1)}\right).\end{align}

Thus, the transfer matrix $E_\mathbb{I}= \sum_{j,k=\pm 1/2}  A_{j}  \otimes \overline{A}_{j}$ is intimately related to correlations—indeed in this limit our MPS is a finitely correlated state <cit.>. We now discuss how the results of Theorem <ref> relate to this framework.

§.§ General discussion

Viewing the transfer matrix $E_\mathbb{I}$ as a $\chi^2 \times \chi^2$ matrix, then it has a Jordan normal form with $\chi^2$ eigenvalues (counting multiplicity).
Let us first consider the case where $E_\mathbb{I}$ is diagonalisable. For $A$ in canonical form, the largest eigenvalue of $E_\mathbb{I}$ (in absolute value) equals one, and moreover, for the purposes of this discussion, let us suppose this largest eigenvalue is unique (see Remark <ref>).
Then (<ref>) simplifies to:

\begin{align}
\bra{\psi}\mathbb{X}_1\mathbb{X}_{N+1}\ket{\psi}= \sum_{k=1}^{\chi^2} \mu_k^{N-1} \bra{l_1} E_\mathbb{X} \ket{r_k}\bra{l_k}E_\mathbb{X} \ket{r_1},\label{MPScorr2}
\end{align}

where $\mu_1=1$, $\lvert \mu_k \rvert < 1$ for $k>1$ and $\{\ket{r_k},\bra{l_j}\}$ are the left and right eigenvectors of $E_\mathbb{I}$, satisfying $\braket{l_j}{r_k}=\delta_{jk}$.

Now, consider the operators $\mathcal{O}_\alpha$ defined in (<ref>). For $\alpha=2k+1$ these are correlators of local operators, for example:

\begin{align}
\langle \mathcal{O}_{1}(1) \mathcal{O}_{1}(N+1) \rangle&= \sum_{k=1}^{\chi^2} \mu_k^{N-1} \bra{l_1} E_{X} \ket{r_k}\bra{l_k}E_{X} \ket{r_1}\\
\langle \mathcal{O}_{-1}(1) \mathcal{O}_{-1}(N+1) \rangle&= \sum_{k=1}^{\chi^2} \mu_k^{N-1} \bra{l_1} E_{Y} \ket{r_k}\bra{l_k}E_{Y} \ket{r_1}\\
\langle \mathcal{O}_{3}(1) \mathcal{O}_{3}(N+1) \rangle&= \sum_{k=1}^{\chi^2} \mu_k^{N-3} \bra{l_1} E_{X}E_{Y}E_{X} \ket{r_k}\bra{l_k}E_{X}E_{Y}E_{X} \ket{r_1}.
\end{align}

Thus, using the results of Theorem <ref>, the terms proportional to $r_M^N$ that appear mean we can identify eigenvalues of the transfer matrix: $\{(-1)^{n_P+1}r_M\} \subseteq \{\mu_k \}$. Moreover, the terms $C_M$ will correspond to the overlaps such as $\bra{l_1} E_{X} \ket{r_k}$.
Note that if $E_\mathbb{I}$ had off-diagonal terms in the Jordan normal form, then these would appear in correlators as, e.g., $n \mu_k^{n-1}$. Such correlations do not appear in Theorem <ref> in the generic case, although they can appear for non-generic cases, for example in (<ref>). In any case, since non-trivial Jordan blocks correspond to degenerate eigenvalues, the dimension of the transfer matrix is always lower bounded by the number of distinct $r_M$.

Now, for $\mathcal{O}_\alpha$ with $\alpha =2k$ we have a non-local operator. However, we can do a similar analysis, using the generalised transfer matrix $E_Z$. For example, we have:

\begin{align}
\langle \mathcal{O}_{2}(1) \mathcal{O}_{2}(N+1) \rangle&=\lim_{L\rightarrow \infty} \tr\left(E_X^{\vphantom{N-1}}E_Y^{\vphantom{N-1}}E_{Z}^{N-2}E_Y^{\vphantom{N-1}}E_X^{\vphantom{N-1}}E_\mathbb{I}^{L-(N+2)}\right)\\&= \sum_{k=1}^{\chi^2} \tilde\mu_k^{N-2} \bra{l_1\vphantom{\tilde l}} E_{X}E_{Y} \ket{\tilde{r\vphantom{l}}_k}\bra{\tilde{l}_k}E_{Y}E_{X} \ket{\vphantom{\tilde{l}}r_1}
\end{align}

where we suppose for simplicity that we can diagonalise
\begin{align}E_Z= \sum_{k=1}^{\chi^2} \tilde{\mu}_k \ket{\tilde{r\vphantom{l}}_k}\bra{\tilde{l}_k}.\end{align}

Because $\prod_j Z_j$ is a symmetry of our system, it can be shown that $E_Z$ and $E_\mathbb{I}$ have the same eigenvalues up to phase factors <cit.>. This follows from the result that a symmetry operator on the physical index corresponds to a transformation on the bond indices of the form $A\rightarrow WUAU^\dagger$ where $U$ is unitary and $W$ is a diagonal matrix of phase factors $\rme^{\rmi\theta_j}$. This is illustrated in Figure <ref>.
As before, we can use Theorem <ref> to identify eigenvalues $\tilde \mu_k$ of $E_Z$. Moreover, we know that for each $k$, $\lvert \tilde \mu_k \rvert = \lvert \mu_j \rvert$ for some $j$. Hence, we can lower bound $\chi^2$ by the number of $r_M$ that give us distinct $\mu_j$ and distinct $\tilde{\mu}_j$ where none of these $r_M$ have the same absolute value.

An MPS is called injective if the transfer matrix has a unique largest eigenvalue[In some literature this is called a pure MPS <cit.>, and injective has a slightly different definition <cit.>.]. In the models considered in this paper, the MPS ground state is not necessarily injective. Indeed, suppose that $\omega=2k+1$ for $k\in \mathbb{Z}$. Then we have a local order parameter $\mathcal{O}_{2k+1}$. Note that $\mathcal{O}_{2k+1}$ is odd under the symmetry $\prod_j Z_j$. Hence, given a unique largest eigenvalue, we must have:

\begin{align}
\langle \mathcal{O}_{\omega}(1)\rangle = \bra{l_1} E_{X}E_{Y}\dots E_{X} \ket{r_1}= 0 \label{eq:mixed}
\end{align}

where for definiteness we fix $\omega>0$. From (<ref>), the two-point correlation function will behave like:

\begin{align}
\lim_{N\rightarrow\infty}\langle \mathcal{O}_{\omega}(1) \mathcal{O}_{\omega}(N+1) \rangle=\bra{l_1} E_{X}E_{Y}\dots E_{X} \ket{r_1}\bra{l_1}E_{X}E_{Y}\dots E_{X} \ket{r_1},\label{eq:mixed2}
\end{align}

where the left-hand-side has a non-zero limit given in Remark <ref>. Equations (<ref>) and (<ref>) are inconsistent, and so for $\omega =2k+1$ we must have further eigenvalues of absolute value one. This means that if the MPS is injective, we must have $\omega =2k$.

It is simple to generalise (<ref>) to the non-injective case, and the relation between $E_Z$ and $E_\mathbb{I}$ illustrated in Figure <ref> does not rely on the MPS being injective. However, if we do have an injective MPS, this means that $W$ is a matrix of the form $\rme^{\rmi \theta} \mathbb{I}$. By applying the symmetry transformation on the physical index twice we have $A \rightarrow \rme^{2\rmi\theta} U AU^\dagger$. Since the transfer matrix $E_{Z^2}=E_\mathbb{I}$, we can conclude that $\rme^{\rmi \theta}=\pm1$. Hence, in the injective case, either the spectrum of $E_Z$ is the same as the spectrum of $E_\mathbb{I}$ or it is the same as the spectrum of $-E_\mathbb{I}$.

Throughout the paper we have considered the generic case to be where the zeros and inverse zeros $\{z_{j_1}^{\vphantom {-1}},z_{j_2}^{-1},Z_{k_1}^{\vphantom {-1}},Z_{k_2}^{-1}\}_{j_1,j_2=1\dots n_z; k_1,k_2=1\dots n_Z}$ are pairwise distinct.
For the purposes of finding a lower bound on the bond dimension, we will now suppose the following `strongly generic' condition.
For $m_z\subseteq\{1,\dots,n_z\}$ and $m_Z\subseteq\{1,\dots,n_Z\}$
define the products:

\begin{align}
R(m_z,m_Z)=\prod_{j \in m_z} z_j \prod_{k \in m_Z} Z_k^{-1}.
\end{align}

The `strongly generic' condition is that $\lvert R(m_z,m_Z) \rvert =\lvert R(m'_z,m'_Z) \rvert$ if and only if $m_z$ and $m_{z}'$ (similarly, $m_Z$ and $m_Z'$) differ only by replacing the index of any zero by the index of its complex-conjugate[Recall that since $f(z)$ has real coefficients, zeros are either real or come in complex-conjugate pairs]. This is a natural condition given the form of $f(z)$, and moreover implies the usual generic case assumed above.

(nodeE) at (2.5,0) [thick,minimum size=1cm,circle,draw] $E_{Z}$;
[rounded corners] (node cs:name=nodeE, angle = 70) – +(0,.5)–+(.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -70) – +(0,-.5)–+(.75,-0.5);
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.5)–+(-.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-.75,-0.5);

at (4,0) $=$;

(nodeA1) at (5.5,1) [thick,minimum size=.15cm,circle,draw] $A$;
[rounded corners] (node cs:name=nodeA1, angle = 180)–+(-.5,0.);
[rounded corners] (node cs:name=nodeA1, angle = 0) –+(.5,0.);
(nodeX) at (5.5,0)  [thick,minimum size=.15cm,circle,draw] $Z$;
(nodeA2) at (5.5,-1) [thick,minimum size=.15cm,circle,draw] $\overline{A}$;
[rounded corners] (node cs:name=nodeA2, angle = 180)–+(-.5,0.);
[rounded corners] (node cs:name=nodeA2, angle = 0) –+(.5,0.);
[rounded corners] (node cs:name=nodeA1, angle = -90) –+(0,-0.25);
[rounded corners] (node cs:name=nodeA2, angle = 90) –+(0,0.2);
at (7,0) $=$;
(nodeU1) at (6.5,1) [thick,minimum size=.15cm,circle,draw] $U^\dagger$;
(nodeU2) at (4.5,1) [thick,minimum size=.15cm,circle,draw] $U^{\vphantom{\dagger}}$;
(nodeU3) at (3.5,1) [thick,minimum size=.15cm,circle,draw] $W$;
(nodeA1) at (5.5,1) [thick,minimum size=.15cm,circle,draw] $A$;
[rounded corners] (node cs:name=nodeA1, angle = 180)–+(-.22,0.);
[rounded corners] (node cs:name=nodeA1, angle = 0) –+(.17,0.);
[rounded corners] (node cs:name=nodeU3, angle = 0) –+(.15,0.);
[rounded corners] (node cs:name=nodeU3, angle = 180) –+(-.25,0.);
[rounded corners] (node cs:name=nodeU1, angle = 0) –+(.25,0.);
(nodeA2) at (5.5,-1) [thick,minimum size=.15cm,circle,draw] $\overline{A}$;
[rounded corners] (node cs:name=nodeA2, angle = 180)–+(-2.3,0.);
[rounded corners] (node cs:name=nodeA2, angle = 0) –+(1.3,0.);
[rounded corners] (node cs:name=nodeA1, angle = -90) – (node cs:name=nodeA2, angle = 90);

at (8,0) $=$;
(nodeU1) at (3.5,1) [thick,minimum size=.15cm,circle,draw] $U^\dagger$;
(nodeU2) at (1.5,1) [thick,minimum size=.15cm,circle,draw] $U^{\vphantom{\dagger}}$;
(nodeU3) at (0.5,1) [thick,minimum size=.15cm,circle,draw] $W$;
[rounded corners] (node cs:name=nodeU3, angle = 0) –+(.15,0.);
[rounded corners] (node cs:name=nodeU1, angle = 0) –+(.25,0.);
[rounded corners] (node cs:name=nodeU3, angle = 180) –+(-.25,0.);
(nodeE) at (2.5,0) [thick,minimum size=1cm,circle,draw] $E_\mathbb{I}$;
[rounded corners] (node cs:name=nodeE, angle = 70) – +(0,.52)–(node cs:name=nodeU1, angle = 180);
[rounded corners] (node cs:name=nodeE, angle = -70) – +(0,-.5)–+(1.55,-0.5);
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.52)–(node cs:name=nodeU2, angle = 0);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-2.5,-0.5);

Graphical representation of the relationship between $E_Z$ and $E_\mathbb{I}$. $U$ is a unitary matrix, and $W$ is a diagonal matrix of phase factors. This follows from Theorem 5 of <cit.>.

§.§ Lower bound on the bond dimension

We will now consider the number of different $r_M$ that appear in Theorem <ref>, and use this to find a lower bound on the bond dimension. For ease of notation, we analyse the case $n_P=n_z+n_Z$, where, for $n_P$ even, we can find a lower bound that coincides with the upper bound in <cit.>. The general lower bound then follows easily. We remark that the $N_{\alpha}$ that appears in Theorem <ref> has a natural explanation for this value of $n_P$. Recall that for $N\geq N_{\alpha}$ the given formula for the correlator applies, where $N_{\alpha} =\max\{\lvert \alpha \rvert,1\}$ (for $N<N_{\alpha}$ we can take a limit to evaluate the correlator). This is exactly the number of sites that the local endpoint operator of $\mathcal{O}_\alpha$ occupies—in particular, this is the value of $N$ for which the product of $E_\mathbb{X}$ at each endpoint is fixed, and as $N$ increases we simply increase the number of $E_\mathbb{I}$ or $E_Z$ that appear in between the endpoints.

Suppose that $f(z)$ satisfies the strongly generic condition given in Remark <ref>. Then

\begin{align}
2^{n_z+n_Z}\leq \chi^2.
\end{align}

In <cit.>, the upper bound  $ \chi^2 \leq 2^{2 \lceil\textrm{range}(H)/2\rceil}$ is derived, where $\textrm{range}(H)$ is defined in Remark <ref>. If $n_P=n_z+n_Z$, then $\textrm{range}(H)= n_z+n_Z$. Hence, when $n_P=n_z+n_Z$ is even (implying also that $\omega$ is even),
combining the two inequalities gives us that $2^{n_z+n_Z}= \chi^2$. When $n_P=n_z+n_Z$ is odd, and for other values of $n_P$, $2 \lceil\textrm{range}(H)/2\rceil \geq n_z+n_Z$. Hence, these two inequalities do not fix the optimal value of $\chi$ in these cases, although we conjecture that the upper bound is optimal;[More precisely, we make this conjecture given the generic condition used throughout the paper: i.e., that  $\{z_{j_1}^{\vphantom {-1}},z_{j_2}^{-1},Z_{k_1}^{\vphantom {-1}},Z_{k_2}^{-1}\}_{j_1,j_2=1\dots n_z; k_1,k_2=1\dots n_Z}$ are pairwise distinct. If this does not hold then there are counter examples. For example, $f(z)=z^{-2}(z-2)^2(z-1/2)^2$ has a product state ground state with bond dimension $\chi=1$, while $\textrm{range}(H)=2$.] i.e., $\log_2(\chi)=\lceil\textrm{range}(H)/2\rceil $.

§.§.§ Proof of Corollary <ref>

Let us fix $n_P=n_z+n_Z$, Theorem <ref> gives non-zero correlators for $-n_Z\leq \alpha\leq n_z$. Note that there are $2^{n_z+n_Z}$ subsets of $S=\{z_1,\dots z_{n_z}, Z_1^{-1},\dots Z_{n_Z}^{-1}\}$. Given any subset, one can take the product of the elements contained within it. Clearly all $r_M$ appearing in Theorem <ref> are of this form. Moreover, we show in Appendix <ref> that for any such product, there is a set $M$ that appears in $\langle \mathcal{O}_{\alpha}(1) \mathcal{O}_{\alpha}(N+1) \rangle =\sum_M C_M^{\vphantom N}r_M^N$ for $-n_Z\leq \alpha\leq n_z$ such that $r_M$ gives this product. The constants $C_M$ are non-zero, so we can identify $r_M$ as eigenvalues $\mu_k$ of $E_\mathbb{I}$ or $\tilde{\mu}_k$ of $E_Z$. The local (non-local) correlation functions are those with $\alpha$ odd (even). Note that in Theorem <ref>, for a fixed $\alpha$ the products $r_M$ that appear always correspond to either an even-size or an odd-size subset of $S$. If $\omega$ is odd, then
the local order parameter contains only even-size subsets, since there must be an $M$ such that $r_M=1$. If $\omega$ is even, then the non-local string order parameter contains only even-size subsets for the same reason. Moreover, if we shift $\alpha$, we alternate between correlators containing only even- and only odd-size subsets. Suppose $\omega$ is odd, then products of terms in all of the even-size subsets of $S$ (and only these products) appear as $r_M$ in correlations of local operators. Using the strongly generic condition, this gives us $2^{n_z+n_Z-1}$ different $\mu_k$. Moreover, products of terms in all of the odd-size subsets of $S$ (and only these products) appear as $r_M$ in correlations of non-local operators. This gives us $2^{n_z+n_Z-1}$ different $\tilde\mu_k$. Since each $\tilde\mu_k$ is a product of an odd-size subset in $S$, by the strongly generic condition, it cannot have the same absolute value as any $\mu_j$ that appears in a local correlator. This means that there must be additional eigenvalues of $E_\mathbb{I}$ that correspond to each of these $\tilde\mu_k$. If $\omega$ is even the same argument goes through exchanging odd-size and even-size. Hence, we can conclude that there are at least $2^{n_z+n_Z}$ different eigenvalues of $E_\mathbb{I}$, i.e.,  $2^{n_z+n_Z}\leq \chi^2$. Moreover, for $n_P= n_z+n_Z+k$,  we simply shift the labels of the $\mathcal{O}_\alpha$ by $\alpha \rightarrow \alpha-k$ in each formula. This means that  $2^{n_z+n_Z}\leq \chi^2$ holds in general.

§.§ Spectrum of the transfer matrix

In the proof of Corollary <ref>, we identify eigenvalues of the transfer matrix, up to complex phase factors, in order to find a lower bound on the bond dimension. In the case that we have $n_P=n_z+n_Z$ is an even integer, and using the upper bound on the bond dimension from <cit.>, we can go further and find the spectrum of the transfer matrix. This is summarised as follows:

Consider a model of the form given in (<ref>), with $\sigma = \pm 1$, that satisfies the strongly generic condition given in Remark <ref> and has $n_P=n_z+n_Z=2n$ for some $n\in\mathbb{Z}$.
We define $S=\{z_1,\dots z_{n_z}, Z_1^{-1},\dots Z_{n_Z}^{-1}\}$.
Then the MPS is injective and the $2^{n_z+n_Z}$ eigenvalues of the transfer matrix are labelled by subsets $s\subseteq S$ and are given by:

\begin{align}
\mu(s) = (-\sigma)^{\lvert s \rvert}\prod_{\tau_j \in s}   \tau_j. \label{eq:spectrum}
\end{align}

In Remark <ref> we show that subject to the strongly generic condition, and with $n_P=n_z+n_Z+k=2n$, then the non-zero eigenvalues are given by (<ref>), and all further eigenvalues of the transfer matrix are zero. Note that if $n_P$ is even, then so is $\omega$. We do not give results for the spectrum of the transfer matrix in cases where $\omega$ is odd.

Let us define $\rme^{\rmi \theta}=-\sigma$ and
\begin{align}
\eps_j &=\begin{cases} -\log{z_j}  \qquad &1\leq j\leq n_z \\
-\log{Z_{j-n_z}^{-1}}&n_z+1\leq j \leq n_z+n_Z
\end{cases},
\end{align}

where $\eps_j \in \mathbb{C}$ in general.
Then given the assumptions in Corollary <ref>, we have that the eigenvalues of the transfer matrix are given by:

\begin{align}
\mu_{n_1,\dots n_{n_z+n_Z}} =\rme^{-\sum_j n_j (\eps_j+\rmi \theta)}.
\end{align}

We can view these eigenvalues as coming from an, in general non-hermitian, free-fermion effective Hamiltonian:

\begin{align}
H_{\mathrm{eff}}&= \sum_{j=1}^{n_z+n_Z}\left(\eps_j+\rmi\theta\right) a^\dagger_j a_j \label{eq:h_eff}
\end{align}

for some fermionic modes $a_j$. Then the spectrum of $E_\mathbb{I}$ coincides with the spectrum of $\rme^{-H_{\mathrm{eff}}}$. It would be interesting to determine whether there exist choices of $a_j$ such that $E_\mathbb{I}$ is equal to
and, if so, whether the form of the $a_j$ can be determined straightforwardly from properties of the function $f(z)$. Such a free-fermion form of $H_{\mathrm{eff}}$ in the XY model[Except on the disorder line $a=b$, this model is outside the subclass analysed here. The MPS and transfer matrix that correspond to $H_\mathrm{eff}$ are exact only in the limit of infinite bond dimension. ] with $f(z) = \frac{1}{z}(z-a)(z-b)$ is analysed analytically for the case of $a,b\in\mathbb{R}$ in <cit.>. The zeros are restricted to be real so that the quantum to classical mapping can be used <cit.>. For the case of complex zeros, the eigenvalues of the transfer matrix were studied numerically in <cit.>, but no formula for $H_\mathrm{eff}$ is given.

As shown in the proof of Corollary <ref>, the phase factor $\rme^{\rmi\theta} = \pm 1$ is the phase difference between $E_\mathbb{I}$ and $E_Z$. In particular, whereas the dominant eigenvalue of $E_\mathbb{I}$ is always $1$ (by normalisation), the dominant eigenvalue of $E_Z$ is $\rme^{\rmi\theta}$. This is a $\mathbb Z_2$ invariant of the gapped phase of matter <cit.>. Its value is independent of the topological invariant $\omega$; indeed, even when $\omega = 0$, one can have $\theta = 0$ or $\theta = \pi$. These are known as two distinct symmetry-protected trivial (SPt) phases of matter <cit.>: whilst they are not accompanied by protected edge modes or degeneracies in the entanglement spectrum, they are nevertheless separated by a quantum phase transition.

Physically, the SPt phases associated to $\theta=0$ or $\theta=\pi$ are distinguished by noting whether or not the string order parameter with long-range order oscillates, i.e., $\langle \mathcal O_\omega(1) \mathcal O_\omega(N+1) \rangle = \textrm{constant} \times \rme^{\rmi\theta N} (1+o(1)) $. Hence, one can interpret $\theta$ as the momentum of the string order parameter. Intuitively, the reason this is quantized (i.e., that it has to take on a discrete value $\theta \in \{0,\pi\}$) is due to the string order parameter being related to a $\mathbb Z_2$ symmetry, implying that $\rme^{2\rmi\theta} = 1$. As discussed in Remark <ref>, this can be made precise in the MPS framework by noting that $E_{Z^2} = E_\mathbb{I}$.

As derived in Remark <ref>, all models with $n_P = n_z +n_Z-k$ for $k>0$ are limiting cases of models with $n_P = n_z' +n_Z$, where $n_z' = n_z+k$ and we have $k$ additional zeros $\eps x_j$ with $x_j$ pairwise distinct. Suppose also that the zeros of these models (including our chosen $x_j$) satisfy the strongly generic condition, and moreover that $n_P$ (and therefore $\omega$) is even along this path. Then all models with $\eps>0$ have a transfer matrix that can be analysed using Corollary <ref>; this means that the MPS is injective and that we can identify the spectrum of the transfer matrix as $2^{n_z+n_Z+k}$ eigenvalues that are products of the zeros and inverse zeros. In the limit $\eps=0$, the limiting MPS is the ground state of a model with $n_P = n_z +n_Z-k$.
The corresponding limiting transfer matrix has $2^{n_z+n_Z}$ known non-zero eigenvalues, and has $2^{n_z+n_Z+k}-2^{n_z+n_Z}$ zero eigenvalues (arising from products of the form (<ref>) for sets that contained any $\eps x_j$).
Note that a priori this limiting MPS need not be in canonical form; bringing it to canonical form can potentially reduce the bond dimension, thereby removing zero eigenvalues of the transfer matrix[However, our conjecture for generic models implies that the limiting MPS does have optimal bond dimension—see the discussion below Corollary <ref>.].
An analogous argument can be made for $k<0$ by considering $f(1/z)$ and noting that the models defined by $f(1/z)$ and $f(z)$ differ by an on-site change of basis and so have the same transfer matrix.

§.§.§ Proof of Corollary <ref>

The proof essentially follows from the proof of Corollary <ref>. In particular, given our assumptions we have that $\omega$ is even and so, as explained above, we have $2^{n_z+n_Z-1}$ distinct eigenvalues of $E_\mathbb{I}$ labelled by odd-size sets $s\subseteq S$:

\begin{align}
\mu(s) = \sigma (-1)^{n_P+1}\prod_{\tau_j \in s}   \tau_j  = -\sigma \prod_{\tau_j \in s}   \tau_j.
\end{align}

(Recall that if $\sigma=-1$ in (<ref>), we have an additional factor of $(-1)^N$ in all correlators given in Theorem <ref>.) Furthermore, we also have $2^{n_z+n_Z-1}$ distinct eigenvalues of $E_Z$, labelled by the even-size subsets:

\begin{align}
\tilde\mu(s) = -\sigma \prod_{\tau_j \in s}   \tau_j.
\end{align}

Moreover, by the strongly generic condition, $\lvert \tilde\mu(s)\rvert \neq \lvert \mu(s')\rvert$ for any $s,s'$.
Now, the upper bound on the bond dimension means that we have identified all eigenvalues of $E_\mathbb{I}$ up to phase factors, and since all $\lvert \tau_j \rvert<1$ we have a unique largest eigenvalue: this is $ \tilde\mu(s) =-\sigma$, where $s$ is the empty set. Hence, the MPS is injective.
As discussed in Remark <ref>, having an injective MPS implies that the spectrum of $E_Z$ is either the same as the spectrum of $E_\mathbb{I}$ or the same as the spectrum of $-E_{\mathbb{I}}$. We can then conclude that $\pm \tilde\mu(s)$ is an eigenvalue of $E_\mathbb{I}$ for all $s$ (where the sign does not depend on $s$). We can identify this sign since we know that $\mu=1$ is an eigenvalue of $E_\mathbb{I}$. Hence, if $\sigma=-1$, we have that $\tilde\mu(s)$ is an eigenvalue of $E_\mathbb{I}$ for all $s$, while if  $\sigma=1$, then $-\tilde\mu(s)$ is an eigenvalue of $E_\mathbb{I}$ for all $s$. We have hence identified $2^{n_z+n_Z}$ eigenvalues of $E_\mathbb{I}$, and by the upper bound on the bond dimension this means we have the full spectrum.

(nodeA1) at (5.5,1) [thick,minimum size=.15cm,circle,draw] $A$;
[rounded corners] (node cs:name=nodeA1, angle = 180)–+(-.5,0.);
[rounded corners] (node cs:name=nodeA1, angle = 0) –+(.5,0.);
(nodeA2) at (5.5,-1) [thick,minimum size=.15cm,circle,draw] $\overline{A}$;
[rounded corners] (node cs:name=nodeA2, angle = 180)–+(-.5,0.);
[rounded corners] (node cs:name=nodeA2, angle = 0) –+(.5,0.);
[rounded corners] (node cs:name=nodeA1, angle = -90) – (node cs:name=nodeA2, angle = 90);

at (7,0) $=$;

(noder) at (8,0) [thick,minimum size=.15cm,circle,draw] $r_1$;
[rounded corners] (node cs:name=noder, angle = 90) – +(0,.5)–+(-.75,0.5);
[rounded corners] (node cs:name=noder, angle = -90) – +(0,-.5)–+(-.75,-0.5);
(nodel) at (9,0) [thick,minimum size=.15cm,circle,draw] $l_1$;
[rounded corners] (node cs:name=nodel, angle = 90) – +(0,.5)–+(.75,0.5);
[rounded corners] (node cs:name=nodel, angle = -90) – +(0,-.5)–+(.75,-0.5);

at (1,0) $-\frac{1}{b}  \times \frac{1}{C_X}$;

(nodeE) at (2.5,0) [thick,minimum size=.15cm,circle,draw] $E_X$;
(nodev) at (3.5,0) [thick,minimum size=.15cm,circle,draw] $r_1$;
[rounded corners] (node cs:name=nodeE, angle = 70) – +(0,.5)–+(.8,0.5)– (node cs:name=nodev, angle = 90);
[rounded corners] (node cs:name=nodeE, angle = -70) – +(0,-.5)–+(.8,-0.5)– (node cs:name=nodev, angle = -90);
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.5)–+(-.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-.75,-0.5);
(nodeE) at (5.5,0) [thick,minimum size=.15cm,circle,draw] $E_X$;
(nodev) at (4.5,0) [thick,minimum size=.15cm,circle,draw] $l_1$;
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.5)–+(-.8,0.5)– (node cs:name=nodev, angle = 90);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-.8,-0.5)– (node cs:name=nodev, angle = -90);
[rounded corners] (node cs:name=nodeE, angle = 90) – +(0,.5)–+(.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -90) – +(0,-.5)–+(.75,-0.5);

at (1,0) $-a\times\frac{1}{C_Y}$;

(nodeE) at (2.5,0) [thick,minimum size=.15cm,circle,draw] $E_Y$;
(nodev) at (3.5,0) [thick,minimum size=.15cm,circle,draw] $r_1$;
[rounded corners] (node cs:name=nodeE, angle = 70) – +(0,.5)–+(.8,0.5)– (node cs:name=nodev, angle = 90);
[rounded corners] (node cs:name=nodeE, angle = -70) – +(0,-.5)–+(.8,-0.5)– (node cs:name=nodev, angle = -90);
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.5)–+(-.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-.75,-0.5);
(nodeE) at (5.5,0) [thick,minimum size=.15cm,circle,draw] $E_Y$;
(nodev) at (4.5,0) [thick,minimum size=.15cm,circle,draw] $l_1$;
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.5)–+(-.8,0.5)– (node cs:name=nodev, angle = 90);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-.8,-0.5)– (node cs:name=nodev, angle = -90);
[rounded corners] (node cs:name=nodeE, angle = 90) – +(0,.5)–+(.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -90) – +(0,-.5)–+(.75,-0.5);

at (1,0) $+\frac{a}{b}\times\frac{1}{C_{\tilde{Z}}}$;

(nodeE) at (2.5,0) [thick,minimum size=.15cm,circle,draw] $E_{\tilde Z}$;
(nodev) at (3.5,0) [thick,minimum size=.15cm,circle,draw] $r_1$;
[rounded corners] (node cs:name=nodeE, angle = 70) – +(0,.5)–+(.8,0.5)– (node cs:name=nodev, angle = 90);
[rounded corners] (node cs:name=nodeE, angle = -70) – +(0,-.5)–+(.8,-0.5)– (node cs:name=nodev, angle = -90);
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.5)–+(-.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-.75,-0.5);
(nodeE) at (5.5,0) [thick,minimum size=.15cm,circle,draw] $E_{\tilde{Z}}$;
(nodev) at (4.5,0) [thick,minimum size=.15cm,circle,draw] $l_1$;
[rounded corners] (node cs:name=nodeE, angle = 110) – +(0,.5)–+(-.8,0.5)– (node cs:name=nodev, angle = 90);
[rounded corners] (node cs:name=nodeE, angle = -110) – +(0,-.5)–+(-.8,-0.5)– (node cs:name=nodev, angle = -90);
[rounded corners] (node cs:name=nodeE, angle = 90) – +(0,.5)–+(.75,0.5);
[rounded corners] (node cs:name=nodeE, angle = -90) – +(0,-.5)–+(.75,-0.5);

Graphical representation of the transfer matrix for $f(z)=z^{-2}(z-a)^2(z-b)^2$ with $a<1$ and $b>1$. We define $\tilde{Z} = Z - \frac{\left(a^2-1\right) b^2-a b+1}{b (a-b)}$. The normalisation constants $C_X$, $C_Y$ and $C_{\tilde{Z}}=-C_XC_Y$ are given in (<ref>).

§.§ Diagonalising the transfer matrix for $f(z) = z^{-2}(z-a)^2(z-b)^2$

We now explore how our results can allow us to constrain the eigenvectors of the transfer matrix as well as the eigenvalues. Let us consider $f(z) = z^{-2}(z-a)^2(z-b)^2$ with $a<1$ and $b>1$. The correlators are non-zero for $-1\leq \alpha \leq 1$. The values of these correlators are given in Example <ref> (after shifting the index appropriately). Let $P_1=\ket{r_1}\bra{l_1}$, then, given the upper bound $\chi^2=4$, in Appendix <ref> we prove the following:

\begin{align}
E_\mathbb{I} &=P_1 - \nonumber
a \left( \frac{(b-a)}{{(1-b^{-2})(1-ab)}}E_YP_1E_Y \right)-\frac{1}{b} \left(\frac{(a-b)}{{(1-a^2)(1-ab)}}E_XP_1E_X\right) \\&\qquad+\frac{a}{b}  \left( \frac{(a-b)^2}{{(1-a^2)(1-b^{-2})(1-ab)^2}}E_{\tilde Z}P_1E_{\tilde Z} \right),\label{eq:diagonaltransfer}
\end{align}

where $\tilde{Z} = Z - \frac{\left(a^2-1\right) b^2-a b+1}{b (a-b)}$. The eigenvalues are $\{1,-a,-1/b,a/b\}$ and the bracketed expressions are the relevant $P_j=\ket{r_j}\bra{l_j}$. This is given graphically in Figure <ref>. Note that we did not need the MPS tensor to derive this result, although it is needed to evaluate this formula. We can furthermore use the formula to identify that, say, $\ket{r_2} \propto E_Y \ket{r_1}$. Given this formula for $E_\mathbb{I}$, further calculations are required to determine whether we can find expressions for fermionic modes $a_j$ in (<ref>) such that $E_\mathbb{I}=\rme^{-H_{\mathrm{eff}}}$.

One can straightforwardly generalise the analysis of Appendix <ref> to derive results about the eigenvectors of $E_\mathbb{I}$ in cases where $f(z)$ has more zeros. However, this will not immediately lead to a diagonal form as in (<ref>). In that example, operators such as $E_YP_1E_Y$ corresponded to a unique eigenvalue, while in the general case such operators can correspond to multiple eigenvalues (i.e., we have a sum over more than one set of zeros in the expressions given in Theorem <ref>).

§ RECOVERING GENERIC MODELS VIA A LIMITING PROCEDURE

The results in this paper are for the special subclass of BDI and AIII models where all non-zero zeros of $f(z)$ have an even multiplicity. This subclass admits closed results for correlation functions and an exact matrix product state representation of the ground state. In this section, we show how a generic BDI or AIII model can be obtained as a limit[While it will be a sequence $f_m(z)$ of finite-range models, the range increases linearly with $m$ (i.e., it increases without bound). Nevertheless, the limit $\lim_{m \to \infty} f_m(z)$ is finite-range.] of such special models. This shows that any ground state in these subclasses can be obtained as a limit of matrix product states. Moreover, such a sequence can be used to derive results about generic models by using results derived in the present work.

We note that sequences approximating certain Toeplitz symbols (see Section <ref>) by rational functions are discussed in <cit.>. The methods of expansion are different in general, although coincide for the case of the quantum Ising model. Moreover, our approach focuses on a sequence of Hamiltonians with MPS ground states, from which one can then derive a sequence of Toeplitz symbols.

§.§ A sequence of models

A general gapped BDI or AIII model corresponds, up to normalisation, to a Laurent polynomial of the form

\begin{align}
f(z) =  \sigma \frac{1}{z^{N_P}} \prod_{j=1}^{N_z} \left(z-z_j\right)\prod_{k=1}^{N_Z} \left(z-Z_k\right). \label{fz_generic}
\end{align}

We have $\lvert z_j\rvert <1$, $\lvert Z_k\rvert >1$, and $\sigma \in \{\pm1\}$ for the BDI class <cit.>, while for AIII, $\sigma=\rme^{\rmi\theta}$.
For a generic model in this class, we have that all zeros have multiplicity one, and moreover that no zeros coincide with inverse zeros (or inverse conjugate zeros in AIII). We now construct a sequence $f_m(z)$ ($m=1,2,\dots$) of models of the form (<ref>) such that in an open region of the complex plane containing the unit circle, we have $\lim_{m \to \infty} f_m(z) = f(z)$.  As discussed in Section <ref>, the reason that the unit circle is special is because it stores all physical information. To construct the sequence, we first rewrite

\begin{align}
f(z) = \sigma \left(\prod_{k=1}^{N_Z} \left(-Z_k\right)\right) \frac{1}{z^{N_P-N_z}} g(z)^2 \quad \textrm{with } g(z) = \prod_{j=1}^{N_z} \sqrt{ 1-\frac{z_j}{z}} \prod_{k=1}^{N_Z} \sqrt{  1-\frac{z}{Z_k} }.
\end{align}

By using the series expansion of the square roots, we obtain a well-defined series expansion for $g(z)$ for $\max_j\{|z_j|\} < |z| < \min_k \{|Z_k| \}$ (i.e., an annulus containing the unit circle). More precisely, if we define the partial sums

\begin{align}
s_m(z) = \sum_{l=0}^m (-1)^l \left( \begin{array}{c} 1/2 \\ l \end{array} \right) z^l
\end{align}

then $\sqrt{1-z} = \lim_{m\to \infty} s_m(z)$ (where we take the principal branch of the square root).
We can thus define a sequence

\begin{align}
g_m(z) = \prod_{j=1}^{N_z} s_m(z_j/z) \prod_{k=1}^{N_Z} s_m(z/Z_k) \label{eq:path}
\end{align}

where $g_m(z)$ converges uniformly to $g(z)$ on the annulus given above.
Hence, $f_m(z) \propto z^{N_z-N_P} g_m(z)^2$ defines a sequence of polynomials which all belong to the special subclass of models studied in the present work, and where the limit of this sequence is a generic gapped BDI or AIII model[One can analyse non-generic models in a similar way, for ease of presentation we focus on the generic case. For example, if $f(z) =(z-z_1)(z-z_2)^2 $, we could take the sequence $g_m(z)=s_m(z_1/z)(z-z_2)$.]. We explain below that every truncated model has the same winding number as the limiting $f(z)$.

This path can then be used to extend certain results of the present work to the more general case. We illustrate this now for the order parameter in the BDI class.

§.§ Order parameter

In a previous work <cit.>, for a general $f(z)$ of the form (<ref>) we derived the following value of the order parameter in class BDI:

\begin{align}
\lim_{N\rightarrow\infty}\lvert\langle\mathcal{O}_\omega(1)\mathcal{O}_\omega(N+1)\rangle \rvert = \left( \frac{ \prod_{j_1,j_2=1}^{N_z} (1-z_{j_1} z_{j_2})
\prod_{k_1,k_2=1}^{N_Z} \left(1-\frac{1}{Z_{k_1} Z_{k_2}}\right)}{\prod_{j=1}^{N_z}\prod_{k=1}^{N_Z} \left (1-\frac{z_j}{ Z_k}\right)^2} \right)^{1/4}. \label{eq:orderparameter_generic}
\end{align}

This result is proved using Szegő's theorem <cit.>. To illustrate the usefulness of approximating a generic $f(z)$ by a sequence of models with degenerate roots, as in (<ref>), here we demonstrate how this expression (<ref>) is naturally obtained by taking a limit of the formula (<ref>) which is derived in this work using Day's formula (Theorem <ref>).

For any $m$, since $s_m(z)$ is a polynomial with $s_m(0) = 1$, we know there exists a set of complex numbers $\left\{ \lambda_l^{(m)}\right\}_{l=1,\cdots,m}$ such that we can write $s_m(z) = \prod_{l=1}^m \left(1-\lambda_l^{(m)} z \right)$. Moreover, since $s_m(z)$ is real-valued on the real-line, the $\lambda_l^{(m)}$ are real or come in complex-conjugate pairs. One can also show[It is straightforward to see that $\lvert 1-s_m(z) \rvert <1$ for $z=\rme^{\rmi k}$, and the result follows using Rouché's theorem <cit.>.] that $\lvert\lambda_l^{(m)}\rvert<1$ for all $m$. This means that the winding number of $g_m(z)$ is zero, and so, for all $m$, $f_m(z)$ has winding number $N_z-N_P$.

We thus have

\begin{align}
\sqrt{1-z} = \lim_{m \to \infty} \prod_{l=1}^m \left(1-\lambda_l^{(m)} z\right) \qquad \lvert z\rvert<1. \label{eq:sqrt_identity}
\end{align}

By taking the square root of both sides of (<ref>) we obtain the useful identity

\begin{align}
(1-z)^{1/4} = \lim_{m \to \infty} \prod_{l_1,l_2=1}^m \left( 1 - \lambda_{l_1}^{(m)} \lambda_{l_2}^{(m)} z \right), \label{eq:sqrt_identity2}
\end{align}

we prove this in Appendix <ref>.

Using the above factorisation of $s_m(z)$,  we have that

\begin{align}
g_m(z) = \prod_{j=1}^{N_z} \prod_{l=1}^m \left(1-\frac{\lambda_l z_j}{z} \right) \times \prod_{k=1}^{N_Z} \prod_{p=1}^m \left(1- \frac{\lambda_p z}{Z_k} \right).
\end{align}

There are $n_z = N_z m$ roots inside the unit circle given by $\lambda_l z_j$ and $n_Z = N_Z m$ roots outside given by $Z_k/\lambda_p$.

Using the result (<ref>) derived in this work, the order parameter for a given $m$ is

\begin{align}
\lim_{N\rightarrow\infty}\lvert\langle\mathcal{O}_\omega(1)\mathcal{O}_\omega(N+1)\rangle \rvert = \frac{ \prod_{j_1,j_2=1}^{N_z} \prod_{l_1,l_2=1}^m  (1-\lambda_{l_1} \lambda_{l_2} z_{j_1} z_{j_2}) \prod_{k_1,k_2=1}^{N_Z} \prod_{l_1,l_2=1}^m  \left(1-\frac{\lambda_{l_1} \lambda_{l_2}}{Z_{k_1} Z_{k_2}}\right)}{\prod_{j=1}^{N_z}\prod_{k=1}^{N_Z} \prod_{ l_1, l_2=1}^m  \left (1-\lambda_{ l_1} \lambda_{ l_2} \frac{z_j}{ Z_k}\right)^2}.
\end{align}

Note that, as discussed in Remark <ref>, this formula applies even in the non-generic case. Using (<ref>) to take the limit, we obtain

\begin{align}
\lim_{m \to \infty} \lim_{N\rightarrow\infty}\lvert\langle\mathcal{O}_\omega(1)\mathcal{O}_\omega(N+1)\rangle \rvert = \left( \frac{ \prod_{j_1,j_2=1}^{N_z}  (1- z_{j_1} z_{j_2} )
\prod_{k_1,k_2=1}^{N_Z} \left(1-\frac{1}{Z_{k_1} Z_{k_2}}\right)}{\prod_{j=1}^{N_z}\prod_{k=1}^{N_Z}  \left (1- \frac{z_j}{ Z_k}\right)^2} \right)^{1/4}. \label{eq:orderparameter_limit}
\end{align}

We see that this coincides with the known formula (<ref>). For this to be a rigorous independent derivation of this formula, one should also prove that the two limits on the left-hand side of (<ref>) commute. This would require bounding the subleading terms which are also given by Theorem <ref>, but we will not pursue this further here.

§ TOEPLITZ DETERMINANTS

An $N\times N$ Toeplitz matrix takes the following form:

\begin{align}
(T_{mn}) = ({t}_{m-n}) = \left(\begin{array}{ccccc}{t}_0 & {t}_{-1} & {t}_{-2} & \hdots & {t}_{-(N-1)} \\{t}_{1} & {t}_0 & {t}_{-1} & \hdots &{t}_{-(N-2)} \\{t}_2 & {t}_{1} & {t}_0 & \hdots&{t}_{-(N-3)} \\\vdots & \vdots & \vdots  & \ddots&\vdots\\ {t}_{N-1} &{t}_{N-2} & {t}_{N-3} & \hdots & {t}_0 \end{array}\right).\label{toeplitzform}
\end{align}

We say the Toeplitz matrix has (scalar) symbol $t(z)$ if the entries of the Toeplitz matrix are the Fourier coefficients of the function $t(z)$, i.e.,

\begin{align}
t_{n} = \frac{1}{2\pi} \int^{2 \pi}_{0}t\left(\rme^{\rmi k}\right) \rme^{- \rmi n k} \rmd k. \label{fourier}
\end{align}

Given this symbol, we can generate an $N \times N$ Toeplitz matrix for any $N$. We define $D_N[t(z)]$ to be the determinant of the $N\times N$ Toeplitz matrix generated by the symbol $t(z)$. Toeplitz determinants have a rich theory; see, for example, the review <cit.> and references therein. We will also need to consider block Toeplitz matrices—generalising the above to symbols that are matrix valued functions. Our cases of interest have a $2\times2$ matrix symbol $\Phi(z)$, where the determinant $D_N[\Phi(z)]$ will be the determinant of the $2N \times 2N$ matrix with form as in (<ref>), with $2\times 2$ matrix Fourier coefficients defined as in (<ref>).

Now, for translation-invariant free-fermion chains in class BDI, we have that the correlator of string operators is the following Toeplitz determinant:

\begin{align}
(-1)^{N(\alpha-1)}\langle\mathcal{O}_\alpha(1) \mathcal{O}_\alpha(N+1)\rangle  &= \det \left( \frac{1}{2\pi} \int_0^{2\pi} \frac{ {f(\rme^{\rmi k})}}{\big\lvert f(\rme^{\rmi k})\big\rvert}\rme^{-\rmi\alpha k} \rme^{-\rmi(m-n) k}\rmd k\right)_{m,n=1}^{N}\nonumber\\\label{detcorr}
&= D_N\left[ \sqrt{\frac{f(z)}{f(1/z)}}z^{-\alpha}\right]= D_N\left[ \sqrt{\frac{f(1/z)}{f(z)}}z^{\alpha}\right],\end{align}

(the final equality is $\det M = \det M^T$). Since the non-trivial fermionic two-point function is given by

\begin{align}
\langle -\rmi\tilde\gamma_n \gamma_m \rangle = \int_0^{2\pi} \frac{ {f(\rme^{\rmi k})}}{\big\lvert f(\rme^{\rmi k})\big\rvert} \rme^{-\rmi(m-n) k}\rmd k,
\end{align}
the formula (<ref>) follows from Wick's theorem—a derivation is given in <cit.>. At $z=1$, $f(z)$ is real and we choose the branch of the square-root in (<ref>) so that $\sqrt{f(z)/f(1/z)}$ is equal to the sign of $f(1)$. For $f(z)$ defined in (<ref>), the scalar symbol can be simplified to a rational function. The Fourier coefficients that give the matrix entries are evaluated in Appendix <ref>, while Day's formula gives us an expression for the determinant in this case <cit.>. Day's formula may also be applied to evaluate string correlations in class AIII, this is derived in Appendix <ref> where we find:

\begin{align}
\langle\tilde{\mathcal{O}}_\alpha(1) \tilde{\mathcal{O}}_\alpha(N+1)\rangle
&= \big\lvert D_N[t(z)]\big\rvert^2=\big\lvert D_N[t(1/z)]\big\rvert^2,\end{align}

where $t(z)= \sqrt{{f(z)}/{\overline{f}(1/z)}}z^{-\alpha}$. Again this can be simplified to a rational function in our case of interest.

Consider symbols of the form:

\begin{align}
t(z)= \rho  \frac{\prod_{j=1}^{s} \left(z-\tau_j\right)}{\prod_{j=1}^{q} \left(1-z/\gamma_j\right)\prod_{j=1}^{p} \left(z-\delta_j\right)}\label{daycanon}
\end{align}

with  $p\geq 1, q \geq 1$, $s \geq p+q$, $|\delta_j|<1$, $|\gamma_j|>1$ and $\{\tau_1, \dots, \tau_s\}$ pairwise distinct. Then we have:

\begin{align}
D_N[t(z)] = \sum_M C_M^{\vphantom N} r_M^N \qquad N\geq 1,
\end{align}

where the sum is taken over all subsets of $M \subseteq \{1, \dots,s\}$ such that $\lvert M \rvert = p$ and $C_M$, $r_M$ are defined as follows. For each $M$ define $M^c =  \{1, \dots, s\} \setminus M$, and also define $P=\{ 1 \dots p\}$ and $Q= \{ 1 \dots q\}$. Then

\begin{align}
r_M&= (-1)^{s-p} \rho \prod_{k \in M^c} \tau_k \\
C_M &= \prod_{k \in M^c,m\in P} (\tau_k - \delta_m) \prod_{l \in Q, j\in M} (\gamma_l - \tau_j) \prod_{l \in Q, m\in P} (\gamma_l - \delta_m)^{-1} \prod_{k \in M^c,j\in M} (\tau_k - \tau_j)^{-1}  .
\end{align}

The condition that the $\tau_j$ are pairwise distinct is necessary because differences $\tau_k - \tau_j$ appear in the denominator of $C_M$. However, we will want to study cases where the $\tau_j$ are degenerate, this is important since the factor $z^k = (z-0)^k$ appears in our analysis of correlators. As pointed out in <cit.>, one can still use Day's formula by taking an appropriate limit. In particular, let us write  $z^k = \lim_{\eps\rightarrow 0} \prod_{j=1}^k (z-\eps x_j)$ for $x_j$ pairwise distinct. For finite $\eps$ we do not have degenerate zeros and can apply Day's formula. Moreover, the Fourier coefficients, and hence the matrix elements of the Toeplitz determinant will depend continuously on $\eps$. Indeed the Fourier coefficients of $\prod_{j=1}^k (z-\eps x_j)t(z) $ will just be sums over shifted Fourier coefficients of $t(z)$. Since we are taking a finite Toeplitz determinant, we can then take the limit of Day's formula as $\eps \rightarrow 0$ and this will be the determinant we are interested in.

We are also interested in the eigenvalues of the correlation matrix. For translation-invariant free-fermion chains in class BDI, the correlation matrix for a subsystem of size $N$ is a $2N\times 2N$ block Toeplitz matrix generated by the symbol $\Phi(z,0)$ <cit.>, where

\begin{align}
\Phi(z,\lambda) = \left(\begin{array}{cc}\rmi \lambda & \sqrt{\frac{f(z)}{f(1/z)}} \\ -\sqrt{\frac{f(1/z)}{f(z)}} & \rmi \lambda \end{array}\right). \label{blocksymbol}
\end{align}

By evaluating the determinant generated by this symbol for general $\lambda$, we can find the eigenvalues of the correlation matrix. Note that changing the branch of the square-root takes $\Phi(z,\lambda)\rightarrow-\Phi(z,-\lambda)$; since the eigenvalues come in $\pm \lambda$ pairs, they are independent of this branch choice. Moreover, notice that $\Phi(1/z,\lambda)=-\Phi(z,-\lambda)$. Again, since the eigenvalues come in $\pm \lambda$ pairs this means that this determinant is unchanged under $f(z) \rightarrow f(1/z)$. Note that for $\lambda=0$ we have, following <cit.>, that:

\begin{align}
D_N[\Phi(z,0)]&=(-1)^N D_N\left[\sqrt{{f(z)}/{f(1/z)}}\right]D_N\left[-\sqrt{{f(1/z)}/{f(z)}}\right]\nonumber\\&=\Big(\langle \mathcal{O}_0(1) \mathcal{O}_0(N+1) \rangle\Big)^2.\label{eq:detAN2}
\end{align}

In our analysis of the entanglement spectrum in class AIII, a closely related block Toeplitz matrix arises, the symbol is given in Appendix <ref>.

As before, in the models that we consider these symbols simplify to a matrix where each element is a rational function. For a certain class of such rational matrix functions, we have Gorodetsky's formula <cit.>. The general statement is too long to include here, but let us give an idea of the result:

Suppose $a(z) = \sum_{j=0}^s a_j z^j $ is a matrix polynomial, where $a_j \in \mathbb{C}^{r \times r}$ and suppose that $a_s$ is invertible. For symbols of the form:

\begin{align}
\Phi(z)=  \frac{1}{\prod_{j=1}^{p} \left(1-z/\gamma_j\right)\prod_{j=1}^{q} \left(z-\delta_j\right)} a(z)
\end{align}

where $|\delta_j|<1$ and $|\gamma_j|>1$, and subject to some conditions involving $p,q,s$ and the matrix $a_0$, then:
\begin{align}
D_N[\Phi(z)] =(\det{a_s})^N \frac{\det{\mathcal{M}[N,\Phi(z)]}}{\det{\mathcal{M}[0,\Phi(z)]}} \qquad N\geq 1,
\end{align}

where $\mathcal{M}[n,\Phi(z)]$ is an $rs \times rs$ matrix such that $\det{\mathcal{M}[0,\Phi(z)]}\neq 0$. The construction of $\mathcal{M}$ can be found in <cit.> and requires the Smith canonical form of $a(z)$.
When we apply this theorem we will give all details for the particular cases we need.

§ STRING CORRELATORS—ANALYSIS

As stated above, the BDI string correlators satisfy

\begin{align}(-1)^{N(\alpha-1)}\langle\mathcal{O}_\alpha(1) \mathcal{O}_\alpha(N+1)\rangle =D_N[t_1(z)]=D_N[t_1(1/z)] \end{align}
where $t_1(z)= \sqrt{{f(z)}/{f(1/z)}}z^{-\alpha}$. Moreover, in AIII the string correlators satisfy

\begin{align}\big\langle\tilde{\mathcal{O}}_\alpha(1) \tilde{\mathcal{O}}_\alpha(N+1)\big\rangle =\big\lvert D_N[t_2(z)]\big\rvert^2=\big\lvert D_N[t_2(1/z)]\big\rvert^2 \end{align}

where $t_2(z)= \sqrt{{f(z)}/{\overline{f}(1/z)}}z^{-\alpha}$. We hence see that we can understand both classes by analysing $D_N[t(z)]$ for $t(z)=t_2(z)$. Then for class BDI we reach Theorem <ref> by noting that $\overline{f}(1/z)=f(1/z)$ for $t_\alpha\in\mathbb{R}$ and, since zeros come in complex-conjugate pairs, we can simplify certain formulae. For class AIII we reach Theorem <ref> by simply taking the absolute value squared of the result.

We will show now that for

\begin{align}
f(z) = \frac{\rme^{\rmi\theta}}{z^{n_P}} \prod_{j=1}^{n_z} \left(z-z_j\right)^2\prod_{k=1}^{n_Z} \left(z-Z_k\right)^2, \label{fzcanon2}
\end{align}

then either $D_N[t(z)]$ or $D_N[t(1/z)]$ can be evaluated using Day's formula.
Some further analysis is then needed to reach our results. First, note that we can simplify by setting $\theta =0$. In the BDI case, $\theta\in\{0,\pi\}$; if $\theta=\pi$ we take the other branch of the square-root and multiply the Toeplitz matrix by $-1$. This sign hence gives an additional factor of $(-1)^N$ when we evaluate the determinant, as mentioned in Remark <ref>. In the AIII case $\theta\in[0,2\pi)$, this would give an additional complex phase when we evaluate the determinant, but since we take the absolute value this drops out.
If $n_z=0$ and $n_Z=0$ then $f(z)=z^{-n_P}$ and the string correlators are trivial. We treat the case $n_z>0$ and $n_Z>0$ first. Then we analyse the remaining cases where either $n_z=0$ or $n_Z=0$.

§.§ Proof of Theorems <ref> and <ref> with zeros inside and outside the unit circle

Suppose first that $n_z>0$ and $n_Z>0$. By inserting (<ref>) into $t(z)$ and rearranging we reach:

\begin{align}
t(z) &=  z^{n_z+n_Z-n_P-\alpha}\left(\prod_{j=1}^{n_Z}(-\overline{Z}_j)\right)^{-1} \frac{\prod_{j=1}^{n_z} \left(z-z_j\right)\prod_{k=1}^{n_Z} \left(z-Z_k\right)}{\prod_{j'=1}^{n_z} \left(1-z/\overline{z}_{j}^{-1}\right)\prod_{k=1}^{n_Z} \left(z-\overline{Z}_{k'}^{-1}\right)}\label{tdef}\\