Suppression of odd-frequency pairing by phase-disorder in a nanowire coupled to Majorana zero modes

Odd-frequency superconductivity is an exotic phase of matter in which Cooper pairing between electrons is entirely dynamical in nature. Majorana zero modes exhibit pure odd-frequency superconducting correlations due to their specific properties. Thus, by tunnel-coupling an array of Majorana zero modes to a spin-polarized wire, it is in principle possible to engineer a bulk onedimensional odd-frequency spinless s-wave superconductor. We here point out that each tunnel coupling element, being dependent on a large number of material-specific parameters, is generically complex with sample variability in both its magnitude and phase. Using this, we demonstrate that, upon averaging over phase-disorder, the induced superconducting, including odd-frequency, correlations in the spin-polarized wire are significantly suppressed. We perform both a rigorous analytical evaluation of the disorder-averaged T -matrix in the wire, as well as numerical calculations based on a tight-binding model, and find that the anomalous, i.e. superconducting, part of the T -matrix is highly suppressed with phase disorder. We also demonstrate that this suppression is concurrent with the filling of the single-particle excitation gap by smearing the near-zero frequency peaks, due to formation of bound states that satisfy phase-matching conditions between spatially separated Majorana zero modes. Our results convey important constraints on the parameter control needed in practical realizations of Majorana zero mode structures and suggest that the achievement of a bulk 1D odd-ω superconductivity from MZMs demand full control of the system parameters.


I. INTRODUCTION
As was originally demonstrated by Berezinskii, 1 the superconducting pair amplitude of fermions may also be odd under exchange of time coordinates, or, equivalently, frequency ω, which extends the usual classification of superconducting (SC) pair symmetries at equal times. 2,3 Part of the importance of odd-frequency (odd-ω) pairing is based on it allowing for highly unconventional SC dynamical correlations, even suggested to be an instance of hidden order. 4 Examples of oddfrequency pairing include long-range proximity effect in superconductor-ferromagnet junctions, 5,6 paramagnetic Meissner effect, 7,8 and Majorana zero modes (MZMs) in topological SCs. 9,10 MZMs are their own antiparticles and emerge isolated as end states in certain topological SCs. 11 They currently hold particular relevance for potential applications in topological quantum computation. [11][12][13][14] Remarkably, MZMs are also a clear example of states that carry pure odd-ω spin-polarized correlations. This is due to their intrinsic nature of being their own antiparticles that implies that their particle-hole propagator is at the same time a particle-particle propagator, which for a zero-energy state automatically acquire an odd-ω dependence. 9,10,15,16 A promising experimentally feasible platform for topological SC, and therefore for MZMs, is based on combining nanowires (NWs) with Rashba spin-orbit coupling, a strong magnetic field, and proximity induced s-wave SC. [17][18][19][20] Other proposals for MZMs include magnetic atoms on spin-orbit coupled SCs, 21,22 and proximitycoupled SC on metallic edge states of topological insula-tors. 23,24 These physical realizations have motivated an enormous interest since the initial experiments. [25][26][27][28][29] Due to the vast experimental progress on MZMs, a promising direction to realize odd-ω SC is using MZMs. In fact, it was recently proposed how to engineer a bulk one-dimensional (1D) odd-ω spinless s-wave SC by coupling an array of MZMs to a spin-polarized wire (SPW). 30 The setup was also shown to exhibit a paramagnetic Meissner effect, whereby an applied magnetic flux is enhanced rather than screened by the induced supercurrent. Additionally, the induced pure odd-ω SC state was demonstrated to be robust against disorder in the coupling coefficients, which, nonetheless, were assumed to be real. However, under more realistic conditions any couplings to MZMs can easily become complex. [31][32][33][34] In fact, in an array of MZMs, it is very unlikely to have equal or even real valued couplings due to the large amount of parameters possible to tune in topological SC wires. [25][26][27][28][29] This then implies that it is crucial to analyze the robustness of the induced odd-ω pairing also under complex and different couplings, and in particular disorder with respect to the complex phase of the coupling terms.
In this work, we consider an array of MZMs coupled to a SPW through generalized complex couplings. In particular, we focus on the role of disorder in the phases of the complex couplings. We find that, apart from local odd-ω pairing, there generically also exist non-local odd-and even-ω correlations. All these amplitudes correspond to the equal-spin triplet (spin-polarized) symmetry as they arise as an induced effect due to MZMs. Most importantly, based on both analytical and numerical approaches, we demonstrate that even moderate phase disorder in the couplings drastically suppresses the induced SC in the SPW, in stark contrast to the effect of real disorder where even large variations of the couplings has been shown to only produce very small changes. 30 Furthermore, we show that the suppression of the induced pair correlations is accompanied by the filling of the energy gap in the density of states (DOS) due to formation of bound states in the SPW. As a consequence, our work shows that experimental realizations of a 1D odd-ω superconductor from an array of MZMs must have full control of the phases of the couplings.
The remainder of this article is organized as follows. In Sec. II we present the model, introduce the distribution of complex couplings, and discuss the Green's function approach for the pair correlations as well as a numerical tight-binding setup. In Sec. III we analytically calculate the pair amplitudes for one and two MZMs coupled to the SPW, and also perform a perturbative analysis for an array of MZMs. Complementary to this, in Sec. IV we carry out a tight-binding numerical study in order to support our analytical findings. Finally, we summarize our work in Sec. V. For completeness, we also provide additional supporting details in four Appendices.

