A Analytic Green’s function for the infinite wire using conformal mapping

Effective theory approach to the Schrödinger-Poisson problem in semiconductor Majorana devices


We propose a method for solving to the Schrödinger-Poisson problem that can be efficiently implemented in realistic 3D tight-binding models of semiconductor-based Majorana devices. The method is based on two key ideas: i) For a given geometry, the Poisson problem is only solved once (for each local orbital) and the results are stored as an interaction tensor; using this Green’s function scheme, the Poisson component of the iteration procedure is reduced to a few simple summations. ii) The 3D problem is mapped into an effective multi-orbital 1D problem with molecular orbitals calculated self-consistently as the transverse modes of an infinite wire with the same electrostatic potential as the local electrostatic potential of the finite 3D device. These two ideas considerably simplify the numerical complexity of the full 3D Schrödinger-Poisson problem for the nanowire, enabling a tractable effective theory with predictive power. To demonstrate the capabilities of our approach, we calculate the response of the system to an external magnetic field, the dependence of the effective chemical potential on the work function difference, and the dependence of the effective semiconductor-superconductor coupling on the applied gate potential. We find that, within a wide range of parameters, different low-energy bands are characterized by similar effective couplings, which results in induced gap features characterized by a single energy scale. We also find that electrostatic effects are responsible for a partial suppression of the Majorana energy splitting oscillations. Finally, we show that a position-dependent work function difference can produce a non-homogeneous effective potential that is not affected by the screening due to the superconductor and is only partially suppressed by the charge inside the wire. In turn, this potential can induce trivial low-energy states that mimic the phenomenology of Majorana zero modes. Thus, any position-dependent work function difference (even at the 1 level) along the nanowire must be avoided through carefully engineered semiconductor-superconductor interfaces.


[name=Supplementary Figure]suppfigure

I Introduction

Motivated by the theoretical model proposed by KitaevKitaev (2001) and the concrete predictionsSau et al. (2010a); Alicea (2010); Oreg et al. (2010); Lutchyn et al. (2010); Sau et al. (2010b) about the existence of zero energy Majorana modes in proximity-coupled semiconductor-superconductor (SM-SC) hybrid structures, a systematic experimental search for Majorana zero modesMajorana (1937); Read and Green (2000) (MZMs) has gained momentum in the past few years.Mourik et al. (2012); Deng et al. (2012); Das et al. (2012); Rokhinson et al. (2012); Churchill et al. (2013); Finck et al. (2013); Chang et al. (2015); Albrecht et al. (2016); Deng et al. (2016); Chen et al. (2017) Recent improvements in materials science and nanofabricationZhang et al. (2017a); Nichele et al. (2017); Zhang et al. (2017b) have led to the observation of stable zero energy subgap states that manifest the predicted quantization of the zero bias tunneling differential conductance at low temperatures.Sengupta et al. (2001); Law et al. (2009); Flensberg (2010) The signatures observed experimentally provide strong indication that MZMs localized at the ends of proximitized semiconductor nanowires may have been realized in the laboratory. However, based on the existing evidence one cannot rule out the possibility that these experimental signatures are, in fact, generated by non-topological Andreev bound states (ABSs), which are ubiquitous in the presence of non-uniform system parameters (e.g., variations of the electrochemical potential) or when the wire is coupled to a quantum dot.Kells et al. (2012); Prada et al. (2012); San-Jose et al. (2013); Ojanen (2013); Stanescu and Tewari (2014); Lee et al. (2014); Cayao et al. (2015); Klinovaja and Loss (2015); Stanescu and Tewari (2016); Liu et al. (2017); Moore et al. (2017) In particular, the possible presence of partially separated ABSs (ps-ABSs) consisting of pairs of Majorana bound states separated by a distance comparable to or larger than the characteristic Majorana length-scale (but less than the length of the wire) should raise serious concern, as one cannot distinguish between these trivial low-energy modes and genuine non-Abelian MZMs using any type of local measurement at the end of the wire.Moore et al. (2017) Considering this rather disturbing state of affairs, in conjunction with the promising proposalsAlicea et al. (2011); Sau et al. (2011); Aasen et al. (2016); Karzig et al. (2017); Litinski et al. (2017) for testing the predicted non-Abelian properties of the MZMs and building topological qubits, which will require exquisite control of the hybrid system, it becomes clear that a major theoretical task is to develop a more detailed modeling of semiconductor-superconductor Majorana devices. The minimal model used extensively so far in the Majorana nanowire literature is simply insufficient for describing the SM-SC structure at the level of essential details necessary to distinguish between MZMs and ABSs, as well as in the elucidation of other basic properties of the hybrid device.

A key component of this task is to account for the electrostatic effects that are naturally induced by the presence of the superconductor-semiconductor interface and external potential gates. Understanding these effects is critical in the context of two important aspects of the modeling of hybrid devices. On the one hand, they control three basic system parameters: the chemical potential, the Rashba spin-orbit coupling, and the induced superconducting pair potential. Typically, these parameters are treated as independent phenomenological parameters. In fact, they are all controlled by the effective electrostatic potential inside the wire generated by the work function difference at the SM-SC interface and by the applied gate potential in the presence of a low (but non-vanishing) electron density. The work function difference and the gate potential determine the number of charge carriers in the wire (hence, the value of the chemical potential relative to the bottom of the conduction band). In addition, the transverse profile of the effective potential is directly linked to the Rashba spin-orbit coupling and determines the amplitudes of the wave functions at the SM-SC interface, which, in turn, control the strength of the proximity-coupling to the superconductor. Understanding the dependence of these system properties on control parameters such as external gate potentials and applied magnetic fields is important for correctly interpreting the experimental data and optimizing the Majorana devices. On the other hand, electrostatic effects are critical ingredients of existing and proposed Majorana devices, ranging from the controllable tunnel barrier in a charge transport measurement, to the electrostatic confinement in two-dimensional SM-SC structures,Suominen et al. (2017); Hell et al. (2017a, b) and electrostatic operations in Majorana nanowire-based topological circuits,Alicea et al. (2011); Sau et al. (2011); Aasen et al. (2016); Karzig et al. (2017); Litinski et al. (2017) while being a major potential source of unwanted inhomogeneity in the active segments of these devices (i.e. those that host the non-Abelian MZMs). From this perspective, understanding in detail the electrostatic effects in semiconductor Majorana devices represents a requirement. Clearly, the minimal model in which all of these crucial parameters (e.g. chemical potential, spin-orbit coupling, proximity-induced pair potential) are assumed to be independent adjustable parameters, is highly inaccurate (and perhaps even incorrect) and non-predictive, since these parameters cannot be freely tuned in any experimental hybrid system by controlling the electrostatic environment (i.e. various gate voltages).

