Topology of critical chiral phases: multiband insulators and superconductors

Recent works have proved the existence of symmetry-protected edge states in certain one-dimensional topological band insulators and superconductors at the gap-closing points which define quantum phase transitions between two topologically nontrivial phases. We show how this picture generalizes to multiband critical models belonging to any of the chiral symmetry classes AIII, BDI, or CII of noninteracting fermions in one dimension.


I. INTRODUCTION
The presence of topological edge states, decoupled from the bulk, is a key characteristic of symmetryprotected topological phases of quantum many-body systems 1 . In one dimension (1D), these states are exponentially localized at the physical boundaries of the system, making their energy vanish identically in the thermodynamic limit (and sometimes also at certain finetuned points in the phase diagram of a finite system). With no interactions present, the possible 1D fermionic topological phases are those of the topological band insulators and mean-field superconductors 2,3 , classified by the "ten-fold way" [4][5][6] . For this class of models the very existence of edge states is a consequence of the topologically nontrivial phase structure of the single-particle bulk states ("bulk-boundary correspondence" 7 ), with their robustness against local perturbations (or uncorrelated disorder) being ensured by the symmetries enforced on the perturbations. Well-known examples include the fractionalized soliton mode of the Su-Schrieffer-Heeger model 8 and the Majorana zero-energy mode of the Kitaev chain 9 .
The existence of edge states has conventionally been thought to require that perturbations do not close the insulating band (or quasiparticle) gap. This assumption was proven wrong in 1D in a work by Verresen, Jones, and Pollmann 10 . These authors showed that exponentially localized edge states may survive at quantum criticality, at the gap-closing quantum phase transition (QPT) between two topologically nontrivial gapped phases, i.e. phases which both support topological edge states. Earlier works [11][12][13][14][15][16][17][18][19][20][21][22][23][24] exploring specific models in 1D had anticipated that distinct topological phases − supporting robust edge states − may in fact form also at quantum critical points (and possibly also in higher dimensions in the presence of additional gapped degrees of freedom 25,26 ). However, why and how this happens was first unveiled in detail by Verresen et al. 10 , providing important intuition.
Their theory is underpinned by a study of QPTs within the BDI symmetry class of the ten-fold way, where the single-particle Hamiltonians H exhibit spinless timereversal symmetry, T HT −1 = H, particle-hole symmetry CHC −1 = −H, and chiral symmetry SHS −1 = −H, where T, C, and S are the corresponding first-quantized sym-metry operators 27 with T 2 = C 2 = 1. The critical point separating the two topologically nontrivial phases is labeled by two numbers: a topological invariant ν, which, if positive, counts the edge states, and the central charge c of the conformal field theory (CFT) which describes the scaling limit of the unperturbed Hamiltonian, with no gapped degrees of freedom present 28 . The key to the analysis is provided by a meromorphic function which encodes the properties of BDI Hamiltonians, the zeros and poles of which control the values of ν and c 10 . Given this, Verresen et al. argue that two critical Hamiltonians in the BDI class can be smoothly connected (by tuning a control parameter) only if they share the same values of ν and c. In a follow-up work it is shown that ν and c can be encoded also by correlation functions of certain nonlocal string order parameters 29 . The theory was subsequently put into a larger framework of "symmetryenriched quantum criticality" 30,31 where the presence of nonlocal symmetry operators implies localized and topologically robust edge modes. Also, a generalization to all symmetry classes where the topological classification is larger than Z 2 , independent of dimension, was put forward by Verresen in Ref. 32.
The case when a lattice system is not strictly translational invariant (with respect to the underlying lattice), but instead has a repeating enlarged unit cell, is briefly touched upon in Ref. 33, but then only for the BDI symmetry class, and with no proof of the existence of the critical edge states. In this paper we wish to add to the picture by addressing the unit-cell problem for all three chiral symmetry classes of noninteracting fermions in 1D, i.e. the symmetry classes in 1D with a topological classification beyond Z 2 : BDI, CII, and AIII, containing models subject only to perturbations which respect chiral invariance. Note that other symmetry classes in 1D exhibit at most a single topologically nontrivial gapped phase, and therefore, in keeping with Refs. 10,32, do not support critical edge states. One should here recall that an enlarged unit cell implies that the spectrum of the model − be it a band insulator or a mean-field superconductor − displays a multiband structure in the Brillouin zone. Aside from the possible relevance for experiments − including studies of multiband topological nanowires [34][35][36][37] , quasi-1D fermionic gases in synthetic gauge fields 38 , and 1D topological quantum phase tran-sitions out of equilibrium 39,40 − an analysis of the multiband problem introduces several new facets which may advance our general understanding of topology at quantum criticality. One may here mention that the extended scenario presented in Ref. 32 builds on a low-energy representation of an appropriate Hamiltonian and is therefore not directly applicable to multiband systems where higher-energy bands may impact the topological classification. This provides yet another motivation for a more thorough study of the multiband problem, here narrowed to critical chiral phases in 1D.
The paper is organized as follows. In the next section we address the problem of 1D gapless chiral phases in symmetry class AIII. This is the most general case supporting such phases in 1D since only chiral invariance is enforced on the allowed perturbations; time-reversal and particle-hole symmetries, if at all present, are considered as accidental symmetries. As a warm-up, in Sec. II.A we construct the topological invariant for the simple twoband case, borrowing some of the machinery from Ref.
10, but slightly adapted so as to be immediately extendable to the multiband case, as shown in Sec. II.B. In Secs. II.C and II.D the critical edge states are constructed explicitly for the two-band and multiband case respectively. While our approach for the two-band case in Sec. II.C again closely follows Ref. 10, the approach in Sec. II.D is new. Section III contains a test of our approach in Sec. II when applied to symmetry class BDI, here treated as a special case of AIII, where, in addition to chiral symmetry, time-reversal and particle-hole symmetries are also enforced. This section aims to thoroughly compare our results to those in Ref. 10 where class BDI is treated explicitly. In Sec. IV we turn to the CII symmetry class, with the same symmetries enforced as for BDI but with T 2 = C 2 = −1, thus comprising spinful free fermion models protected by all three symmetries; chiral, timereversal, and particle-hole symmetries. For transparency − and also to chisel out the similarities and differences between the treatment of critical BDI models in Ref. 10 − we here focus on the case of a CII spinful Majorana chain with four bands. Section V contains a numerical check of the robustness of critical edge states against uncorrelated disorder, employing the spinful Majorana chain from the previous section as benchmark model. Section VI, finally, briefly summarizes our work.

