Universal phase diagram of topological superconductors subjected to magnetic flux

We perform a theoretical study of the orbital effect of a magnetic field on a proximity-coupled islands array of px+ipy topological superconductors. To describe the system, we generalize the tightbinding model of the Hofstadter butterfly to include the effect of the superconducting islands. The quantum Hall topological phases, appearing in the absence of superconductivity, are characterized by integer fermionic Chern numbers corresponding to the number of occupied bulk Landau levels. As the strength of the superconducting pairing increases a series of transitions occurs, with one less chiral Majorana edge mode at each consecutive phase, leading to a reduction of the fermionic Chern number by a half. When the pairing potential exceeds the tight-binding model bandwidth, Cooper pairs are localized in the islands, the Chern number is zero, and there are no low-energy edge modes. We identify domains in the model’s parameter space for which the system is topological and supports an odd number of chiral Majorana edge modes. While the precise shape of the domains depends on the details of the model, the general structure of the phase diagram is robust, and it is obtained numerically and in several simplified traceable analytical models. We discuss the relevance of this study to recent experimental research of two-dimensional superconductor arrays on semiconductor systems.


I. INTRODUCTION
Low-dimensional topological superconductors (SCs) have been the subject of vast interest recently 1,2 . That is because they constitue a novel form of quantum matter, and also because they may support Majorana zero modes, which are germane for non-Abelian topological quantum processing 3 .
A time-reversal symmetry breaking topological superconductor is realized when an odd number of energy bands, which are not spin degenerate, are unstable to the creation of Cooper pairs, while the other bands are gapped.
This may be obtained by the proper combination of materials that have strong spin-orbit coupling, are proximity coupled to a superconductor, and break timereversal symmetry by internal magnetic phenomena 4,5 or external application of a magnetic field 6,7 . Several realization schemes for topological superconductivity have been put forward in two spatial dimensions 4,[8][9][10][11] .
Application of an external magnetic field affects electrons both through the orbital degrees of freedom and through the Zeeman term coupled to the spin. For systems with a large g factor the latter is more significant, especially when the magnetic field is parallel to the wire's axis in the one-dimensional case 12 or to the system's plane in the two-dimensional case. However, when the magnetic field is applied perpendicular to a twodimensional system, we expect an interesting interplay between the (quantum) Hall effect [13][14][15] and the emergence of topological superconductivity, decorated by vortices induced by the perpendicular field.
To access this interplay, we study a tight-binding model which gives rise to Hofstadter-butterfly physics 16,17 , with proximity coupling to an array of SC islands (see Fig. 1). Whereas the parts that are not covered by the SC are expected to be susceptible to an external gate potential, sites that are in proximity to the SC islands are less sensitive to the gate potential but are coupled via a pairing potential. This setup, inspired by a recent experiment 18 , enables proximity effects together with tunability of the chemical potential.
The phase of the SC order parameter is assumed to be constant inside each island. This is a good approximation as long as the flux per plaquette Φ is much smaller than the flux quantum Φ 0 = h/e, since in this regime the SC phase does not wind inside an island. To establish the regime of validity of this approximation, we note that the flux piercing an island of size L × L when applying a perpendicular magnetic field B is Φ/Φ 0 = B[T]·(L[nm]) 2 · 2.4 × 10 −4 . Typically a magnetic field of the order of 1T is required to get a substantial Zeeman splitting and tune into the topological SC phase, so L must be smaller than ∼ 64nm. We note that in such small sizes the Coulomb energy of an isolated island may be large; however, if the island is thick and the effective coupling between the islands is strong, this effect may be neglected 19 , as is done throughout this paper.
The system we investigate may be viewed as the superconducting generalization of Hofstadter's model. In addition to the known commensurability effects which are related to the magnetic unit cell, here the geometry of the superlattice plays a role. Namely, in addition to the unit plaquette there is also a SC unit plaquette (see Fig. 1). The orbital effect of the magnetic field modifies the state of the SC, and therefore the ground-state phase configuration of the SC order parameter is not uniform. We find this phase-configuration numerically and also by drawing an analogy to the simpler, exactly-solvable frustrated XY model 20 .
Our model describes a p x +ip y SC, with an application of a non-zero perpendicular magnetic field. Besides being interesting on its own merit, this model has experimental relevance as the orbital effects are almost inevitable 21 .
The magnetic field may induce vortices in the SC (or between the SC islands); since the p x +ip y SC can be tuned into a topologically non-trivial phase (i.e. one that supports chiral Majorana edge modes), a Majorana bound state is expected to reside at the core of each vortex. In the effective low-energy description of the system, the matrix elements between the Majorana bound states at close vortices are non-zero. The coupling between them may give rise to a chiral Majorana mode at the edge of the whole sample.
In this manuscript, we show that Majorana edge modes indeed appear, and the system can be driven into a global topological phase, by the application of an external perpendicular field. This is formally shown in Sec. II by calculating the topological invariant of symmetry class D in two dimensions -the Chern number [22][23][24] N , which is by definition equal to the number of chiral Majorana edge modes. When this Chern number is even, two Majorana modes can pair up to form one fermion edge mode, and the Chern number N /2 is similar to that of the integer quantum Hall effect. However when the Chern number is odd, one unpaired Majorana edge mode exists, and the system enters the topological superconducting state.
Our main result is the phase diagram shown in Fig. 2, where the parameter space includes substantial regions with odd Chern number. This diagram exhibits a series of quantum Hall transitions, where the Chern number N jumps in steps of two, as we sweep the chemical potential keeping the pairing potential ∆ = 0. At finite ∆ the jumps split to two "Majorana transitions" in which the Chern number changes by one. At large ∆ the system becomes topologically trivial, which we attribute to the localization of Cooper pairs in islands.
We then turn to examine the universal features of the phase diagram in Sec. III. We show that even though the system under study is very rich and complex, and depends on many parameters, its phase diagram shares many generic features with the phase diagrams of much simpler models. The fact that the SC pairing potential is staggered is found in Subsec. III A to be important to yield zero Chern number for large pairing potential. The universal properties are understood by considering the symmetries of such models, which depend (for a fixed external field) on two parameters, e.g., p 1 and p 2 , the analogs of ∆ and µ. When p 1 = 0, a sweep of p 2 induces transitions between even Chern numbers in steps of two. Relieving the constraint on p 1 splits the transitions, and exposes regions with odd Chern numbers. As p 1 increases the Chern number decreases in unit steps until it vanishes for large p 1 , making the system topologically trivial. To elucidate these features, we study a simplified version of the 2D model in Subsec. III B, and a generalization of Kitaev's chain model 25 in Subsec. III C. These two models yield phase diagrams which are similar to that of the full 2D model (see Figs. 7,8).
We conclude with a discussion of the main results of the paper. We comment on possible realizations of the effective p-wave pairing we have assumed here in systems x y FIG. 1. Illustration of the islands array: a 3D view is shown in the left panel and a top view (of a part of the system) is depicted in the right panel. Superconductivity exists only in the red squares. In the right panel blue circles denote lattice sites of a tight-binding model of the system. Nearest-neighbor hopping is allowed and uniform throughout the entire lattice. Inside each island, the phase of the superconducting order parameter is approximately constant and its direction is illustrated by an arrow. The magnetic flux per plaquette is Φ, and the flux per SC plaquette is ΦSC, as illustrated by the gray squares. We define the parameters g = Φ h/e which is the flux per unit plaquette and f = Φ SC h/2e which is the flux per SC plaquette. In the geometry depicted here, f = 8g (see Appendix B).
with spin-orbit coupling, and ways to detect the topological phase using heat transport.
Details of the calculations are relegated to the following appendices. In Appendix A we provide additional details on the numerical approach used to find the groundstate phase configuration. In Appendix B we elaborate on the geometrical properties of the superlattice. In Appendix C we review the method used to calculate the Chern number in real space. In Appendix D we present several additional numerical results in the 2D system to support our main results.