In general, accounting for electrostatic effects requires solving a Schrödinger-Poisson problem self-consistently. The Schrödinger-Poisson problem is important in understanding the properties of low-dimensional semiconductor structures and, indeed, over the years many self-consistent treatments have been carried out in semiconductor inversion and accumulation layers,Stern (1972); Ando (1976); Das Sarma and Vinter (1982) semiconductor heterojunctionsStern and Das Sarma (1984) and quantum wells,Bastard et al. (1983); Bauer and Ando (1986) semiconductor nanowires,Lai and Das Sarma (1986); Laux and Stern (1986) and semiconductor quantum dots.Stopa (1996) Most of these self-consistent theories are carried out within the continuum effective mass approximation (sometimes with additional approximations to simplify the numerics) where the self-consistency is limited to the electrons in the semiconductor itself, thus motivating our work. In general, these theories capture the electronic structure of the low-dimensional semiconductor systems extremely well,Ando et al. (1982); Bastard (1988) and have become a standard tool in the semiconductor industry. Our goal here is to develop a similar self-consistent tool in hybrid structures with SM-SC interfaces, whereas by contrast the standard low-dimensional semiconductor systems have typically SM-SM (or SM-insulator, as in Si MOSFETs) interfaces. The presence of superconductivity, spin-orbit coupling, and magnetic field makes our problem much richer (and more difficult technically) than the above-mentioned pure semiconductor low-dimensional systems.

The most relevant components that determine the electrostatic effects in SM-SC structures are the applied gate potentials, the work function difference at the SM-SC interface, and the screening due to the presence of the superconductor and the finite charge in the wire. The topological superconducting phase and the emerging MZMs have been been found to be relatively stable against disorder and weak interaction.Gangadharaiah et al. (2011); Stoudenmire et al. (2011); Sela et al. (2011); Lutchyn et al. (2011); Lobos et al. (2012); Crépin et al. (2014); Manolescu et al. (2014) Considering the properties of the semiconductor materials used in the fabrication of Majorana devices and the strong screening by the superconductor, it is reasonable to assume that the main effects of electron-electron interaction are faithfully captured at the mean-field level, i.e. within the Hartree approximation. Exchange-correlation effects may have some small quantitative effects, but given that the typical semiconductor materials used in Majorana nanowires (e.g. InSb and InAs) have very small electron effective masses and rather large lattice dielectric constants, we expect such exchange-correlation corrections to be rather negligible since the relevant dimensionless interaction coupling constant (the so-called value) is very small. Therefore, the task at hand is to find a self-consistent solution of a three-dimensional (3D) Schrödinger-Poisson problem associated with a given semiconductor-superconductor Majorana device. This task, however, poses a significant challenge due to the enormous number of relevant degrees of freedom that have to be taken into account. A possible path would be the brute force approach to the 3D Schrödinger-Poisson problem. This could be helpful in the engineering process of a specific device, but has two major disadvantages: it is an extremely costly numerical scheme and it provides virtually no additional understanding of the relevant physics and has limited predictive power. An additional (and rather serious) empirical problem associated with a brute-force 3D Schrödinger-Poisson approach is that the relevant experimental parameters are simply not known at the level of accuracy necessary for such a method to provide reliable results at the eV energy scale operational for the MZM problem of interest here.

In this work we propose and develop an alternative approach involving an effective theory of the 3D Schrödinger-Poisson problem that can be efficiently implemented numerically and can provide insight into the low-energy physics of the SM-SC device, particularly in terms of the dependence of key low-energy features on the SM-SC materials parameters and the applied gate voltages. Our method is based on two key ideas: i) We split the actual 3D problem into a 2D problem corresponding to an infinite (uniform) wire and an effective multi-orbital 1D problem with “molecular” orbitals calculated (self-consistently) using the infinite 2D system. ii) For a given geometry, the Poisson problem is solved once for each lattice site and the results are stored; using this Green’s function scheme, the Poisson component of the iteration procedure becomes trivial. More specifically, we first consider an infinite nanowire-superconductor system in the presence of an external gate potential that is translation invariant (along the wire) and calculate the transverse profiles of the wave functions associated with each confinement-induced band by solving self-consistently the corresponding 2D Schrödinger-Poisson problem. Next, we construct an effective multi-orbital 1D model of the 3D device by dividing the system into “slices” and associating to each “slice” molecular orbitals given by the transverse profiles of the confinement-induced bands corresponding to an infinite wire with the same electrostatic potential as the local electrostatic potential of the “slice,” which is obtained by solving a 3D Laplace equation. Of course, including all the bands would simply imply a change of basis. The point is that the subspace spanned by a relatively small number of bands calculated self-consistently by solving the (2D) infinite wire problem provides a good approximation for the low-energy Hilbert space of the 3D system. The projection reduces the numerical complexity of the problem enormously, since it eliminates a large number of (irrelevant) high-energy degrees of freedom that have to be considered when using the brute force approach to the full 3D problem. We note that both the 2D problem and the 1D effective model are solved self-consistently. The first self-consistency condition ensures that the calculated transverse profiles (hence, the “molecular” orbitals) accurately include interaction effects (at the Hartree level of Coulomb energetics), while the second condition ensures that the charge is correctly distributed along the wire (within the same approximation). This effective approach is both computationally efficient and physically substantive, as demonstrated explicitly in the current work, being characterized by numerical tractability and predictive power.

This work focuses on a method to effectively solve the Schrödinger-Poisson problem in semiconductor Majorana devices, elucidating the implicit approximations as well as additional possible simplifications and refinements of the proposed approach. In addition, we provide specific examples of how one can use this method to address important questions regarding the low-energy physics of proximity-coupled SM-SC structures. We first consider the case of an infinite semiconductor wire in the presence of an external gate potential and a work function difference at the interface between the wire and the superconductor. We calculate the response to an external magnetic field and compare the predictions based on first order perturbation theoryVuik et al. (2016) with the fully self-consistent results. We also calculate the dependence of the “effective chemical potential” (in fact, the energies of the interacting semiconductor bands) on the work function difference and show that the corresponding linear coefficient is of order unity. By contrast, the dependence on the applied gate potential is strongly suppressed due to the screening provided by the superconductor. We also investigate the dependence of the band-dependent induced pair potential on the work function difference and the applied potential and find that the low-energy bands are characterized by similar values of this parameter, in sharp contrast with predictions based on simple noninteracting models. Next, we consider a finite wire and investigate the energy splitting oscillations of the Majorana modes arising from the overlap of the MZMs at the two wire ends of the wire. We find that interaction partially suppresses these oscillations, which is an effect arising from the self-consistency in the Schrödinger-Poisson solution. We then consider a finite system with a nonuniform work function difference at the SM-SC interface. This non-uniformity in the work function could arise, for example, from physical structural fluctuations at the interface, which are invariable at the few mono-layer level even in the best epitaxial interfaces. We find that small variations of the work function difference (of the order of ) can generate variations of the effective electrostatic potential larger than the induced gap. The screening by the superconductor plays no role in suppressing the emergence of this inhomogeneous potential, while the screening by the charge inside the wire is only effective at high occupancies. This calculation provides concrete support to the possibility of long-range potential inhomogeneities in proximitized nanowires, which are predictedStanescu and Tewari (2016); Liu et al. (2017); Moore et al. (2017) to induce trivial low-energy states that mimic the (local) signatures of non-Abelian MZMs. We note that the typical absolute work function at the SM-SC interface is of the order of hundreds on meVs whereas the relevant low-energy energy scale (e.g. the induced gap in the nanowire) is only eV, making the homogeneous control of the work function along the whole SM-SC interface a rather formidable materials science, fabrication, and engineering challenge, which must be solved for future progress in the field. We mention as an aside that the work function inhomogeneity issue discovered in the current work is quite distinct from the short-range disorder problem associated with the SM-SC interface discussed earlier in the literature within the minimal model.Takei et al. (2013)

