Skeleton of Matrix-Product-State-Solvable Models Connecting Topological Phases of Matter

Models whose ground states can be written as an exact matrix product state (MPS) provide valuable insights into phases of matter. While MPS-solvable models are typically studied as isolated points in a phase diagram, they can belong to a connected network of MPS-solvable models, which we call the MPS skeleton. As a case study where we can completely unearth this skeleton, we focus on the one-dimensional BDI class -- non-interacting spinless fermions with time-reversal symmetry. This class, labelled by a topological winding number, contains the Kitaev chain and is Jordan-Wigner-dual to various symmetry-breaking and symmetry-protected topological (SPT) spin chains. We show that one can read off from the Hamiltonian whether its ground state is an MPS: defining a polynomial whose coefficients are the Hamiltonian parameters, MPS-solvability corresponds to this polynomial being a perfect square. We provide an explicit construction of the ground state MPS, its bond dimension growing exponentially with the range of the Hamiltonian. This complete characterization of the MPS skeleton in parameter space has three significant consequences: (i) any two topologically distinct phases in this class admit a path of MPS-solvable models between them, including the phase transition which obeys an area law for its entanglement entropy; (ii) we illustrate that the subset of MPS-solvable models is dense in this class by constructing a sequence of MPS-solvable models which converge to the Kitaev chain (equivalently, the quantum Ising chain in a transverse field); (iii) a subset of these MPS states can be particularly efficiently processed on a noisy intermediate-scale quantum computer.

Models whose ground states can be written as an exact matrix product state (MPS) provide valuable insights into phases of matter. While MPS-solvable models are typically studied as isolated points in a phase diagram, they can belong to a connected network of MPS-solvable models, which we call the MPS skeleton. As a case study where we can completely unearth this skeleton, we focus on the one-dimensional BDI class-non-interacting spinless fermions with time-reversal symmetry. This class, labelled by a topological winding number, contains the Kitaev chain and is Jordan-Wigner-dual to various symmetry-breaking and symmetry-protected topological (SPT) spin chains. We show that one can read off from the Hamiltonian whether its ground state is an MPS: defining a polynomial whose coefficients are the Hamiltonian parameters, MPS-solvability corresponds to this polynomial being a perfect square. We provide an explicit construction of the ground state MPS, its bond dimension growing exponentially with the range of the Hamiltonian. This complete characterization of the MPS skeleton in parameter space has three significant consequences: (i) any two topologically distinct phases in this class admit a path of MPS-solvable models between them, including the phase transition which obeys an area law for its entanglement entropy; (ii) we illustrate that the subset of MPS-solvable models is dense in this class by constructing a sequence of MPS-solvable models which converge to the Kitaev chain (equivalently, the quantum Ising chain in a transverse field); (iii) a subset of these MPS states can be particularly efficiently processed on a noisy intermediate-scale quantum computer. The realization that the entanglement of gapped manybody ground states obeys an area law was a breakthrough for condensed matter physics [1]. It justifies the use of tensor network states as a description of the wavefunction, having become a key analytic and numerical tool [2][3][4][5][6][7][8][9]. These tools are most refined for the case of matrix product states (MPS) describing one-dimensional systems. In most scenarios, such MPS are approximations to the true ground states. However, a wide variety of Hamiltonians are known where the ground state is an exact MPS-i.e., with a finite bond dimension in the thermodynamic limit [10][11][12][13][14][15][16][17][18][19][20][21][22][23].
Less explored are continuous paths of MPS-solvable models connecting distinct phases of matter through a quantum phase transition. One option is to simply define paths in the manifold of exact MPS states in Hilbert space. Indeed, using the well-established parent Hamiltonian construction, this gives a path of MPS-solvable models [9,34,35]. While MPS cannot capture conformal critical points which have diverging entanglement entropy, they can describe certain multicritical points where the gap closes [14]. Paths approaching such points can exhibit a diverging correlation length in a finite bond dimension MPS. Explicit discussions in the literature of such instances seem to be rare, an example being the disorder line in the spin-1/2 XY chain [38][39][40][41][42][43]; this interpolates between two distinct ferromagnets by passing through a multicritical point 1 with dynamic critical exponent z = 2. A reincarnation of this example-related by a Kramers-Wannier transformation-is the MPS path connecting the trivial phase to the Haldane phase as realized by the cluster state [14,44]. Let us note that the aforementioned parent Hamiltonian construction is not unique and can give rise to unwieldy Hamiltonians which are not necessarily in a class of interest.
In this work, we do not start from a path of MPS: instead, we specify the Hamiltonian class and ask which models have an MPS ground state. This leads to MPSsolvable paths forming the skeleton around which the rest of the phase diagram is structured. For a particular class of non-interacting symmetric Hamiltonians (BDI class [45]), we develop a general understanding of paths of MPS-solvable models, connecting the distinct SPT phases. Different paths connect at joints where the system is multicritical and still has an MPS ground state. We refer to this network as the MPS skeleton. Remarkably, this skeleton is dense in this class (similar to how rational numbers are dense on the real line): any gapped ground state can be obtained as a sequence of Hamiltonians whose ground state is an exact MPS.
We note that while the idea of the MPS skeleton is by no means particular to non-interacting systems, this setting is an interesting case study. Despite free-fermion Hamiltonians and MPS-solvable systems both being pinnacles of solubility, they have a rich interplay: one cannot typically write the ground state of a free-fermion system as an exact MPS due to its entanglement spectrum having infinite rank 2 , and there is no analytic handle on truncating this to a particular bond dimension. This truncation has been investigated numerically: an approach for the XY model is given in Ref. [47]; and more generally there are approaches based on truncating the free-fermion correlation matrix [48,49], the 'MPO-MPS method' [50] and through Schmidt decomposition [51]. The MPS description of free-fermion states has been explored before in the context of Gaussian MPS [49,52,53]. Using this framework, Ref. [53] showed that free-fermion states admitting an exact (Gaussian) MPS representation have a correlation matrix that satisfies a certain property [9], readily applying to arbitrary dimensions. We will see that, for the BDI class, this property coincides with our characterization of the MPS skeleton. Indeed, our analysis shows that the implication also works the other way: this property is sufficient for MPS-solvability, and moreover we give an explicit construction of the ground state. We do not appeal to the formalism of Gaussian MPS, and it would be interesting to translate our results into that language, giving a concrete subclass of Gaussian MPS for which we have an explicit construction.
To briefly outline the paper, in Section II we introduce BDI Hamiltonians and summarize our main results, followed by explicit examples in Section III. We provide the derivation of our results in Section IV and elaborate on special cases in Section V.

A. Model
We consider a chain of Majorana operators γ n andγ n , which are respectively real and imaginary under timereversal 3 . The most general Hamiltonian term in the free-fermion BDI class is h n,α = iγ n γ n+α , which has the convenient property h 2 n,α = 1. Any translation-invariant BDI Hamiltonian can be written as where t α ∈ R due to hermiticity. The special cases H ω where t α = δ α,ω are stabilizer code Hamiltonians (all terms commute), and are fixed-point Hamiltonians in the phase with winding number ω ∈ Z. Note that ω = 1 gives the Kitaev chain [54]. It is convenient to encode the information of the Hamiltonian parameters in a Laurent 4 polynomial 3 More precisely, if cn is a complex fermionic operator, the Majorana operators are γn = cn + c † n andγn = i c † n − cn . Then {γn,γm} = 0 and {γn, γm} = 2δnm. 4 Note that it can contain negative powers of z if tα = 0 for some negative α.
Previous work has already shown that a multitude of physical information can be readily extracted from f (z). E.g., the single-particle spectrum is given by ε k = |f (e ik )|; the correlation length is ξ = max i {1/| ln |ζ i ||} where ζ i are the roots of f (z); the topological invariant ω = N z − N p , where N z is the number of roots inside the unit circle and N p is the degree of the pole at z = 0 [55,56]. We will refer to the ground state of f (z): this means the ground state of the related Hamiltonian. Under the usual Jordan-Wigner transformation (see Eq. (37)), the Hamiltonian terms become In particular, the fixed-point Hamiltonians H ω correspond to symmetry-breaking or SPT 5 Hamiltonians such as the Ising model H 1 = − 1 2 n X n X n+1 and cluster model H 2 = − 1 2 n X n−1 Z n X n+1 . More generally, H α is a generalized cluster model, or α-chain [57][58][59][60][61].