II. SYMMETRY CLASS AIII
A. Topological invariant for two-band gapless AIII systems in 1D A noninteracting fermion system with first-quantized single-particle Hamiltonian H is said to possess chiral symmetry if there exists a unitary operator S such that SHS −1 = −H. For a lattice system one can then define two effective sublattices associated with this symmetry by using projectors P A = (1 + S) 2 and P B = (1 − S) 2 where P A B projects states onto the sublattices labelled by A and B, respectively. Note that by using this definition, A and B act indeed as sublattices: the terms present in H cannot couple states from the same sublattice due to the restriction imposed by SHS −1 = −H. Now, let us consider a two-band fermionic lattice system from the AIII symmetry class, defined by requiring that all perturbations, and also the Hamiltonian, respect chiral symmetry. By this, any AIII Hamiltonian H can connect only sites from different sublattices, implying the generic expression where j (n) runs from 1 to N (from −N to N ), with N the number of lattice sites. The two sublattices are denoted by A and B, and t n are the corresponding hopping amplitudes, here allowed to be complex but restricted to a finite range, i.e. t n = 0 for large enough n . Note that A and B denote states which correspond to a pair of internal degrees of freedom. This can be spin, the two sites of a unit cell, or Nambu degrees of freedom. The Hamiltonian in Eq. (1) is easily diagonalized by a Fourier transformation: with k = 0, 2π N, ..., 2π(N − 1) N . Or simply with f (k) = ∑ n t n exp (ikn) in the basis spanned by A, k⟩, B, k⟩. Importantly, the function f (k) is seen to be in one-to-one correspondence with H(k). It follows from Eq. (3) that the eigenenergies are given by k = ± f (k) (which means that f (k) = k e iϕ k for some ϕ k ) and therefore the zeros of k coincide with the zeros of f (k). If nondegenerate, such a zero, k = k 0 , implies that k ∼ k − k 0 , bringing about a massless relativistic excitation of the critical theory (with the bulk energy gap being closed). When the internal states A and B in Eq. (1) are Nambu degrees of freedom, allowing for a nontrivial Majorana representation of the Hamiltonian in Eq. (1), the excitations are those of Majoranas (as in Ref. 10), else they are ordinary fermions. Provided that all zeros of f (k) are nondegenerate, precluding the appearance of dispersions k ∼ (k − k 0 ) m with m ≠ 1 a dynamical exponent, one infers that the effective field theory which describes the critical phase is that of a conformal field theory 10 . For gapped systems, with k ≠ 0, f (k) = k e iϕ k is a well-defined function on the unit circle with the phase factor e iϕ k prescribing a mapping S 1 → S 1 . By this, the winding number ν which defines the topological invariant for the 1D AIII symmetry class can be identified with the number of times that f (k) winds around the origin in the complex plane as k is swept through the Brillouin zone This standard formula for a winding number topological invariant breaks down for gapless AIII phases, since now k vanishes for at least one value of k. However, as realized by Verresen et al. for the BDI symmetry class 10 , one may circumvent this difficulty by performing an analytic continuation of f (k) to the entire complex plane, i.e., taking f (k) → f (z) in such a way that The winding number ν can now be calculated using Cauchy's argument principle as with N z the number of zeros of f (z) inside the unit disk, z < 1, and where N p is the number of poles (counting multiplicities) also inside the unit disk. Importantly, the right-hand side of (6) remains well defined also in the gapless case. The only difference to the BDI case is that t n are now allowed to be complex-valued constants since time-reversal invariance is not enforced for AIII. The quantity on the right-hand side of Eq. (6)  We next consider the extension to the case with 2n distinct states per unit cell, with n > 1. The corresponding single-particle Hamiltonian H(k) in k-space will now be represented by a 2n × 2n matrix, implying 2n bands in the Brillouin zone.
To set the stage, let us choose a basis for the unit cell that diagonalizes the chiral operator S, and call the corresponding basis states X, i⟩ (with eigenvalue +1) and Y, i⟩ (with eigenvalue −1) where the labels X, Y run over the 2n different states within a cell and i runs over the cells. We can always go to this basis by a unitary transformation. To satisfy the chiral symmetry condition SHS −1 = −H we can have couplings only between states of opposite eigenvalues, in other words, X, i⟩ can couple only to Y, j⟩. Therefore, the most general single-particle multiband Hamiltonian can be written as with i, j running over all cell indices. , finally, are hopping amplitudes, which, as for the two-band case, are allowed to be complex.
We now diagonalize H by performing a Fourier transformation, writing with k = 0, 2π N, ..., 2π(N −1) N . Reading off the Hamiltonian in k-space: with the n × n matrix F (k) composed of elements f XY (k) = ∑ j t XY j exp (ikj). For example, with four bands we have Clearly, the eigenenergies are zero if and only if det F (k) is zero. Therefore, the zeros of H(k) coincide with the zeros of d(k) = det F (k), with F (k) in one-to-one correspondence with H(k). The winding number ν for a multiband gapped system is now calculated as For the same reason as for the two-band case, this formula becomes inapplicable for gapless phases. To get around this we again use analytic continuation, d(k) → d(z), with where T j is the corresponding hopping matrix constructed out of the t XY j coefficients. By employing Cauchy's argument principle, we obtain the same expression for the winding number ν as in the two-band case, ν = N z − N p , now with N z the number of zeros and N p the number of poles (including multiplicity) of d(z) inside the unit disk. By construction, this expression for the winding number is valid also for a gapless multiband system in the AIII symmetry class and hence can be used to label its distinct critical phases.