The rest of the paper is organized as follows. In Sec. II we present our approach to the Schrödinger-Poisson problem in proximitized semiconductor nanowires. We describe the Green’s function scheme, its implementation in the case of infinite nanowires, and the scheme for constructing and solving the effective 1D problem corresponding to finite systems. In Sec. III we apply our method to infinite Majorana nanowires and investigate the response to an external magnetic field, the dependence of the effective chemical potential on the work function difference at the SM-SC interface, and dependence of the proximity-induced pair potential on the relevant parameters. Sec. IV is dedicated to finite hybrid structures of experimental relevance. We discuss the suppression of the Majorana splitting oscillations due to interaction and the emergence of inhomogeneous potentials in systems with a non-uniform work function difference. We conclude in Sec. V with a summary of the results and a discussion of the relevance of this work to future studies of Majorana systems.

Ii Theoretical methods

In this section we describe our approach to the Schrödinger-Poisson problem in proximitized semiconductor nanowires. We discuss A) the Green’s function scheme, B) the infinite wire case, and C) the effective 1D problem. We restrict ourselves to the weak coupling regime, i.e. we assume that the low-energy wave functions have almost all their weight inside the semiconductor nanowire (with an exponentially-small tail penetrating inside the superconductor). Within this approximation, the superconductor can be treated as a “boundary condition” for the electrostatic potential. The strong/intermediate coupling regime, which is expected to exhibit interesting new physics at low energies,Stanescu and Das Sarma (2017a) will be discussed elsewhere.

ii.1 The Green’s function scheme

Consider a -dimensional semiconductor system described by a multi-orbital tight-binding Hamiltonian of the form


where is a non-interacting Hamiltonian, which includes hopping terms, spin-orbit coupling, and external field contributions, and accounts for the electron-electron interaction. At the mean-field level, the interaction term has the form


where and label the lattice on which the tight-binding model is defined, and are combined orbital and spin indices, and are matrix elements of the Hartree potential with the basis states of the tight-binding model. The operator creates an electron in a single particle state with orbital/spin index centered at site (i.e. the state ). The Hartree (or Coulomb) potential satisfies the Poisson equation


where is the background dielectric constant of the semiconductor and is the charge density. In turn, the charge density can be expressed in terms of the eigenstates of Hamiltonian (1) as a sum over the occupied states,


Equations (1-4) define a Schrödinger-Poisson problem that has to be solved self-consistently. The self-consistency arises from the fact that the eigenstates , which define the charge density through Eq. (4), are in turn determined by the charge density through Eqs. (2) and (3). We note that having a unique solution of the Poisson equation (3) requires specified boundary conditions. Also, in general, the non-interacting Hamiltonian contains an external electrostatic potential generated, for example, by an applied gate voltage. Finding the spatial dependence of this external potential may require solving an additional Laplace equation, which involves knowledge of various geometrical and materials details characterizing each given device. It is convenient to solve the Poisson equation (3) with homogeneous boundary conditions and incorporate all non-homogeneous contributions (e.g., a non-vanishing gate voltage) into the boundary conditions of the Laplace equation. Note that the Laplace equation has to be solved once (for a given external potential configuration), while the Poisson equation has to be solved self-consistently, together with the Schrödinger problem defined by Hamiltonian (1), within an iterative scheme, which can be computationally expensive. For example, having to solve the Poisson equation numerically at every iteration represents a serious practical obstacle when exploring the large parameter space that typically characterize the heterostructure model. In addition, numerical accuracy demands very precise solutions of the Poisson equation, making this the essential roadblock in the efficiency of the computational scheme.

To address this challenge, we reformulate the problem so that the Poisson component of each iteration becomes trivial. First, we write the eigenstates in terms of the localized basis states as


Defining , we can write the charge density in the form


where are local orbitals. Note that the second term in Eq. (6) is due to orbital overlap and can be neglected in single-band models (see below).

Next, we introduce the Green’s function defined by the equation


with homogeneous boundary conditions. Note that represents the electrostatic potential generated by an electron occupying the orbital at site . Finally, we define the following “interaction tensor”:


The element represents the interaction energy between two electrons occupying the orbitals at site and at site , respectively. Note that, in general, the Green’s function defined by Eq. (7) and the interaction tensor defined by Eq. (8) are complex quantities. Using these quantities, we can write the matrix elements of the Hartree potential in the form


Our strategy is to solve Eq. (7) for every lattice site in the system, which can be done numerically or, in some cases, even analytically (see, for example, Appendix A), perform the integration in Eq. (8), and store the interaction tensor. The Poisson component of the iterative scheme reduces to tensor contraction in Eq. (9). We note that, in practice, many elements of the interaction tensor are small and can be safely neglected. Also, the calculation of the interaction tensor using Eq. (8) requires knowledge of the basis states , which can be found using ab-intio techniques. In the applications discussed in this work we only consider single-orbital models and we assume that has spherical symmetry and is strongly localized near site (i.e. we neglect the overlap with neighboring orbitals).

The general scheme described above simplifies significantly in the case of single-orbital tight-binding models. Since the only internal degree of freedom is spin, we have , where is the spin index. Furthermore, the spatial profile of the local orbital is spin-independent, so that we have , and we neglect the overlap between neighboring orbitals, . With these simplifications the relevant Green’s function that has to be calculated (for each lattice site ) becomes


and the interaction tensor (8) reduces to an interaction matrix,


Note that is simply the effective Coulomb interaction energy between two electrons at sites and , respectively. Finally, the interaction term from Hamiltonian (1) becomes local and can be expressed in terms of the matrix elements of the Hartree potential as


The usefulness of this method becomes clear if we consider exploring a large parameter space within a given device geometry. As long as the geometry of the system remains fixed, we can change various system parameters, such as back gate potentials, magnetic fields, and spin-orbit couplings, while using the same interaction matrix, which is determined once at the beginning of the calculation. Moreover, since finite element computational methods can automatically handle unconventional and complicated device geometries, this method can be applied to devices having arbitrary shape, with any number of gates, different dielectric materials, and arbitrary spatial dimension . Thus, the method described above is of wide applicability to actual systems of experimental relevance.

ii.2 Schrödinger-Poisson scheme for infinite nanowires

While using the Green’s function method makes a 3D Schrödinger-Poisson problem significantly more manageable (this technique being clearly preferable to the pure brute force self-consistent approach), a direct 3D calculation may still be prohibitively costly due to the large number of (relevant) degrees of freedom. To overcome this challenge, we split the 3D problem into a 2D problem corresponding to an infinite uniform wire and an effective 1D problem associated with the actual finite structure. In this section we describe the self-consistent procedure for solving the 2D Schrödinger-Poisson problem using the general framework discussed above.

Figure 1: (Color online) Typical transverse profile of a Majorana SM-SC heterostructure. The SM nanowire (yellow) is partially covered by an s-wave SC (blue) and placed on an insulating substrate (light red). A back gate (black) creates a controllable electrostatic potential.

