Majorana bound states in vortex lattices on iron-based superconductors

Here we briefly review the model of surface states of a strong three-dimensional topological insulator (TI) that is used throughout this work [19]. We start by defining a simple two-dimensional square lattice model,

where $\sigma^{i}$ are Pauli matrices in spin space and the Fermi velocity $\sim\lambda$ re-appears as a hopping parameter in the tight-binding description below. With four Dirac cones in the Brillouin zone, Eq. (1) cannot represent a TI which is characterised by an odd number of surface Dirac cones. This is in accord with the Nielsen-Ninomyia theorem [28, 29], which states that it is impossible to construct a purely 2D, TR invariant lattice model with an odd number of massless Dirac fermions. The surface states of a strong TI, with an odd number of Dirac cones, are holographic in the sense that they cannot exist without the topologically non-trivial 3D bulk. A way to resolve this issue is to introduce a mass term proportional to $\sigma_{z}$ [19], parametrized as $M_{k}=m[(2-\textrm{cos}k_{x}-\textrm{cos}k_{y})-\frac{1}{4}(2-\textrm{cos}2k_{x}-\textrm{cos}2k_{y})]$. With the extra mass term, we define the Hamiltonian

where $\mu$ is the chemical potential. The mass term gaps out all the Dirac cones except the one at $\Gamma\equiv(0,0)$. We anticipate that adding such a TR symmetry breaking term will not compromise the physics of the TI-SC interface that we want to study because the presence of vortices will ultimately break TR symmetry anyway [20]. Moreover the effect of TR breaking is small around the $\Gamma$ point, where $M_{k}\propto k^{4}$.

To incorporate the proximity-induced SC in the surface states, we use the Bogoliubov-de Gennes Hamiltonian given by

where $\Delta$ is the induced SC order parameter which enters with the identity matrix in spin space. The term $h_{\textrm{0}}(\bm{k})=h_{\bm{k}}(\lambda,m,\mu)$ represents the particles and its time-reversed counterpart, $-\sigma^{y}h_{\textrm{0}}^{*}(-\bm{k})\sigma^{y}=-h_{\bm{k}}(\lambda,-m,\mu)$, represents the holes in this basis. Figures 1a-b depict the spectrum of the normal state Hamiltonian in Eq. (3) when the mass term is zero and non-zero, respectively. Figure 1c shows the effect of a finite SC pairing $\Delta$. Figure 1d depicts the bands with finite $\mu$ and $\Delta$, resulting in the formation of Dirac cones at energy $\pm\mu$ at the $\Gamma$ point. Here and in the remainder of our work we focus only on the low energy physics of the model, and choose parameters for which the non-superconducting Hamiltonian Eq. (2) faithfully represents a single Dirac fermion at the $\Gamma$ point.

Upon Fourier transforming the Hamiltonian in Eq. (2) leads to the second-quantized real-space lattice model for the TI surface state that reads [20]

Here $\alpha=x,~{}y$ labels both the hopping direction and Pauli matrices, and $\psi_{{\bm{r}}}=(c_{{\bm{r}}\uparrow},c_{{\bm{r}}\downarrow})^{T}$ denotes the two-component spinor at lattice site ${\bm{r}}$. The hopping $\lambda$, mass $m$, and chemical potential $\mu$ are as above except for an absorbed lattice constant. Upon proximity-coupling the TI to a $s$-wave SC, Cooper pairs tunnel across their interface, inducing an effective $s$-wave pairing amplitude in the surface states of the TI. The resulting TI-SC model is written in terms of a standard BdG Hamiltonian [20],

where $\Delta_{{\bm{r}}}$ is the induced SC pairing at site ${\bm{r}}$. The on-site basis now expands to include both particle and hole states, and is given by the Nambu spinor $\Psi_{\bm{r}}=(\psi_{\bm{r}},i\sigma_{y}\psi_{\bm{r}}^{\dagger})^{T}=(c_{{\bm{r}}\uparrow},c_{{\bm{r}}\downarrow},c^{{\dagger}}_{{\bm{r}}\downarrow},-c^{{\dagger}}_{{\bm{r}}\uparrow})^{T}$. In this basis the pairing matrix $\Delta({\bm{r}})$ is proportional to the identity in spin space. Thus, we arrive at the eigenvalue problem

where $H_{\rm 0}(\lambda,\mu,m)$ is the Hamiltonian of particles given by Eq. (4), and its time-reversed counterpart $-H_{\rm 0}(\lambda,\mu,-m)$ describes holes. The spin-space vectors $\bm{u}_{n}({\bm{r}})$ and $\bm{v}_{n}({\bm{r}})$ encode particle and hole eigenstates, respectively, and $E_{n}$ are the corresponding energy eigenvalues. The lattice BdG Hamiltonian can likewise be viewed as stemming from a direct Fourier transform of the momentum-space version in Eq. (3). Finally, let us recall that the particle-hole symmetry of the BdG Hamiltonian is a built-in constraint of the theory that must be obeyed. We discuss apparent violations of particle-hole symmetry stemming from inconsistent gauge choices, as well as their resolution, in Sec. 2.2 and A.

We now discuss how one can impose arbitrary vortex lattices in generic SC systems, within the BdG formalism and while using periodic boundary conditions (PBC). To this end, note that for a translation-invariant system we can define a uniform superconducting order parameter, i.e. $\Delta({\bm{r}})=\Delta_{0}$. In the presence of vortices, however, the order parameter acquires a spatially varying phase, $\Delta({\bm{r}})=\Delta_{0}e^{i\phi({\bm{r}})}$, where $\phi({\bm{r}})$ winds by $2\pi$ around each vortex. Although one expects physical observables to be perfectly periodic in the presence of periodic vortex lattice the factor $e^{i\phi({\bm{r}})}$ lacks this periodicity thus precluding the naive inclusion of PBC. As first pointed out by Anderson [30] and by Gorkov and Schrieffer [31] the situation can be remedied when the external magnetic field is included in the theory to model a realistic, thermodynamically stable vortex lattice. The latter is incorporated in Eq. (6) via the Peierls substitution, according to which all hopping terms acquire a phase $\varphi_{{\bm{r}},{\bm{r}}+\alpha}=\frac{e}{\hbar c}\int^{{\bm{r}}+\alpha}_{{\bm{r}}}\bm{A}\cdot\bm{dr}$, with the magnetic vector potential $\bm{A}$. The key insight is that the lack of translation invariance in the naive BdG Hamiltonian is a gauge artifact. Correspondingly, the solution consists of finding the correct gauge in which the Hamiltonian recovers the invariance. To achieve this we shall use a modified version of the FT singular gauge transformation [21, 22], designed to circumvent the problem of branch cuts. This allows us to impose PBC, and thus to describe large vortex lattices without boundary effects.

Let us consider a vortex lattice of $2n$ vortices that we arbitrarily divide into groups A and B of $n$ vortices each. Figure 2 shows the A and B vortex sub-lattices for systems with two or four vortices per magnetic unit cell. The total phase contributions to the SC order parameter are noted as $\phi_{A}({\bm{r}})$ and $\phi_{B}({\bm{r}})$, respectively, and obey

where ${\bm{r}}^{j}_{g}$ is the position of $j^{th}$ vortex of type $g=A,B$. With this, we can define a unitary transformation [21, 22]