II. MODEL AND METHOD
To create a 1D odd-ω superconductor we consider an SPW that is tunnel coupled to an array of MZMs. The MZMs are modeled as end modes of semiconducting NWs that are in a topologically non-trivial state, similar to Ref. 31. Each NW carries one MZMs at each of its two ends, but only one of them is coupled to the SPW, as illustrated in Fig. 1. The Hamiltonian H of the whole system is therefore composed of three parts describes the SPW with c † x (c x ) creating (destroying) a spin-polarized fermion in the SPW andξ(x) = −( 2 /2m)∂ 2 x − µ w is a ordinary parabolic dispersion relation, with µ w being the chemical potential that determines the filling of the SPW.
The tunnel coupling of the SPW with MZMs at the topological NW end points is described by a tunneling Hamiltonian: 31 where Γ n represents the coupling strength between the MZM γ n and the fermionic mode in the SPW c x at position x n = n b, where b is the separation between two neighboring MZMs. The MZM γ n = f n + f † n = γ † n is one of the two MZM formed from a zero-energy Bogoliubov excitation f † n of the topological NW. For topological NWs of finite length the two MZMs split from zero-energy due to hybridization, 35-38 which is commonly described by: 11 where δ n is the energy splitting in the n-th NW. 39 We consider a general situation in which the coupling strengths Γ n are allowed to be all different from each other and also complex quantities, which describe a realistic situation for an array of MZMs coupled to a SPW. [31][32][33] Here, the SC phases in the topological NWs can easily be one of the sources for the complex Γ. Moreover, under realistic conditions, complex couplings Γ might also appear just due to different parameters in the topological NWs. For example, recent studies in systems with spin-orbit coupling and magnetic field, have demonstrated that the coupling between a MZM and normal systems (quantum dots) gives rise to a tunneling element that is both spin-dependent and complex. 34 Thus, since the couplings depend on the parameters of each NW and it is very hard to the imagine the NWs in an array all having exactly the same parameters, the couplings will generically be complex and vary randomly. We are therefore interested in analyzing the robustness of the induced SC correlations into the SPW against disorder in the couplings, and particularly phase disorder, not previously considered. For real and equal Γ n , Eq. (2) reduces to the model used in Ref. 30, where real disorder was also studied.
More specifically, we model the complex couplings Γ n = |Γ n | e i θn as random variables with the logarithm of the magnitude ln(|Γ n |/Γ 0 ) being uniformly distributed and the scaled phase θ n /π being independently distributed according to a weighted uniform distribution around 0 with weight (1 + s)/2 and around 1 and −1 with weight (1 − s)/4, with in total the width D φ of the 2π phase interval covered. We thus characterize phase disorder by D φ , where D φ = 1 means full phase disorder. Physically, a finite value of D φ indicates the lack of control over the phases in the complex couplings, an issue impossible to avoid experimentally. We also introduce a sign bias, i.e. weight of positive values minus weight of negative values, equal to s, −1 ≤ s ≤ 1. There is nothing fundamental in this choice of distribution, and, in fact, other distributions that model phase disorder would produce similar conclusions. However, this choice has the advantage of allowing the moments |Γ| 2 − |Γ| 2 q , and Γ q to be expressed analytically, exhibiting also a high suppression with increase in the power q and width of the phase distribution D φ . In particular, the expectation value: where sinc(x) = sin(π x) π x . Another advantage of this distribution is that the case of real couplings strengths is achieved as a continuous limit by simply taking D φ → 0. In that case, Eq. (4) reduces to [(1 + s) + (1 − s) cos(π q)] /2, as expected for a distribution of real numbers with a sign bias s. More details on the parameters of the distribution are given in Appendix A.