B. The MPS skeleton
Here, we add to this body of knowledge by characterizing when the ground state of Eq. (1) is an MPS. For our purposes, this means that we have an explicit finitedepth circuit representation for the ground state (note that the gates of this circuit will not necessarily be unitary). In Section IV C we will make the connection to the usual definition of an MPS [9] as a tensor network where the ground state is a contraction of tensors with virtual indices that have bond dimension χ. We have that the ground state is an exact MPS if f (z) is a square; more precisely: Result 1 (Existence of MPS.) If f (z) = z p g(z) 2 for p ∈ Z and g(z) = d k=0 s k z k (with s k ∈ R), then the Hamiltonian given in Eq. (1) is frustration-free. Moreover, its ground state can be exactly represented as an MPS with finite bond dimension χ. If we ensure that s 0 = 0 = s d (which one can always do by appropriately choosing p and d ≥ 0), then where range(H) is defined as the largest power of either z or 1/z in f (z) and · : R → Z is the ceiling function.
The formula for χ is for MPS which are symmetric under fermion parity (or, equivalently, spin-flip symmetry); in the case of spontaneous symmetry breaking (for the spin chain), this formula applies to the cat state, whereas the symmetry-broken state has log 2 χ = range(H)/2 . We believe that this formula for χ is generically optimal, and in certain cases this can be proved-this is discussed in Section IV C.
This gives a complete characterization of all BDI Hamiltonians with an exact MPS ground state (with finite χ in the thermodynamic limit). More precisely, the above result holds for the broader class of Hamiltonians f (z) = ±z p g(z) 2 h(z), where h(z) is any Laurent polynomial that satisfies h(1/z) = h(z), has no roots on the unit circle and has positive constant term. However, (i) the sign can easily be toggled (see Eq. (43)) and so we will henceforth consider a positive sign, and (ii) the ground state is independent of h(z). Any f (z) not of this form has correlation functions with asymptotics containing terms of the form N α exp(−N/ξ) for α / ∈ N + and ξ ≤ ∞, which cannot be captured by an MPS with a finite bond dimension 6 . These claims are proved in Appendix A.
We point out that a similar necessary condition for MPS-solvability was obtained in Ref. [52] in the context of lattice models for free (bosonic) oscillators; the argument straightforwardly extends to the fermionic case 7 . A relation was also given between the bond dimension and the interaction range of the Hamiltonian. The fermionic analogue of these Gaussian MPS are discussed in Refs. [9,53]; a characteristic of such states is that rational functions of z = e ik generate correlations by Fourier transform. Indeed, for the case under discussion, correlations are generated by f (z)/f (1/z) which, on the MPS skeleton, reduces to the rational function z p g(z)/g(1/z). Our results thus show that this characterization is sufficient as well as necessary. Moreover, our proof is constructive-we will turn to this now.

C. Construction of MPS
Excluding a measure zero set, we have a closed form for the MPS wavefunction. This is most easily described in terms of d real parameters b k=1,··· ,d that are obtained from the following recursion. Here d is the degree of g(z) and (s 0 , · · · , s d ) are its coefficients, as defined in Result 1; writing s = (s 0 , s 1 , · · · , s d ) and flip( s) = (s d , s d−1 , · · · , s 0 ) we have: 6 The case ξ = ∞ corresponds to a critical system. 7 We are grateful to N. Schuch for clarifying this point.
drop last entry of s; The outcome of this algorithm will only be used if |b k | = 1 for all k (here we thus exclude a measure zero case). Given this condition, one can show that s 0 = 0 at each step, ensuring that the ratio s k /s 0 is well-defined.
As an example of the above algorithm, consider g(z) = 1 + 4z + 2z 2 , then s = (1,4,2). From the first recursion, we obtain b 2 = 2 and s = (−3, −4). From the second recursion, we have b 1 = 4 3 and s = 7 3 . As we will now see, the values for b 1 and b 2 directly give us the ground state as a quantum circuit.
In Section IV B, we derive the following, using the same conditions as listed in Result 1 and the b k obtained through Algorithm 1: Result 2 (Construction of MPS.) If |b k | = 1 for k = 1, . . . , d, the ground state of Eq. (1) can be constructed with d layers of circuits: where |ψ p is the ground state of the fixed-point Hamiltonian H p = 1 2 n h n,p . Each layer is generated by a fixed-point Hamiltonian as follows: This circuit can be rewritten as an MPS with the bond dimension claimed in Result 1; see Section IV C.
The gates M (k) appearing in this result are generically not unitary (so one still has to normalize the wave function). However, |b k | = 1 implies they are invertible. In fact, if |b k | < 1 for k = 1, . . . d − 1 then we can give an especially efficient unitary circuit representation for the MPS whose unit element scales logarithmically with bond dimension-see Section V A.
Note that since H k is a sum of commuting terms, M (k) can be written as a product M (k) = n M (k) n as follows: This local form leads to the MPS description and will be important in our analysis below 8 .
If |b k | < 1 then β k ∈ R and so M (k) can be seen as an imaginary time evolution generated by the fixed-point Hamiltonian of the phase with winding number ω = k+p. If |b k | > 1 then M (k) can be written as an imaginary time evolution exp(−arctanh(1/b k )H p+k ) followed by a unitary SPT entangler 9 W p+k , where W a = exp i π 2 H a . (It is straightforward to then show that b k ∈ R in Eq. (5) is equivalent to the wavefunction being real, as required for the BDI class.) While the imaginary time evolution cannot change the winding number, these SPT entanglers permute the fixed-point Hamiltonian as follows: Using this identity, we can move all SPT entanglers so that they act on the initial state, shifting it to the fixed-point ground state with the same winding number as the model under consideration. Result 2 can therefore be paraphrased as follows: Generic states on the MPS skeleton of translationinvariant BDI models with winding number ω are equivalent to sequences of imaginary time evolutions with fixed point Hamiltonians H k applied to the fixed-point ground state of H ω . We specify generic states to exclude cases with |b k | = 1, and note that unlike Result 2, these imaginary time evolutions are not necessarily applied in order of increasing range k due to their transformation under the SPT entanglers. Result 2 implies Result 1 through a continuity argument, taking into account cases with |b k | = 1, given in Section V D.
An alternative understanding of our construction is as follows. In terms of the polynomial g(z), for each step k in Algorithm 1 we can define with g d (z) = g(z). The coefficients of g k−1 (z) are the entries of s after the k-th step. The algorithm decreases the degree step by step until we have g 0 (z) ∝ z 0 . If we consider a sequence of models by f k (z) ∝ z p g k (z) 2 , we derive Result 2 by showing that applying M (k) to the ground state of f k−1 (z) gives us the ground state of f k (z).
The ground state of f 0 (z) ∝ z p is the fixed-point state |ψ p . We note that the finite-depth circuit representation of the ground state in Eq. (6) holds for both infinite as well as finite periodic chains (for fermionic and spin chain representations). However, when rewriting this circuit as a translation-invariant MPS in Section IV C, the treatment will be most natural for the case of an infinitely-long chain.

D. Consequences
Beyond constructing the MPS ground state on the MPS skeleton, the above results have some interesting consequences.
Firstly, we can construct a path of MPS-solvable models between any two gapped phases, labelled by wind- 9 This follows from the identity arctanh(b k ) = arctanh(1/b k ) − isign(b k )π/2. Note that W † a = WaP where P is the fermion parity. Hence for our purposes we can always use Wa.
ing numbers ω 1 and ω 2 , as follows. Let us first consider the case ω 1 − ω 2 = 2k for some k ∈ N. Then define the path: f (z) = z ω2 (z + a) 2k , where a ∈ R. For a = 0 we have f (z) = z ω1 , while for a → ∞ we have f (z) = z ω2 . At a = 1 we have a phase transition; at that point, using the results of Section V C, the ground state is an MPS with the fixed-point ground state of f (z) = z ω2+k . If ω 1 − ω 2 = 2k + 1 then first take f (z) = z ω1−2 (z+a) 2 for 0 ≤ a ≤ 1. At the point a = 1 we are connected to the following path (at the point A = 1): f (z) = z ω1−1 (Az + 2 + A/z) = z ω1−1 h(z). Then taking A → 0 we have a path (with constant ground state independent of h(z)) connecting to f (z) = z ω1−1 . We then can use the previous path to connect to z ω2 . This construction is simply an example, and we will encounter other paths in the next section.
Secondly, one can show that any model in the BDI class arises as a limit of a sequence of models with an MPS ground state. Indeed, this follows from being able to obtain a generic polynomial as a limit of polynomials which are squares. We showcase this phenomenon explicitly in Section III C by constructing a path of MPS-solvable models which converge toward the quantum Ising chain in a transverse field. We demonstrate how this can be used to extract the scaling dimension ∆ = 1/8 associated to the Ising universality class [62]. The general claim is proved in the concurrent work Ref. [63].
Thirdly, if we are in the case where we have a unitary circuit representation for the MPS, then we can use this representation to derive a formula for the (string) order parameter. For p = 0 and |b k | < 1, such a unitary representation exists and we obtain: This result is derived in Section V A 2 and is noteworthy since it does not rely on Wick's theorem or the Toeplitz determinant theory that appears in standard approaches (see, for example, [38,64,65]; Ref. [56] gives results in the notation of this paper for the general BDI class). In fact, using Toeplitz determinant theory on the MPS skeleton leads to a number of interesting exact results-this is explored in the concurrent work [63]. Finally, in Section V B, we apply our results to particlenumber-conserving models protected by a sublattice symmetry (class AIII), containing deformations of the Su-Schrieffer-Heeger (SSH) chain [66].