where $\phi_{A}({\bm{r}})+\phi_{B}({\bm{r}})=\phi({\bm{r}})$ gives the total SC phase. Because of the delta functions in Eq. (7), the above transformation has to be regularized on the lattice [22]. Note that the A and B blocks in Eq. (8) refer to the Nambu structure in the BdG Hamiltonian, cf. Eq. (3). Hence $U$ in Eq. (8) is not a pure gauge transformation, and the effective magnetic fields that are experienced by electrons and holes are changed.

Applying the unitary Eq. (8) on the BdG Hamiltonian, $\mathcal{H}\rightarrow U^{-1}\mathcal{H}U$, removes the phase factors from the SC pairing terms in the off-diagonal blocks of Eq. (5). The transformed BdG Hamiltonian matrix then is given by

where $\mathcal{V}_{A}$ and $\mathcal{V}_{B}$ are new phase factors appearing in the hopping terms of electron and hole component. They are

with $g=A,B$, and $\Phi_{0}=hc/e$ the flux quantum. The resulting total phases in Eq. (10) are manifestly gauge invariant, making the unitary transformation Eq. (8) a singular gauge transformation. Further, the integrands in Eq. (10) can be viewed as an effective vector potential, leading to the effective magnetic field

One may interpret this as quasi-particles experiencing an effective magnetic field $\mathcal{\vec{B}}_{A}$, which comprises a uniform physical magnetic field $\vec{B}$ and delta function spikes of opposing polarity at vortices of type A. Similarly, quasi-holes experience an effective magnetic field $\mathcal{\vec{B}}_{B}$. The reason for choosing an equal number $n$ of vortices of type A and B is so that $\mathcal{\vec{B}}_{g}$ is zero when averaged over the lattice. Alternatively, there is one flux spike arising due to each vortex of type A or type B that compensates for each flux quantum of the physical magnetic field. Now the condition $\langle\mathcal{\vec{B}}_{g}\rangle=0$ allows us to incorporate PBC in the transformed BdG Hamiltonian Eq. (9), which had only allowed open boundary conditions prior to the application of the unitary transformation Eq. (8). Therefore, we can create a periodic vortex lattice with arbitrarily-chosen sub-lattices of type A and B.

As a practical matter, it would seem that we just need to calculate the phase factors in Eq. (10) in a convenient manner, followed by diagonalization of the Hamiltonian Eq. (9) to find the spectrum and eigenstates. The phase factors in Eq. (10) can be constructed in terms of the reciprocal lattice vectors $\bm{G}$ of the magnetic unit cell in Fig. 2, as described in Refs. [22, 24] and A.1. This results in

where ${\bm{r}}^{j}_{g}$ are the vortex positions and $\lambda_{s}$ is the London penetration depth. However, there are some subtleties in applying the FT transformation in systems with PBC.

First, the application of the FT transformation must ultimately yield the same physical observables independent of the choice of A-B sub-lattice degrees of freedom. For instance, the four-vortex configurations, ABAB and AABB, illustrated in Fig. 2, have different assignment of A-B vortices but must be equivalent in all physical aspects. Furthermore, the particle-hole symmetry must be obeyed for any vortex configuration. These two symmetry constraints that must be respected by any physical model are broken if one directly applies the FT recipe [21, 23, 22] under PBC using the phase factors in Eq. (12) (see A.2).

The resolution of this apparent symmetry breaking lies in understanding that the electrons and holes on the underlying physical lattice must experience the same magnetic field. This translates into the condition that flux threaded through any plaquette of the lattice can only differ by multiples of $2\pi$, as dictated by the FT transformation. However, upon imposing PBC on the FT-transformed BdG Hamiltonian, we need to ensure explicitly that the fluxes through the two additional, non-contractible loops in the torus geometry (2D lattice with PBC) obey the same principle. This leads to a modified prescription for the phase factor in Eq. (12), which thus are given as

where $\alpha=\delta x,~{}\delta y$ are the hopping directions on the lattice. The correction “wave-vector” $\bm{k}$ here follows from

where $\Delta{\bm{r}}=\sum_{j\in A}{\bm{r}}_{A}^{j}-\sum_{j\in B}{\bm{r}}_{B}^{j}$ with $j$ running over the vortices of type A and B, respectively. A detailed derivation of this result is given in A.3.

We here show how to fit the physical triangular vortex lattice into the square lattice of the underlying materials model, cf. Sec. 2.1. We then discuss how one can represent the model parameters in such a way that the mean inter-vortex distance $d_{v}$ and the disorder strength, Sec. 3.4, enter as tunable parameters. Material parameters other than the mean vortex distance (set by the applied field) are chosen fixed and correspond to sensible experimental values [13, 12, 27], see Table 1. The mean vortex distance $d_{v}$ and disorder parameters, see Sec. 3.4, will also be used for the simpler Majorana lattice model in Sec. 4.

In our simulations, we choose square lattice model sizes $L_{x}\times L_{y}$ that allow to closely fit a triangular vortex pattern with the vortex positions fixed to centers of the square lattice plaquettes. Such vortex arrangements are generated from simple geometric considerations, with two example configurations depicted in Fig. 3. We here fix $x$ and $y$ directions to have the same number of square lattice cells, $L_{x,y}=L$, and further demand that they exhibit the same pattern of vortices, related by a mirror about the $x$-$y$ axis. This naturally keeps the two lattice directions equivalent, and defines a convenient magnetic unit cell for triangular lattice simulation. Let us denote the two Bravais lattice vectors of the desired, ideal triangular lattice that is rotated like in Fig. 3, as

with $c_{15}=\cos(15^{\circ})$, $s_{15}=\sin(15^{\circ})$. We have omitted an overall scale factor here. The goal now is to find a lattice arrangement of vortices that allows us to approximate these vectors, but where the entries naturally are integers (viz., integer multiples of the lattice constant). Simply put, we want to find pairs of integers $(n,m)$ that approximate

and to define lattice vectors ${\bm{\delta}}_{1}=(n,m)$ and ${\bm{\delta}}_{2}=(m,n)$. Integer pairs that realize gradually improving rational approximations of this number are

Furthermore, the triangular lattice generated by ${\bm{\delta}}_{1,2}$ has to fit into our square lattice model. Going to larger $(n,m)$ naturally demands a larger underlying lattice, which thus becomes increasingly expensive to simulate. To fix the required lattice size $L$, consider the corners $(x,y)=(0,0)$ and $(L,0)$ of the magnetic unit cell, cf. Fig. 3. Since these are identified with each other, we have to be able to connect them by a sum of lattice vectors ${\bm{\delta}}_{1,2}$. Hence

with some integers $k_{1,2}$. For the pairs $(n,m)$ in Eq. (16), we then find coefficients

The possible lattice sizes hence are $L=8,~{}15,~{}112,~{}209,$ or multiples thereof. Once $L$ is fixed, one can read off $k_{1,2}$ above. Finally, recall that our model necessitates an even number of vortices in the magnetic unit cell. For the above pairs $(n,m)$ and quoted lattice sizes $L$, the respective number of vortices is $n_{v}=8,~{}15,~{}112,~{}224$. Hence the configuration $(n,m)=(4,1)$ cannot be used as magnetic unit cell, however one may repeat it to obtain a larger, admissible unit cell with - say - $L=30$ and $n_{v}=60$. In Fig. 3 we have shown two representative examples, namely $(n,m)=(6,2)$ with $L=16$ and $n_{v}=8$ and $(n,m)=(11,3)$ with $L=112$ and $n_{v}=112$. Using sparse-matrix methods for the Hamiltonian construction and diagonalization here allows us to consider large lattice sizes $L$, both to improve the resolution and to put more vortices that also more closely match a triangular lattice. Finally, having a sizable number of vortices ($n_{v}=112$) will be important once we introduce positional disorder in the vortex lattice, see Sec. 3.4.