Consider an infinite quasi-1D semiconductor (SM) nanowire proximity coupled to an s-wave superconductor (SC). The axis of the wire is oriented along the -direction, while the finite cross section has a geometry similar to that shown in Fig. 1, which is the typical experimental setup for Majorana nanowires. The semiconductor nanowire (e.g., InSb or InAs) is partially covered by an s-wave superconductor (e.g., Al or NbTiN) and placed on an insulating substrate. A controllable back gate allows one to change the electrostatic potential across the wire. For clarity and to avoid cumbersome notations, we restrict ourselves to single-orbital tight-binding models and we neglect the overlap between neighboring orbitals, which allows us to use the simplified version of the Green’s function scheme described above. However, we emphasize that the approach is generic and can be directly generalized to the multi-orbital case. The non-interacting part of the Hamiltonian describing the nanowire has the form


where are position labels in the transverse (-) plane (i.e. normal to the nanowire direction taken to be the -direction throughout) and creates an electron at position with longitudinal wave vector and spin . Note that the lattice is only defined inside the SM nanowire. In Eq. (13) are matrix elements for hopping across the wire, (with being the effective mass) is the longitudinal component of the kinetic energy, represents the external potential at site arising from the back gate and the work function difference at the SC-SM interface, and is a reference energy (determined by the value of the SM band gap and the possible presence of dopants) that controls the minimum of the (noninteracting) spectrum for an isolated SM wire. In the last two terms, represents the (half) Zeeman splitting due to a magnetic field applied parallel to the wire, is the Rashba spin-orbit coupling, and (with ) are Pauli matrices associated with the spin degree of freedom. Note that the (infinite) wire has translational invariance in the -direction and, therefore, is a good quantum number. Also, we assume that the SM-SC coupling is weak, which means that the SC can be treated as i) a source of Cooper pairs for the wire (with pairing potential ) and ii) a boundary condition for the electrostatic problem. The weak-coupling assumption, used extensively in the Majorana nanowire literature, enables one to integrate out all the complications of the underlying superconductor in terms of a single pairing potential parameter characterizing the induced proximity effect.

The electrostatic potential has to be calculated by solving a Laplace equation with boundary conditions determined by the geometry of the problem and by two key parameters: the gate voltage and the work function difference at the interface, (see Fig. 1). We emphasize that, for a given SM model, the parameters , , and completely determine the carrier concentration in the nanowire and the transverse profiles of the wave functions and effective electrostatic potential (which includes the interaction effects at the mean-field level). Hence, the chemical potential of the wire (relative to, e.g., the bottom of the spectrum), the Rashba coefficient , and the induced pairing potential are not independent parameters (as implicitly assumed in the extensively used minimal model), but rather functions of , , and , the actual independent parameters of the microscopic theory.

The interaction effects are incorporated at the mean field (Hartree) level by adding to Hamiltonian (13) the term


where are the matrix elements of the Hartree potential. These matrix elements are determined by the interaction matrix (11) and by eigenstates of the full Hamiltonian , where labels the confinement-induced transverse modes. Explicitly, we can write the matrix elements of the Hartree potential in the form


Solving the Schrödinger-Poisson problem for the infinite wire implies solving Eq. (10) with homogeneous boundary conditions (once) for each lattice site corresponding to a transverse section of the wire, calculating and storing the interaction matrix given by Eq. (11), then solving self-consistently the Schrödinger problem for with the matrix elements of the Hartree potential being given by Eq. (15).

Figure 2: (Color online) Left: Confinement-induced bands for an infinite wire with in the absence of an applied Zeeman field (i.e. for ). Right: Applying a magnetic field () splits the bands into pairs of spin sub-bands (solid lines). The dashed lines represent the (fictitious) energy dispersion corresponding to that is used in the definition of the effective chemical potential . For each band is defined with respect to the minimum (at ) of the corresponding dashed line. The effective chemical potential of the second band, , is shown as an example.

A few comments regarding the practical implementation of this scheme are warranted. First, we note that the basis states of the tight-binding model are typically unspecified. Moreover, we often deal with effective tight-binding models defined on a lattice having a unit cell much larger then the atomic unit cell of the semiconductor. Hence, should not necessarily be regarded as atomic-type orbitals. In such cases, a reasonable approximation that can be easily implemented numerically is based on the assumption that the charge associated with is uniformly distributed throughout the unit cell. Second, we note that, imposing only minor additional restrictions, we can find an analytic solution of Eq. (10). The main idea is to solve the Poisson problem in a cylindrical geometry, then use a conformal mapping to obtain the results for, e.g., a hexagonal wire. The details of this calculation are provided in Appendix A.

Finally, let us discuss qualitatively the effect of the (mean-field) electron-electron interaction on the energy spectrum of the infinite SM wire. A quantitative analysis will follow in Sec. III. The transverse confinement of the nanowire gives rise to confinement-induced one-dimensional sub-bands (henceforth referred to as “bands”), as shown in Fig. 2. Here, for simplicity, we take . We define the effective chemical potential measured relative to the bottom of a given band as , where is the n energy eigenvalue of corresponding to (i.e. no Zeeman field) and . In the presence of a Zeeman field, the effective chemical potential is defined as , where is the energy (at ) of the corresponding spin-split sub-band (see Fig. 2). Defining such a quantity can be useful in the context of Majorana physics, for example when discussing the “topological condition”, , where is the topmost occupied band. Note that is positive for occupied bands and negative for empty bands. Neglecting interactions results in a chemical potential that is independent of the applied Zeeman field, However, due to interaction effects, the dependence of on control parameters such as the Zeeman field becomes nontrivial. Indeed, turning on splits each band into two spin sub-bands, as shown in Fig. 2. With increasing , the higher energy spin sub-band “loses” occupied states, while its lower energy partner gains occupied states. The net gain (or loss) is, in general, nonzero, which implies that the occupation of each band will change and, consequently, the Hartree potential (15) will change. In turn, this shifts the effective chemical potential of each band by an amount that has to be determined self-consistently. We conclude that applying a magnetic field does not simply split the bands. Instead, due to interactions, the Zeeman effect has to be supplemented by band-dependent energy shifts that can only be determined by solving the Schrödinger-Poisson problem self-consistently. Hence, the effective chemical potential varies with , leading to important consequences regarding the dependence of various low-energy features on the applied magnetic field.

ii.3 The effective 1D problem for finite systems

Consider a finite nanowire oriented along the x-direction and having a certain transverse profile. We divide the wire into layers (or slices), each containing sites, as shown schematically in Fig. 3. The corresponding 3D Hamiltonian has the form


where creates an electron with spin localized near the site of layer , is the number operator, and are intra- and inter-layer nearest neighbor hopping matrix elements, respectively, is the (half) Zeeman splitting, and is the Rashba spin-orbit coefficient. The electrostatic effects are described by the external potential and by the mean field contribution , which will be determined self consistently using the generic Green’s function method discussed in Sec. II.1 and the procedure described below. Note that a reference energy that controls the minimum of the (noninteracting) spectrum for an isolated SM wire [see Eq. (13)] can be incorporated into .

First, for each layer we define the following auxiliary Hamiltonian:


