Topological properties of the long-range Kitaev chain with Aubry-Andr\'e-Harper modulation

We present a detailed study of the topological properties of the Kitaev chain with long-range pairing terms and in the presence of an Aubry-Andr\'e-Harper on-site potential. Specifically, we consider algebraically decaying superconducting pairing amplitudes; the exponent of this decay is found to determine a critical pairing strength, below which the chain remains topologically trivial. Above the critical pairing, topological edge modes are observed in the central gap. For sufficiently fast decay of the pairing, these modes are identified as Majorana zero-modes. However, if the pairing term decays slowly, the modes become massive Dirac modes. Interestingly, these massive modes still exhibit a true level crossing at zero energy, which points towards an initimate relation to Majorana physics. We also observe a clear lack of bulk-boundary correspondence in the long-range system, where bulk topological invariants remain constant, while dramatic changes appear in the behavior at the edge of the system. In addition to the central gap around zero energy, the Aubry-Andr\'e-Harper potential also leads to other energy gaps at non-zero energy. As for the analogous short-range model, the edge modes in these gaps can be characterized through a 2D Chern invariant. However, in contrast to the short-range model, this topological invariant does not correspond to the number of edge mode crossings anymore. This provides another example for the weakening of the bulk-boundary correspondence occurring in this model. Finally, we discuss possible realizations of the model with ultracold atoms and condensed matter systems.