III. EXAMPLES
Here we will discuss three examples for which our results can be applied. As a first example we take the simplest case for Result 1. This leads us to a model introduced by M. Wolf et al. [14].
The second example introduces a two-parameter model, and we find several MPS-solvable paths that make up the MPS skeleton. Certain special cases appear where our results do not strictly apply, however, we can still find the ground state wave functions (these special cases are analyzed in Sections V C and V D).
In a third example we discuss how the quantum Ising chain can be approximated by a series of MPS-solvable parent Hamiltonians. This is illustrative of how any model within the BDI class can arise as a sequence of Hamiltonians with an exact MPS ground state.
A. Transition from ω = 0 to ω = 2 Based on Result 1, the simplest example of an MPSsolvable model in the BDI class that one might come up with is one for which f (z) is the square of a first order polynomial (i.e, d = 1, p = 0). Let us take: with the parameter λ ∈ R. This gives the polynomial that parameterizes a family of models within the space of the three-parameter Hamiltonian where we have written out the fixed-point models H α defined in Section II A. Note that this model does not have the Z 2 × Z 2 symmetry that conventionally protects the cluster SPT phase, but it has the anti-unitary symmetry ( n Z n ) K (where K is complex conjugation) which also protects the cluster model [57]. This parameterized family of Hamiltonians is the same path of Hamiltonians that is introduced by M. Wolf et al. in Ref. [14] as an example for a quantum phase transition within the MPS-framework, and was later used by A. Smith et al. in Ref. [44] to simulate a quantum phase transition on a noisy intermediate-scale quantum computer 10 .
A global, positive prefactor in front of the Hamiltonian in Eq. (11) does not change its ground state, so we can normalize the Hamiltonian such that t 0 + t 1 + t 2 = 1. This makes it possible to draw the phase diagram, as shown in Fig In the fermionic representation, these correspond to Kitaev chains with distinct winding numbers. In the dual spin chain formulations, these phases are the trivial paramagnet (ω = 0), the Ising magnet (ω = 1) and the symmetry-protected topological cluster phase (ω = 2). The solid red and blue lines show the parameterized paths along which we can find the MPS representation of the ground state, see Eqs. (10) and (12), the dashed gray lines show the lines where one parameter equals zero. See also Refs. [14,44].
through λ from −∞ to ∞ means traversing the red curve from left to right. For λ = 0 the corresponding Hamiltonian is H 0 , at λ = 1 2 the phase transition occurs and at λ = 1 the corresponding Hamiltonian is H 2 . Note that this is Kramers-Wannier dual to the disorder-line of the XY chain [38,57].
Within the phase diagram, there is actually another line that corresponds to MPS-solvable Hamiltonians, shown in blue. On this line the polynomial describing the Hamiltonian is withλ ≤ 1/4. Note that h(1/z) = h(z), has no zeros on the unit circle and has a positive constant term for allλ < 1/4-hence, we have that the ground state is the same as f (z) = z along this entire path (see the discussion following Result 1). Forλ = 1/4 we are at the multicritical point (λ = 1/2 on the red curve), while for λ > 1/4 we are on a critical line. On this critical line, we still have the form f (z) = zh(z) but now h(z) has zeros on the unit circle and the low energy physics is described . . . The repeating unit element of the circuit forms a tensor A j αβ , as defined in Eq. (14), which shows the equivalence of the circuit to a matrix product state with bond dimension χ = 2.
Turning back to the question of finding the MPS representation of the ground state of the model in Eq. (10), we can apply Result 2. For b 1 we simply find Note that for λ = 1 2 and λ → ±∞ we have that b 1 = ±1 and so Result 2 does not apply. These special points-which also happen to be the phase transitions-will be discussed as special cases in Section V C.
From the circuit construction of the state we can then obtain the usual MPS tensors (see Section IV C for the definition). The circuit construction and equivalence to an MPS is illustrated in Fig. 2.
We get the MPS tensor by interpreting the circuit gate as a four-legged tensor, where one ingoing leg acts on a spin |↓ , one outgoing one corresponds to the physical index, and the two legs connecting the ladder structure can be interpreted as the virtual legs-this is illustrated in Fig. 2. Therefore we have which in this case gives the MPS tensor We can compare our solution to the MPS given in Ref. [14] for this path. After rotating into the basis of Eq. (15) and inserting g = 2λ−1, the MPS from Ref. [14] is It can be checked that the matrix where s denotes the sign of 2λ − 1, relates the two MPS tensors as A j ∝ V −1 M j V . Therefore the two MPS representations are equivalent.
For all values of λ, we can use the results of Section V A to find a unitary circuit representation. This means that, using Eq. (8), we have the following expression for the order parameter. For the ω = 0 phase (λ < 1/2) we have: and for the ω = 2 phase (λ > 1/2): for further details see Section V A 2.
B. Transitions between ω = 0, ω = 2 and ω = 4 As a second example, we take a look at a path containing a phase transition between phases with the topological invariant ω = 0 and ω = 2, as well as a transition between phases with ω = 2 and ω = 4. We note that in the spin representation, H 0 , H 2 and H 4 are all in distinct interacting SPT phases protected by the Z 2 ×Z T 2 symmetry generated by P = n Z n and complex conjugation T = K [57].
Let us consider the model described by the polynomial with µ, ν ∈ R to ensure hermiticity of the Hamiltonian. By varying the parameters µ and ν we can explore the different phases of the model. The phase diagram is shown in Fig. 3. There, the differently colored regions correspond to the different phases, which are labelled by the topological invariant ω.
There are two choices for ν in Eq. (20) for which we can express f (z) as the square of a function g(z); if we choose ν = µ we find and if we choose ν = µ µ+1 we find For these particular choices of ν we can then apply Result 2 to find the MPS representation of the ground state.
The two options are also plotted as lines in the phase diagram in Fig. 3, where the blue line indicates the first case and the red line indicates the second case.
There is a third line, shown in the phase diagram in orange, that corresponds to a family of MPS-solvable Hamiltonians, but is not a square of a function g(z). If we choose ν = µ − 1, we find Here h(z) is a Laurent polynomial that satisfies h(z) = h(1/z), has a positive constant term and no roots on the unit circle for µ / ∈ {−1, 0, 1, 2}. Note that for 0 < µ < 1 we need to factor out a minus sign, so that the constant term of h(z) is positive. This means that the ground state along the orange line in Fig. 3 is that of H 2 for µ < 0 and µ > 1, and that of −H 2 for 0 < µ < 1.
Let us now return to the two cases of f (z) = g(z) 2 in Eqs. (21) and (22). As the first case (the blue line in Fig. 3) is essentially the previous example up to replacing p = 0 by p = 2, we will focus on the second case herethis is the red line in Fig. 3.
In order to apply Result 2, we first need to expand the function g(z) in Eq. (22) and calculate the set of b k . Expanded, g(z) becomes Using Algorithm 1 we can calculate b 2 = µ+1 , which fully specify the gates M with |ψ 0 being the ground state of H 0 . This result holds as long as |b k | = 1 holds for k = 1, 2 and for µ ∈ R. One can check that the cases where this does not hold are µ ∈ −1, − 1 2 , 1 where |b 1 | = 1, and The values of µ where |b 1 | = 1 happen to be the phase transitions (see Fig. 3). At these points the usual procedure fails but it turns out Parameter µ that we can find the gates that construct the ground state in all cases. We make some general points about these exceptional cases in Sections V C and V D. In particular, the ground state of the model with µ ∈ −1, − 1 2 , 1 is constructed in Section V C and the ground state of the model with µ = 1 2 1 ± √ 5 is constructed in Section V D. While the case µ = 1 2 1 ± √ 5 is actually already included in the discussion of the model in Eq. (23), Section V D gives a more general discussion of cases where |b k | = 1.

