Effects of gate-induced electric fields on semiconductor Majorana nanowires

We study the effect of gate-induced electric fields on the properties of semiconductor-superconductor hybrid nanowires which represent a promising platform for realizing topological superconductivity and Majorana zero modes. Using a self-consistent Schr\"odinger-Poisson approach that describes the semiconductor and the superconductor on equal footing, we are able to access the strong tunneling regime and identify the impact of an applied gate voltage on the coupling between semiconductor and superconductor. We discuss how physical parameters such as the induced superconducting gap and Land\'e g-factor in the semiconductor are modified by redistributing the density of states across the interface upon application of an external gate voltage. Finally, we map out the topological phase diagram as a function of magnetic field and gate voltage for InAs/Al nanowires.


I. INTRODUCTION
Composite heterostructures provide an opportunity to realize exotic phases of matter by exploiting the properties of individual components. A particularly interesting example involves semiconductor-superconductor (SM-SC) hybrid structures which represent a promising platform for the realization of topological superconductivity [1][2][3][4][5][6][7][8][9]. Topological superconductors support exotic neutral excitations consisting of an equal superposition of an electron and a hole-Majorana zero-energy modes (MZMs) [10][11][12]. Because of the particle-hole symmetry in a superconductor, such modes appear at zero energy and, thus, there is no cost to occupy these states. This leads to a growing degeneracy of the ground state as the number of MZMs is increased, a hallmark of topological superconductors. Theory predicts that exchanging the position of MZMs [10,13] or performing certain nonlocal measurements of the charge encoded in a collection of MZMs [14] leads to a nontrivial transformation within the degenerate ground-state manifold, and represents a non-Abelian operation which is independent of the details of its execution. This property of topological superconductors has generated a lot of excitement in the condensed matter physics, quantum information, and material science communities [15][16][17][18] as it opens up the possibility of Majorana-based topological quantum computing [6,9,19,20].
Realizing topological superconductivity in the laboratory is not an easy task since the originally proposed models [10,12] involved spinless p-wave superconductivity. Electrons in solids have spin 1=2 and most of the common superconductors have s-wave pairing which involves electrons with opposite spins. Therefore, quenching spin degeneracy and preserving superconducting pairing is quite nontrivial. One way to overcome the problem is to use materials with a strong spin-orbit interaction which couples spin and orbital degrees of freedom. A number of platforms for realizing MZMs in the laboratory have been recently proposed . The most promising proposal for realizing MZMs is based on one-dimensional SM-SC hybrid structures [27,28] and involves a semiconductor with strong spin-orbit coupling (such as InAs or InSb) and an s-wave superconductor (such as Al). In this proposal, a magnetic field or another time-reversal breaking perturbation is needed to drive the system into the spinless topological regime [27,28]. This proposal has triggered significant experimental activity , and there is a compelling body of experimental evidence that MZMs have been realized in these systems. For a very recent example, see Ref. [64], which reports a robust quantized 2e 2 =h zero-bias conductance consistent with the Majorana scenario.
Much of the progress in realizing MZMs with proximitized nanowires is attributed to the material science advance in fabricating semiconductor-superconductor heterostructures. In the first generation of experiments [49][50][51][52][53][54] the superconductor was deposited ex situ, which required removing the native oxide forming on the semiconductor's surface due to air exposure. In the second generation of experiments the thin aluminum shell [56] is deposited epitaxially and is thus grown on pristine SM facets without breaking the vacuum; see Fig. 1. Tunneling spectroscopy measurements of the induced superconducting gap [55,60,62,64,68] in such samples exhibit a large induced gap (i.e., close to the bulk gap of the superconductor), which indicates that the improved epitaxial interfaces are characterized by a strong hybridization of the states in the semiconductor and superconductor. In this strong tunneling regime, many physical parameters such as the g factor and spin-orbit coupling are strongly renormalized due to the hybridization. In order to quantitatively understand the hybridization and its implications on the band structure as well as other physical properties, one has to consider the band offset at the superconductor-semiconductor interface. Depending on the sign of the band offset, one can have either a Schottky barrier or an accumulation layer [70][71][72][73]. Based on preliminary angle-resolved photoemission spectroscopy studies [74], one finds that the band offset for epitaxially grown InAs=Al heterostructures is −ð200-300Þ meV, supporting the accumulation layer scenario.
Proper theoretical treatment of the strong-coupling regime is also necessary to understand how external gates affect the electronic state, and in particular the topological nature, of SM-SC heterostructures. Furthermore, recent proposals for realizing scalable architectures for topological quantum computation with MZMs rely on fine electrostatic control [75][76][77][78][79]. Thus, understanding the effect of electric fields on the low-energy properties of the proximitized nanowires is critical both for the interpretation of the existing Majorana experiments [57,60,62,64,68] as well as for the optimization of proposed Majorana devices [9].
In order to understand the physical properties of the proximitized nanowires, one needs to solve the electrostatic and quantum-mechanical problems self-consistently, i.e., perform Schrödinger-Poisson (SP) calculations. Compared to the case of purely semiconducting heterostructures [80][81][82], the problem at hand is much more challenging technically because it involves disparate materials with very different effective masses, Fermi energies, g factors, etc. (see Table I). In other words, the standard numerical tools based on the continuum mass approximation cannot be applied to semiconductor-superconductor hybrid systems. Therefore, modeling of the semiconductor-superconductor hybrid structures requires developing numerical techniques which can effectively take into account different length scales in the semiconductor and superconductor.
Previous effective models for superconductorsemiconductor hybrids [87][88][89][90][91][92] do not properly describe the experimental system and provide only qualitative predictions for the electric field dependence. These models rely on independent phenomenological parameters such effective masses, spin-orbit couplings, g factors, as well as tunneling strength between semiconductor and superconductor. While this approach may be suitable for the weak tunneling regime, naive extensions of such models to the strong-coupling limit are inadequate. This is because the electric field applied to the semiconductor can drastically change the electrons' confinement, i.e., push or pull electron density in the semiconductor to or away from the interface. This in turn strongly affects physical parameters of the system, including, as we will see, the tunneling rate, effective spin-orbit coupling, g factor, as well as induced superconducting gap.
More advanced models have been introduced recently [93][94][95][96][97] which treat the effects of an electric field within some effective models where the superconductor is taken into account via boundary conditions. This approach, while being computationally advantageous, does not take into account the effects arising from the redistribution of the wave function between the semiconductor and the superconductor. In this work, we treat the superconducting and semiconducting degrees of freedom explicitly on the same footing. Using an adaptive discretization algorithm for the SM and SC components, we develop an effective model which is computationally tractable and allows us to adequately capture the effect of the gate-induced electric field on the heterostructure.  L z (nm) 60 10 L y (nm) 52 52 Our results allow one to understand and interpret recent experiments investigating the electric field and disorder dependence of the effective parameters [57,60,62,64,68]. They have motivated the recent systematic experimental study of the tunability of the superconductor-semiconductor coupling in Majorana nanowires [98]. A better understanding of the effect of the external electric fields on the phase diagram of superconductorsemiconductor nanowires is important to reinforce the confidence that the observed zero-bias peaks are due to Majorana zero modes. Our results describe the conditions necessary to achieve a strong tunneling regime between superconductor and semiconductor and show how external electric fields and the presence of disorder affect the topological phase diagram of SC-SM nanowires. They provide an important contribution for the realization of the next-generation experiments designed to verify the non-Abelian nature of the modes responsible for the zero-bias peak observed so far in transport measurements and to give an indisputable proof of the realization of Majoranas.
The paper is organized as follows. We begin with a discussion of the setup and methods in Sec. II, where we provide technical details of the Schrödinger-Poisson approach. In Sec. III we present our results. We first focus on the limit of zero magnetic field and then discuss the behavior at finite magnetic fields. We then discuss the resulting topological phase diagram. We conclude the section with studying the effects of disorder that change the strength of SM-SC coupling. We summarize our results in Sec. IV and discuss their relevance for current and future experiments.