A. Dyson's equation for the SPW
The focus of this work is on the pair correlations induced in the SPW, and thus it is convenient to work with a Nambu-Gor'kov Green's function in imaginary time: where ψ(x) = c x , c † x , and the imaginary-time Heisenberg operators are ψ(x, τ ) = e τ H ψ(x) e −τ H and ψ(x , τ ) = e τ H ψ † (x ) e −τ H . With only the relative time τ − τ being relevant we perform a Fourier series expansion with respect to the relative time into Matsubara frequencies ω m = (2n + 1) π T . In this way, the Green's function is a 2 × 2 matrix in particle-hole space: The equation of motion forĜ(x, x ; i ω m ) requires the evaluation of the commutator of H = H w +H Γ +H δ , given by Eqs. (1)-(3), with c x (and c † x ). When this is evaluated, terms proportional to c x (or c † x ) and γ n = f n + f † n are obtained. Thus, the equation of motion for G contains the correlators f n |c † x ω , f n |c x ω , f † n |c † x ω , and f † n |c x ω . In order to eliminate these cross-correlators, another set of commutators, of H with f n (or f † n ), is used, which generates terms proportional to f n (or f † n ) and i 2 δ(x − n b) Γ n c † x + Γ * n c x . In this way eliminating the cross-correlators, while using the Green's function for isolated NW modes: where δ mn is a Kronecker-delta in the NW position, we finally arrive at the Green's function for a SPW coupled to an array of MZMs, dependent on two spatial coordinates and Matsubara frequency and satisfying the Dyson equation: Here all the effects of the n-th topological NW are contained within the vertex: In Eq. (7)Ĝ 0 is the propagator in the SPW given by: where the first diagonal element correspond to the free electron propagator in the SPW given by g e (x − x ; i ω m ) = −i sgn(ω m )(m/k 0 )e i sgn(ωm) k0 |x−x | with k 0 = 2 m (µ w + i ω m ) and the second term to the hole propagator in the SPW g h (x − x ; i ω m ) = −g e (x − x; −i ω m ).
An equation forV n ·Ĝ(n b, x ) is obtained by multiplying Eq. (7) byV n from the left, and setting x = n b. The solution of this equation is then expressed as: where the matrixT nl (i ω m ) plays the role of a kernel of the equation: Here we introduce the shorthand notation (ĝ 0 ) nl (i ω m ), which is the Green's function for a bare SPW evaluated at the discrete subset of positions of the MZMs, i.e.: Finally, the solution of Eq. (7), given Eqs. (10) and (11), may be written as: (13) making it clear thatT nm (i ω m ) plays the role of a Tmatrix for the SPW with respect to the effect of the MZM array. The Green's function given by Eq. (13) allows the calculation of the induced pair correlations in the SPW and is used for the analytical calculations in Sec. III.

B. Numerical tight-binding setup
Beyond studying the MZM array using Green's functions we also adopt a tight-binding version of the model described by Eqs. (1)-(3) to facilitate numerical results going beyond the perturbative analytical treatment.
Here we place M MZMs on a ring with z − 1 lattice sites of the SPW inbetween two MZMs and one SPW site at each MZM site, with the periodic boundary conditions c n+z M = c n for the SPW modes and f n+M = f n for the Bogoliubov modes in the n-th NW through which the MZM is expressed. The tight-binding version of Eq. (1) is then given by: For the clean, non-disordered system we can further impose periodic boundary conditions with period z, corresponding to a lattice constant equal to the MZM distance b. This leads to folding of the original SPW dispersion relation ξ(p) = −2t cos(p/z) − µ, z times as p goes over from −z π to z π onto the new unit cell with Brillouin zone from −π to π. Thus, specifying the band index l (0 ≤ l ≤ z − 1) of the conduction band cut by the chemical potential, as well as the (dimensionless) Fermi wavevector −π < p F ≤ π, fixes the chemical potential µ to: The meaning of Eqs. (15) is sketched in Fig. 2.
The tight-binding version of the tunnel Hamiltonian Eq. (2) becomes: where the factor 1/ √ z is necessary if we require |Γ| to be the energy gap in the single-particle excitation spectrum in case of a periodic non-disordered MZM array Γ n = Γ, just as in the continuum case.
We can now define a large 2 M (z + 1)-dimensional Nambu spinor Ψ = c n , c † n |f m , f † m , such that the tightbinding Hamiltonian may be written as a bilinear form correlators are then in principle calculable. In particular, the tight-binding version of Eq.
, which using the eigenvalues and eigenvectors can be calculated as: From this we can directly extract the even-and odd-ω anomalous components G E,O eh as well as the LDOS N i (ε) using: where η is a small positive smearing factor of a Dirac delta function.

III. ANALYTICAL RESULTS
In this part we demonstrate that disorder in the phases of the complex couplings are detrimental for the induced pair amplitudes. Moreover, we show that when there is a finite relative phase between the couplings of at least two MZMs, bound states are induced in the SPW. In order to show these arguments, we consider three simple, yet useful, cases that provide a clear visualization of the role of complex couplings on the pair amplitudes.