Extending the unit cell
Any lattice system invariant under a translation by a unit cell has an equivalent description in which the unit cell has been enlarged, but now with more internal degrees of freedom. In particular, chiral symmetry is preserved under an extension of the unit cell since the new Hamiltonian, written in a basis with the enlarged unit cell, still only couples sites belonging to distinct sublattices. We can enlarge the unit cell until we end up with hopping only between nearest-neighbor cells, i.e. only the matrices T −1 , T 0 , and T 1 in Eq. (12) are nonzero (with the row-and column-indices of these hopping matrices now running over all internal states in the extended unit cell). This drastically simplifies the analysis since we are left with only three terms in the summation over j in Eq. (12). As a result, the winding number ν = N z − N p can now be calculated using The poles in this expression appear only in the prefactor z −n and it follows that we can write ν =N z − n whereN z is the number of zeros inside the unit disk of detF (z), whereF Note that the counting of zeros include possible zeros at the origin that will cancel one or several of the poles of z −n in Eq. (13). This is a useful result that we shall exploit when proving the existence of critical edge states for a multiband system. Also important to note here that the enlarging of the unit cell does not change the generalized winding number ν = N z − N p 33 . In passing, let us stress that the overbars in the formulas above and in the rest of the paper are used to denote different variables, not complex conjugation.

Short comment
Any chiral critical noninteracting fermion system in 1D − i.e. any system exhibiting a critical phase and belonging to symmetry class AIII, BDI, or CII − can be classified using the approach above. While we have here exploited only the presence of chiral invariance, assuming symmetry class AIII, the additional symmetries enforced by CII and BDI will only add restrictions on F (z), and hence on d(z) in Eq. (13). It may also be worth emphasizing that we have not assumed any specific form for the chiral symmetry operator S, but only used that S is unitary and that hence there exists a basis which diagonalizes S.

C. Edge states in two-band gapless AIII systems in 1D
We now turn to the demonstration of AIII critical edge states. Our plan of attack for the case of two bands is very similar to that of Verresen et al. for the BDI symmetry class 10 (with substantial changes when we later turn to multiband systems). The goal is to construct ν linearly independent states (per edge) "by hand", with the properties that (i) their energies vanish identically for a semi-infinite chain, and (ii) their wave functions decay exponentially as one moves away from the edge, by this establishing a bulk-boundary correspondence 2 for 1D critical two-band systems in the AIII symmetry class.
For clarity, let us again write the Hamiltonian in Eq.
(1), but now explicitly for a semi-infinite chain: Here ±Λ are cutoffs beyond which the range of hopping vanishes. The state B, i + n⟩ is a null vector for i + n ≤ 0. Introducing a state this state will represent a zero mode of H if with the sum over unit cells constrained by i+n≥1. This gives us the following conditions on the coefficients which multiply the A, i⟩ and B, i⟩ states in Eq. (17), call them C A,i and D B,i respectively: with i ≥ 1. The idea is now to use the zeros of f (z), Eq. (5), to prove that the coefficients in Eqs. (18a) and (18b) can be chosen so as to yield precisely ν edge states ψ α ⟩ at edge edge of the chain, α = 1, ..., ν, with the properties (i) and (ii) above. We denote by z α the largest ν zeros within the unit disk, with the rest of the zeros denotedz s : f (z α ) = 0 with α = 1, ..., ν, and f (z s ) = 0 with s = 1, ..., N p .
Case N p = 0 : Let us first consider the case when there are no poles in f (z). This means that t n<0 = 0 in Eq. (5) since the corresponding terms are the ones that create poles in f (z). In this case Eqs. (18a) and (18b) contain all hopping amplitudes t n present in f (z). It is therefore easy to construct the zero modes by taking b with the second line obtained by summing the geometric series, given that z α < 1. Note that we here use a notation where (z α ) 0 = 1 also for z α = 0. The ν states obtained by inserting b .., ν, into Eq. (16) have zero energy by construction, and moreover, they decay with the unit cell index i since z α is inside the unit disk, implying an exponential decay ∼ exp(i ξ α ) with localization length ξ α = −1 ln z α . Thus, both conditions (i) and (ii) above are satisfied. Let us note in passing that when a zero z α approaches the unit circle, the localization length is seen to diverge, signaling criticality, with the corresponding edge mode hybridizing with the bulk spectrum, leaving room for a massless bulk excitation at z α = 1, in accord with the discussion in Sec. II.A.
All cases with distinct z α , α = 1, ..., N z , yield N z linear independent edge states. But what if there are m degenerate zeros of f (z) (meaning that for α 1 , ..., α m the zeros are z α1 = ... = z αm = z α )? Following Ref. 10 we can here take and one verifies that this implies thus producing m linear independent edge states as required. Summarizing, when N p = 0 we can readily construct N z edge states with the required properties (i) and (ii), valid also when the system is critical. Note that the edge states thus obtained are nonzero only on sublattice B since we have chosen a (α) m = 0. One should here note that given a basis in which the chiral symmetry operator S is diagonal, zero-energy edge states of any chiral-symmetric model necessarily have support on only one sublattice of a bipartite lattice 41 .
Case N p ≠ 0 : In this case we look for b s } complex numbers to be determined. As before we take all a (α) m = 0. We know that t n<−Np = 0 since otherwise the multiplicity of the pole at the origin would be larger than N p . It follows from Eq. (18a) that the condition C A,i>Np = 0 is trivially satisfied. It remains to find a set of constants λ (α) s that will make the rest of the expansion coefficients C A,1≤i≤Np vanish identically. Again invoking Eq. (18a), these constants must satisfy the equations where s and summation over s is implied. These equations are equivalent to a matrix equation of the type Aλ = b, with A an N p × N p matrix with elements A is . As follows from the discussion in Ref. 33, provided that z α is nondegenerate, a row reduction turns A into an invertible Vandermonde matrix, implying the existence of unique solutions for the λ (α) s . If there are degenerate zeros of f (z) one can employ the same device as in the case with degeneracies when N p = 0, working with derivatives instead; cf. Eq. (20). However, we shall not go into details here.
Turning to the case of negative winding numbers, we cannot directly employ the elegant inversion argument used in Ref. 10 since in our case, with complex hopping amplitudes allowed in the AIII symmetry class, inversion symmetry is in general broken. However, the necessary modification of the inversion argument is minor, and essentially comes down to doing some "relabeling" in the preceding equations for positive winding numbers.
We start off with the same Hamiltonian H as before, Eq. (15), and interchange A with B, and n with −n. We thus obtain with t ′ n = (t −n ) * and with A, i−n⟩ a null vector for i−n ≤ 0. Note that this is just a rewriting of the Hamiltonian and hence the edge states of H and H ′ must be the same. Let us focus on H ′ , in one-to-one correspondence with the number of zeros on (outside) the unit circle. We can then write with, as before, N P the multiplicity of the pole of f (z) at the origin. It follows that the multiplicity of the pole of f ′ (z) is N ′ P = N − N P . Moreover, by the inversion z → 1 z, f ′ (z) must have the same number of zeros inside the unit circle as f (z) has outside the unit circle, that is, N ′ Z = N O . The winding number ν ′ for the inverted system is now easily calculated as Given this result we can construct the edge states for ν < 0 in exact analogy to the case ν > 0, using the zeros , the condition that this state is a zero mode of H ′ takes the form with the sum over unit cells constrained by i − n ≥ 1. It follows that the coefficients C ′ A,i and C ′ B,i which multiply the A, i⟩ and B, i⟩ states respectively must satisfy These are the same equations as for the case with ν > 0, Eqs. (18a) and (18b), and hence we can construct the edge states in exactly the same manner, but now using the zeros of f ′ (z) that will give us ν ′ = −ν −N C edge states. It may be worth pointing out that these states live on sublattice A, reflecting the fact that we interchanged A and B when rewriting the Hamiltonian, replacing H by H ′ .