where . The auxiliary model, which describes an infinite wire, is defined on a lattice with a transverse profile that matches the lattice of layer , i.e. the local transverse profile of the original 3D system. Note that, Hamiltonian (17) represents a specific case of the infinite wire problem considered in Sec. II.2 corresponding to an external potential and no Zeeman field, i.e. . In other words, the auxiliary Hamiltonian describes an infinite system in the presence of a translation-invariant external potential that matches the local external potential of the actual 3D wire on layer .

Figure 3: (Color online) Schematic representation of the layers (slices) used for constructing the effective 1D model. A generic site of the 3D lattice is labeled , where is the layer index and indicates the transverse position within the layer.

The -independent transverse components of the single-particle eigenstates of the auxiliary Hamiltonian have the form


where is the local orbital at site and is a band index. Note that the label for the spin degree of freedom has been suppressed. By convention, the Roman letters label (transverse) positions within the wire, as well as the corresponding local orbitals. On the other hand, the Greek letters will be used to designate confinement-induced bands and the “molecular orbitals” associated with the transverse profile of the corresponding band.

Next, we perform a change of basis in the tight-binding Hamiltonian , from the local orbitals to the molecular orbitals given (for each layer) by the eigenstates of the auxiliary problem (17). For convenience, we introduce the “vector” operator with components labeled by . Here, we have , , and or , so that the total number of degrees of freedom (which gives the size of ) is . Similarly, we label the molecular orbital basis with . Using these notation, we rewrite the Hamiltonian (16) in a more compact (and generic) form as


where the nonzero matrix elements match the corresponding quantities from Eq. (16). The structure of these matrices is discussed in Appendix B. Now let be the transformation matrix that generates the desired change of basis. The element of the transformation matrix corresponding to and is given by the coefficient in Eq. (18), . Inserting the identity in Eq. (19) and defining the annihilation operator for the molecular orbital, , leads to


where for all matrices from Eq. (19). The potential is the mean field contribution determined self-consistently by solving the auxiliary problem (17) for each layer , i.e. for we have . The additional term represents the difference between the mean-field potential calculated self-consistently for the original 3D problem and . Noticing that the quantity between the square brackets in the first term of Eq. (20) is nothing but an eigenvalue of the auxiliary Hamiltonian (17) for , we can write the 3D Hamiltonian in the form


So far, we have made no approximation; the physics described by Eq. (21) is exactly the same as that described by the original Hamiltonian (16). However, the key point of this construction is that the low-energy sub-space of the original problem (which is the relevant sub-space for understanding Majorana physics) is well approximated by the low-energy subspace spanned by a relatively small number of molecular orbitals. In other words, we can project the 3D Hamiltonian onto the low-energy sub-space spanned by the molecular orbitals with . The projection generates the following effective 1D Hamiltonian


where and label the sites of the (finite) 1D lattice, and designate the molecular orbitals, and the summations marked by a symbol are restricted to the lowest energy orbitals, i.e, . The hopping matrix elements can be written in terms of the hopping matrix between layers and as


Starting with nearest-neighbor hopping in Eq. (16) results in an effective 1D model with nearest-neighbor hopping . Note, however, that the hopping matrix elements of the effective Hamiltonian are, in general, orbital- and position-dependent. The position dependence and orbital mixing can be particularly strong at the ends of the wire or inside the transition regions between a segment of the wire that is covered by a superconductor and a segment that is not covered (e.g., a tunnel barrier region). This behavior is generated by the transverse profiles (i.e. molecular orbitals) being position-dependent inside the transition region. Similar considerations also apply to the spin-orbit coupling term. However, for a quantitative description of position-dependent spin-orbit coupling one should start with a more detailed model of the 3D wire, e.g., using an 8-band Kane-type Hamiltonian, rather than the simple phenomenological term discussed here. This is certainly doable, but unnecessary at this stage, as we focus on the basic ideas of the effective theory. Nonetheless, it is important to emphasize that, based on the present analysis, we can conclude that accurate modeling of inhomogeneous regions such as, for example, the tunnel barrier region at the end of a proximitized wire, using effective 1D Hamiltonians should necessarily involve position-dependent hopping/spin-orbit coupling and orbital mixing terms, in addition to the potential barriers that are typically considered in the literature. This physics of the position dependence is not accounted for in the usual minimal model of Majorana nanowires.

Calculating the matrix elements of the mean-field potential is a straightforward extension of the Green’s function method discussed in Sec. II.1. Let be an eigenstate of the effective Hamiltonian (22). Expanding it in terms of molecular orbitals, , then in terms of local orbitals, , we have


The interaction matrix of the original 3D problem is determined by solving equations (10) and (11) for the corresponding system. This encodes the interaction energy between two electrons occupying the local orbitals and , respectively. It is convenient to define the molecular orbital interaction tensor given by


where . Note that has the same structure as the interaction tensor (8), with and . In particular, the element represents the interaction energy between two electrons occupying the molecular orbitals on site and on site , respectively. Finally, using the results of Sec. II.1, one finds that the matrix elements of the mean-field potential are given by


where and we have subtracted the matrix elements of the mean-field potential associated with the auxiliary problem (17).

We conclude this section with a summary of our approach to the Schrödinger-Poisson problem in semiconductor Majorana devices. Assume that a finite nanowire described by a 3D tight-binding model, e.g., the Hamiltonian given by Eq. (16), is weakly coupled to a superconductor. The first step is to calculate the external electrostatic potential by solving a Laplace equation with appropriate boundary conditions. The result will depend on the geometry of the system, as well as the applied gate potential (or, more generally, in a system with multiple gates) and the work function difference at the SM-SC interface, . Second, we divide the nanowire into layers and solve the auxiliary (infinite wire) problem (17) for each layer, following the self-consistent procedure described in Sec. II.2. The third and final step involves solving the effective 1D problem (22) self-consistently using the matrix elements (26) of the mean-field potential. The essence of the approximation involved in this procedure is the ansatz that the transverse profiles of the low energy states at a given location along the wire are similar to the profiles of the low-energy confinement-induced bands of an infinite wire under the same electrostatic conditions. The theory includes mode mixing due to off diagonal terms in the effective Hamiltonian, which allows for corrections to these profiles. If one includes enough molecular orbitals into the basis of the effective model, the low-energy physics of the system is accurately described. One can systematically check if enough orbitals have been included by increasing and monitoring the convergence of the results. Finally, we emphasize that both the auxiliary (infinite wire) problem and the effective 1D problem are solved self-consistently. Using the Green’s function approach reduces the Poisson components of these problems to the summations in Eqns. (15) and (26), respectively.

Iii Electrostatic effects in infinite wires

In this section we illustrate the implementation of the general scheme described above focusing on the infinite wire case. We address three basic questions: i) how are the spectral features (in particular the effective chemical potential) modified by the presence of an external Zeeman field? ii) what is the dependence on the work function difference ? and iii) how does the effective SM-SC coupling depend on the back gate voltage ? Throughout this section we consider an infinitely long wire of radius nm (see Fig. 1) described by a Hamiltonian given by Eqns. (13) and (14). The parameters of the model correspond to an InSb nanowire and we have , where is the bare electron mass, the nearest-neighbor hopping matrix element eV, and the relative permittivity . The total number of lattice sites corresponding to the hexagonal cross section of the wire is . For simplicity, we ignore the spin-orbit coupling (i.e., we set ) and use the analytical solution of the Green’s function described in appendix A. The self-consistent Schrödinger-Poisson scheme that we use is discussed in Sec. II.2.