After fixing the lattice setup we now discuss parameter choices and their relation to experiments, including how one can simulate different mean inter-vortex distances $d_{v}$ in fixed-size lattices by rescaling. First, the experimentally determined parameters in the FeTe_0.55Se_0.45 materials platform are the SC gap $\Delta_{0}$, the Fermi velocity $v_{F}$, and the Fermi wave-vector $k_{F}$. Their values are quoted in Table 1. In our simulations, we choose the SC gap $\Delta_{0}$ as the unit of energy, and the Fermi wave-length $\lambda_{F}=k_{F}^{-1}$ as basic length scale. From these, the chemical potential $\mu$ and Majorana coherence length $\xi$ follow via well-known relations for the model in Sec. 2.1. Additional externally set parameters are the inter-vortex distance $d_{v}$ and the STM resolution $\eta_{\rm STM}$, see Table 1. For the triangular lattice, the inter-vortex distance follows by simple geometric consideration as $d_{v}=(2\Phi_{0}/\sqrt{3}B)^{1/2}$, with the magnetic flux quantum $\Phi_{0}$ and applied magnetic field $B$. We take typical values $B=1$T$-8$T, i.e. $d_{v}=70$nm$-17$nm. Now for $n_{v}$ vortices in a lattice model of size $L\times L$, the vortex distance is given by $d_{v}=d_{v}^{\rm lattice}a_{0}$, with

the vortex distance in lattice units. Since we want to set $d_{v}$ as an external parameter while $d_{v}^{\rm lattice}$ is already fixed by the lattice setup $(L,n_{v})$, we will thus have to re-scale the lattice constant $a_{0}$ in our model. The scaled lattice constant $a_{0}$ also translates to rescaled model parameters $\lambda$ and $m$. In particular, the hopping parameter $\lambda$ is determined via Fermi velocity and lattice spacing as

where we inserted natural units, cf. Table 1. The mass $m$ then is chosen as $m\simeq 0.5\lambda-1.0\lambda$, cf. Sec. 2.1. Equivalent arguments apply to the square lattice, where for the $n_{v}=2$ unit cell in Fig. 2 one has $d_{v}^{\rm lattice}=L/\sqrt{2}$.

It may seem unconventional that the lattice spacing $a_{0}$ and model parameters $\lambda$ and $m$ are varied together with the vortex distance $d_{v}$, and hence with the applied magnetic field $B$. To this end, it helps to remember that the only crucial feature of the TI surface state model in Sec. 2.1 was the presence of an isolated Dirac cone at the $\Gamma$ point. This fact is not altered by the above parameter scaling procedure. The square crystal lattice should be viewed as a means to regularize the model at short distances and thus allow for straightforward numerical simulation. An intuitive physical picture is that our simulation keeps the number of vortices $n_{v}$ in the investigated “field of view” constant, i.e. while increasing the magnetic field we also zoom in on a smaller patch of the material. As discussed above, this approach allows us to consider a fixed number of vortices for which the triangular vortex lattice fits well into the fixed-size materials model.

In figures below, cf. Fig. 4, we show the local density of states (LDOS) integrated over a window around a vortex position ${\bm{r}}_{v}$, and the spatially resolved DOS at zero energy or at energies that correspond to peaks in the LDOS. To make this concrete, let us denote the position and energy-resolved DOS as $\rho({\bm{r}},\omega)$. For site ${\bm{r}}$ of the lattice and for each eigenstate $j$ of the Hamiltonian, we have a four-component spin-Nambu vector $\Psi_{j}({\bm{r}})=[\bm{u}_{j}({\bm{r}}),\bm{v}_{j}({\bm{r}})]^{T}$ with electron- and hole-components $\bm{u}_{j}({\bm{r}})$ and $\bm{v}_{j}({\bm{r}})$. The spatially and frequency-resolved DOS then follows as

Obviously this quantity contains too much information to show in a simple plot. Instead, the LDOS at a vortex position ${\bm{r}}_{v}$ can be calculated as

with integration window of size $2W\times 2W$. The spatially resolved DOS at energy $E$ is ${\rm DOS}_{E}({\bm{r}})=\rho({\bm{r}},\omega=E)$, where we omit the integration over an energy window. Rather the spectral broadening $\eta$ (STM resolution $\eta_{\rm STM}$) is included in the calculation of spectral densities $\rho({\bm{r}},\omega)$, where we put Lorentzian line shapes of width $\eta$. Alternatively, one could also consider a thermally broadened peak at temperature $T$, i.e. $\delta_{T}(\omega)\sim 1/\cosh^{2}(\omega/T)$. Note that the LDOS as defined here is not necessarily particle-hole symmetric. While one may start with a PH symmetric $\rho({\bm{r}},\omega)$ in Eq. (20), the quantity in Eq. (21) reflects what is measured as LDOS in experiments. Indeed, we find that the LDOS plots shown in Fig. 4 and Secs. 3.3 and 3.5 generally are PH asymmetric even for non-disordered vortex lattices, though the effects are more pronounced when disorder is included. This happens even though the total, spatially averaged DOS is necessarily PH symmetric, and indicates that pairs of particle and hole eigenstates can have spatially distinct spectral weight distributions. By picking a fixed-size spatial window to represent the “vortex LDOS”, we then implicitly choose how much spectral weight of individual eigenstates is included.

To check that the scaled model introduced in Secs. 2.1 and 3.1 correctly reproduces the vortex physics of proximitized TI surface states, we simulate the simplest square lattice with two-vortex unit cell, cf. Fig. 2(a). The Majorana vortex mode (MVM) hybridization is expected to follow the analytical prediction [32]

with Fermi momentum $k_{F}$, phase $\theta$, and coherence length $\xi$. Following Refs. [27, 32] and Table 1, we take $k_{F}^{-1}=5{\rm nm}$, $\theta=\pi/4$ (or $\theta\approx 1.4$), and $\xi=13.9{\rm nm}$. As shown in Fig. 4, our model indeed closely matches the analytical formula for any $d_{v}\gtrsim 20\,{\rm nm}$. We here set $t_{0}=4\Delta_{0}$ and $\theta=\pi/4$ to fit our numerical results. For smaller vortex separations, the hybridization of MVMs and CdGM modes becomes so strong that higher-energy features of our model become important. These are non-universal, i.e. they do not reflect an isolated TI surface Dirac cone, and hence this breakdown at very high vortex densities (magnetic fields) is expected. Going forward, we thus are confident that our model faithfully reproduces the low-energy physics even for triangular or disordered vortex lattices, so long as the typical vortex separations are sufficiently large, $d_{v}^{\rm typ}\gtrsim 20\,{\rm nm}$.