D. Edge states in 2n-band gapless AIII systems in 1D
We now proceed to generalize our finding from the previous section and establish a bulk-boundary correspondence for 1D critical multiband systems in the AIII symmetry class, i.e. systems with 2n bands (with n ≥ 1, treating two bands as a special case). As for the simple two-band case, our objective is to construct ν linearly independent states per edge with the properties that (i) their energies vanish identically for a semi-infinite chain, and (ii) their wave functions decay exponentially as one moves away from the edge. Here we shall take the unit cell large enough so that only the hopping matrices T 0 , T ±1 in Eq. (12) are nonzero. This is a convenient choice which simplifies the analysis, and which does not impact the result for the number of edge states.
For transparency, let us begin by explicitly writing down the Hamiltonian in Eq. (7) on a semi-infinite lat-tice, having chosen a sufficiently large unit cell so that there are only three sets of nonzero hopping amplitudes t XY j , with j = 0, ±1: where Y, 0⟩ is a null vector. As spelled out in Sec. II.B, the labels X and Y run over the 2n different states within a unit cell, X ∈ {A, B, ...} on one sublattice and Y ∈ {C, D, ...} on the other, with X and Y implicitly summed over in Eq. (28). From now on, all occurrences of repeated indices X and Y are summed over. Next, we write a general expression for a multiband state, This state is a zero mode of H if Similarly to the two-band case in Sec. II.C, Eqs. (18a) and (18b), this gives us the following constraints on the coefficients which multiply the X, i⟩ and Y, i⟩ states: with i ≥ 1. Using matrix notation, where T m−i is the hopping matrix constructed out of the amplitudes t XY i (and similarly for T * i−m ), and where a Case T −1 = 0: For this case, with only T 0 and T 1 being nonzero, it is straightforward to construct the zero modes. Recall from Sec. II.B that the winding number ν can be calculated from is the number of zeros (poles) of det F (z). When T −1 = 0 there are no poles, and the expression for the winding number simplifies to ν = N z . Let us take a more pragmatic route here and simplify the computations by assuming that every zero of det F (z) is nondegenerate. This is motivated by the fact that since a degeneracy will split unless the hopping amplitudes are fine-tuned, any generic experimental uncertainty will wash away the degeneracies, making the degenerate case irrelevant for applications.
Any zero of det F (z), denoted by z α with α = 1, ..., N z , guarantees that there exists a nonzero eigenvectorb Choosing a (α) m = 0 in Eq. (32b) one infers the existence of ν zero-energy states, living on one of the sublattices.
In exact analogy to the two-band case, the states obtained via the construction above − inserting b − decay exponentially as exp(i ξ α ) with localization length ξ α = −1 ln z α , z α < 1. Thus, also condition (ii) above is satisfied.
Case T −1 ≠ 0: Let us denote by z β (with z β < 1) and byb (β) all solutions tō Assuming that all zeros are nondegenerate (cf. the discussion above when T −1 = 0), their total number isN z .
that are all linearly-independent among each other because they decay at distinct rates z m−1 β . Given this, we now aim to construct independent edge states out of b (β) m . Thus, we look at superpositions b (α) m = ∑N z β=1 λ αβ b (β) m and notice that for each such superposition the coefficients in Eq. (32a) vanish identically. The existence of zeroenergy states then hinges on the possibility to construct independent vectors b (α) m such that with some coefficients x δγ . It follows that for each δ we can then construct a state from b (α) m which trivially fulfills the condition C 1 = 0. Importantly, the states thus obtained are linearly independent since they are each tied to a particular component b (δ) m , with each b (δ) m having a unique decay rate z m−1 δ . Following the protocol above, the vectors b (α) m are seen to define at least ν =N z −n = N z −N p independent edge states. This concludes the analysis for this case.
We approach this case in the same way as when we analyzed negative winding numbers for two bands in Sec. II.B, essentially relying only on a "relabelling" of the equations derived for positive winding numbers.
By interchanging the sublattice indices X and Y and flipping the sign of the index j which labels the hopping amplitudes t XY j in Eq. (28), the multiband Hamiltonian H in the same equation takes the form with (t ′ j ) XY = (t XY −j ) * and where X, 0⟩ is a null state. Recall that repeated indices X and Y are always summed over. The function F (z) = T −1 z −1 + T 0 + T 1 z, associated to H, now gets replaced by a new function F ′ (z), with the property that det F ′ (z) = (det F (1 z * )) * . By exactly the same argument as for the two-band case, it follows that if det F (z) calls for ν < 0 then det F ′ (z) entails ν ′ = −ν −N C with N C the number of zeros of det F (z) on the unit circle.
The construction of the edge states is simple, and can be carried out with the case ν > 0 as a template. To see how, we act with H ′ on the multiband state ψ α ⟩ in Eq. (29), reading off the condition that ψ α ⟩ is a zero-energy state: where a 0 = 0. This equation implies the following constraints on the coefficients (for i ≥ 1): m,X = 0 (35b) or, using matrix notation, These are the same equations as (32a) and (32b) for the case with ν > 0, but with T replaced by T ′ . Thus, we can construct the zero-energy edge states in exactly the same manner as for ν > 0, but now using the zeros of det F ′ (z) that will give us ν ′ = −ν − N C such states.