A. Model
We introduce the following tight-binding Hamiltonian for single-species fermions on a two-dimensional square lattice, with superimposed p x + ip y SC islands: Here c † m,n , c m,n are creation and annihilation operators of fermions in the lattice site labelled by (m, n), whose location is a (mx + nŷ) where a is the lattice constant; µ is the chemical potential; t is the hopping amplitude between nearest-neighboring sites; Φ is the flux per unit plaquette, which enters the Hamiltonian via the Peierls substitution 26 using the Landau gauge A = Φxŷ/a 2 ; ∆ is the magnitude of the induced superconducting pairing potential; and θ j is the phase of the superconducting order parameter in the j th island, which is approximately constant inside the island.
The Hamiltonian can be brought to the Bogoliuobvde-Gennes (BdG) form by defining the Nambu spinor where H 0 corresponds to the normal µ, t terms and H ∆ = −H T ∆ corresponds to the SC terms. H BdG has particlehole symmetry {H BdG , Λ} = 0 with Λ = τ x K, where τ x is a Pauli matrix acting in particle-hole space and K is complex conjugation. The magnetic field and the chiral p-wave pairing terms break time-reversal symmetry [H BdG , K] = 0, placing H BdG in symmetry class D with a Z topological invariant [22][23][24] .
In the Hamiltonian Eq. (1), the phases {θ j } are treated as parameters. In reality, due to the interaction between the SC and the magnetic field, the ground state configuration of these phases should be determined in a selfconsistent way. This can be done numerically, as elaborated in Appendix A. Due to the computational complexity of this method, we employed a simpler, approximate way to obtain the ground state phase configuration. We drawing an analogy between the current system and the two-dimensional frustrated XY model 20 -a classical theory of U (1) spins with a nearest-neighbor exchange interaction in the presence of a gauge field. Our model is composed of an array of Josephson junctions, and therefore, upon integrating out the fermionic degrees of freedom, the low-energy effective theory for the bosonic SC phases {θ j } has the same symmetry properties as the frustrated XY model. The effective flux per SC plaquette f = ΦSC h/2e is related to the actual flux per plaquette g = Φ h/e by geometric factors (see Appendix B for elaboration); concretely, f = 8g for the geometry of Fig. 1. We therefore substitute the known ground state of the frustrated XY model 27 with flux f as the phase configuration of our system.