Figure 4: (Color online) Normalized potential profiles corresponding to the Green’s function generated by an infinite line charge placed inside the nanowire at a position given by the lattice site . Placing the charge in the vicinity of the superconductor [panel (b)] results in a strongly screened potential.

Before addressing the main questions, we make two general remarks. First, we note that the potential created by the charge inside the semiconductor is strongly screened by the superconductor and the back gate. To illustrate this point and to show the structure of the Green’s function, we calculate the potential profile created by an infinite line charge placed inside the nanowire at a position corresponding to the lattice site , i.e., we calculate the Green’s function . The results are shown in Fig.(4). Note that a charge placed near the middle of the wire, i.e. far from the SM-SC interface and the back gate [panel (a)], generates a potential characterized by a spatial extent much larger than that of a potential created by a charge in the vicinity of the SM-SC interface [panel (b)]. This implies that the effect of Coulomb interactions (at the mean-field level) is significantly reduced due to screening by the SC. The back gate has a similar effect. Consequently, the spatial profile of the (occupied) transverse modes is expected to determine the strength of interaction effects: the effects will be strong if the charge is located away from the SM-SC and back-gate interfaces and weak if (most of) the charge is localized in the vicinity of an interface.

Second, we would like to estimate the importance of self-consistency in solving the Schrödinger-Poisson problem. Are fully self-consistent calculations really necessary? This is obviously important from a practical viewpoint since the self-consistent procedure is computationally costly even within our effective theory approach (and hopelessly complicated in a brute-force direct 3D approach). To address this question, we compare fully self-consistent calculations with results obtained by treating electronic interactions within first order perturbation theory. We note that the first order perturbation theory relies on the assumption that the wave functions associated with different transverse modes are not affected by interactions (i.e., that they are solely determined by the external fields). Therefore, any discrepancy between the two methods is a result of the electronic interactions changing the wave function profiles. Details concerning the perturbative calculations are given in Appendix C. A comparison between self-consistent calculations and perturbative results for two different sets of parameters is shown in Fig. 5. We plot the energy difference between the eigenstates calculated using the two methods (for the lowest four bands) as function of the applied gate voltage. To understand the behavior illustrated in Fig. 5, we note that positive values of , as well as negative gate voltages , result in the electrons being pushed toward the SM-SC interface, where the screening by the superconductor reduces interaction effects. We emphasize that, even in this situation, the energies of the eigenstates are significantly renormalized by interactions, but the profiles of the wave functions are barely affected, as demonstrated by the low values of in Fig. 5 corresponding to this regime. Applying a positive gate potential moves the charge distribution toward the center of the wire, where the interaction effects are stronger. In addition, choosing a negative reference energy [see panel (b)] corresponds to the isolated nanowire being electron-doped, i.e. having more charge carriers. Increasing the charge density enhances the strength of interaction effects, including the interaction-induced change of the wave function profiles. A second factor that contributes to the enhancement of in panel (b) is a lower value of [as compared to that used in panel (a)], which diminishes the attraction of electrons toward the SM-SC interface and reduces screening. As a final comment, we note that the energy differences in Fig. 5 can be large on the scale relevant for Majorana physics. Thus, a perturbation theoretic treatment of Coluomb interaction may be quantitatively completely unreliable since the sub-band energy scale is large compared with the delicate energy scale associated with the near-zero-energy Majorana physics. Also, if one is interested in Majorana devices that contain segments of the wire that are not covered by a superconductor (e.g., a tunnel barrier region), one should expect strong interaction effects, which requires a fully self-consistent treatment.

Figure 5: (Color online) Energy difference between eigenvalues calculated i) fully self-consistently and ii) using perturbation theory as function of the back gate voltage, . Only the lowest four bands are shown. The parameters used in the calculations are: (a) mV, meV, and (b) = 50 mV, meV. Note that the energy scales in the two panels differ by an order of magnitude. The color code for the bands is: black (), green (), gray (), orange (). We note that is a measure of how strongly the wave function profiles are affected by interactions.

iii.1 Electrostatic response to an applied magnetic field

We investigate the response of the system to an applied Zeeman field focusing on the field dependence of the effective chemical potential. In Sec. II.2 we have defined the chemical potential measured relative to the bottom of a given band as , where is the energy (at ) of the corresponding spin-split sub-band (see Fig. 2). The dependence of on the Zeeman field has been studied in Ref. Vuik et al., 2016 based on a perturbative scheme. Here, we systematically compare the perturbation theory results with the fully self-consistent calculation. This has a double purpose: on the one hand it serves as a test ground for the numerical implementation of our self-consistent scheme and, on the other hand, it provides a systematic evaluation of the accuracy of the perturbation theory approach.

Figure 6: (Color online) Dependence of the effective chemical potential on to applied magnetic field for a system with single-band occupancy. The solid blue line corresponds to the analytic solution given by Eq. (51), while the orange dots are the numerical results of the fully from self-consistent calculation. The parameters that control the electrostatic properties of the system are: meV, meV, and meV.

We start with a comparison between the chemical potential calculated fully self-consistently for a system with single band occupancy and the low magnetic field analytical solution obtained in Appendix C. We note that Eq. (51) is valid in the low field (high chemical potential) regime, . The results are shown in Fig. 6. Note that the two methods are in excellent agreement, suggesting that the transverse profile of the lowest-energy band is practically independent of the applied magnetic field. In the light of the general comments made at the beginning of this section, these results are not surprising. Indeed, the relatively large (positive) and the negative gate voltage strongly push the charge toward the SM-SC interface. Increasing the Zeeman field changes the occupation of the lowest band, but the effect is too weak to modify the transverse profile. Note, however, that the energy of the band (i.e. the effective potential ) changes significantly with the applied Zeeman field.

Next, we consider several cases characterized by different values of the parameters that control the electrostatic properties of the system, , , and . In Fig. 7 we fix the intrinsic system parameters and and tune the device from a regime characterized by single-band occupancy (top panel) to a regime characterized by three occupied bands (bottom panel) by changing the applied gate voltage. The self-consistent and the perturbative results are shown as points and solid lines, respectively. The wave function profile corresponding to the highest energy occupied bands are shown in the lower left insets. Note that the effective chemical potential initially increases with the Zeeman field, until the highest energy spin sub-band is completely depleted. At higher fields, the effective chemical potential decreases to reduce the amount of the charge that is added to the system as the low-energy spin sub-band “sinks” with increasing .

Figure 7: (Color online) Dependence of the effective chemical potential on the applied Zeeman field field. The difference is defined as . Top: System with single band occupancy corresponding to the parameters meV, meV, and meV. Bottom: Same system (i.e. meV, meV), but with three occupied bands, which corresponds to applying a positive gate potential meV. The dotted and solid lines are obtained using the self-consistent approach and the perturbation method, respectively. The gray, green, and black lines represent the highest energy, middle, and lowest energy bands, respectively. The corresponding spectra are shown in the upper right insets, while the wave function profiles of the highest occupied bands are shown in the lower left insets.