A. Two-band gapless BDI systems in 1D: topological invariant and edge states
The problem of critical edge states in the BDI symmetry class of 1D models was solved by Verresen et al. 10 for the case of two bands. Since any model which belongs to the BDI symmetry class is also contained in AIII − the symmetry class which we analyzed in the previous section, including the multiband case − one may think that there is not much to add. Still, it is instructive to revisit the problem to see how the BDI formalism in Ref. 10 fits into the scheme presented here, taking into account the symmetries which mark out the BDI class. Beyond chiral invariance, intrinsic also to AIII, these are time-reversal symmetry, T HT −1 = H, and particle-hole symmetry, CHC −1 = −H, with T 2 = C 2 = 1. To avoid any misconception, it is important to stress that these symmetries are to be considered as accidental when a Hamiltonian H is placed in the AIII symmetry class, while being enforced on any perturbation of the same Hamiltonian when placed in BDI. In particular, this implies that the symmetry-protection of edge states is stronger if H is placed in AIII rather than in BDI.
Verresen et al. 10 performed their analysis on a representative BDI Hamiltonian, H BDI , describing a spinless superconducting wire. In real space, assuming translational invariance on a chain with unit cells labeled by j, the second-quantized Hamiltonian can be written as H BDI = ∑ n t n H n , where Here γ j = 1 2 (c † j + c j ) andγ j = i 2 (c † j − c j ) are Majorana modes, with c j and c † j fermion operators. Note that timereversal symmetry precludes terms of the form iγ j γ and iγ jγ as well as complex hopping amplitudes in the decomposition of H BDI . As implied by the notation, each unit cell with index j contains two Majorana modes, γ j andγ j .
Particle-hole symmetry is built into H BdG , with C = τ x K (where τ x is the Pauli x-matrix and K the complexconjugation operator), as is also time-reversal symmetry with T = K, with both operators here written in a k-space representation. Here T 2 = C 2 = 1 as must be since the model is spinless. The chiral symmetry operator S is given by τ x in the chosen basis. We can make a rotation to a "chiral basis" where S is diagonal by simply acting with U = (1 + τ y ) √ 2 on the basis states. As a result, τ z → τ x and τ y → τ y and therefore H BdG → − 1 4 ∑ j,n t n (τ x ⊗ j⟩⟨j + n + iτ y ⊗ j⟩⟨j + n ) + H.c., now with chiral symmetry operator S = τ z . Rewriting H BdG in terms of the eigenstates of τ z , call them A⟩ and B⟩, one obtains H BdG = − 1 2 ∑ j,n t n A, j⟩⟨B, j + n ) + H.c., of the very same form as our Eq. (1) (up to an immaterial prefactor of −1 2). It follows that the f (z) function derived for the BDI symmetry class in Ref. 10 is given by the same expression as in our Eq. (5), with the crucial difference that the parameters t n are now constrained to be real. According to the fundamental theorem of algebra, this restricts the zeros of f (z) to be real or to come in complex-conjugate pairs. Apart from this, the overall picture does not change. As worked out by Verresen et al. 10 , the number of topological edge states is determined by the zeros and poles inside the unit disk of the f (z) function, just as for AIII. presents a formal proof that the topological invariant for a critical translational invariant system in the BDI symmetry class is independent of the choice of unit cell. While intuitively clear that this must be the case, the proof primarily serves as a consistency check of the construction of an f (z) function when the unit cell is enlarged. If we restrict the generic AIII multiband Hamiltonian in Eq. (1) to BDI by requiring invariance under time-reversal and particle-hole symmetry where T = 1 ⊗ 1K and C = τ x ⊗ 1K in a k-space representation our d(z) function in Eq. (12), now restricted by these symmetries, can be shown to be identical to that of Ref. 33, denoted by f (z) in Eq. (S8) of that same reference. Note that the unit matrix multiplying K in the expressions for T and C acts in the n-dimensional space of internal states contained in the unit cell (not counting the Nambu degree of freedom with one particle-and one hole-state). The necessary calculation is straightforward but long-winded, and for this reason we omit it here. Moreover, the result is entirely expected: Our construction of critical edge states for multiband AIII models ensures the existence of exponentially localized zero-energy edge states also for critical multiband models in the BDI symmetry class, with the same counting as before, ν = N z − N p .

IV. SYMMETRY CLASS CII
A model which is placed in symmetry class CII is protected by all three symmetries S, T, and C of the tenfold way 1 , analogous to BDI, but now with T 2 = C 2 = −1. Thus, this is the appropriate symmetry class to use for chiral-invariant spinful fermion systems with topological edge states protected by time-reversal and particle-hole symmetry. We address this problem similar to how we approached the BDI class in the previous section, where we analyzed how the functions f (z) (for two bands, Eq. (5)) and d(z) (for the general multiband case, Eq. (12)) get restricted by the additional symmetries T and C. Here, we have to proceed with some care.
First, we cannot satisfy the constraints T 2 = C 2 = −1 for a two-band model. The reason is simple. For these constraints to be satisfied, the symmetry operators T and C, both being anti-unitary, have to be represented by a product of a Hermitian (and unitary) matrix with purely imaginary elements and the complex conjugation operator K. There is only one such 2 × 2 matrix, the Pauli y-matrix, and hence, with T and C being distinct operators, there is no faithful representation adequate for two bands.
Secondly, turning to more than two bands, one still has to be wary: There are multiple choices for the representation of the symmetries and they may influence the d(z) function in Eq. (12) differently. Here, for ease and clarity, we limit our discussion to the simplest multiband case, 2n = 4, with conventional realizations of the timereversal and particle-hole symmetries; see below.