We now investigate the triangular vortex lattice, using the setup of $L=112$ and $n_{v}=112$ vortices laid out in Sec. 3.1 and Fig. 3. Vortex lattice disorder will be introduced and simulated in Secs. 3.4 and 3.5. Before repeating the comprehensive analysis of Fig. 4 for the triangular lattice, let us consider a few typical vortex separations $d_{v}\approx 49\,{\rm nm},~{}35\,{\rm nm},~{}28\,{\rm nm},~{}24\,{\rm nm}$ (magnetic fields $B=1\,{\rm T},~{}2\,{\rm T},~{}3\,{\rm T},~{}4\,{\rm T}$) and review the simulation results in more detail. Results for both the spatially-resolved DOS at various energies and the frequency-resolved LDOS at a vortex position are shown in Fig. 5. For vortex separation $d_{v}\approx 49\,{\rm nm}$ we clearly recover the expected low-energy Majorana vortex modes and higher-energy CdGM modes, both resolved in real space (DOS) and in the local vortex spectrum (LDOS). In fact, at such a large vortex separation $d_{v}$, the spectral peaks are relatively sharp and not split in energy. From Sec. 3.2 and Fig. 4, recall that the MVM hybridization and thus the low-energy peak splitting closely followed the analytical prediction in Eq. (22). For the above quoted vortex separations $d_{v}$, we would thus obtain

In the calculation of spectral densities we set the broadening $\eta=0.005\Delta_{0}$, roughly reflecting the experimental STM resolution, cf. Table 1. This explains why for $d_{v}\approx 49\,{\rm nm}$ and $35\,{\rm nm}$ the MVM hybridization and splitting is not yet resolved. However note that higher-energy modes generally show a stronger hybridization and splitting, simply because they are less strongly localized and thus overlap even for larger vortex separations. This can clearly be observed in the LDOS for the second CdGM mode (third spectral peak) for $d_{v}\approx 49\,{\rm nm}$, and for both first and second CdGM modes (second and third spectral peaks) at $d_{v}\approx 35\,{\rm nm}$, see Fig. 5.

At higher magnetic fields, for $d_{v}\approx 28\,{\rm nm}$ and $24\,{\rm nm}$, it generally becomes harder to distinguish individual CdGM modes. Due to their strong spatial overlap and hybridization, the CdGM modes here form broad bands rather than individual localized vortex modes. Also the lowest-energy MVMs become split, reflecting the above sizable MVM hybridization $\sim 0.1\Delta_{0}$, and eventually broaden into bands with a number of higher-energy but lower-intensity side peaks. This is in contrast to what we found for the simple square vortex lattice in Fig. 4, where the vortex LDOS showed a sharp pair of low-energy peaks (split by the MVM hybridization) down to small vortex separations $d_{v}=20\,{\rm nm}$. To better understand this, we now analyze the LDOS and peak splitting behavior for vortex separations $d_{v}=20{\rm nm},...,50{\rm nm}$, see Fig. 6.

The LDOS at a vortex and the MVM hybridization, for the triangular lattice with varying vortex separation $d_{v}$, are shown in Fig. 6. While the overall qualitative behavior of the vortex LDOS is similar as for the square vortex lattice, cf. Fig. 4, it differs in important details. Markedly, the MVM energy splittings do not quantitatively follow the analytical prediction [32] in Eq. (22) that was obeyed in the square lattice case. Further, as already pointed out in Fig. 5, both high-energy CdGM modes but also the low-energy Majorana vortex modes for small vortex separations $d_{v}\lesssim 30\,{\rm nm}$ develop multi-peak structures with several lower-intensity “side bands”. It thus may appear that the low-energy physics of our model at high vortex density (large field) should not be interpreted in terms of a simple triangular Majorana lattice.

We here identify two main reasons for the breakdown of the naive analytical picture that would predict MVM peak splittings with a gap given by Eq. (22). First, the Majorana hybridization formula Eq. (22) was derived for the case of two vortices at fixed distance $d$ [32]. In that case, both the tunnel coupling and thus the energy splitting of the MVMs is given by $\varepsilon_{\rm MVM}$. In a triangular lattice, however, even if one directly simulates a Majorana-only model with tunnel couplings given by Eq. (22), cf. Sec. 4 below, one does not find a single “Majorana band” but rather a band with multiple side-peaks [33, 34]. Roughly speaking, this is because it is not possible to define local pairs of MVMs that hybridize and split in a triangular lattice, because the latter is frustrated. This is in contrast with the square vortex lattice, cf. Fig. 4, which is not frustrated and where the observed peak splitting closely follows the prediction of Eq. (22). Having multiple side-peaks for energetically split MVMs thus is not a flaw of our model, but rather a feature that correctly reproduces the expected triangular Majorana lattice physics.

Second, as introduced in Sec. 3.1, the triangular vortex arrays are embedded in a microscopic model that lives on a square lattice. Hence, there will always be small deviations from an ideal Majorana-only model, cf. Sec. 4, for which one can freely choose the arrangement of vortices. Since the MVM tunnel couplings are strongly dependent on distance, this effect should not be underestimated. To be more quantitative, let us revisit the embedded vortex lattice construction in Sec. 3.1. For the $L=112$ and $n_{v}=112$ vortex case, the lattice vectors were ${\bm{\delta}}_{1}=(11,3)$ and ${\bm{\delta}}_{2}=(3,11)$ (in units $a_{0}$). The base triangle spanning the lattice in Fig. 3 thus has distinct side lengths $|{\bm{\delta}}_{1}|=|{\bm{\delta}}_{2}|\approx 11.4$, and $|{\bm{\delta}}_{1}-{\bm{\delta}}_{2}|\approx 11.3$. Entering this into the MVM hybridization formula (22), one finds couplings $t_{\rm long}$ and $t_{\rm short}$ with the ratio

for vortex separations $d_{v}\approx 49\,{\rm nm},~{}35\,{\rm nm},~{}28\,{\rm nm},~{}24\,{\rm nm}$, respectively. Even if we directly simulate the Majorana-only lattice model for the near-triangular lattices used in this section, we thus expect fairly significant deviations from the perfect triangular lattice case. While from a technical standpoint the results our model produces are correct, they do not quite reflect the triangular lattice configurations we were after. In the Majorana lattice model below, see Sec. 4, we thus choose to work based on ideal (non-embedded) triangular lattices. To validate our above assessment, we simulated the ideal triangular Majorana lattice model for a small number of modes $(n_{v}=120)$, and extracted the low-energy vortex LDOS and spectral gap. We here recovered the “fanning out” of the lowest-energy peak as in Fig. 6, while the gap more closely traces the analytical result in Eq. (22), with $t_{0}=2\Delta_{0}$ and $\theta=\pi/4$. (Note a factor two in the gap energy scale compared to the square lattice, cf. Sec. 3.2 and Fig. 4.)

We now compare our results with experimental observations [13, 14, 15, 16, 18]. Due to translational symmetry, a perfect triangular lattice shows a uniform tunneling spectrum (LDOS) for all vortices, which does not happen in experiment. Setting aside minor deviations from the perfect triangular lattice, our results also show that the splitting of low-energy peaks in the LDOS should exhibit periodic oscillations, cf. Fig. 6. While the general trend of obtaining fewer vortices with zero-bias peaks for higher fields (smaller vortex separations) is consistent with a growing Majorana hybridization, oscillations were never seen [13, 14, 15, 16, 18]. Lastly, we observed a fanning out of the low-energy spectral peaks which gives a cluster of side-peaks with lower spectral weight toward higher energies. This might be a first indicator for the development of “pyramid-shape” distributions in the LDOS peak statistics [14, 27], which we will discuss in detail for the disordered vortex lattice.

As seen here, hybridization effects alone cannot explain all experimental observations [14, 27]. In the simulations up till now, vortices necessarily are identical due to translational invariance. Instead in actual experiments there are prominent variations between the individual vortex spectra [14], implying that the translational symmetry of the vortex lattice is broken and disorder effects need to be taken into account. These may include positional disorder (vortices shifted from regular lattice sites), as well as fluctuations of various parameters such as the SC pairing amplitude and chemical potential. We note that Refs. [14, 27] have explored various possible sources of disorder, and concluded that most likely positional disorder is the culprit behind the strong vortex-to-vortex variations. We thus focus on the latter mechanism, and compare our results to their conclusions where appropriate.