C. A path of MPS for the quantum Ising chain
We now show how our results can also tell us something about general models in this class. We will take as an instructive example of such a model that does not have an exact MPS ground state (for J = 0). The fermionic chain with this Hamiltonian interpolates between the trivial and the Kitaev chain (with critical point at |J| = 1); while the corresponding spin chain is the transverse field Ising model. Even though this model is not MPS-solvable, we will construct a sequence of MPSsolvable models which converges towards it. Note that the idea used here to find the path can be generalized to give a path of MPS-solvable models converging towards any Hamiltonian in the BDI class-the proof can be found in Ref. [63]. The related Laurent polynomial is f (z) = 1 + Jz, from which we can read off the topological invariant ω = 0 for |J| < 1 and ω = 1 for |J| > 1. We will focus on the ω = 0 phase 11 , as well as the critical point. Of course, 1 + Jz is not a square, and thus does not have an exact MPS ground state. However, we can write f (z) = g(z) 2 with g(z) = √ 1 + Jz. We can then use the series expansion to expand g(z). This expansion is valid if it converges to f (z) on the unit circle-indeed it is on the unit circle where we connect f (z) to our Hamiltonian (as discussed in Section II, the absolute value of f (e ik ) gives the energy spectrum and, moreover, its phase encodes the singleparticle modes [57]). Hence, if |J| ≤ 1, we can define (28) This converges to the quantum Ising chain: corresponds to a Hamiltonian on the MPS skeleton, and has an exact MPS ground state with bond dimension χ = 2 m . Note that this path can be used even for |J| = 1: for all m, all roots of f m (z) lie strictly outside the unit disk 12 ; hence, f m (z) gives a path of gapped Hamiltonians that approximate a critical Hamiltonian.
More explicitly, for any m ∈ N + , the perturbed Ising chain (with |J| ≤ 1) which has an exact MPS ground state with χ = 2 m is given by The perturbation δH is obtained by calculating the coefficients of f m (z) = g m (z) 2 and using binomial identities to simplify the double sum, resulting in where H α = − 1 2 n X n Z n+1 · · · Z n+α−1 X n+α are (fixedpoint) generalized cluster models.
One can prove that the MPS path in Eq. (28) is the optimal path of MPS approximations in the space of polynomials g m (z) which do not have roots inside the unit disk; see Appendix B 1 for a proof. Moreover, the same derivation also tells us that the energy density E m = ϕ m | (Z n − JX n X n+1 ) |ϕ m /2, where |ϕ m is the ground state of f m (z), is given by: Indeed, this converges to E ∞ := lim m→∞ E m which is equal to − 1 2π 2π 0 √ 1 + J 2 + 2J cos kdk, the ground state energy density of the quantum Ising chain. From Eq. (31) we also learn that the deviation for a given truncation m is In particular, for |J| < 1, the energy deviation decreases exponentially in m. Do note that since m = log 2 χ, this is only a polynomial decay in χ.
We thus obtain a path of MPS with ever-increasing bond dimension that converges to the ground state of the quantum Ising chain. Moreover, the above mechanism (using series expansions) can be applied to any generic model in the BDI class; this is worked out in more detail in the concurrent work [63] where it is used to analytically derive results about generic models. Note that sequences of free-fermion MPS converging to the ground state of the XY model are investigated in Ref. [47]. The approach there is valid for a particular region of the phase diagram that includes the quantum Ising model; however, in contrast to our path, performing the truncation requires numerical calculations.
Given this path of MPS ground states approximating the Ising model, it is interesting to see what we can derive about the critical model this way. Let us recall that the Ising CFT has two non-identity local scaling operators and σ (corresponding to the energy term and order parameter, respectively) as well as two nonlocal ones, µ and ψ (corresponding to the disorder operator and fermion, respectively) [62]. On the lattice, · · · Z n−2 Z n−1 Z n ∼ µ(x). We will now show how to extract the scaling dimension ∆ µ of µ using the above path of MPS. (Note that due to Kramers-Wannier duality, this also immediately gives us ∆ σ = ∆ µ .) We will use the path in Eq. (29) (with J = 1) where δH in Eq. (30) detunes it into the trivial paramagnetic phase for any finite m ∈ N. This detuning gives a finite energy gap ε π (at k = π) and long-range order to the disorder parameter, lim N →∞ Z 1 Z 2 · · · Z N −1 Z N =: µ 2 = 0. We will obtain both quantities. From their relative scaling µ ∼ ε ∆µ/∆ π and noting that it is a simple argument 13 13 For the quantum Ising chain, f (z) = 1 + Jz, hence the gap scales to derive that ∆ = 1, we thus extract ∆ µ .
First, at J = 1, the gap is given by Second, in Section V A 2, we derive a formula for the order parameter in the ω = 0 phase that is applicable to our path, from which we obtain µ = Note that {b k } k=1,··· ,m (which implicitly depend on m) can be efficiently obtained from the coefficients of g m (z) by order m multiplications and additions 14 . We plot the result in Fig. 4 where we find µ ∼ 1/m 1/8 . More precisely, by fitting the exponent (black dashed line), we find the critical exponent associated to the disorder operator of the Ising CFT: of approximations to this determinant, based on expanding the square-root as in Eq. (27). We point out that our treatment does not require this theory. An advantage of the above is that it gives an analytic expression for a path of MPS-solvable parent Hamiltonians which is optimal in some respect, as explained above. However, it is not optimal in the space of all MPSsolvable Hamiltonians: by allowing roots of g m (z) inside the unit disk, the variational energy can be decreased. As a χ = 2 example, consider f var (z) = 1 where |z 1 | < 1 and |Z 1 | > 1, which has winding number ω = 0. For any such choice of roots, the resulting MPS will have χ = 2. One can optimize these roots to minimize the variational energy with respect to the quantum Ising chain Hamiltonian given above. The variational energy is given by the negative real root of greatest absolute value of the equation: see Appendix B 2 for details. For J = 0, Eq. (35) gives the exact ground state energy density, and the deviation from the exact result increases with 0 ≤ |J| ≤ 1. For the critical quantum Ising chain (i.e., |J| = 1), the best variational energy for χ = 2 (using f var (z)) is which can be compared to the exact result − 2 π ≈ −0.63662. We thus see that the χ = 2 free-fermion MPS can reproduce the correct energy density within 0.3%. For |J| = 1, this is almost an order of magnitude better than the variational energy given by Eq. (31) for m = 1, namely −5/8 = −0.625, which is within 2% of the true energy density.

IV. ANALYSIS
To derive the results stated in Section II, we first show that the ground state is frustration-free if the associated polynomial is of the form f (z) = z p g(z) 2 . In Section IV B, using Witten's conjugation method [70,71], we then derive an explicit (non-unitary) circuit mapping its ground state to a fixed point wavefunction. This circuit can then be explicitly rewritten as a matrix product state with the claimed bond dimension, this is derived in Section IV C.
The results in Section II were stated for both the fermionic chain and the Jordan-Wigner dual spin chain simultaneously; however, certain sections below are more straightforwardly presented with one or the other picture. In particular, Sections IV A-IV B make use of the fermionic notation for ease of presentation. In Section IV C we explain how to construct the MPS tensor. While fermionic MPS are well understood [9,72], the spin chain representation is more conventional. Despite working with the Hilbert space of the spin chain, it will still be useful to present certain formulas using fermionic operators throughout this work. The underlying Jordan-Wigner transformation is given by: While care usually has to be taken about the precise meaning of the product over sites j < n, dependent on boundary conditions, in Section IV C we will implicitly be working in case of an infinitely-long chain, such that these subtleties do not arise. The Jordan-Wigner dual expressions for h n,α = iγ n γ n+α were already given in Eq. (3), which includes the identity Z n = iγ n γ n . We also have: The ground state of H 0 is denoted by |ψ 0 . This corresponds to the completely filled fermionic state c † n c n |ψ 0 = |ψ 0 ; while for spins the corresponding state is |↓ · · · ↓ . Using these identities, all formulas below can be transformed from fermions to spins and vice-versa.

A. The Hamiltonians are frustration-free
A frustration-free model is one where the Hamiltonian can be written as a sum of terms such that each term is individually minimized in the ground state [71,73].
Here we derive the frustration-free property of the above systems. This will also form the starting point of the wavefunction construction in Section IV B.
If f (z) = z p g(z) 2 with g(z) = s α z α , then Eq. (1) can be written as which is confirmed by expanding it out. Note that each term indexed by n in Eq. (39) has eigenvalues ±| s| 2 . We now show that the ground state minimizes each term in Eq. (39), i.e., that the energy density is e 0 = − 1 2 | s| 2 . This is the defining property of a frustration-free model.
For any f (z), the ground state energy density can be expressed as where the contour integral is along the unit circle. If f (z) = z p g(z) 2 , this simplifies to As explained above, this shows that the model is frustration-free.
For what follows, it will prove to be useful to define Γ n = α s α γ n+α+p − iγ n−α and then to rewrite Eq. (39) as By expanding Eq. (42) and observing that terms of the form iγ n γ m and iγ nγm do not survive (either by explicit computation or by noting the complex-conjugation symmetry of the model), one verifies that it equals Eq. (39). Similarly, one sees that the frustration-free property of Eq. (39) is equivalent to the ground state |ψ of Eq. (42) satisfying Γ † n Γ n |ψ = 0 for all n. For the fermionic case under consideration 15 , this is in turn equivalent 16 to Γ n |ψ = 0.