A. Topological invariant for four-band gapless CII systems in 1D
For the purpose of connecting the representative fourband CII models to possible physical realizations (for examples, see e.g. Refs. 42,43), we shall begin by identifying the possible varieties of spinful superconducting chains that satisfy the required symmetries. The case of four bands is easy to picture, thinking in terms of a Bogoliubov-de-Gennes construction using Nambu spinors as in our discussion of the BDI symmetry class, but now with an added spin-1/2 degree of freedom which doubles the number of bands. We take the time-reversal symmetry operator to be of conventional type for spin-1/2 systems, T = 1 ⊗ σ y K (in a k-space representation), with the first (second) factor acting in the particle-hole (spin) space, and with K acting in both, yielding T 2 = −1. The particle-hole symmetry is decoupled from the spin degree of freedom and therefore its symmetry operator has to be C = τ y ⊗ 1K (in the same representation) in order to fulfill the requirement C 2 = −1. Here τ y is also a Pauli y-matrix, but acting in the particle-hole space (cf. the text after Eq. (38)). The chiral symmetry operator S is fixed by S = T −1 C. In the following we shall explore how the restrictions implied by the T and C symmetries map out the possible spinful chains belonging to the CII class. With an eye to parallels with the BDI symmetry class as studied in Ref. 10, we will work in a representation with Majorana modes. Different from Ref. 10, the Majoranas will now carry spin-1/2.

Spinful Majorana chains
To construct the possible CII Majorana chains we must analyze what types of spinful Majorana bilinears that survive the restrictions from time-reversal and particlehole symmetry. Our strategy is to first construct all real-space first-quantized Hamiltonian matrices H allowed by these symmetries and then use a Nambu representation to go to second-quantization, and from there, to a representation in spinful Majoranas, defined by , with c † j,σ and c j,σ fermion operators at lattice site j with spin σ =↑, ↓.
As an outcome we obtain the spinful Majorana chains with serving as bases for all spin-preserving and spin-flipping Majorana chains, respectively. The amplitudes t n and t n are real, with [−Λ, Λ] the range of couplings between the Majoranas. The associated hopping matrix in spin space, constructed as in Sec. II.B, is readily read off from Eqs. (39) - (41): For details, see the Appendix. By going to k-space via a Fourier transform and performing an analytic continuation to the entire complex plane, one can now use Eq. (42) to define an F -matrix by F (z) = ∑ n T n z n , in exact analogy to how it was done in Sec. II.B for the AIII symmetry class. We can again use the construction d(z) = det F (z) for extracting the winding number ν = N z − N p , with d(z) well defined also at criticality when the gap is closed (and d(z) has one or several zeros on the unit circle).
Any first-quantized Hamiltonian which can be put in symmetry class CII can also be put in AIII (however, with a stronger protection of its topological edge states since now there are many more types of perturbations against which the states are protected, having removed the constraints of time-reversal and particle-hole symmetry). One thus expects that the bulk-boundary correspondence that we derived in Sec. II for AIII should still be valid. In other words, one anticipates that the winding number ν = N z − N p obtained for a 1D model in symmetry class CII correctly counts its number of topological edge states, also at criticality. In fact, when discussing the multiband problem for symmetry class BDI in Sec III we relied precisely on this line of argument, grounded in the work by Verresen et al. 10 on two-band models. While the argument is expected to be valid also for CII, the absence of results for two-band models in CII − since none exists in this symmetry class! − may call for a closer look. To this we turn next. However, rather than proving the existence of edge states by an explicit construction for CII, we shall again exploit the bulk-boundary correspondence for critical edge states in AIII. To make the argument formally sound we will make the connection between the CII and AIII winding numbers mathematically manifest, taking off from the four-band Majorana chain in Eq. (41). This will clarify why the CII winding numbers can take only even integer values, whereas there is no such restriction for AIII 1 . Let us here point out that the generalization to 2n bands with n > 2 becomes cumbersome to handle in the second-quantized Majorana representation. Instead, it is preferable to start directly with a first-quantized single-particle CII Hamiltonian and proceed as for the multiband BDI systems in Sec. III.B, exploiting the results in Sec. II.D for multiband singleparticle models in the AIII symmetry class. Again, the calculation is straightforward but long-winded. Since the existence of edge states is fully expected given the analysis of the four-band case, we omit it here.

Connection between CII and AIII winding numbers
We begin by writing the general CII Hamiltonian in Eq. (41) on first-quantized form, reversing the Nambu procedure from above and by this extracting the corresponding BdG Hamiltonian H CII BdG , call it simply H: We then perform a sequence of unitary transformations acting in the particle-hole and spin spaces, with the transformations given by The chiral symmetry operator S becomes diagonal in the rotated basis and takes the form S = ∑ j τ z ⊗ 1 j⟩⟨j (as is easily verified by checking the identity SH ′ S −1 = −H ′ ). Moreover, H ′ is diagonal in spin space and can be decomposed as H ′ = j,n (t n + it n ) A, ↑, n⟩⟨B, ↑, j + n + (t n − it n ) A, ↓, j⟩⟨B, ↓, j + n + H.c., (45) where A, σ, j⟩ and B, σ, j⟩ are particle and hole states, respectively with spin σ =↑, ↓, attached to the unit cell with index j. But this is nothing but two copies, labeled by ↑ and ↓, of the general two-band AIII Hamiltonian in Eq. (1)! It follows immediately that the winding number, is derivable from d(z) = det F (z) = det ∑ n T ′ n z n , with We can explicitly rewrite d(z) = det F (z) = f 2 (z) + g 2 (z) with f (z) = ∑ n t n z n and g(z) = ∑ ntn z n . This result is in perfect agreement with the one obtained above using a Majorana representation, but now we can directly refer to the bulk-edge correspondence proved for class AIII. This establishes that the winding number ν = N z − N p for symmetry class CII in 1D correctly counts the number of topological edge states, also at criticality. As mandated by the tenfold way for the 1D CII symmetry class 1 , ν is restricted to even integers: Factorizing d(z) = (f (z) + ig(z))(f (z) − ig(z)), complex zeros are seen to come in conjugate pairs, with real zeros coming with even multiplicities since t n andt n are real. (The poles trivially come with even multiplicities.)