II. SETUP AND METHODS
We consider the system shown in Fig. 1. The nanowires used in current experiments typically have a hexagonal shape as shown in Fig. 1(a). The cross section of the wire, which we take to be the ðy; zÞ plane, consists of a 10-nmthick Al film (blue) covering two facets of InAs nanowire (orange). The electrostatic environment is controlled by a back gate (gray). For practical reasons we do not explicitly treat this gate and the separating dielectric medium in our calculations, but rather take the gate into account only as a boundary condition for the potential in the wire. In order to convert this into the actual voltage applied to the gate (which is sample dependent), the distance to the gate and the dielectric constant have to be taken into account. For the devices of interest, the length of the wire L x is much larger than its transverse dimensions L y;z .
The presence of the Al layer breaks the hexagonal symmetry of the nanowire cross section and, as shown in Fig. 2, causes the formation of an electrostatic potential that strongly confines the electrons close to the SM/Al interface. For this reason the hexagonal cross section of the wire can be well approximated by an effective rectangular cross section, as shown in Fig. 1(b). We henceforth refer to the effective wire with rectangular cross section as the slab model. By choosing L y for the slab model to be such that the number of cross-sectional modes is the same as for the hexagonal cross section wire, the use of the slab model does not cause any significant loss of accuracy and significantly simplifies the numerical implementation and solution of the SP problem.
The Hamiltonian for the heterostructure in the normal state can be written as (ℏ ¼ 1) where the spatially dependent effective mass m Ã ðzÞ, Fermi energy ε F ðzÞ, spin-orbit coupling strength αðzÞ, and g factor gðzÞ are equal to m Ã ðzÞ ¼ m SM [m Ã ðzÞ ¼ m SC ] for z < 60 nm (z > 60 nm) and similarly for ε F ðzÞ, αðzÞ, and gðzÞ; k x ,k y are the momentum operators in the x and y direction, respectively; ϕðzÞ is the electrostatic potential; σ x;y;z are the Pauli matrices in spin space; μ B and B are the Bohr magneton and the external magnetic field, respectively. The values for the material parameters used henceforth are given in Table I. In this work we investigate bulk properties of the heterostructure. Therefore, we assume henceforth that the nanowire is infinitely long and translationally invariant along the x direction. This allows one to use the basis of plane waves along the x direction and to replace the operator k x in Eq. (1) by its eigenvalue. In the clean limit considered here, due to the finite-size quantization in the y and z directions, the spectrum of the system consists of effectively 1D subbands. We obtain the eigenvalues and eigenstates of FIG. 2. The electrostatic calculation uses Dirichlet boundary conditions at z ¼ 0 and z ¼ L SM z , i.e., the top and bottom of the semiconducting wire. At z ¼ 0, the boundary condition is given by the gate voltage V g , while at the interface to the aluminum it is given by the band offset W (see also Table I). This leads to an accumulation layer at the interface. the resulting HamiltonianĤ n ðk x → k x Þ corresponding to these subbands via a mode decomposition in the y direction and by replacing the derivatives with respect to z with finite differences using a nonuniform grid [99] with two different spacings corresponding to the semiconducting and superconducting components, respectively. The spacings are chosen such that dz < π=k F in order to minimize discretization errors. Using a nonuniform spacing significantly alleviates the computational cost and allows us to systematically study the phase diagram of the problem.
In the absence of spin-orbit coupling and disorder, the discrete modes along the y direction are ψ α¼0 n y ðyÞ ¼ with the different n y ∈ N modes being decoupled. The spin-orbit coupling term hybridizes them [87]. The corresponding matrix elements are In the fψ k x ;n y ;z ¼ e ik x x Ψ α¼0 n y ðyÞδðzÞg basis the Hamiltonian matrix takes the form where A is the N y × N y matrix with elements given by Eq. (3), N y is the number of discrete modes along the y direction considered, and H α¼0 n is the matrix obtained by projecting the Hamiltonian Eq. (1) for α ¼ 0 on this basis.
We treat the s-wave superconductor at the BCS meanfield level. The Bogoliubov-de Gennes (BdG) Hamiltonian for the system can be written as where τ 0;x;y;z are Pauli matrices in Nambu (particle-hole) space. We include the superconducting pairing only in the superconductor, i.e., ΔðzÞ ¼ Δ 0 for z > 60 nm, where Δ 0 is the SC gap of Al (see Table I), and ΔðzÞ ¼ 0 for z < 60 nm. In a finite magnetic field, the superconducting gap in the Al shell is suppressed due to the inclusion of a finite g factor for the Al (see Table I). Given that the Al film is very thin, see Fig. 1, we neglect orbital effects due to the magnetic field.