B. Constructing the circuit (Result 2)
The MPS parent Hamiltonian construction leads to a frustration-free Hamiltonian where a given MPS is the ground state. The converse typically holds, although not all frustration-free models have MPS ground states 17 : an example is given in Ref. [77]. Here we give a direct proof that this is the case for translation invariant BDI models with f (z) = z p g(z) 2 by explicitly deriving the quantum circuit that constructs the ground state.
Firstly, note that it is sufficient to prove this for p = 0. Indeed, one can shift f (z) → z q f (z) by shifting all h n,α → h n,α+q (see Eqs. (1) and (2)). Hence, we see that starting from the p = 0 result one recovers p = 0 by shifting h n,α → h n,α+p in all formulas.
Secondly, we assume a positive overall sign of f (z). Different global signs of f (z) are related by the unitary transformation S, given by This global sign generically does not affect the analysis so we can account for it by applying S to the final state.
(The only exceptions to this are cases with zeros on the unit circle, discussed in Section V C.) Under conjugation by S we invert the gates M (k) 15 Note that the Jordan-Wigner transformation of Eq. (42) for periodic boundary conditions will give nonlocal 'boundary' terms involving phase factors. However, using the translation-invariance of the state, a frustration-free local spin Hamiltonian is obtained by simply dropping the nonlocal phase factor. 16 E.g., writing the singular value decomposition Γn = U SV , we see that Γ † n Γn = V † S 2 V . Hence, they have the same zero eigenvalues/eigenvectors. 17 References [74][75][76]  To reach Result 2, we derive the following stronger result:

Result 3 (Relating circuits and polynomials.)
If |ψ i is the initial ground state associated to some polynomial f i (z) = g i (z) 2 , then for any k ∈ Z and B k ∈ R with |B k | = 1, the transformed state is the ground state for (46)

Analysis: Construction of MPS.
Using Result 3, one can start from the trivial case g(z) = 1 (with trivial ground state |ψ 0 ) and successively apply layers of gates to obtain the ground state of the desired g(z). Let us first take an example. Starting with g i (z) = 1 and applying a layer generated by the Kitaev or Ising chain, M (1) , we obtain g f (z) = 1 + B 1 z. With a second layer, M (2) , applied to g i (z) = 1 + B 1 z, we obtain the ground state of g f (z) = 1 + B 1 z + B 2 z 2 (1 + B 1 /z) = B 2 z 2 + B 1 (1 + B 2 )z + 1. If we set B 1 = 4 3 and B 2 = 2, we recover the polynomial g(z) = 2z 2 + 4z + 1 that we discussed in Section II.
More generally, let us show that the recursion in Algorithm 1 leads to the desired ground state. Define b k for k = 1, . . . , d using this recursion. Recall that we start from the Hamiltonian corresponding to f (z) = g d (z) 2 and then the recursion defining b k amounts to Eq. (7), which is, for each k, given by: This is equivalent to: As long as |b k | = 1 we can thus invert the recursion; the term (1 − b 2 k ) = 0 is an unimportant constant. Hence, setting B k = b k , we can work up from the ground state of f 0 (z) = 1 to f (z) = g d (z) 2 using Result 3.
We now have Result 2 for all values of p, i.e., the ground state of f (z) = z p g d (z) 2 is given as a product as long as |b k | = 1 for any k. We now note that the fixed-point wave function |ψ p can itself be written down as a circuit, the precise form depending on whether p is odd or even: |ψ 2q = W q |ψ 0 |ψ 2q+1 = W q+1 |ψ 1 = W q+1 n P n |ψ 0 (50) where W q are the SPT entanglers introduced above and are commuting projectors. Due to the two cases in Eq. (50), it is helpful to define q ∈ Z and r ∈ {0, 1} such that p = 2q + r. By writing the SPT entangler W q as a product of unitary and mutually commuting operators we have the following representation for the fixed-point wave function |ψ p : where U n = 1 Note that Result 3 is the fundamental statement, and can be used to transform between ground states of any two Hamiltonians in our class that are related by Eq. (45) (including relating Hamiltonians that are cases with |b k | = 1 and where we do not have a construction using Result 2); for a particular choice of transformations we derive Result 2 from Result 3.

Analysis: Relating circuits and polynomials
To derive Result 3, we start with a ground state |ψ i corresponding to a polynomial f i (z) = g i (z) 2 . Writing g i (z) = α s α z α , we know that |ψ i is annihilated by Γ n in Eq. (42) (with p = 0). Hence, M (k) |ψ i is annihilated byΓ n := M (k) Γ n M (k) −1 , implying that it is the ground state of H = 1 2 nΓ † nΓn . All that remains is to calcu-lateΓ n to confirm that it corresponds to the polynomial defined in Eq. (45).
Since M (k) n = 1 − A k iγ n γ n+k (remember that we have set p = 0), its inverse is well-defined because A 2 k = 1 (equivalently, |B k | = 1, see Eq. (46)). Hence, up to an irrelevant global factor, we have where in the last step we divided by (1 + A 2 k ). Similarly, Taken together, we see that Γ n = α s α (γ n+α − iγ n−α ) gets mapped tõ Hence, the transformed state is the ground state corre- This completes the proof. Note that if instead B k = ±1, the gate in Eq. (46) becomes a projector onto the ground state of ±H k . Hence in that case, |ψ f is zero if |ψ i is the ground state of f (z) = ∓z k , otherwise |ψ f in Eq. (44) is the ground state of f (z) = ±z k . See Section V D for further discussion of these cases.
We note that this way of constructing frustration-free Hamiltonians-i.e., where Γ n is conjugated by an invertible operator M -is known as the Witten conjugation method [70,71]. Usually, one simply chooses M and considers the resulting frustration-free Hamiltonian. What is special to our case is that we have an explicit formula for a set of M that are quadratic fermionic gates and can be used to generate any frustration-free model in the BDI class.

C. MPS representation
Given the explicit circuit construction of the ground state from the previous section, namely Eqs. (49) and (53), we will now show that this corresponds exactly to a finite bond dimension translation-invariant MPS. This means that the ground state can be written as the contraction of a translation-invariant tensor network, where the virtual indices between unit cells have finite dimension. For a spin chain, a translation invariant MPS, with MPS tensor A, is a state of the form For fixed j, A j is a χ × χ matrix, where χ is the bond dimension [9]. The fermionic case is similar, for details see Ref. [72]. For the purposes of defining the MPS tensor, A, it is simplest to contract an index of the circuit with the state |↓ n , and so we will work in the spin chain picture. The algebraic steps involved apply (by definition) in the same way to spin or fermionic operators; and for ease of presentation we will continue to use fermionic notation for the gates. These fermionic operators can be interpreted as short-hand for the spin operators as given in Eqs. (3) and (38).
To better understand the product of operators in Eq. (49), we introduce a graphical notation, defined in Fig. 5. All of the operators appearing in this product have the same form, namely, 1 + ch n,α = 1 + icγ n γ n+α for some c ∈ C. We represent the sites by black "wires", similar to quantum circuit diagrams. The operators are then denoted by a colored "gate" with ends labelled with either circles denoting γ n operators, or squares forγ n operators. Using this notation, an example of a product of the form in Eq. (49) is shown in Fig. 6.
The equivalence to MPS is established by grouping the gates in this product into a repeating unit element, illustrated by the gray box in Fig. 6. This repeating element has one wire that corresponds to a site and several wires that connect to different unit cells, corresponding to the virtual indices of an MPS tensor. The bond dimension of the corresponding MPS is χ = 2 N , where N is the number of wires connecting the unit cells, see Figs. 7(b-c).
While we have put the circuit in MPS form, this provides a very loose bound for the bond dimension χ = 2 N , where N is quadratic in d and linear in p. However, by using the commutation properties of the gates we can provide a much tighter bound on the bond dimensionin fact, we conjecture that this bound gives the optimal 6. Example of the mapping from circuit to MPS for d = 2, p = 3 using the graphical notation in Fig. 5.
The different colored gates correspond to the different layers in the state construction. The gray box indicates the repeating circuit element equivalent to an MPS tensor.
bond dimension, and this can be proved in certain cases. These gates are able to commute past each other except when the same symbol in our graphical notation is acting on the same spin, i.e., we cannot bring a circle past a circle or a square past a square, see Fig. 5(b). Furthermore, we are able to bring the unitary gate U n (defined in Eq. (53)) past those appearing in M (k) , which results in the algebraic relation shown in Fig. 5(c). Let us consider here p = 2q + r > 0, (see Appendix C for the more general expression), then whereM (k) n+q+r = 1 − ia k γ n+q+r γ n+k+p . To see this, This is shown schematically in Fig. 5. Using Eq. (38), we have the corresponding spin operator. These algebraic relations allow us to drastically reduce the bond dimension, as shown for the particular example of d = 2, p = 3 in Fig. 7(a). This more compact form follows from the non-trivial application of the commutation relations for the gates and the algebraic identities using the U n unitary gates and is explained in detail in Appendix C. In general the bond dimension is given by χ, where We show this in Appendix C using the methods introduced here. The bond dimensions for different values of d and p are shown in Tab. I. The bond dimension in Eq. (61) is an upper bound for the bond dimension required for an exact MPS representation of the ground state. That is, this bond dimension is sufficient for an exact representation. In the case that we have a gapped model on the MPS skeleton with no mutually inverse zeros 18 , we believe that this bond dimension is also necessary; i.e., Eq. (61) gives the optimal bond dimension. (For gapless points in the MPS skeleton, we show in Section V C that the ground state can be found by considering a related gapped model.) In Ref. [63], it is explicitly shown in the spin chain representation that when p = −d and p is even, we have a lower bound on χ that coincides with Eq. (61), and hence this proves the optimality in this case. This is proved by analyzing ground state correlation functions in these models 19 . To test the formula for the optimal bond dimension more generally, we compare the analytical upper bound with the bond dimension obtained from finding the ground state numerically using the density matrix renormalization group (DMRG) [2,6] with explicitly conserved global Z 2 symmetry. We find that the numerically and analytically obtained bond dimension perfectly coincide for all cases tested, suggesting that this bond dimension is indeed optimal beyond the cases where we 18 As seen in Section III A, there are cases with mutually inverse roots where this upper bound does not give the optimal bond dimension. These are models defined by f (z) =f (z)h(z) where h(z) = h(1/z), h(z) has a positive constant term and has no zeros on the unit circle. In such cases, we conjecture that the bond dimension will be determined byf (z) according to Eq. (61). 19 Denote zeros of g(z) inside the unit circle by z j and zeros outside the unit circle by Z k and let S = {z 1 , . . . , zn z , Z −1 1 , . . . Z −1 n Z }. The result in [63] also assumes that given any subset of S, if we take the product of zeros in that subset then the absolute value of that product is different to any other subset except for subsets containing any conjugate roots-this condition holds generically for gapped models with f (z) = z p g(z) 2 .  have an analytic proof. It would be of interest to see if the methods of Ref. [78] could be used to find a lower bound that coincides with our upper bound.

