Effects of the electrostatic environment on the Majorana nanowire devices

One of the promising platforms for creating Majorana bound states is a hybrid nanostructure consisting of a semiconducting nanowire covered by a superconductor. We analyze the previously disregarded role of electrostatic interaction in these devices. Our main result is that Coulomb interaction causes the chemical potential to respond to an applied magnetic field, while spin-orbit interaction and screening by the superconducting lead suppress this response. Consequently, the electrostatic environment influences two properties of Majorana devices: the shape of the topological phase boundary and the oscillations of the Majorana splitting energy. We demonstrate that both properties show a non-universal behavior, and depend on the details of the electrostatic environment. We show that when the wire only contains a single electron mode, the experimentally accessible inverse self-capacitance of this mode fully captures the interplay between electrostatics and Zeeman field. This offers a way to compare theoretical predictions with experiments.


I. INTRODUCTION
Majorana zero modes are non-Abelian anyons that emerge in condensed-matter systems as zero-energy excitations in superconductors. [1][2][3] They exhibit non-Abelian braiding statistics 4 and form a building block for topological quantum computation. 5 Following theoretical proposals, 6,7 experiments in semiconducting nanowires with proximitized superconductivity report appearance of Majorana zero modes signatures. [8][9][10][11][12] These "Majorana devices" are expected to switch from a trivial to a topological state when a magnetic field closes the induced superconducting gap. A further increase of the magnetic field reopens the bulk gap again with Majorana zero modes remaining at the edges of the topological phase.
Inducing superconductivity requires close proximity of the nanowire to a superconductor, which screens the electric field created by gate voltages. Another source of screening is the charge in the nanowire itself that counteracts the applied electric field. Therefore, a natural concern in device design is whether these screening effects prevent effective gating of the device. Besides this, screening effects and work function differences between the superconductor and the nanowire affect the spatial structure of the electron density in the wire. The magnitude of the induced superconducting gap reduces when charge localizes far away from the superconductor, restricting the parameter range for the observation of Majorana modes.
To quantitatively assess these phenomena, we study the influence of the electrostatic environment on the properties of Majorana devices. We investigate the effect of screening by the superconductor as a function of the work function difference between the superconductor and the nanowire, and we study screening effects due to charge. We focus on the influence of screening on the behavior of the chemical potential in the presence of a magnetic field in particular, because the chemical potential directly impacts the Majorana signatures.
The zero-bias peak, measured experimentally in Refs. 8-12, is a non-specific signature of Majoranas, since similar features arise due to Kondo physics or weak antilocalization. 13,14 To help distinguishing Majorana signatures from these alternatives, we focus on the parametric dependence of two Majorana properties: the shape of the topological phase boundary 15,16 and the oscillations in the coupling energy of two Majorana modes. [17][18][19][20][21] Both phenomena depend on the response of the chemical potential to a magnetic field, and hence on electrostatic effects. Majorana oscillations were analyzed theoretically in two extreme limits for the electrostatic effects: constant chemical potential [19][20][21] and constant density 20 (see App. A for a summary of these two limits). In particular, Ref. 20 found different behavior of Majorana oscillations in these two extreme limits. We show that the actual behavior of the nanowire is somewhere in between, and depends strongly on the electrostatics.