The trends revealed by Figs. 6 and 7 can be naturally interpreted as corresponding to the intermediate regime between the constant chemical potential and constant density limits. Indeed, in the absence of electronic interactions the effective chemical potential is independent of the Zeeman field. In the opposite limit, which corresponds to strong interactions, the effective chemical potential for a system with a single occupied band will decrease in a (nearly) one-to-one correspondence with the (half) Zeeman splitting to maintain a constant charge density (so as to minimize the Coulomb energy cost). The situation is slightly more complicated in the case of multiple occupied bands. Nonetheless, the results in Fig. 7 show clearly that the rates of change of the effective chemical potentials with respect to are significantly lower than the expected behavior in the constant density limit.

An important feature in Fig. 7 is the good agreement between the perturbative results and the the self-consistent solutions. The agreement is slightly better in the case of a single occupied band (top panel) primarily due to the wave function profile, which is very localized near the SM-SC interface, resulting in a nearly complete screening of the electronic interactions. By contrast, in the lower panel (i.e. for a system with three occupied bands) the top band has a significant portion of its wavefunction near the center of the wire, where screening is incomplete. As a result, the wave function profile acquires a dependence on the applied Zeeman field and the self-consistency starts to matter. To clearly see why a discrepancy between the two methods implies a change in the wave function profile, we recall that within the perturbative method (see Appendix C) the response to an applied magnetic field can be obtained by solving the equation


where , and the matrix elements of the reciprocal capacitance, are calculated using the fully self-consistent wave functions (for arbitrary and ) at . Hence, good agreement between the two methods implies that is, practically, -independent, while discrepancies reveal the change of the wave function profile with the Zeeman field.

Figure 8: (Color online) Dependence of the effective chemical potential on the applied Zeeman field field, , for a doped wire with meV, meV, and meV. The dotted and solid lines are obtained using the self-consistent approach and the perturbation method, respectively. The green, black, and gray lines represent the 3, 4, and 5 energy bands, respectively.The wave function profile of the highest energy band (lower left inset) shows that most of the charge is localized away from the interfaces with the SC and the back gate. Consequently, the agreement between the fully self-consistent calculation and the perturbative result is significantly weaker than in Fig. 7.

To further test these findings, we consider a doped wire (i.e. ) with five occupied bands and a wave function profile heavily peaked in the middle of the wire. the results are shown in Fig. 8. The self-consistency is clearly more important in this case, although for low Zeeman fields the perturbation theory still provides a reasonably good approximation. In addition, we note that having two nearly degenerate top bands results in a second increase of with (in the low-field regime) associated with the depletion of a spin-split sub-band. We conclude that using perturbation theory with a reciprocal capacitance matrix calculated (fully self-consistently) at reference field (e.g., ) provides a very good approximation over a wide regime of parameters. This result can be understood in the light of our discussion of the results shown in Fig. 5. Indeed, the typical values of the Zeeman splitting are small on the energy scale corresponding to the variation of the gate voltage in Fig. 5. Hence, the wave function profiles are largely determined by the electrostatic parameters , , and and have a very weak dependence on . However, we emphasize again that the perturbative approach itself starts with a self-consistent calculation of the wave function profiles at a reference field, e.g., , and then treats the field dependence perturbatively.

iii.2 Dependence on the work function difference

A key parameter that controls the electrostatic properties of the system is the work function difference . Unfortunately, this parameter is not uniquely determined by the materials of the heterostructure (i.e. the SM and the SC), as it depends on certain details of the SM-SC interface that, in turn, are determined by the fabrication procedure, e.g., the exact procedure used for treating the SM wire surface before depositing the superconductor.Gül et al. (2017) In fact, it is rather difficult to obtain the interface work function difference experimentally, particularly at the level of accuracy (better than 1) relevant for Majorana physics in nanowires. In particular, there could very easily be sample-to-sample work function differences for the same type of SM-SC hybrid structures depending on the fabrication details. In fact, even within a single sample, there could be local position dependent variations in along the nanowire length. Here, we treat as an unknown phenomenological parameter and determine the dependence of the low-energy spectrum on this parameter by solving the Schrödinger-Poisson problem self-consistently.

Figure 9: (Color online) Dependence of the band energies (for arbitrary ) on the work function difference for a doped system with meV and two different gate potentials: (a) meV and (b) mV. Note that the slope is of the order one (more specifically, approximately ) over a wide range of parameters. When , most of the charge is located away from the SM-SC interface and the bands depend weakly on [panel (b)].

The results for a doped nanowire with meV and two different values of the gate potential are shown in Fig. 9. First, we note that in the presence of a negative gate voltage [panel (a)] the slope is approximately (i.e. of order one) for all low-energy bands and for a wide range of values. This behavior can be understood in terms of the charge being pushed toward the interface, i.e. being localized in a region where the effective potential is of the order of . Changing the sign of the applied voltage [panel (b)] results in wave functions that are more spread over the cross section of the wire. However, quite remarkably, for most of the low-energy modes still exhibit a strong dependence on . This dependence becomes weaker when . Nonetheless, for positive gate potentials and arbitrary values of there are many modes that are predominantly localized near the SM-SC interface and show a strong dependence on the work function difference. In addition, there are some modes that are localized away from the interface, which exhibit a significantly weaker dependence on . These modes are also expected to have weaker proximity-induced superconductivity, as discussed below. A major consequence of the strong dependence on illustrated in Fig. 9 is that weak inhomogeneities in the work function difference (e.g., due to the surface treatment of the SM wire) could result in significant inhomogeneities of the effective potential along the wire. For example, considering the system from Fig. 9, a 2 variation of may result in a variation of the effective potential of the order of meV, which is large (typically, by a factor of ) when compared with the induced gap. We emphasize that the screening by the superconductor plays no role in reducing these potential variations. On the other hand, screening by the charge inside the wire may suppress the inhomogeneity. We will address this problem in Sec. IV in the context of finite wires.

Figure 10: (Color online) Top: Energy eigenvalues as a function of the applied gate potential for a system with meV and meV. Bottom: Dependence of the normalized effective SM-SC coupling matrix, , on the gate voltage, , for the four lowest energy bands (). The solid and dashed lines represent diagonal matrix elements and off-diagonal elements of the form , respectively. Note that all diagonal elements have similar magnitudes, while the off-diagonal elements are negligible, until becomes comparable to .

iii.3 Effective semiconductor-superconductor coupling

The last property that we investigate in the context of infinite wires is the dependence of the effective SM-SC coupling on the applied gate potential. In general, we can define the effective SM-SC coupling asStanescu et al. (2010); Stanescu and Tewari (2013)


where are matrix elements for hopping across the SM-SC interface and is the surface Green’s function of the parent superconductor. Working within a local approximation,Stanescu et al. (2010); Stanescu and Tewari (2013) we have , with being nonzero if labels a site at the SM-SC interface and zero otherwise. As evident from Eq. (28), the effective coupling is determined by the hopping across the SM-SC interface and by the surface density of states of the parent SC. Note that the position-dependent quantity is only defined at the interface and does not contain all the information necessary for evaluating the strength of the superconducting proximity effect. Indeed, quantities such as the induced gap or phenomena such as the proximity-induced low-energy renormalizationStanescu and Das Sarma (2017a) are controlled by the band-dependent effective coupling