A. Electrostatics
In order to obtain the electrostatic profile, one has to solve the SP equations self-consistently. Given that the BCS meanfield approximation breaks charge conservation, this is a nontrivial task; see, e.g., discussion in Ref. [100]. However, electrostatic screening of a metal is only weakly modified by the superconductivity with the small parameter being Δ 0 =ϵ F ≪ 1. As a consequence, to obtain the electrostatic potential within this accuracy the charge density entering the Poisson equation can be calculated neglecting the superconducting pairing, i.e., using the Hamiltonian H n instead of the full Hamiltonian H. The effects of the spin-orbit coupling and Zeeman terms [94] on the total electron density profile nðzÞ are also very small and can be neglected. Thus, to solve the full problem we follow a two-step approach. We first solve the SP problem in the normal state taking α ¼ 0 and B ¼ 0 to obtain the electrostatic profile. We then use the obtained electrostatic profile to find the eigenvalues and the eigenstates of the system for Δ 0 ≠ 0, α ≠ 0, and different values of B.
The first step consists in solving self-consistently the Schrödinger equation H n jΨi ¼ EjΨi, requiring Ψ to vanish at the boundaries of the system [101], and the Poisson equation, where nðzÞ ¼ ð2πL y Þ −1 P n y ;E½Ψ<¼0 R dk x jΨ k x ;n y ðzÞj 2 , ϵ r is the relative dielectric constant of the SM, see Table I, and ϵ 0 is the vacuum dielectric constant. The setup for the Poisson equation is shown in Fig. 2. At z ¼ L SM z the boundary condition for ϕðzÞ is given by the band offset W between the SM and the SC. The boundary condition at z ¼ 0 is set by the back gate. The coupled Schrödinger-Poisson equations are solved iteratively until convergence is achieved, using Anderson's mixing algorithm [102].

B. Band structure
The calculated electrostatic profile ϕðzÞ is inserted into the full Hamiltonian H to obtain the band structure fε ðnÞ ðk x Þg and the corresponding eigenstates of the nanowire. Since the chemical potential is included in the Hamiltonian, the effective Fermi energy for each band is set simply by the bottom of the band. We can find the Fermi momentum in each band k ðnÞ F by solving ε ðnÞ ðk ðnÞ F Þ ¼ 0. The Fermi velocity is given by v n F ¼ ½ðdε ðnÞ Þ=ðdk ðnÞ Þj k ðnÞ ¼k ðnÞ F . In addition, from the eigenstates at k ¼ k F we extract how strongly different subbands are coupled to the superconductor, which we define through the weight of the corresponding state in the superconductor: We define the gap as the minimum of the energy of the first excited state: BdG ðk x Þj. At zero magnetic field, E g gives an estimate for the induced gap Δ ¼ E g ðB ¼ 0Þ.

III. RESULTS
In this section, we discuss the results of our numerical simulations. We first discuss only the electrostatic problem for both a model of a hexagonal wire and the slab model introduced above. We then investigate the nature of the electronic states in a limit of strong coupling between the semiconductor and superconductor and discuss their superconducting properties at zero magnetic field. Then we study properties of the hybrid nanowires in a finite magnetic field and obtain estimates for the effective g factor in the hybrid structure. We present the topological phase diagram and compare it with previous results [87,103]. Finally, we present the results for the wires with the disorder potential present in the superconductor and show its impact on the induced gap and the phase diagram.