For the simulation of disordered vortex lattices, a key ingredient is the generation of “good” vortex lattice configurations that qualitatively agree with those observed in experiment [14, 27]. Once generated, we can easily compare and qualify them by inspecting the reciprocal space structure function $S(\bm{k})=N_{v}^{-1}\sum_{j,k}e^{-i\bm{k}\cdot({\bm{r}}_{j}-{\bm{r}}_{k})}$ or by computing pair-correlation functions. Let us denote the vortex positions in the disordered lattice as

with $\bm{R}_{j}$ the near-triangular Bravais lattice sites, see the discussion in Sec. 3.1 and Fig. 3, and $\bm{w}_{j}$ the deviations.
Inspecting the experimental data in Ref. [14], cf. Ref. [27], the vortices appear to be disordered globally while locally repellent and in a roughly triangular arrangement. This is in agreement with the basic understanding of vortex physics in superconductors. As an input to the initialization of $\bm{w}_{j}$, we hence consider two parameters: the “displacement variance” $\sigma_{d}$, and “disorder correlation length” $\lambda_{d}$. Here $\sigma_{d}$ signifies the typical allowed displacement of each vortex from its triangular lattice position, independent of its neighbors, while $\lambda_{d}$ describes the typical length scale on which vortices tend to be displaced in a similar way - thus locally retaining a triangular lattice arrangement.

An operational means by which we initialize the displacements $\bm{w}_{j}$, with the above prescription for a lattice of $n_{v}$ vortices, is as follows. Consider the set of all $\bm{w}_{j}({\bm{r}})$ as $n_{v}$ separate Gaussian random variables with a fixed covariance matrix $K$. The covariance matrix can be specified by the expectation values of products of pairs of different $\bm{w}_{j}$, and for cross-correlated disorder (as we want to consider here) will have finite off-diagonal elements. The displacements $\bm{w}_{j}$ then are chosen as Gaussian random variables with zero mean, $\langle\bm{w_{j}}\rangle=0$, and (co)variance

For simplicity we here take independent displacements in $x$ and $y$ directions with identical covariance matrix, i.e. $K^{xx}=K^{yy}=K$ and $K^{xy}_{jk}=\langle w_{j}^{x}w_{k}^{y}\rangle=0$. Note that the separations $R_{jk}=|\bm{R}_{j}-\bm{R}_{k}|$ entering Eq. (24) refer to the Bravais lattice sites $\bm{R}_{j}$ in Eq. (23) and Fig. 3. To generate random samples of displacements $\bm{w}_{j}$, one thus has to specify the covariance matrix $K$ for the desired vortex lattice (with fixed $n_{v}$ and $L$) and for the pair of parameters $\sigma_{d}$ and $\lambda_{d}$. Independent randomly displaced vortices are obtained for $\lambda_{d}\to 0$, while values larger than the typical inter-vortex separation ($\lambda_{d}\gg d$) lead to “quasi-ordered” lattices.

From the above prescription we generated a few sample lattices with different $\sigma_{d}$ and $\lambda_{d}$, and inspected the corresponding structure functions as done in Ref. [14]. Since overall scales do not matter at this moment, we set the vortex lattice spacing $d_{v}=1$. We found that parameter ranges $\sigma_{d}\in[0,0.25]$ and $\lambda_{d}\in[0,3]$ are reasonable for reproducing experiment-like vortex configurations [14, 27]. Here one may increase $\lambda_{d}$ as $\sigma_{d}$ grows larger, such that a local approximate triangular lattice is retained. In our disordered vortex lattice simulations below, we fix $\lambda_{d}=2$ and $\sigma_{d}\in[0,0.25]$. The corresponding real- and reciprocal-space lattices are shown in Fig. 7.

We now have to address one additional issue, which already affected the setup of vortex lattices in Sec. 3.1. Below we want to simulate disordered vortex lattices in the FeTe_0.55Se_0.45 platform, extending the results of Sec. 3.3. As with the embedding of regular vortex lattices, cf. Sec. 3.1, we here have to amend the disordered vortex configurations such that they fit into the underlying square lattice of the materials model. Let us start with the non-disordered but slightly distorted triangular lattice in Fig. 3, where we fixed vortices to sit in the centers of the materials lattice plaquettes. Now upon adding disorder, the vortex positions are displaced from the lattice plaquette centers, or even shifted to neighboring plaquettes. We then need to make sure that vortices keep a “safe” minimum distance from materials lattice sites, since otherwise we have to include more and more reciprocal lattice vectors in the calculation of FT phase factors. This is because the cutoff for sums over reciprocal lattice vectors, see Eq. (12) and A, determines how finely we are able to resolve modulations of the FT phase factors (and hence the SC phase) in real space. If vortices are very close to materials lattice sites, these modulations in FT or SC phases become very strong. Therefore, if one includes too few reciprocal lattice vectors in the FT phase calculation, an erroneous particle-hole symmetry breaking of the obtained spectrum occurs. Performing tests similar to the ones in A, we thus determine the maximum allowed deviations of vortex positions from plaquette centers for a given fixed and numerically feasible cutoff of reciprocal lattice summations.

This leads us to restrict the disordered vortex lattice positions to the vicinity of plaquette centers as follows: for each vortex position ${\bm{r}}_{j}$ in Eq. (23), we determine which plaquette of the materials lattice it is in. We may call this the “plaquette location” $\bm{R}_{j}^{\prime}$ of the vortex. Note that the latter can be distinct from the original (non-disordered) vortex location $\bm{R}_{j}$ if the displacement $\bm{w}_{j}$ was sufficiently large. We then calculate the embedded vortex position

with $0\leq\alpha\leq 1$. For $\alpha=1$, the embedded vortex position is simply the bare disordered position ${\bm{r}}_{j}$, while for $\alpha=0$ we force the disordered lattice positions to also be centered at plaquettes, ${\bm{r}}_{j}^{\prime}=\bm{R}_{j}^{\prime}$. In the latter case, positional disorder hence only comes in multiples of the materials lattice constant. We find that numerically feasible reciprocal lattice summation cutoffs give us sufficient accuracy in the FT phase calculation to set $\alpha=0.6$, while keeping erroneous PH symmetry breaking effects small enough to not adversely affect our results or conclusions. All vortex positions and lattices in the following section hence will be embedded in the sense that vortex positions are at most $\alpha\cdot 0.5=0.3$ linear distance away from plaquette centers, i.e. each vortex keeps a linear distance of at least $0.2$ from materials lattice sites and edges. In fact the vortex positions and reciprocal lattices in Fig. 7 were already generated by this prescription, and reflect the embedded vortex lattices used henceforth. Effects both due to the use of periodic boundary conditions (PBC) in a finite lattice and due to the vortex embedding are visible in Fig. 7. Besides the overall hexagonal pattern in reciprocal space, one can see faint horizontal and vertical ripples which are due to the square lattice.

Finally we note that in Section 4 we will introduce a Majorana-only model that effectively simulates the lowest-energy modes of the vortex lattice. That model does not explicitly simulate the underlying material, and thus none of the above embedding issues apply, discussed here and for the non-disordered lattice case in Sec. 3.1. For that reason the Majorana-only model is much easier to set up and analyze, both conceptually and numerically.

We here simulate near-triangular vortex lattices with positional disorder parametrized by the disorder variance $\sigma_{d}$ and correlation length $\lambda_{d}$ (quoted in units $d_{v}$), cf. Sec. 3.4. As for the regular triangular lattice in Sec. 3.3 and Fig. 5, we first investigate the DOS and LDOS for a few typical vortex separations $d_{v}\approx 49\,{\rm nm},~{}35\,{\rm nm},~{}28\,{\rm nm},~{}24\,{\rm nm}$. For illustration purposes, we here use strong disorder with $\sigma_{d}=0.2$ and $\lambda_{d}=2$, cf. Fig. 7 or the spatial DOS plots in Fig. 8 from which one can easily identify the vortex positions. Also note that though PH symmetry on the level of individual vortices generally is broken, below we focus on the positive-frequency part of the LDOS and spatially resolved DOS. The behavior at negative frequencies is qualitatively similar.

The main observations of these simulations are discussed in Fig. 8. Briefly, the disorder effects in the spatial DOS and vortex LDOS become more dramatic as the vortex separation $d_{v}$ is lowered. This is expected, since the hybridization energy scales increase and fluctuations in the latter can be resolved once they surpass the spectral broadening $\eta$. Also at low vortex separations the disorder conspires with the lattice frustration effects that lead to side-peaks for Majorana and CdGM modes, and produces a more complex multi-peak LDOS. Finally, the low-energy peaks in LDOS and the zero-energy spatial DOS are most sensitive to disorder when the corresponding clean system is near a node of the Majorana hybridization oscillation. This is the case for $d_{v}\approx 35\,{\rm nm}$ in Fig. 8, where some of the vortices show zero-bias peaks in the LDOS, while others do not. We attribute this to the combination of strong fluctuations of coupling strengths under weak variation of the inter-vortex distance, together with the overall sizable mode hybridizations at small distances (overcoming the spectral broadening $\eta$). As discussed in detail in Sec. 4, the fact that the sign of the effective Majorana hybridizations fluctuates between distinct modes then is the most crucial ingredient to obtain a strong LDOS variation [33]. Qualitatively, the results in Fig. 8 look more similar to the experimental results of Machida et al. [14] as compared to the non-disordered case, Fig. 5. This motivates us to further investigate the effects of positional disorder.

To make more qualified observations, we next consider disorder-averaged spectral functions. We here focus on the lowest-energy (most localized) modes and solve for “only” the $5\cdot n_{v}$ eigenstates closest to zero energy. Hence the overall spectral weight at high energies will be suppressed or altogether absent. In contrast, in Fig. 8 and to properly resolve higher-energy modes, we solved for significantly more eigenstates of the Hamiltonian which is correspondingly more expensive. The results for three vortex separations $d_{v}\approx 49{\rm nm},~{}35{\rm nm},~{}28{\rm nm}$ and various disorder strengths are shown in Fig. 9. As for a single disorder realization in Fig. 8, for a large vortex spacing $d_{v}\approx 49{\rm nm}$ the disorder effect is hardly resolved, since the fluctuations of the vortex mode hybridizations are on the scale of the spectral broadening $\eta$. When the average vortex spacing is close to a node in the MVM hybridization, $d_{v}\approx 35{\rm nm}$, the initially split low-energy peaks merge and broaden as the disorder strength is increased, but a zero-energy “hump” remains. Hence, on average we would expect to find a peak at or close to zero energy more likely than not when scanning over an ensemble of vortices. Higher-energy CdGM modes that initially are split likewise merge into a smoothed out hump of spectral weight. Finally, at a small average vortex spacing $d_{v}\approx 28{\rm nm}$ where the low-energy modes are strongly split, disorder becomes ineffective in closing the gap. While some spectral weight is shifted towards zero energy, a distinct hump is left where the zero-disorder clusters of spectral peaks were located, cf. Fig. 5. We again note that very similar behavior for the low-energy (Majorana) modes of the vortex lattice will be recovered in the simpler Majorana-only model, see Sec. 4 below.

Motivated by the experimental analysis of Refs. [14, 27], we finally consider the spectral peak statistics in the LDOS of vortex modes. I.e. we here inquire about an ensemble of vortices ($n_{v}=112$) that are in the “field of view” of our simulations (viz., the STM in experiment), and deduce the presence or absence of zero- and low-energy modes based on the locations of peaks in the vortex’ LDOS. In Fig. 10 we consider a single vortex lattice and disorder realization, to illustrate how the peak statistics are obtained. First, we calculate the vortex LDOS for all $n_{v}$ vortices in the lattice, cf. Figs. 8 and 9. The peaks in the vortex LDOS are identified simply by finding its local maxima, and peak heights and locations are indicated for a subset of vortex spectra in the inset of Fig. 10. We then histogram the peak positions $\omega_{n}$, with the result shown in the main panel of Fig. 10. Overall, after considering several disorder realizations, we find that the peak position distributions tend to be “pyramid-shaped”, with the lowest spectral peaks clustering towards zero energy. One quantity that often is quoted in STM studies of the vortex LDOS is the zero-bias peak rate (ZBPR), i.e., the fraction of vortices that shows a zero-energy peak. For the parameters considered here, we find ${\rm ZBPR}\approx 0.61$.

Beyond a single disorder realization, we consider the disorder-averaged LDOS peak statistics shown in Fig. 11. Here we find clear evidence for the pyramid-shaped distribution of LDOS peak positions. As mentioned in the discussion of our non-disordered vortex lattice results, cf. Sec. 3.3 and Fig. 5, we can partially attribute them to the fanning out of the Majorana and CdGM bands in the triangular lattice model. Ultimately what is behind this is not the specific form of model or disorder, but rather the generic fact that (Majorana) fermions on a triangular lattice are frustrated. As expected based on the behavior of the averaged LDOS, cf. Fig. 9, the ZBPR increases significantly only when the non-disordered system already is close to a gapless state (for vortex spacing close to a node of the MVM hybridization oscilllation). If the system is initially strongly gapped, it takes rather strong disorder to generate any vortices that show zero-bias peaks in their LDOS. For $d_{v}\approx 28{\rm nm}$, this means a disorder variance as large as $\sigma_{d}=0.2$, see Fig. 11.
We thus conclude this section by emphasizing that a key feature of vortex modes in experiments [14, 27], i.e. pyramid-shaped low-energy peak distributions, are recovered in our model. In agreement with the analysis of Refs. [14, 27], it was sufficient to include simple positional disorder of vortices as opposed to chemical or other disorder effects.

Here we show how to set up the basic Majorana lattice model, to which one can introduce a tunable degree of positional disorder as discussed in Sec. 3.4. A simple way to implement a triangular lattice is to choose a two-site unit cell of “$a$” and “$b$” Majorana vortex modes (MVMs). Their relative positions in the unit cell are given as

with $dx=1$ and $dy=\sqrt{3}$. The lattice translation vectors $\bm{R}_{x}=(dx,0),\bm{R}_{y}=(0,dy)$ span a rectangular lattice of this two-site unit cell. MVM positions are thus labelled by two integers, $n_{x}$ and $n_{y}$, and their “flavor” $a$ or $b$:

Aside from the form and choice of disorder, cf. Sec. 3.4, there is only one ingredient to the model. The Majorana hybridization $t_{jk}$ between MVMs localized at ${\bm{r}}_{j}$ and ${\bm{r}}_{k}$ follows from the well-established analytical form [32, 38, 24, 27]