B. Chern number
Having established the ground state configuration, we shall henceforth treat the phases {θ j } in Eq. (1) as fixed parameters and study the properties of the resulting Hamiltonian. A direct way to classify the topological properties of the system is calculating the Z topological invariant -the Chern number [22][23][24] . Upon casting the Hamiltonian in BdG form, the Chern number N counts the number of chiral Majorana modes at the edge of the sample; an odd Chern number therefore corresponds to an unpaired chiral Majorana edge mode, which is the hallmark of topological superconductivity.
The Chern number may be calculated in momentum space 28 or in real space 29,30 . In this work we used a realspace calculation, since it is computationally more efficient and numerically robust (the details of the real-space calculation are explained in Appendix C). In particular, since the SC phases vary in space, going to momentum space is not very beneficial. All numerical calculation schemes are prone to finite size or resolution effects; we minimized these by using large enough systems such that in the absence of superconductivity (i.e. ∆ = 0) we obtained the known Chern numbers for the original Hofstadter model 31 . We then scanned the parameter space by varying µ and ∆ (t is held constant, fixing the energy scale), and calculated the Chern number for each set of parameters. The resulting phase diagram for f = 1 2 is shown in Fig. 2.
Several important features appear in the phase diagram At large ∆ the Chern number tends to zero. The reason for this is the islands structure, namely that the SC covering is partial. For large ∆, such a structure gives rise to localized pairs which only realize a topologically trivial phase. This argument is further supported in Subsec. III A, where we study the effect of SC staggering analytically.
In order to complement the Chern number phase diagram we examined the energy gap as a function of µ and ∆ (see Fig. 10 in Appendix D), and found that near µ = 0 and at finite ∆ the system becomes gapless. This is because for µ = 0 and ∆ = 0 the system undergoes a transition of a large Chern number difference, resulting in a very small gap. Finite ∆ then mixes the states near µ = 0, resulting in a gapless state. Notice that the Chern number is ill-defined in such a gapless state, so the area of the phase diagram close to µ = 0 is unreliable. The method we used to calculate the Chern number always yields and integer, as the system is finite so a gap always exists; however, this integer number has no real physical meaning. Therefore, in Fig. 2 we shade regions where the energy gap is smaller than the finite-size level spac- ing, estimated by the tight-binding bandwidth t divided by the total number of sites. For completeness we show the Chern numbers of these regions in the phase diagram, although the true behavior there is metallic.