A. SPW coupled to a single MZM
We start by considering the simplest case, a single MZM coupled to the SPW. In this case the system of equations Eq. (11) is actually a single equation involving 2 × 2 matrices in Nambu space, which is readily solved: where 2 |Γ| 2 is an even function of frequency ω m . The only difference between Eq. (19) and Eq. (8) is thus the substitution of the denominator D. This clarifies the meaning of the Tmatrix as involving a series of coherent events between the MZM and the SPW of scattering as an electron and back-scattering as a hole. Then, by using Eq. (13), the Green's function for the SPW is easily obtained, with the entries specified in Eq. Eq. (6) given by: where T ij are the elements of theT matrix given by Eq. (19).
The diagonal (normal) elements of the Green's function allow the calculation of the LDOS, while from the off-diagonal (anomalous) terms we extract information about the pair correlations induced in the SPW. Since we are interested in the induced pair correlations, we now discuss G eh which written out is: where we have used the expressions for the free particle propagators in the SPW given in Eq. (9), T eh = −(i ω m /D)Γ 2 , and k 0 is the electron wave vector defined after Eq. 9. The pair amplitude G eh thus exhibits translational invariance breaking through the exponent k 0 |x| − k * 0 |x |, which mixes electron (k 0 ) and hole (k * 0 ) wavevectors with different spatial coordinates. This is similar to the breaking of spatial parity in other superconducting junctions. [40][41][42] In order to visualize this point it is useful to write the wave vector k 0 , defined after Eq. (9), in the Rearranging the exponent of Eq. (21), results in k F (|x|−|x |)+i κ (|x|+|x |) with κ = (ω m k F )/(2µ) and k F = √ 2m µ (we choose a system of units in which = 1). Plugging this expression into Eq. (21), leads to an exponential decay term with decay length given by the inverse of κ and an oscillatory term governed by the Fermi wave-vector. This oscillatory term can further be written in terms of cosines and sines, which are even and odd-functions in space, respectively, and thus leads to the coexistence of even-and odd-ω components of the pair amplitude. 41,42 Interestingly, locally in space, i.e. x = x , G eh is purely odd in frequency, an induced effect that comes directly from the MZM. The odd-ω dependence is here captured in the element T eh , since its denominator D is an even function of ω m . We therefore conclude that the induced local pairing due to a single MZM is odd in frequency, in agreement with previously reported the odd-ω nature of Majorana pair amplitudes. 9,15,30,41,[43][44][45][46] Furthermore, from Eqs. (20) we observe that the anomalous elements are proportional to Γ 2 . This stems from Eq. (19), where the diagonal elements T ee,hh are proportional to |Γ| 2 , while the off-diagonal T eh is proportional to Γ 2 . This is in notable contrast to the findings of Ref. 30, where the pair amplitude was assumed to always be proportional to real couplings. If Γ is real, it does not present any issues for the pair amplitudes, however, it does when Γ is complex. In fact, by performing a disorder average over the phase of Γ, the expression for G eh might give zero or be exponentially suppressed, depending on the disorder distribution. Of course, there is no physical meaning of performing disorder average over the phase of Γ when only a single MZM is coupled to the SPW. However, this simple example still clearly shows the potentially detrimental effect of disorder averaging over the phases of Γ on the pairing amplitudes. As we will see in the following subsections, this seemingly innocent effect at this level has profound consequences when considering more MZMs coupled to the SPW.

