Bulk-boundary correspondence and singularity-filling in long-range free-fermion chains

The bulk-boundary correspondence relates topologically-protected edge modes to bulk topological invariants, and is well-understood for short-range free-fermion chains. Although case studies have considered long-range Hamiltonians whose couplings decay with a power-law exponent $\alpha$, there has been no systematic study for a free-fermion symmetry class. We introduce a technique for solving gapped, translationally invariant models in the 1D BDI and AIII symmetry classes with $\alpha>1$, linking together the quantized winding invariant, bulk topological string-order parameters and a complete solution of the edge modes. The physics of these chains is elucidated by studying a complex function determined by the couplings of the Hamiltonian: in contrast to the short-range case where edge modes are associated to roots of this function, we find that they are now associated to singularities. A remarkable consequence is that the finite-size splitting of the edge modes depends on the topological winding number, which can be used as a probe of the latter. We furthermore generalise these results by (i) identifying a family of BDI chains with $\alpha<1$ where our results still hold, and (ii) showing that gapless symmetry-protected topological chains can have topological invariants and edge modes when $\alpha -1$ exceeds the dynamical critical exponent.

Introduction.The bulk-boundary correspondence is a central concept in the study of topological phases of matter [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].This relates topologically stable edge effects with topological features of the bulk Hamiltonian.A simple manifestation of this is in certain translationinvariant quantum chains with time-reversal symmetry, where the Hamiltonian on a periodic chain can be used to define a winding number which counts the number of topologically protected Majorana zero modes localised at the edge [1,4,[18][19][20][21].Research on this topic has predominantly focused on the short-range case where lattice Hamiltonians couple sites up to some finite range.In the past decade there has been significant interest in quantum systems with long-range interactions [22,23].This has been motivated by proposals for, and progress in, experimental systems, such as Ref. [24] for effective freefermion chains.Here long-range typically means that couplings decay as a power of the distance [i.e., Hamiltonian terms acting between sites at distance r are O(r −α )].Interesting physical effects have been observed including algebraically localised edge modes and the breakdown of the entanglement area law [25] and conformal symmetry at criticality [26].
Regarding topological edge modes in such long-range chains, most results in the literature concern the canonical Kitaev chain [27] with additional long-range hopping or pairing terms [22,[28][29][30][31][32][33][34][35][36].(For interacting studies see Refs.[37,38].)The long-range Kitaev chain sits in the BDI symmetry class of free-fermion Hamiltonians [4,8,18,39], and it is straightforward to see that for α > 1 the bulk winding number remains well defined [30].Very recently, Ref. [40] treated the free-fermionic phase diagram in great generality and gave a proof that the short-range phase classification is preserved in the longrange case with α > d (in general dimension and symme-try class).Work on the long-range Kitaev chain showed that topological edge modes exist, but only in particular models.This leaves open important questions for topological Majorana zero modes in long-range chains: when do they exist, what is their connection to the bulk invariant, and what are their localisation properties at the edge?
Here, we present the first systematic study of a whole symmetry class, giving rise to a detailed bulk-boundary correspondence in long-range chains.We focus on the exemplary BDI class as mentioned above, although the results carry over for the AIII class [41] which famously includes the Su-Schrieffer-Heeger chain [42].
We show that the bulk invariant corresponds exactly to the number of topological edge modes and give a rigorous method to find the edge mode wavefunctions.Additionally, we find that the bulk string-order parameters for the short-range case continue to reveal the bulk topology.We complement these results by outlining a principle for calculating the finite-size energy splittings for the zero modes in long-range chains, that we call singularity filling.Together with our analysis of the localisation properties of the edge modes, this brings a number of disparate results in the literature into a coherent picture.
The methods we use are from the mathematical theory of Toeplitz determinants (see, e.g., [43]), a key technique in the analysis of the two-dimensional Ising model [44].We expect this approach to long-range chains to be fruitful more generally.
The model.Consider the BDI class of translation invariant spinless free-fermions with time-reversal sym-arXiv:2211.15690v3[cond-mat.str-el]14 Apr 2023 metry: Here are the real [imaginary] Majorana fermions constructed from spinless complex fermionic modes c n on each site.The real coupling coefficients t n are called α-decaying [40] if t n ≤ const(1 + |n|) −α .Assuming absolute-summability of the t n (implied by α > 1) we can solve the closed chain by a Fourier transformation and Bogoliubov rotation (see Appendix A).This information is summarised by the continuous complex function: The eigenmode with momentum k is defined by the phase of f (e ik ) and has energy ε k = |f (e ik )|.Thus, the Hamiltonian ( 1) is gapped when f (z) = 0 on the unit circle.In that case, the argument of f (z) is well-defined, and we have the winding number This is the bulk topological invariant, which cannot change without a gap-closing if we enforce the absolutesummability condition.
Bulk-boundary correspondence and edge mode wavefunction.We now consider the Hamiltonian (1) with open boundary conditions (we keep only the couplings that do not cross the boundary).We first consider the limit of a half-infinite chain, where edge modes have zero energy (later we study finite-size splitting).
In this limit, the edge mode wavefunctions are zeroeigenvectors of a Toeplitz operator, which can be solved using the Wiener-Hopf method.More directly, define a real Majorana zero mode as γ L = ∞ n=0 g n γ n that satisfies [γ L , H BDI ] = 0. Evaluating the commutator gives us a Wiener-Hopf sum equation, which is straightforwardly solved [45] using results of McCoy and Wu [44], leading to: Theorem 1 (Bulk-boundary correspondence) Take a half-infinite open chain H BDI , where the related bulk Hamiltonian has winding number ω and absolutely-summable couplings, then there exist exactly |ω| zero-energy edge modes.
More constructively, writing f (z) = z ω b + (z)b − (z) (here b ± (z) are the Wiener-Hopf factors defined below), then for ω > 0 we have ω linearly independent normalisable real edge modes given by γ For ω < 0 the same results hold upon substituting Here and throughout we use the notation that (h(z) n+1) dz is the nth Fourier coefficient of a function h(z).Key to our result is a canonical form called the Wiener-Hopf decomposition.First define f 0 (z) = z −ω f (z), which is non-vanishing on the unit circle and has a continuous logarithm log(f 0 )(z).We fix the normalisation of H BDI such that the zeroth Fourier coefficient (log(f 0 )) 0 = 0. Then we can always write: where the Wiener-Hopf factor given by b ± (z) = e ∞ n=1 (log(f0))±nz ±n is analytic strictly inside (outside) the unit disk.We note that z ω encodes the winding around the unit circle and hence the topological invariant of the system.Multiplying f (z) by z m shifts [46] the hopping t n → t n−m , such that f 0 (z) defines a topologically trivial 'version' of the system.This is analogous to the trivial insulator and the Kitaev chain being related by a shift.
Theorem 1 extends the bulk-boundary correspondence from the short-range to the long-range case: the bulk winding number counts edge modes everywhere in the space of Hamiltonians with absolutely-summable couplings [(α > 1)-decay implies absolute-summability, but examples like the Weierstrass function [47,48] can be used to construct families with 0 < α ≤ 1].Our result is also constructive: we have the edge mode wavefunction in terms of Fourier coefficients of a particular function.To construct the exact edge mode, one needs to first calculate the Wiener-Hopf decomposition.However, we will see below that this can often be bypassed if one is interested only in the asymptotic edge-mode profile.
In short-range models we expect exponentiallylocalised edge modes, corresponding to roots of f (z) [21,49] (see Appendix C 4).Based on Theorem 1 we see that the localisation follows from analytic properties of the Wiener-Hopf factors.If (b ± (z ±1 )) −1 is analytic to some distance outside the unit circle, we will see exponential decay (this appears consistent with previous such observations in the long-range Kitaev chain at fine-tuned points [33]).Exponential localisation was also observed in Ref. [38], but for a different reason-there the shortrange (parity-odd) edge modes cannot couple to the longrange density-density interactions in perturbation theory due to fermion parity symmetry.In our long-range case, the edge modes are generically algebraically-decaying and guaranteed to be normalisable due to the Wiener-Lévy theorem [44,47].
One can read off b + (z) = Li α (z)/z and b − (z) = zLi α (1/z).Suppose ω = 1, then we have one edge mode with: the second equality is derived using contour integration and known asymptotics for Li α (z) on the real line (assuming α / ∈ N) [50].The analysis is given in Appendix FIG. 1. Finite-size splitting from singularities.(a) As an example of our general results, we consider a long-range chain whose hopping coefficients define the complex function f (z) [Eq.( 2)] with singularities of f (1/z) D and further terms in the asymptotic expansion can be found using the same methods.
For ω = 2, we see we have two edge modes, with the same leading order behaviour.This means we can take the difference n −α − (n − 1) −α = Θ(n −α−1 ), and have a faster decaying strictly localised mode (see Theorem 2).
Singularity-filling for wavefunctions.While the bulk-boundary correspondence of Theorem 1 is our most general result, we can give additional results in a broad class of (α > 1)-decaying models.We say that 1/f (1/z) has singularities at {k s } 1≤s≤r if it has asymptotic Fourier coefficients (1/f (1/z)) n = r s=1 e inks n −Ω ks (a s + o(1)) as n → +∞.We call Ω ks > 1 the order of the singularity at k s , and assume the o(1) term is 'nice', i.e., can be expressed as a sum of inverse powers of n [as is the case in Eq. ( 5)].We also assume that Ω min = min s {Ω ks } / ∈ Z.This implies that 1/f (1/z) has δ 0 = Ω min − 1 continuous derivatives [48,51].
We can find a basis of mutually anticommuting edge modes γ(p) For ω < 0 analogous results hold where we now take γ n → γn and f (1/z) → f (z).
The idea of the proof is as in the ω = 2 example following (5): we take linear combinations of edge modes that cancel the dominant asymptotic term(s), and then use the Gram-Schmidt process (with respect to the anticommutator) to construct anticommuting modes [21].We note that if the Fourier coefficients of the Wiener-Hopf factors themselves have a 'nice' expansion, then singularityfilling will hold with no limiting ν (see Appendix C 2).
Example.The long-range Kitaev chain corresponds to: This model was studied for various choices of couplings in Refs.[22,29,[31][32][33]36].Computing (1/f (1/z)) n gives the asymptotic behaviour of the edge mode wavefunction in the ω = 1 case: agreeing with results in the literature (see Appendix E).
In analogy with the singularity-filling for edge-mode wavefunctions above, we have a conjecture for the finitesize splittings for the edge modes.In this case the levels associated to singularities go up in steps of two.

Conjecture 1 (Splitting from singularity-filling)
Take an open chain H BDI of size L, where the related bulk Hamiltonian has winding number ω > 0 and 1/f (1/z) has singularities as defined above.
For ω < 0 analogous results hold where we replace This conjecture is based on numerical experiments (see Fig. 1) and theoretical results (see below).The underlying theory indicates that for a family f (z) = z ω f 0 (z), there may exist an ω max such that this holds only for ω < ω max .In fact, given Ω min > 5, and an assumption on the spectrum, we can prove the conjecture up to ω max = 3.However, empirically we expect the conjecture to hold more generally, as observed in Fig. 1.
The conjecture allows us to understand how finite-size effects hybridise the edge modes.For ω = 1 we see that the predicted splitting comes from the dominant singularity ε 1 = Θ(L −Ωmin ).Since this has the same asymptotics as the edge-mode wavefunction, this agrees with an intuitive connection between the spatial profile of the wavefunction and the induced splitting from the boundaries (see Appendix F 1) that does not generically hold for the higher-winding case.For ω = 2 we expect to have two edge modes, one with ε 1 = Θ(L −Ωmin ) and one with either ε 2 = Θ(L −(Ωmin+2) ) or ε 2 = Θ(L −Ωnext ), depending on which has the slower decay.In the case of higher winding numbers, our conjecture predicts the hybridisation of the boundary modes, which is not in direct correspondence to the maximally localised basis identified in Theorem 2.
To justify the conjecture, consider models f (z) = z ω f 0 (z) with open boundary conditions; each such model has a corresponding single-particle (block Toeplitz) matrix, with determinant equal to L j=1 (−ε 2 j ), where ε j are single-particle energies.Assuming (α > 1)-decay, it can be shown, using Toeplitz determinants, that for the trivial model f 0 (z) this product is finite in the limit L → ∞, while for ω = 0, the corresponding determinant decays to zero with L (with power depending on ω and Fourier coefficients of 1/f (z)); see also Appendix F. Our method is to use the scaling of this determinant to predict the edge mode splitting.E.g., for ω = 1 we interpret: as predicting a single edge mode with finite-size splitting ε 1 = Θ(L −ν ).For multiple edge modes (and ω > 0), we further assume inductively that the ω − 1 edge modes shared between the models z ω f 0 (z) and z ω−1 f 0 (z) have the same energy splitting power-law in each model, and hence the additional decay in the determinant for z ω f 0 (z) comes from the ωth edge mode [52].This is plausible since for periodic boundaries the models defined by f (z) have spectrum independent of ω, and we expect the system with open boundaries to differ from the bulk only 'near the edge'.With finite-range interactions we believe this could be proved using results about eigenvalues of banded block Toeplitz matrices [53], for long-range chains we take it as an assumption that the scaling to zero with L comes only from edge modes rather than the bulk band.In an earlier work the idea appeared in reverse: utilising the existence of exponentiallylocalised edge modes in short-range chains to predict asymptotics of block Toeplitz determinants [54].
We thus convert the question of finite-size edge mode splitting to a question about asymptotics of Toeplitz determinants.While there are several assumptions required to connect this theory to the edge mode splittings, the underlying singularity-filling picture for Toeplitz determinant asymptotics is in many cases fully rigorous.We outline some of these results in Appendix F; important references are [43,51,55,56].
Novel topological probe.A remarkable consequence is that the finite-size splitting of the lowest energy mode depends on the total number of edge modes.In fact, we can turn this into a probe of ω: by perturbing a short-range chain f s (z) (with winding ω) by a long-range test function, its finite-size splitting exponent will allow us to find ω (note that this is the scaling of the lowest one-particle energy, no further information about the spectrum is required).An example test function would be f LRK (z), with ∆ = 0. Then for the function f (z) = f s (z) + f LRK (z), for small, our picture gives a finite-size splitting L −(α+2(|ω|−1)) .
String-order parameters.We now consider the periodic chain.Define the finite fermion parity string by It is know that the set of O κ (n) form order parameters for the gapped phases in the short-range case [57].In the long-range case we have: Theorem 3 (String order) Consider a gapped (α>1)decaying H BDI , in the thermodynamic limit with periodic boundaries, and write f (z)/|f (z)| = z ω e W (z) .Then: Thus the O κ act as order parameters in the long-range case.The idea of the proof is as follows: the stringcorrelation functions O κ (1)O κ (N ) are Toeplitz determinants generated by z −κ f (z)/|f (z)|.The function f (z)/|f (z)| generates the correlation matrix of the chain, and it was proved in Ref. [40] that for an α-decaying chain with α > 1, the correlation matrix is (α − ε)-decaying for any ε > 0. This is sufficient regularity for us to use the results of Ref. [55] to prove Theorem 3 (see Appendix G).
Gap-closing and edge modes at critical points.For H BDI with finite-range couplings, topological edge modes can persist at critical points [21,58].We give some results in this direction for the long-range case.
Suppose we have a gapless bulk mode with dynamical critical exponent z dyn .In the continuum limit, the dimension of the long-range term in the action δS ∼ ψ(x)ψ(y)(x − y) −α dtdxdy is (z dyn + 1 − α), which is irrelevant for α > z dyn + 1.On the lattice, we hence expect that for gapless models of the form f crit (z) = (z − 1) z dyn f gap (z) (which has the aforementioned lowenergy description if f gap (z) is non-vanishing on the unit circle), the edge modes will be stable as long as f (z) is (α > z dyn + 1)-decaying.Indeed, our Theorem 1 can be adapted to show that this f crit (z) has ω localised edge modes where ω is the winding number of f gap (z).This follows from expanding (z − 1) z dyn in f crit (z), and interpreting this as a sum of (z dyn + 1) gapped Hamiltonians, all sharing the same ω edge modes as per Theorem 1.
The above functional form can arise by interpolating between topologically distinct gapped Hamiltonians.For instance, between two phases with winding numbers ω = 1 and ω = 2, there will generically be a single gap-closing with a linearly-dispersing mode if α > 2.More precisely, if this occurs at momentum k = 0, then f gap (z) := f (z) z−1 should define a gapped model with ω = 1.We can then apply the above discussion to infer the existence of the localised edge mode at criticality.We have confirmed this for an explicit example in Appendix H.
Outlook.We have shown how general analytic methods can be used to establish the bulk-boundary correspondence in a class of long-range chains, and give insights into edge mode localisation and finite-size splitting.This included examples with α < 1 and certain gapless models.Key questions remain within this class: what happens in the general case when α < 1 and the integer winding classification breaks down?Can we establish general stability results in critical lattice models, and do these coincide with our field-theoretic analysis?We expect extensions of analytic techniques used above to provide further insights.Moreover, it is worth exploring how broadly our results can be generalised, including to other freefermion classes (beyond BDI and AIII) [4,8,18] and higher-dimensional models.
The extension to long-range multi-band cases would be interesting, likely requiring block Toeplitz operators.In the short-range BDI and AIII classes, edge modes were constructed in Ref. [49], where the bulk topological index is the winding of the determinant of a chiral block of the Hamiltonian.
From the mathematical side, it would be most interesting to find a proof of the singularity-filling conjecture.It would be interesting to see if this picture generalises beyond the studied cases, perhaps even to interacting models with algebraically decaying edge modes, and whether their finite-size splitting also depends on the value of the topological invariant.
Then we can use the usual method of solving such chains via Fourier transformation and Bogoliubov transformation (see e.g., Ref. [21]), leading to the complex function f (z).As a consequence of the absolute-summability, in the thermodynamic limit we can write continuous functions ε k and ϕ k such that on the unit circle f (e ik ) = ε k e iϕ k .Then the Hamiltonian has the diagonal form: In this expression, c k = n e −ikn c n are the Fourier transformed spinless fermion operators c n = (γ n + iγ n )/2.
Appendix B: Singularities of f (z) and related functions

Defining singularities via Fourier coefficients
In the main text we introduce the idea of singularity-filling, for singularities of certain functions defined on the unit circle in the complex plane.These singularities correspond to particular momenta 0 ≤ k s < 2π, or, equivalently, to points on the unit circle e iks .In general, we say that a function h(z) on the unit circle has singularities at {k s } 1≤s≤r if it has an asymptotic Fourier expansion of the form h n = r s=1 e inks n −Ω ks (a s + o(1)).In the main text we suppose that the o(1) terms are all of the form n −b for b > 0. One may consider generalisations, such as allowing terms of the form log(n) a n −b , for some a ∈ Z and b > 0. We explain below that the proof of Theorem 2 can accommodate this particular generalisation.
This definition of singularity, based on Fourier coefficients, can be related to other notions of analytical singularity.For example, suppose h(z) has branch point(s) on the unit circle at e iks , and that we can analytically continue the function to the complex plane up to some branch cuts.When we compute asymptotic Fourier coefficients, we are dominated by integrals near the branch points.Then, supposing an appropriate expansion at the singularity, using Watson's lemma [63,64] we find a dominant contribution (at each singularity) of the form e iksn n −Ω ks (a particular example of this is the calculation following Eq.D3 below).Another relevant definition of singularity is a discontinuity in some derivative of the function h(z).More precisely, we can characterise the smoothness of the function h(z) according to the number of continuous derivatives.Then we have well-known results relating this smoothness to the asymptotic decay of the Fourier coefficients; see, for example, Ref. [48].
Since our results in the main text depend directly on certain Fourier coefficients, we choose to use this definition of singularity for clarity.In analysing a particular problem with a chain corresponding to a function f (z) that has some analytical singularity, one needs to then justify how this is reflected in the asymptotic Fourier expansion.Whether this is straightforward depends upon the particular choice of f (z), but there are many general results available [43,48].

Relationship between singularities of f (z) and the Wiener-Hopf factors
Our main results depend on several different, but closely related, functions.The function f (z) = ∞ n=−∞ t n z n corresponds directly to the Hamiltonian, and the dominant asymptotic decay of the Fourier coefficients of f (z) tells us the algebraic decay of the coupling coefficients.Note that f (z) and z k f (z) necessarily have the same singularities, since we simply shift t n → t n−k and this does not change the asymptotic Fourier coefficients.
For ω > 0 [ω < 0] the edge mode wavefunctions depend on Fourier coefficients of the inverse Wiener-Hopf factor Moreover, based on the Toeplitz determinant theory that underlies Conjecture 1, the edge-mode splittings depend on the asymptotic Fourier coefficients of m(z (this is explained in greater detail in Appendix F).These functions are clearly closely related, and this can be made quantitative.
Let us then consider an (α > 1)-decaying Hamiltonian, with corresponding f (z).These functions f (z) are a subset of the class considered in Ref. [55], and we can make a corresponding analysis.Following Ref. [55], denote the Fourier coefficients of b + (z) by r n , and the Fourier coefficients of b − (1/z) −1 by q n .Then, from the Wiener-Hopf decomposition, we have that r 0 = q 0 = 1 and r −m = q −m = 0 for any m ∈ N.Moreover, r n and q n are absolutely summable and we have: Similarly denote the Fourier coefficients of f 0 (1/z) −1 by s n (which is in general doubly-infinite).Then we have: Note that this calculation is exact for 1/f 0 (1/z), but if we instead take 1/f (1/z) = z −ω /f 0 (1/z) the Fourier coefficients are simply shifted.
Finally, b + (z) 2 has the same properties as b + (z), with Fourier coefficients rn .We can then write Let us analyse (B2), with analogous conclusions holding in the other cases.Suppose, as in the main text, that we have an asymptotic expansion for s n with certain singularities where we say the o(1) term is 'nice', i.e., that each term is an inverse power of n (we may also have further subdominant terms that decay faster than any power of n, we will usually suppress them below).Now, using the absolute summability of the r n , we see that q n has an asymptotic expansion with identical singularities {k s } (and corresponding orders {Ω ks }) to s n , we simply renormalise the coefficients in the expansion.
To justify that the orders of the singularities are the same, note that since f (z) corresponds to a gapped Hamiltonian, b + (z) cannot vanish on the unit circle.
For our purposes in Theorem 2 and Conjecture 1, we also want the o(1) terms here to have a nice dependence on n.We now show that the expansion is in inverse powers of n up to some cut off that depends on Ω min , the dominant singularity.In particular, we will now show that for some known constants A j , and δ 0 = Ω min − 1 .The same conclusion will hold for m n , by analogous calculations.
To prove this, we first recall that functions in the class C β have n = β continuous derivatives.Moreover, our assumption on the singularities of 1/f (1/z), where Ω min is the order of the dominant singularity, implies the following.[51], and in particular have δ 0 continuous derivatives.Let us now revisit the crucial term in the expansion of q n : ∞ j=0 e iksj r j n n + j This is in a form amenable to Watson's lemma [63,65], leading to the following asymptotic expansion: for some constant coefficients A k that depend on derivatives of b + (z).This establishes (B6) above.Consider the BDI Hamiltonian with open boundaries.In general [66], we can diagonalise by finding raising and lowering operators χ ε = L−1 n=0 a n γ n + b n γn that satisfy [H BDI , χ ε ] = 2εχ ε .Evaluating the commutator reduces to the mathematical problem of finding the eigenvectors of a block Toeplitz matrix, for which analytical results are available only in special cases.Note that Toeplitz matrices are matrices that are constant along diagonals.These constants are determined as Fourier coefficients t n of a generating function t(z).A block Toeplitz matrix has the same structure as a (scalar) Toeplitz matrix, but the scalar constants on each diagonal are replaced by constant matrices of fixed size.
One such solvable case is that of exact zero modes; then ε = 0, and the problem reduces to finding eigenvectors in the kernel of scalar Toeplitz matrices.There do exist a variety of results for asymptotic behaviour of eigenvectors of scalar Toeplitz matrices [67,68]; for a review of the field see [69].However, for topological Majorana zero modes, the splitting is generically exactly zero only in the infinite system size limit.

a. Wiener-Hopf sum equations
In their textbook on the Ising model [44], McCoy and Wu solve the following Wiener-Hopf sum equation: subject to the condition n∈Z |c n | < ∞, and solutions are sought with bounded norm, i.e., n∈Z |x n | < ∞.
For our application, y n = 0, and we give the results for that case.Define c(z) = n∈Z c n z n .Then, assuming c(z) does not vanish on unit circle, we have the Wiener-Hopf decomposition c(z) = z ν β + (z)β 0 β − (z), and let us fix the overall normalisation so that β 0 = 1.Note this decomposition exists and each function has an absolutely convergent Fourier series due to the Weiner-Lévy theorem [44,47].
The general solutions of (C1) for y = 0 are as follows: We thus see that for ν ≥ 0 there are no non-trivial solutions, while for ν < 0 we have |ν| solutions.Note the fixed chirality of this problem, it is always negative winding allowing solutions.

b. Application to edge modes
Consider the half-infinite OBC Hamiltonian H BDI = i n≥0,m≥0 t m−n γn γ m , assuming that For real chiral edge modes we have γ L = n≥0 α n γ n .These satisfy [H BDI , γ L ] = 0. Calculating this commutator we find: where tα = t −α .We thus see that if this commutator vanishes, then the α m must be solutions to (C1), for the choice By considering dependence on z and 1/z we have b ∓ (1/z) = β ± (z), we reach the first part of Theorem 1.Note that to prove the independence of the edge modes, we use b − (1/z) −1 = e ∞ k=1 V −k z k has all negative Fourier coefficients equal to zero [70].The second part of Theorem 1 follows straightforwardly by considering the inverted chain.Then we take f (z) → f (1/z) and switch γ and γ.Alternatively one can repeat the above Wiener-Hopf calculation after inserting the ansatz for an imaginary chiral edge mode.
As an aside, suppose we did not restrict to chiral zero modes, and take the following ansatz χ L = n (α n γ n +β n γn ).Then calculating the commutator gives two independent problems of the form (C1), one for f (z) and one for f (1/z).Hence, depending on winding number, at least one of them will have no non-trivial solutions and we are back in the chiral case.
Note that a Majorana edge mode is normalisable if n≥0 |g n | 2 < ∞.The proof of Theorem 1 leads to a stronger conclusion than stated: in fact the edge modes given are the only edge modes that exist satisfying the condition n≥0 |g n | < ∞.Hence, the discussion above based on results of [44] does not immediately exclude 'accidental' (non-topological) localised edged modes that are sufficiently delocalised that n≥0 |g n | → ∞.However, appealing to general results [43,59,60] on invertibility of Toeplitz operators (over the sequence space l 2 ) leads to the conclusion that Theorem 1 does indeed give us all of the Majorana edge modes.

Proof of Theorem 2
Here we prove a stronger form of Theorem 2 given in the main text.The version in the main text is simpler to state, and follows from: where P is a finite set of non-negative reals that contains zero.Define ν 1 , . . ., ν ω by the ω lowest levels E s (n) = Ω ks +n over all singularities s and n ∈ Z ≥0 ('singularity filling').Define also ν = min s {Ω s + p s }.
We can find a basis of mutually anticommuting edge modes γ(r) νr ); here νr = min{ν r , ν }; while q r is equal to q s for the corresponding singularity.
Let us first do an analysis of linear combinations of asymptotic expansions.Suppose we have an expansion: Then we have that: e inks e −imks a s,p n −Ω ks −p (C6) where a s,0 = a s,0 and the other terms can in principle be computed from the expansion of (n − m) −Ω ks −p .By taking a linear combination gn = g n − Ag n−m we can cancel the leading term from one (and only one) of the singularities.Indeed, we simply choose A = e imks to cancel the leading term of the series about the sth singularity.Now we show we can cancel terms inductively according to singularity-filling. Suppose we have a set {g n , . . ., g n } constructed from {g n , g n−1 , . . ., g n−m+1 }, such that each of the terms g Here p 0 (s) is the 'filling' of the singularity s.For j = 0, p 0 (s) = 0 for all s, for j = 1, p 0 (s) = δ st where Ω kt is the minimum of all of the Ω ks and so on.Now, we take g n−m to add to our set, this will be of the form (C6).We can then take a linear combination of g n−m and g n will cancel the next leading term (according to the singularity-filling prescription).Continuing in this way we cancel dominant terms until we reach g n .Since we always cancel the dominant term we are in accordance with singularity-filling.Now, suppose the g n−m are the wavefunction coefficients of our linearly independent zero modes γ (m) L as given in Theorem 1.We can take the linear combinations prescribed above and it is clear that we maintain linear independence.For us to have a good basis of edge modes we also need them to mutually anti-commute.This can be achieved by a Gram-Schmidt process [21]-if we do this in order of fastest decaying to slowest decaying we will preserve the asymptotic decay rates found above.
Note that it may be the case that when taking linear combinations some a s,p vanishes accidentally, this simply means we can find even faster decaying modes.We also need to deal with more general asymptotic expansions that have discrete sets of powers as well as logarithmic terms.
First consider an expansion of the form where p is discrete.An example would be the expansion with a single singularity: with α > 0 (i.e., P = {0, α} in Theorem 4).If we take a linear combination to cancel a 0 it will necessarily also cancel b 0 , so we indeed restrict to the series Ω + n and ignore α, consistent with the claim in Theorem 4.
Consider now logarithmic terms in the asymptotic expansion.Suppose then: ∞ q=0 e inks a s,p,q log(n) qs−q n −Ω ks −p .(C9) Note that now g n = Θ(log(n) qs n −Ω ks ), where s minimises Ω ks .Using that log(n we have that: e inks e −imks a s,p,q log(n) qs−q n −Ω ks −p , (C10) where a s,0,q = a s,0,q .We can then fill singularities inductively as above, where the decay associated to each singularity will be O(log(n To complete the proof, we need to consider the restriction on the sums by p s and q s .The key point is the error term o(log(n) qs n −Ω ks −p s ), where we do not know the explicit n dependence, and hence the behaviour on taking linear combinations.We can apply singularity filling as described above up to the point this term is no longer subdominant.This leads to the ν in Theorem 4.
Having Theorem 4, we can deduce Theorem 2 of the main text using the connection between singularities of 1/f (1/z) and b − (1/z) −1 .
In particular, suppose that 1/f (1/z) has an expansion of the form

Weierstrass chains
Here we present an example of a chain where we establish the bulk-boundary correspondence for t n that decay slower than 1/n.This is perhaps surprising based on previous literature on long-range chains [30], but is straightforward given our condition of absolute-summability.
Take an integer b > 1, and a real number β > 0. The corresponding Weierstrass chain has couplings is gapped and has winding number zero.The series converges absolutely [47,48], and one can see that t n ≤ C/(1 + |n|) β ; i.e., the t n are β-decaying.Hence, we have examples of β-decaying chains for any β > 0. We can then use Theorem 1 to find edge modes for the shifted Weierstrass chains f (z) → z ω f (z).Note that for 0 < β < 1, f (z) is nowhere differentiable [47].So long as we maintain the gap, we can add couplings t n corresponding to another absolutely-summable chain; this means we can find further (β > 0)-decaying chains that are not restricted to the special case studied here where many t n = 0.

Short-range chains
Here we connect Theorem 1 to the known results in the finite-range case [1,20,21].Indeed, in this case, the results reduce to those given in Ref. [21].
If the t n are non-zero for only a finite range, then we have that: The z j are inside the unit circle, and Z k are outside the unit circle.We can read off ω = N z − N p , and for ω > 0 we have that there are ω edge modes and there is a basis where the localisation lengths are set by the ω zeros closest to the unit circle.A proof is given in Ref. [21].
We can get the same result using our Theorem 1. First, up to a factor e V0 that we fix by rescaling the Hamiltonian, we have: For ω > 0 we then use Theorem 1 to identify ω edge modes with wavefunctions given by the Fourier coefficients: for n sufficiently large and where As in the proof of Theorem 2 (and in corresponding analysis in [21]) we can then take appropriate linear combinations to get the claimed localisation lengths.
An analogous discussion holds for ω < 0 and zeros outside the unit circle appearing in b + (z).Moreover, we can use the analysis of gapless models given in the main text (see also below) to see that short-range gapless models have edge modes with localisation lengths determined by zeros of f (z).
One may consider long-range chains as a limiting case of short-range chains, where the interaction range tends to infinity.Then, the degree of the pole and/or the number of zeros on the unit circle increases without bound.The results in this paper apply to fixed Hamiltonians.This means that short-range chains (even with arbitrarily large but finite range) have edge-modes with exponentially-decaying wavefunctions for sufficiently large site index.On the other hand, long-range chains, even with very weak long-range couplings (e.g., α-decaying models with arbitrarily large α), will typically have an algebraically-decaying wavefunction for sufficiently large site-index.Physically we expect that these cases should behave similarly; the starkly different behaviours for a fixed Hamiltonian are a consequence of the particular sequence of limits that we are working in.For short-range chains with a large finite range, the edge mode will have a wavefunction corresponding to (C14), and for a certain values of n (depending on the range of the Hamiltonian) this will approximate the algebraic decay of a long-range chain.Hence, by considering a sequence of finite-range Hamiltonians converging to a long-range model, we expect to see agreement.This is comparable to the approximation of the ground-states of critical spin chains using a sequence of matrix-product states of increasing bond dimension [71].
Aside from certain special cases, we do not have a closed-form for the Wiener-Hopf decomposition (see the next subsection for an example of such a special case), as would be needed to find the exact edge-mode as given in Theorem 1.However, as discussed above, the asymptotic Fourier coefficents of 1/f LRK (1/z) and the corresponding 1/b − (1/z) have the same singularities.Hence, for the topological phase, we can find the asymptotic decay of the wavefunction by calculating large Fourier coefficients, G n , of 1/f LRK (1/z) (one could also use the statement of Theorem 2 directly).
The calculation of these large Fourier coefficients is similar to that given in the previous section and goes as follows.By definition, The integrand is analytic in the same cut-plane as f LRK (z), excluding isolated poles at z = 0 and z = 1/z j where z j are the zeros of f LRK .We deform the contour out to infinity, and the contour gets snagged on the branch cut and at the poles outside the unit circle.The contribution from these poles will decay exponentially as Θ(z n j n k−1 ) for some k ∈ N, corresponding to the degree of the pole.
Using that Li α (e −x ) is real for x > 0 [50, Eq. 25.12.11], and continuous across the branch cut, the integral along the branch cut is: Using the expansion (D4), and that all other contributions are exponentially decaying, we have that: where z 0 is the zero of f LRK inside of and closest to the unit circle, k ∈ N, and to evaluate the branch cut contribution we also use the expansion [50, Eq. 25.12.12]: We hence see that in general, the edge mode decays as n − min(α,β) .For the case where J = ∆ and α = β, the branch cut integral vanishes, and we have exponentially localised modes, with localisation length ξ −1 = − log(|z 0 |).We can take further terms in the expansion of the denominator in (E3) to derive subdominant terms in the asymptotic expansion, noting that they decay as inverse powers of n and so we have a nice expansion for these Fourier coefficients.With different justifications, two closely related integrals to (E2) were analysed in Refs.[33,36], leading to the same conclusion: a single edge mode decays as n − min(α,β) .This result was moreover in agreement with previous numerical results [28,29,31].
Our method not only gives an analytic approach to finding the asymptotic decay of the single edge mode in this model, it also allows us to predict the decay of the edge modes for higher winding numbers, as discussed in the main text (although we note that finite-size effects can hybridise the edge modes and so the energy eigenbasis in a finite chain may not have the same form as the basis given in Theorem 2).

Finite-size splittings
Let us consider a case of the long-range Kitaev chain (E1) when α = β and J = ∆, and consider |a| = |2J/µ| < ζ(α) −1 , we are then in the trivial phase.We will use the rigorous methods explained in the next section (Rigorous Underpinnings for Conjecture 1) to calculate the relevant Toeplitz determinant that we believe gives us the edge mode splittings.After an overall renormalisation: In this case, we have that 1/f 0 (z) has no negative Fourier coefficients.We can then use Theorem 6 (given below) to see that the splitting is exactly zero if we take f (z) → z ω f 0 (z) for ω > 0. This is trivial to observe at the Hamiltonian level, we have decoupled Majorana modes for this choice of f (z).More interesting is that for ω < 0 the same theorem This is a weak form of Conjecture 1, for the case ω = 1 (and for ω = −1 with the usual replacements), since we conjecture the splitting is Θ(L −ν ), and the proposition is restricted to (α > 2)-decay.The result applies also for |ω| > 1, giving an even weaker form of Conjecture 1 in such cases, since we conjecture the splitting is Θ(L −ν ) where ν ≥ ν, where the inequality is strict if any singularity is filled more than once.
To prove this result, we split the half-infinite edge mode into a mode supported on a region A, consisting of the sites 0 up to L, and the region B consisting of the remaining sites.
We choose the normalisation Z L so that γ 2 A = 1.In the large L limit, Z L tends to a constant; we suppress this from the notation below.Similarly we split the Hamiltonian H BDI = i n≥0,m≥0 t m−n γn γ m on a half-infinite chain as: Since γ A has support only on γ n for n ∈ A, and |ψ has no correlations between A and B, this reduces to: Now: Since we assume that the t n are (α > 2)-decaying, this last expression is summable, and so we have that ε = O(L −ν ).
c. Some results on shifted determinants Our method for calculating the edge mode splittings requires the asymptotics of D N [z ω f 0 (z)].These determinants are related to the functions l(z) and m(z); analysis and a general result can be found in Ref. [55].Roughly speaking, for ω = 1 the edge mode splitting decays like m N , while for ω = 2 the product of edge mode splittings behaves like m 2 N − m N −1 m N +1 , this has some cancellations and behaves like a discrete derivative, leading to the singularity-filling picture.In general this statement holds only up to some error terms that depend on the analytic properties of f 0 (z).
where M (N ) is an |ω| × |ω| matrix with matrix elements: Functions in the class C β have n = β continuous derivatives and the nth derivative satisfies a Hölder condition where β = n + β 0 .As we will see below, the asymptotics of det(M (N )) will decay faster as |ω| increases.Hence, Theorem 5 is limited when looking at large values of ω, where det(M (N )) can be of the same order as the unspecified error term (this motivates the ω max in Conjecture 1).We can evaluate the asymptotics of det(M (N )) using Proposition 2 (corresponding to singularity-filling, see below) if we have an appropriate asymptotic expansion for l(z) or m(z).
For the models considered in the main text, we have that 1/f (1/z) ∈ C β−1 for β = Ω min − ε where ε > 0. Using Proposition 2 below, we have that the decay of det(M (N )) is at most N −γ where ωΩ min ≤ γ ≤ ωΩ min + ω(ω − 1).We are then justified in using the singularity filling picture [73] as long as γ < 3(β − 1) < 3(Ω min − 1).This inequality is violated for ω = 3, while it is satisfied for ω = 2 as long as Ω min > 5. (Note: we also need a nice expansion for m n (or l n ), and as proved above we can use the nice expansion for 1/f (1/z) to infer this up to the first subleading term whenever Ω min > 3.) Note that this is a conservative estimate for the applicability of Conjecture 1, since it is based on the possibility of the subleading term in Theorem 5 becoming relevant.Numerics such as Figure 1 in the main text and Figure 2 indicate that, in those models, the singularity filling continues to apply for higher winding numbers.
The following result from [43] is useful in certain special cases, including when we have f (z) depending only on z: Theorem 6 (Boettcher and Silbermann) Suppose f 0 (z) satisfies the conditions f 0 (z) = 0 on the unit circle, has winding number zero and has absolutely convergent Fourier series [59,74].Furthermore, suppose that the nth Fourier coefficient of 1/f 0 (z) is zero for n < −n 0 ≤ 0. Then for N ≥ n 0 : Now suppose that the nth Fourier coefficient of 1/f 0 (z) is zero for n > n 0 ≥ 0. Then for N ≥ n 0 : This exact formula on the right-hand-side means we can analyse the asymptotics without limits on the winding.We do this analysis in the next subsection.The conclusion is that: Proposition 2 Suppose that l(z) has an asymptotic expansion of the form: then for ω < 0: where c k1,...km > 0, and the sum over k 1 , . . .k m is over all choices where this minimum is achieved.If the minimum is unique then we are guaranteed that this is the dominant term for all N , otherwise the sum may contain cancellations.
The proof relies on a truncation of (F24), so we can also consider cases where we have a nice expansion of l n up to some power.Note that an identical proposition can be written for ω > 0 and where the parameters correspond to the asymptotic expansion of m N (note that these will in general be different to the parameters corresponding to l N ).Using this proposition, we can evaluate the asymptotics of Theorem 6.The decay of this determinant N −e is consistent with the singularity-filling picture [indeed, one can think of the formula (F26) for e as singularity-filling, it is in this sense that we say singularity-filling is rigorous for certain Toeplitz determinants; however our Conjecture 1 also supposes that we can identify the individual eigenvalues that go to zero, going beyond the determinant].The proof of this proposition follows closely parts of the proof of Widom's theorem that we turn to now.

d. Widom's theorem
In Ref. [56], the asymptotics of Toeplitz determinants of z ω f 0 (z) where f 0 (z) is continuous and piecewise C ∞ but not C ∞ are analysed.This means that there are finitely many points (singularities) z h = e iθ h , and at each such point there is a finite α h ∈ N where the α h th derivative is discontinuous (and for all integers k < α h the derivatives are continuous).This is a different definition of singularity compared to the one used in the main text (based on certain asymptotic expansions), but is related, and indeed the same picture emerges.
Theorem 7 (Widom) Suppose we have f (z) = z ω f 0 (z), where f 0 (z) is continous, non-zero and piecewise C ∞ but not C ∞ on the unit circle, and has winding number zero.Suppose that f 0 (z) has m singularities at z 1 = e iθ1 , . . ., z m = e iθm with corresponding α 1 , . . ., α m .Then: Defining Ω θ h = α h + 1, the formula for e is consistent with the singularity-filling picture as discussed in the main text; and in the case the minimum is unique will give the leading term (with non-zero coefficient for all N ) as a rigorous result for the Toeplitz determinant.

Determinants of asymptotic expansions-Proof of Proposition 2
In this section we find the asymptotics of Toeplitz determinants D n [z −N l(z)] based on the asymptotic expansion of l N .The key ideas are all based on Widom's proof of Theorem 7. Let us first recall a lemma [75]: Lemma 1 (Widom 1990) Suppose we have a finite set of measures dµ n (s, t) and functions ϕ n (s), ψ n (t) such that ϕ n (s) j ψ n (t) k ∈ L 1 (dµ n (s, t)).Define the matrix M by: dµ n l (s l , t l ), (F30) where the sum is over all r-tuples (n 0 , . . .n r−1 ).As before we rescale t j → t j /N .Note that for n j = n k we have: . (F49) where we define x to be the number of pairs (j, k) such that j > k and n j = n k , and the integral c({n l }) is positive.Now, let us find the dominant term(s) among these choices of {n j }.Recall that n 0 , . . .n r−1 will correspond to r choices of the m singularities.Let us suppose that the the singularity corresponding to α n has filling k n , then overall we have n k n = r.We also have that the number of (j, k) such that n j = n k = n is given by k n (k n − 1)/2.Hence, 2x = n k n (k n − 1), and moreover: Note that we can write c({n l }) = c k1,...km since we can order the t k as we like.As in Widom's result, it is possible that the different dominant contributions could cancel, so that the dominant asymptotic term is not N −e .We have the correct asymptotics if the minimum is unique, for example.This analysis for l(z) could be repeated identically for m(z), now the behaviour will depend on α n that characterise the singularities of m(z).If either l(z) or m(z) is analytic on an annulus containing the unit circle, then we should include in the expansion those terms coming from the nearest singularity to the unit circle.Finally, putting Ω θn = α n +1 we have Proposition 2, and this agrees with the singularity-filling picture of Conjecture 1.

Appendix C: Analysis of edge modes 1 .
Proof of Theorem 1

Theorem 4
Consider a model corresponding to f (z) = z ω b + (z)b − (z) with open boundary conditions, and suppose that the Fourier coefficients of b − (1/z) −1 have an expansion

n
to cancel the leading term, and get a new expansion g n−m .A linear combination of g n−m and g (1)

r s=1 e
inks n −Ω ks p≥0 a s,p n −p .Then we can use the discussion in the previous section to see b − (1/z) has an expansion r s=1 e inks n −Ω ks Ωmin −2 p=0 a s,p n −p + O(n − Ωmin −1 ) ; thus ν = Ω min + Ω min − 2, recovering Theorem 2 of the main text.
(b)We illustrate this for ω = 4, where we show the numerically-obtained splittings for system size L. Their power-law decays ∼ 1/L ν i are accurately predicted by the 'singularity-filling' of Conjecture 1.For ω > 0

)
Notice that H AA is the Hamiltonian for a finite chain of size L + 1 with open boundaries.Let |ψ be the ground state for the half-infinite chain with Hamiltonian H AA + ∞ n=L+1 c † n c n .The state γ A |ψ is orthogonal to the ground state, and we can consider the variational energy in this parity sector relative to the ground state: