{'fraction_non_alphanumeric': 0.04611732790776534, 'fraction_numerical': 0.011631061376737877, 'mean_word_length': 4.0987206085753805, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Low-energy in-gap states of vortices in superconductor-semiconductor heterostructures

The recent interest in the low-energy states in vortices of semiconductor-superconductor heterostructures are mainly fueled by the prospects of using Majorana zero modes for quantum computation. The knowledge of low-lying states in the vortex core is essential as they pose a limitation on the topological computation with these states. Recently, the low-energy spectra of clean heterostructures, for superconducting-pairing profiles that vary slowly on the scale of the Fermi wavelength of the semiconductor, have been analytically calculated. In this work, we formulate an alternative method based on perturbation theory to obtain concise analytical formulas to predict the low-energy states including explicit magnetic-field and gap profiles. We provide results for both a topological insulator (with a linear spectrum) as well as for a conventional electron gas (with a quadratic spectrum). We discuss the spectra for a wide range of parameters, including both the size of the vortex and the chemical potential of the semiconductor, and thereby provide a tool to guide future experimental efforts. We compare these findings to numerical results.

The recent interest in the low-energy states in vortices of semiconductor-superconductor heterostructures are mainly fueled by the prospects of using Majorana zero modes for quantum computation. The knowledge of low-lying states in the vortex core is essential as they pose a limitation on the topological computation with these states. Recently, the low-energy spectra of clean heterostructures, for superconducting-pairing profiles that vary slowly on the scale of the Fermi wavelength of the semiconductor, have been analytically calculated. In this work, we formulate an alternative method based on perturbation theory to obtain concise analytical formulas to predict the low-energy states including explicit magnetic-field and gap profiles. We provide results for both a topological insulator (with a linear spectrum) as well as for a conventional electron gas (with a quadratic spectrum). We discuss the spectra for a wide range of parameters, including both the size of the vortex and the chemical potential of the semiconductor, and thereby provide a tool to guide future experimental efforts. We compare these findings to numerical results.

I. INTRODUCTION
The research of in-gap states of vortices started in the 1960s when approximate solutions to the Bogoliubov-de Gennes (BdG) equations for a vortex in pure type-II swave superconductors (SC-s) were found and the scaling of the states bound to the vortex core was determined as ∆ 2 SC /E F ; with ∆ SC the bulk superconducting gap and E F the Fermi energy [1,2]. Those states arise due to the phase winding of the superconducting order parameter around the normal core [3] and their scaling was verified numerically in Ref. [4]. The study of in-gap states became even more exciting when it was realized that vortices of chiral p-wave superconductors exhibit zero-energy Majorana modes [5] that braid as non-Abelian particles [6,7]. As there are no superconductors known that unambiguously show a chiral p-wave pairing, these ideas remained academic for some time.
Recently, it was realized that the combination of a conventional s-wave superconductor with the surface states of a topological insulator (TI) in a heterostructure emulates the same physics [8]. An s-wave superconductor is deposited on the surface of a three dimensional topological insulator gapping out its surface states except for a closed inner region that is left uncovered. Within this unproximitized region, a superconducting flux quantum is trapped which introduces the desired phase winding and confines a Majorana mode at zero energy. It was immediately realized that the design flexibilities of the heterostructure approach also allow to design the minigap, i.e., the gap of the Majorana mode to the first excitation [9]. In fact, typical values of ∆ 2 SC /E F in natural superconductors are smaller than 1 mK, rendering the level spectrum inside a vortex core akin to a continuum. For heterostructures, the Fermi energy E F of the superconductor is replaced by the much smaller value inside * alexander.ziesen@rwth-aachen.de the semiconductor, thus allowing to reach values for the minigap in the order 1 K [9]. As such spacings are resolvable by tunneling spectroscopy [10][11][12], this insight spurred a renewed interest in understanding the in-gap spectrum of heterostructures in details [9,13,14].
For superconducting-pairing profiles varying slowly on the scale of the Fermi wavelength of the semiconductor, Refs. [15] analytically calculated the low-energy spectrum of heterostructures. They obtained the spacing of the states to be proportional to ∆ 2 ∞ /µ, where ∆ ∞ denotes the bulk value of the superconducting gap in the proximitized semiconductor and µ is the chemical potential measured with respect to the charge neutrality point in the semiconductor.
In this work, we employ an alternative method based on conventional perturbation theory that allows to calculate the energy of the low-lying Andreev states of vortices given the pairing-potential and the magnetic-field profile as input. Our results are valid for semiconductorsuperconductor heterostructures where the semiconductor is either a TI or a conventional two-dimensional electron gas (2DEG) with quadratic spectrum. The formulas allow to determine the effects of the many design parameters. Different from a natural vortex, in a heterostructure with an etched superconductor, the size of the vortex, the effective coherence length, and the chemical potential are all individually tunable.
The paper is organized as follows. In Sec. II, we introduce the model for the SC-TI system. In Sec. III, we provide results for the minigap in terms of the different design parameters. In Sec. IV, the SC-2DEG heterostructure is introduced with the results for the minigap in Sec. V. Section VI discusses oscillatory corrections in the vortex spectra in the case of a sharp gap profile. We conclude in Sec. VII. Details of the numerical calculations performed in this paper are given in appendix A. Schematics of the topological insulatorsuperconductor heterostructure. Panel (a) is the side view, panel (b) the top view. The superconductor (gray) with pairing potential ∆SC is deposited on the surface of a 3D topological insulator (blue) with chemical potential µ. A superconducting flux quantum Φ0 is pinned to a circular region (light gray) with reduced pairing potential. Figure 1 displays a schematic of the heterostructure under consideration, see also Refs. [9,14,15] for related work. An s-wave superconductor with a bulk pairing strength ∆ SC is deposited on the surface of a 3D topological insulator with chemical potential µ and a superconducting flux quantum Φ 0 = π c/e is trapped on the surface where e > 0 denotes the elementary charge. The closeness of the two materials results in a proximity induced pairing potential ∆ TI (r) on the surface of the TI [9,14]. Given the bulk value ∆ SC in the superconductor, the exact function ∆ TI (r) is dependent on the transparency between the two materials [9]. There, it has been found that for good transparency the ∆ TI (r) is proportional to ∆ SC (r). We drop the subscript TI in the following, since we always refer to the proximitized pairing function. We write the Hamiltonian of the proximitized TI surface in the BdG formalism as
The vortex at position r = 0 introduces a phase winding with ∆(r) = ∆(r)e −iϕ , ∆(r) > 0. Note that we have assumed that the setup is rotationally symmetric for simplicity. The vortex traps a superconducting flux quantum Φ 0 = 2π ∞ 0 dr rB(r) with B(r) the z-component of the magnetic field at the radius r [16]. The vector potential A has to be in the London gauge [17] and we find A(r) = (1/r) r 0 ds sB(s)e ϕ = A(r)e ϕ . The Hamiltonian for the surface states of the topological insulator is given by with v D the Dirac velocity and σ j the Pauli matrices in spin-space.
Due to the rotational symmetry, the Hamiltonian H BdG commutes with the total-angular momentum where η j are the Pauli matrices in particle-hole space. The states are thus classified by the eigenvalues m of J z which leads to the ansatz [14,18,19] Ψ(r) = e i(m− 1 with F : [0, ∞) → C 4 . Note that the factor r −1/2 ensures the Hermiticity of the radial equation below [20]. Uniqueness of the wavefunction requires that J z − 1 2 (η z +σ z ) ∈ Z which yields m ∈ Z. We split the radial Hamiltonian H R + H A into two parts with the magnetic term H A = (π v D /Φ 0 )A(r)σ x η z and the remainder The eigenvalues E of H R + H A with |E| < ∆ ∞ form the in-gap spectrum which is particle-hole symmetric. This symmetry is implemented in the radial equation by the mapping (m, E) → (−m, −E). Note that ∆ ∞ denotes the superconducting bulk gap away from the vortex, i.e., we assume lim r→∞ ∆(r) = ∆ ∞ and lim r→0 ∆(r) = 0 .
In the following, we analyze the pairing-profiles in details, see Fig. 2. Different from the discussion of the in-gap spectrum in superconductors, for heterostructures, we keep the radius R of the vortex as a free parameter independent of the coherence length ξ = v D /∆ ∞ in the TI. The term 'natural' vortex will refer to the fact that the radius of the vortex coincides with the coherence length [21]. We use ∆ a with the sharp pairing profile as a model for an etched hole in the superconductor. The profile ∆ c on the other hand corresponds to a magnetic vortex that is trapped inside the SC with R = ξ SC the coherence length of the superconductor. We show in Fig. 3 that ∆ b approximates the latter well. As the ramp profile ∆ b allows for explicit analytical expressions, we will concentrate on ∆ a and ∆ b in the following.

A. Results and comparison to numerics
Due to the presence of the superconductor, gating capabilities of the TI surface are limited and we therefore FIG. 2. Panel (a) displays the gap profiles ∆x(r), we have considered. The profile ∆a corresponds to a hole that is etched into the superconductor. On the other hand, ∆ b and ∆c model a magnetic vortex where the pairing amplitude is suppressed due to the presence of the magnetic field and the resulting screening currents. All the profiles vanish for small r/R and approach the bulk value for r R. Panel (b) shows the weighting functions w(r), defined below Eq. (11). Note that by decreasing X = R/ξ the weight w(r) is spread out to larger values of r. focus on the regime where the chemical potential is larger than the Thouless energy, µ v D /R. We first present the results and move the derivation to the next subsection. We find that the energies of the in-gap Andreev states are given by E m = m δ, m ∈ Z, with the level spacing here and below, we defined the average with the weighting function w(r) = e − 2 ξ r 0 ∆(r )/∆∞dr , see Fig. 2(b). The weighting function is close to one in-side the vortex core and decays for larger radii. Note that the weighting function for ∆ b approximates the one for ∆ c rather well. Note that the results (valid in the limit of a small Fermi wavelength) coincide with the findings of Ref. [15] derived by a different method.
Given a paring profile ∆(r) and magnetic field strength B(r), Eq. (10) immediately yields the in-gap spectrum. As such, it constitutes an important guidance for on-going experimental efforts. Note that in the selfconsistent description of superconductors these two profiles are coupled and require the solution of, e.g., the Ginzburg-Landau equation for a specific setup [22]. The contribution of the magnetic term δ A depends on the type of superconductor and the radius R. In order to obtain a result, we need the magnetic field B(r) in the region of suppressed ∆, where the weighting function w(r) is largest. A vortex traps the flux Φ 0 in an area of roughly π(R + λ) 2 centered at the normal core with λ the London penetration depth. For type-II superconductors (ξ λ) , we assume a homogeneous magnetic field of size B(r) = Φ 0 /π(R + λ) 2 throughout the region of suppressed ∆ and find the correction We have to distinguish two cases. For λ R ξ, the magnetic correction is negligible [1]. On the other hand, for λ R, the magnitude of δ A scales as Therefore, this scaling selfconsistently satisfies the condition for the perturbation theory presented in the subsequent section. For a type-I superconductor (λ ξ R) the correction needs to be evaluated separately, since the full flux quantum penetrates the region occupied by the low-energy states. Considering a strong type-I superconductor with the homogeneous field B(r) = Φ 0 /πR 2 for r < R and zero outside, we obtain Again, the result self-consistently satisfies the perturbation theory and we observe that Eq. (12) also holds for type-I, where we have R λ. We expect this expression to work rather well for the profile ∆ a , where no screening currents are present in the relevant region r R and give a reasonable estimate for the arbitrary gap profiles. In the following, we will be focusing on the second term δ ∆ .
Evaluating the contribution of the pairing term for the profiles ∆ a and ∆ b yields which is valid for µ ∆ ∞ v D /R. Note that for a sharp step, the level spacing that determines the minigap agrees in scaling with quantization modes in a box  14) and (15). Panel (x), x ∈ {a, b}, shows the energy levels for a pairing profile ∆x(r). The gray, red, blue lines correspond to m ∈ {0, 1, 2} for a 'natural' radius R = ξ. The black line depicts the result for R = 2ξ and m = 1. Note that the numerically obtained spectra agree well with the predictions in (b). For the step potential in (a), additional oscillations appear due to the sharp features in ∆a(r). and decreases faster than for the ramp profile. Because of this, it is preferable not to etch the superconductor but to have magnetic flux line trapped inside the SC [23]. Irrespective of the gap profile, gating capabilities are needed to lower µ and ensure that the energy spacing δ is experimentally resolvable. This agrees with the findings of Refs. [9] and [15] for natural vortices.
In Fig. 4, we compare the numerical solution of H R to the approximate expressions from Eq. (10). A multidimensional shooting method is employed to calculate the exact solutions throughout this paper; for details see App. A. In the case of ∆ b , the analytic solution agrees with the numerical solution for different values of m and R. Using the step profile ∆ a , however, the approximation for E m only captures the average value of the energy, omitting a significant oscillating contribution. The latter can be attributed to interference effects of the radially oscillating wavefunctions within the extended core region with ∆ = 0. Note that breaking the rotational symmetry or introducing disorder in the core region should destroy these interference effects (cf. Sec. VI).

B. Derivation
In this section, we derive the expression Eq. (10) which was used in the last section to discuss our main findings. Different from previous semi-classical methods [1,15], we employ standard perturbation theory for this task. We split the Hamiltonian into the unperturbed part H (0) with the perturbation H 1, 1, 0). The reason for this choice is that we can find an exact zero mode for H (0) , see below. Note that the full radial Hamiltonian is given by A = H A also acts as a perturbative term.
We find that the unperturbed Hamiltonian H (0) has a solution to the energy E = 0 for each angular momentum m ∈ Z; where J m is the m-the Bessel function. Thus, the unperturbed Hamiltonian H (0) features a zero mode in each angular momentum sector. The particle-hole symmetry requires that under the perturbation the state at m = 0 stays at E = 0 while the condition for the rest of the spectrum simply demands E m = −E −m . To illustrate the existence of the exact zero mode at m = 0, we solve the full Hamiltonian (1) for a strong type-I SC in App. B. This mode at m = 0 and E = 0 is a MZM.
In order for the perturbation theory to work, we need to make sure that the perturbative terms H (1) are small compared to the first excited state of H (0) . We have numerically checked that the first excited state of H (0) is at energy ∆ ∞ for µ v D /R irrespective of the pairing profile for radii of the order of the coherence length. To obtain the approximate spectrum for H R + H A , we determine the energy in first order as Before investigating the integrals explicitly, we give an estimate of the order of magnitude for H ∆ . Due to the exponential term in the weight function w(r), the probability distribution |F m (r)| 2 has the highest weight for r < R. In contrast, the gap profile in H (1) ∆ leads to dominant contributions for r > R. Therefore, we obtain an estimate of the energy by evaluating H ∆ at r = R, which yields This reproduces the scaling from Refs. [1,15] for the case of a natural vortex with R = ξ.

The explicit form of Eq. (18) is given by
If we assume that ∆(r) and A(r) are slowly varying on the scale µ/ v D , we can replace J m (x) 2 with 1/πx by averaging over the fast oscillating component of the Bessel function, see App. C. As a result, we obtain E m ≈ m δ with the minigap δ given by Eq. (10) above.

IV. THE ELECTRON GAS HETEROSTRUCTURE
For comparison, we consider in this section the case of a conventional 2DEG such as GaAs with a quadratic spectrum that is proximity coupled to a SC. In particular, this allows to connect to earlier results that have been obtained for in-gap spectra of vortex lines in pure type-II superconductors [1]. The Hamiltonian is given as H = d 2 r Ψ † (r)H BdG (r)Ψ(r) where we choose Ψ(r) = (ψ ↑ (r), ψ † ↓ (r)) T . The Bogoliubov Hamiltonian is given by (m * is the effective mass) (21) where we have neglected the Zeeman energy. The calculation proceeds analogously to Sec. II where, here, we employ the ansatz Uniqueness of the wave function now requires m ∈ Z + 1 2 . The application of Eq. (22) leads to the real Hamiltonian for the radial motion H = H R + H A , which we split in the magnetic term and the remaining parts Note that an important difference to Sec. II is the fact that in this case m is half-integer. Because of this, all the states will be gapped and no zero modes remain.

V. ELECTRON GAS ENERGY GAP
A. Results and comparison to numerics As above, we first present our main findings. The ingap spectrum is approximately given by E m = m δ . For µ 2 /m * R 2 , the level spacing assumes the form Note that comparing (25) to (10), we see that with the correspondence v F = 2µ/m * ↔ v D , we find δ = 1 2 δ. Rewriting the first summand δ ∆ in Eq. (25) in terms of k F , it coincides with the formula derived in Refs. [1,15].
As for the TI heterostructure, we find that the minigap decreases with increasing R. Analogous to section III, approximating the magnetic field to be homogeneous yields The strength of this term again depends strongly on the type of superconductor. In the following, we focus on the correction due to the pairing potential. The comparison between the R ξ approximation to δ ∆ in Eq. (25) and the numerical solution to Eq. (24) is displayed in Fig. 5.
As for the TI, we evaluate the expression explicitly for natural vortices with R = ξ and find δ ∆x = β x ∆ 2 ∞ /µ with β a ≈ 0.12, β b ≈ 0.47, β c ≈ 0.43. Thus, the results agree with the ∆ 2 SC /E F of Ref. [1] and reproduce the factor 1 2 obtained in Ref. [15] for case b.

B. Derivation
Analogously to Sec. III B, we want to set up a perturbative treatment of H R valid in the limit µ 2 /m * R 2 .
We write As before, we have a zero mode G m (r) of H (0) in each sector m given by The spectra of H (0) are gapped by ∆ ∞ for µ 2 /m * R 2 , irrespective of the pairing profile for radii of the order of the coherence length. For slowly changing ∆(r) and A(r), we obtain E m ≈ m δ with δ given by Eq. (25) analogous to section III B (details see App. D). Notably, the correction terms proportional to σ z drop out.