C. Robustness to flux
We now turn to test the dependence of the above results on the magnetic flux threaded through the system. Up to this point the flux was f = 1 2 flux quantum per SC plaquette; we shall now examine two additional test cases, f = 2 5 and f = 1 3 (corresponding to g = 1 20 and g = 1 24 respectively), which give rise to simple enough unit cells that allow numerical investigations. The Chern number phase diagrams for these two values of the flux are shown in Fig. 3.
These phase diagrams are similar to the one obtained for f = 1 2 (see Fig. 2). The effect of small ∆ -splitting of the two-fold transitions -persists. The effect of very large ∆ -destroying the topological phase -persists as well. The specifics do depend on the flux; in particular, the number of different phases at ∆ = 0 varies. Nevertheless, the general behavior appears to be flux-independent. . Both diagrams share their generic features with these of Fig. 2 which was obtained for f = 1 2 . Substantial areas of odd Chern number, indicating a topological SC phase, can be seen at intermediate values of ∆. Regions of small energy gap (smaller than the bandwidth t divided by the total number of sites) are shaded. For display purposes, we used linear interpolation between the sampled points.
Remarkably, it seems that regardless of f , substantial areas of odd Chern number are formed in the parameter space at intermediate values of ∆. We have shown this behavior for three values of f . The generic pattern we observed suggests that by tuning µ and ∆, a topological phase can be stabilized for various values of f .

III. UNIVERSALITY OF THE PHASE DIAGRAM
Having calculated the phase diagrams in the previous section, we next wish to understand their generic properties. Some details vary between different values of f , for example the shape of the boundary between different topological phases, the value of the gaps and the available Chern numbers of the phase diagram are not identical (see Sec. II and Figs. 2, 3). However, we find three features appearing in all cases: (i) For large pairing potential ∆ the system always becomes trivial. (ii) For ∆ = 0 the Chern numbers are even; as we change µ, twofold jumps of the Chern number occur at the transitions between the phases. (iii) As we increase ∆ for fixed µ, successive transitions occur with mostly a reduction, and sometimes an increase, of the Chern number by one in each transition, until we reach the trivial phase with zero Chern number.
Motivated by these observations, we argue in this section that these features are universal and are shared also by other models with staggered ∆. To demonstrate this, we analyze simplified models in both one and two dimensions that can be treated semi-analytically.
We first show that when turning off the magnetic flux in the islands model, the system becomes topologically trivial for large ∆. We attribute this feature to localization of the Cooper pairs in regions with large pairing potential. We show that in a simplified model of the 2D system that includes the magnetic flux, where the SC phase configuration can be approximately derived, the universal features appear in the phase diagram.
The observation regarding localization of Cooper pairs applies also in one dimension, as we show analytically. This is demonstrated by analyzing a generalized onedimensional Kitaev model with staggered pairing potential. To obtain a system with large Chern numbers in 1D, we cascade several of these staggered Kitaev chains. The resulting system is still one dimensional, as the number of chains is kept finite. The system belongs to symmetry class BDI, which means it possesses a Z topological invariant [22][23][24]32 , just like the two-dimensional class D system we studied in Sec. II. We use this analogy to engineer a system with the same pattern of topological transitions as the full 2D system, by judiciously controlling the coupling between the Kitaev chains.