A. Electrostatics and density distribution 1. Hexagonal cross section
In order to obtain the correct number of subbands for a given gate configuration for the wire with the hexagonal cross section, it is sufficient to solve the SP problem using the Thomas-Fermi approximation and simply requiring the wave function to vanish at the boundaries of the cross section. The solution of the full SP problem is computationally expensive due to the shape of the cross section and unnecessary for the purpose of simply estimating the number of cross-sectional modes. We perform this calculation in COMSOL and obtain eigenstates using the KWANT package [104].
Our results are summarized in Fig. 3, where we show the density for all occupied modes below the Fermi energy. This calculation does not explicitly treat the aluminum shell; instead, it assumes that the only effect of the presence of the Al layer is to induce a band offset. We set this band offset to W ¼ −0.25 eV [74]; see Table I. The approximations used to obtain the results of Fig. 3 cause quantitative inaccuracies for the local density of states (LDOS) and the carrier density profile. However, these results are sufficiently accurate to estimate the number of electronic cross-sectional modes below the Fermi energy for a given V g . In addition, the results of Fig. 3(a) show the qualitatively correct result that for V g ≤ 0 most of the charge density is localized at the semiconductor-superconductor interface due to the strong band offset between the InAs and Al. This fact means that for the slab model, the thickness of the SM wire in the z direction does not affect the electronic properties in a significant way as long as it is few times larger than the confinement length in the z direction (∼20 nm). The effective width L y of the slab model can then be fixed by requiring the number of subbands to be equal to the number of cross-sectional modes obtained from the hexagonal calculation, as long as L y is also larger than the confinement length in the z direction. For V g ¼ 0, the hexagonal cross section results show that there are six modes; see Fig. 3(b). From this we obtain that for the slab model, L y ¼ 52 nm, larger than the confinement length for V g ¼ 0.
In the remainder, all of the results are obtained using the effective slab model with L y ¼ 52 nm width and L SM z ¼ 60 nm thickness for the SM and L ðSCÞ z ¼ 10 nm for Al, as shown in Table I.

Slab model
We now switch to the slab model, which explicitly treats the superconducting Al shell. We self-consistently solve the coupled Schrödinger-Poisson equations for three different values of V g to obtain the electrostatic potential ϕðzÞ and the density nðzÞ, respectively, shown in Figs. 4(a) and 4(b). Since the Al shell is taken to be metallic with an extremely short screening length, the electrostatic potential is assumed to be constant throughout the Al. The dashed  line in Fig. 4(a) shows the Fermi level in Al. It is worth pointing out that because Δ 0 ≪ ϵ F , including the pairing term for the Al makes only a negligible difference to the electrostatic profile. For V g ≤ 0, the electrostatic potential confines the carrier density in a layer about 20 nm wide close to the SM/Al interface, as shown in Fig. 4(b). For V g > 0, the electrostatic potential is below the Fermi energy also on the gate side. This allows the accumulation of charges also near the gate, as shown by the result in Fig. 4 B. Nature of electronic states in strong-coupling limit We now discuss the nature of the electronic states in the electrostatic environment determined by the gate as well as the band offset between the semiconducting wire and the metallic shell. In particular, we investigate how strongly states are hybridized between the two materials depending on the gate voltage.
In Figs. 5(a)-5(c), we show the electrostatic profile (cf. Fig. 4) for three values of the gate voltage. For the case of V g ¼ 0 [Figs. 5(b) and 5(e)], we find nine hybridized subbands, some of which are mostly localized in the SM, whereas the others have large weight in the superconductor [105].
For V g < 0, the electrostatic potential confines the wave function in the SM to a very narrow region close to the SM/Al interface. Such confinement favors a strong hybridization of the SM and Al eigenstates, thus giving rise to states which have large weight in both the SM and Al. Such large hybridization is prevented in the absence of the confining electrostatic potential due to the large mismatch between the Fermi velocities of the two materials. The strong confining potential due to the band offset is therefore critical for the hybridization of the SM and Al states.
For V g > 0, we see in Fig. 5(c) that a number of subbands closer to the Fermi energy appear which are not confined to the interface, and instead have appreciable weight throughout the SM. These states have very small W SC . Their contribution to the density can also be seen in Fig. 4(b) in the peak of the density near the gate.
While one might naively expect that the lowest bands are most confined to the interface and thus hybridize most strongly, this is not reflected in the data shown in Fig. 5. To further elucidate which bands most strongly hybridize with the superconductor, we show the full band structure at Fig. 6. Here, color again indicates the weight of the state in the SC; however, in contrast to Fig. 5, we do not just consider k x ¼ 0. We observe that hybridization with the superconductor may depend strongly on k x , and in this case is generally strongest at large k x .

C. Superconducting properties at B = 0
We now turn our attention to the situation where the Al shell is in the superconducting state. The value of W SC at k ðnÞ F correlates well with the magnitude of the induced superconducting gap Δ ind for a given subband. From the discussion of the previous section, we can then immediately conclude that different subbands will have different values of Δ ind . This is illustrated in the lower panel of Fig. 6, in which the subbands are shown for the case when Δ 0 ≠ 0, for energies of the order of Δ 0 . We see that bands with smaller W SC (shown as more blue) also have smaller Δ ind . The smallest value of Δ ind is what fixes the superconducting gap for the SM-SC hybrid nanowire. This again emphasizes the importance of the strong confining potential to increase the hybridization between the two materials and thus a large induced gap. For V g < −0.6 V, all subbands in the SM hybridize very strongly with the Al subbands and, thus, have a large induced SC gap. As a result, there are no subgap states below ε ≈ 0.6Δ 0 . As V g increases and the electrostatic potential becomes less confining, additional SM subbands become occupied for certain threshold values of V g . As shown in Figure 7(a), the number of subbands jumps at V g ≈ ð−0.6; −0.45; −0.15; −0.06; 0.03Þ V. In some cases the additional subbands have a smaller value of W SC resulting in a decrease of Δ ind . From Fig. 7(b) we can see that this happens for the V g threshold values of −0.6 and −0.06 V. For V g > 0, as shown in Fig. 5(d), some of the subbands have states that are not localized close to the SM-SC interface and for which W SC is negligible. In this situation Δ ind → 0, and the system becomes gapless.
The evolution of Δ ind with V g is shown in Fig. 7(b), together with the evolution of W SC . From this figure we see that for V g < −0.6 V, Δ ind ≈ 0.75Δ 0 . Furthermore, these results indicate that the evolution of the nanowire's superconducting gap with V g can be quite nontrivial and is closely related to W SC . In order to have Δ ind ∼ Δ 0 , strong confining potentials (V g < −0.6 V) are necessary. Conversely, in the case of a positive gate voltage, there are occupied states in the SM [Figs. 5(c)-5(f)] which are far away from the SC and, as a result, are weakly proximitized. Figures 7(c)-7(f) show the evolution of ϵ F , k F , v F , and ξ with V g for the subband with the smallest induced superconducting gap, which determines Δ ind for the system. For a fixed number of subbands, as V g increases ϵ F , k F , and v F grow; see Figs. 7(c)-7(e). Using the values of v F and Δ ind one can estimate the coherence length ξ ¼ ℏv F =ðπΔ ind Þ. From Fig. 7(b), we see that change in V g preserving the number of occupied subbands leads to small changes in Δ ind . Thus, the variations of ξ are mostly due to the changes in v F ; see Fig. 7(f). We see that, as long as the number of subbands is constant, ξ grows with V g and follows v F . The discontinuities in ϵ F , k F , v F , and ξ appear when the number of occupied subbands changes; see Figs. 7(a)-7(f).

