Interaction-induced zero-energy pinning and quantum dot formation in Majorana nanowires

Majorana modes emerge in non-trivial topological phases at the edges of specific materials such as proximitized semiconducting nanowires under an external magnetic field. Ideally, they are non-local states that are charge-neutral superpositions of electrons and holes. However, in nanowires of realistic length their wave functions overlap and acquire a finite charge that makes them susceptible to interactions, specifically with the image charges that arise in the electrostatic environment. Considering a realistic three-dimensional model of the dielectric surroundings, here we show that, under certain circumstances, these interactions lead to a suppression of the Majorana oscillations predicted by simpler theoretical models, and to the formation of low-energy quantum-dot states that interact with the Majorana modes. Both features are observed in recent experiments on the detection of Majoranas and could thus help to properly characterize them.


I. INTRODUCTION
Semiconducting nanowires with strong spin-orbit interaction, like InAs or InSb, are becoming ideal systems for the artificial generation of topological superconductivity [1][2][3].In addition to its fundamental interest, such nanowires, that may host Majorana bound states (MBSs) at their ends or interfaces [4,5], constitute promising platforms for Majorana-based quantum computing devices [6][7][8][9].Progress in fabrication techniques have allowed to induce a hard superconducting gap in InAs [10] or InSb [11] nanowires with epitaxially deposited Al layer.Moreover, last generation devices exhibit a very low degree of disorder which allows them to almost reach the ballistic limit [12][13][14].
In spite of these advances, the experimental signatures of MBSs in the nanowire devices deviate significantly in several aspects from the theoretical predictions of minimal models.This is the case, for instance, regarding the behavior of the subgap conductance through the proximitized nanowire, which has been addressed in several experiments [10,[12][13][14][15][16][17][18][19].In a long wire (whose length is much greater than the induced coherence length) the presence of MBSs manifests in the appearance of a zero bias conductance peak whose width is controlled by the normal state conductance [20].However, for the typical wire lengths explored in actual experiments (which are of the order of a few µm) it is expected that the overlap between MBSs located at both ends of the wire gives rise to conventional Andreev bound states which deviate from zero energy, leading to an oscillatory pattern in the low bias conductance as a function of Zeeman field, chemical potential or wire length [21][22][23].Conspicuously, in most A nanowire of rectangular cross-section (green) lying on an insulating substrate (grey) and in contact to a thin metallic layer in one of its faces (light blue), corresponding to the parent superconductor, and two normalmetal leads at its ends (orange) separated by tunnel barriers (brown).Typical values for the dielectric constants for each region are indicated.(b) Low energy spectrum versus chemical potential µ for a wire of thickness W = 100nm and length L = 1µm.Other parameters are the spin-orbit coupling α = 20nm•meV, the induced pairing energy ∆ = 0.3meV and the Zeeman energy VZ = 2meV.Electrostatic environmentinduced zero energy pinned regions between Majorana oscillations are indicated in red.Quantum dot levels (in blue), originated at the wire's edges due to the interaction with the bulk contacts, anticross with Majorana levels and remove their zero-energy pinning.
of the available experimental data the emergence of a robust zero bias conductance peak is observed above some critical Zeeman field without the expected oscillatory pattern [12,19,24,25].Several mechanisms have been proposed to account for the reduction or lack of oscillations, such as smooth confinement [21,[26][27][28], strong spin-orbit coupling [29], position dependent pairing [30], orbital magnetic effects [31], Coulomb repulsion among the carriers in the nanowire [22], or the presence of the normal drain lead connected to the hybrid wire [32].
Another source of Majorana oscillation suppression was put forward by some of us in a recent work [33].The key realization is that MBSs in a finite-length wire posses a finite charge, typically distributed uniformly along the wire [34], which can be susceptible to electrostatic interactions with the surrounding medium.We considered the case of a grounded parent superconductor, thus avoiding the effect of a charging energy associated to the Cooper pairs, and showed that, in such case, a residual effect of interactions may arise from the image charges induced in the electrostatic environment of the nanowire.Using a simple model for the induced potential we concluded that, in typical experimental setups, interactions would lead to a pinning of the MBSs to zero energy around parity crossings and, thus, to more robust zero bias conductance peaks than predicted by the non-interacting models.
The aim of the present work is to test the validity of the predictions of Ref. [33] for the case of more realistic calculations of the induced electrostatic potential, taking into account the actual three dimensional (3D) geometry as well as the effect of nearby metallic leads.We consider the geometry depicted in Fig. 1(a), where a nanowire of rectangular cross section lies on an insulating substrate (typically SiO 2 ) and is contacted to a thin superconducting (SC) layer on one of its faces and to two bulk normal leads at both ends, separated by thin insulating barriers.In Fig. 1(a) we indicate the characteristic dielectric constants of each region, which are relevant for the calculation of the induced potential through the Poisson's equation (discussed bellow).Our aim is to solve this equation together with the Bogoliubov de Gennes equation for determining self-consistently the charge density ρ(x) along the nanowire.For this purpose we derive a generalized method of image charges which allows us to calculate the induced potential in rather general conditions, taking into account a 3D electrostatic environment as the one shown in Fig. 1(a).
We find two main effects coming from this interaction and exemplified in Fig. 1(b).One is, as stated before, the suppression of Majorana oscillations around parity crossings (zero energy crossings where the total fermion parity of the wire changes), both as a function of the Zeeman energy V Z and the wire's chemical potential µ.This effect is produced because, at each parity crossing, a finite Majorana charge Q M enters the wire from the reservoir in an abrupt fashion.If the electrostatic screening is smaller inside the wire than in the reservoirs, a repul-sive interaction is produced between the incoming charge and its images, preventing its entrance.This translates into finite regions in parameter space [in red in Fig. 1(b)] where Majorana modes are pinned to zero energy within a finite range of V Z or µ proportional to the Majorana charge Q M and the strength of the interaction.This was already shown in Ref. [33] but for a simplified dielectric profile where the presence of the superconducting shell had been ignored.We here include it and find that the size of the pinned regions decreases but the pinning effect is still present under certain conditions that we discuss in detail below.Moreover, we explain the incompressible behavior of the electron liquid within these pinned regions in terms of the Majorana wave functions and their charge.
Another important effect of the electrostatic environment unexplored before is the creation of deep potential wells at the ends of the wire close to the bulk metallic electrodes.These wells, obtained explicitly here through the self-consistent calculation, are similar to the confinement potentials typical of quantum dots.Localized quantum dot-like energy levels in these regions disperse with magnetic field (or chemical potential) and appear below the induced gap in the wire spectrum [in blue in Fig. 1(b)].In the topological regime, dot-like levels interact with Majorana states, anticrossing them when they approach zero energy.Similar phenomenology can be observed in some experiments [14,19], and has been likely found in other occasions but discarded by experimentalists looking for the simpler picture.Interestingly, it has been shown that the shape of these anticrossings can be used to quantify the degree of Majorana wave function non-locality [35,36], a prediction that has been experimentally demonstrated recently [25].Here, we show that if the dot-Majorana levels interaction occurs in a pinning region, Majorana levels are forced to depart from zero energy, revealing the existence of a finite wave function overlap between them in spite of their zero energy.We analyze this behavior again in terms of the Majorana (and dot) wave functions and their charge.
The paper is organized as follows: in the following section (Sec.II) we provide insight into the theoretical model used to treat interactions.In the next section (Sec.III A) we analyze the case in which the influence of the bulk normal leads can be neglected, recovering the pinning effect found in Ref. [33] for a repulsive electrostatic environment.However, we focus here on the electrostatic environment effects on the Majorana wave function, rather than on its spectral properties.In the next section (Sec.III B) we study the effect of including the bulk normal leads of Fig. 1(a), finding that they give rise to the formation of quantum dot-like bound states.We further analyze the interplay of such states with the MBSs.Finally, we present the conclusions of our work at the end (Sec.IV).The robustness of the pinning effect is analyzed in detail in Supporting Information (SI) 4.

II. MODEL AND THEORETICAL APPROACH
We model the electronic states along the proximitized Rashba nanowire of length L using the following single channel Hamiltonian [4,5] where are electron annihilation operators, and σ and τ are the Pauli matrices in spin and Nambu space, respectively.The model is defined by setting the parameters m * , µ, α, V Z and ∆, corresponding to the effective mass, the chemical potential, the spin-orbit coupling, the Zeeman energy caused by an external magnetic field, and the induced SC pairing potential.
In Eq. ( 1) we also include the electrostatic potential φ(x) felt by charges in the nanowire, which can be decomposed as φ(x) = φ int (x) + φ b (x), where φ int is the potential that arises from the free charges inside the nanowire, while φ b corresponds to the potential created by bound charges that emerge in the electrostatic environment.We compute the electrostatic potential using the Poisson equation where ( r) is the non-homogeneous dielectrical permittivity of the entire system and ρ( r) is the quantum and thermal average of the charge density of the nanowire obtained with Eq. ( 1).The intrinsic part φ int (x) of the potential satisfies an analogous equation with a uniform equal to that of the nanowire.The geometry depicted in Fig. 1(a) is taken into account through a piecewise ( r) function where each material is characterized by a different dielectric constant, so that ( r) changes abruptly at the interfaces.Then, assuming that the charge density in the nanowire is located along its symmetry axis (x-axis), we obtain the electrostatic potential φ b (x) using the method of image charges, as explained in detail in SI 1.More precisely, φ b is given by φ b (x) = dxV b (x, x ) ρ(x ) , where V b (x, x ) is a kernel determined in order to satisfy the proper boundary conditions.We find analytical expressions for V b (x, x ).They are simple but rather lengthy and are given in the SI for two different cases: neglecting the effect of the bulk normal leads at the wire ends and including it.The results for these two cases are analyzed in the following sections.The obtained potential φ(x) on the nanowire axis should be plugged back into Eq.( 1).The combined Poisson-Schrödinger problem must then be iterated until it achieves self-consistency.As shown in Ref. [33], the φ int (x) part of the electrostatic solution (i.e. the intrinsic electron-electron interaction part of the problem), treated at the Hartree-Fock level, has a negligible effect on the low energy spectrum in the topological regime.We may therefore concentrate only on the self-consistency with φ b (x).In SI 2 we explain in detail the self-consistent numerical method used to compute the electrostatic potential profile as well as the eigenvalues and eigenvectors of Eq. (1).For completeness, in SI 3 we also show the effect of including the intrinsic interaction from φ int (x), proving that its effect is small and that the main contribution stems from φ b .
In the following calculations, we consider the dielectric constants shown in Fig. 1(a): for the dielectrics materials (the wire, the substrate and the surrounding medium), we use typical values [37] ( = 17.7, d = 3.9 and a 1, respectively); while for the metallic leads we assume that, because they are bulky, they screen external electric fields perfectly (i.e.M → ∞).This may not be the case for the SC shell, whose capability for screening external electric fields may be weaker due to its small thickness and unavoidable presence of possible disorder [38].If this is the case, it is then characterized by a finite effective dielectric permittivity which depends on the SC shell width as well as its composition, as we show in SI 1.Some experiments [39] have reported that for ultrathin metallic layers (∼ 5 − 10nm) it is of the order of SC 100.For these values, as we show in the next section, we find a repulsive environment, i.e., an environment whose effective permittivity is smaller than that of the wire so that the bound charges that arise at the interfaces have in average the same sign as the free charges.We consider in SI 4 the generality of our results versus the SC dielectric constant and the location of the charge density within the nanowire section.In Fig. 4(c) we show that, when the charge density is fixed at the center of the wire, as SC becomes larger the dielectric environment turns into an attractive one and the pinning effect is eventually lost.This however strongly depends on the location of the charge density.If, as pointed out in Ref. [40], it happens to be close to the SC shell, the screening effect is larger and the pinning is suppressed.Nevertheless, as we analyze in Fig. 4(e), even if SC , the pinning effect survives when the wave function is located further away from the SC.

A. Results without bulk normal leads
It is convenient to start by analyzing the simpler case in which we neglect the effect of the bulk normal leads in the induced potential φ b .As an example we consider a nanowire of width W = 100nm, length L = 1µm and the following choice of realistic parameters: m * = 0.015m e , α = 20nm•meV, ∆ = 0.3meV, µ = 0.5meV and T = 10mK.These could correspond, for example, to a InSb  nanowire in contact to an Al superconducting shell [14], but similar results are obtained for InAs wire parameters [19].For an infinite wire, a schematic representation of the energy bands is shown in Fig. 2(a) in the absence and in the presence of the Zeeman field.At zero temperature, the occupied states below the Fermi level are those between the horizontal dashed line and the band bottom.Apart from a small contribution coming from the spin-orbit energy, the position of the band bottom is controlled by the wire's chemical potential µ, the Zeeman energy V Z and the induced potential energy eφ b .The magnetic field lowers the band bottom, charging the wire, whereas the induced potential energy, coming from electrostatic repulsion, tends to compensate that trend.
In the finite-length wire, the evolution of the induced potential profile along the nanowire length (x-axis) for different Zeeman fields is shown in Fig. 2(b).As can be observed, the induced potential tends to expel charge from the center of the wire, where it is positive, while it bends downwards at its ends.On the other hand, the evolution of the potential with Zeeman field exhibits a step-like behavior with regions where it increases linearly with V Z (red curves), screening the magnetic field effects, and regions where it remains almost constant as V Z increases (grey curves).This causes the electron fluid to behave in an incompressible or compressible manner, respectively.This different behavior can be clearly seen in Fig. 2(c) where the electrochemical potential at the center of the wire, given by V Z + µ − eφ b (L/2), is plotted against Zeeman splitting, both in the presence and absence of interactions.The effect of this peculiar evolution of the electrostatic potential has direct consequences on the wire spectral properties, as we analyze in Fig. 4, but for comparison, let us first see what happens in the non-interacting case.The wire's spectrum is shown in Fig. 3(a).There we can observe the emergence of low energy subgap states for V Z > ∆ 2 + µ 2 , corresponding roughly to the critical field for the bulk topological transition.We also obtain the typical energy oscillations produced by overlapping Majorana wave functions due to the finite length of the wire [21][22][23].More insight can be obtained by analyzing the evolution of the total wire's charge Q tot = L 0 dx ρ(x) as well as the Majorana charge Q M , whose absolute value is given by Here, Q ±1 are the charge corresponding to the even/odd lowest energy eigenstates ψ ±1 , and u L,R are the electron components of the Majorana wavefunctions . The total charge increases in general with magnetic field but, for finite length wires, it does so by jumping abruptly a quantity equal or smaller than e at each parity crossing (where the Majorana oscillations cross zero energy and the electron parity of the wire changes from even to odd or viceversa), as shown in Fig. 2(d), dashed curve.This abrupt change in charge is actually injected into the fermion state created by the two overlapping Majoranas and is given by |Q M | at the parity crossings.The (oscillatory) evolution of |Q M | with magnetic field is given in Fig. 3(b).Strikingly, |Q M | is maximum at the parity crossing, where the energy is zero, and goes to zero at the oscillation cusps.
As the length of the wire tends to infinity, Q M tends to zero (not shown).Indeed, the finite value of Q M at the parity crossings is a direct measurement of the Majorana overlap, as shown in Ref. [33].Note that the Majorana overlap is defined similarly to the right-hand side of Eq. ( 3), but with the absolute value inside the integral.The behavior of the Majorana wavefunctions is illustrated in Figs.3(c-f).The probability density for the left and right Majorana wave functions exhibits an overall decay towards the center of the wire controlled by the length ξ ∼ v F /∆ and an oscillatory pattern controlled by λ F [41,42].Moreover, the number of oscillations that fit in L increases by one with Zeeman field at each parity crossing.Interestingly, we observe that the left-right oscillatory patterns are out of phase for the cases where the splitting of the MBSs is maximum, panels (c) and (e).This minimizes the left-right wave function overlap and the Majorana charge goes to zero.On the other hand, the oscillations are in phase (d) when the energy splitting is zero, at the parity crossings, producing a maximum in |Q M | and overlap.Although the Majorana wave functions are more strongly located at the wire edges, we note that the charge density of this fermionic state is uniform across the wire [34] and, thus, it is uniformly affected by the interaction with the environment when this is present.
In the presence of interactions with the image charges, the single-point parity crossings versus V Z in the spectrum are replaced by extended regions where the subgap states remain pinned at zero energy, indicated by the red lines in Fig. 4(a).The abrupt jumps in Q tot in the noninteracting case are replaced by a linear increase along the V z -values where the zero-energy pinning occurs, see Fig. 2(d).This is a consequence of the repulsive environment that inhibits the entrance of charge in the wire where the electron liquid behaves in an incompressible manner.On the other hand, the Majorana charge remains basically constant at the pinning plateaux, as shown in Fig. 4(b).The finite value of Q M in these regions indicates that zero-energy does not imply absence of overlap between the left and right Majorana states.This is actually a common misconception that we would like to point out here.The Majorana overlap, which is a measurement of the degree of non-locality of the two Majorana wave functions, mostly depends on the length of the nanowire (and to a lesser extent on other parameters, such as the induced superconductor gap and the Rashba coupling), but it is not necessarily correlated to the Majorana energy splitting.Different mechanisms can reduce this splitting, such as interactions with the environment as studied here, smooth potential or gap profiles [21,[26][27][28]30], or orbital magnetic effects [31], and still leave the Majorana overlap unaffected.The behavior of the Majorana wavefunctions in this case is illustrated in Figs.4(c-f).In the pinning regions the Majorana wave functions remain practically frozen and in phase.This in turn explains why |Q M | is maximum in these regions.
The generality of these results with wire parameters is analyzed in the SI 4.There we show how the width of the pinning plateau evolves with V Z when we change the chemical potential, the dielectric permittivity or the width of the SC shell, and the aspect ratio of the nanowire section.We find that pinning survives for any chemical potential, while it vanishes when the attractive contribution of the SC shell becomes dominant over the dielectric repulsion.

B. Effect of bulk normal leads
In this section we analyze the effect of including the bulk normal leads in the calculation of the induced potential φ b .Fig. 5(a) illustrates the evolution of φ b with increasing Zeeman field for the same choice of parameters as in the previous section but including the normal contacts.As can be observed, while in the central region of the wire a similar repulsive step-like evolution with V z is found (corresponding to compressible/incompressible electron fluid behavior), significant attractive regions appear at the wire ends produced by the metallic character ( M → ∞) of the adjacent leads.As we discuss below, these attractive regions give rise to the formation of quantum-dot (QD) like bound states that may interact with the low energy subgap states of the Majorana wire.The evolution of the spectral properties and of the total charge Q tot in this case are shown in Fig. 5(c) and (e), respectively.On the one hand, we observe that the pinning plateaus around each parity crossing (in red) are still present although with a smaller width.On the other hand, the main effect of the presence of the attractive potential regions is the appearance of an additional four energy levels (two per contact, in blue) that approach zero energy for V Z around 2.5 − 3meV.At the same time we observe a rather abrupt decrease in the total wire charge (of roughly 2e), see Fig. 5(e).We can associate these additional levels with QD-like bound states arising in the attractive regions of the induced potential that anticross with the Majorana levels when their energies are on resonance [35,36,43,44].
To demonstrate the validity of this interpretation we show in Fig. 5(d) and (f) the spectral properties and the total charge evolution for an isolated wire with a simple double potential well taken to mimic the effect of the electrostatic environment, shown in (b).Notice that in this case we do not attempt a self-consistent calculation but rather include the Zeeman field as a rigid shift of the two spin bands (like in the non-interacting case but with an inhomogeneous potential profile).Although the zero-energy pinning is not captured by this model, one can clearly observe the presence of four levels coming down towards zero-energy for V Z around 2.5 − 3meV, as in the interacting case.The presence of these states is a consequence of the renormalization of the topological phase transition due to the electrostatic potential (either φ b or φ f ixed ) which is not constant along the wire because φ b (or φ f ixed ) depends on x.For the shown values of V Z , only the central part of the wire is in the topological regime (V Z > V c Z ), corresponding to an effectively shorter Majorana wire, whereas the outer parts are trivial (V Z < V c Z ), corresponding to two effective QDs attached to it.Specific details of how QD-levels interact with Majorana nanowire ones can be found in Ref. [35,36,[43][44][45].
Further evidence on the nature of the low energy states at V Z ∼ 3meV is provided in Fig. 6 where we plot the wave-function probability profiles (in the Majorana basis) of the low energy states around the QD-Majorana levels anticrossing.For simplicity, we consider only the case of the potential barrier model.At the anticrossing, the Majorana and dot states merge and cannot be really told apart, but we will refer to the two lowest in energy as Majorana levels and to the other two as dot levels.As can be observed, at the anticrossing the Majorana levels (green circle) leak into the QD regions leaving the central (topological) part practically void.Conversely, the two dot-like states (immediately above in energy, orange squares) penetrate and delocalize along the wire.When the Zeeman field increases and the QD and Majorana levels are detuned, the dot states depart from low energy (pink rhombus) and from the topological part of the wire, Majorana nanowire subject to interactions from the electrostatic environment (including the influence of the bulk normal leads at its ends).(a) Self-consistent induced potential energy eφ b (x) along the wire's length for increasing values of the Zeeman splitting.Same wire parameters as in Fig. 2. Note that the main effect of the bulk normal leads is to create confining potential wells at the wire edges.(b) Barrier-like potential energy profile used to mimic the self-consistent solution.Spectra of the Majorana nanowire as a function of VZ in the (c) interacting case and in the (d) non-interacting case but using the potential profile model of (b).(e,f) Evolution of the total charge Qtot with Zeeman splitting for the two previous cases, respectively.Red color indicates incompressible electron fluid behavior as before, while blue color indicates QD-like behavior due to the metallic contacts.
whereas the usual overlapping behavior of the MBSs is recovered but with the Majoranas bound to the effective topological edges (yellow triangle).The absolute value of the Majorana charge vs. V Z is shown in Fig. 6(b), calculated considering only the two lowest energy states (as before).At the anticrossing the Majorana charge oscillation is distorted, see blue region, but the area below the curve is conserved.The missing charge in Fig. 5 (b) does not come from the Majorana states, but from the dot ones.At the anticrossing region, the two QD states (one per potential well) that were occupied (below the Fermi level) move upwards in energy as the Zeeman field increases and cross the Fermi level, emptying themselves.This is why in the blue regions of Fig. 5 (e) and (f) the wire's total charge does not increase at the corresponding parity crossing, but instead decreases loosing effectively twice the charge of the electron e.
Finally, we would like to point out that, when the dot levels anticross the Majorana ones in a pinning region, the Majorana states detach from zero energy.This can be seen in Fig. 5(c) and Fig. 1(b).The reason is that, although in the pinning regions the Majorana energy is zero, their wave function overlap is not.It is actually maximum, as explained when discussing Fig. 4.Each QD acts as a local probe (one couples to the left topological region of the wire, the other to the right).If the wire's length were large (much bigger than the coherence length), left and right Majoranas would be disconnected from each other, and a local probe coupled to one of them would not be able to change its energy or perturb it.This is actually the core manifestation of their topological protection.However, when the wire's length is finite and the Majoranas overlap, each QD couples to both Majoranas at either end and their energies are modified.The typical shapes of the anticrossing were recently analyzed and can be used to quantify the degree of Majorana non-locality [35,36].

IV. CONCLUSIONS
In this work we have studied the low energy characteristics of Majorana nanowires when including its interaction with a realistic 3D electrostatic environment.This is done by solving self-consistently the Bogoliubovde Gennes equation together with the Poisson's equation.Typically, the total charge of the wire in equilibrium with the reservoirs increases with magnetic field (or the wire's chemical potential).However, if the electrostatic screening is smaller inside the wire than at the contacts, a repulsive interaction arises that leads to zero energy pinning around parity crossings in the wire's spectrum.While the screening due to the parent SC shell tends, in general, to reduce this pinning effect, we find that it still persists depending on the quality of the SC layer and the location of the charge density within the nanowire.The pinning mechanism could help explain the precise shape of the Majorana oscillations (or lack thereof) observed in some dI/dV experiments, which exhibit substantial deviations from the predictions of simple models for finite length wires.
On the other hand, and more importantly, the self consistent solution of the electrostatic potential varies non-homogeneously along the wire.It is relatively flat in the central region but, due to the screening from the left/right metallic contacts, it becomes strongly negative at the edges.This creates potential wells that confine QD-like states at the ends of the wire, which appear in the spectrum as discrete states within the induced gap that disperse with Zeeman energy or chemical potential.These QD levels interact with the Majorana states in a specific way which is strongly dependent on the Majorana wavefunction, and particularly on its degree of spatial non-locality.The pinning mechanism and the coupling to QD-like states compete against each other, so that the pinned zero-energy plateaux may become lifted at resonance with the dot states, thus revealing their electrostatic origin (as opposed to true wavefunction nonlocality).

Introduction
The electrostatic potential φ b produced by the environment is calculated from the interaction between the nanowire charge density and bound charges it creates at the surrounding medium shown in Fig. (1) of the main text.To do this we use the method of the image charges.We model the nanowire as a semiconductor rod of square section of relative permittivity , length L and rectangular section of width W = 2R, being R the half-width.We assume that the charge density ρ ( r) of the nanowire is located along its symmetry axis (xaxis) as a linear charge density.The nanowire faces are in contact with different dielectric materials of permittivities 1 (substrate), 2 (superconducting shell), 3 and 4 (surrounding medium); while two metal leads of permittivities M1 and M2 are placed at both ends of the nanowire.In order to give more insight on the solution of this problem, we solve first some simpler cases.In the first section we obtain the electrostatic potential created by a linear charge density placed before one, between two, and before two infinite planes.These problems can be understood as 1D problems.In the second section we obtain the potential with all four surrounding media but without the bulk leads, which can be treated as a 2D problem.Finally, in the third section we obtain the full model including also the interaction with the leads.

A. One infinite plane
The solution of this problem can be found in textbooks on electromagnetism [46].When one charge q is placed (at the origin of coordinates) in a medium of dielectric constant and at a distance R from the interface between this and another medium of permittivity 1 , a bound charge of magnitude κ 1 q appears spread at the interface, where The effect of this bound charge is equivalent to that of a point charge of the same magnitude and located at a specular distance with respect to the plane from the original charge, see Fig. 7(a).This is why it is called an image charge.The classical electrostatic potential due to the interaction between the original and the image charge through the Coulomb's law takes the form where 0 is the vacuum permittivity.Because φ b is linear in q, this result can be generalized to an arbitrary 1D density charge ρ (x): As argued in Ref. [33], since the bound charges are distinguishable from the nanowire charges (cannot tunnel in and out), this potential can be directly transformed into a quantum operator as a Hartree interaction without any Fock correction.Assuming a purely local polarizability (Thomas-Fermi limit), the density operator may be transformed as ρ (x ) → ρ (x ) , which is perfectly equivalent to the above classical equation.Then, the potential takes the form where the kernel of the interaction V b (x, x ) ≡ 2 encodes the geometrical information of the interaction.We note that Eq. ( 8) is general for any geometry and charge density.For this reason, in the following sections we first obtain the kernel function V b for a given geometry using a single original charge, and then we generalize our results to an arbitrary density ρ (x) and to its corresponding quantum expression using Eq. ( 8).

B. Two opposite infinite planes
We now consider a charge between two infinite and parallel planes separating the nanowire's dielectric constant from two other media of permittivities 1 and 3 [see Fig. 7(b)].In order to try to satisfy the boundary conditions imposed by the Gauss' law, in a first step two image charges κ 1 q and κ 3 q are required in each dielectric medium at the same distance R from the interface [see red dots in Fig. 7(b)].However, each image charge does not satisfy the boundary condition with respect to the opposite interface.For this reason, in a second step, two (accidentally equals) additional image charges of magnitude κ 1 κ 3 q are required at a distance 3R from each interface [see purple dots in Fig. 7(b)].But, again, these additional image charges do not satisfy the boundary conditions with respect to the opposite interfaces and another step must be taken.It is possible to see that for each image charge q (n) α in the dielectric α created at the n-th step of the image charge method, another image charge is required in the opposite β dielectric at a distance 2 (n + 1) R from the original charge q in order to satisfy the boundary conditions, with q (0) α,β = 1.Thus, the interaction kernel takes the form of an infinite series Note that each term of the sum decreases with n as V (n) b ∼ κ n /n.Since |κ| < 1 for a dielectric medium and κ M = −1 for a metallic one (where M → ∞), then the kernel converges to V b ∼ ln (1 − κ) in spite of the infinite summation.

C. Two consecutive infinite planes
Finally, we consider the case depicted in Fig. 8(a) where a (real) charge q is placed in a dielectric material characterized by a permittivity A .This charge is at a distance R from a first interface with a material of dielectric constant B and width W , and at a distance R + W from a second (parallel) interface with another medium of permittivity C .As we know, the real charge q creates an image charge κ B q at a distance 2R from it inside the B medium.However, the potential created by both charges is only valid inside the A medium.The potential created inside the B medium can be found using the image method [46]: it is the same potential than that created by an image charge of magnitude q[2 A /( A + B )] located at the same position of the real charge.From this point, the same steps explained in the previous subsection can be followed, and once the infinite series is obtained, the potential can be transformed back to the A medium.Thus, one can prove that the bound charges potential is given by Since this expression is rather complex, it is convenient to replace the effect of both media B and C by just one characterized by an effective permittivity eff which, from an electrostatic point of view, is equivalent.Hence, the bound charges potential would be given by Eq. 6 with → A and 1 → eff .Comparing both equations, the effective permittivity (at x = 0) is where is in a medium of dielectric constant A (the nanowire) and placed at a distance R from a thin material of permittivity B and thickness W (the SC shell).Above the thin film there is another semi-infinite space of permittivity C (the vacuum).When all the permittivities are finite, bound charges arise in both surrounding mediums.(b) Effective SC permittivity SC calculated using Eq. 12 vs the SC thin film permittivity tf for two different film thicknesses WSC .Parameters are the same as in the main text.
We note that this system corresponds to the nanowire-SC shell-vacuum double interface of 1(a).There, the effect of the shell finite width and the vacuum on top has been condensed in an effective SC permittivity.This means that, in Eq. 12, eff → SC , A → , C → a and B is the true SC thin film permittivity tf .
In Fig. 8(b) we show the effective SC permittivity SC considered in the main text as a function of tf for two different shell widths W SC .Notice that, as the SC shell becomes thinner (red curve corresponds to 8nm), the effective permittivity gets more renormalized.A value SC ∼ 100 corresponds to tf ∼ 4000.
Section II: Two dimensions

A. One rectangular corner
This is also a textbook problem [46].A charge q inside a dielectric with permittivity is placed at a distance R 1 from dielectric 1 and at a different distance R 2 form dielectric 2 , which are perpendicular to one another [see Fig. 7(c)].To satisfy the boundary conditions, two image charges κ 1 q and κ 2 q are required in each dielectric at (x, −2R 1 , 0) and (x, 0, −2R 2 ), see red dots.Because of these, another image charge of magnitude κ 1 κ 2 q is required at (x, −2R 1 , −2R 2 ), purple dot.In this case, and due to the closed geometry of the problem, three image charges are enough to satisfy the boundary conditions.
The kernel function takes thus the form B. Four rectangular corners: interaction without the leads Now a charge q is placed inside an infinite wire of rectangular section and of permittivity , see Fig. 7(d).This charge is at a distance R 1 from two parallel flat dielectrics with permittivities 1 and 3 , and at a distance R 2 from another two parallel dielectrics with permittivities 2 and 4 which are perpendicular to the previous ones.Combining what we have learned in the previous Sections, to satisfy the boundary conditions an infinite ensemble of image charges have to be placed in the different dielectrics as shown in Fig. 7(d).For each image charge in one dielectric, another one appears at a specular distance from the opposite interface, while for each two image charges placed in two perpendicular media, just one more appears at the corner.The electrostatic potential due to all these charge is given by where: If the wire's section is square (R 1 = R 2 ), then the kernel function can be rewritten in a more compact way: − δ m,0 Note that the number of charges increases as 4n at each n-th step of the image charge method, whereas the rest of the expression inside the brackets decreases as κ n /n as before.Thus, each term of the sum goes as so the convergence of the kernel is ensured in this case as well.
Section III: The full-model Finally, we solve the full system of Fig. (1) of the main text.Now, apart from the four dielectric media in each of the four faces of the square section, there are another two faces in the x-direction in contact to metallic regions.We consider that the nanowire has a square section of semiwidth R. First, we assume that the charge q is placed at the coordinates origin and at the same distance R from each metallic region M 1 and M 2 .Following the same procedure as before we obtain where If now the charge q is placed at an arbitrary position x inside the nanowire, and the metal M 1 interface is at x = 0 and the M 2 interface is at x = L, then the kernel function is given by If L is large enough compared to 2R, we can take into account only the lowest order of the image charges at the metals, q (k=0,1) Mi , and the number of charges at each n-th step of the image charge method increases only as ∼ 12n.Then, the system follows the same convergence criterion as in the previous Section.If L ∼ 2R, Eq. ( 20) converges as well, but the demonstration is longer.

Mean field approximation to treat electron-electron interactions
We want to solve the energy spectrum of the nanowire Hamiltonian Ĥ when the interaction between electrons φ is included.In general, this interaction can be written in second quantization as where č † α , čα are defined as the Nambufied vector of operators Here, c † iσ and c iσ are electron creator/annihilation operators with quantum numbers α.V αβ above encodes the electronic interaction.Therefore, Greek indexes α, β encode both particle/hole character, spin (↑, ↓), and any other indexes the electron might have, such as site index i = 1, ..., N in a tight-binding description.
To treat the quartic interaction we resort to a mean field approach called the Hartree-Fock-Bogoliubov (HFB) approximation.Using the Wick's theorem and neglecting fluctuations and constant terms we can write φeff = α,β The first two terms are known as Hartree terms, which include information about direct (repulsive/attractive) interaction between electrons.The second and third are the Fock terms, which include the exchange interaction due to the electron indistinguishable properties.These terms ensure that non-physical self-interactions introduced by the first terms are cancelled.The last two are known as Bogoliubov terms, which include possible pairing correlations between electrons.We want to rewrite Eq. ( 23) in a more compact manner.In order to do that, we define two matrices: • The lambda matrix: where I space is the identity matrix in real space (for a one-dimensional tight-binding model with N sites, this matrix is the N × N identity), I 2×2 is the identity matrix in spin space, and σ x is the Pauli x-matrix in Nambu space.This matrix satisfies the property č † = Λč.
• The density matrix: We can express this matrix in terms of the eigenvectors γ n of the Hamiltonian Ĥ +e φeff , which are related to the čα 's through a unitary transformation Ψ as γ n = Ψ nα čα .Then where Using these two matrices, Eq. ( 23) can be rewritten as where we assume a symmetric potential V αβ = V βα .Here we have used the notation D [v] as the diagonal matrix with vector v in its diagonal, d {A} as a column vector whose elements are the diagonal elements of matrix A, the dot product A • B as a matrix product, and the star product A B as an element wise product between matrices.However, Eq. ( 27) does not have Nambu structure because in general V iστ,iστ = V iσ τ ,iσ τ , so Bogoliubov-de Gennes formalism cannot be applied.We symmetrize this expression by doing where we have used the anticommutation relation c, c † = 1, the property Λ t = Λ and we neglect constant terms again.Thus, the general interaction between electrons in the HFB approximation can be written as where φ eff is given by Eq. ( 27).

Inclusion of the intrinsic interaction
The interaction between the electrons inside the nanowire (intrinsic interaction) is given, in principle, by the bare Coulomb potential in one dimension.Taking into account the finite radius (half-width) R of the wire, a more precise form for the interaction is [47] V where x is the distance between electrons.When we consider this potential only at the Hartree level, we find zero energy pinning around parity crossings, just like we did with the extrinsic interaction.However, this pinning is unphysical since it comes from spurious self-interaction terms introduced by the Hartree approximation [33].For the intrinsic interaction it is thus necessary to include the Fock correction due to the indistinguishability of electrons in the nanowire.
If we consider the bare interaction, Eq. ( 30), in the Fock terms, an overcompensation of the pinning effect is found and unphysical jumps appear at each parity crossing.To cure this problem, we introduce screening in the quasi-static Thomas-Fermi limit so that the potential in the Fock terms acquires an additional exponential decay e −|x|/λ T F that depends on the Thomas-Fermi length λ T F (which should be of the order of the Fermi wavelength).
In this case, we find that the parity crossings induced by the intrinsic interaction are suppressed, in agreement with the self-interaction argument.
If the nanowire is discretized using a thigh-binding model, the interaction can be written as where the indexes α = {i, σ, τ } and β = {j, σ , τ } include all the quantum numbers of the electrons, a is the distance between two neighboring sites (lattice constant), and the term [1 − δ α,β ] ensures that an electron cannot interact with itself.As stated before, the above equation is only valid for the Fock terms, while for the Hartree terms it is the bare interaction (same expression with λ T F → ∞).Then, one can obtain the potential in the HFB approximation using Eq. ( 29).
The electron-electron interaction between the nanowire and the bound charges of the dielectric environment (extrinsic interaction) is found in SI 1, and it may be implemented following the same procedure by substituting x → (i − j) /a.We note that now the term [1 − δ α,β ] should not be included since electron α is always inside the nanowire while β is outside, in the surrounding medium (or the other way around).Finally, the potential in the HFB approximation can be computed using Eq. ( 29), but now the Fock and the Bogoliubov terms have to be ignored as we argued in the SI 1, so that the last four terms of Eq. ( 27) are not considered.

Numerical self-consistent method used to solve the eigenspectrum
We note that Ĥ + e φHFB depends on its own eigenvectors (see Eq.( 26)), and thus it has to be solved selfconsistently.We solve this problem numerically using the following procedure: in the first iteration of the selfconsistent method, we find the density-matrix ρ using the eigenvectors of the Hamiltonian Ĥ as a test solution.In the next step we obtain a new ρ diagonalizing Ĥ + e φHFB where the interaction has been obtained using the previous density matrix.In the following steps, the density-matrix is found using a linear combination of the eigenvectors of Ĥ + e φHFB in the two previous iterations.This is done to introduce damping in the iteration procedure in order to ensure the convergence of the self-consistent method.In each step, we compare the eigen-energies of Ĥ + e φHFB with those of the previous step.We repeat this procedure until convergence.We consider the iteration has converged when the difference between both energy-spectra is much smaller than the main energy scale of our problem (i.e. the superconductor gap ∆).Here we show that the features studied in the main text (zero-energy pinning and QD formation) remain qualitatively unaltered when including electron-electron interactions φ int inside the nanowire.The intrinsic interaction introduces small quantitative changes in the spectrum, but the qualitative behavior stays the same.Following our previous work [33], we treat this interaction at the mean field level, within the Hartree-Fock-Bogoliubov approximation, and assume a bare Coulomb interaction for the Hartree terms and a screened Coulomb interaction in the quasi-static Thomas-Fermi limit for the Fock terms (see SI 2 for more details).
In Fig. 9 we show the bound charges electrostatic potential along the nanowire and the energy spectrum versus the Zeeman field ignoring (a and c) and including (b and d) the intrinsic interaction.Note that here we do not include the Bogoliubov correction since this term just renormalizes the gap as shown in Ref. [33].In general, both spectra are qualitatively similar and thus we can conclude that the intrinsic interactions do not alter the features studied in the main text.However, there are some quantitative differences.First, the dispersive QD levels approach zero energy at a slightly smaller Zeeman energy as a result of small changes in the induced potential φ b , as can be seen in Fig. 9(d).Second, the position of the gap closing and thus, the topological phase transition, shifts to a different magnetic field.This is also a consequence of small changes in φ b , as well as small changes in the Zeeman energy induced by the Fock terms, that modify the value of the critical Zeeman field through Eq. ( 4) of the main text.Finally, the energy splitting of the Majoranas is larger due to the renormalization of the Fermi momentum induced by the Fock terms as well.

SI 4: ROBUSTNESS OF THE PINNING EFFECT
We want to test the validity of our results when varying different parameters of the electrostatic environment.Fig. 10 provides various phase-diagrams indicating the occurrence of Majorana bound states zero-energy pinning (in red) as a function of the different parameters.Although we have taken µ = 0.5meV for the simulations of the main text, pinning is general for all chemical potentials (within the topological phase), as can be seen in the upper panels of Fig. 10.As a function of µ and V Z , the non-interacting lines of (a) corresponding to pointlike parity crossings transform into incompressible finite width stripes in the interacting case (b).Pinning regions are bigger for lower chemical potentials and for higher magnetic fields, since the repulsive interaction is larger too.It can also be observed that the onset of the topological phase is different in the interacting system than in the non-interacting one (black dashed line), at least for positive µ.This is because the electrostatic potential renormalizes the chemical potential [48] and thus it modifies the value of the bulk critical magnetic field, given in Eq. ( 4) of the main text.
However, pinning is not general for all kind of environments.Figure 10(c) shows the zero energy regions across the V Z − SC space (the µ − SC diagram exhibits a similar behavior).For SC 300 the pinning plateau width shrinks into points because the electrostatic environment turns into an attractive one.This means that bound charges of the opposite sign arise in the dielectric medium at these large permittivities, so that the entrance of charge at each parity crossing is no longer suppressed.Note that SC represents the effective SC permittivity of the system composed by a SC thin film (epitaxially grown over the nanowire) of intrinsic permittivity tm , finite width W SC and covered by vacuum, as we argue in SI 1(IC).For a film width of 8nm, an effective SC 300 corresponds to a SC permittivity of SC 15000, i.e., basically a perfect metal.
The width of the nanowire also plays a role in the appearance of pinned regions.Fig. 10(d) shows the incompressible regions as a function of V Z and r yz , where r yz = W y /W z is the aspect ratio of the nanowire section.When the distance between the SC shell and the opposite side is large (large r yz ), pinning is bigger.This is because the relative coverage of the wire by the SC shell decreases and so does its attractive contribution.
Finally, if we consider perfect metallic screening by the SC shell, i.e., SC → ∞, we can also study the appearance of the pinned regions depending on the distance between the nanowire charge density and the SC shell.In Fig. 4(e) we study the phase diagram as a function of V Z and the position of the charge density across the nanowire section, y/R, where R is the wire's half width.In Fig. 4(f) the position is fixed to the center of the wire, but we vary the wire's aspect ratio.In both cases we observe that, as the Majorana wave functions approaches the SC, the screening effect increases and the pinning disappears [40].However, if the wave function separates from the SC shell, the pinning survives.This may happen when the bottom gate attracts the charge density away from the SC or when the wave function is more spread throughout the wire's section, as for example for higher sub-bands.

FIG. 1 .
FIG. 1.(a) Schematic representation of the setup analyzed in the present work.A nanowire of rectangular cross-section (green) lying on an insulating substrate (grey) and in contact to a thin metallic layer in one of its faces (light blue), corresponding to the parent superconductor, and two normalmetal leads at its ends (orange) separated by tunnel barriers (brown).Typical values for the dielectric constants for each region are indicated.(b) Low energy spectrum versus chemical potential µ for a wire of thickness W = 100nm and length L = 1µm.Other parameters are the spin-orbit coupling α = 20nm•meV, the induced pairing energy ∆ = 0.3meV and the Zeeman energy VZ = 2meV.Electrostatic environmentinduced zero energy pinned regions between Majorana oscillations are indicated in red.Quantum dot levels (in blue), originated at the wire's edges due to the interaction with the bulk contacts, anticross with Majorana levels and remove their zero-energy pinning.

FIG. 2 .
FIG.2.Majorana nanowire subject to interactions from the electrostatic environment (ignoring the influence of the bulk normal leads at its ends).(a) Schematic of the nanowire's dispersion relation in the absence and in the presence of the Zeeman field.(b) Self-consistent induced potential energy eφ b (x) along the wire's length for increasing values of the Zeeman splitting.Wire parameters as in Fig.1(b) and with µ = 0.5meV.(c) Energy difference between the Fermi level and the band bottom at the center of the nanowire, VZ + µ − eφ b (L/2), and (d) total charge Qtot of the nanowire as a function of VZ for the non-interacting (dashed) and interacting (solid line) cases.Red curves highlight parameter regions for which there is interaction-induced zero-energy pinning in the spectrum.

FIG. 3 .
FIG.3.Majorana wave functions in the non-interacting case: Energy levels (a) and the absolute value of the Majorana charge QM (b) vs Zeeman energy.Panels (c-f) show the wave-function probability profiles of the two lowest energy states in the Majorana basis at selected values of the Zeeman field within the topological region.When the splitting is maximum (green circles and yellow triangles) the left and right Majorana wave-function oscillations are out of phase, whereas when the splitting is zero (orange square) they are in phase.

FIG. 4 .
FIG.4.Same as Fig.3but for the interacting case (without leads).In the pinned regions the Majorana wave functions remain in-phase as a function of Zeeman field and the Majorana charge (b) freezes at its local maximum value (in red), instead of continuing the oscillation like in the non-interacting case (dashed curve).

FIG. 6 .
FIG. 6. Evolution with Zeeman field of the spectrum (a) and the absolute value of the Majorana charge QM (b) for the barrier-like potential model of Fig. 5(b).Panels (c) and (e) show the wave-function probability profile (in the Majorana basis) of the two lowest energy states at the VZ values indicated in (a).Panels (d) and (f) show the same but for the second and third energy states (QD-like states).At the dot-Majorana levels anticrossing, the Majorana wave function leaks into the dot regions leaving the topological region of the wire practically void.This is manifested in |QM | by two consecutive zeros, one per dot level (around VZ = 2.5meV).

FIG. 7 .
FIG.7.Configuration of the image charges produced by the original charge (in black) studied in different Sections.The original charge is in a medium of dielectric constant (green color in the plots) that extends in the third direction.(a) Sec.IA: a charge in front of an infinite semi-space of permittivity 1 at a distance R. (b) Sec.IB: a charge between two semiinfinite and parallel regions of permittivities 1 and 3. (c) Sec.IIA: a charge close to the intersection between two media of permittivities 1 and 2. (d) Sec.IIB: a charge between four different media.As before, the black point depicts the real charge, while the red, purple and orange points depict the first, second and third image charges found in the first, second and third steps of the image method procedure, respectively.For (a) and (c) all the image charges are shown, whereas for (b) and (d) there is an infinite number of them and only the first ones are shown.Black solid lines are interfaces, and black dashed lines are image interfaces.

FIG. 8 .
FIG.8.(a) System studied in Sec.IC.A charge (black dot) is in a medium of dielectric constant A (the nanowire) and placed at a distance R from a thin material of permittivity B and thickness W (the SC shell).Above the thin film there is another semi-infinite space of permittivity C (the vacuum).When all the permittivities are finite, bound charges arise in both surrounding mediums.(b) Effective SC permittivity SC calculated using Eq. 12 vs the SC thin film permittivity tf for two different film thicknesses WSC .Parameters are the same as in the main text.

FIG. 9 .
FIG. 9. Majorana nanowire in the presence of interactions (including the influence of the bulk normal leads at its ends).Self-consistent induced potential energy eφ b (x) along the wire's length for increasing values of the Zeeman splitting ignoring (a) and including (b) the electron-electron interactions inside the nanowire.(c) and (d) are their corresponding energy spectra.Wire parameters are the same as in the main text, and the Thomas-Fermi length is λT F = L/3.

FIG. 10 .
FIG. 10.Phase diagrams indicating the parameter regionswhere the Majorana bound states are pinned to zero-energy (in red).In the upper panels the phase diagram is calculated as a function of Zeeman field and chemical potential for the non-interacting case (a) and interacting case(without leads) with SC = 100 (b).The central panels correspond both to the interacting case, but (c) considers variations in the effective dielectric constant of the thin superconducting layer, SC , whereas (d) explores different aspect ratios of the nanowire's section ryz = Wy/Wz, where Wy,z are the y and z widths of the nanowire faces.In the lower panels we consider perfect screening by the SC shell, SC → ∞, but we vary the distance between the transversal Majorana charge density and the SC shell.In (e) y/R is the position of the Majorana wave function across the nanowire section.When y = 0 the charge density is at the center of the nanowire, whereas when y = R it is a the opposite face of the SC.In (f) the wave function is fixed at the center of the nanowire section, but we vary its aspect ratio.Here µ = 0.5meV and Wz = 100nm, as in the main text.