where is an eigenstate of the system associated with the n confinement-induced band. In general, one would expect the coupling matrix to be non-diagonal and the diagonal terms to be strongly band-dependent. Note that the key ingredient in Eq. (29) is the amplitude of the wave function at the SM-SC interface. In turn, this amplitude is determined by the electrostatic properties of the system, in particular by the parameters , , and . Consequently, the strength of the superconducting proximity effect is expected to be strongly affected by these parameters, in particular by the applied gate voltage.

Figure 11: (Color online) Top: Energy eigenvalues as a function of the applied gate potential for a system with meV and meV. Bottom: Dependence of the normalized effective SM-SC coupling matrix, , on the gate voltage, , for the energy bands closest to the Fermi energy (). The solid and dashed lines represent diagonal matrix elements and off-diagonal elements of the form , respectively. Note that all diagonal elements have similar magnitudes, while the off-diagonal elements are negligible, until becomes comparable to .

To evaluate the dependence of the effective SM-SC coupling on the applied gate potential we consider two systems characterized by the same position-dependent coupling, (i.e. independent of ) if is at the SM-SC interface and otherwise, same work function difference, meV, and different reference energies, meV and meV, respectively. The dependence of the corresponding effective coupling on the applied gate potential is shown in Fig. 10 and Fig. 11, respectively. Note that the wire with meV has occupied bands (top panel of Fig. 10), while the doped wire with meV is characterized by a significantly higher occupancy ( occupied bands). Only the values of corresponding to the bands closest to the Fermi energy are shown.

The results in Figs. 10 and 11 reveal three important features. First, we note that the off-diagonal components of are significantly smaller than the diagonal components, except in the regime characterized by large positive values of . Second, there is a clear trend: the strength of the effective SM-SC coupling decreases with increasing , i.e. as the electrons are attracted toward the back gate and away from the SM-SC interface. The trend is less clear in the system characterized by high occupancy (see Fig. 11) for . This is due to the presence of different types of modes, some localized predominantly near the SM-SC interface and some away from the interface, as discussed in the context of Fig. 9. Finally, we note that the magnitude of is about the same for several low-energy bands within a significant range of parameters. This result is somehow unexpected, considering predictions based on simple noninteracting models, and has direct experimental implications. Specifically, in a system with multi-band occupancy (a class which probably includes all experimental nanowires), the existence of (significantly) different values of the band-dependent coupling should lead to the observation of different proximity-induced gaps. By contrast, similar band-dependent couplings will lead to the observation of a single proximity-induced gap, unless a high-resolution measurement of the induced-gap features is possible. The results in Figs. 10 and 11 are consistent with the second (i.e. single-gap) scenario. Of course, a more detailed study of the specific experimental setup is necessary in order to gain complete understanding of any given device. In particular, the possibility of distinct band-dependent proximity gaps in SM-SC hybrid systems cannot be ruled out apriori.

Iv Electrostatic effects in finite wires

In this section we illustrate the implementation of our general scheme for solving the Schrödinger-Poisson problem for finite systems by constructing the effective 1D model described in Sec. II.3. As proof-of-concept examples, we consider two problems that play a major role in understanding the significance of recent experimental observations on SM-SC Majorana structures: i) the Majorana energy splitting oscillations, and ii) the emergence of trivial low-energy states (i.e. generic low-energy non-topological Andreev bound states) in inhomogeneous Majorana wires. Throughout this section we consider a finite wire of radius nm (see Fig. 1), length m, and unit cell length in the direction parallel to the wire nm, which corresponds to dividing the wire into layers. Each layer contains sites. The parameters are again taken to correspond to an InSb nanowire with and relative permittivity , while the Rashba spin-orbit coupling coefficient is set to meVÅ. We note that an analytical solution for the Green’s functions is not possible due to the broken translation symmetry. For this reason, the Green’s functions are calculated numerically using the finite element analysis software package FEniCS.Alnaes et al. (2015)

The results presented in this section are obtained using the method described in Sec.II.3 with a few additional simplifications. First, we neglect all the terms in the effective Hamiltonian (22) that are off-diagonal in molecular orbitals. Explicitly, we set for . Although naively one may expect different bands to be decoupled, in general they are not, due to mixing terms introduced by the mean field fluctuations , by the proximity-induced band coupling , and by intrinsic inhomogeneities (e.g., the presence of a barrier region at the end of the wire), which lead to spatial variations of the transverse profiles and generate off-diagonal hopping, and spin-orbit coupling, . The last two sources of band mixing can be neglected assuming weak coupling to the parent superconductor and weak inhomogeneity. Understanding the role of the off-diagonal mean field fluctuations deserves a separate study. Neglecting band mixing leads to a simplification of Eq. (26), which becomes


where we use the notation .

Figure 12: (Color online) Decay of the molecular orbital interaction tensor as function of distance, , for a system with meV, , and meV. The solid orange and blue lines correspond to the elements and , respectively. The dashed gray line represents the element . Note that the interaction tensor decays nearly exponentially, with similar characteristic length scales for all elements. The elements that are not shown in the figure have a similar behavior.

Another simplifying assumption is that the interaction matrix between two layers dependents only on the distance between them, . In other words, the 3D Green’s function for a given site inside layer is assumed to be the same as the Green’s function of corresponding site in layer , up to an overall shift by . This approximation neglects the edge effects, but reduces the number of Green’s functions that need to be calculated by a factor of . However, we expect the edge effects to be small because of the strong screening provided by the SC and back gate. Indeed, the tensor elements decay rapidly as a function of , as shown in Fig. (12). Notice the nearly exponential decay, which implies that the interaction tensor elements become negligible over distances corresponding to a few layers (i.e. lattice sites of the effective 1D model). This demonstrates that approximation is accurate everywhere, except the very edge of the wire. Also note that the diagonal and off-diagonal elements in band space are of the same order. Consequently, charge fluctuations in one band will have a large effect on the other bands through the mean field interaction term. We note that the simplifying assumption will manifestly break down if we are interested in the electrostatic properties of a tunnel barrier region (or any other type of strong inhomogeneity). In this case the full interaction tensor has to be calculated for the barrier region; the simplifying assumption is still valid inside the (homogeneous) proximitized segment of the wire.

The final simplification involves the construction of the auxiliary Hamiltonian (17) for systems with inhomogeneous electrostatic boundary conditions along the wire. We assume that the inhomogeneity is weak (in practice we consider a variation of ) and we construct the auxiliary Hamiltonian using the local boundary conditions, instead of the local electrostatic potential. More specifically, we construct as the potential of an infinite wire problem with boundary conditions given by the local boundary conditions, i.e. and , of the full 3D device. Note that, ideally, one has to solve the Laplace equation for the whole 3D devices and obtain the electrostatic potential , then construct the auxiliary Hamiltonians for each layer using . However, if the variations in the boundary conditions are small, we expect the two constructions to produce similar results, the only significant difference being the presence of spurious discontinuities in the approximate construction in regions where the boundary conditions change abruptly. Taking into account all the simplifications discussed above, the effective 1D Hamiltonian (22) becomes


where , while the off-diagonal elements are