where $t_{\rm mat}(r_{jk})$ denotes a material-dependent Majorana hybridization, and the $\sin(\omega_{jk})$ factor captures geometric vortex lattice phases. Note that $t_{\rm mat}$ depends only on the relative distance $r_{jk}=|{\bm{r}}_{j}-{\bm{r}}_{k}|$ between the MVMs. We first specify the materials-dependent part, cf. Eq. (22),

with Fermi momentum $k_{F}$, a phase $\theta$, and coherence length $\xi$. The bare tunneling amplitude $t_{0}$ sets the overall energy scale of the model; we henceforth put $t_{0}=2\Delta_{0}$ ($=3.8$ meV), cf. the discussion in Sec. 3.3. This will allow us to make quantitative comparisons with results in Sec. 3. According to Table 1, we further take $k_{F}^{-1}=5{\rm nm}$, $\theta=\pi/4$ (or $\theta\approx 1.4$), and $\xi=13.9{\rm nm}$. Below we measure all distances in nanometer (${\rm nm}$), and scale the lattices with the characteristic inter-vortex distance $d$. Note that $d=15-80{\rm nm}$ in experiments [27, 14].

Now consider the geometric vortex lattice phases, entering via $\omega_{jk}$ in Eq. (28). The simplest way to evaluate $\omega_{jk}$ is by following the Grosfeld-Stern (GS) rule [35], which strictly only holds for regular lattices like the triangular or square lattice cases considered in Refs. [24, 38]. Here we assume that the GS rule also applies in moderately disordered lattices, and discuss possible deviations. Taking the results of Refs. [24, 38] and including next-neighbor hopping, we obtain the phase assignments shown in Fig. 12.

For disordered lattices and without including phase factors $\omega_{jk}$ beyond the GS rule, we have to argue whether we can trust our numerical results at least qualitatively. Inspecting Eq. (28), the main effect of phases $\omega_{jk}$ different than those in the regular lattice will be to further modulate and randomize couplings $t_{jk}$. Also local fluctuations of physical parameters $k_{F}$, $\theta$, and $\xi$ will tend to randomize the couplings $t_{jk}$. Introducing just one source of randomness and one way in which it enters, namely the positions of MVMs and their separations entering $t_{jk}$, might underestimate the effect of disorder. Conversely, to capture the experimentally observed effects, we might have to increase the disorder strength to larger values than naively expected. Still this allows for a significant computational simplification since calculating the tunneling amplitudes only requires knowing the distances between MVMs, Eq. (29), while to calculate $\omega_{jk}$ exactly one has to consider the full MVM lattice arrangement. We have also repeated the analysis of the disordered triangular lattice as in Sec. 3.4 and Fig. 7, but now using the base triangular lattice defined by Eq. (27). While we here do not provide an additional figure, note that since for the Majorana-only model we are able to go to significantly larger vortex lattices, the effects of employing PBC on a torus are much less visible. Further the embedding problem is completely avoided in the Majorana-only model, since there is no underlying materials lattice.

Let us also note that in earlier work, Kraus and Stern [33] and Laumann et al. [34] considered the case of disordered tunneling amplitudes and phases in Majorana lattices; these were simply picked as random variables. Naturally it would be easy for us to do the same, and pick the phase values $\omega_{jk}$ according to some random distribution. However we do not expect to gain much additional insight from this. Rather, here we take into account fluctuating MVM separations as motivated from and observed in the recent experiments [14, 27]. Still it will be instructive to compare our results with those obtained in Refs. [33, 34].

With the setup of the Majorana lattice in Secs. 4.1 and 3.4, Figs. 12 and 7, we now calculate GS phase factors and tunneling amplitudes according to Eq. (28). It is useful to rewrite the Majorana tight-binding Hamiltonian (26) as

with Majorana vectors $\hat{\gamma}^{T}=(\gamma_{1},...,\gamma_{2N})$, and $\hat{h}$ a real anti-symmetric matrix with entries $t_{jk}$ in Eq. (28). After a Schur decomposition of the matrix $\hat{h}$ into real, anti-symmetric $2\times 2$ blocks on its diagonal [8], one obtains

where $\tilde{h}=O\hat{h}O^{T}$ and $\tilde{\gamma}=O\hat{\gamma}$ with orthogonal matrix $O$. The pairs $\tilde{\gamma}_{2m},~{}\tilde{\gamma}_{2m+1}$ then define $N$ complex fermions $c_{m}=\tilde{\gamma}_{2m}+i\tilde{\gamma}_{2m+1}$, and $H=\sum_{m}\varepsilon_{m}(c_{m}^{\dagger}c_{m}-\frac{1}{2})$. One can now bin and plot the energy eigenvalues $\varepsilon_{m}$, averaged over disorder realizations or with varying $\sigma_{d}$ and $\lambda_{d}$, to compare with the earlier results of Refs. [33, 34].

For us more interestingly, we can obtain the Greens and spectral functions for the original Majoranas $\gamma_{j}$ by using the above transformation $O$. The local spectral density of $\gamma_{j}$ is what is measured by an STM tunneling experiment [14, 15, 16] into the MVM and vortex at position ${\bm{r}}_{j}$. Let us consider the Lehmann representation of the retarded Majorana Greens function (GF), at zero temperature,

where $|n\rangle$ is the many-body (MB) eigenstate at energy $E_{n}$, and $|0\rangle$ is the ground state. To relate to eigenvalues $\varepsilon_{m}$, we write $\gamma_{j}$ via $\tilde{\gamma_{j}}$ in terms of complex fermions $c_{m}$,

Here $\psi_{mj}=\frac{1}{2}(O_{2m,j}-iO_{2m+1,j})$ is derived from the above orthogonal transformation $O$. Since we are dealing with a single-particle problem, the MB states are product states $|n\rangle=|n_{0}\cdots n_{N-1}\rangle$. For the MB ground state $|0\rangle$, each fermion is filled or empty depending on the sign of the corresponding energy $\varepsilon_{m}$. Applying $\gamma_{j}$ to $|0\rangle$, either $c_{m}$ or $c_{m}^{\dagger}$ thus annihilate the state. For the $m$th term in the sum, matrix elements in Eq. (32) hence evaluate as

If we remove or add an electron from/to the ground state (first/second term), the energy is lowered/raised by $\varepsilon_{m}$, giving $E_{0}-E_{n}=\mp\varepsilon_{m}$. Since we consider a particle-hole symmetric form for the GF, adding $E_{0}\leftrightarrow E_{n}$ in Eq. (32), energy denominators with either sign are present in the final result. Hence the retarded Majorana GF for modes $\gamma_{j}$ reads

We thus obtain a symmetrized sum of retarded GFs of the systems’ eigenmodes $c_{m}$, weighted by their amplitudes $|\psi_{mj}|^{2}$ at the MVM lattice site ${\bm{r}}_{j}$. The associated local spectral function then follows as $\rho_{j}(\omega)=-\frac{1}{\pi}{\rm Im}G_{j}^{R}(\omega)$. Note that in constrast to the LDOS and spatial DOS in the TI-SC model, cf. Sec. 3.2, the Majorana LDOS as defined here is necessarily PH symmetric.