B. Edge states in four-band gapless CII systems in 1D
Let us finally address the character of the zero-energy edge states, showing that they are Majoranas.
Zero-energy states have support on only one sublattice in a basis where the chiral symmetry operator S is represented by a diagonal matrix. This transpired from our analysis in Sec II, and holds quite general 41 . Knowing that S is indeed diagonal in the rotated spin-diagonal basis of H ′ , we can therefore consider, without loss of generality, a zero-energy state, call it φ⟩, with support on the A sublattice only. When writing Eq. (45) we associated A with a particle state, (1 0) T [with sublattice B being associated with a hole state, (0 1) T ]. In this notation, and leaving out the spatial part of φ⟩, we thus write φ⟩ = (1 0) T ⊗ (a b) T = (a b 0 0) T , with (a b) T a spinor in spin space with amplitudes a and b. Particle-hole symmetry implies that C φ⟩ is also a zero-energy state, where, in the rotated basis, C = τ z ⊗σ x K. We can always construct two new zero-energy states by adding and subtracting φ and C φ⟩, in this way obtaining φ⟩ ± = (u ± 1 u ± 2 0 0) T with u ± 1 = a ± b * and u ± 2 = b ± a * . We now go back to the original basis by carrying out the inverse transformations A acting in the composite particle-hole and spin space in any unit cell on the A sublattice), we then express the ψ⟩ ± states in second-quantization as with the composite particle-hole and spin space spanned by c † ↑ 0⟩ = (1000) T , c † ↓ 0⟩ = (0100) T , c ↑ 0⟩ = (0010) T , and c ↓ 0⟩ = (0001) T . Reading off from Eq. (47), one immediately identifies the mode operators which correspond to these states: Using that u ± 1 = ±(u ± 2 ) * , it follows that d ± = ±id † ± , implying that the zero-energy states thus constructed are Majorana states.

V. ROBUSTNESS OF THE EDGE STATES: A NUMERICAL TEST
In the previous sections we showed that the boundary states in unperturbed critical 1D multiband chains are connected to a topological number ν, generalizing the ordinary winding number for gapped systems to critical systems in chiral symmetry classes AIII, BDI, and CII. In gapped systems the edge states are robust to any perturbation from uncorrelated disorder which respects the relevant symmetries and leaves the bulk gap open. Here we numerically verify on a case study that the edge states of a CII critical chain also exhibit robustness to disorder although the gap is now closed. By a symmetry analysis we identify the chiral symmetry as solely responsible for the protection of the states. A possible mechanism of such protection was proposed in Ref. 32 for the case of two-band models, but its generalization to the multiband case has remained somewhat unclear and requires a more thorough investigation.
For a numerical test we take a topologically nontriv- rewritten in first-quantization as with N + 1⟩ and N + 2⟩ null states. Having placed H {1,2} in symmetry class CII, its relevant symmetries are identified by the time-reversal operator T = ∑ N j=1 1 ⊗ σ y K ⊗ j⟩⟨j , the particle-hole operator C = ∑ N j=1 τ y ⊗1 K⊗ j⟩⟨j , and the chiral symmetry operator S = ∑ N j=1 τ y ⊗σ y ⊗ j⟩⟨j , where K is the complex-conjugation operator.
In analogy to gapped symmetry-protected topological systems we are interested to see if the edge states of H {1,2} in Eq. (49) remain localized under symmetrypreserving locally uncorrelated perturbations and, if so, if they stay pinned exactly at zero energy. The localization of the edge states can be quantitatively addressed by calculating the participation ratio (PR) that is defined by R = 1 ∑ i p 2 i , where p i is the occupation of the Bogoliubov quasiparticle at site i, explicitly obtained by summing up all probability amplitudes of the corresponding Nambu spinor. The PR quantifies localization of an eigenstate: A completely localized state has R = 1 while a completely delocalized state, such as a plane wave, has R = N (where N measures the extent of the chain, count- ing the total number of unit cells). The robustness of an edge state is established if both of the following conditions are satisfied: its energy level remains at zero and the PR is of order unity. In Fig. 1 we have perturbed the edge states with random uncorrelated on-site disorder proportional to τ x ⊗ 1 (T and C preserving), τ y ⊗ σ y (C breaking), and τ y ⊗ 1 (T breaking). The edge states display robustness when both symmetries are preserved (implying that also chiral symmetry, S = T −1 C, is preserved) but get destroyed once one of the symmetries, T or C, is broken (implying that also chiral symmetry gets broken 1 ). Thus, the critical edge states are seen to be well protected in the CII symmetry class.
It is interesting to inquire what happens if allowing disorder which breaks both T and C symmetries but leaves chiral symmetry unbroken. This is tantamount to put the Hamiltonian H {1,2} in the AIII symmetry class where any perturbation which respects only chiral symmetry is allowed. One may maybe object and say that this assignment of symmetry class is impossible: "Particle-hole symmetry is a 'built-in' feature of any second-quantized Hamiltonian expressed in a Nambu spinor basis, implying that the corresponding first-quantized BdG Hamiltonian (like H {1,2} ), as well as any perturbation thereof, is also particle-hole symmetric." However, whereas the use of the Nambu basis does preclude particle-hole symmetry-breaking perturbations in the second-quantized theory, there is no such constraint on perturbations of the BdG Hamiltonian. In other words, once the BdG Hamiltonian matrix has been extracted from the underlying second-quantized theory using a Nambu basis, this matrix defines a single-particle theory which can be put into the AIII symmetry class, allowing for perturbations that break particle-hole symmetry. True, such perturbations do not represent physical symmetry-breaking perturbations of the original unperturbed theory of a mean-field superconductor since they cannot be traced back to a second quantized formulation using the Nambu basis. However, the maneuver allows us to formally address the question whether critical edge states of a first-quantized Hamiltonian (like H {1,2} ) are robust against perturbations which respect only chiral symmetry. This is an important issue, independent of whether such perturbations can be realized in an underlying second-quantized theory or not.
With this as a backdrop, we now consider perturbations which break both T and C but, different from the cases of Fig. 1, preserve S. The results are displayed in Fig. 2 where we have applied random uncorrelated onsite disorder proportional to τ x ⊗ σ y and τ z ⊗ σ y (T and C breaking, S preserving), showing that the critical edge states of H {1,2} do survive perturbations which respect only chiral symmetry.
Given these results we conjecture that it is precisely the chiral symmetry S which protects the topological edge states in disordered 1D critical phases. Numerical examinations of other disorder types and critical chains support this conjecture 44 . One should here note that while chiral symmetry by definition is indeed the only possible protecting symmetry for critical phases of AIII, it is a priori not evident that it actually fulfills this role. However, our numerical results provide strong evidence that it does.