A. The effect of staggered SC
We consider a generalization of Kitaev's chain model 25 , which is described by the Hamiltonian where c † j creates a single-species fermion at the j th site of the 1D lattice, µ is the chemical potential, t is the (real) nearest-neighbor hopping amplitude, and ∆ j is the SC pwave pairing potential between sites j and j + 1. When ∆ j is uniform, i.e. ∆ j = ∆ for all j, the model 25 supports a topological phase with Majorana end modes for |µ| < 2t for any ∆ = 0.
Consider now a staggered SC version of this model, where the pairing potential takes the values ∆ j = ∆ for odd j and ∆ j = 0 for even j, as illustrated in Fig. 4(a). In momentum space, this modulation corresponds to an enlarged unit cell: (we take the lattice constant to be unity). By analyzing the Pfaffian 25 of the Hamitlonian H st (k) at k = 0, π we find that the topological phase appears for µ 2 + ∆ 2 2 < (2t) 2 , see the ellipse in Fig. 4(b). This means that SC staggering limits the range of ∆ values for which the topological phase survives, reminiscent of the result obtained for the original 2D system. Indeed, this staggered SC configuration can be thought of as the 1D analog of the SC islands configuration we studied in the 2D case. The above result is in sharp contrast with the uniform SC case, where the phase boundaries do not depend on ∆, provided it is non-zero. In order to further establish the effect of a staggered SC, we also studied a simple 2D model where a similar phenomenon occurs. Consider the p x + ip y model on a square lattice where the SC only exists in islands as illustrated in Fig. 5(a). This is a simplification of the model we started with to the case of zero magnetic flux. For simplicity we choose each island to contain four sites (two along the x direction and two along the y direction). The topological properties of the system may then be analyzed by the Pfaffian of the Hamiltonian at the time-reversal invariant momenta (k x , k y ) = (0, 0), (0, π), (π, 0), (π, π). The phase diagram for this model is shown in Fig. 5(b), and it exhibits similar features to those of the 1D model. In particular, the range of ∆ supporting the topological phase is limited, unlike the uniform p x + ip y case. We have further verified that in order to get this effect, it is enough for ∆ to be staggered in any way -islands, stripes, or any other staggered patterns.

B. Semi-analytical 2D model
Having established the importance of a staggered ∆, we now introduce a 2D model with magnetic flux The phase diagram shows the sign of the Pfaffian (which equals +1 for the trivial phase and −1 for a topological phase, with a Chern number +1 or −1) multiplied by the energy gap as a function of the chemical potential, µ and the SC pairing potential ∆. The phase boundary is marked in a black dashed line for clarity, and we set a cutoff scale for the energy gap so the important details are clear. The staggering of the SC limits the range of ∆ in which the system is in the topological phase. The phase diagram shows the sign of the Pfaffian (which equals +1 for the trivial phase and −1 for the topological phase) multiplied by the energy gap, as a function of the chemical potential µ and the SC pairing potential ∆. As in the 1D case (see Fig. 4), the staggering of the SC limits the range of ∆ in which the system is in the topological phase.
which can be understood almost entirely from analytical considerations. We consider a system of coupled SC stripes, with a constant perpendicular magnetic field, see Fig. 6. The virtue of this model is that in the strongly anisotropic limit, where the stripes can be treated independently, we can analytically determine the ground- Illustration of the stripes system (top view) described in Eq. (7). Blue circles denote lattice sites in the tight-binding model, and p-wave superconductivity is introduced only on the red areas. A magnetic field B is applied perpendicular to the system's plane. The ground-state phase configuration of the SC in this model may be analytically approximated, see Eqs. (9a)-(9b).