We can now calculate the local spectra $\rho_{j}(\omega)$ for large lattice sizes ($N=L_{x}L_{y}$ unit cells, $2N$ MVMs, PBC), varying the mean vortex distance $d$ (magnetic field $B$) and disorder parameters $\sigma_{d}$ and $\lambda_{d}$. The broadening $\eta$ here mimics the finite spectral resolution of the STM measurement in experiments, as in Sec. 3. Depending on parameters, we then check how many of the MVMs have zero-bias peaks (zero-bias peak rate, ZBPR), what is the distribution of the lowest and higher-energy peak positions across all vortices, and so on. Having access to large lattices with thousands of MVMs allows us to obtain better statistics than is possible with the more detailed models of Sec. 3 and Ref. [27]. This approach thus should enable us to reproduce some of the experimental observations [14, 15, 16, 18], and provide valuable information on whether the low-energy physics in iron-based SCs can be understood in terms of Majorana vortex modes.

Here we discuss some results from the MVM lattice simulations outlined in Secs. 4A-B. The single-particle density of states (DOS) of the model is shown in Fig. 13. We obtain it by binning the energies $\varepsilon_{m}$ in Eq. (31), with an average over disorder realizations for fixed parameters $\sigma_{d}$ and $\lambda_{d}$, cf. Sec. 3.4, and fixed spacing $d$.

Comparing to results of Kraus and Stern [33], we find that the DOS in Fig. 13 are qualitatively similar to their random-hopping models with nearest-neighbor terms. While we have a more regular structure of couplings, derived from the underlying near-triangular lattice, this difference does not seem to affect the density of states too much. However we do expect the vortex spacing $d$ to play an important role, especially for weak disorder when the cosine-modulation of couplings in Eq. (29) is prominently visible. To this end, note that parameters $d$ and $\sigma_{d}$ translate into typical vortex separations $r\in(d-\sigma_{d},d+\sigma_{d})$. Entering these into Eqs. (28)-(29), we notice that if the average vortex spacing $d$ is tuned close to a node of the oscillating MVM hybridization, couplings $t_{jk}$ will be random both in sign and amplitude. Instead if $d$ is close to a local maximum (or minimum), introducing a small variance $\sigma_{d}$ only scrambles the amplitude but not the sign of $t_{jk}$’s. As shown in Ref. [33], the two cases are quite different: random-sign models are gapless and show a peak of the DOS at zero energy, while with random amplitude (but fixed sign) the model stays gapped. We hence expect a gapped system, without zero-bias peaks for the local DOS at vortex cores (see below), whenever $d$ is far from a node of the MVM hybridization and disorder is weak (small $\sigma_{d}$ and/or large $\lambda_{d}$). When $d$ is close to a node, even weak disorder will randomize the signs of couplings and lead to a gapless system with zero-modes in many vortices. Finally, for strong disorder (large $\sigma_{d}$), the precise value of the spacing $d$ does not matter anymore and typical vortex separations $r\in(d-\sigma_{d},d+\sigma_{d})$ lead to random-sign couplings and hence a gapless system. Indeed the results shown in Fig. 13 reflect these trends. Comparing to results in the TI-SC model of Sec. 3.5, cf. Fig. 9, we find that its low-energy modes show very similar behavior as the Majorana-only model in Fig. 13.

Next we consider the LDOS and spectral peak statistics of MVMs, cf. Sec. 3.5 and Refs. [14, 27]. As discussed in Sec. 4.2 and shown in Fig. 14, calculating the LDOS allows to extract spectral peak frequency histograms and the zero-bias peak rate (ZBPR). When the system becomes gapless and consecutively develops a zero-bias peak in the total DOS of Fig. 13, it also shows a non-zero ZBPR for the LDOS of MVMs. Moreover, we find that the lowest and second-lowest energy peak positions $\omega_{1,2}$ generically are distributed in a “pyramid-shape” way in Fig. 14, in accordance with observations in experiment [14] and the more complex models of Sec. 3 and Ref. [27].

Beyond a single disorder realization in Fig. 14, we again take disorder averages to analyze trends in the ZBPR and peak statistics upon varying the disorder parameters $\sigma_{d}$, $\lambda_{d}$ and spacing $d$. This data is shown in Fig. 15.

The MVM peak statistics observed in Fig. 15 closely agree with our intuition gained from Fig. 13, and also with the more complex model in Sec. 3.5 and Fig. 11. Generally, disorder broadens spectral features and smears out the distributions of spectral peak positions in Fig. 15. The overall trends depend on whether the system initially, in the absence of disorder, is strongly gapped or close to a gapless state (spacing $d$ far from or close to a node of the oscillating MVM hybridization, Eq. (29)). If the initial system is near-gapless, disorder quickly closes the gap and leads to an accumulation of spectral peaks at zero energy. Ramping up $\sigma_{d}$ then spreads out the spectral weight and actually decreases the ZBPR, which is different from the full TI-SC model in Fig. 11. This may be due to the absence of higher-energy CdGM modes in the Majorana-only model, which otherwise would provide a source from which spectral weight also is transferred to lower energies. For a strongly gapped system, the threshold to generate a finite ZBPR is larger, but the ZBPR increases even for strong disorder. This is more similar to the behavior in the full TI-SC model. Finally at very strong disorder, say $\sigma_{d}>0.2d$, the cosine-modulation of the MVM hybridization Eq. (29) becomes inessential. The main trend that determines the ZBPR then is the decay $\langle t_{jk}\rangle_{\sigma_{d},\lambda_{d}}\sim e^{-d/\xi}/\sqrt{d}$ of the hybridizations $t_{jk}$ with increasing vortex spacing $d$.

Last, we analyze how the ZBPR varies as function of the spacing $d$ and disorder strength $\sigma_{d}$. The extraction of the ZBPR here proceeds analogous to Figs. 14 and 15, and results are shown in Fig. 16. Naturally the ZBPR also depends on the disorder correlation length $\lambda_{d}$ and the peak broadening $\eta$ (“STM resolution”), and we here pick values that are convenient for illustration purposes.

The ZBPR in Fig. 16 for small to intermediate $\sigma_{d}$ shows clear signatures of the hybridization oscillations, Eq. (29). Whenever the MVM hybridization is strong the ZBPR is suppressed, and vice versa. For stronger disorder $\sigma_{d}\gtrsim 0.15d$ and large spacing $d\simeq 30-60{\rm nm}$ that suppresses the overall energy scales in the system, the ZBPR quickly averages to around $0.4-0.6$. The trend of increasing ZBPR with increasing spacing $d$, at strong disorder, is clearly observed in Fig. 16.

Finally, let us relate back to experimental findings. Machida et al. [14] observe an increase in vortex lattice disorder as the magnetic field is ramped up. In terms of our model, we hence should consider an increasing disorder variance $\sigma_{d}$ with fixed or decreasing disorder correlation length $\lambda_{d}$ as the vortex spacing $d$ is reduced. In Fig. 16, this means taking some line cut that starts at large $d$ and small $\sigma_{d}$, and ends at small $d$ and large $\sigma_{d}$. Likewise one could investigate the effect of a changing broadening $\eta$, reflecting the STM resolution in experiment [14, 13, 18, 15, 16, 12]. In this section we set $\eta=0.002t_{0}$, which according to our discussion in Figs. 4 and 6 translates to $\eta=0.004\Delta_{0}\approx 7.2\,\mu eV$. This is close to the experimentally quoted $\eta_{\rm STM}\approx 20\mu eV$. Also note that Ref. [14] employs a more sophisticated peak-fitting procedure than simply checking for maxima in the local spectral density, which effectively enhances their spectral peak resolution to a larger value. Last and more importantly, our simulations are implemented in the zero-temperature limit; here one could also include finite temperature, in form of an additional Gaussian broadening of the spectral peaks.