B. SPW coupled to two MZMs
The simplest system where the relative phase of the coupling strengths Γ has a physical manifestation is a SPW coupled to 2 MZMs with distance b and complex couplings Γ 1,2 . Then, Eq. (11) is a 4×4 matrix in site and Nambu spaces, but still analytically solvable. Following the same method as for a single MZM we obtain the evenand odd-ω pair amplitudes where g e,h correspond to the free electron and hole propagators in the SPW defined in Eq. (9). Moreover, the de-nominatorD is given byD . Based on these results we first conclude from Eqs. (22) that both even-and odd-ω pair correlations are induced in the SPW due to its coupling to MZMs. Both pair amplitudes also acquire a dependence on the phases of the complex couplings Γ 1,2 , as evident in Eqs. (23). Moreover, since these amplitudes correspond to equal-spin spin-triplet pairing, the even-frequency component is finite non-locally but vanishes locally, which is similar to the even-ω pairing discussed for the single MZM in the previous subsection. On the other hand, we find that odd-ω pairing exists both locally and non-locally. Both even-ω and odd-ω pair amplitudes exponentially decay from the position of the MZMs with a decay length determined by the inverse of κ, specified in the previous subsection. Second, we conclude that the pair amplitudes given by Eqs. (22) emerge proportional to Γ 2 1(2) or Γ 1 Γ 2 through the numerator of the T ij , N eh ij , given by Eqs. (23). This, therefore, implies that any average over the complex phases of Γ 1,2 leads to a suppression of the induced pair correlations in the SPW, demonstrating a detrimental role of complex couplings.
Before concluding this part, we explicitly demonstrate the effect of a finite relative phase θ between the couplings of the two MZMs Γ 1,2 . This is illustrated in Fig. 3, where we plot the LDOS in the SPW as a function of frequency for nine values of the phase θ uniformly distributed on the segment [0, 2π]. We calculate the LDOS from the standard expression LDOS(x, ω) = (−1/π)ImG ee (x, x; ω + iη), with ω real frequency, η an infinitesimal positive number, and G ee the electron-electron Green's function obtained after solving Eq. (13). The system is symmetric with respect to the midpoint x = b/2 between the two MZMs and we therefore plot the LDOS at three representative points on the side closer to the first MZM, where we measure the MZM distance b in units of the Fermi wavelength λ F = 2π/k F in the SPW.
The most prominent feature in Fig. 3 is the presence of resonance-like peaks at ±δ, corresponding to the leaking of the MZM states into the SPW. The weight of these peaks is highly dependent on position, having the highest weight (anti node) at b/4, while the lowest weight (node) is directly under the MZM and also in the middle between the two MZMs at b/2. These correspond to delocalized states in the SPW with a half-wavelength (nodeto-node distance) b/2, i.e. a wavelength equal to b. Since in this plot we have chosen b = λ F , it is not clear from these results which of these length scales determines the wavelength of the bound state like feature. We therefore investigate other MZM separations in App. B and find that the wavelength is determined by λ F . In addition, it is clear that the spectrum is symmetric in θ relative to θ = π. The peaks are sharpest for θ = 0, π, and 2π, and have the largest smearing for θ = π/2, 3π/2. The smearing, associated with a finite lifetime of these resonance states, is highly phase-dependent, which is reminiscent of the conditions for constructive and destructive interference. To summarize, these results suggest that delocalized bound states appear naturally systems with multiple MZMs due to the complex couplings, but that their behavior is also naturally expected to be less controllable as the number of Majorana NWs increases.

C. SPW coupled to a MZM array
Having studied one and two MZM coupled to the SPW, we finally turn to the case of an infinite array of MZMs and for this system calculate the induced pair amplitudes in the SPW using disorder-averaging to capture the intrinsic variability in the complex phases of the MZM couplings.
For an infinite array, it is convenient to express Eq. (11) in terms of the Fourier series coefficients: where p, p are continuous variables on the interval (−π, π]. The functions defined in Eqs. (24) are all periodic in each p with a period 2π. With this Eq. (11) acquires the form: This is an integral equation forT (p, p ; i ω m ) as a function of continuous variables, which reflects the now infinite dimension of the matrix equation Eq. (11). Nevertheless, it is a convenient formulation for three reasons. First, the periodicity with a lattice constant b imposed on the continuum wire by the MZM array is reflected in the foldedĝ 0 (p; i ω m ). Expressing it in terms of the continuum Green's functionĜ 0 (k, i ω m ), the following relation holds: In fact, due to the simple parabolic dispersion relation, the sum over m in Eq. (26) may be performed analytically with the result: where k 0 is the electron wavevector defined after Eq. (9). Second, the Fourier-transformed vertex coefficientŝ V (p − p ; i ω m ) become approximately normally distributed when disorder is taken into account, regardless of the exact distribution of the real-spaceV n (i ω m ) according to the Central Limit Theorem. This allows us to apply Wick's theorem for an average of a product of several such terms in an expansion.
Third, after disorder averaging, the T -matrix recov- Using the Fourier series formulation we can now, from Eqs. (13), (24a), (24c), and (28), express the disorderaveraged Green's function for the SPW as: Here we note that the second term incorporates the effects of BZ folding through the matricesĝ 0 (p; i ω m ), with integration over p only in the interval from −π to π. As a direct consequence, the effects of disorder averaging are now completely encompassed within the disorderaveraged T -matrix. Finding a general expression for the disorder-averaged T -matrix is straightforward but lengthy and we therefore provide the details as Appendix C.
We are here primarily interested in the induced anomalous correlations in the SPW, characterized by the electron-hole component of the Green's function. To obtain them, according to Eq. (29) and the fact that bothĜ 0 (k; i ω m ) andĝ 0 (p; i ω m ) are diagonal in Nambu space, only the electron-hole component of the disorderaveraged T eh is needed. The detailed derivation of expressions for T within the second Born approximation (i.e. keeping only the n = 1 term in Eq. (C8), which is of second order in δV ) are provided in Appendix D, and here we produce schematic, but, nevertheless sufficiently detailed expressions to draw the most important conclusions. In particular, for the electron-hole component we obtain whereg e,h are given by Eqs. (27). The denominator D T is here an even function of ω m , but its explicit form is not necessary for our discussion. Similarly, the exact form of φ 1 (i ω m ), given by Eq. (D2b), andĝ 1 (i ω m ), given by Eqs. (C4b), (D3), is not important for our discussion. Although the expression for the disorder-averaged pair amplitude G eh that follows from combining Eqs. (29) and (30) is not simple, some conclusions can be drawn already from analyzing the components of T eh . For example, we obtain that A 2 = |Γ| 4 , B 2 = |Γ| 2 Γ 2 , and C 2 = Γ 4 , while A 1 = |Γ| 2 and B 1 = Γ 2 , see Appendix D for details on their evaluation. In the generic situation, we expect D φ to be close to 1. Phase disorderaveraging thus forces C 2 B 2 A 2 . This directly results in the A 2 -term being the dominant term in the square brackets expression in Eq. (30). However, this term also contains a factor of g eh 1 that, according to Eqs. (C4b) and (D2) is also proportional to B 1 . Thus, even this dominating term is in fact highly suppressed after phase-disorder averaging, because B 1 A 1 . As a consequence, the disorder-averaged pair amplitude in Eq. (30) is small, and we can conclude that both the evenand odd-ω induced pair correlations is highly suppressed in an MZM array.
Additionally, we have verified that the formation of bound states due to the coupling between the SPW and the array of MZMs follows the case discussed in previous section for two MZMs. In fact, we find bound states emerging at energies around the Majorana splittings δ n . However, unlike the case of two MZMs, where the two isolated modes are not able to open a gap in the continuum of the SPW spectrum, in the case of an array such a gap is opened, and the in-gap states also do not acquire smearing, but, are shifted in energy instead. We briefly discuss this feature in the following section.
We conclude this analytical part of our work by sum-marizing our findings: A MZM coupled to a SPW induces even-and odd-ω correlations in the SPW, whose amplitude is proportional to Γ 2 . Locally in space only oddω pairing exists. When two MZMs are coupled to the SPW with different complex couplings, delocalized bound states emerge in the SPW between the two MZMs, and both pair amplitudes become dependent on the phase difference and proportional to Γ 2 1(2) or Γ 1 Γ 2 . This already on the level of two MZMs indicate a suppression of the induced pairing in the presence of phase disorder in couplings Γ i . Finally, for an infinite array, we verify that the induced pair amplitudes in the SPW are necessarily considerably suppressed under disorder-averaging over the phases of the complex couplings.