INTRODUCTION
Topological phases of matter have now been established as an intriguing research area in condensed matter physics [1][2][3][4][5]. Non-interacting symmetry protected topological (SPT) phases are either band insulators or superconductors: they have a gap in their energy spectrum of the bulk and are characterized by a global order parameter, called topological invariant. These topological invariants are integers. As a consequence, they cannot be smoothly deformed, within the same symmetry class of the Hamiltonian, to another phase labelled by a different value of the invariant. Therefore, topological phase transitions are always accompanied by the closing and subsequent reopening of the bulk gap. A classification of SPT phases in terms of their symmetries (time-reversal, particle-hole and chiral) and their dimensionality has been performed in the celebrated "periodic table of topological insulators and superconductors" [6][7][8][9][10]. Each class is characterized by a different topological invariant which can take integer Z, even-integer 2Z or binary Z 2 values. Finally, SPT phases possess gapless boundary states protected against local perturbations and the number of these protected surface states is proportional to the topological invariant.
The superconducting Kitaev chain [11] is a prime example of a 1D non-interacting topological system. The main constituents of the p-wave Kitaev chain [11] are spinless fermions that undergo short-range hoppings, primarily between nearest neighbors, in a 1D lattice. They can also get paired together into superconducting Cooper pairs with opposite momenta. The Kitaev chain exhibits time-reversal, chiral and particle-hole symmetry; thus it belongs to the BDI class of the topological classification [6]. This model can host non-local Majorana modes (MMs), i.e. zero-energy modes localized at the two boundaries of the chain. Such system is also characterized by a topological invariant, called winding number, and the number of MMs at each end of the chain is equal to this winding number. Therefore, these MMs are topologically protected against local perturbations and, hence, cannot be removed without a topological phase transition, i.e. a transition from a topological superconducting phase to the trivial superconducting phase when all MMs become paired.
The robustness of the MMs makes them promising candidates to build topological quantum computers. They could be used as qubits that can store and manipulate quantum information in a decoherence free manner [11]. Recently, there have been many proposals for the realisation of MMs in different systems like heterostructures of topological insulators and s-wave superconductors [12] or cold fermion systems near p-wave Feshbach resonance [13], or with attractive s-wave interaction in the presence of Rashba spin-orbit coupling and Zeeman field [14,15]. Other proposals include heterostructures of spin-orbit-coupled semiconductor thin films [16,17], or nanowires [17][18][19] coupled via proximity effect with s-wave superconductors and a Zeeman field. Majorana fermions may also emerge from magnetic nanoparticles on a superconductor without spin-orbit coupling [20,21].
The short-range Kitaev chain is a version of the inte-grable p-wave superconducting chain of fermions. The primary focus of the present work is a long-range Kitaev chain, in which the superconducting pairing term decays with distance as a power-law [22][23][24][25][26]. While the topological phases of the short-range Kitaev chain are characterized by Z topological invariants, the present model posseses a new unconventional topological phase characterized by half-integer Z/2 winding numbers, even though the system belongs to the same BDI class and respects the same discrete symmetries of the short-range model. For slow power law decay, i.e. for decay exponent α given by α 1, the MMs of the short-range model (α 1) coalesce to form massive nonlocal edge states called massive Dirac modes (MDM). These new edge states lie within the bulk gap and are topologically robust against local perturbations that do not violate fermionic parity and particle-hole symmetry. Thus, like the MMs, MDMs may also find novel applications in the area of topological quantum computations.
One must note that the quantum critical behavior of systems with long-range interactions/pairings can be very different from that of the short-range ones, due to the non-local nature of the interactions. Therefore, over the recent years, much effort has been devoted to understanding the static [27,28] and dynamic properties [29] of long-range systems through quantum quenches [30][31][32], periodic drives [26], through establishment of the Lieb-Robinson bounds [33], and through studies of the growth of entanglement entropy [34], Bell inequalities [35], and thermalization [36,37] (for a detailed review, see [38]).
Interestingly, although there have been numerous recent claims regarding the detection of Majorana fermions in superconducting/semiconducting heterostructures [39][40][41][42][43][44][45][46], the quest for Majorana fermions has been quite challenging and even controversial. A major issue that needs to be addressed in any system hosting MMs or MDMs is the presence of correlated or uncorrelated local inhomogeneities in the onsite potential, which may alter the physics of such systems.
In recent years, new experimental techniques in the field of ultracold atoms have emerged, where a superposition of two periodic lattices with commensurate or incommensurate periods can be created [47][48][49][50][51][52]. In particular, this allows for the experimental implementation of theoretical models such as the Aubry-André-Harper (AAH) model [53,54] and opens the road for the study of the role of spatial inhomogeneities in close connection to the experiments discussed above. The AAH model raised great interest as, for incommensurate lattices, the system is quasiperiodic and exhibits a localization-delocalization transition due to its self-dual nature [53,[55][56][57]. More recently, the AAH model and similar quasiperiodic models have also been extensively studied for their topological properties [58][59][60][61][62]. In fact, the AAH Hamiltonians can be seen as the dimensional reduction of a 2D Hofstadter model [63], which describes electrons in a 2D lattice subjected to a perpendicular magnetic field akin to that of a quantum Hall system. Both models, the AAH chain and the 2D Hofstadter lattice, are described by a set of Harper equations [54], and the two momentum variables of the 2D model map onto one momentum variable and a 2π-peridodic Hamiltonian parameter in the 1D model. Accordingly, the AAH model also exhibits a fractal Hofstadter butterfly energy spectrum, and can host topologically protected boundary modes reminiscent of the 2D system; it is in turn characterized by a 2D topological invariant called Chern number.
The interplay between nearest neighbor p-wave superconducting pairing and the AAH potential has been studied in Refs. [64][65][66], for both the commensurate and the incommensurate cases. One of the key observations in such systems has been that a finite amount of superconductivity is required for unpaired Majorana modes to be present in the system [67]. Although comprehensive studies [66,68,69] have been performed on unifying the topological phase diagram for a range of periodic, incommensurate and disordered potentials in the short-range Kitaev chain, the combined effect of long-range pairings and AAH onsite potential on the topological phase diagram of the Kitaev chain has never been explored. Here, we present a comprehensive study of the 1D long-range Kitaev Hamiltonian with AAH type spatially inhomogeneous potential. Through extensive numerical calculations of the topological invariants and edge spectra, we reveal how the critical behavior changes throughout the topological phase diagram. Specifically, we show that, for sufficiently slowly decaying pairing, the central bands are characterized by the same half-integer topological invariant as in the long-range Kitaev chain with homogeneous chemical potential. The edge modes appear at finite energy, which suggests to interpret them as massive Dirac modes (MDMs). However, the presence of the AAH potential induces a non-monotonous dispersion of these modes when tuning the chemical potential or the superconducting pairing strength. This leads to a crossing of these levels in isolated points at zero energy, suggesting an intrinsic relation to Majorana physics. At the same time, none of this interesting behavior at the system edge is reflected by the topological invariant, suggesting that the system exhibits a weakened bulk-boundary correspondence. This conclusion is corroborated by studying the behavior in higher band gaps induced by the AAH potential. Specifically, we find that long-range pairing removes edge state crossings by turning them into avoided crossings, whereas the topological invariants of the corresponding bands remain unchanged. Plan of the paper. Section II introduces the model and discusses its general properties. Section III is devoted to the discussion of the winding number and challenges related to the numerical computation of such winding number in the presence of long-range pairing. In Section IV, we analyze the effects of the AAH modulation on the central gap behavior of the long-range Kitaev chain. This includes characterizing the MMs and MDMs and their corresponding winding numbers, studying the existence of a critical superconducting pairing, and a detailed anal-ysis of the level crossings at zero energy. Section V is devoted to the effect of the long-range pairings on higher energy gaps, which contain edge states arising from the AAH potential. Finally, in Section VI we discuss possible experimental realizations of the model in condensed matter and in systems of ultracold atoms or molecules.

II. MODEL HAMILTONIAN
We consider a chain of spinless fermions with an onsite AAH modulation of the chemical potential and a longrange p-wave superconducting pairing. Its Hamiltonian can be written as where t is the hopping amplitude, µ is the chemical potential, ∆ is the superconducting pairing amplitude and c i (c † i ) are the annihilation (creation) operators at the i-th site of the chain. The superconducting pairing is taken to decay as a power law l −α , where l denotes the distance between the sites and the scaling exponent α ∈ R [70]. For a constant chemical potential f (i) = 1, we recover the long-range Kitaev chain with homogeneous onsite potential. As mentioned above, this model is known to be topologically equivalent to the short-range Kitaev chain with nearest-neighbor pairing terms for α 1. In contrast, for α 1 the model can host a long-range topological phase with MDMs, characterized by a half integer topological invariant.
We also consider an AAH onsite chemical potential where the modulation frequency β gives the periodicity of the potential and the parameter φ shifts the origin of the modulation. In particular, for irrational β, the system becomes incommensurate with respect to the lattice periodicity. The latter can lead to interesting effects such as gap openings or localization-delocalization transitions [53,56,57]. In this paper, we consider a system with commensurate values of β = p/q, and characterize its topology. Both the AAH model and the Kitaev model are interesting due to their topological properties separately from each other. Here, we study the interplay between the two. The AAH potential has periodicity of q sites and we therefore consider a system of N sites with L supercells, where N = qL. First, we write the Hamiltonian in the Bogoliubov-de-Gennes (BdG) basis to properly treat the pairing term. The Hamiltonian of Eq. (1) then reads where T is the basis in real space within a supercell denoted by u. H local stands for the elements of the Hamiltonian involving operators within the supercell. H hop contains the hopping terms between supercells. Finally, H l includes the longrange pairing terms. The explicit form of each contribution can be found in Appendix A. Moreover, we can either impose open boundary conditions (OBC) to study the edge states, or anti-periodic boundary conditions (APBC) assuming χ u+L = −χ u to study the bulk properties of the system (see Appendix B for details).
To focus on the bulk properties of the system, let us impose APBC and compute the Hamiltonian in the Fourier space, which reads where χ k denotes the Fourier transform of the supercell vector, with the Fourier sum being performed over the quasimomenta k = (2m + 1)π/L, 0 ≤ m ≤ L − 1. Finally, for the thermodynamic limit (L → ∞), the Hamiltonian reads as In order to better understand the contribution of the long-range pairing, it is worth looking into the definition of H inf , which takes the form (see App. A for more details on the derivation) where C x,y = ig xy σ y . The function g xy is defined as follows (7) where HLP stands for the Hurwitz-Lerch-Phi function, also called the Lerch transcendent function [22,71]. The HLP function has a singularity at k = 0 for α < 1. This singularity encodes the long-range character of the Hamiltonian and contributes both to the edge dynamics and to the bulk topology of the system, creating a longrange unconventional topological phase. Our goal will be to characterize it.

III. WINDING NUMBER
As mentioned in Section I, the model under study lies in the BDI symmetry class of topological insulators and superconductors, that is, it is particle-hole, time-reversal and chiral symmetric. Therefore, in order to characterize the bulk topology of the model, we compute the winding number considering three methods. The first one is the Fukui-Hatsugai-Suzuki algorithm [72], which, as we will see, fails to characterize the unconventional long-range topological phase. Secondly, we consider the real-space winding number for finite systems, which converges to the expected value of the winding number for large systems. Finally, we implement an algorithm which relies on the chiral symmetry of the Hamiltonian. For this approach, we make use of the analytical expression for the infinite system from Eq. (5).
Before getting into the details of each method, we would like to stop and briefly analyze the long-range case α 1 [22][23][24][25][26]. On the one hand, the long-range character of the Hamiltonian affects the edge states of the system: the MMs hybridize and become massive and form MDMs. On the other hand, the bulk topology of the model hosts two different unconventional topological phases, characterized by a half-integer Z/2 winding numbers. It is important to note that the winding number for short-range systems can only take integer Z values [6][7][8][9][10]. Nevertheless, the long-range character of the Hamiltonian allows for a new type of characterization of the bulk topology. We will see that the bulk topology of our system is well-defined, but not directly connected to the edge dynamics. In fact, the bulk-edge correspondence gets weakened for long-range systems [25].

A. Fukui-Hatsugai-Suzuki algorithm
A common approach to compute the winding number is the algorithm introduced by Fukui, Hatsugai and Suzuki [72]. The main idea is to discretize the Brillouin zone of the system {k 1 , ..., k s }, and then to construct the tensor U (i) m,n = ψ m (k i+1 )|ψ n (k i ) for every k i . Here ψ m(n) is the m(n)-th eigenvector of the Hamiltonian. The winding number ν is then defined as where | · | denotes the determinant of the tensor U . To compute this quantity, we work with periodic boundary conditions. This means that we need to close the loop by computing the scalar product between |ψ m (k s ) and |ψ n (k 1 ) .
Recall that the long-range character of the Hamiltonian is encoded in the HLP function in Eq. (7), and this function is singular at k = 0 for α < 1. That is why, for the case of α < 1, one needs to consider a discretization of the Brillouin zone without an explicit inclusion of k = 0 in order to avoid the singularity. Nevertheless, even if one considers an appropriate discretization of the Brillouin zone, this algorithm always fails to characterize the long-range topological phases through half-integer Z/2 winding numbers. In fact, it can be proven that this algorithm always leads to an integer Z winding number. In order to understand why, let us consider the case of f (i) = cst. This case is straightforward as the model has only two bands. However, this argument can be generalized to the multi-band (non-Abelian) case, as required, for investigating the properties that we are interested in [see Eq. (8)]. The Zak phase of a single band n is given by the sum of the relative phases over the Brillouin zone where Note that the function F 1 (k l ) is periodic along the first Brillouin zone, meaning that F 1 (k l +sδ k ) = F 1 (k l ), where δ k is the spacing of the discretization and s are the number of points considered. Moreover, F 1 takes values within the principal branch of the logarithm, such that −π < F 1 (k l )/i < π. This means that F 1 (k l ) can be written in terms of forward differences and an integer-valued field n 1 If one closes the loop around the Brillouin zone, the finite differences cancel out, and the Chern number is always integer valued The only way to obtain non-integer winding numbers would be either to neglect the contribution of F 1 (k s ), therefore not closing the loop, or to include the point k = 0. Neither option is reliable, since the first one results in gauge dependent quantities and the second one relies on a singularity point.

B. Real-space winding number
A second approach to compute the winding number is based on the real-space Hamiltonian of a finite system with N = qL sites [73]. In this case, the winding number is defined as follows where Tr is the trace per unit cell and The quantity Q is defined as Q = P + − P − , where P + and P − are the projectors on the positive or negative eigenstates respectively P ± = ± |ψ ± ψ ± |. Γ A and Γ B are the projector operators onto the A and B sublattices of the system [74], namely Γ A = diag(1, 0, 1, 0 · · · ) and Γ B = diag(0, 1, 0, 1 · · · ) and X is the position operator.
We emphasize that computing [X, Q AB ] is not straightforward for a system with periodic or anti-periodic boundary conditions as the position operator is not periodic. We here use the approach introduced in Ref. [75], for which we consider a system with a domain −N/2 ≤ x i ≤ N/2, where x i 's are the eigenvalues of the position operator X. Then, the commutator can be written as where λ are the solutions to z 2N +1 = 1 and c λ = λ N +1 1−λ . For this approximation to be effective, the system has to be sufficiently large.

C. Infinite system winding number
Finally, a third approach to compute the winding number relies on the definition of the winding number for the infinite system. Since the system is chiral, one can always find a basis in which the Hamiltonian is block-offdiagonal (see Appendix C) For such Hamiltonian, the chiral operator is diagonal and the winding number can be computed with the help of the expression [74] where h is the upper-right block of the Hamiltonian in Eq. (15). In our case, we are interested in using the discrete version of Eq. (17), which takes the form Finally, note that h contains the HLP function defined in Eq. (7). In order to compute the winding number, we need the analytical derivative of this function We will use the winding number computation explained in this section to study the central gap of the model in the following section.

IV. CENTRAL GAP OF THE KITAEV CHAIN WITH AAH POTENTIAL: A QUANTITATIVE STUDY
The AAH potential splits the energy spectrum of the Kitaev chain into several bands. In this Section, we focus on the behavior of such system around its central gap, that is, around E = 0. Specifically, we compare the behavior of short-ranged and long-ranged Kitaev chains in the presence or absence of the AAH modulation.
In this context, we first study the energy spectra as a function of the chemical potential µ, which reveal the existence of MMs and MDMs for different values of the decay exponent α. Then, we characterize these different phases using the real-space III B and the infinite winding numbers III C. Here, we also compare the accuracy of the two methods. We then investigate how energy spectra and winding numbers depend on the superconducting pairing term ∆, and find that the presence of an AAH modulation leads to a critical superconducting pairing ∆ C : below this critical value, we observe neither MMs nor MDMs in the system. The value of ∆ C depends on the decay exponent α as well as the amplitude of the chemical potential µ. Finally, we discuss a peculiar behavior found in the presence of both long-range pairing and AAH potential: while the long-range nature of the pairing leads to a degeneracy splitting of the MMs and turns them into MDMs, the AAH is responsible for a true level crossing of these modes which allows them to re-combine at isolated points, and which yields an adiabatic connection between MM phase and MDM phases.
A. Effect of the AAH potential on the energy spectrum We first study the energy spectrum as a function of the chemical potential µ, comparing the case of homogeneous chemical potential, f (i) = cst = 1, with the case of a modulated AAH potential. We consider a system of 233 sites, and for the AAH potential, we take the modulation frequency β to be p/q = 144/233. We expect, however, that the results qualitatively hold also for other values of β, if q is sufficiently large. Figures 1 (a) and (c) depict the energy spectrum of the system for homogeneous chemical potential with APBC (blue) and OBC (green). For the short-range system (a), we see the appearance of zero-energy MMs in the central energy gap for |µ/t| < 1. For |µ/t| > 1, the system Energy spectrum for the system with OBC (green) and APBC (blue) for q = 233, L = 1, ∆ = 0.5t and φ = 0. We consider f (i) = 1 in (a) and (c) and AAH chemical potential in (b) and (d). Panels (a) and (b) depict the energy spectrum for the short-range system and thus, exhibit MMs. Panels (c) and (d) depict the energy spectrum for the long-range system, which harbours MDMs in the central gap.
is topologically trivial, and no edge states exist within the bandgap. For the long-range system, shown in Figure 1 (c), we only see one gap closing at µ/t = 1. For µ/t < 1, the system exhibits MDMs, which are clearly separated from the bulk and only exist for OBC. For µ/t > 1, in contrast, the system is topologically trivial and does not exhibit any edge states. Figures 1 (b) and (d) depict the energy spectrum for a system with an AAH potential. The AAH potential leads to a splitting of the spectrum into many bands. Band gaps at non-zero energy and the massive edge states hosted by these gaps will be analyzed in Section V. In the present Section, we focus on the central gap around E = 0. For the short-range system, shown in Figure 1(b), we observe the existence of the zero-energy MMs in the central region around µ ≈ 0. As compared to the case of homogeneous potential, the AAH modulation significantly increases the parameter range supporting MMs, which now extends to |µ/t| 2. For the long-range system, shown in Figure 1(d), the MMs hybridize to form MDMs. In this case, the central region for which MDMs exist is |µ/t| < 1. Interestingly, the MDMs exhibit a crossing at zero energy which will be further discussed in Sec. IV E. Moreover, at |µ/t| ≈ 1.5, we see another opening of weakly gapped regions which do not host MDMs or MMs.
We conclude that the Kitaev chain with AAH potential exhibits a variety of different phases which are distinguished through the presence or absence of MMs and MDMs. In the following, it will be our goal to further characterize these different phases through the winding numbers introduced in Sec. III.
Real-space winding number for L = {1, 3, 5} and infinite system winding number. In both cases, we considered q = 233, ∆ = 0.5t and φ = 0. We see that, when increasing L, the real-space winding number converges to the one for the infinite system. This convergence is slower for long-range systems (c-d). Note that the computation of winding number is not valid for the values of µ for which there is no gap at zero energy (shaded regions).

B. Comparison between real-space and infinite system winding numbers
Before performing an in-depth analysis of the winding numbers and the edge states of the Kitaev chain in the presence of an AAH potential, we compare the two algorithms to compute the winding number introduced in the previous Section. Figure 2 shows the comparison between the real-space winding number for different values of L and the infinite system winding number. For the short-range system (a)-(b), where we consider α = 5, the winding number takes the discrete values 0 and 1. Note that for the short-range case, the FHS algorithm from Section III A would lead to valid results. The non-trivial central region for which ν = 0, which is |µ/t| < 1 for f (i) = cst (a) and |µ/t| < 2 for AAH onsite potential (b) coincides with the existence of MMs in Figures 1(a)-(b). The real-space winding number here is computed for a system with L = {1, 3, 5}, while the infinite system winding number takes L → ∞. For the short-range system, both algorithms show good qualitative agreement even for L = 1. Thus, here N = q = 233 is large enough for the approximation in Eq. (14) to converge.
Figures 2 (c)-(d) depict the winding number for the long-range system with α = 0.5. Here again, the realspace winding number is computed for a system with q = 233 and L = {1, 3, 5}. The winding number takes half-integer values, ±1/2. In the case f (i) = cst, shown in panel (c), the positive winding number corresponds to the region for which µ/t < 1, which is also the region that corresponds to the presence of MDMs [See Figure 1 (c)]. For µ/t > 1, the system does not exhibit any edge states In the latter case, the combined effect of AAH onsite potential and longrange pairing leads to a succession of topologically differentiated regions separated by gap closings, see panel (d). As for the case of constant potential, the winding numbers in presence of AAH potential are integer in the short-ranged system, and half-integer in the long-ranged system. All quantities are computed using the real-space Hamiltonian for a system with q = 233 and L = 1.
but it is neither trivial, since the winding number takes a non-zero value: this system exhibits weak bulk-edge correspondence [25]. For the AAH onsite potential, shown in Figure 2 (d), the positive winding number corresponds to the region where |µ/t| < 1. Again this is the region for which the system presents MDMs in the energy spectrum, cf. Figure 1(d). The two neighboring regions, 1.5 < |µ/t| < 2, show a negative value of the winding number and do not exhibit MDMs. Note that the gap closes for |µ/t| > 2, and correspondingly, the winding number is no longer well defined, as expected.
For the long-range case, the real-space winding number does not agree exactly with the infinite winding number for the values of L that we have considered. The reason for this is that for long-range systems, the approximation in Eq. (14) converges more slowly, which leads to small finite-size effects at L = 1.
C. Critical superconducting pairing in short-ranged and long-ranged systems We now investigate how the presence of MMs and MDMs depends on the value of the superconducting pairing ∆. To this end, we compute the energy spectra of the finite Hamiltonian with OBC as a function of ∆, as well as the real space winding number for the same system size with APBC. The results are depicted in Figure 3, showing energies in blue and winding number in red. Figures 3 (a)-(b) show the results for a system with rapidly decaying pairing (α = 5). Panel (a) considers a system with homogeneous chemical potential (µ/t = 0.5) and panel (b) considers a system with AAH chemical potential (µ/t = 1.2). For the homogeneous potential, the system exhibits topological edge modes even in the limit ∆ → 0. In contrast, the presence of the AAH potential leads to a critical pairing ∆ C > 0 which is required for the opening of the gap and the appearance of MMs. The real space winding number also becomes non-zero only for ∆ > ∆ C . This consistency between number of edge modes and topological invariant is an example of the celebrated bulk-edge correspondence.
Figures 3 (c)-(d) depict energy spectrum and winding number for a slow power law decay (α = 0.5), again distinguishing between homogeneous potential (c), and AAH potential (d). The critical pairing also remains ∆ C = 0 for homogeneous chemical potentials. However, the combination of an AAH potential and long-range pairing is found to yield a very rich scenario: the system is gapless for pairing terms up to ∆ = 0.25t. At this value, a gap opens and closes at ∆ C = 0.55t. There are no edge states within this gap. For ∆ > 0.55t, another gap opens, but now containing MDMs. As already seen in the energy spectrum as a function of µ, the MDMs cross at zero energy at ∆ = 0.85t. A discussion of this interesting behavior will be given in Section IV E. The gapped regimes are characterized by non-zero winding numbers: The winding number takes the value ν = −0.5 in the first gap (the one without edge states), and ν = 0.5 in the second gap (the one with MDMs). Thus, the appearance of MDMs is accompanied by a jump of the winding number by 1, but we emphasize that in general, for α < 1, the system has a weak bulk-boundary correspondence [25]. This means that a non-zero value of the winding number does not always correspond to MDMs. In Appendix D, we show that the bulk-boundary correspondence is weak by looking at the system when varying the superconducting pairing ∆ from negative to positive values. In this case, the winding number also jumps from −0.5 to 0.5, but the edge physics remains the same. Nevertheless, the value of the winding number has physical consequences in the distribution of the Schmidt eigenvalues and the violation of the area law for the von-Neumann entropy [25] or the quantization or half quantization of the multipartite entanglement [76]. Real-space winding number for the system with f (i) = 1 and µ/t ∈ {0.5, 1, 1.5, 2}. We observe two regions, corresponding to α < 1 (long-range) and α > 1 (short-range). At µ/t = 1 (b) there is a closing of the gap, for which the winding number is not well-defined. When we go from µ/t < 1 (a) to µ/t > 1 (c-d) we see a topological phase transition, in which the winding number for short-range goes from 1 to 0 and for long-range from 0.5 to −0.5. Here, ∆C = 0. The realspace winding number is computed for q = 233 and L = 1. We consider that the central gap is closed when ∆E < 0.035t (black regions).

D. Phase diagram
We now study the phase diagram of the model by fixing the chemical potential to a discrete value µ/t ∈ {0.5, 1, 1.5, 2}, while continuously varying α and ∆. The phase diagrams are expressed in terms of the winding number, shown in Figure 4 for homogeneous chemical potentials, and in Figure 5 for AAH potentials.
For the homogeneous chemical potential, Figure 4, we observe that the critical superconducting pairing is always at ∆ C = 0. We see two regions, one corresponding to the short-range Kitaev chain (α > 1), and one corresponding to long-range system (α < 1). In the shortrange case, the winding number is 1 for µ/t < 1, and 0 for µ/t > 1. The long-range scenario has a winding number of 0.5 for µ/t < 1, and −0.5 for µ/t > 1. At µ/t = 1, the system is gapless and the winding number is ill-defined. We also observe in Figure 4 (a) that, for small values of ∆ and for 0.5 < α < 1, the long-range system exhibits MMs, which later hybridizes to form MDMs [24] when ∆ is increased.
For the AAH chemical potential, Figure 5, we make the following observations: For µ/t < 1 (a), the system behaves as in the homogeneous case, exhibiting two different regions, corresponding to long-range and shortrange behavior. From this we can assume that, for very small chemical potential, the AAH modulation does not affect the topology of the system. However, for µ/t ≥ 1, the situation changes, see Real-space winding number for the system with AAH onsite potential and µ/t ∈ {0.5, 1, 1.5, 2}. When µ/t < 1 (a), we only observe two regions, corresponding to α < 1 (longrange) and α > 1 (short-range). Here, the system behaves like in the homogeneous case in Figure 4, and thus ∆C = 0. When µ/t ≥ 1, the critical superconducting pairing ∆C = 0 both for the short-range and the long-range cases. For the short range case there is simply a region where the gap is closed for ∆ < ∆C and a region with winding number 1 for ∆ > ∆C . For the long range case, we see three different regions: first the gap is closed, then it opens with winding number −0.5, then it closes again and it re-opens with winding number 0.5. The real-space winding number is obtained for a system with size q = 233 and L = 1. We consider that the central gap is closed when ∆E < 0.025t (black regions).
∆ C is now shifted to a non-zero value. The latter seems to increase continuously for decreasing values of α. For ∆ < ∆ C , the central gap of the system is closed and the winding number is not defined. If we now focus on the region α < 1, we clearly see a region in which the sign of the winding number changes. This phase does not host MDMs and it is separated from the other phases through a gap closing. Finally, the value of ∆ C depends drastically on the value of the chemical potential amplitude µ.

E. Edge modes in the slow power law decay regime
As observed in the previous subsections, the central gap of the (sufficiently) long-range Kitaev chain with AAH potential hosts an edge mode at positive energy and an edge mode at negative energy, but very interestingly, these modes exhibit crossings at zero energy [see Figure 1 (d) and Figure 3 (d)]. In the present subsection, we focus on this crossing.
First, we note that the chiral symmetry of the system dictates that a crossing of the two modes can only happen at zero energy. However, we are not aware of a symmetry which would prevent the two modes from Energy spectrum (blue) and quantum fidelity F (red) for a system with size q = 233 and L = 1. Here we considered ∆ = 0.5t, α = 0.5 and φ = π. The quantum fidelity F takes the value 1 everywhere except for the closings of the gap at |µ/t| = 1.2 and the crossing at µ/t = −0.85. The quantum fidelity is computed following the MDM indicated in green.
mixing, and thus, they could also perform an avoided crossing. To distinguish a true crossing from an avoided one, we study the quantum fidelity between eigenstates as we adiabatically change the chemical potential. The fidelity is defined as where ψ µ is the single-particle wave function describing the edge mode with positive energy. The fidelity together with the energy spectrum is shown in Figure 6 for a finite size system with OBC and α = 0.5. The positive mode ψ µ is plotted in green, and the crossing happens at µ/t = −0.85. At this crossing point, the quantum fidelity drops to zero. This demonstrates that, at the crossing point, the two orthogonal modes are exchanged, and thus, that the crossing is indeed a true crossing and not an avoided one. Moreover, we also compute the winding number to see if the crossing corresponds to a change in the bulk topology. If one interprets the vicinity of the zero-energy crossing as a reappearance of Majorana modes, one would expect that the winding number becomes 1. However, the crossing is found to have no effect on the winding number, which remains to be 0.5 throughout the whole gapped region (i.e. for |µ/t| < 1.2). On the other hand, as mentioned already in Sec. IV C, long-range interactions can weaken or destroy the bulk-edge correspondence, and therefore, a half-integer winding number may not necessarily exclude the existence of Majorana edge modes.
For a better understanding of the relation between Majorana physics and the observed zero-energy crossings, we shall ask how the zero-energy crossing connects to the regime in which pairing is sufficiently short-range to support an extended phase with zero-energy modes. To this end, we plot in Figure 7 (a-c) the energy splitting ∆E between the two edge modes, the quantum fidelity F of the positive mode, and the winding number ν as a function of both µ and α. For α > 1, the energy splitting becomes very small, and, as we demonstrate below, in this regime the finite value of the splitting is a mere finite-size effect. Therefore, the edge modes in this regime can be interpreted as Majorana zero modes, and such interpretation is corroborated by the winding number which takes the value ν = 1 for α > 1. Nevertheless, the finite-size splitting still allows to energetically differentiate between the two modes, and also for α > 1 we observe a true level crossing, as indicated by the black line of zero fidelity in Figure 7 (b). Interestingly, this curve seamlessly extends to the level crossing at α < 1. In this regime, however, the splitting between the two modes takes significantly larger values [see Figure 7 (a)], and cannot any longer be interpreted as a finite-size effect. This demonstrates that the isolated zero-energy modes at α < 1 are smoothly connected to the MMs at α > 1. In contrast to this, the bulk topology is different in the two regimes: the winding number drops from 1 to 0.5 when α < 1 [see Figure 7 (c)]. We also mention that, for α < 1 and µ/t < −1, the central gap closes [marked in black colors in Figure 7 (c) and yielding the second black line in Figure 7 (b)]. As can be seen from Figure 7(b), the zero-energy crossing of the edge modes (right black line) and the gap closing (left black line) approach each other, when α is reduced. For completeness, we also mention that the gap re-opens at even smaller value of µ. In this regime, no edge modes are observed, yet the winding number takes a non-zero value, ν = −0.5.
To conclude this section, we support our claim regarding the qualitative changes at α = 1 by a finite-size scaling of the degeneracy splitting of the modes, and of the winding number. Therefore, we focus on two values of α = 0.8 and α = 1.1, and fix µ/t = −1 to a value on the left side of the zero-energy crossing. For this choice, we plot the mode splitting ∆E against the inverse of the system size, 1/L, in Figure 8 (a). For the finite-size scaling, we assume an exponential closing of the splitting, ∆E = ae −bL , or a polynomial closing of the splitting, ∆E = a(1/L) b , in both cases with two fit parameters a and b. The exponential fit matches the data neither at α > 1 nor at α < 1. In contrast, the polynomial fit perfectly matches the data for α > 1. However, for the largest system sizes considered (L = 25 supercells), the observed energy splitting at α = 0.8 remains above the polynomial fit. This discrepancy becomes more pronounced and affects more data points for smaller values of α. Thus, we conclude that the energy splitting remains finite for α < 1, while it vanishes polynomially with the system size for α > 1.
The qualitative change at α = 1 is also backed by the behavior of the winding number: While sufficiently above α = 1, it takes the value one, the winding number is 0.5 when α is sufficiently below one. In the vicinity of α ≈ 1, the winding number takes intermediate values between 0.5 and 1, no matter which algorithm we use to compute the winding number. However, we find that in the realspace algorithm the intermediate values are due to the finite size of the system, whereas for the chiral algorithm All the computations consider a system with q = 233, L = 1, and φ = π. The winding number is computed using the real-space approach. Gapless regions (defined by a threshold ∆E < 0.02t), where the winding number becomes ill-defined, are colored black in panel (c). Note that in (b) the region that corresponds to α > 1 is expected to host MMs, but even in this regime, there is a finite-size splitting of the two modes and a true energy crossing at zero energy occurs.
(considering an infinite system) the finite discretization of the momentum grid is the reason for the intermediate values. Scaling the winding number vs. system size or discretization, we observe that the invariant approaches the value one, if α > 1, or the value 1/2, if α < 1. This is illustrated in Figure 8 (b), for the scaling of the winding number vs. discretization D at α = 0.8 and at α = 1.1 This observation suggests a sharp change in the bulk topology at (or in the vicinity of) α = 1.

V. AUBRY-ANDRÉ-HARPER EDGE STATES
In the presence of an AAH potential, the system exhibits topological edge states in the higher energy gaps. These AAH edge states arise from the AAH potential and can be understood through the dimensional reduction of the 2D Hofstadter model. The latter describes electrons in a 2D lattice subjected to a perpendicular magnetic field [63] and its dimensional reduction onto a family of 1D chains corresponds to the AAH model. In particular, different instances of the AAH model are characterized by the 2π-periodic dephasing parameter φ, which can be interpreted as one of the quasi-momenta of the original 2D system, i.e. k x = k and k y = φ.
In Figure 9, we compare the edge states for the system with α < 1 (long-range) and with α > 1 (shortrange) and focus on the energy gap around E = 0.5t. In the normal (short-range) scenario [panel (a)], the edge states cross within the energy gap a given number of times. In the gap considered here, they cross twice (at φ = {1.95, 5.07}). Figure 9(c) shows the energy spectrum around the second crossing and analyzes the quantum fidelity F at the crossing. It drops to zero at φ = 5.07. This must be contrasted to the situation in the long-range system. Figure 9(b) depicts the energy spectrum for the long-range system with OBC, and two avoided crossings are observed at the same values of φ = {1.95, 5.07}. However, the quantum fidelity F in Fig. 9(d) stays at 1 around φ = 5.07. This indicates that the long-range superconducting pairing hybridizes the two edge states, such that the total charge transported in one period of the dephasing now becomes zero.
We now focus on the computation of the 2D Chern numbers. They characterize the bulk topology associated to these edge states, and for a system with bulkedge correspondence, the value of the two-dimensional Chern number should coincide with number of crossings, since it indicates the charge that is pumped at the edges of a 1D cylinder. Here, we use the the Thouless-Kohmoto-Nightingale-Nijs (TKNN) algorithm correspond to the short-range system (α = 5). In (a) we see the energy spectrum for the full period of the dephasing, while in (c) we zoom into one of the avoided crossings. Here, we also compute the quantum fidelity F (dashed red line). Panels (b) and (d) correspond to the long-range scenario (α = 0.5). Again, the upper panel (b) shows the energy spectrum for the full period of the dephasing, and the lower panel (d) zooms into the crossing and also considers the quantum fidelity F (dashed red line). To compute the quantum fidelity, we follow the edge state highlighted in green. We considered a system with OBC and q = 233, L = 1, µ/t = 1 and ∆ = 0.3t. introduced in Ref. [77] where they define the Berry curvature as follows where k x and k y are quasi-momenta of the Brillouin zone. Recall that here, k x = k and k y = φ. We choose the value of the Fermi energy E F to be located inside the gap that we want to characterize. Note that where we consider the analytical derivative of the infinite Hamiltonian in Eq. (5). Then, the Chern number is defined as We here discretize the integral as a Riemann sum where ∆ kx and ∆ ky are the two discretization lengths. We see that at the energy gap around E = 0.5t the 2D Chern number takes the value 2 for any value of the decay exponent α and of the superconducting pairing ∆. This does not agree with the number of crossings within the gap, which we extract from Figures 10 and 11. In particular, in Figure 10 we study the crossing at φ = 5.07. We identify the edge states in the energy spectrum of the system with OBC and study them while changing the parameters of the Hamiltonian. We see that the crossing splits for α < 1 and increasing values of ∆. In Figure 11 (a), we show the energy gap ∆E between the AAH edge states at the crossing point, which clearly increases for α < 1 and for any value of ∆. Accordingly, in Figure 11 (b) the quantum fidelity min(F ) shows a minimum at zero for short-range systems but it does not drop to zero when the system is long-range, indicating that there is no crossing.
We can conclude that the bulk topology arising from the AAH potential does not match the edge dynamics in the presence of the long-range superconducting pairing, and from this we can infer that the weakening of the bulk-boundary correspondence extends, not only to the MDMs but also to the AAH topology.

VI. EXPERIMENTAL REALIZATION
Before concluding we discuss shortly the experimental feasibility of realizing the 1D long-range Kitaev model with Aubry-André-Harper potential. Below we consider condensed matter systems and ultracold atomic systems.
In case (a) under the strong coupling limit, it is theoretically assumed that the magnitude of the proximityinduced superconductivity is the same as that of the host superconductor. In experiments, however, the magnitude of the induced superconductivity is actually found to be much smaller. To take this into account, it is necessary to assume a weak-coupling between the surface system and the host superconductor, where it can be shown that exponentially decaying oscillating long-range pairing and hopping can arise in the chain with a characteristic length scale of the order of the coherence length of the Cooper pairs. Hence, the coherence length is many times larger than the lattice constant, and the long-range pairings are practically non-local [78].
In case (b) of a chain of magnetic impurities placed in a conventional superconductor, each individual adatom can create localized, individual sub-gap Shiba states. These states hybridize with each other and with the bulk superconducting condensate and, through multiple Andreev reflections, form energy bands. This band structure resembles that of the Kitaev chain for a helical spin structure of the adatoms and can be probed by scanning tunneling spectroscopy. In the particular case of deep-lying Shiba states, it has been shown in Ref. [89] that the effective tight-binding Bogoliubov-de Gennes Hamiltonian representing the Kitaev chain has hoppings and pairings that are long-range, with a 1/r-power-law decay within the coherence length of the host superconductor. Beyond this length scale, the decay is exponential.
In both cases, the semiconducting nanowire in a magnetic field, and the nanowire-like chain of magnetic impurities are plagued by defects, either originating from the spatial variance of the magnetic field, or through the distribution of the impurities constituting the chain. This can be well modeled via the Aubrey-Andre model which includes the effect of spatial inhomogeneity. Recently, it has also been proposed that an effective Kitaev chain with long-range pairing and hopping terms can also be formed from the combination of a solid-state Majorana platform made out of a planar Josephson junction in proximity to a 2D electron gas with Rashba spin-orbit coupling placed in an external magnetic field [90].

B. Utracold atoms
We now discuss the building blocks necessary to realize such a system in a cold-atom quantum simulator. On the one hand, realizing a setup that generates the Aubry-Andre modulation is possible. In fact, it has been demonstrated in numerous experiments using super lattice techniques, pioneered in Refs. [91,92]. On the other hand, there has also been some proposals to realize the Kitaev chain in cold-atomic setups [93,94]. The main difficulty for our Hamiltonian is to propose a cold-atomic setup that includes long-range power law superconducting pairing with a controlled decay rate. One possibility would be to engineer the long range pairing with the help of an attractive long-ranged interaction. Here we discuss three classes of cold atom setups that could lead to such interactions: (i)Fermionic atoms/molecules with static dipole moments. Ultracold dipolar gases are in center of interests of the ultracold atoms quantum matter community since the end of 1990s (for reviews/books see Refs. [95][96][97][98]). Static dipoles can be magnetic, as it happens in "magnetic atoms" like Chromium (condensed by T. Pfau group [99]), Erbium (condensed by F. Ferlaino's group [100]), or Dysprossium (condensed by B. Lev's group [101]). Much stronger static electric dipoles can be induced in hetero-nuclear molecules (first achieved in D. Jin's-J.Ye's groups [102]). Strong dipole-dipole interactions lead to non-standard Hubbard models, with many ingredients necessary for realization of the long-range Kitaev chain [103]. The progress in the field of ultracold dipolar gases, and in particular in quantum engineering of chemistry/quantum matter with ultracold molecules is spectacular [104]. For instance, a degenerate Fermi gas of polar molecules [105], or a dipolar quantum gas with metastable supersolid properties [106] were observed recently. The drawback of the static dipole-dipole interac-tions is that they decay as 1/d 3 , that is with α = 3 the system is not very different from the short-range chain.
(ii)Fermionic atoms/molecules with laser-induced dipole moments. Already in 2000, Kurizki and collaborators propose to induce dipole-dipole interactions using off-resonant laser excitation. Such interactions are subject to retardation effects, and as such, in addition to 1/d 3 term, they include 1/d 2 and 1/d 1 terms. O'Dell et al. [107][108][109] showed that particular configurations of intense off-resonant laser beams can give rise to an attractive 1/d interatomic potential. Such a "gravitational-like" interaction leads to stable Bose-Einstein condensates that are self-bound (without an additional trap) with very interesting properties. Light-induced dipole-dipole interactions were intensively studied in experiments [110] and theory [111], leading to more recent observations of long-range one-dimensional gravitational-like interactions in a neutral atomic cold gas [112], or light induced inverse-square law interactions between nanoparticles [113]. The drawback of this approach is that light induced dipole interactions contain not only a conservative part, but a dissipative part too.
(iii) Fermionic atoms/molecules with synthetic Coulomb interactions. Finally, there is a recent proposal for analogue quantum chemistry simulation [114], allowing in particular for synthetic Coulomb interactions between fermionic particles playing role of electrons. This proposal does not belong to the mainstream of quantum computational chemistry [115], which maps fermionic operators to qubits (spins 1/2) via Jordan-Wigner transformations, and uses quantum devices operating on qubits. Still it has generated a lot of interest and excitement in the community, and first experiments in this direction are on their way [116]. The challenge here would be to find a regime with attractive interactions.

VII. CONCLUSIONS
In this paper, we studied in detail the bulk topology and the edge states of the long-range Kitaev chain with an AAH potential. In this context, we discussed several algorithms to characterize the bulk topology, and we derived the phase diagram of the system at half filling. We also studied the appearance of MMs and MDMs in the system: in both short-range and long-range systems with an AAH potential, the appearance of MMs and MDMs requires a non-zero critical pairing. The critical values of the pairing depend on the long-range decay exponent α, and as such these values are very different for the short and the long-range scenario. We also discovered a peculiar behavior: in the long-range case with AAH potential, MDMs exhibit true crossings at zero energy. We found that these crossings have no effect on the winding number. Nevertheless, we showed that they are adibatically connected to the extended phase with MMs for alpha 1, and thus have a direct relation to Majorana physics. Finally, we studied the energy gaps at higher fillings. They are characterized by a 2D Chern invariant, which can efficiently be computed numerically. An interesting behavior appeared in the study of the edge states: the AAH edge modes were dramatically affected by the long-range superconducting pairing, and the bulkboundary correspondence appeared to be weakened.
As an outlook, it would be interesting to study the quasi-periodic limit of the system. Quasiperiodic systems are known to exhibit a multifractal energy spectra [117,118] and critical, extended or localized wavefunctions depending on the chemical potential amplitude [119]. The multifractal and localization properties of the AAH model have been extensively studied for different generalizations of the model [64,66,[119][120][121][122], but the interplay between long-range superconducting pairing and incommensurability has not yet been explored. where: with A j = −µf (j)σ z , B = t 2 σ z −∆iσ y and C l = − ∆ d α l iσ y . Here, we have replaced the l of Eq. (1) by d l = min(l, N − l) to take into account the APBC.
From Eq. (2), we know that the system has a periodicity of q sites, therefore we can rearrange the Hamiltonian into three contributions: H local , which contains the terms within each supercell of q sites, H hop which contains the hoppings between the adjacent supercells and H l , which connects all the supercells of the system through the longrange superconducting pairing where χ u = c qu , c † qu , ..., c qu+(q−1) , c † qu+(q−1) T and L is the number of supercells of the system, namely L = N/q. Each contribution to the Hamiltonian is defined as follows H hop =      0 0 · · · 0 0 0 · · · 0 . . . . . . . . . . . . and, C l,0,0 C l,0,1 · · · C l,0,q−1 C l,1,0 C l,1,1 · · · C l,1,q−1 C l,2,0 C l,2,1 · · · C l,2,q−1 . . . . . . . . . . . . . . . C l,q−2,0 C l,q−2,1 · · · C l,q−2,q−1 C l,q−1,0 C l,q−1,1 · · · C l,q−1,q−1 where B = t 2 σ z and C l,x,y = − ∆ 2d α l,x,y iσ y . For the system with APBC, d l,x,y = min (lq − (x − y), N − (lq − (x − y))) . (A7) Moreover, we explicitly impose APBC by assuming that In order to obtain the momentum space Hamiltonian, we need to use the Fourier transformation of the spinor χ u where χ k = c k,0 , c † −k,0 , ..., c k,q−1 , c † −k,q−1 T . Note that u ∈ {0, ..., L − 1} denotes the supercell index. Combining Eq. (A9) with Eq. (A8), we find which implies that the momentum k is defined as where m ∈ {0, 1, 2, ...L − 1}. The momentum representation of the Hamiltonian takes the form where, in order to simplify the expression, we used If we consider a system with β = p/q, it has a periodicity of q sites. Then, the total length of the system will be N = Lq, where L is the number of supercells. We can take L to be a finite number or consider the system for the limit L → ∞. Either way, we can still consider different types of boundary conditions. Usually, in order to construct the momentum space representation of the Hamiltonian, we need to consider a system with either periodic (PBC) or anti-periodic (APBC) boundary conditions. The long-range pairing terms of the Hamiltonian get canceled if we impose PBC, and that is why we consider APBC when needed in our calculations. It is important to note that, in the infinite limit L → ∞, the boundary conditions can be neglected. Then, there is no technical difference between considering APBC or PBC.

Appendix C: Block-off diagonal infinite Hamiltonian
We transform the infinite Hamiltonian H inf of Eq. (5) as follows where R y = e iσy π 2 ⊗ 1 q . Moreover, we perform a rotation in order to rearrange the basis to the following form χ k = c k,0 , ..., c k,q−1 , c † −k,0 , ..., c † −k,q−1 T .

(C2)
We are left with where H inf is block-off diagonal. Energy spectrum with OBC and real-space winding number for a system with µ/t = 0.5, α = 0.5 and φ = 0. We see that the real-space winding number changes sign at ∆ = 0, while the edge states exist at both sides of the gap closing. The plot is obtained for a system with size of q = 233 and L = 1. Figure 12 shows the comparison between the behavior of the edge dynamics and the bulk topology of the longrange system. In particular, by looking at the energy spectrum for OBC we notice that the system exhibits MDMs for any value of ∆/t ∈ [−1, 1] but ∆ = 0, which corresponds to a gap closing, and the region around ∆ = 0, for which we find MMs. The existence of MMs for small values of ∆ was already adressed in Figure 4. Here, instead, we are interested in comparing the edge dynamics with the bulk topology. If we look at the real-space winding number for the same system, we see that it takes a positive value +0.5 for ∆/t > 0 and a negative value −0.5 for ∆/t < 0. This means that at ∆ = 0 we have a topological phase transition, but this does not affect the edge dynamics. Therefore, we clearly show how the bulk-boundary correspondence is broken when the sign of the superconducting pairing ∆ is changed.