A. Unitary version
The circuit construction described in the previous sections is generally made from non-unitary gates M (k) and projectors P . However, in certain cases we can instead use a circuit made entirely of unitary gates. Note that while any MPS can be written as a sequential unitary circuit [79], the depth of the repeating unit element would scale with bond dimension χ, whereas our construction scales as log χ. Such a unitary circuit representation may be useful for processing quantum states on quantum computing platforms. We first explain how this works for p = 0, 1, and a k ∈ R for all k (equivalently |b k | < 1 for all k, recall that a k is defined in Eq. (6)). We then show how this can be extended to all p ∈ Z and a d . We show in Section V A 2 that the condition that |b k | < 1 for k < d means that the zeros of g(z) are either all inside or all outside the unit circle.

Circuit construction
To put our circuits into unitary form we need to use a substitution that is demonstrated in Fig. 8(a). Schematically, we are able to turn a square symbol acting directly on the initial state into a circle (see Fig. 5 for the definition of these symbols). Explicitly, we have the relation where X(n) is any operator that is not supported on site n. This follows from iγ n γ n |ψ 0 = Z n |↓ · · · ↓ = − |↓ · · · ↓ . On the right-hand side of Eq. (62) we have a gate proportional tõ which is unitary for a k ∈ R. In the case that p = 0 or p = 1, from Eq. (5) (and the discussion following it), our circuit is of the form for some A k ∈ R. Indeed, if p = 0, then A k = a k whereas if p = 1, then A 1 = 1 (this is the gate P n in Eq. (51) that projects the initial state into the ground state of H 1 ) and A k>1 = a k−1 . We can then make the substitution in Eq. (62) for each of the (1 − iA 1γn γ n+1 ). The resulting V (1) n will commute past the other gates and so we can bring it down, as illustrated in Figs. 8(c-d). We then have the gate (1 − iA 2γn γ n+2 ) acting on the state X(n) |ψ 0 , where X(n) is some operator not supported on site n. Repeatedly using the substitution and commuting gates, we end up with a version of the circuit consisting of only unitary gates, as shown in Fig. 8(f).
By applying extra layers of unitary gates we can extend the set of states that we can construct with unitary circuits to include all p ∈ Z and no restriction on b d (i.e. a d can be real or on the unit circle). This is achieved using the SPT entanglers W α , recall that these are given by: Conjugating by W α corresponds to the mapping t n → t 2α−n in the Hamiltonian in Eq. (1), or, equivalently, in the polynomial f (z) in Eq. (2). By using W d+p we can remove the constraint on a d 20 . Furthermore, combining 20 Note that applying this transformation takes g(z) → z d g(1/z).
two of these entanglers results in an even shift, that is W β W α : t n → t n+2(β−α) . Starting from p = 0, 1 this allows us to transform to any value of p. Since all of the gates appearing in the product in Eq. (65) commute, we can similarly collect gates into a repeating unit cell.
It is important to note that this alternative construction of states using W α does not generally correspond to the optimal bond dimension found in previous sections. However, the unitary circuit representation may be of practical use, for instance in processing these states on a quantum computer [44].

A formula for the order parameter
As a result of the unitary circuit representation given above, we can derive a formula for the order parameter in those models. In this section we will work in the spin chain picture. Let us consider g(z) = d k=0 s k z k , such that the corresponding |b k | < 1. Then, we have: where the left-hand side is the string order parameter for ω = 0.
To prove this, consider the ground-state correlation function for a half-infinite string ∞ j=m Z j , which is This takes b d → 1/b d and b k → b k for k < d. Hence, we still require |b k | < 1 for k < d.
given by: Note that nŨ n is such thatŨ n acts beforeŨ n+1 on the string-see also Fig. 8(f). In terms of spins, we have: Then, for (m − k) ≤ n < m we havẽ where For all other values of n,Ṽ (k) n commutes with ∞ j=m Z j . Moreover,Ṽ (k ) n commutes with O n,k for k > k. One can then see that, since ↓| n Y n |↓ n = 0, the second term in Eq. (69) does not contribute. Then by replacing: implying the result. Now, within the BDI class, lim N →∞ | Z 1 . . . Z N | = 0 is equivalent to being in the gapped phase with ω = 0 [56] (and in particular, lim N →∞ | Z 1 . . . Z N | → 0 implies that the gap closes). This means that, since p = 0, Eq. (66) tells us that if |b k | < 1 then all zeros of g(z) are outside the unit circle. In fact these conditions are equivalent: if all zeros of g(z) are outside the unit circle, then all |b k | < 1. To see this, consider g 0 (z) = 1 + z d for < 1. This has all zeros outside the unit circle, and b d = , b k<d = 0. Let g(z) have all zeros outside the unit circle, we can tune the zeros of g 0 (z) to the zeros of g(z) along paths outside the unit circle, and this corresponds to a path of gapped Hamiltonians. Moreover, the b k vary continuously along this path, and at no point along the path can we have |b k | → 1 as this would contradict the fact that the path is gapped, hence, g(z) also has all |b k | < 1. As explained above, the case |b d | > 1 and |b k | < 1 for k < d can be analyzed by applying the SPT entangler W d . This transformation maps zeros of g(z) to inverse zeros of g(z), and so must correspond to the case that g(z) has all zeros inside the unit circle.
Recall that g(z) = d k=0 s k z k , and let us fix s 0 = 1 for convenience. As an illustration, we can evaluate Eq. (66) for d = 1: where |s 1 | < 1 and for d = 2: where |s 2 | < 1 and |s 1 | < |1 + s 2 |. If we denote the zeros of g(z), which lie outside the unit circle, by Z 1 , . . . , Z d , then we also have: This is a special case 21 of a more general formula for the order parameter in terms of zeros of g(z), given in Ref. [63]. This means that: Note that for d = 1, s 1 = 1/Z 1 , so using Eq. (73) we see immediately that this equality is satisfied. For d = 2 one can show directly that Eq. (74) and (75) are equal. Eq. (66) applies for the case p = 0. By applying SPT entanglers, it is an immediate consequence that if we allow general p ∈ Z with the condition |b k | < 1 then: where O p is the (string) order parameter for the phase with winding number p. In the spin chain representation O p is local (non-local) for p odd (even)-general definitions are given, for example, in Ref. [56]. We Thus far we have focused on the BDI class which contains (superconducting) pairing terms. The Kitaev chain [54] is the generating SPT of this class, in the sense that all topological phases are obtained by considering stacks.
In this section, we discuss the AIII class which preserves particle number and has the SSH chain [66] as its generator. We will show that it can be embedded into the BDI class, thereby offering a reinterpretation of some of our results. Similar to Eq. (1), we can define the following fermionic hopping model: As before, we define the range of this Hamiltonian to be the largest |α| of all non-zero τ α . For τ α ∈ C, this is a general translation-invariant model in the AIII class. Let us first discuss the special case τ α ∈ R. As explained in Appendix D 1, Eq. (78) can be rewritten as a translationinvariant Majorana chain where the range has been doubled; more precisely, it has the form of the BDI Hamiltonian in Eq. (1) with t 2α = τ α and t 2α−1 = 0. Hence for τ α ∈ R our results above apply directly to these models. In particular, there exists an MPS with bond dimension log 2 χ = range(H); and if all |b k | = 1 then we have a construction of this MPS.
To state the analogous construction, let us define the fixed point Hamiltonians For instance, the case k = 1 is the SSH chain [66]. We define the associated polynomial f (z) = α τ α z α = z p g(z) 2 . (As before, the quasiparticle dispersion relation is ε k = |f (e ik )|. Recent work has established the relation between f (z) and topological edge modes [80].) If we calculate the b k (and corresponding β k ) exactly as above, then we have that the ground state is given bŷ As before, theseM (k) correspond to imaginary time evolutions with fixed point Hamiltonians, followed by SPT entanglers for |b k | > 1.
In Appendix D 1 we show that the general Hamiltonian given in Eq. (78) with complex hopping τ α ∈ C is equivalent to a Majorana chain with a two-site unit cell, for which we have not given an explicit construction of the ground state in the present work. In the concurrent work [63], the existence of an exact MPS ground state is proved in the case that f (z) = α τ α z α = z p g(z) 2 , but no construction or upper bound on the bond dimension is given. We conjecture that the same bond dimension formula holds in this case (i.e., that log 2 χ = range(H)).

