Phase tunable second-order topological superconductor

TL;DR

This paper proposes a phase-tunable two-dimensional nanowire array to realize second-order topological superconductors with Majorana corner states, demonstrating control over topological phases and robustness against disorder.

Contribution

It introduces a novel, tunable 2D nanowire array platform for second-order topological superconductivity with detailed characterization and disorder resilience.

Findings

Coupling control via phase-biasing enables access to various topological phases.

System remains robust against disorder in wires and junctions.

Identification of a topological Anderson phase induced by disorder.

Abstract

Two-dimensional second-order topological superconductors (SOTSCs) have gapped bulk and edge states, with zero-energy Majorana bound states localized at corners. Motivated by recent advances in Majorana nanowire experiments, we propose to realize a tunable SOTSC as a two-dimensional nanowire array. We show that the coupling between the Majorana modes of adjacent wires can be controlled by phase-biasing the device, allowing to access a variety of topological phases. We characterize the system using scattering theory, which provides access to its transport properties and its topological invariants. The setup is robust against disorder, both in the nanowires themselves and in the Josephson junctions formed between adjacent wires. Further, we identify a parameter regime in which an initially trivial system is rendered topological upon adding disorder, providing an example of a second-order topological Anderson phase.

Equations

$$\begin{array}{rl}& {\mathcal{H}}_{\mathrm{wire}}(k_x)=[2t_x(1-\cos k_x)-\mu ]{\tau }_z{\sigma }_0+\\ & V_x{\tau }_0{\sigma }_x+V_z{\tau }_0{\sigma }_z+\Delta {\tau }_x{\sigma }_0+\alpha \sin k_x{\tau }_z{\sigma }_y,\end{array}$$

$$\begin{array}{rl}& E^2={ϵ}_{k_x}^2+V_x^2+V_z^2+{\Delta }^2+{\alpha }^2{\sin }^2{k}_x\\ & \pm 2\sqrt{(V_x^2+V_z^2){\Delta }^2+(V_x^2+V_z^2+{\alpha }^2{\sin }^2{k}_x){ϵ}_{k_x}^2},\end{array}$$

$$\begin{array}{rl}\mathcal{H}(k_x)=& [2t_x(1-\cos k_x)+2t_y-\mu ]{\eta }_0{\tau }_z{\sigma }_0+\\ & V_x{\eta }_0{\tau }_0{\sigma }_x+V_z{\eta }_0{\tau }_0{\sigma }_z+\\ & \Delta \cos \frac{\varphi }{2}{\eta }_0{\tau }_x{\sigma }_0+\Delta \sin \frac{\varphi }{2}{\eta }_z{\tau }_y{\sigma }_0+\\ & \alpha \sin k_x{\eta }_0{\tau }_z{\sigma }_y-2t_y{\eta }_x{\tau }_z{\sigma }_0+\beta {\eta }_y{\tau }_z{\sigma }_x,\end{array}$$

$$\begin{array}{rl}{H}_{\mathrm{TB}}& =\sum _{x,y}{\Psi }_{x,y}^{†}[(2t_x+2t_y-\mu ){\tau }_z{\sigma }_0+V_x{\tau }_0{\sigma }_x+\\ & V_z{\tau }_0{\sigma }_z+\Delta \cos \varphi (y){\tau }_x{\sigma }_0+\Delta \sin \varphi (y){\tau }_y{\sigma }_0]{\Psi }_{x,y}\\ & +\{{\Psi }_{x,y}^{†}[-t_x{\tau }_z{\sigma }_0+i\frac{\alpha }{2}{\tau }_z{\sigma }_y]{\Psi }_{x+1,y}+h.c.\}\\ & +\{{\Psi }_{x,y}^{†}[-t_y{\tau }_z{\sigma }_0+i\frac{\beta }{2}{\tau }_z{\sigma }_x]{\Psi }_{x,y+1}+h.c.\},\end{array}$$