state SC phases.
The Hamiltonian takes the form The j th stripe lies at m = 2j −1, 2j. The Peierls phase in the Landau gauge is θ P (m, n) = −2π Φ Φ0 m. In the ground state, the gauge-invariant phase difference between every two SC bonds is zero 33 , yielding The θ x SC phases reside at half-integer x and integer y (they connect two adjacent sites along x), and vice versa for θ y SC . Therefore, in our gauge choice, a reasonable configuration is given by Unlike the more sophisticated model studied in Sec. II, here the SC phase configuration is known. It is indeed an approximation -each stripe is treated independently of the others -but at least in the limit t y t x it is sensible. Determining the phase configuration correctly is of paramount importance: The phase diagram, and in particular the topological regions, are extremely sensitive to the SC phases (this is also demonstrated in Appendix D for the full model). Armed with the phase configuration, we are in a position to analyze the Chern number, just as was done in Sec. II. The topological phase diagram for this model, as a function of the chemical potential µ and the pairing potential ∆, is shown in Fig. 7. The phase diagram shares the important features mentioned in Subsec. II B with the phase diagrams of the full model (see Figs. 2, 3). In particular, it displays the crossover from the Hofstadter transitions at ∆ = 0 to the localized superconducting state at large ∆, passing through a series of topological phases at intermediate ∆.
The stripes structure, chosen here for its simplicity, is not special. We tested two additional variants of this semi-analytical 2D model: one where the SC stripes are replaced by SC islands, and one where the flux is only threaded in the normal regions between the SC stripes. We have also corroborated the results by applying different magnetic fluxes. The phase diagrams in all of these variants are similar, and most importantly they all share the features we refer to as universal.

C. Stacking several staggered Kitaev chains
Let us now turn our attention to 1D. Using the staggered Kitaev chain Eq. (5) with ∆ = µ as a building block, we present a way to construct a phase diagram similar to the ones of the full 2D model shown in Sec. II. This is achieved by stacking M staggered Kitaev chains, each having a different hopping parameter t. In the ab-sence of coupling between the chains, we get independent Chern number transitions at the µ values appropriate for each t j (j = 1, . . . , M ) and ∆. Since we want to mimic the behavior of the full 2D system, we take an even number of chains and choose their hopping parameters to be pairwise equal, i.e. t 1 , t 1 , t 2 , t 2 , . . . , t M/2 , t M/2 . This choice guarantees that in the absence of inter-chain coupling, the Chern number may only change by an even number at each transition.
Next, we introduce normal (non-SC) coupling of strength w between neighboring chains. The role of this coupling is similar to that of ∆ in the 2D case: when turned on, it splits the two-fold Chern number transitions into single transitions, where the Chern number changes by ±1. At large µ or large w, the Chern number vanishes due to the SC staggering. Qualitatively, we infer that this construction reproduces the universal features we mentioned before.
In order to make the phase diagram of this model even more similar to that of the original 2D system, we perform a slight modification to w. A single staggered chain with µ = ∆ undergoes a transition of the Chern number from −1 to 1 at µ = 0. Thus, our coupled system undergoes a transition from −2M to 2M at µ = 0 without passing at zero Chern number. Let us then choose w = w 0 / √ µ + µ 0 , where w 0 is the bare coupling strength and µ 0 is a small constant. This way, at µ → 0 the actual coupling w becomes strong, thus driving the Chern number to zero. This transformation, which can be seen as a redefinition of the axes, gives rise to the phase diagram shown in Fig. 8. The resulting phase diagram is similar to the 2D phase diagrams shown in Figs. 2, 3, 7. We therefore conclude that this simple 1D model captures most of the generic features of the complicated 2D models. This observation supports the notion of universality in the phase diagrams, which applies to systems related by dimensionality and symmetry properties: class D in 2D and class BDI in 1D.