C. Multicritical points
Since |f (e ik )| gives the single-particle energy spectrum, if the polynomial f (z) has zeros on the unit circle, at points z = e ikn , then the system is gapless. Let us consider gapless models on the MPS skeleton: the starting point is, as above, that f (z) = z p g(z) 2 h(z), where h(z) has no zeros on the unit circle, and so any zeros on the unit circle must have even multiplicity. In the phase diagram of translation-invariant BDI models, these are multicritical points. Rather than applying our algorithm immediately, we first simplify the problem by reducing to an equivalent gapped model. This is possible due to the even multiplicity of all zeros on the unit circle. In Appendix D 2 we show that the ground state of a gapless model given by f gl (z) = z p g gl (z) 2 h(z) is the same as the ground state of a closely related gapped model f g (z) = (−1) m0/2 z p+Nc/2 g g (z) 2 . Here g g (z) is the polynomial g gl (z) after dividing by (z − e ikn ) for all zeros on the unit circle, N c /2 ∈ Z is the number of zeros of g gl (z) on the unit circle counting multiplicity and m 0 is the multiplicity of the zero at z = 1 (if there is no zero there, m 0 = 0). This means that we can find the ground state at the multicritical point by applying our methods to this related gapped model.
To see how this works in practice, we can refer back to our earlier examples. In the first example, see Eq. (10), we have a gapless point at λ = 1/2. In that case f gl (z) = (z + 1) 2 and this implies f g (z) = z. Hence, this ground state is simply given by the Ising ferromagnet. Similarly for λ → ±∞ (normalizing f (z) appropriately) we have f gl (z) = (z − 1) 2 , leading to f g (z) = −z.
In the second example, see Eq. (22), we have gapless points when µ ∈ {−1/2, ±1}. Then the ground state can be obtained from Note that the case µ = −1 requires taking a well-behaved limit for the ratio f g (z)/|f g (z)|.
D. Cases where |b k | = 1 Above we impose the condition that |b k | = 1. This ensures that both the gates M (k) defined in Eq. (6) and the recursion in Eq. (7) are invertible, and moreover that Algorithm 1 is unambiguous. Here we discuss what happens when we relax this condition.
First, take a model defined by f (z) = g(z) 2 and recall the approach in Section IV B 1. By inverting the recursion in Eq. (7) and then fixing B k = b k in Result 3, we could transform a fixed-point ground state into the ground state of f (z). We could equally have set B k = −b k in Result 3, giving a transformation from the ground state of f (z) to a fixed-point ground state. These are equivalent since b k → −b k is the same as M (k) → (M (k) ) −1 .
In the case that |b k | = 1 for k > k 0 and |b k | = 1 for k 0 , the second point of view is helpful. We can then use Result 3 for each k > k 0 to write: where |ψ is the ground state of f (z) = g(z) 2 and |ψ is the ground state of f (z) = g k0 (z) 2 . Our method breaks down at the next step because, as explained in Section IV B 2, if |b k0 | = 1 then applying amounts to applying the projector P (k0) (note we have a +b k0 here because we set B k = −b k ). There is a special case where this can work-when the state |ψ that the projector acts on is an eigenstate of the projector.
In particular, our method is still able to construct the ground state if the projector P (k0) annihilates this state 22 , i.e., Due to translation symmetry, we can conclude that |ψ is the ground state of ±H k0 . In Appendix D 3 we show that this case applies if and only if g k0−1 (z) = 0 (equivalently: in the application of Algorithm 1 the vector s for the iteration step k = k 0 − 1 vanishes). Moreover, we show that given b k0 = ±1, |ψ is the ground state of b k0 H k0 , i.e., |ψ = S (1−b k 0 )/2 |ψ k0 . Thus, from Eq. (82), we can construct the ground state of our initial model by: In the above discussion, we took p = 0. As in Section IV B, the result for general p follows by applying SPT entanglers that shift h n,α → h n,α+p and |ψ k0 → |ψ k0+p . These results are relevant to our second example, defined in Eq. (22).
In particular, at the points µ = 1 2 (1 ± √ 5), our algorithm leads to b 2 = 1. Applying Eq. (7) gives g 1 (z) = 0. Then using the above, we have that the ground state is |ψ 2 -the cluster state. Alternatively, as shown in Eq. (23), we have mutually inverse zeros at these values of µ-this means that the ground state is given by the ground state of f (z) = z 2 . It is easy to generalize this observation to any g(z) of degree two which has s 0 = s 2 , implying that b 2 = 1.
There are nevertheless models with |b k | = 1 where our approach does not work. Following on from the immediately preceding example, consider g(z) with degree two and with b 2 = −1, i.e., s 0 = −s 2 . This cannot be simplified unless s 1 = 0-then g(z) = (z − 1)(z + 1) which is a multicritical point with the same ground state as f (z) = −z 2 . For other values of s 1 we have a gapped model where our approach fails. One can argue that for these gapped models, we can define a perturbed model with s 0 = −s 2 + , where both the ground state wave function and the corresponding Hamiltonian continuously depend on the parameter such that the limit → 0 is well-defined. However, the question remains how to take this limit. Note that despite not having an explicit MPS in the limit, we can argue that the upper bound on χ remains valid at the limiting point. Indeed, the optimal χ 2 is the number of non-zero eigenvalues of the reduced density matrix of a subsystem, and since the state is continuous, the number of non-zero eigenvalues cannot be greater at the limit point.
To see explicitly that there still exists an MPS representation with an appropriate bound on the bond dimension we can consider the positive eigenvalues of the correlation matrix [81,82]. For a subsystem of size N there will be N of these eigenvalues, {ν 1 , . . . ν N }, where any ν j = 1 is a trivial eigenvalue. The eigenvalues of the reduced density matrix can be derived from these {ν j }, and the number of non-zero eigenvalues of the reduced density matrix is 2 x where x is the number of non-trivial ν j . Now, consider the model with where |z 1 | < 1 and |Z 1 | > 1. In [63], the eigenvalues of the correlation matrix in this model are found for any subsystem size. For a subsystem of size N → ∞, there are two generically non-trivial eigenvalues given by: Now, the model in Eq. (86) has b 2 = 1 when Z 1 = z −1 1 . Then ν 2 1 = ν 2 2 = 1: all correlation matrix eigenvalues are trivial. This is consistent with what we proved above, k 0 + p = 0, and thus our system has the ground state |ψ 0 .
The model in Eq. (86) has b 2 = −1 when Z 1 = −z −1 1 . Then This corresponds to a bond dimension χ = 2 which is the same as the general case Z 1 = ±z −1 1 where our construction above applies. Hence, although we do not have a construction of the MPS in the case b 2 = −1, the limiting point is an MPS with a bond dimension that is upper bounded by the path approaching it, as we expected by the general continuity argument.