IV. NUMERICAL RESULTS
In order to demonstrate the effects of phase-disorder, without relying on the approximations that were essential in obtaining closed-form expressions in the analytical approach in Sec. III C, such as the Born approximation for the infinite array, and also going beyond a qualitative interpretation of disorder effects, we next perform a numerical diagonalization of a tight-binding model described in Sec. II B and add disorder explicitly.
We choose a ring of M = 40 MZMs, with z = 10 atomic sites between two neighboring MZMs. Since the folded bands alternate from electron-like to hole-like around momentum p = 0, we center the chemical potential to cut the third band (l = 2), as it is the first higher electronlike band, and choose a k F = 0.2 π. We then also scan six equidistant values of the chemical potential, ranging from the third to the fifth band with the same k F , whose numerical values are given in Fig. 4(d). This assures that our results are not fine-tuned but can be interpreted as general. As the band index increases, the bandwidth also increases (see e.g. Fig. 2). We therefore choose Γ 0 to be equal to 0.4 of the smaller of the distances from the Fermi level to the band edge for the lowest scanned band, which turns out to be Γ 0 = 0.0815 t. Since the effects of disorder in the coupling coefficients are investigated, we fix the energy splitting in each NW to a small finite value δ = 0.1 Γ 0 .
First of all, we verify that keeping the disorder real (D φ = 0) reproduces the robustness of the induced SC correlations, as also found in Ref. 30, and that regardless of the value of s (the sign bias of the disorder), and also for relatively large magnitude disorder, up to D r = 1. Then, with the value of Γ 0 fixed, we generate 20 different random disorder configurations according to the distribution described in the introduction of Sec. II, using the parameters s = 1, D r = 0, and with D φ ranging through the values listed in Fig. 4. Here D r is a dimensionless measure of the disorder in the coupling given by the relative error of |Γ| 2 (see Appendix A). The results of these calculations are summarized in Figs. 4 and 5.
In Fig. 4 we plot the magnitude of on-site, and thus dominating, odd-ω component of the sample-averaged anomalous Green's function G O eh i,i (i ω m ) as a function of Matsubara frequency ω m , site location i, and D φ . The overall ω m -dependence exhibits a drastic suppression with phase-disorder strength D φ , as evident in the logscale plot in Fig. 4(a). We stress that this behavior is thus very different for the effect of real disorder, where the pair amplitudes are robust. 30 The data in Fig. 4(a) are the odd-ω component directly under one of the MZMs. In Figs. 4(b) and (c), we plot the site dependence over 6 out of the 40 periods for both the lowest positive and a high frequency component. For low frequencies we see very little site variations. This can be understood from our analytical results for the single MZM, Eq. (21), where the decay length is ∼ 1/ω m . 30,47 At higher frequencies, here represented by Fig. 4(c) evaluated at half the bulk gap Γ 0 /2, we find more variations with position. The anomalous Green's function G O eh i,i (i Γ 0 /2) actually has peaks (approximately) at a quarter distance from a MZM with dips both under a MZM, as well as at a mid-point between them. 48 Finally, in the last panel Fig. 4(d) we show how the low-frequency on-site component is suppressed with disorder for several different values of the chemical potential. Despite the doping-dependence of the DOS, we find a very similar monotonic suppression with phasedisorder strength, which shows that our results are not sensitive to the specifics of the model situation. Beyond the dominant on-site odd-ω contribution to the anomalous Green's function plotted in Fig. 4, there are also sub-dominant non-local contributions, in analogy to the analytical results in Sec. III B, which have both evenω, and odd-ω dependence and are peaked at different positions along the chain depending on the particular choice of model parameters. However, they are equiva- lently suppressed with the dominant on-site contribution as the phase disorder increases.
In Fig. 5 we focus on the DOS in the SPW and plot the the site-averaged DOS for several values of the phase disorder D φ , for the same parameters as in Fig. 4(a-c). With increasing disorder we see that the near-zero-energy peaks (split off zero energy because a finite δ) is diminishing but most notable is that the filling of the energy gap at finite energies. A crucial difference here compared to conventional effects of pair-breaking disorder is that the gap gets filled by smearing of the near-zero-energy peaks and by states near the gap edge. Thus the in-gap states are naturally associated with the MZMs, and we interpret their shift in energy as the formation of non-local bound states corresponding to phase-matching conditions, similarly to the simplest example of two MZMs considered in Sec. III B.
Another feature that is prominent in Fig. 5 is the large peak that occurs near −Γ 0 . This peak is more robust to the disorder than the energy gap. By carefully inspecting the site-resolved LDOS, we identify this as a contribution coming from sites that are at a quarter of the MZM distance from a particular MZM. Under an MZM, and at the midpoint between two MZMs, this peak has actually the smallest weight. Note here that we plot the electron DOS thus the DOS is not needed to be symmetric around zero energy. The missing DOS for creating a symmetric spectrum are hole-like quasiparticles.

