Gaussian matrix product states cannot efficiently describe critical systems

Gaussian fermionic matrix product states (GfMPS) form a class of ansatz quantum states for 1d systems of noninteracting fermions. We show, for a simple critical model of free hopping fermions, that: (i) any GfMPS approximation to its ground state must have bond dimension scaling superpolynomially with the system size, whereas (ii) there exists a non-Gaussian fermionic MPS approximation to this state with polynomial bond dimension. This proves that, in general, imposing Gaussianity at the level of the tensor network may significantly alter its capability to efficiently approximate critical Gaussian states. We also provide numerical evidence that the required bond dimension is subexponential, and thus can still be simulated with moderate resources.

The growing complexity of quantum many-body wavefunctions with increasing system sizes has motivated the development of variational classes of states. By exploiting simplifying features of a given problem, ansatz states can help optimize numerical resources, as well as provide an insightful new perspective into the inner workings of quantum correlations in these systems. For instance, Gaussian states and their correlation matrix formalism greatly facilitate computations involving non-interacting particles. On the other hand, tensor network states have become an essential theoretical framework and numerical toolbox for quantum many-body physics, excelling at the representation of area law and similarly low-entangled states [1].
For systems of free (or weakly interacting) fermions, both classes can be combined to give rise to Gaussian fermionic tensor network states [2]. Relevant examples include the ground state of the Kitaev-Majorana chain, and models of topological insulators and superconductors [3][4][5][6]. In the 1d case, the resulting tensor network is the Gaussian fermionic matrix product state or GfMPS, which has been shown to outperform non-tensor network based methods in free fermion computations for very large systems [7]. This motivates the question about the expressivity of GfMPS, namely what kind of free fermionic states can be efficiently described by them.
In the case of general MPS, it was proved in [8] that an efficient approximation (i.e. one with bond dimension growing at most polynomially with the system size N ) exists whenever a certain Rényi entropy is bounded by O(log N ). This established the usefulness of MPS to approximate states with at most a logarithmic violation of the area law, including the ground states of both gapped and gapless (critical) local Hamiltonians. In the setting of fermionic chains, it also applies to fermionic MPS (fMPS). However, it is not known whether an analogous result holds for GfMPS whenever the state being approximated is Gaussian with similarly bounded entropies.
Here we answer this question in the negative: we provide a simple counterexample in the form of a critical hopping fermion Hamiltonian, whose Gaussian ground state can be efficiently approximated by fMPS but not by GfMPS. The proof of this last fact combines a rigorous bound on the error incurred by a fixed rank Gaussian approximation of a Gaussian state with specific knowledge of the entanglement structure of the target ground state, obtained from asymptotic Toeplitz determinant theory. Furthermore, we provide evidence, both from conformal field theory arguments and numerical results, that the required bond dimension scaling for a good GfMPS approximation is nevertheless subexponential. This makes the question about the existence of an efficient GfMPS approximation a hard one to settle numerically, which motivated our pursuit of an analytical proof.
Our result disproves the somewhat intuitive assumption that the most bond dimension-efficient approximation to a Gaussian state would come from a Gaussian tensor network. This can be relevant when optimizing resources for computational applications, though it should be noted that the savings due to having access to the correlation matrix formalism may well compensate the extra bond dimension derived from Gaussianity. Additionally, there are other situations where Gaussian tensor networks have faced difficulty approximating Gaussian states, as is the case for ground states of local, gapped, quadratic Hamiltonians displaying chiral topological features [3,4]. In this context, our findings leave the door open to the existence of better, non-Gaussian tensor network approximations that bypass the no-go results.
Model -We consider a periodic chain of length N with a single fermionic mode a i , a † i per site, satisfying the usual canonical anticommutation relations, and study the free hopping Hamiltonian at half filling, where we have defined the momentum modes in the usual form, a k ≡ 1 √ N N j=1 e ikj a j with k ∈ 2π N Z ∩ (−π, π]. In this basis H is diagonal and its ground state, which is Gaussian due to H being quadratic, can be determined by filling the negative energy modes below the Fermi mo-arXiv:2204.02478v1 [quant-ph] 5 Apr 2022 mentum, k F = π/2 [28]. This is encoded in the momentum space correlation matrix, where Θ denotes the Heaviside step function. The position space correlation matrix C ij ≡ a † i a j can then be obtained by an inverse Fourier transform. It exhibits power-law decays as befits a gapless model (its explicit form can be seen in Appendix C). Note that due to particle number conservation, a k a q = a i a j = 0.
The entanglement structure of this state, which will be key to the results presented next, can be obtained from its correlation matrix. Given a bipartition of a pure Gaussian state into complementary regions R,R, we can find a basis of modes on each subsystem such that the state decomposes as the tensor product of entangled fermion pairs [9]. How entangled these pairs are is given by the spectrum of the correlation matrix of either subsystem, which is nothing but the corresponding submatrix of the global correlation matrix, For convenience and notational unity we will work with the eigenvalues of 1], and call the Gaussian entanglement spectrum [29]. The Rényi entropy S α then splits as a sum of contributions from each entangled pair, where so that the entanglement decreases with |λ| from λ = 0 (maximally entangled state) to λ = ±1 (product state). For the ground state of H, the leading scaling of the Rényi entropy of an interval of size L can be seen to be logarithmic [10][11][12], which is consistent with it lying in the universality class of the free boson conformal field theory (CFT) with central charge c = 1 [13,14]. Efficient approximation with fMPS -A fermionic matrix product state (fMPS) [2] is defined in terms of a series of so-called fiducial states of f physical fermions and 2χ virtual fermions. The state represented by the fMPS is obtained by contracting the virtual fermions, i.e. projecting them onto maximally entangled pairs (see Theorem 1. The family of ground states of (2) is efficiently approximable by fMPS.
Proof. Using the Jordan-Wigner transformation, we can map our system to a spin chain model, the XX model. Since we have a logarithmic bound (7) for the Rényi entropies, in particular for α < 1, it follows from Lemmas 1 and 2 in [8] that the ground state of the XX model is efficiently approximable by MPS. Then, undoing the Jordan-Wigner transformation, we can find an efficient fMPS approximation for the fermionic model (this is done in detail in Appendix A).
In general, understanding the entanglement structure of quantum states is key to obtain both approximability and inapproximability results, since the bond dimension of an MPS has a clear interpretation as the maximum Schmidt rank for a bipartition of an MPS into connected subsystems. In order to compare it with an analogous result in the next section, we cite here the following Lemma 1 (Low Schmidt rank approximation). Let |Ψ by a bipartite quantum state with Schmidt spectrum {s j } n j=1 in descending order. Then for any bipartite state |Ψ of Schmidt rank at most r, and the bound is tight: the optimal |Ψ can be found by truncating the Schmidt decomposition of |Ψ .
Lemma 1 is a consequence of the Eckart-Young-Mirsky theorem [15,16], which states that the optimal low rank approximation to a given matrix comes from truncating its singular value decomposition. It leads to a lower bound for the error of an MPS approximation [31] (used for instance in some of the inapproximability proofs in [17]).
No efficient approximation with Gaussian fMPS -A Gaussian fMPS (GfMPS) is an fMPS for which all fiducial states are Gaussian. Since the contraction operation projects each pair of virtual modes onto a Gaussian state, the maximally entangled pair, the global state after contraction is necessarily Gaussian. The contraction of GfMPS tensors can be done at the level of correlation matrices via Schur complements [7,18] (see Appendix D).
An efficient approximation in terms of GfMPS can be defined analogously to the fMPS case. Then, we have our main result as Theorem 2. The family of ground states of (2) is not efficiently approximable by GfMPS.
The proof of Theorem 2 will follow from two lemmas. The first one is the Gaussian version of Lemma 1. Given the Gaussian entanglement spectrum {λ j } of a bipartite state, we call the number of eigenvalues λ j = ±1 its Gaussian rank. We then have Lemma 2 (Low Gaussian rank approximation). Let |Ψ be a bipartite Gaussian state, with Gaussian entanglement spectrum {λ i } n i=1 , ordered so that |λ i | ≤ |λ i+1 |. Then for any Gaussian state |Ψ of Gaussian rank at most r, where S trunc , and the bond is tight: the optimal |Ψ can be found by truncating the Gaussian singular value decomposition of |Ψ .
We were not able to find a proof in the literature so we provide one together with related results in Appendix B. Lemma 2 can be used to lower bound the error of a GfMPS approximation to a given state, since the Gaussian rank of any GfMPS divided into two connected subsystems is upper bounded by χ. We then need information on the entanglement spectrum of our target state, which is provided by Lemma 3. For the ground state of (2), let I L,N (µ) be the number of eigenvalues λ from the Gaussian entanglement spectrum of an interval of size L in a chain of N sites that satisfy |λ| < µ, and let c > 0. Then there exists µ < 1 such that The proof of Lemma 3 can be found in Appendix C. It starts by proving the equivalent property for the Gaussian entanglement spectra of the infinite chain: in the thermodynamic limit, we can exploit the theory of Toeplitz determinants to lower bound the corresponding I L,∞ (µ) function in the asymptotic regime. Then we use standard inequalities to show that the difference between the finite and infinite chain correlation matrices is bounded in trace norm. This ensures that their respective spectra are distributed similarly enough so that Lemma 3 follows.
Proof (of Thm. 2). Suppose there exists a GfMPS approximation with polynomial b.d. D(N ). Then we can find c > 0 such that χ(N ) = log 2 D(N ) ≤ c log N . Thanks to Lemma 3, we know there exists some µ < 1 such that I L,N (µ) > (c + 1) log N , and we have which diverges as N → ∞. Since |Ψ N has Gaussian rank bounded by χ(N ) across the bipartition, Lemma 2 implies that the overlap between the ground state and its GfMPS approximation goes to zero as the system size increases. Thus, by contradicion, any approximation with bounded error must have χ(N ) growing faster than logarithmically, and consequently D(N ) grows superpolynomially.
CFT argument -The techniques used in the proof of Theorem 2 cannot be applied to obtain better lower bounds, or upper bounds on the required b.d., for which more accurate knowledge of the Gaussian entanglement spectra of finite chains would be needed. Here we provide evidence that this b.d. is subexponential. First, a heuristic argument is made, based on conformal field theory (CFT). In the next section, we present some numerical results.
The low-lying entanglement spectrum of a critical model is know to behave universally according to the underlying CFT. With our notation, for an interval of L sites in an chain of N sites, we have [19][20][21][22] |λ n | tanh π 2 2 where is the effective length of our interval in units of some UV cutoff a, and the ε n are fixed by the CFT. In our model, this is the free compactified boson CFT, and (More precisely, these numbers characterize the spectrum of scaling dimensions of an associated boundary CFT.) The CFT spectrum (13) not only satisfies the condition in Lemma 3, but it also allows us to estimate the tail contribution to the ∞-Rényi entropy, Thus, if the CFT prediction were exact for the whole spectrum, by Lemma 2 the required scaling for χ is for some proportionality constant = ηN . Here we used the fidelity error attached to a single bipartition of the system. Our experience from MPS is that we should account for all bipartitions, with different L/N and thus different η [8]. We do this coarsely by replacing → /N . This results in which is subexponential.
Numerics -We have also performed some numerical studies of this problem which seem to partially confirm the scaling from (16). We introduce them briefly here, and elaborate on them in Appendix D. Essentially, we defined a subclass of translation invariant GfMPS (which we call ladder GfMPS), which exactly represent states with the following occupation number in momentum space: for p, q arbitrary real odd monic polynomials. Clearly, to reproduce (3), we need to have p (resp. q) supported mostly inside (resp. outside) the Fermi surface. We tried several choices, the best of which is represented in Fig. 2.
There we have used δ ≡ Ψ MPS |H fb |Ψ MPS − E 0 , with H fb the flat band Hamiltonian with the same ground state as H, and E 0 its ground state energy, as a proxy for the fidelity error , which it upper bounds. Further energy optimization (as done in [23]), which we did not pursue, could improve the results. Fig. 2 also shows the estimation (17) for η ∼ 1.3 (resulting from optimization), which gives an idea of the scaling of our numerical results.
Toy model -Finally, we present a simple toy model of a family of Gaussian states, with no reference to a Hamiltonian, that are efficiently approximable by fMPS but not by GfMPS. We begin by explaining the intuition behind it. Lemmas 1 and 2 highlight the differences between Gaussian and non-Gaussian truncation of a bipartite state [32]. This dichotomy is also reflected in different bounds for the truncation errors, G and N G resp., in terms of the α < 1 Rényi entropy. Assuming G , N G 1, we have, For the non-Gaussian case this is Lemma 2 in [8]. The Gaussian case can be proved similarly by minimizing the entropy over all possible Gaussian entanglement spectra with a fixed Gaussian truncation error. Since χ ∼ log D, we see that the bond for G is much weaker, i.e. there exist states with small S α and N G but high G for the same bond dimension [33]. These states are usually characterized by Gaussian entanglement spectra involving many low-entangled pairs. We exploit this in our toy example, leveraging the addition of entangled pairs as we grow the system by making their entanglement increasingly weaker. Consider a family of states |ψ N on rings of N sites obtained by distributing ν(N ) entangled pairs between opposite sites in the chain, each with strength λ N , where for some β > 0. The α < 1 Rényi entropy is maximal for a bipartition of the ring into two equal halves, since all pairs contribute to give which is upper bounded by O(log N ) for α ∈ β 1+β , 1 . Therefore, for any β > 0 an efficient fMPS approximation exists. However, for any GfMPS approximation |ψ N , the error bound (10) reads Thus if | ψ N |ψ N | ≥ 1 − , then for large enough N and small enough we have so that D(N ) > ∼ N (log N ) β . This scaling can be easily seen to be sufficient (since χ(N ) = ν(N ) allows for an exact representation), so that the required b.d. to approximate |ψ N is superpolynomial but still subexponential.
Acknowledgments -A.F.R. would like to thank F. Ares for introducing him to Toeplitz asymptotics a few years ago. A.F.R. is supported by the Alexander von Humboldt Foundation. I.C. acknowledges funding from the ERC Grant QUENOCOBA (No 742102) and the DFG through the DACH Lead Agency Agreement Beyond C (No 414325145). Following [2], a fermionic MPS with f = 1 physical fermions and 2χ = 2 virtual fermions per site can be defined by means of a series of "projectors" where ⊕ denotes the sum in Z 2 , j is the position index along the chain, [A j ] k lm is a coefficient tensor, and a j , α j , β j are the physical, left virtual and right virtual fermionic modes at site j, respectively. The restriction on the indices ensures a well-defined fermionic parity p for Q j . One also needs to define the operators which generate entangled pairs of virtual fermions from the vacuum. The fMPS state is then defined by the action of both sets of operators on the global vacuum, followed by projecting out the virtual fermions: The generalization to larger physical and/or bond dimensions is straightforward. If we now interpret the fermionic Fock space as the Hilbert space of a spin 1 2 chain (as is done implicitly by the Jordan-Wigner transformation), one can see that |Ψ can be obtained from a standard "spin" MPS, with local tensors that coincide with [A j ] k lm , possibly up to signs. Consequently, the set of 1d fMPS states coincides with those obtained from MPS that satisfy: (i) their bond dimension is a power of 2, and (ii) each tensor has a well-defined parity.
In the case we are interested in, the Jordan-Wigner transformation maps the fermionic Hamiltonian H from (2) to the XX model spin chain, for which the results from [8] apply, since all its Rényi entropies grow logarithmically (see Eq. (7)). The existence of an MPS approximation with polynomial bond dimension then follows. Additionally, this spin MPS has a global parity symmetry represented by the action of the product of Z operators on each physical spin. To see that this implies the existence of an fMPS approximation to the ground state with polynomial b.d., we show that we can make the spin MPS satisfy conditions (i) and (ii) without a substantial increase in bond dimension. Let D be the bond dimension of the spin MPS, and [B j ] k lm denote its tensors. Then there is q such that D ≤ 2 q < 2D, and we can embed the D × D × 2-dimensional MPS tensors into 2 q × 2 q × 2-dimensional ones without changing the state just by padding each tensor with zeros, i.e. by extending the range of l, m to 2 q , and defining the additional tensor elements as 0. Thus we can satisfy (i) with less than twice the original bond dimension. To get condition (ii), we modify the tensors by adding an additional pair of indices l , m ∈ {0, 1} such that where |l|, |m| denotes the parity of the corresponding index. The new tensors are individually parity symmetric, so that (ii) holds, and they can be seen to generate the same state as the original ones (which is only possible because said state has global parity symmetry). Thus it follows that there exists an fMPS approximation to the ground state of H, with less than four times the bond dimension of the spin MPS approximation to the XX model ground state, which therefore grows at most polynomially.

Appendix B: Gaussian bipartite state overlaps and Gaussian entanglement spectrum
In this Appendix we prove Lemma 2 as a corollary to a more general result. For convenience, we work here in the Majorana representation, introducing Majorana operators Finally, we introduce the Gaussian singular value decomposition (Gaussian SVD) [9], which states that, given Γ the correlation matrix of a pure bipartite Gaussian state on two subsystems of 2n Majorana fermions each [34], we can find O, Q ∈ O(2n) that block diagonalize Γ, where the 4 × 4 blocks are given by and the θ j can all be chosen to lie on the first quadrant, 0 ≤ θ j ≤ π 2 , in which case we have cos θ j = |λ j |. The θ j are another possible way to write the Gaussian entanglement spectrum. Indeed, W (θ) is the correlation matrix of a pair of fermionic modes, which is in a product state for sin θ = 0 (|λ| = 1) and maximally entangled whenever cos θ = 0 (λ = 0). The number r of entangled pairs (sin θ = 0) is what we called in the main text the Gaussian rank. Now we are ready to state and prove the following Theorem 3. Let |Γ , |Γ be pure bipartite Gaussian states on two subsystems of 2n Majorana fermions. Let their correlation matrices be Γ,Γ and their Gaussian entanglement spectra be given by {θ j } n j=1 , {θ j } n j=1 respectively. Then, and the bound is tight (it is reached for some Γ,Γ).
Proof. We follow a strategy inspired by Thm. VI.7.1 in [24]. Γ,Γ will be of the form for some O,Õ, Q,Q ∈ O(2n). To begin with, we shall assume that We are seeking to upper bound where we have used (B4) and the purity condition Γ −1 = −Γ. In other words, our problem consists in determining We know the maximum exists since we are optimizing over a closed manifold. Further, we can assume O, Q = 1, which amounts to fixing the mode basis on which we express our states and does not affect their overlap.
Assume that (Γ,Γ) constitute an extreme point of the target function. This implies that no infinitesimal change in the matrixΓ will change the overlap, that is, (B15) By differentiating, and then using det(Γ +Γ) > 0 (since we are looking for maxima) and the cyclicity of the trace, we arrive at , which is skew-symmetric, have the following block structure (according to the bipartition of our states), where Condition (B22) can be seen to imply which, thanks to our assumptions about the θ i , is enough to force for some b i ∈ R. Now we go back to (B19) and write hence, where the inverse of 1 − B exists since det(1 − B) = det(1 + B T B) > 0. Defining β i ≡ 2 arctan b i , the expression above yieldsΓ which can be checked to be consistent with all the conditions derived before, in particular (B31) In conclusion, the pairs Γ,Γ with maximal overlap for fixed spectra satisfy that Γ andΓ are simultaneously singular-value-decomposable, by which we mean there exists a basis in which for some permutation σ. The statement of the theorem then follows from simply computing the overlap of these states, and extends to the case of general {θ i ,θ i } by a continuity argument.
The case described in Lemma 2 follows as a corollary to the previous theorem by forcing all but r of theθ i to be equal to 0. It can then be seen that the optimal choice for the remaining ones is for them to equal the r largest θ i (the most entangled pairs) and for the permutation σ to match them accordingly, so that the maximum overlap is given by (10), once we express the Gaussian entanglement spectrum back in terms of λ j . The bound is tight since an approximation with such an overlap can be obtained by performing the Gaussian SVD of the target state and setting all but the r largest θ i to 0 (i.e., all but the r smallest |λ j | to 1).
Appendix C: Proof of Lemma 3 As we advanced in the main text, we begin by proving a corresponding result in the thermodynamic limit: Lemma 4. For the ground state of (2) on an infinite chain, let I L,∞ (µ) be the number of eigenvalues λ from the Gaussian entanglement spectrum of an interval of size L that satisfy |λ| < µ, and let c > 0. Then there exists µ < 1 such that as L → ∞.
Proof. Let C L,∞ be the correlation matrix of the interval of length L, and V L, , and let f (z) be a holomorphic function on a domain that includes the interval [−1, 1] where all the eigenvalues {λ i } of V L,∞ lie. Since we have it follows from Cauchy's integral formula that where C is a contour within the domain of holomorphicity of f encircling the interval [−1, 1]. Since V L,∞ is a Toeplitz matrix with an adequate symbol, the asymptotic value of D L (z) as L → ∞ is given to us by the Fisher-Harwig conjecture, in particular by a subcase thereof which was proven by Basor [25]. This property has been exploited for various computations in the XX model [11]. In our particular case, it tells us where by ∼ we mean both sides are equal up to O(1) terms that do not grow with L. The right hand side of (C3) then reads Using complex variable techniques, this finally results in which is a strong constraint on the distribution of eigenvalues. It hints at the fact that they are asymptotically distributed with a density 2 log L/(π 2 (1 − λ 2 )) along the interval [−1, 1], with the rest of them (a number of order L) eventually clumping at the endpoints. We are now in a position to bound the function I L,∞ (µ). It can be written as a sum over eigenvalues, with f the indicator function of the interval [−µ, µ], which of course cannot be extended to a holomorphic function. Still, to get intuition, the result would be (C9) and since the coefficient of log L diverges as µ → 1, we would have the result. To make a proper statement, we use the functions which lower bound the indicator function Θ(µ − |λ|) and are holomorphic on a disk containing [−1, 1]. Thus we can assure, and since the coefficient of log L on the rhs still diverges as µ → 1, the result follows. Now we will show that the spectra of the correlation matrices for the finite and infinite chains are close enough that Lemmas 3 and 4 imply each other. Denote the Frobenius norm by · 2 and the Schatten 1-norm (or trace norm) by · 1 . We then have, Lemma 5. Let C L,N , C L,∞ be the correlation matrices for an interval of L sites of a finite chain of N sites and an infinite chain, respectively, and let L/N = ϕ stay constant as we increase N . Then C L,N − C L,∞ 1 is bounded by a constant.
Proof. Both C L,N and C L,∞ are Toeplitz matrices. Their matrix elements read where m = 2, 1, 0, −1 whenever N ≡ 0, 1, 2, 3 mod 4 respectively. Define L × L Toeplitz matrices T even j , T odd j with elements T even j i,i+r ≡ cos By expanding and collecting terms cautiously in (C12), it can be seen that where a j , b j are the coefficients in the series expansion of the holomorphic functions which are absolutely summable within their disc of convergence (the unit disc). The trace norm of T even j , T odd j can be bounded by using together with the inequality to find T even j 1 ≤ 4j + 2 T even and Going back to (C16) this yields which converges and is thus bounded as N → ∞ with constant L/N .

Finally, we have
Proof (of Lemma 3). We will argue by contradiction. Assume therefore that there is c > 0 such that for all µ < 1, I L,N (µ) ≤ c log L. Since Lemma 4 holds, we can choose µ < µ < 1 such that Recall now the following inequality for the trace norm of the difference of Hermitian matrices, where α i , β i are the eigenvalues of A, B in descending order [24]. We choose In our situation, there are asymptotically at least log L eigenvalues of A that are at least µ − µ away from their associated eigenvalues of B, thus the left hand side of the inequality grows with L while the right hand side is bounded by Lemma 5, a contradiction.
This is useful in the translation invariant case, for which different momenta decouple. At the level of correlation matrices, this implies that the Majorana correlation matrix Γ, is block diagonalized by the Fourier transform F, Our translation invariant GfMPS will be determined by a fiducial state of 2 Majorana fermions, χ left virtual Majorana fermions and χ right virtual Majorana fermions. Note that in this Appendix χ differs by a factor of 2 from χ in the main text, since it counts the number of virtual Majorana fermions. Because we are working with periodic boundary conditions, we can allow χ to be odd. In fact, the parity of χ can have significant consequences for the parity structure of the states in the variational class [23], and in our case, odd χ is actually preferable. We denote the correlation matrix of the fiducial state by where the block structure comes from separating the two physical fermions (A is a 2×2 submatrix) from the virtual fermions (D is a 2χ × 2χ submatrix). The correlation matrix for the GfMPS state is obtained by projecting the virtual Majorana fermions onto entangled pairs, which in momentum space reads [2] Next we define a rail GfMPS, which is characterized by where O 11 is f × f and O 22 is χ × χ. The correlation matrix for the fiducial state of 2f physical fermions and 2χ virtual fermions for the rail GfMPS is defined to be (D8) Therefore, the 2f × 2f correlation matrix G(k) for the resulting GfMPS state is is a unitary matrix, or, in our f = 1 case, a complex phase. In fact, readers familiar with the theory of discrete linear time invariant (LTI) systems may recognize T (z) as the transfer function of a lossless system whose state space representation is given by the matrix O with the blocking from (D7). This analogy could be exploited to import techniques from the LTI system literature to the GfMPS setting. Here we will not pursue it further. It is nevertheless known that T (z) will be a finite Blaschke product, i.e. a unimodular rational function [26], of the form where |η| = 1 and α j are the eigenvalues of O 22 , which can be any conjugation invariant set of complex numbers inside the unit disk [35]. An f = 1 ladder GfMPS is made from juxtaposing two f = 1 rail GfMPS and projecting their respective second physical Majorana fermions onto maximally entangled pairs, to form the "rungs" of the ladder (see Figure 3). The resulting GfMPS has again f = 1 and a bond dimension that equals the sum of those of the rails plus one (from the rungs), χ = χ 1 + χ 2 + 1. Its correlation matrix is given by G(k) = 0 e ik T 1 (e ik )T 2 (e ik ) * −e −ik T 1 (e ik ) * T 2 (e ik ) 0 (D12) What was gained from this construction is that the product T 1 (e ik )T 2 (e ik ) * is now a unimodular rational function with poles no longer confined to the unit disk, and thus more general. It can be written in terms of an arbitrary polynomial and its reciprocal polynomial, and a few additional manipulations lead to the general form of n k we showed in the main text, for p, q arbitrary real odd monic polynomials of degree χ (assuming χ is odd), or equivalently, n k = (1 + cos k) π(cos k) 2 (1 + cos k) π(cos k) 2 + (1 − cos k) θ(cos k) 2 , (D14) for π, θ arbitrary real monic polynomials of degree χ−1 2 . We can then try to guess adequate families of polynomials that make n k close to its exact value from Eq. (3) on the allowed momenta k ∈ 2π N Z∩(−π, π]. We tried expressions based on Fourier expansions of the exact n k and on Chebyshev polynomials, which nevertheless displayed a clearly exponential b.d. Our best results came from picking p (resp. q) so that its zeros are exactly a subset of the allowed momenta that are outside (resp. inside) the Fermi surface. For those selected values, the GfMPS approximation reproduces the target state exactly. Several approaches can be followed to choose which precise momenta we make exact. Choosing all of them next to the Fermi points still leads to exponential b.d., but spreading them logarithmically (so that we still pick exponentially more momenta that are close to the Fermi points), leads to a very well-behaved ansatz family that gives rise to the results shown in the main text.