VI. OUTLOOK
We have introduced the idea of the MPS skeleton underlying the phase diagram of one-dimensional models. To illustrate this concept, we have given a simple characterization of this skeleton for translation-invariant models in the free-fermion BDI class, as well as a construction of the MPS ground state for every model on the skeleton up to a measure-zero set of exceptions. Hamiltonians on this measure-zero set are limits of cases where we have the MPS construction; it would be interesting to see explicitly how to construct the MPS ground state in these cases. It would also be of interest to find a unitary circuit representation that applies more generally than the subset of cases we discuss above.
A natural problem is to extend our work to models in the BDI class with a larger unit cell, as well as other free-fermion classes. As discussed in Section V B, translation-invariant models in class AIII are equivalent to BDI models with a two-site unit cell, so these problems are related. The characterization of the MPS skeleton for the translation-invariant AIII class is the same as in the translation-invariant BDI class, and we expect it would be relatively straightforward to adapt the construction of the MPS given in this work to this class. A more challenging generalization would be to consider Majorana chains in class D where complex hopping and pairing are allowed, or even to two-dimensional free-fermion systems, where one can analogously consider free-fermion circuits and projected entangled pair state (PEPS) skeletons. The results of Refs. [52,53] on Gaussian MPS apply for higher spatial dimensions, and make clear the following important property of MPS-solvable models: the correlation matrix in Fourier space has entries that are rational functions. It remains to be shown that this is a sufficient condition for MPS-solvability. Perhaps this can be shown in a constructive manner; as done in the present work for the one-dimensional BDI class. Moreover, it would be helpful to clarify the exact relationship between the construction of the ground state MPS given in this paper with the Gaussian MPS and PEPS appearing in Refs. [52,53]. This could indicate how to generalize our construction to higher dimensions.
Finding the MPS skeleton for class of models beyond free-fermions would be very interesting. In particular, as discussed in Section II, we found that any MPS state in the BDI class can be obtained by starting with a fixedpoint wavefunction and applying a finite number of imaginary time-evolutions generated by fixed-point Hamilto- nians. An open question is whether this characterization remains true for interacting MPS.
Let us point out that the unitary circuit representa-tion of parts of the MPS skeleton presents a powerful approach for processing quantum states on near-term quantum computers. While these states are specified by a number of parameters b k proportional to the range of the Hamiltonian, the classical processing and extraction of observables requires the contraction of matrices with bond dimension exponential in the number of parameters. The quantum circuit, however, has a unit element with its depth directly proportional to the number of parameters, leading to a possible exponential speed improvement for processing the states. This approach has been demonstrated in Ref. [44] for the example discussed in Section III A. Whether this favourable scaling can be extended to the entire MPS skeleton is an interesting and open question; the methods in Ref. [48] could potentially be useful here. While this advantage is perhaps artificial for the exactly solvable and one-dimensional systems considered in this paper, possible generalizations may be of practical value for studying non-trivial topological phases on quantum computers. Finally, these exactly-solvable MPS states could serve as useful initial states which could be dressed with nonintegrable perturbations. For instance, the states we construct could potentially be useful initializations for Gutzwiller-projected DMRG [50,51,[83][84][85].

FIG. 9. Contour integral for the Fourier transform. (a)
The contour is the unit circle, integrating over this contour defines the Fourier coefficients of (A7). (b) The deformed contour gives the same integral and is snagged at poles and branch cuts (indicated by blue wavy lines).
text with f (z) = σz p g(z) 2 h(z), hence, any translationinvariant BDI Hamiltonian outside of this class has at least one odd m. Now, the ground state correlation function − h n,α is the nth Fourier coefficient of (A7) [56]. Other correlation functions can be derived from this one using Wick's theorem. The function (A7) is analytic on the unit circle and we can compute the asymptotics of − h n,α by analytic continuation. For n > 0, by deforming the contour out to infinity we pick up the dominant contributions wherever the contour gets snagged: at poles and branch points. This is the Darboux principle [86] and is illustrated in Figure 9. The function (A7) has poles at Z k (z −1 j ) for m Z k (m zj ) even, and branch points at Z −1 k and Z k (z −1 no roots inside the unit disk (i.e., q max (z) ∈ P α,β ), then q max gives the global maximum of E[q; p] for q ∈ P α,β . Moreover, E[q max ; p] = 1 2 β k=0 a 2 k .
Note that if we take such a p(z) with no roots inside the unit circle, then it converges uniformly and so we must have that q max (z) has no roots inside the unit circle for sufficiently large β. Moreover, when we apply this result for p(z) 2 = 1 + Jz, this holds for all β.
Proof: For any q(z) ∈ P α,β , letα ∈ N be the largest integer such that q −α = 0, and then letβ ∈ Z be the largest integer such that qβ = 0 (note that 0 ≤α ≤ α and −α ≤β ≤ β). Then we can write q(z) = z −αq (z) withq(z) ∈ P 0,α+β . We will now show that E does not depend on a k if k >β −α. First note that ∂ a k p(z) = z k , thus Sinceq(z) is a polynomial with no roots inside the unit disk, we know that there exists an expansion 1 q(1/z) = ∞ r=0 c r z −r which converges on the unit circle. Hence, the largest power appearing in the integrand of Eq. (B4) is z −k−2α−1+(α+β) = zβ −α−k−1 (note that we use that a 0 = 0, which indeed follows from the assumption that the partial sum of p(z) has no root inside the unit disk). If k >β −α, we thus see that there is no term proportional to z −1 in any of the integrals in Eq. (B4)-then, by the residue theorem, the derivative ∂E/∂a k is always zero.
Due to this independence, we can without loss of generality set a k = 0 for k >β −α, i.e., we truncate p(z). If we chooseα = 0, it is then possible to setq(z) equal to this truncated p(z), denoted pβ(z). To see that this is indeed the optimal choice, writeq(e ik ) = ρ k e iφ k and pβ(e ik ) = r k e iθ k (note that due to the real coefficients, we have thatq(1/z) = ρ k e −iφ k and similarly for pβ(1/z)), then This is clearly maximal if and only if φ k = θ k , which is achieved by settingq(z) = pβ(z). Finally, the value is then given by Although we have a local maximum of E for every allowed choice ofβ, Eq. (B6) is clearly globally maximized if we chooseβ as large as possible, i.e.,β = β.

Variational energy
We now explain how to derive Eq. (35). If we take the ground state of f var (z) = 1 z 2 (z − z 1 ) 2 (z − Z 1 ) 2 where |z 1 | < 1 and |Z 1 | > 1, then the energy density for the Ising Hamiltonian is given by: As explained in Appendix A, these expectation values are Fourier coefficients of (A7). These are calculated in [63] for models on the MPS skeleton. Using the result for our case, we have that E[z 1 , Z 1 ] is equal to: We then reach Eq. (35) by minimizing this expression subject to |z 1 | < 1 and |Z 1 | > 1 [89].
a bond dimension of χ = 2 N , where N is the number of wires connecting unit cells. Hence, the naive grouping of gates into a repeating unit element would give log 2 χ = 12. The first step in simplifying the MPS representation is shown by the arrows in Fig. 11(a), which leads to the circuit in Fig. 11(b). This consists of commuting gates into the following form |ψ = n M (2) n M (1) n U n P n | · · · ↓↓↓ · · · . (C7) This is possible using the graphical rules in Fig. 5(b), namely, that gates commute past each other so long as symbols of the same kind are not acting on the same wire. After this first step, we have already reduced the bond dimension, and by grouping with respect to the repeating element in Fig. 11(b) we find that log 2 χ = 5 is one less than the support of M (2) n . This bond dimension can be further reduced by making use of the algebraic relations for U n , described in Appendix C 1. We use these relations to pull the unitary gate U n past M  (1) n+2 P n | · · · ↓↓↓ · · · , where in this caseM (k) n = 1 − ia k γ n γ n+k+1 . This is shown by the arrows in Fig. 11(b) leading to the circuit in Fig. 11(c).
Finally, by using the commutation relations of the gates again, we bring the gates into the order shown in Fig. 11(d), which corresponds to |ψ = n M (2)
General case: For the general case with p ≥ 0 we can follow the same steps. First, we are able to commute gates into the form |ψ = n M (d) n . . . M (2) n M (1) n U n (P n ) | · · · ↓↓↓ · · · .
(C10) We place P n in brackets to indicate that it is only there in the case that p is odd. Then, we use algebraic relations to pull U n past the layers of M (k) n , reducing the support of these gates. Finally we simply commute gates so that we can group gates by wire, where the left-most operator of all gates except P n appears on a single wire, preceded by the right-most operator of P n (if it is there).
Let us now argue how the bond dimension changes when we do this. Firstly, when p ≥ 0 we have range(H) = 2d + p. Consider the initial grouping of gates Eq. (C10), with general d. The support of M (d) n is d + p + 1, and we see that the bond dimension of Eq. (C10) satisfies log 2 χ = d + p. The support of U n is p/2 , and so by using the algebraic relations for U n as described, we end up withM (d) with support d + p + 1 − p/2 . Finally, if p is odd then we have to include the projectors P n which increase the support of the unit circuit element by 1. Therefore the total support of the reduced circuit is d+p+1− p/2 . The bond dimension is one less than this since one of the wires is physical (or projected on |↓ ). In conclusion, we have an MPS construction with log 2 χ = d + p/2 = range(H)/2 . (C11) b. Case (ii): p < 0, d ≤ |p| Example: Let us now consider d = 3 and p = −6, that is, a polynomial of the form f (z) = z −6 (s 0 + s 1 z + s 2 z 2 + s 3 z 3 ) 2 , which corresponds to a skeleton through the phase diagram of the Hamiltonians of the form The circuit construction for the ground state is shown in Fig. 12(a). In this case, range(H) = 6, and so Eq. (C3) states that the corresponding MPS has a bond dimension with log 2 χ = 6/2 = 3. However, the naive grouping of gates into a repeating unit element would give log 2 χ = 15.
Similarly to the previous example, the first step is to commute the gates past each other to get the grouping in Fig. 12(b), reducing to log 2 χ = 5. Then we use the algebraic relations for the unitary gate U n to bring it past the corresponding M (1) n , as shown going between Fig. 12(b) and (c). (To be precise, n is the site such that the circles or γ operators of U n and M (1) n agree.) Bringing U n past M