D. Superconducting properties at finite magnetic fields
We now study how the properties of the SM-SC nanowire depend on the presence of the external magnetic field B aligned along the longitudinal direction of the wire. As discussed in Sec. II, in our treatment the magnetic field enters only via the Zeeman term. For B ∼ 1 T, the orbital effect of the applied magnetic field is small since the SC is only 10 nm thick and in the regime of interest the wave functions in the SM are confined to the SM-SC interface within a 20 nm range.
We start by investigating Zeeman splitting for the nanowire with multisubband occupancy. The corresponding band structure is shown in Fig. 8(a). Let us consider the gate voltage such that the highest-occupied subband (shown in blue) has small Fermi energy. The application of a magnetic field splits the subband and, at some critical field B c , drives the minority subband across the Fermi level (provided B c is less than the critical field of the superconductor). This is illustrated in Fig. 8(b). At this point, the majority subband becomes the highest-occupied band, and, thus, many properties such as the Fermi energy, Fermi velocity, and Fermi momentum change discontinuously.
Another scenario corresponds to Fig. 8(c), where a band is just above the chemical potential. In this case, the gap of the system is determined by the lower-occupied subband (shown in green). An increasing magnetic field splits the lowest-unoccupied band (shown in red) and eventually it becomes occupied.
In both these cases, we end up with an odd number of occupied subbands at large enough magnetic fields and, thus, the nanowire can be driven into the topological regime. However, the evolution of the gap with the magnetic field is drastically different in these two cases. This can seen in Fig. 9, which shows the evolution with magnetic field of the spectral gap E g , effective Fermi energy for the highest-occupied subband, and corresponding k F , v F , and ξ for different gate voltages close to the threshold value V g;t ¼ −0.427 V. This threshold value, corresponding to the gray curves in Fig. 9, corresponds to the gate voltage at which the relevant subband has exactly zero effective Fermi energy. As B increases from zero the gap E g decreases and eventually vanishes at B ¼ B c , corresponding to the topological phase transition. For B > B c , the nanowire is driven into the spinless regime with p-wave pairing potential. The p-wave gap exhibits a nonmonotonic dependence on the magnetic field, and eventually vanishes because the s-wave superconductor becomes normal. For our parameters this occurs at B SC ¼ 5.8 T. Note that we do not take into account orbital effects here, so in practice Δ 0 may vanish before that.
The blue family of curves in Fig. 9(a) correspond to the case when a highest-occupied subband has a small Fermi energy at zero field [Figs. 8(a) and 8(b)]. In that case, the gap at zero field is already set by the band that will eventually be split to give rise to topological superconductivity, and, thus, the gap evolves as a smooth function for B < B c . At B ¼ B c , the minority subband crosses the Fermi level, and the topological gap is opened in the majority subband. As discussed above, the properties of the Fermi points evolve discontinuously across the transition (panels b,c,d in Fig. 9), and the gap increases rapidly into the topological phase (panel a in Fig. 9).
At more negative gate voltages, the situation shown in Figs. 8(c) and 8(d) is realized, corresponding to the red lines in Fig. 9. Here, the gap at B ¼ 0 is determined by the next occupied subband. Upon the application of a magnetic field, the distance between the majority subband and the chemical potential eventually becomes smaller than the gap induced in the next-highest subband. This distance thus sets the spectral gap. The discontinuity of the gap function can be seen in Fig. 9(a). At B > B c , the topological gap is opened in the majority subband. In this case, we plot in Fig. 9 the properties of the Fermi points only for the subbands that eventually become topological, and thus plot no values below the topological phase transition.
It is very interesting to notice that the size of the induced superconducting gap for B ¼ 0 does not necessarily correlate with the size of the topological gap. This can be understood from the fact that the topological gap for B > B c is always opened in the same band, whereas the Δ ind at B ¼ 0 is opened in a different band when V g becomes smaller than V g ¼ −0.427 V. As can be seen from Fig. 9(c), the Fermi momentum k F for B > B c , which corresponds always to the same band, increases with V g . The topological gap increases with k F since the effective Rashba field is stronger at higher momentum, allowing the s-wave pairing to induce a larger gap. This is a very important result because Δ ind at B ≳ B c sets a crucial scale for the robustness of a topological qubit against error sources such as thermal fluctuations, diabatic corrections, disorder [100,[106][107][108], etc.

E. Effective g factor
A crucial quantity to characterize the semiconductorsuperconductor system is the effective g factor of the hybrid system. Because of the drastically different g factors in the two materials, this will depend intricately on the wave function hybridization between them. Furthermore, the g factor is crucial in enabling a large and robust topological phase, since a large g factor is necessary for the topological phase transition to occur at a magnetic field well below the critical value at which the Al shell is driven normal. A large effective g is thus very helpful in achieving a sufficient separation between these scales.
The effective g factor can be obtained from studying the Zeeman splitting of bands at k x ¼ 0, as illustrated in Fig. 8. In particular, since at k x ¼ 0 the spin-orbit terms in the Hamiltonian Eq. (4) vanish, the spin splitting of the bands at k x ¼ 0 is entirely determined by the Zeeman term. As the change of the energy levels εðk x ¼ 0Þ is linear with the magnetic field, the absolute value of the g factor can be extracted as jgj ¼ 2f½dεðk x ¼ 0Þ=ðdμ B BÞg. This linear fit for B < B c is illustrated in Fig. 9(b) with dashed lines. Note that when the gate voltage is such that the closest subband to the Fermi level is unoccupied [ Fig. 8(c)], the slope near the gap closing is the same as the one at εðk x ¼ 0Þ, allowing for a reliable extraction of the g factor from the tunneling conductance measurements [68]. This is shown in Fig. 9(a) by the dashed lines.
In Fig. 10 we study the dependence of the extracted g factor on the applied gate voltage in the vicinity of the threshold values at which the number of subbands change (see Fig. 7). As expected, we find every subband to be characterized by an almost constant g factor, with significant changes occurring only at transitions between bands. When the hybridization between semi-and superconductor is weak, the g factor is close to the bare semiconductor value jg bare SM j ¼ 15. Conversely, when the voltage is very negative (the value from Fig. 7 is written in every panel) and the hybridization between semi-and superconductor is strong, the g factor is almost as small as the bare superconducting g factor jg bare SM j ¼ 2. Additional information can be extracted from the ratio of the induced gap to the critical field, shown in Figs. 10(b), 10(d), 10(f), 10(h), and 10(j). This quantity is easily accessible in experiments, and has been used in the experimental literature as a proxy for the g factor [68].
Our results clearly show that unlike the g factor, this quantity has a strong dependence on gate voltage over relatively small gate voltage variations. In particular, a resonant structure appears with a peak that corresponds to the gate voltage being tuned to the threshold value at which the subband crosses the effective chemical potential. Only at this point does this quantity reach the values of the effective g factor shown in Figs. 10(a), 10(c), 10(e), 10(g), and 10(i). Figure 11(a) shows the value of g for different topological regions (see also the discussion in Sec. III F). We see that as V g becomes more negative the g factor becomes smaller and approaches the value of g in the SC. As stated above, this is due to the fact that as V g becomes more negative the hybridization between SM and SC states becomes stronger, as clearly shown by the evolution of W SC ; see Fig. 11(b). Larger negative values of V g create an electrostatic potential that more strongly confines the SM states at the SM-SC interface. The tighter confinement results in a stronger hybridization between SM and SC states. Figure 11(c) summarizes the important relation between strength of the hybridization between SM and SC states and the g factor by showing the dependence of jgj on W SC . We see that qualitatively g scales linearly with W SC .  external magnetic field-rather than more abstract quantities such as the Fermi energy of the subbands and the Zeeman splitting, which are dependent on applied electric field. As discussed above, the relation between V g and the parameters characterizing the nanowire band structure, such as the subbands' Fermi energy and induced superconducting gap, is highly nontrivial given the nonlinear nature of the SP problem and the presence of multiple subbands. For this reason, simplified models in which the subband chemical potential is assumed to be directly proportional to V g in general cannot be used to obtain a reliable phase diagram in the ðV g ; BÞ plane. Similarly, we have shown that for the g factor we cannot take the bare value of the SM. One qualitative feature that emerges from the results shown in Fig. 12(a) is how the shape and size of the topological regions depend on V g . We see that for very large negative V g the critical magnetic field is higher than that for small negative values of V g . This is due to the fact that the hybridization of the SM's and SC's states is stronger for larger negative V g and therefore the effective magnitude of g is much smaller than g bare SM , causing an increase of the critical field.

F. Topological phase diagram
The topological phase diagram can be obtained by calculating the topological index M ∈ Z 2 (Majorana number) [12,87,103]: where Bðk x Þ is an antisymmetric matrix which defines the Hamiltonian of the system in the Majorana basis. The negative (positive) sign of M corresponds to a topologically trivial (nontrivial) phase. The latter supports Majorana zero modes at the ends of the nanowire. In the continuum limit, the lattice spacing a → 0, and the sign of Pf½Bðk x → ∞Þ is fixed. Thus, the topological quantum phase transition corresponds to a change of sgnfPf½Bðk x ¼ 0Þg. Note that the topological phase transition in this case is accompanied by a vanishing quasiparticle gap at k x ¼ 0, which provides another way of determining the phase boundary.
It is illuminating to compare the phase diagram of Fig. 12(a) with previous studies [87,103]. Adapting the results of Refs. [27,28,103], the critical magnetic field for the topological transition is given by where ϵðk x ¼ 0Þ defines the position of the band bottom at k x ¼ 0 relative to the Fermi energy in the superconductor and g is the effective g factor. As follows from the discussion above, the effective g factor can be obtained by calculating Zeeman splitting at k x ¼ 0.
Having obtained the dependence of ϵðk x ¼ 0Þ, g, and Δ on V g , one can draw the boundaries in the ðV g ; BÞ plane of the topological phase diagram using Eq. (9). These boundaries are shown in green in Figs. 13(a), 13(c), and 13(e). We see that they exactly match the boundaries obtained by identifying the value of B, B c , where the gap is closing, Δ ind ¼ 0. In particular, Eq. (9) gives the correct boundaries if the renormalization of the g factor is taken into account. On the other hand, if in Eq. (9) we use the bare SM g factor, Eq. (9) gives incorrect boundaries, shown in red in Figs. 13(a), 13(c), and 13(e). The boundaries obtained using the bare SM g factor overestimate the size of the topological region, especially when V g is strongly negative, as shown in Fig. 13(e). As discussed above, this is due to the fact that the value of the g factor, in the strong-coupling regime, is strongly renormalized by the hybridization between the SM and SC states. Figures 13(b), 13(d), and 13(f) show the band structure of the SM-SC nanowire close to the Fermi energy when Δ 0 → 0 for the appropriate values of V g . We see that for very negative values of V g , Fig. 13(f), the SM states are very strongly hybridized with the SC states. This result is consistent with the fact that for this case g is much smaller than g SM and, therefore, the topological region is much smaller that we would have obtained assuming g ¼ g SM .
Figures 12(b)-12(d) show the dependence of the gap [ Fig. 12(b)], Fermi velocity [ Fig. 12(c)], and coherence length [ Fig. 12(d)] on the magnetic field, for B > B c , for different domes in the topological phase diagram shown in Fig. 12(a). To obtain these figures we use the representative gate voltages indicated by the white dashed lines in Fig. 12(a). For fixed magnetic field, the gap decreases with the gate voltage while the Fermi velocity increases. As a consequence, the coherence length increases. The results of Fig. 12(d) show that when the topological gap is maximal, ξ has values in the 100-300 nm range, and that ξ grows linearly with B for B > B c . The growth of ξ with B is slower for more negative gate voltages, a reflection of the fact that for larger negative gate voltages the effective g factor is smaller due to the stronger hybridization of the SM's and SC's states. The results of Fig. 12(d) are important for the design of Majorana-based qubits since their topological protection relies on the exponentially small splitting of Majorana zero modes which, in turn, strongly depends on the coherence length.

G. Effect of disorder in SC-SM heterostructures
Disorder is ubiquitous in solid-state systems and has a strong impact on the physical properties of proximitized nanowires. It has been shown that disorder in the SM and at the SM-SC interface is detrimental to the topological phase [109][110][111][112][113][114][115][116][117][118]. However, in the MBE-grown SC-SM heterostructures [9,56] disorder effects are minimized resulting in high-quality SMs as well as SM-SC interfaces. The SC (i.e., Al) is also nominally of high quality but its outer surface is covered by an amorphous oxide layer; see Fig. 14. Therefore, the scattering from the outer boundary randomizes the motion of quasiparticles in the SC. Because of the large effective mass mismatch and the conservation of the momentum at SM-SC interface, disorder in the SC is not expected to interfere with the observation and manipulation of MZMs [92,112,118]. Nevertheless, as we show it can strongly affect key quantities of the proximitized nanowire, such as Δ ind and the critical field B c , as well as their dependence on the external gate voltage.
We model disorder by adding the random potential V D to the Hamiltonian Eq. (1). Because of the computational complexity of the problem, we consider here the disorder potential V D ðz; yÞ (i.e., homogeneous along the wire). Such potential hybridizes different y and z subbands, which effectively increases density of states in the SC. In contrast to Eq. (4) the modes in the y direction are not separable anymore, and we have to resort to a numerical solution of the full Hamiltonian, leading to a much greater numerical complexity. We calculate results for a given disorder realization and then average physical observables over approximately 25 disorder realizations.
We model the amorphous oxide layer in the superconductor by adding a disorder potential within l z ¼ 2 nm from the outer surface of the SC; see  the disorder potential to have zero average and to be spatially uncorrelated: Here, ⟪ Á Á Á ⟫ denotes averaging over disorder realizations, K 0 parametrizes the disorder strength, and n imp corresponds to the density of impurities. Considering the finite spatial resolution in the numerical calculation and the uniform box distribution of the disorder with the amplitude U D , K 0 is related to U D as δrU 2 D =3 ¼ K 0 n imp , where δr is the volume of a single cell of the spatial discretization (see Table I), that we set to be uniform in the region where the disorder is located. We vary U D in the parameter range 0-1 eV, which corresponds to the effective mean-free length larger than 10 nm.
One of the main effects of the disorder in the SC is to break the conservation of the momentum and to induce broadening of the SC subbands, which effectively increases the number of superconducting subbands hybridizing with a given SM mode. This effect is shown in Fig. 14(b) where we plot the wave function probability density: the wave function probability in the SC is random, which corresponds to chaotic motion of quasiparticles, whereas the one in the SM preserves periodicity in the y direction. One may notice that disorder leads to an enhancement of W SC , as shown in Fig. 15(a). As U D increases, more of the SM's subbands couple to the SC's subbands with comparable strength and therefore changes of V g that "push" different SM's subbands to the Fermi level do not cause sudden jumps of W SC in contrast to the clean case. The behavior of W SC versus V g correlates with the dependence of Δ ind on V g , as shown in Fig. 15(b). Similar to the clean case (see Fig. 7) there is a one-to-one correspondence between the weight in SC and the induced gap. As follows from Fig. 15, the disorder in the SC increases the range of values of V g for which the SC-SM heterostructure is in the strong tunneling regime, which agrees qualitatively with recent experiments [68,98].
As discussed in the previous sections, a large W SC implies not only a large Δ ind but also a reduced g factor. Since the presence of disorder increases W SC (all other parameters being equal), the reduction of the g factor is also enhanced. From Eq. (9), we see that both the increase of Δ ind and the reduction of jgj will cause an increase of the minimal critical magnetic field for the topological phase transition. Considering that the disorder in the SC does not affect directly the SM's subbands, in particular their energy at k x ¼ 0, we conclude that the dominant effect of the disorder in the SC on the topological phase diagram is to cause an increase of B c . Therefore, in the presence of disorder the topological regions of Fig. 12(a) move towards larger values of B. The shift in V g of the boundaries of the topological phase appears to be negligible.
To study this more quantitatively, we determine the minimal critical field for a topological "dome" as previously shown in Fig. 12; we denote this minimal field as B Ã c .   Fig. 16 represent the values of B Ã c for different disorder realizations whereas the solid symbols correspond to the values of B Ã c averaged over 25 disorder realizations. These results clearly demonstrate that as U D increases B Ã c increases as well, mostly due to a decrease of g, and results in a reduction of the area of the topological regions in the ðB; V g Þ plane. For very large disorder strengths and negative gate voltages, the effective g factor becomes so small that B Ã c becomes larger than the critical field of Al. We conclude this section by summarizing that the main effect of disorder in SC is to increase the coupling between SM and SC, resulting in larger induced gap and critical fields to cross the topological phase transition.

IV. SUMMARY AND CONCLUSIONS
We study properties of SM-SC nanowires in the presence of strong external electric fields. Our method is based on self-consistent Schrödinger-Poisson calculations which treat the semiconductor and the superconductor on equal footing. This approach allows one to take into account several semiconductor subbands which, we believe, are present in current experimental devices. We find that the treatment of the SM and SC at the same level is necessary to describe the strong-coupling regime characteristic to the high-quality epitaxial nanowires [56]. Such hybrid nanowires are very promising for the topological quantum computing applications as they exhibit large proximity-induced gaps and very low subgap conductance [55,60,64].
One of the most important results of our work is to provide an insight regarding the necessary conditions for achieving the strong-coupling regime in proximitized nanowires. We find that one of the key ingredients is the presence of an accumulation layer at the interface between the SM and the SC. The presence of an accumulation layer implies a strong confinement of the semiconductor wave function close to the SM-SC interface. Without such confinement, the significant mismatch between the Fermi velocities of SM and SC would significantly reduce the induced gap. This conclusion has recently been supported by angle-resolved photoemission spectroscopy experiments that have shown that in epitaxial InAs/Al systems the band offset W is negative and therefore an electron accumulation layer is present.
We investigate the effect of an external electric field which can be used to modify the confining potential and, thus, modify properties of electronic states in SM-SC devices. We find that external electric field can be used to change the number of subbands in the semiconductor, tunneling rate, induced gap, magnitude of the effective g factor, coherence length ξ, etc. Our results show that the relation between V g and the quantities characterizing the electronic state of SM-SC quasi-1D nanowires is not trivial. The understanding of the interplay of V g , number of subbands, and electronic properties is one of the most important results of our work.
We obtain the topological phase diagram as a function of the gate voltage and magnetic field B. Previous works calculated the topological phase diagram in terms of phenomenological parameters such as effective chemical potential. Our work is the first to present a phase diagram in terms of V g , the experimentally relevant and tunable quantity, instead of the chemical potential. We find that in the strong-coupling regime the renormalization of the g factor due to the strong hybridization between the SM's and the SC's states can significantly modify the topological phase diagram. For typical setups, the g of the SC is smaller (in absolute value) than the SM's g factor, and so the strong hybridization reduces the g factor causing a decrease in the ðV g ; BÞ plane of the region where the system is in the topological phase.
Finally, we take into account the effect of disorder. There is a large body of papers investigating the effect of the disorder in the semiconductor [31,87,111,[119][120][121][122][123] and at the interface [124] concluding that disorder leads to the subgap density of states (i.e., states below the induced gap). However, given the observation of a very small subgap density of states in recent experiments on highquality proximitized nanowires [55,60,62,64], we believe that the semiconductor, as well as the SC-SM interface, is quite clean. The situation with Al is less clear since the presence of the native oxide covering Al may lead to significant impurity scattering. The disorder in the superconductor relaxes the constraint on momentum conservation and enhances the induced SC gap and the critical field. Once again, the optimization of the tunneling rate between SM-SC is very important [112,118].
Our work has important implications for current and future experiments aiming to realize Majorana-based topological qubits using SM-SC heterostructures as it allows one to optimize Majorana devices by tuning key parameters, Δ, g, and ξ with gates. Our results show that in the strong coupling the renormalization of g can be quite significant, increasing the minimal magnetic field necessary to drive the system into the topological phase. Thus, there is a sweet spot, and it is beneficial to operate in the intermediate-coupling regime. This is critical information to design experiments aimed at realizing MZMs.
Finally, we discuss some limitations of our model. First, the disorder is two dimensional, leading to a qualitatively similar picture as in the disorder-free case. Impurities in the SC may induce subgap states [112,[116][117][118] that can be captured in a fully three-dimensional simulation. Another limitation of our model is the lack of orbital effects due to the magnetic field. Because of the strong geometry dependence of the orbital effect [97,125], however, a careful treatment of it needs to go beyond the slab model discussed here [98,126].
Finally, we emphasize that, although in this work we focus on InAs/Al hybrid nanowires, the Schrödinger-Poisson approach proposed in this work can be used to study other heterostructures such as InSb/Al nanowires, two-dimensional SM-SC heterostructure, and the quasi-1D channels created by electrostatic confinement in such structures.
Recently, three other papers on a similar topic appeared [127][128][129]. These papers, as ours, present numerical approaches aimed at a more quantitative description of semiconductor-superconductor heterostructures. Reference [127] studies the electrostatic environment of nanowires in the weak-coupling regime. Reference [128] discusses the electrostatic environment in the presence of metallic Al and focuses on the normalstate band structure of the quasi-one-dimensional heterostructure. Reference [129] studies the renormalization of the semiconductor band structure by the proximity to the bulk superconductor in the strong-coupling regime neglecting electrostatic effects.