VI. SUMMARY
Building on the work by Verresen et al. 10 on critical two-band BDI models in 1D, we have carried out a study of critical multiband models in any of the 1D chiral symmetry classes AIII, BDI, and CII. We use an approach where the enlarged unit cell responsible for the multiband structure of a model is further extended until one is left with intercell hopping of fermions only between nearest-neighbor auxiliary cells. This allows for a transparent and rigorous analysis of the problem, enabling us to prove the existence of critical edge states for any 1D multiband model belonging to one of the chiral symmetry classes. As in the original work in Ref. 10, the number of such edge states is coded by a topological invariant generalizing the winding number of gapped 1D models, now providing a bulk-boundary correspondence for all chiral critical phases in 1D.
A numerical case study of a four-band model in the CII symmetry class − perturbing its Hamiltonian by uncorrelated disorder distributions of different symmetry contents − shows that the robustness of its critical edge states is protected solely by chiral invariance, with timereversal and particle-hole symmetries playing no role. We conjecture that this picture is general, with chiral symmetry being the sole protecting symmetry of 1D critical topological edge states not only in the AIII symmetry class (where any perturbation which respects only chiral symmetry is allowed), but also for any 1D model belonging to the CII or BDI symmetry class. Put differently, we expect that the subsets of the CII and BDI symmetry classes composed of models that support topological critical phases are entirely contained within the AIII symmetry class, demoting T and C to accidental symmetries when it comes to protection of these phases.
Our work may open a path towards a more comprehensive study of symmetry protection of multiband topological phases at criticality, including those of nontranslational invariant models in higher dimensions and artificially generated phases from Floquet topological engineering 45 . The classification of periodically driven (Floquet) systems at criticality is here a particularly promising direction for exploration: Our multiband analysis should be quite useful for treating critical timeperiodic systems within Floquet theory -a theory that is intrinsically multiband due to the repetition of frequency zones 46,47 . In fact, the Floquet systems represented within frequency domain and truncated at sufficiently large frequency index are mathematically equivalent to time-independent multiband systems [46][47][48] . Any Floquet chiral system is then anticipated to satisfy the time-independent chiral relation after the truncation, and combined with our results this shows the existence of the topological edge states. Importantly, these arguments are expected to hold also at the closed induced gap (known as the 'anomalous gap') corresponding to the existence of anomalous edge states 48,49 , having no analog in static systems, also at criticality.
Another important topic to explore is the effect of interactions in multiband critical models. There are already some results 10,30,32,33 on the robustness of critical topological edge modes against interactions, using density matrix renormalization group (DMRG) and effective field theory methods, but so far only for two bands. The extension to more bands is technically challenging, but we expect that our present work will be of use also here.

APPENDIX: CII spinful Majorana chains
To obtain all possible CII spinful Majorana chains we first construct all real-space first-quantized Hamiltonian matrices H allowed by the proper time-reversal and particle-hole symmetries and then use a Nambu representation to go to second quantization with complex fermions and, from there, to a representation in spinful Majoranas.
Recall that the symmetry conditions are T HT −1 = H and CHC −1 = −H, with T 2 = C 2 = −1. In a real-space representation (with H defined on a chain with N unit cells) we have that T = N j=1 1⊗σ y K ⊗ j⟩⟨j , C = N j=1 τ y ⊗1 K ⊗ j⟩⟨j , (A1) with j running over all unit cells, and where τ α and σ α are Pauli matrices when α = x, y, z and equal to the 2×2 unit matrix 1 when α = 0. The operator K effects complex conjugation. For the four-band case, with the real-space Hamiltonian given by a Hermitian 4N × 4N matrix, we make the observation that any such matrix can be written as a linear combination over real numbers of matrices H j,j+n = (i) m τ α ⊗ σ β ⊗ j⟩⟨j +n + H.c., m = 0, 1, (A2) with j constrained by j + n ≤ N . Examining all 32 such matrices (for any fixed j and n), we find that 8 of them respect the T and C symmetries. Writing out only the parts which act in the spin-and particle-hole spaces (suppressing the common j⟩⟨j+n factors and the Hermitian conjugate to save space), these are: To construct the symmetry-respecting secondquantized Hamiltonians H j,j+n corresponding to (a) -(h), we introduce four-component real-space Nambu spinors for the N unit cells: Ψ † j = (c † j,↑ c † j,↓ c j,↓ − c j,↑ ) and Ψ j = (c j,↑ c j,↓ c † j,↓ − c † j,↑ ) T , j = 1, 2, ..., N , yielding the 4N -component Nambu spinors for the full lattice, Ψ † = (Ψ † 1 ...Ψ † j ....Ψ † j+n ...Ψ † N ) and Ψ = (Ψ 1 ...Ψ j ....Ψ j+n ...Ψ N ) T respectively. Inserting the expressions for H j,j+n implied by (a) -(h), cf. Eq. (A2),