V. CONCLUSIONS
In this work we have investigated the robustness of the superconducting pair correlations induced into a spinpolarized wire (SPW) from an array of Majorana zero modes (MZMs) when the couplings between these two systems are complex and acquire different phases. This corresponds to a realistic situation, as complex couplings generically appear for varying system parameters and disorder in the phases of the complex couplings is experimentally unavoidable.
First, we have demonstrated that, in general, the pair correlations induced into the SPW exhibit both even-and odd-frequency (ω) dependence, which exponentially decay from the position of the MZMs. We have shown this effect for a single MZM, a pair of MZMs, and for an array of MZMs. Interestingly, we find that the phases of the complex couplings get transferred to the pair amplitudes.
Second, we have shown that the induced pair correlations, including the odd-ω component, suffer a considerable suppression due to phase-disorder averaging, exactly as a consequence of the transferred complex phases. This is in stark contrast to the effect of real disorder where the pair amplitudes remain robust. 30 We have found this strong disorder dependence both by analytically evaluating the T -matrix within the second Born approximation and performing numerical tight-binding calculations for an array of MZMs.
Third, we find that the suppression of the pair correlations in the SPW occurs concurrently with the filling of the energy gap by in-gap bound states appearing between spatially-separated MZMs. Carefully analyzing the situation of a pair of MZM coupled to a SPW, we have found bound states between the two MZMs that are highly sensitive to the relative complex phase between the two couplings. We can thus conclude that it is the complex phases and their disorder that is also causing the filling of the energy gap.
In summary, our work demonstrates that the conditions favorable for practical realization of a bulk 1D oddω superconductivity from MZMs requires full control of the system parameters.
treat the magnitude and phase of Γ as independent. As for |Γ|, being a non-negative quantity, it is convenient to use the (natural) logarithm ln(|Γ|/Γ 0 ), which spans the whole real line. We choose the distribution of the logarithm to be uniform U(f min , f max ), with the bounds f min and f max chosen such that: which leads to the the following expressions: x coth(x) = 1 + D 2 r .
Here, Γ 0 is to be considered as a constant effective coupling that would reproduce the same T -matrix in Eqs. (8) and (11), D r is a dimensionless measure of the disorder in the coupling given by the relative error of |Γ| 2 . In the case of real couplings, the signs of Γ do not enter Eq. (8), and, thus, play no effect on the SPW. Nevertheless, we model a sign bias by a parameter s, −1 ≤ s ≤ 1, being the difference in the weight of positive and negative signs.
For creating a distribution of complex couplings Γ = |Γ| e i θ , we choose, strictly by convenience, that the phase follows a weighted uniform distribution around 0 with weight (1 + s)/2, and symmetrically around π and −π with weight (1 − s)/4: where W (x) = (1/2) Θ(1−|x|) is a window function from −1 to 1 with unit weight. Here D φ is the total fraction of the 2π phase interval that is covered. We use the distribution P (θ) to numerically model phase disorder in the couplings in Section IV of the main text.