A. The Schrödinger-Poisson problem
We discuss electrostatic effects in a device design as used by Mourik et al, 8 however our methods are straightforward to adapt to similar layouts (see App. B for a calculation using a different geometry). Since we are interested in the bulk properties, we require that the potential and the Hamiltonian terms are translationally invariant along the wire axis and we consider a 2D cross section, shown in Fig. 1. The device consists of a nanowire with a hexagonal cross section of diameter W = 100 nm on a dielectric layer with thickness d dielectric = 30 nm. A superconductor with thickness d SC = 187 nm covers half of the wire. The nanowire has a dielectric constant r = 17.7 (InSb), the dielectric layer has a dielectric constant r = 8 (Si 3 N 4 ). The device has two electrostatic boundary conditions: a fixed gate potential V G set by the gate electrode along the the lower edge of the dielectric layer and a fixed potential V SC in the superconductor, which we model as a grounded metallic gate. We set this potential to either V SC = 0 V, disregarding a work function difference between the NbTiN superconductor and the nanowire, or we assume a small work function difference 22,23 We model the electrostatics of this setup using the Schrödinger-Poisson equation. We split the Hamiltonian into transverse and longitudinal parts. The transverse Hamiltonian H T reads with x, y the transverse directions, m * = 0.014m e the effective electron mass in InSb (with m e the electron mass), −e the electron charge, and φ the electrostatic potential. We assume that in the absence of electric field the Fermi level E F in the nanowire is in the middle of the semiconducting gap E gap , with E gap = 0.2 eV for InSb (see Fig. 2 with z the direction along the wire axis, α the spin-orbit coupling strength, E Z the Zeeman energy and σ the Pauli matrices. The orientation of the magnetic field is along the wire in the z direction. In this separation, we have assumed that the spin-orbit length l SO = 2 /(m * α) is larger or comparable to the wire diameter, l SO W . 24,25 Furthermore, we neglect the explicit dependence of the spin-orbit strength α on the electric field. We ignore orbital effects of the magnetic field, 26  of the transverse wave functions is much smaller than the wire cross section due to screening by the superconductor, as we show in Sec. III. Since the Hamiltonian is separable in the limit we are using, the charge density in the transverse direction ρ(x, y) is: with ψ i the transverse wave function and E i the subband energy of the i-th electron mode defined by H T ψ i = E i ψ i . Further, n(E i , E Z , α) is the 1D electron density, which we calculate in closed form from the Fermi momenta of different bands in App. C. The subband energies E i depend on the electrostatic potential φ(x, y), and individual subbands are occupied by "lowering" subbands below E F (shown schematically in Fig. 2(b)). 27 The Poisson equation that determines the electrostatic potential φ(x, y) has the general form: with the dielectric permittivity. Since the charge density of Eq. (3) depends on the eigenstates of Eq. (1), the Schrödinger and the Poisson equations have a nonlinear coupling. We calculate the eigenstates and eigenenergies of the Hamiltonian of Eq. (1) in tight-binding approximation on a rectangular grid using the Kwant package. 28 We then discretize the geometry of Fig. 1 using a finite element mesh, and solve Eq. (4) numerically using the FEniCS package. 29 Eqs. (1) and (3) together define a functionalρ[φ], yielding a charge density from a given electrostatic potential φ. Additionally, Eq. (4) defines a functionalφ[ρ], giving the electrostatic potential produced by a charge density ρ. The Schrödinger-Poisson equation is self-consistent when We solve Eq. (5) using an iterative nonlinear Anderson mixing method. 30 We find that this method prevents the iteration process from oscillations and leads to a significant speedup in computation times compared to other nonlinear solver methods (see App. E). We search for the root of Eq. (5) rather than for the root of since we found Eq. (5) to be better conditioned than Eq. (6). The scripts with the source code as well as resulting data are available online as ancillary files for this manuscript.

B. Majorana zero modes in superconducting nanowires
Having solved the electrostatic problem for the normal system, i.e. taking into account only the electrostatic effects of the superconductor, we then use the electrostatic potential φ(x, y) in the superconducting problem. To this end, we obtain the Bogoliubov-de Gennes Hamiltonian H BdG by summing H T and H L and adding an induced superconducting pairing term: with τ the Pauli matrices in electron-hole space and ∆ the superconducting gap. The three-dimensional BdG equation (7) is still separable and reduces for every subband with transverse wave function ψ i to an effective one-dimensional BdG Hamiltonian: where p = −i ∂/∂z and we defined µ i = −E i (see Fig. 2(b). Since the different subbands are independent, µ i can be interpreted as the chemical potential determining the occupation of the i-th subband. While the Fermi level is kept constant by the metallic contacts, the chemical potential µ i of each subband does depend on the system parameters: µ i = µ i (V G , E Z ). Most of the model Hamiltonians for Majorana nanowires used in the literature are of the form of Eq. (8) (or a twodimensional generalization) using one chemical potential µ. To make the connection to our work, µ should be identified with µ i , and not be confused with the constant Fermi level E F . For example, the constant chemical potential limit of Ref. 20 refers to the special case that µ i is independent of E Z , and it is not related to E F being always constant. 31 Properties of Majorana modes formed in the i-th subband only depend on the value of µ i (or equivalently E i ). In the following we thus determine the effect of the electrostatics on µ i before we finally turn to Majorana bound states.

III. SCREENING EFFECTS ON CHARGE DENSITY AND ENERGY LEVELS
We begin by investigating the electrostatic effects in absence of Zeeman field and a spin-orbit strength with l SO = 233 nm, negligible for the electrostatic effects. We solve the Schrödinger-Poisson equation for a superconductor with V SC = 0 V and a superconductor with V SC = 0.2 V, and compare the solutions to two benchmarks: a nanowire without a superconducting lead, and a nanowire in which we ignore screening by charge. Specifically, we compute the influence of screening by the superconductor and by charge on the field effect on the lowest energy levels and charge densities. To evaluate the role of screening by charges in the wire, we compare the full solutions of the Poisson equation (4) to its solution with the right-hand side set to zero. Our results are summarized in Fig. 3 showing the dispersion of µ i and Fig. 4 showing the charge density for the same situations and the values of V G marked in Fig. 3.
The approximate rotational symmetry of the wire leads to almost doubly degenerate bands with opposite angular momenta when electric field is negligible-a situation realized either in absence of the superconductor [ Fig. 3 . However in most cases, presence of the superconductor leads to a large V G required to induce a finite charge density in the wire, and the degeneracy is strongly lifted.
The lever arm of the gate voltage on the energies E i , reduces from the optimal value of 1, at V G < 0 by approximately a factor of 4 due to charge screening alone [ Fig. 3(a)]. Screening by the superconductor leads to an additional comparable suppression of the lever arm, however its effect is nonlinear in V G due to the transverse wave functions being pulled closer to the gate at positive V G . Comparing panels (b) and (c) of Fig. 3 we see that screening by the superconductor does not lead to a strong suppression of screening by charge when V SC = 0: the field effect strongly reduces as soon as charge enters the wire when we take charge screening into account. This lack of interplay between the screening by superconductor and by charge can be understood by looking at the  The full self-consistent solution of the Schrödinger-Poisson equation is computationally expensive and also hard to interpret due to a high dimensionality of the space of unknown variables. We find a simpler form of the solution at a finite Zeeman field relying on the large level spacing ∼10 meV in typical nanowires. It ensures that the transverse wave functions stay approximately constant, i.e. | ψ(E Z )|ψ(0) | ≈ 1 up to magnetic fields of ∼ 7 T. In this limit we may apply perturbation theory to compute corrections to the chemical potential for varying E Z .
We write the potential distribution for a given E z in the form where φ b.c. is the potential obeying the boundary conditions set by the gate and the superconducting lead, and solves the Laplace equation The corrections φ i to this potential due to the charge contributed by the i-th mode out of the N modes below the Fermi level then obeys a Poisson equation with Dirichlet boundary conditions (zero voltage on the gates): where we write the chemical potential at a finite value of E Z as µ i (E z ) = µ i +δµ i where µ i is the chemical potential in the absence of a field. We now define a magnetic field-independent reciprocal capacitance as which solves the Poisson equation Having solved the Schrödinger-Poisson problem numer- . The correction δE i to the subband energy E i is then given in first order perturbation as Using Eqs. (12), (14) and δµ i = −δE i we then arrive at: with the elements of the reciprocal capacitance matrix P given by Solving the Eq. (15) self-consistently, we compute corrections to the initial chemical potentials µ i . The Eq. (15) has a much lower dimensionality than Eq. (5)  We start by computing the electrostatic response to changes in the magnetic field when the Fermi level is close to the band bottom for a single band (N = 1, and we write the index µ 1 ≡ µ for brevity). We study the influence of the electrostatic environment and assess whether the device is closer to a constant charge density or constant chemical potential situation (using the nomenclature of Ref. 20 explained in App. A).
The top panel of Fig. 5 shows the chemical potential response to Zeeman field. Without a superconducting contact, the electron-electron interactions in the nanowire are screened the least, and the Coulomb effects are the strongest, counteracting density changes in the wire. In agreement with this observation, we find the change in chemical potential µ comparable to the change in E Z . Hence, in this case the system is close to a constantdensity regime.
A superconducting contact close to the nanowire screens the electron-electron interaction in the wire due to image charges. The chemical potential is then less sensitive to changes in magnetic field. We find that this effect is most pronounced for a positive work function difference with the superconductor V SC = 0.2 V, when most of the electrons are pulled close to the superconducting contact. Then, the image charges are close to the electrons and strongly reduce the Coulomb interactions. In this case the system is close to a constant chemical potential regime. For V SC = 0 V screening from the superconducting contact is less effective, since electric charges are further away from the interface with the superconductor. Therefore in this case, we find a behavior intermediate between constant density and constant chemical potential.
Besides the dependence on the electrostatic surrounding, the magnetic field response of the chemical potential depends on the spin-orbit strength. Specifically, the chemical potential stays constant over a longer field range when the spin-orbit interaction is stronger. 32 The bottom panel of Fig. 5 explains this: when the spin-orbit energy E SO E Z , the lower band has a W-shape (bottom left). A magnetic-field increase initially transforms the lower band back from a W-shape to a parabolic band, as indicated by the dashed red lines. During this transition, the Fermi wavelength is almost constant. Since the electron density is proportional to the Fermi wavelength, this means that both the density and the chemical potential change very little in this regime. We thus identify the spin-orbit interaction as another phenomenon driving the system closer to the constant chemical potential regime, similar to the screening of the Coulomb interaction by the superconductor.
At large Zeeman energies E Z E SO , the spin-down band becomes parabolic (bottom right of Fig. 5). This results in the slope of µ(E Z ) becoming independent of the spin-orbit coupling strength, as seen in the top panel of  Close to the band bottom and when spin-orbit interaction is negligible, we study the asymptotic behavior of µ and n by combining the appropriate density expression Eq. C6 with the corrections in the chemical potential Eq. 15. In that case, the chemical potential becomes We associate an energy scale E P with the reciprocal capacitance P , given by and study the two limits E P E Z and E P E Z . In the strong screening limit E P E Z we find the asymptotic behavior µ ≈ −E Z , corresponding to a constantdensity regime. The opposite limit E P E Z yields µ ≈ − √ E P E Z , close to a constant chemical potential regime. We computed E P explicitly for the chemical potential variations as shown in the top panel of Fig. 5. For a nanowire without a superconducting lead, we find an energy E P ≈ 42 meV E Z , indicating a constantdensity regime. Using the classical approximation of a metallic cylinder above a metallic plate, we find an energy of the same order of magnitude. For a nanowire with an attached superconducting lead at V SC = 0 V, we get E P ≈ 7 meV ∼ E Z , intermediate between constant density and constant chemical potential. Finally, a superconducting lead at V SC = 0.2 V yields E P ≈ 0.5 meV E Z , indicating a system close to the constant chemical potential regime.
Since integrating over density-of-states measurements yields δn 0 , the inverse self-capacitance −e ψ 0 |P 0 |ψ 0 can be inferred from experimental data by fitting the density variation curves to the theoretical dependence µ(E Z ). This allows to experimentally measure the effect of the electrostatic environment, knowing the remaining Hamiltonian parameters. We compare the response to Zeeman field in the multiband case for N = 3 and N = 10 to the single band behavior in Fig. 6. We observe that presence of extra charges further reduces the sensitivity of the chemical potential to the magnetic field. We interpret the nonmonotonous behavior of the chemical potential (most pronounced for N = 10 in Fig. 6, but in principle present for all N ) as being due to a combination of the Van Hove singularities in the density of states and screening by charges. For a fixed chemical potential, the upper band, moving up in energy due to the magnetic field, loses more states than the lower band acquires, since it approaches the Van Hove singularity in its density of states. To keep the overall density fixed, the chemical potential increases. Once the density in the lower band equals the initial density, the upper band is empty and the chemical potential starts dropping again. In the limit of constant density and a single mode the magnetic field dependence of the chemical potential can be solved analytically, reproducing the non-monotonicity and kinks (see App. D).
Relating the variation in µ i to density measurements is experimentally inaccessible for N > 1, since corrections to µ i depend on the density changes of each individual mode, as expressed in Eq. (15).

A. Shape of the Majorana phase boundary
The nanowire enters the topological phase when the bulk energy gap closes at a Zeeman energy of E Z = µ 2 + ∆ 2 . The electrostatic effects affect the shape of the topological phase boundary through the dependence of µ on E Z . To find the topological phase boundary as a function of both experimentally controllable parameters V G and E Z , we perform a full self-consistent simulation at E Z = 0. We then compute corrections to the resulting chemical potential at arbitrary E Z using Eq. (15), and find topological phase boundary E Z = µ 2 + ∆ 2 by recursive bisection. Figure 7 shows the resulting phase boundary corresponding to ∆ = 0.5 meV. The phase boundary has a non-universal shape due to the interplay between electrostatics and magnetic field. In agreement with our previous conclusions, the electrostatic effects are the strongest with absent work function difference V SC = 0 V (top panel of Fig. 7) when the nanowire is intermediate between constant density and constant chemical potential. 33 Close to the band bottom, the charge screening reduces changes in density, and thus lowers the chemical potential by an amount that is similar to E Z . Hence, the lower phase boundary (at smaller V G ) has a weaker slope than the upper phase boundary (at larger V G ). Note that in the limit of constant density, the lower phase boundary would be a constant independent of E Z (see App. D).
For a work function difference V SC = 0.2 V, the system is closer to the constant chemical potential regime. In this regime, µ changes linearly with V G , yielding a hyperbolic phase boundary with symmetric upper and lower arms and its vertex at E Z = ∆. When spin-orbit interaction is strong, a transition in the lower arm of the phase boundary from constant chemical potential (hyperbolic phase boundary) to constant density (more horizontal lower arm) occurs, resulting in a 'wiggle' which is most pronounced for V SC = 0 V and l SO = 60 nm. This feature is less pronounced for V SC = 0.2 V due to the screening by the superconductor suppressing the Coulomb interactions.  Fig. 8.

B. Oscillations of Majorana coupling energy
The wave functions of the two Majorana modes at the endpoints of a finite-length nanowire have a finite overlap that results in a finite nonzero energy splitting ∆E of the lowest Hamiltonian eigenstates. [17][18][19][20][21] This splitting oscillates as a function of the effective Fermi wave vector k F,eff as cos(k F,eff L). 20 We investigate the dependency of the oscillation frequency, or the oscillation peak spacing on magnetic field and the electrostatic environment.
A peak in the Majorana splitting energy occurs when Majorana wave functions constructively interfere, or when the Fermi momentum equals qπ/L, with q the peak number and L the nanowire length. The momentum difference between two peaks is where E Z,q is the Zeeman energy corresponding to the q-th oscillation peak. In the limit of small peak spacing, yielding the peak spacing Since k F,eff = k F,eff (E Z , µ(E Z )), we substitute We obtain the values of ∂k F /∂E Z and ∂k F /∂µ from the analytic expression for k F , presented in App. C. The value dµ/dE Z results from the dependence µ(E Z ) shown in Fig. 5. Fig. 8 shows the peak spacing as a function of E Z for a nanowire of length L = 2 µm. Stronger screening reduces the peak spacing (i.e. increases the oscillation frequency) by reducing the sensitivity of the chemical potential to the magnetic field, as discussed in Sec. IV. In addition, spin-orbit strength has a strong influence on the peak spacing, since for E Z E SO the density and thus k F,eff stay constant. This results in a lower oscillation frequency and hence a larger peak spacing. Correspondingly, we find that the peak spacing may increase, decrease, or roughly stay constant as a function of the magnetic field.
Similarly to the shape of the Majorana transition boundary, Fig. 8 shows that the peak spacing does not follow a universal law, in contrast to earlier predictions. 21 In particular, our findings may explain the zero-bias oscillations measured in Ref. 11, exhibiting a roughly constant peak spacing. Fig. 9 shows Majorana energy oscillations as a function of both gate voltage and magnetic field strength for V SC = 0.2 V, with L = 1000 nm to increase the Majorana coupling. The diagonal ridges are lines of constant chemical potential. The difference in slope between the ridges of both plots indicates a difference in the equilibrium situation: closer to constant density for weak spin-orbit coupling, closer to constant chemical potential for strong spin-orbit coupling. The bending of the constant chemical potential lines in the lower panel indicates a transition from the latter mechanism to the former mechanism, due to the increase of magnetic field, as explained in Sec. IV.

VI. SUMMARY
We have studied the effects of the electrostatic environment on the field control of Majorana devices and their properties. Screening by charge and by the superconductor strongly reduce the field effect of the gates. Furthermore, screening by the superconductor localizes the charge and induces a large internal electric field. When we assume the superconductor to have a zero work function difference with the nanowire, charge localizes at the bottom of the wire, which reduces the induced superconducting gap.
Coulomb interaction causes the chemical potential to respond to an applied magnetic field, while screening by the superconductor and spin-orbit interaction suppress this effect. If a superconductor is attached, the equilibrium regime is no longer close to constant density, but either intermediate between constant density and constant chemical potential for a superconductor with zero work function difference, or close to constant chemical potential for a superconductor with a positive work function difference. An increasing spin-orbit interaction also reduces the response of the chemical potential.
Due to this transition in equilibrium regime for increasing screening and spin-orbit interaction, the shape of the Majorana phase boundary and the oscillations of Majorana splitting energy, depend on device parameters instead of following a universal law.
We have shown how to relate the measurement of density variations to the chemical potential response. Since the Majorana signatures directly depend on this response, our work offers a way to compare direct experimental observations of both signatures with theoretical predictions, and to remove the uncertainty caused by the electrostatic environment.
Our Schrödinger-Poison solver, available in the supplementary files for this manuscript, can be used to compute lever arms and capacities for different device dimensions and geometries, providing practical help for the design and control of experimental devices.
In Ref. 20 Das Sarma et al. considered Majorana oscillations as a function of magnetic field. The authors considered there two extreme electrostatic situations that they refer to as constant chemical potential and constant density.
In particular, Ref. 20 considers a one-dimensional nanowire BdG Hamiltonian as in Eq. (8), with µ 1 being denoted as µ. In this model, the subband energy E b is fixed and set to 0. The electron density is changed by adjusting µ (shown for the E Z = 0 case in Fig. 10(a).
For fixed µ in Eq. (8) electron density will change upon changing E Z . For example, if E Z µ, E so , electron density will increase monotonically as E Z is increased (see Fig. 10(b). This constant chemical potential situation is Ref. 20 also considered the opposite case of infinitely strong Coulomb interaction. In this case the electron density is fixed, and consequently µ must change as E Z changes. This constant density situation is schematically shown in Fig. 10(c).
Appendix B: Lever arms in an InAs-Al nanowire Another promising set of devices for the creation of Majorana zero modes is an epitaxial InAs-Al semiconductor-superconductor nanowire. These systems exhibit a hard superconducting gap and a high interface quality due to the epitaxial growth of the Al superconductor shell. 34 Figure 11 shows a cross section of the device. The r = 14.6 nanowire (InAs) lies on an r = 4 dielectric layer (SiO 2 ) of thickness d dielectric = 200 nm and is connected on one side to an Al superconducting shell. The device has two gates: a global back-gate with a gate potential V BG , and a side gate with a potential V SG , separated by a vacuum gap of width d gap . We model the superconductor again as a metal with a fixed potential V SC . These three potentials form the boundary conditions of the system.
We estimate the dependence of the lever arm of the side date dE/dV SG on d gap using the self-consistent Schrödinger-Poisson simulations. We set the back gate to V BG = −3.5 V, and choose the work function difference of the Al shell equal to 0.26 eV, such that one electron mode is present at a side gate voltage of V SG = −2 V, with d gap = 145 nm, as was observed in experiments. 35 We use the band gap 0.36 eV for InAs.
Our results are shown in Fig. 12, and allow to translate the gate voltages into the nanowire chemical potential. The work for the InAs-Al device shows that our numerical algorithm is easily adjusted to different device geometries, as long as the nanowire stays translationally invariant.

Appendix C: Electron density in a nanowire
Integration over the 1D density of states yields the electron density n(E, E Z , α), related to the charge density by Eq. (3). To derive the density of states, we start from the nanowire Hamiltonian, consisting of the transverse Hamiltonian of Eq. (1) and the longitudinal Hamiltonian of Eq. (2): Assuming that the wave function has the form of a plane wave ∝ e ikz in the longitudinal direction, and quantized transverse modes ψ i with corresponding energies E i in the transverse direction (where i denotes the transverse mode number), the energies of the Hamiltonian are of energy yields where α, E, E Z , and E i are in units of 2 /2m * . The relation between the density of states D(E) and k is We obtain the density n(E i , E Z , α) by integrating the density of states up to the Fermi level E F . The Zeeman field opens a gap of size 2E Z between the upper and the lower spin band. Due to the W-shape of the lower spin band, induced by the spin-orbit interaction, we distinguish three energy regions in integrating up to E F . If E F > E Z , both spin bands are occupied and the integration yields If −E Z < E F < E Z , only the lower band is occupied, and the dispersion has two crossings with the Fermi level, yielding a density n(E i , E Z , α) = 1 π k + (E F , E i , E Z , α).
For a nonzero spin-orbit strength, we have four crossings of the lower spin band with E F if E F < −E Z (see Fig. 5, bottom panel). Since only the interval k − ≤ k ≤ k + contributes to the density, integration of the density of states yields Eqs. (C5), (C6), and (C7) provide analytic expressions for the electron density. We use these equations to calculate the charge density of Eq. (3). Optimization methods find the root of the functional form of Eq. (E1), as given in Eq. (5). As opposed to other methods, the Anderson method uses the output of last M rounds as an input to the next iteration step instead of only the output of the last round. 30 The memory of the Anderson method prevents the iteration scheme from oscillations and causes a significant speedup in computation times in comparison to other methods, and in particular the simple under-relaxation method often used in nanowire simulations. 36,37 As a test system, we take a global back gate device, consisting of a hexagonal InSb nanowire on an r = 4 dielectric layer (SiO 2 ) of thickness 285 nm, without a superconducting lead. Due to the thick dielectric layer in comparison to the Majorana device, this device is more sensitive for charge oscillations (a different number of electron modes in the system between two adjacent iteration steps). This makes the device well-suited for a performance benchmark. We compare the Anderson method to three other nonlinear optimization methods: Broyden's First and Second method 38 and a method implementing a Newton-Krylov solver (BiCG-stab). 39 Fig . 13 shows the results. In this plot, we show the cumulative minimum of the error. Plateaus in the plot correspond to regions of error oscillations. The figure shows that the Anderson method generally converges quickly and is not affected by error oscillations. However, the three other methods show oscillatory behavior of the error over a large range of iterations. Both Broyden's methods perform worse than the Anderson method, but generally converge within ∼ 10 3 iterations. The Newton-Krylov method performs the worst, having a large region of oscillations up to ∼ 10 3 − 10 4 iterations. Due to its robustness against error oscillations, the Anderson method is the most suited optimization method for the Schrödinger-Poisson problem. For a much thinner dielectric layer, such as the 30 nm layer in the Majorana device, the iteration number is typically ∼ 10 1 for all four tested optimization methods.
In our approach, we choose not to use a predictor-corrector approach 40,41 that can also be used together with a more advanced nonlinear solver such as the Anderson method. 42 The advantage of the direct approach used here is its simplicity, without a significant compromise in stability and efficiency.