$$\phi(y)=\left\{\ \begin{array}[]{ll}0$ , $y=1,2\\ \phi$ , $y=3,4\\ 2\phi$ , $y=5,6\\ \vdots\\ (\frac{N_{\rm wires}}{2}-1)\phi$ , $y=N_{\rm wires}-1,N_{\rm wires}.\end{array}\right.$$

$$\phi(y)=\left\{\ \begin{array}[]{ll}0$ , $y=1\\ \phi$ , $y=2,3\\ 2\phi$ , $y=4,5\\ \vdots\\ \frac{N_{\rm wires}}{2}\phi$ , $y=N_{\rm wires}.\end{array}\right.$$

$$S={\tau }_y{\sigma }_y{S}^{\ast }{\tau }_y{\sigma }_y.$$

$$Q_l=\mathrm{sign}(\det {\mathrm{r}}_{\mathrm{l}}),$$

$$r_l=\left(\begin{array}{cc}{r}_l^{ee}& r_l^{eh}\\ r_l^{he}& r_l^{hh}\end{array}\right),$$

$$G_A=G_0[N-\mathrm{tr}(r_{ee}{r}_{ee}^{†})+\mathrm{tr}(r_{eh}{r}_{eh}^{†})],$$

$$\begin{array}{rl}& \mathcal{H}=[\frac{-{\hslash }^2{\partial }_x^2}{2m}-\mu ]{\eta }_0{\tau }_z{\sigma }_0+V_z{\eta }_0{\tau }_0{\sigma }_z+\\ & \Delta \cos \frac{\varphi }{2}{\eta }_0{\tau }_x{\sigma }_0+\Delta \sin \frac{\varphi }{2}{\eta }_z{\tau }_y{\sigma }_0-i\alpha {\partial }_x{\eta }_0{\tau }_z{\sigma }_y,\end{array}$$

$$\begin{array}{r}{\mathcal{H}}_{\mathrm{SNW}}=[\frac{-{\hslash }^2{\partial }_x^2}{2m}-\mu ]{\tau }_z+V_z{\sigma }_z+\Delta {\tau }_x-i\alpha {\partial }_x{\tau }_z{\sigma }_y,\end{array}$$

$$\begin{array}{r}{\mathcal{H}}_{\mathrm{SNW}}^{\prime }=\frac{-{\hslash }^2{\partial }_x^2}{2m}{\tau }_z-V_z{\sigma }_y+\Delta {\tau }_x-i\alpha {\partial }_x{\tau }_z{\sigma }_z,\end{array}$$

$$\begin{array}{rl}{\Psi }^{{}^{\prime }†}& =\frac{1}{\sqrt{2}}({\psi }_{\uparrow }^{†}+i{\psi }_{\downarrow }^{†},{\psi }_{\downarrow }^{†}+i{\psi }_{\uparrow }^{†},{\psi }_{\downarrow }-i{\psi }_{\uparrow },-{\psi }_{\uparrow }+i{\psi }_{\downarrow })\\ & \equiv ({\psi }_{\uparrow }^{{}^{\prime }†},{\psi }_{\downarrow }^{{}^{\prime }†},{\psi }_{\downarrow }^{{}^{\prime }},-{\psi }_{\uparrow }^{{}^{\prime }}).\end{array}$$

$${\mathcal{H}}_{\mathrm{T}+\mathrm{SOC}}^{\prime }=-\frac{{\hslash }^2{\partial }_x^2}{2m}{\sigma }_0-i\alpha {\partial }_x{\sigma }_z,$$

$${\tilde{H}}_{\mathrm{T}+\mathrm{SOC}}^{{}^{\prime }}=\sum _{\sigma ,{\sigma }^{\prime }}\int dx{\tilde{\psi }}_{\sigma }^{{}^{\prime }†}(x)e^{i\sigma k_{\mathrm{SO}}x}{\mathcal{H}}_{\mathrm{T}+\mathrm{SOC}}^{\prime }{e}^{-i{\sigma }^{\prime }{k}_{\mathrm{SO}}x}{\tilde{\psi }}_{{\sigma }^{\prime }}^{{}^{\prime }}(x),$$

$${\tilde{\mathcal{H}}}_{\mathrm{T}+\mathrm{SOC}}^{\prime }=-\frac{{\hslash }^2}{2m}({\partial }_x^2+k_{\mathrm{SO}}^2),$$

$${\tilde{\psi }}_{\sigma }^{{}^{\prime }}(x)={\tilde{R}}_{\sigma }^{{}^{\prime }}{e}^{i{k}_{\mathrm{SO}}x}+{\tilde{L}}_{\sigma }^{{}^{\prime }}{e}^{-i{k}_{\mathrm{SO}}x},$$

$$\begin{array}{rl}{\tilde{H}}_{\mathrm{T}+\mathrm{SOC}}^{{}^{\prime }}=-i\alpha \int dx& [{\tilde{R}}_{\uparrow }^{{}^{\prime }†}{\partial }_x{\tilde{R}}_{\uparrow }^{{}^{\prime }}+{\tilde{R}}_{\downarrow }^{{}^{\prime }†}{\partial }_x{\tilde{R}}_{\downarrow }^{{}^{\prime }}\\ & -{\tilde{L}}_{\uparrow }^{{}^{\prime }†}{\partial }_x{\tilde{L}}_{\uparrow }^{{}^{\prime }}-{\tilde{L}}_{\downarrow }^{{}^{\prime }†}{\partial }_x{\tilde{L}}_{\downarrow }^{{}^{\prime }}]\end{array}$$

$$\begin{array}{rl}{\tilde{H}}_z^{\prime }& =-V_z\sum _{\sigma ,{\sigma }^{\prime }}\int dx{\tilde{\psi }}_{\sigma }^{{}^{\prime }†}{e}^{i\sigma k_{\mathrm{SO}}x}({\sigma }_y{)}_{\sigma ,{\sigma }^{\prime }}{e}^{-i{\sigma }^{\prime }{k}_{\mathrm{SO}}x}{\tilde{\psi }}_{{\sigma }^{\prime }}^{{}^{\prime }}\\ & =-V_z\int dx[-i{\tilde{\psi }}_{\uparrow }^{{}^{\prime }†}{e}^{2i{k}_{\mathrm{SO}}x}{\tilde{\psi }}_{\downarrow }^{{}^{\prime }}+i{\tilde{\psi }}_{\downarrow }^{{}^{\prime }†}{e}^{-2i{k}_{\mathrm{SO}}x}{\tilde{\psi }}_{\uparrow }^{{}^{\prime }}]\\ & =-V_z\int dx[-i{\tilde{R}}_{\uparrow }^{{}^{\prime }†}{\tilde{L}}_{\downarrow }^{{}^{\prime }}+i{\tilde{L}}_{\downarrow }^{{}^{\prime }†}{\tilde{R}}_{\uparrow }^{{}^{\prime }}],\end{array}$$

$$H_{\mathrm{SC}}^{{}^{\prime }}=-\Delta \int dx[{\psi }_{\uparrow }^{{}^{\prime }}{\psi }_{\downarrow }^{{}^{\prime }}-{\psi }_{\downarrow }^{{}^{\prime }}{\psi }_{\uparrow }^{{}^{\prime }}+h.c.],$$

$$\begin{array}{rl}{\tilde{H}}_{\mathrm{SC}}^{{}^{\prime }}=-\Delta \int dx& [{\tilde{R}}_{\uparrow }^{{}^{\prime }}{\tilde{L}}_{\downarrow }^{{}^{\prime }}+{\tilde{L}}_{\uparrow }^{{}^{\prime }}{\tilde{R}}_{\downarrow }^{{}^{\prime }}-\\ & {\tilde{R}}_{\downarrow }^{{}^{\prime }}{\tilde{L}}_{\uparrow }^{{}^{\prime }}-{\tilde{L}}_{\downarrow }^{{}^{\prime }}{\tilde{R}}_{\uparrow }^{{}^{\prime }}+h.c.].\end{array}$$

$${\tilde{H}}^{{}^{\prime }l}=\frac{1}{2}\int dx({\psi }^{{}^{\prime }l}{)}^{†}{\tilde{\mathcal{H}}}^{{}^{\prime }l}{\psi }^{{}^{\prime }l},$$

$$\begin{array}{rl}& {\tilde{\mathcal{H}}}^{{}^{\prime }i}=-i\alpha {\tau }_z{\sigma }_z{\partial }_x+\Delta {\tau }_x-V_z{\sigma }_y\\ & {\tilde{\mathcal{H}}}^{{}^{\prime }e}=i\alpha {\tau }_z{\sigma }_z{\partial }_x+\Delta {\tau }_x.\end{array}$$

$$\begin{array}{rl}& {\tilde{\mathcal{H}}}^i=-i\alpha {\tau }_z{\sigma }_y{\partial }_x+\Delta {\tau }_x+V_z{\sigma }_z\\ & {\tilde{\mathcal{H}}}^e=i\alpha {\tau }_z{\sigma }_y{\partial }_x+\Delta {\tau }_x\end{array}$$

$$\begin{array}{rl}& ({\tilde{E}}^i{)}^2=-{\alpha }^2{\kappa }_i^2+(\Delta \pm V_z{)}^2\\ & ({\tilde{E}}^e{)}^2=-{\alpha }^2{\kappa }_e^2+{\Delta }^2.\end{array}$$

$${\tilde{\Psi }}_1^i(x)$$

$${\tilde{\Psi }}_2^i(x)$$

$${\tilde{\Psi }}_1^e(x)$$

$${\tilde{\Psi }}_2^e(x)$$

Full text

I Introduction

The hallmark characteristic of topological phases of matter is the presence of robust zero-energy modes at the interfaces between topologically trivial and nontrivial regions. The existence of these modes is the consequence of a nonzero, quantized value of the topological invariant, a relationship known as the bulk-boundary correspondence. Until recently, this correspondence only predicted the existence of gapless states with dimension lower by one than that of the bulk. The stability of these gapless modes is ensured by the presence of local symmetries such as the particle-hole, time-reversal or chiral symmetry.Hasan and Kane (2010); Qi and Zhang (2011); Bernevig (2013) In the presence of crystalline symmetries, this correspondence holds only for the surfaces left invariant by these symmetries. Fu (2011); Hsieh et al. (2012); Xu et al. (2012)

Recently, research focus has shifted towards higher-order topological phases, Khalaf (2018); Geier et al. (2018); Trifunovic and Brouwer (2019); Benalcazar et al. (2017a, b); Song et al. (2017); Schindler et al. (2018a, b); Langbehn et al. (2017); Wang et al. (2018a); Zhu (2018); Ezawa (2018a, b, c); Dwivedi et al. (2018); van Miert and Ortix (2018); Ezawa (2018d); Hsu et al. (2018); Yan et al. (2018); Wu et al. (2019); Zhu (2019); Wang et al. (2018b); Liu et al. (2018); Peng and Xu (2019); Serra-Garcia et al. (2019); Xie et al. (2018); Peterson et al. (2018); Serra-Garcia et al. (2018); Imhof et al. (2018); Zhang et al. (2019); Franca et al. (2018); Călugăru et al. (2019); Wang et al. (2018c); Benalcazar et al. (2019); Attig et al. (2018); Queiroz et al. (2018); Liu et al. (2019a); Araki et al. (2019); Ezawa (2018e, f, 2019); Ahn and Yang (2019); Bultinck et al. (2019); You et al. (2018); Liu et al. (2019b); Bomantara et al. (2019); Rodriguez-Vega et al. (2018); Huang and Liu (2018); Volpez et al. (2019); Hayashi (2018); Noh et al. (2018); Agarwala et al. (2019); Kheirkhah et al. (2019); Slager et al. (2015); Ghorashi et al. (2019) in which both the D-dimensional bulk and the (D-1)-dimensional boundary are gapped, whereas topological zero modes have a dimension of (D-2) or less. A D-dimensional, Nth order topological phase is defined as hosting boundary modes of dimension (D-N), localized, for instance, at the corners or hinges of a sample. To date, higher-order phases are divided into two groups, labeled intrinsic and extrinsic.Geier et al. (2018); Trifunovic and Brouwer (2019) The former ones require lattice symmetries and are termination independent, as long as the termination is compatible with these symmetries. The lattice symmetries are essential, as they impose unique and opposite mass terms on the neighboring surfaces, thus localizing gapless modes at the corresponding hinges or corners. This is quantified by a nontrivial bulk invariant. On the other hand, extrinsic phases do not require any spatial symmetries. In essence, the bulk is trivial and boundary modes exist as a consequence of a nontrivial surface.Trifunovic and Brouwer (2019)

A majority of proposals for the realization of higher-order topological insulators rely on the presence of crystalline symmetries. For instance, Refs. Benalcazar et al., 2017a, b discussed a two-dimensional/three-dimensional (2D/3D) Su-Schrieffer-Heeger model,Su et al. (1979) where a topological invariant, protected by mirror symmetries, was interpreted as a bulk quadrupole/octupole moment. So far, this model has been realized in meta-materials,Serra-Garcia et al. (2019); Xie et al. (2018) as well as in microwave and electrical circuits.Peterson et al. (2018); Serra-Garcia et al. (2018); Imhof et al. (2018) Furthermore, higher-order topology was used to explain gapless modes along the step edges of a bismuth crystal,Schindler et al. (2018b) protected by the combination of time-reversal symmetry and spatial symmetries like three-fold rotation and inversion.

Topological superconductors (TSCs) also represent a promising venue to achieve higher-order phases. Some proposals rely on $p$-wave superconductivity,Langbehn et al. (2017); Wang et al. (2018a); Khalaf (2018); Zhu (2018) where a bilayer of $p_{x}+ip_{y}/p_{x}-ip_{y}$ superconductors has gapped edges due to the application of a time-reversal symmetry breaking term,Khalaf (2018); Zhu (2018) or the coupling between the layers is such that certain edges are gapped.Langbehn et al. (2017) Other approaches are based on mixed pairing (usually $p+id$ or $s+id$),Wang et al. (2018a); Wu et al. (2019); Zhu (2019) and several realizations have been proposed, including $\rm Sr_{2}RuO_{4}$ under certain conditions. Furthermore, SOTSC phases have been proposed in engineered systems like a 2D topological insulator placed in the proximity of a high-temperature superconductor,Wang et al. (2018b); Liu et al. (2018) or a combination of a 3D anti-ferromagnetic topological insulator and a conventional s-wave superconductor.Peng and Xu (2019) An additional advantage of higher-order TSCs is their potential use as generators of quantum codes.You et al. (2018) Finally, a recent work of Volpez et al. discusses a SOTSC phase induced once the in-plane Zeeman field is present in $\pi$-junction Rashba layers.Volpez et al. (2019)

The experimental realization of topological superconductivity has been an active field of research in the last decade.Qi and Zhang (2011); Beenakker (2013) Initially, it was predicted in $p$-wave superconductors, leading to proposals based on the proximity introduced superconductivity in topological insulatorsFu and Kane (2008); Qi et al. (2010) and semiconductors,Oreg et al. (2010); Lutchyn et al. (2010) both with strong spin-orbit coupling. Among these, the Majorana nanowire, made of a superconductor deposited or epitaxially grown on a nanowire has been the focus of intense experimental efforts. Mourik et al. (2012); Deng et al. (2016); Zhang et al. (2018); Albrecht et al. (2016); Krogstrup et al. (2015); Gazibegovic et al. (2017); Casparis et al. (2018); Vaitiekėnas et al. (2018)

Motivated by the recent advances in the fabrication techniques and characterization methods of Majorana nanowires, we use them as a platform to create an extrinsic SOTSC. The 2D array of coupled Majorana nanowires is already studied in the context of weak-TSCs,Seroussi et al. (2014); Fu et al. (2007); Fu and Kane (2007); Diez et al. (2014) where a set of equally spaced nanowires is proximitized by a single superconductor. Here, we propose a simple design by which the array realizes a SOTSC once the inter-wire couplings are dimerized. Such a dimerization could be the result of an alternating wire spacing, which would remain fixed once the device has been fabricated. Alternatively, we show that a continuously tunable setup is realized if the wires are equally spaced, but proximitized by different superconductors. By controlling the superconducting phase differences between adjacent wires, it is possible to tune the coupling between adjacent Majorana modes such that the array forms a SOTSC. We characterize the system using the scattering matrix approach, which is discussed in the context of second-order topological phases by Geier et al.Geier et al. (2018) Finally, our numerical simulations show that this phase is robust to disorder, and can, in fact, expand in the presence of it. Therefore, our model realizes an example of a secoSzczecinnd-order “topological Anderson phase”.Li et al. (2009); Groth et al. (2009); Song et al. (2012); Lv et al. (2013); Adagideli et al. (2014)

The rest of the paper is organized as follows. In Section II, we do a step-by-step construction of the SOTSC, starting from the model of a single Majorana nanowire. In Section III, we characterize the system using scattering matrices, determining its topological invariant and transport properties. Section IV is devoted to the effects of disorder on the system. Finally, in section V, we summarize our results and draw comparisons to other proposals.

II Models

In this Section, we show how the coupling between Majorana bound states (MBS) belonging to adjacent wires may be tuned by varying the superconducting phase difference, and then exploit this tunability to generate a SOTSC phase in a nanowire array.

II.1 The model of a Majorana nanowire

We begin by reviewing the model of a single nanowire oriented along the x-direction and proximity coupled to an s-wave superconductor. Only the nanowires with strong Rashba spin-orbit coupling (SOC) are relevant as the SOC, once a magnetic field is applied, creates an effectively spinless state, thus mimicking the physics of the Kitaev chain.Kitaev (2001) Then, the Bogoliubov-de Gennes (BdG) Hamiltonian $H_{\rm wire}=\int dk_{x}\Psi^{\dagger}\mathcal{H}_{\rm wire}\Psi$ is defined through the Hamiltonian density

[TABLE]

in the basis $\Psi^{\dagger}=(\psi_{k_{x},\uparrow}^{\dagger},\psi_{-k_{x}\downarrow}^{\dagger},\psi_{k_{x}\downarrow},-\psi_{-k_{x}\uparrow})$, where $\psi_{k_{x},\uparrow}^{\dagger}$ denotes the creation operator of an electron with spin up and momentum $k_{x}$. As this basis is assumed hereafter, in the following we specify only the Hamiltonian densities.

Furthermore, $t_{x}$ and $\alpha$ are the nearest-neighbor hopping strength and the SOC in the x-direction (along the wire), $\mu$ is the chemical potential, $V_{x}/V_{z}$ is the Zeeman energy in the x-/z-direction, and $\Delta$ is the superconducting pairing strength. Pauli matrices $\tau$ and $\sigma$ represent particle-hole and spin degrees of freedom, respectively. This system belongs to class BDI of the Altland-Zirnbauer classification,Altland and Zirnbauer (1997) with particle-hole symmetry $\mathcal{P}=\tau_{y}\sigma_{y}\mathcal{K}$, chiral symmetry $\mathcal{C}=\tau_{y}\sigma_{y}$ and time-reversal symmetry $\mathcal{T}=\mathcal{K}$, where $\mathcal{K}$ denotes complex conjugation.

From the eigenvalue equation, one obtains

[TABLE]

where $\epsilon_{k_{x}}=2t_{x}(1-\cos{k_{x}})-\mu$.

To achieve the topological phase, it is sufficient that only one of the Zeeman terms is nonzero. The external magnetic field applied in the x-direction gives rise to a Zeeman term $V_{x}$, while $V_{z}$ is more experimentally relevant for ferromagnetic atoms deposited on an s-wave superconductor,Nadj-Perge et al. (2013, 2014) where it arises due to the local magnetic moments of atoms.

The topological criterion, $V_{x}^{2}+V_{z}^{2}>\Delta^{2}+\mu^{2}$, that distinguishes a topologically trivial and a nontrivial phase, is obtained by setting $k_{x}=0$. In the case of only one wire, regimes $V_{x}\neq 0;V_{z}=0$ and $V_{z}\neq 0;V_{x}=0$ give the same probability distribution of zero-energy modes. This is due to the fact that the two terms are mapped onto each other ($V_{x}\to V_{z}$ and $V_{z}\to-V_{x}$) by the unitary transformation $(\sigma_{0}+i\sigma_{y})/\sqrt{2}$. This is no longer the case for a 2D system consisting of coupled nanowires with different superconducting phases, where we find that the direction of the Zeeman field affects the number and the probability distribution of gapless modes.

II.2 Two-wire system

We now proceed with the model of two Majorana nanowires with opposite superconducting phases (Fig. 1) describing how the energy of the mid-gap states can be tuned by changing the phase difference.Haim (2019)

The Hamiltonian density takes the form

[TABLE]

where $t_{y}$ is a spin-conserving hopping term between wires and $\beta$ is a term which describes a hopping which flips the electron spin, a consequence of the Rashba SOC in the y-direction. The superconducting phase difference between the wires is $\phi$ and the Pauli matrix $\eta$ denotes the wire space, while all other symbols retain their meaning. Particle-hole symmetry is preserved with the same operator, $\mathcal{P}=\tau_{y}\sigma_{y}\mathcal{K}$, and we can define a new chiral symmetry operator, $\mathcal{C}=\eta_{x}\tau_{y}\sigma_{y}$, and, consequently, $\mathcal{T}=\eta_{x}\mathcal{K}$.

From now on, unless otherwise specified, the parameters used in the calculations are $\mu=0$, $t_{x}=1.7$, $\alpha=2.5$, $\Delta=2.5$, $t_{y}=0.4$ and $\beta=0.8$. Furthermore, we will consider two separate cases: $V_{x}=4,V_{z}=0$ or $V_{z}=4,V_{x}=0$.

Throughout the following, the key insight we will use to generate a SOTSC is that, by altering the phase difference between the wires, for instance by applying a supercurrent, it is possible to tune the splitting between their MBS. This is shown in Fig. 2, where the spectrum of the two-wire system is plotted as a function of $\phi$.111All numerical simulations are performed using Kwant,Groth et al. (2014) and the code is included as part of the Supplementary Material. The splitting depends both on the type of inter-wire coupling (only $t_{y}$, only $\beta$, or both of them), as well as on the direction of the Zeeman field.

If $V_{z}=0$ and $V_{x}\neq 0$ (top panels of Fig. 2), the interaction strength between the pair of MBS on the left and right side of the system is identical. This is due to a mirror symmetry, ${\cal M}_{x}=\sigma_{x}$, such that $\sigma_{x}H(-k_{x})\sigma_{x}=H(k_{x})$. When only the spin-conserving hopping is present, the phase dependence of the mid-gap energy levels takes the form $\pm\sin{\frac{\phi}{2}}$, and the two levels cross at $\phi=0$ (Fig. 2a). On the other hand, when only the spin-flip hopping is included, the crossing occurs at $\phi=\pi$, with end-state energies behaving as $\pm\cos{\frac{\phi}{2}}$ (Fig. 2b). These numerical results are confirmed in Appendix A using first order perturbation theory in $t_{y}$ and $\beta$. When both coupling terms are nonzero, the crossing point will happen at a generic value of $\phi$, which depends on $t_{y}$ and $\beta$ (Fig. 2c). In this case, by tuning the phase difference, the MBS coupling can be made either stronger or weaker than if the wires were proximitized by the same superconductor (meaning $\phi=0$).

Setting $V_{x}=0$ and $V_{z}\neq 0$ results in similar MBS couplings (bottom panels of Fig. 2). The principal difference between the effect of $V_{z}$ and $V_{x}$ is the fact that, for $V_{z}\neq 0$, there is no unique mirror symmetry operator in the x-direction. Instead, one can define two different $\mathcal{M}_{x}$ for each of the couplings $t_{y}$ ($\mathcal{M}_{x}=\sigma_{z}$, Fig. 2d) and $\beta$ ($\mathcal{M}_{x}=\eta_{z}\sigma_{z}$, Fig. 2e). As such, a broken mirror-symmetry phase is realized once both couplings are present, one in which the MBS splitting on the right boundary of the system is different from that on the left boundary. This is shown in Fig. 2f, where the pair of end modes on the right and left is shown in red and blue, respectively. Due to the symmetry $\eta_{x}\sigma_{z}H(k_{x},\phi)\eta_{x}\sigma_{z}=H(-k_{x},-\phi)$, left-end modes at phase difference $\phi$ have the same energy as right-end modes at $-\phi$.

II.3 Many-wires system

Having shown how superconducting phase differences can be used to engineer end-state energies, we now use this control parameter to generate a SOTSC phase. The setup, consisting of an array of equally spaced nanowires, is shown in Fig. 3. We introduce superconducting bridges which connect neighboring wires into different patterns, corresponding to panels a, b and panels c, d of Fig. 3. We assume that these superconducting “connectors” are placed far from the ends of the wires, such that they do not alter the MBS wavefunctions. However, due to these bridges, pairs of nanowires are effectively proximitized by the same superconductor, leading to a vanishing phase difference between them. Disconnected pairs of wires, on the other hand, may still show a nonzero $\phi$, for instance when a supercurrent is passed through the device.

The real-space Hamiltonian written in the basis $\Psi_{x,y}^{\dagger}=(\psi_{x,y,\uparrow}^{\dagger},\psi_{x,y,\downarrow}^{\dagger},\psi_{x,y,\downarrow},-\psi_{x,y,\uparrow})$, where the subscripts $x,y$ denote the positions of sites in the tight-binding model, reads

[TABLE]

where $\phi(y)$ is the superconducting phase that depends on the coordinate $y$. For the setup shown in Fig. 3a, b it takes values

[TABLE]

where, $N_{\rm wires}$ is the (even) number of Majorana nanowires in the y-direction. The other possibility, Fig. 3c, d, has

[TABLE]

In both arrangements, the superconducting phase difference between adjacent superconductors is the same throughout the system and equals $\phi$.

We first discuss the case when only the Zeeman field in the x-direction is present. Given current nanowire fabrication techniques, we expect this to be the most relevant scenario, since the superconducting shell covering the wires is typically very thin, such that a magnetic field perpendicular to the wires could destroy superconductivity.

For the arrangement of nanowires depicted in Fig. 3a, the SOTSC phase is obtained when the phase difference is tuned such that the energy of MBS that originate from different superconductors is increased compared to the case $\phi=0$. On the contrary, the setup of Fig. 3c requires a stronger coupling between the MBS of wires proximitized by the same superconductor. This produces a dimerization of the MBS couplings in adjacent wires, such that the edge perpendicular to the wires (the $y$-edge) becomes a nontrivial Kitaev chain.Kitaev (2001) The result is a phase with a gapped bulk, gapped edges, but Majorana modes localized at the corners. The latter occur due to the fact that both edges are topologically nontrivial: the edge parallel to $x$ is the last topological nanowire of the array, whereas the edge parallel to $y$ is a nontrivial Kitaev chain.

Given the spectrum shown in Fig. 2c, the appropriate range to observe the nontrivial phase in the setup of Fig. 3a is $\phi\in(\sim 2.26,2\pi)$. We choose a value $\phi=4.5$ and calculate the energy spectrum and the probability distribution of zero-energy modes (see Figs. 4a and 4b). As the mirror symmetry in the x-direction is not broken, we expect four zero-energy modes that are pinned to the corners, which is confirmed by our calculations. For the other arrangement of nanowires (Fig. 3c), one should tune the phase difference in the regime where the energy of MBS decreases compared to $\phi=0$ to obtain the SOTSC phase.

For the sake of completeness, we also examine the case of a Zeeman field pointing in the z-direction. Note that an externally applied magnetic field which is perpendicular to the nanowires will lead to a complicated phase relation across the Josephson junctions, in addition to the afore mentioned possibility of destroying superconductivity altogether. It is however possible to generate such a Zeeman term in the absence of external fields, by taking advantage of the ferromagnetic ordering of atoms deposited on a superconductor.Nadj-Perge et al. (2013, 2014) In this case, the behavior of the Josephson junctions may be treated in the same way as before.

For the system in which $V_{z}\neq 0$, Fig. 2f suggests that any nonzero $\phi$ produces a SOTSC phase with only two gapless corner modes, regardless of a particular realization of the system. This is due to the fact that MBS on the left and right edges of the array will be dimerized in an opposite fashion once a small phase difference is applied. For instance, setting $\phi=4.5$, the configuration given by Eq. (4), supports two zero-energy modes (Fig. 4c) located on two right corners (Fig. 4d). The same phase difference for the other configuration would produce two MBS located on the left corners, as this realization requires smaller energy splittings compared to their value at $\phi=0$ (Fig. 2f). The phase difference $\phi=\pi$ is special, however, as there are four corner modes for a configuration as in Fig. 3b and no corner modes for the setup of Fig. 3d.

III Characterization

In order to show the topological nature of the setup shown in Fig. 3, we calculate the topological invariant using scattering theory.Akhmerov et al. (2011); Fulga et al. (2011, 2012a); Geier et al. (2018) This invariant is based on the parity of the number of MBS at one corner and can be used to characterize both first-order Akhmerov et al. (2011); Fulga et al. (2011, 2012a) and higher-order phases.Geier et al. (2018)

To determine the scattering matrix, a unitary matrix that connects incoming to outgoing modes, we consider a four-terminal transport geometry, consisting of a 2D system, given by the Hamiltonian Eq. (3), and four normal-metal leads oriented along the x-direction, each attached to one corner. Each lead $l$ supports four incoming modes $\Psi^{\rm in}_{l}$ and four outgoing mode $\Psi^{\rm out}_{l}$ at the Fermi level, corresponding to spin and particle-hole degrees of freedom.

In this setup, the scattering matrix $S$ has a $4\times 4$ block structure, whose diagonal blocks are reflection matrices $r_{l}$, $l=1,...,4$, and whose non-diagonal blocks are transmission matrices $t_{ln}$ ($l,n=1,..,4$ and $l\neq n$) from one lead to another. In class D,Altland and Zirnbauer (1997) with the particle-hole symmetry $\mathcal{P}=\tau_{y}\sigma_{y}\mathcal{K}$, the scattering matrix obeys

[TABLE]

The procedure to obtain this relation is detailed in Appendix B. Given the relation Eq. (6), it follows that the determinant of the full scattering matrix, $\det S$, as well as that of each reflection block $\det r_{l}$ are real. This allows to define a topological invariant: the sign of the determinant of the reflection matrixAkhmerov et al. (2011); Geier et al. (2018)

[TABLE]

In total, four invariants can be defined, one for each corner, and they are not mutually independent, as MBS always appear in pairs.

The reflection matrix is a matrix block of the form

[TABLE]

where $r^{ee}$ and $r^{hh}$ are $2\times 2$ matrices that represent normal reflections (scatterings of electrons into electrons and holes into holes), whilSzczecine matrices $r^{eh}$ and $r^{he}$ describe Andreev reflection processes.

Once the bulk and the edges of a 2D system are gapped, the incoming mode can only be reflected. Then, reflection matrices $r_{l}$ become unitary due to current conservation and their determinants satisfy $|\det{r_{l}}|=1$. If the incoming mode is reflected on a Majorana bound state, $\det{r_{l}}=-1$ due to a $\pi$ phase shift in the scattered states.Akhmerov et al. (2011) This contrasts the case of no Majorana corner states, for which $\det{r_{l}}=1$.

The phase transition occurs once the y-edge gap is closed and two counter-propagating edge states appear, protected by the translation symmetry. At this point, the topological invariant has to change sign $(Q=0)$, which is possible as the reflection matrix is not unitary anymore. This change of sign is accompanied by the appearance of the quantized value of thermal conductance, defined as G=G_{\rm th}\rm{tr}$$(t^{\dagger}t).Akhmerov et al. (2011) Here, $t$ is the transmission matrix between the two relevant leads and $G_{\rm th}=\frac{\pi^{2}k_{B}^{2}T}{6h}$ is the thermal conductance quantum. Note that this phase transition point separating the trivial and SOTSC phase is precisely a weak-TSC, hosting gapless Majorana modes on the y-edges, modes which are protected by translation symmetry.

In Fig. 5a, we show the value of $\det{r}$ (blue line) as the phase difference between the nanowires is changed for a system, represented in Fig. 3a, with the Zeeman field in the x-direction. Here, we have dropped the index $l$ as all gapless modes have the same dependence on $\phi$. Furthermore, we use $\det{r}$ as it is a continuous function whose values match $Q$ far away from the phase transition. The boundaries of the topological phase agree with the values of $\phi=0$ and $\phi\simeq 2.26$ that produce uniform energy splittings of MBS throughout the system (see Fig. 2c). These transitions, are accompanied by quantized thermal conductance peaks, represented with red lines in Fig. 5a.

We further study the nature of the corner modes by calculating the Andreev conductance, one of the main experimental tools to confirm the presence of MBS.Fu and Kane (2009); Akhmerov et al. (2009); Mourik et al. (2012); Deng et al. (2016); Zhang et al. (2018) The latter lead to resonant Andreev reflection, giving rise to a quantized value of the Andreev conductance $G_{A}=2G_{0}$ according to the relationBlonder et al. (1982)

[TABLE]

where $G_{0}=\frac{e^{2}}{h}$, $\mathop{\mathrm{tr}}$ denotes the trace, and $N$ is the number of incoming electron modes in the lead ($N=2$ in our case). The numerical results of Fig. 5b show that in the SOTSC phase there is a quantized zero bias conductance peak, whereas this peak is absent in Fig. 5c, which represents a trivial phase.

We now discuss the case of the Zeeman field in the z-direction and for the configuration given in Fig. 3b. We expect that the leads attached on left/right corners show different dependence of the topological invariant and the thermal conductance on $\phi$. We plot the topological invariant/thermal conductance calculated using left leads in Fig. 6a and the same quantities for right leads are given in Fig. 6b. On the left side of the system, the phase transition occurs for $\phi=0$ and $\phi\simeq 4.02$. On the right side of the system, the nontrivial phase happens when $\phi\in(\sim 2.26,2\pi)$. As expected for this arrangement of wires, all four leads measure a nontrivial phase at $\phi=\pi$.

IV Disorder

To model realistic systems, we have to take into account the effects of disorder. Although the breaking of the translational invariance does not affect the validity of our topological invariant,Akhmerov et al. (2011) it makes invariants at the corners mutually independent. Various kinds of disorder act differently on the topological phase in Majorana nanowires.Lutchyn et al. (2010); Stanescu et al. (2011); Akhmerov et al. (2011); Brouwer et al. (2011); Adagideli et al. (2014) In particular, it was found that MBS are robust against short-ranged disorder in chemical potential in single-band nanowiresAdagideli et al. (2014) and their multi-band generalizations.Lutchyn et al. (2010); Stanescu et al. (2011) In the following, we will study the effect of chemical potential fluctuations in the nanowires on the SOTSC phase. Furthermore, we anticipate that the distance between the superconductors that form Josephson junctions will not be uniform throughout the realistic system. Such imperfections create fluctuations of the critical currents $I_{C,j}$ ($j$ denotes the junction),Ashcroft and Mermin (1976); Abrikosov (1988) which in turn lead to fluctuations of the phase differences acquired across the junctions, once an external supercurrent $I$ is imposed. Therefore, on the Hamiltonian level, we simulate this type of disorder by randomizing the phase differences of neighboring superconductors.

Disorder in the chemical potential can originate from impurities in the nanowire. In our simulations, we model these fluctuations by replacing $\mu\to\mu+v_{\mu}$, where $v_{\mu}$ is drawn independently for each lattice site from the uniform distribution $[-\chi_{\mu},\chi_{\mu}]$, with $\chi_{\mu}$ the strength of disorder. The phase diagram, as a function of the phase $\phi$ and disorder strength, is represented in Fig. 7a for the arrangement of nanowires given in Fig. 3a and in Fig. 7b for the arrangement of nanowires as in Fig. 3c. In both cases, we consider a SOTSC phase with four gapless modes in the clean limit, while the case of a non-vanishing $V_{z}$ will be addressed in Appendix C. The diagram represents an averaged value of four topological invariants, each calculated for one corner, over $301$ independent disorder realizations. These results are accompanied by the calculations of x-edge/y-edge/bulk conductance, represented in Figs. 7e, 7f and 7g, respectively. The x-edge/y-edge conductance is calculated using transmission matrices between the leads that share the same horizontal/vertical edges, while the bulk conductance is calculated using a system with periodic boundary conditions in one direction. The phase boundaries look the same for the two different system configurations, so we plot only once the dependence of these conductances on $\phi$ and disorder strength.

We notice that the topological phase is stable up to moderate strengths ($0<\chi_{\mu}<4$) of disorder. Furthermore, as smaller values of $\chi_{\mu}$ do not strongly break the translational invariance in the y-direction, the thermal conductance along the y-edges remains close to unity (bright red line in Fig. 7f). Stronger disorder leads to a metallic phase, characterized by the appearance of x-edge and bulk conductance.Medvedyeva et al. (2010); Fulga et al. (2012b) Finally, the system is in the trivial Anderson insulator phase for very large disorder strengths ($\chi_{\mu}>8$), with $Q=1$ for any value of $\phi$. In Fig. 7a, we notice the appearance of the topological Anderson insulator phase,Li et al. (2009); Groth et al. (2009); Song et al. (2012); Lv et al. (2013); Adagideli et al. (2014) as the boundaries of the SOTSC phase around $\phi=2.26$ are extended compared to the clean limit. In effect, this means that a system which it topologically trivial in the clean limit can be rendered nontrivial once disorder is introduced.

As already pointed out, experimental realizations of our system might have different phase differences between the adjacent superconductors. Therefore, we also model disorder in the phases of superconducting pairing. Phase diagrams of setups shown in Figs. 3a and 3c are presented in Figs. 7c and 7d, respectively. More precisely, we take $\phi_{\rm{dis}}(y)=\phi(y)+v_{\phi}$, where $v_{\phi}$ is uniformly distributed in the range $[-\chi_{\phi},\chi_{\phi}]$, such that $\braket{\phi_{\rm{dis}}(y)}=\phi(y)$. As before, we consider $301$ independent disorder realizations. For the configuration proposed in Fig. 3a, we observe that disorder in $\phi$ reduces the extent of the topological phase in $\phi$ as compared to the clean limit. The opposite happens for the arrangement of nanowires like in Fig. 3c, which represents another example of a higher-order topological Anderson insulator. As the two configurations share the same phase boundaries, the y-edge conductance, represented in Fig. 7h, has the same behavior in both cases. For an almost-clean system, this conductance will be close to one, while larger fluctuations in the phases diminish its value. On the other hand, phase disorder does not cause any x-edge or bulk conductance, since for our parameters the bulk states remain gapped for any $\phi$ (see Fig. 2c). We see that these configurations have different tendencies for large disorder strength, i.e. the SOTSC phase either increases or decreases. However, in limit of small phase deviations, the boundary between the two phases remains close to the one in the clean limit, such that a small value of $\phi$ is still sufficient to cause a topological phase transition.

V Conclusion

Higher-order topological phases present an opportunity to explore novel phenomena. In this work, we constructed a 2D model of a second-order topological superconductor controlled by means of superconducting phase differences. The latter lead to dimerized MBS couplings in an array of Majorana nanowires, allowing to access both the SOTSC phase as well as its transition to a trivial phase. Phase-biasing the device may be achieved by means of a supercurrent (as shown in Fig. 3), or alternatively by using flux loops.

From an experimental point of view, our model offers several advantages. First, it only requires conventional, s-wave superconductors, as well as the fabrication techniques recently developed in the context of Majorana nanowires,Krogstrup et al. (2015); Gazibegovic et al. (2017); Casparis et al. (2018); Vaitiekėnas et al. (2018) including the deposition of a superconductor on a nanowire using in situ methods. Second, it allows for a topological phase to be reached for a range of phase differences, without the latter having to be set to a particular value. In fact, an arbitrarily small value of $\phi$ is sufficient to dimerize the MBS couplings of the nanowire array and reach a topological phase, though it may be that corner modes have a very large localization length in this case. Furthermore, if the SOTSC phase is controlled by a Josephson current, the latter could be made time-dependent, opening the possibility to study periodically-driven SOTSC phases.Bomantara et al. (2019); Rodriguez-Vega et al. (2018); Huang and Liu (2018)

As the nontrivial phase is protected only by particle-hole symmetry, our model is robust to fluctuations in both the chemical potential and phase differences between the nanowires. Depending on the configuration of the system, we have shown that the topological phase can expand in the presence of disorder in the chemical potential, leading to a second-order topological Anderson phase. Furthermore, simulations reveal that randomness of the superconducting phases has a similar effect. Thus, depending on the way we couple the superconductors, the topological phase can either enlarge or decrease.

VI Acknowledgments

We thank Ulrike Nitzsche for technical assistance. This work was supported by the DFG through the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147, project-id 39085490).

Appendix A Corrections to energies of gapless modes

For very small inter-wire coupling, the correction to the energies of the non-perturbed Hamiltonian can be found using perturbation theory. Throughout this Section, we closely follow the derivation presented in Ref. Klinovaja and Loss, 2012. This procedure requires knowing the initial wave-functions, and for this reason we start from a system of two uncoupled nanowires with opposite superconducting phases. Then, the BdG Hamiltonian $H=\int dx\Psi^{\dagger}\mathcal{H}\Psi$ is defined through the Hamiltonian density

[TABLE]

in the basis

$\Psi^{\dagger}=(\psi_{1,\uparrow}^{\dagger},\psi_{1,\downarrow}^{\dagger},\psi_{1,\downarrow},-\psi_{1,\uparrow},\psi_{2,\uparrow}^{\dagger},\psi_{2,\downarrow}^{\dagger},\psi_{2,\downarrow},-\psi_{2,\uparrow})$. Here, $\psi^{\dagger}_{j,\sigma}$ is a fermionic creation operator on wire $j$, with spin $\sigma$.

Finding zero-energy solutions of the above Hamiltonian density can be simplified as it is a block-diagonal matrix. In fact, it is sufficient to solve a $4\times 4$ matrix differential equation with real pairing strength and then restore phases $\exp(\pm i\phi/2)$ into the solution [$\Psi_{\phi}(x)=U\Psi_{\phi=0}(x)$] of Eq. (10) with the unitary transformation $U_{\phi}=\exp{\left(-\frac{i}{2}\frac{\phi}{2}\eta_{z}\tau_{z}\sigma_{0}\right)}$.Jiang et al. (2011) Thus, the relevant Hamiltonian density is

[TABLE]

in the basis $\Psi^{\dagger}=(\psi_{\uparrow}^{\dagger},\psi_{\downarrow}^{\dagger},\psi_{\downarrow},-\psi_{\uparrow})$ and we have dropped the wire index. We will work in the limit where the SOC energy is the largest energy scale of the problem and for simplicity, we assume $\mu=0$. As we use KwantGroth et al. (2014) for numerical calculations, where hard-wall boundary conditions are assumed, our analytical solution in semi-infinite geometry has to obey $\Psi(x=0)=0$ and $\Psi(x=\infty)\rightarrow 0$. Thus, we concentrate on the zero-energy state localized around $x=0$. The use of semi-infinite geometry to calculate this wave-function is a valid approximation as MBS at two ends of the system should have a negligible overlap in order to avoid hybridization.Klinovaja and Loss (2012)

Now, we perform a unitary transformation $U_{\rm{rot}}\mathcal{H}_{\rm{SNW}}U_{\rm{rot}}^{\dagger}$ where $U_{\rm{rot}}=\tau_{0}\frac{\sigma_{0}-i\sigma_{x}}{\sqrt{2}}$ to bring the SOC term into a diagonal form

[TABLE]

and

[TABLE]

The eigenenergies of the system are obtained using the Fourier transformation and are $E^{2}_{\pm}=\frac{\hbar^{4}}{4m^{2}}(k_{x}^{4}+4k_{\rm{SO}}^{2}k_{x}^{2})+\Delta^{2}+V_{z}^{2}\pm 2\sqrt{\Delta^{2}V_{z}^{2}+\frac{\hbar^{4}k_{x}^{4}}{4m^{2}}(V_{z}^{2}+\frac{\hbar^{4}}{m^{2}}k_{\rm{SO}}^{2}k_{x}^{2})}$, where $k_{\rm{SO}}=\frac{m\alpha}{\hbar^{2}}$ and $k_{x}$ denotes the momentum. The gap $E_{\rm{gap}}=2E_{-}$ is reduced for two values of $k_{x}$; at $k_{x}=0$ this gap equals $\Delta-V_{z}$ and its sign determines whether the system is in the topological phase $(V_{z}>\Delta)$ or not, while the gap at $k_{x}=k_{F}=2k_{\rm{SO}}$ is always nonzero as it is induced by superconductivity. Thus, all states close to $k_{x}=0$ form the interior branch while the states close to $k_{x}=k_{F}$ belong to the exterior branch of the spectrum.

In the regime of strong SOC, the magnetic field and proximity-induced superconductivity are considered as perturbations. To transfer the problem in the rotating frame, we use a spin-dependent gauge transformationKlinovaja and Loss (2012); Braunecker and Simon (2015) $[\psi^{\prime}_{\sigma}(x)=\exp(-i\sigma k_{\rm{SO}}x)\tilde{\psi}^{\prime}_{\sigma}(x)]$, where $\sigma=\pm 1$ denotes spin-up/spin-down.

Before the gauge transformation, kinetic energy and SOC of the electrons reads

[TABLE]

and after the transformation

[TABLE]

it becomes

[TABLE]

as the first-order derivatives in $x$ cancel exactly. Thus, this transformation eliminates the SOC term and now, the spectrum of $\tilde{\mathcal{H}}^{{}^{\prime}}_{\rm{T+SOC}}$ consists of two parabolas centered at $k_{x}=0$ $[E=\frac{\hbar^{2}}{2m}(k_{x}^{2}-k_{\rm{SO}}^{2})]$. Then, around the Fermi points $k=\pm k_{\rm SO}$, the spectrum can be linearized, i.e. the electron operator becomes

[TABLE]

where $\tilde{R}^{{}^{\prime}}_{\sigma}/\tilde{L}^{{}^{\prime}}_{\sigma}$ denote right/left movers.

In the linearized system, we ignore the second derivatives of $\tilde{R}^{{}^{\prime}}_{\sigma}/\tilde{L}^{{}^{\prime}}_{\sigma}$. Further, in the rotating frame, we also neglect the fast-oscillating terms $[\exp(\pm ink_{\rm{SO}}x)$, where $n>1]$. Then, the term $\tilde{H}^{{}^{\prime}}_{\rm{T+SOC}}$ becomes

[TABLE]

The spin-dependent gauge transformation also impacts the Zeeman field term, as it becomes a helical field that rotates in the plane that is perpendicular to the spin-orbit coupling vector in order to minimize the energy.Braunecker and Simon (2015) The form of this field $\mathbf{V}_{z}=V_{z}(\hat{y}\cos{2k_{\rm SO}x}+\hat{x}\sin{2k_{\rm SO}x})$, where $\hat{x}$/$\hat{y}$ are unit vectors in the x-/y-direction, is inherited from the Zeeman Hamiltonian

[TABLE]

where we have dropped the fast oscillating terms to obtain the last row.

The initial superconducting Hamiltonian was

[TABLE]

and after the transformation, it becomes

[TABLE]

Finally, we define two basis vectors $\phi^{{}^{\prime}i\dagger}=(\tilde{R}_{\uparrow}^{{}^{\prime}\dagger},\tilde{L}_{\downarrow}^{{}^{\prime}\dagger},\tilde{L}_{\downarrow}^{{}^{\prime}},-\tilde{R}_{\uparrow}^{{}^{\prime}})$ and $\phi^{{}^{\prime}e\dagger}=(\tilde{L}_{\uparrow}^{{}^{\prime}\dagger},\tilde{R}_{\downarrow}^{{}^{\prime}\dagger},\tilde{R}_{\downarrow}^{{}^{\prime}},-\tilde{L}_{\uparrow}^{{}^{\prime}})$. Then, we are able to construct Hamiltonians corresponding to interior ($k\sim 0$) and exterior ($k\sim 2k_{\rm{SO}}$) branches in the laboratory frame from terms in Eqs. (14), (15) and (16):

[TABLE]

where $l=i,e$ and

[TABLE]

The unitary transformation $U_{\rm{rot}}$ applied on Hamiltonians in Eq. (17) rotates the SOC/Zeeman energy into the form as in Eq. (11). Thus, these Hamiltonians, in the basis $\phi^{i\dagger}=(\tilde{R}_{\uparrow}^{\dagger},\tilde{L}_{\downarrow}^{\dagger},\tilde{L}_{\downarrow},-\tilde{R}_{\uparrow})/\phi^{e\dagger}=(\tilde{L}_{\uparrow}^{\dagger},\tilde{R}_{\downarrow}^{\dagger},\tilde{R}_{\downarrow},-\tilde{L}_{\uparrow})$, become

[TABLE]

Zero-energy solutions obey $H^{l}\Psi^{l}=0$ and to solve this set of equations, we assume the ansatz $\Psi^{l}=\exp(\kappa_{l}x)\psi^{l}$, where $\psi^{l}$ is a spinor. Then, the eigenenergies are

[TABLE]

For zero-energy modes, we obtain $\kappa_{i}=\pm\frac{\Delta\pm V_{z}}{\alpha}$ and $\kappa_{e}=\pm\frac{\Delta}{\alpha}$. As we look for a normalizable solution in the segment $[0,\infty)$, this leaves us with $\kappa_{i,1}=\frac{\Delta-V_{z}}{\alpha}$,$\kappa_{i,2}=-\frac{\Delta+V_{z}}{\alpha}$ as well as $\kappa_{e}=-\frac{\Delta}{\alpha}$. Then, to obtain these eigenfunctions in the basis $\tilde{\Psi}^{\dagger}=(\tilde{\psi}_{\uparrow}^{\dagger},\tilde{\psi}_{\downarrow}^{\dagger},\tilde{\psi}_{\downarrow},-\tilde{\psi}_{\uparrow})$, we add phases $\exp{(\pm ik_{\rm{SO}}x})$ associated with right/left movers to spinor parts of the ansatz $\Psi^{l}$. Finally, we obtain four eigenfunctions

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, coefficients $C_{i,1},C_{i,2},C_{e}$ are normalization coefficients that should be obtained from $\int_{0}^{\infty}dx|\Psi^{l}(x)|^{2}=1$.

The wave-function of a Majorana bound state should be a linear superposition of these four eigenfunctions. However, upon imposing the boundary condition $\tilde{\Psi}(x=0)=0$, we see that the only possible combination is $\tilde{\Psi}=\tilde{\Psi}_{1}^{i}-\tilde{\Psi}_{1}^{e}$. By using the relation $\Psi_{\sigma}(x)=\exp(-i\sigma k_{\rm{SO}}x)\tilde{\Psi}_{\sigma}(x)$, we obtain the wave-function in the laboratory frame

[TABLE]

where $C=\sqrt{\frac{2\Delta(V_{z}-\Delta)(V_{z}^{2}+4\frac{m^{2}\alpha^{4}}{\hbar^{4}})}{V_{z}\alpha[(V_{z}-2\Delta)^{2}+4\frac{m^{2}\alpha^{4}}{\hbar^{4}}]}}$. We find a good overlap of the analytical solution with the numerical one (see Fig. 8a) in the limit where SOC energy is the largest energy scale.

The wave-functions that are solutions of Eq. (10) for zero energy [$\Psi_{1}(x)/\Psi_{2}(x)$] are obtained after applying a unitary transformation $U_{\phi}=\exp{\left(-\frac{i}{2}\frac{\phi}{2}\eta_{z}\tau_{z}\sigma_{0}\right)}$ to $(1\;0)^{\dagger}\otimes\Psi^{\dagger}(x)/(0\;1)^{\dagger}\otimes\Psi^{\dagger}(x)$ and are equal to

[TABLE]

and

[TABLE]

Having two zero-energy states requires a degenerate perturbation theory to calculate the corrections to energy of these modes once the inter-wire coupling is present. Thus, we look at the matrix

[TABLE]

where $H_{\rm coupling}$ denotes the coupling between the nanowires.

We extract the coupling terms from Eq. (2), i.e. $H_{t_{y}}=-t_{y}\eta_{x}\tau_{z}$ and $H_{\beta}=\beta\eta_{y}\tau_{z}\sigma_{x}$. Thus, only the off-diagonal terms in Eq. (21) are non-zero and complex-conjugate to each other.

For the normal hopping, one of these off-diagonal terms is

[TABLE]

therefore, giving the matrix

[TABLE]

whose eigenvalues are $\pm t_{y}\sin{(\frac{\phi}{2})}$, which is also confirmed by a numerical simulation presented in Fig. 8b.

Once only the spin-orbit coupling in the y-direction is present, the off-diagonal term becomes

[TABLE]

and in the last line we have used the fact that SOC energy is the largest energy scale of the problem. This produces a Hamiltonian

[TABLE]

whose eigenenergies are $\pm\frac{\Delta}{V_{z}}\beta\cos{(\frac{\phi}{2})}$ and they accurately match numerical results presented in Fig. 8c.

Appendix B Constraints to the scattering matrix due to particle-hole symmetry

We have characterized our system using the determinant of the scattering matrix as the topological invariant. As the scattering matrix is unitary, the absolute value of its determinant has to be one. However, as a result of the constraint the charge-conjugation symmetry imposes on this matrix, this determinant is real and takes two opposite values across the phase transition.

The charge-conjugation symmetry is an anti-unitary symmetry that anti-commutes with the Hamiltonian, such that

[TABLE]

Thus, this symmetry can be formally written as $\mathcal{P}=\mathcal{U}\mathcal{K}$, where $\mathcal{U}$ is a unitary operator and $\mathcal{K}$ denotes the complex-conjugation. Furthermore, this symmetry also acts on the incoming/outgoing modes present in the scattering process. The incoming modes $\Psi_{n}^{\rm in}$ are all plane waves with velocity vector pointed into the scattering region, while this vector points outwards for the outgoing modes $\Psi_{n}^{\rm out}$. As the particle-hole symmetry flips energy as well as the momentum vector, the velocity defined as $\mathbf{v}=\frac{1}{\hbar}\frac{\partial\epsilon}{\partial\mathbf{k}}$ is unaffected and therefore

[TABLE]

where $V^{\mathcal{P}}$ and $Q^{\mathcal{P}}$ are unitary matrices and the summation over the repeated indices in assumed here and in the following.

The scattering matrix can be obtained by solving the Schrödinger equation

[TABLE]

where $\psi^{\rm loc}$ is a wave-function localized near the scattering region. We further simplify the proof by assuming $\Psi^{\rm loc}=0$, and for the particle-hole symmetric system, $\epsilon=0$. Then, by applying the charge-conjugation operator on the Eq. (23):

[TABLE]

Additionally, if $\Psi_{n}$ satisfies $H\Psi_{n}=0$ for any $n$, then the linear combination $HA_{nm}\Psi_{m}=0$. By taking $A=V^{-1}$ and inserting it into the last line of Eq. (24)

[TABLE]

Finally, by comparing Eqs. (23) and 25, the transformation law for the scattering matrix can be deduced as

[TABLE]

By transposing the last equation, one obtains

[TABLE]

Since, in our model $V^{\mathcal{P}}=Q^{\mathcal{P}}=\tau_{y}\sigma_{y}$, the scattering matrix obeys $S=\tau_{y}\sigma_{y}S^{*}\tau_{y}\sigma_{y}$.

Appendix C Effects of disorder

In this section, we study the effects of disorder in the chemical potential and phase difference on the system with the Zeeman field in the z-direction. As this field orientation produces a phase with only two gapless corner modes, we expect that the resulting phase diagrams calculated from right/left leads differ. Furthermore, the simulations of disorder effects on the system with only $V_{x}$ showed that different system configurations have opposite phase diagrams. Thus, here we study only the first setup (Fig. 3b), as we expect that the phase diagrams obtained from right/left leads of this setup are similar to the phase diagrams obtained from left/right leads of the other setup.

As before, fluctuations in the chemical potential are simulated with a variable $v_{\mu}$ that is uniformly distributed in the range $[-\chi_{\mu},\chi_{\mu}]$ and drawn independently for each lattice site. The phase diagrams calculated from left/right leads are presented in Figs. 9a and 9b, respectively. The x-edge conductance is shown in Fig. 9e, while the conductance along left vertical edge is given in Fig. 9f. Finally, the bulk conductance is plotted in Fig. 9g. Overall, the system behaves similarly to the case of $V_{x}$ only.

Next, we consider the effects of disorder in the phases of superconductors. This is simulated with a variable $v_{\phi}$ that is uniformly distributed in the range $[-\chi_{\phi},\chi_{\phi}]$ such that $\phi_{\rm{dis}}(y)=\phi(y)+v_{\phi}$. The resulting phase diagrams are presented in Figs. 9c and 9d. The only non-vanishing conductance (on the y-edge), calculated from the left leads, is presented in Fig. 9h. As for the other field orientation, disorder in $\phi$ can act detrimentally to the topological phase for one setup while the opposite happens for the other setup. Like before, in the region of small disorder strengths, the topological phase does not expand/reduce significantly.