VI. OSCILLATION CORRECTIONS FOR STEP PAIRING POTENTIAL
In this section, we would like to understand the oscillations of the in-gap spectrum that appear for the pairing profile ∆ a , see insets of Figs. 4 and 5. In particular, they arise due to the fact that ∆ a (r) has a sharp feature close to r = R invalidating the replacement J 2 α (x) → 1/πx. A better approximation is given by (20) and (D1). In case of Eq. (20), the incorporation of the oscillation is straight forward and we obtain ( This result matches well with the numerical solution of H R , cf. Fig. 6(a). Note that the oscillations are suppressed by a factor of v D /Rµ 1. For the 2DEG, the oscillating correction is analogously given by (µ ∆ ∞ 2 /m * R 2 ) The perturbative result is compared to the exact solution of the full problem in Eq. (24) in Fig. 6(b). Note that Note that while the oscillations decrease faster than the mean in (a), for panel (b) they both scale as µ −1/2 such that we observe many energy level crossings as well as sign changes in this case.
in this case, the oscillation amplitudes are of the same magnitude as the leading term for m = 1 2 . Because of this, the oscillations for small m lead to crossings at E = 0.
While our theory accurately predicts the oscillating behavior of the in-gap energies for a sharp gap profile, we do not expect these oscillations to be relevant for the experiments. The reasons is that they are easily washed out by either disorder or imperfect radial symmetry. The latter is most relevant for the 2DEG as many levels (of different angular momentum) cross at low energies, see Fig. 6(b).

VII. CONCLUSION
In this work, we considered in-gap states of vortices in clean heterostructures. Using Rayleigh-Schrödinger perturbation theory, we rederived the concise formula Eq. (10) of Ref. [15] to predict the energy gap above the Majorana zero mode in vortices in proximitized topological insulators. With this formula, the effects of the magnetic field profile on the energy difference of the ingap states is predicted, in addition to its pairing profile dependence. We have compared different design princi-ple given by different gap profiles. We have found that it is unfavorable to etch the superconductor if the radius is larger than the coherence length ξ. We have additionally discussed the orbital effects due to the magnetic field for type-I superconductors. We have found that the magnetic field contribution is of the same order as the pairing gap contribution. Our results can guide experimental efforts aiming at increasing the energy separation of the Majorana zero mode to finite energy excitations. We have checked the approximations employed by comparing our results with numerics. In this section, we discuss the multidimensional shooting method that we used to numerically solve the systems of radial equations (5) and (24). The idea is to match the boundary conditions at r = 0, propagate the solution in a shooting step to a final position R R and see for which energies the wave functions can satisfy the boundary condition at the final point.
For small r, we can neglect the pairing and obtain the pair of solutions [x ± = (µ ± E)r/ v D ] F a (r) of Eq. (5) (with ∆ = 0). Using a standard ODE solver, we find solutions F a and F b of the differential equation (5) (for finite ∆ and r ∈ [0, R ]) that coincide withF a andF b for r → 0, see Fig. 7. The task then is to find the energy E for which there exists a linear superposition For this it is necessary that the 2 × 2 matrix has a vanishing determinant for R → ∞. For our numerical results, we have evaluated det(M 12 ) for R = R + 20ξ as a function of E and found its roots. We have checked that the results are the same choosing M 13 , M 24 and M 34 . We cannot choose the orbitals that are related by particle hole symmetry (1 ↔ 4, 2 ↔ 3).
In this section, we discuss the zero mode of the full Hamiltonian Eq. (1) for pairing profile ∆ a and strong type-I SC for which the magnetic field B(r) = Φ 0 /πR 2 is homogeneous within the region of suppressed pairing profile and zero outside. In the chosen gauge, the respective vector field follows to We solve the Schrödinger equation for E = 0 in the inside (r < R) and outside (r > R) separately and match the solutions at r = R. In the inside, the wave function results in is the Kummer function of the first kind with (a) k = a(a + 1)(a + 2) · · · (a + k − 1) and N a normalization constant [24]. Note that the wave function (B2) is the same as for a graphene quantum dot in an applied magnetic field [25]. For the outside region, we obtain (B4) with c 1 , c 2 arbitrary constants. For the solution in Eq. (17), we see that the angular momentum quantum numbers of the orbitals at m = 0 are given by ± 1 2 ± 1 2 , where one summand stems from the Berry phase of the TI and the other is a result of the phase winding of the pairing potential. In contrast to this, the quantum numbers in Eq. (B4) follow to ± 1 2 since the full superconducting flux quantum contained in the region r < R cancels the effect of the phase winding. This is a manifestation of the Meissner effect arising due to the self-consistency equations supplementing the mean-field equations.
Matching the wave functions at r = R determines c 1 and c 2 . In general these are straightforward to obtain but rather lengthy. In the limit µ v D /R, we find the approximate expressions To visualize the exact zero mode, we plot r|Ψ| 2 , the probability to find the mode at a distance r from the vortex core, as a function of r in Fig. 8. Both, the zero mode in Eqs. (B2) and (B4) and the perturbative wave function (17), where the magnetic field is not yet taken into account, converge for µ ∆ ∞ to a constant probability for r < R, which then decays exponentially in the proximitized region. This means that irrespective of the specific FIG. 9. Distance between the exact zero energy state from Eqs. (B2) and (B4) and the perturbative zero mode from Eq. (17) as a function of the chemical potential µ. For the latter, the magnetic field is neglected. The fidelity F is defined according to Eq. (B6). In case of a natural vortex R = ξ the distance between the states decreases with µ such that for µ ∆∞ the effects of the magnetic field on the zero mode are negligible. This tendency only gets stronger upon increasing R.
magnetic field, the zero mode is equally likely found at any arbitrary distance r < R for large chemical potential. Furthermore, we observe that as µ increases, the perturbative zero mode seems to approximate the exact zero mode more closely. To verify this, we use the fidelity to measure the closeness of two wave functions Ψ 1 and Ψ 2 . In Fig. 9 we show the distance 1 − F between the state in Eqs. (B2) and (B4) and the one from (17) as a function of the chemical potential µ. For µ > ∆ ∞ both states converge and the magnetic field effects are thus negligible. This tendency is enhanced for larger R.
Subsequently, the integral in Eq. (20) is explicitly solved for the step profile ∆ a (r) = ∆ ∞ Θ(r − R). Inserting the averaged Bessel functions into Eq. (20) results in where we defined X = R/ξ and approximated X 1, meaning ∆ ∞ v D /R in the second step. Additionally, we used the exponential integral If we do not average the oscillations and insert Eq. (C1), we instead find where we approximated µ ∆ ∞ v D /R in the second step. Besides the average energy contribution, there is an oscillatory term that is suppressed by an additional prefactor v D /Rµ. Analogously, the corrections of the different pairing profiles can be calculated. Note that using the approximation Eq. (C1) does not lead to significant oscillating corrections for the two remaining profiles ∆ b and ∆ c discussed here.
Appendix D: 2DEG integral solutions for µ 2 /m * R 2 In this section, we analyze the first order correction for the spectrum in the 2DEG heterostructure and explicitly evaluate it for ∆ a . Assuming µ 2 /m * R 2 and setting A(r) = 0, the correction for Eq. (27) takes the integral expression (k F = √ 2m * µ/ ) .

(D1)
Analogous to the previous section, the contribution to the integrals for k F r 1 vanishes and we approximate the Bessel functions for the dominant contribution with large arguments x = k F r 1.
Note that the contributions of the second and third term in the perturbation (27) vanish, when we make use of Eq. (C2), but not for Eq. (C1). For the latter case, we analyze the relative strengths of the perturbation terms. Proceeding analogously to Eq. (19), with ∆(R) = ∆ ∞ , we obtain as characteristic scales where the first and third term resemble the Caroli-de-Gennes-Matricon scaling [1] for R = ξ. The slope of ∆(r) around position R determines the relative strength of the remaining contribution. For ∆ b (r) and ∆ c (r), the average slope is sub-linear, such that the contribution can be neglected. For ∆ a (r) = ∆ ∞ Θ(r − R) however, this term is expected to have a sizable effect due to the sharp feature. Subsequently, the integral in Eq. (D1) is solved for ∆ a (r). Inserting Eq. (C2) into Eq. (D1) results in E m = 1 2 E m , but we need to consider X = R/ξ being a µdependent ratio. Note that in contrast to equation (C3), m takes half-integer values. If we do not average the oscillations and insert Eq. (C1), we instead find under the assumptions µ ∆ ∞ 2 /m * R 2 . Notably, the amplitude of the oscillation matches the average value and is not suppressed.