Appendix B: Bound states for two MZMs
In this Appendix we present complementary information for the delocalized bound states in the case of 2 MZMs coupled to a SPW, originally discussed in Fig. 3. In Figs. 6 and 7 we present the same spectrum as Fig. 3, but now for b = 0.5 λ F , and b = 2 λ F . Clearly, the relative position of the nodes and antinodes changes with this change of MZM-MZM distance. In Fig. 6, the nodes occur at x = 0 (and x = b), while there are peaks or antinodes both at x = b/2 and at x = b/4 (and x = 3b/4), but clearly highest at x = b/2. This suggests that the wavelength of the standing wave, λ s , is λ s /2 = b, or λ s = 2 b.
Because b = λ F /2 in this case, λ s = λ F . On the other hand, in Fig. 7 it there are nodes on all panels x = 0, x = b/4, x = b/2 (and, by symmetry at x = 3b/4 and x = b as well). This implies 2 λ s = b, or λ s = b/2. Because b = 2 λ F in this case, λ s = λ F . This shows how the periodicity is set by the Fermi wavelength as stated in the main text.  In this Appendix, we derive a working expression for T (p; i ω m ), the disordered-average of the matrix T acting as an impurity T -matrix. The starting point is Eq. (25) in the main text, which acquires the shorthand form: Here, every symbol is a matrix in Nambu and momentum space, and matrix multiplication involves an integration over a dummy momentum variable ( π −π dp 2π ), and the 1 is an identity matrix in both Nambu and momentum space, with the Dirac delta in momentum space containing an extra factor of 2π (1 → 2π δ(p − p )1). With this, all the rules of matrix algebra are readily applicable to Eq. (C1). In particular, Eq. (C1) can be rewritten in an expanded form: Expressing both T and V as their expectation values plus a deviation from the expectation value we get: Then, taking the expectation value, and keeping in mind that δV = δT = 0 by definition, we arrive at the following equation for the expectation value: Subtracting this equation from the full equation, we arrive at the equation for the deviations: We see here that to find T we need to evaluate the higher moment δV · g 0 · δT in Eq. (C2). To calculate that average, we reformulate Eq. (C3) in a more convenient form in terms of g 0 · δT as: The newly introduced matrix g 1 incorporates the averaged out effects of the MZM array and is diagonal in momentum space. Using Eq. (C4b), Eq. (C2) may be rewritten as: Now by multiplying Eq. (C4a) by δV from the left and taking an average, we see that δV · g 0 · δT is expressed in terms of the new averages δV · g 1 · δV , and (δV · g 1 ) · δV · g 0 · δT : δV · g 0 · δT = δV · g 1 · δV + + δV · g 1 · δV · g 0 · T + + (δV · g 1 ) · δV · g 0 · δT .
As a consequence it means that the following recursive relation for the averages X n = (δV · g 1 ) 2n · δV · g 0 · δT , involving the expectation values Π * n = (δV · g 1 ) 2n−1 · δV , holds: where we have also used Eq. (C5). Summing from n = 0 to ∞ in Eq. (C6), and defining Π * = ∞ n=1 Π * n , we get the following equation for X 0 : Plugging this equation into Eq. (C5), and keeping in mind that g 1 · V equals g 1 · g −1 0 − 1 according to Eq. (C4b), we have: Eq. (C7) simultaneously defines a "dirty" discrete Green's function g for the SPW, incorporating the effects of any disorder in the matrices V n on top of an averaged-out homogeneous coupling V , as well as an implicit equation for the "dirty" discrete T -matrix T . At this point it is customary to introduce the oneparticle-irreducible self-energies: where with the Wick contractions we retain only the ones that cannot be disconnected by "cutting" a g 1 line. It is obvious that this self energy satisfies: which implies a Dyson equation for g : By using Eqs. (C4b) and (C10) we may write: which defines the self-energy. Then, using the second equality in Eq. (C7) for T , after some algebra, we obtain: If this equation is compared with an alternative form of Eq. (C1), namely T −1 = V −1 − g 0 , we see that the role of V −1 after averaging over disorder is played by Σ −1 , or equivalently, the inverse of the self-energy is the inverse of V .