IV. DISCUSSION
In this paper, we studied the orbital effects of a magnetic field on the topological properties of the p x + ip y SC. Motivated by relevant experiments, we focused on a model of SC islands arranged in a square lattice. After determining the ground state configuration of the SC, we derived the topological phase diagram for different values of the magnetic flux (see Figs. 2, 3).
The tunable parameters in our model are the magnetic flux Φ, the chemical potential µ, and the strength of the induced SC pairing potential ∆. Our results suggest that for general values of the flux, it is possible to tune into a topological phase supporting Majorana edge modes by varying µ and ∆. The regions in µ, ∆ space supporting the topological phase are substantial, so the parameters do not have to be extremely fine-tuned.
Experimental detection of the Majorana edge modes Topological phase diagram for the coupled staggered Kitaev chains, as a function of the chemical potential µ and the coupling strength w0. This specific system is composed of six coupled chains, with hopping parameters t = {1, 1, 2, 2, 3, 3} and µ0 = 0.01. When the chains are decoupled, i.e. w0 = 0, the Chern number is even, since the hopping parameters were chosen to be pairwise equal. Finite w0 splits these transitions, much like ∆ does in the full 2D case (cf. Figs. 2, 3, 7). Regions of small energy gap are shaded.
can be done by interference [34][35][36] . In addition, these Majorana modes can be detected by measuring the heat conductivity, which is expected to be a half-integer multiple of π 2 k 2 B T /3h in the presence of a chiral Majorana edge mode 37 . We also note that in the model we analyzed, the islands are connected by few sites, without any normal regions between them. The universality of the phase diagram we found suggests that the inclusion of normal regions between the islands, which exist in several experimental realizations, may not alter the main features of the phase diagram.
Our model should be understood as an effective description of a more complicated physical system. We assumed the existence of an induced p-wave pairing potential and single-species fermions. One can take a more microscopic approach and study spinful electrons with spinorbit coupling, proximity coupled to an s-wave SC with an applied Zeeman field, which can give rise to induced p-wave pairing. It is also possible (at least numerically) to account for disorder, which is not expected to have a drastic effect on the results provided it is smaller than the energy gap. Likewise, our zero-temperature study may be generalized to finite temperatures, but as long as the SC gap is larger than the temperature, our results are expected to hold qualitatively.
We also studied simplified, analytically solvable models in 1D and 2D which share common features with the original 2D model. We found that the islands structure (and more generally the staggering of the SC pairing potential) has a profound implication on the topological phase diagram: Pair localization at large ∆ drives the system away from the topological phase. We showed that a model consisting of coupled staggered 1D Kitaev-like chains exhibits a topological phase diagram which greatly resembles those of the original 2D model (see Fig. 8). These findings suggest that several features of the phase diagram are shared by many models. The results shed light on those obtained for the original model, which relied mainly on numerics.
This notion of universality may also be viewed from a general mathematical perspective. Consider a model with two tunable parameters p 1 and p 2 , such that at p 1 = 0 the Chern number has to be even. Assume further that at large |p 1 | and |p 2 | the Chern number has to vanish, and that at p 1 = 0 the generic behavior is onefold Chern number transitions. Under such settings, it is almost impossible to construct a phase diagram which is topologically distinct from those shown in Figs. 2, 3, 7, 8. We conclude that the phase diagrams we found have universal features and may appear under very general settings. for {θ j }. We compared this numerical method and the known frustrated XY ground states for several fluxes, and found that the configurations are qualitatively similar. We note that the optimization is sensitive to the initial conditions, and many local minima exist, so multiple runs are necessary in order to get a good agreement with the results of the frustrated XY model. The numerical optimization shows that the ground-state configurations of the frustrated XY model are at least local minima of our model. Let us denote the number of lattice sites per dimension per island by S i where i = x, y is the spatial dimension, and the number of plaquettes separating adjacent islands per dimension by V i (see Fig. 9). Then, the area of a superlattice unit cell is A SL = (S x + V x − 1) (S y + V y − 1), and the area of a magnetic unit cell is A M = (S x − 1) (S y − 1). The flux per magnetic plaquette is g = Φ h/e and the flux per superlattice unit cell is f = ΦSC h/2e (see Fig. 1). The factor of 2e comes from the fact that the relevant flux quantum for the superconductor is that of Cooper pairs. The relation between f and g is thus In the current study we used S x = S y = 2 and V x = V y = 1, yielding the relation f = 8g. The largest value of f we investigated is 1 2 and therefore the largest value of g is 1 16 , so the approximation of a constant phase inside each island is reasonable. Here we review the method devised in Refs. 29,30 of calculating the Chern number from the real space Hamiltonian. Consider a 2D lattice with N = L x L y unit cells, and denote their positions by r = (x, y) where x, y are integers. Let us define twisted periodic boundary conditions by where ϕ m θ (x, y) are the single-particle wavefunctions (m = 0, . . . , M − 1 where M is the number of electrons), which are vectors in the space of inner degrees of freedom. The many-body ground-state wavefunction Ψ θ ({ r i }) is the Slater determinant of the single-particle wavefunctions of the occupied states, and the Chern number is given by where T θ is the torus 0 ≤ θ x , θ y ≤ 2π. Going to momentum space, we denote by F m θ ( k) the Fourier components of ϕ m θ (x, y). The twisted boundary conditions dictate the allowed momenta, Denoting F m q ( k (0) ) ≡ F m θ ( k), we note that the manybody wavefunction in momentum space Φ q k (0) i is just the Slater determinant of F m q ( k (0) ). Thus, upon substituting ∂ θµ = L −1 µ ∂ qµ (here µ = x, y) we obtain an equation for the Chern number: where R q is the rectangle 0, 2π Lx × 0, 2π Ly . Using Stokes' theorem, the above expression can be written as a winding number along the boundary ∂ Rq , Next, we divide ∂ Rq into N q small line segments and replace the derivatives and integral by their discrete counterparts to obtain where C m,n α,α+1 = F m qα |F n qα+1 are the coupling matrices and q α are the endpoints of the small line segments.
If L x , L y 1, i.e., we consider a sufficiently large system, it is enough to take N q = 4, corresponding to θ = 0 (i.e. just the four corners of ∂ Rq ). In this case we only need to deal with quantities calculated for periodic boundary conditions. Transforming the coupling matrices back to real space for θ = 0 we obtain C m,n α,α+1 = ϕ m θ=0 |e i( qα− qα+1)· r |ϕ n θ=0 .
We may thus define the matrixC = C 0,1 C 1,2 C 2,3 C 3,0 , diagonalize it to obtain the eigenvalues {λ m }, and calculate the Chern number: If the Hamiltonian is written in the BdG form, this method yields the BdG Chern number N .
Appendix D: Additional numerical results in the 2D system We now present several additional results in the full 2D system, described by the Hamiltonian Eq. (3). These results were obtained numerically, and they support the arguments given in Sec. II of the main text: the vanishing of the energy gap near µ = 0 and the strong dependence of the topological phase diagram on the SC phases. Fig. 10 shows the energy gap as a function of the model's parameters for f = 1 2 . Compared with the phase diagram Fig. 2, it is evident that the topological phases are gapped. In addition, this map makes the transitions of the Chern number transparent -they are accompanied by a closing of the energy gap. Fig. 11 shows the Chern number phase diagram for f = 1 2 , in the case where the SC phases are uniform, i.e. θ j = 0 for all j. This phase configuration is not the ground state due to the presence of the magnetic field. Though it resembles the true phase diagram Fig. 2 in its general shape, this phase diagram does not support odd Chern number phases. This observation emphasizes the role of the SC phases and the vortices induced by the magnetic field. FIG. 11. Topological phase diagram for f = 1 2 flux quantum per superconducting plaquette, as a function of the chemical potential µ and the induced SC pairing potential ∆, for the case where the SC phases are uniform. The phase diagram has the same "domes" structure of the ground-state phase diagram (see Fig. 2), but it almost does not support odd Chern number phases.