Manifestation of chirality in the vortex lattice in a twodimensional topological superconductor
Abstract
We study the vortex lattice in a twodimensional wave topological superconductor with Rashba spinorbit coupling and Zeeman field by solving the Bogoliubovde Gennes equations selfconsistently for the superconducting order parameter. We find that when spinorbit coupling is relatively weak, one of the two underlying chiralities in the topological superconducting state can be strongly manifest in the vortex core structure and govern the response of the system to vorticity and a nonmagnetic impurity where the vortex is pinned. The Majorana zero mode in the vortex core is found to be robust against the nonmagnetic impurity in that it remains effectively a zeroenergy bound state regardless of the impurity potential strength and the major chirality. The spin polarization of the Majorana bound state depends on the major chirality for weak spinorbit coupling, while it is determined simply by the vorticity when spinorbit coupling is relatively strong.
pacs:
74.20.z, 74.20.Rp, 74.25.Uv, 75.70.TjI Introduction
One of the most promising models proposed so far for platforms to realise topological quantum computation is the twodimensional (2D) wave topological superconductivity (TSC) model with Rashba spinorbit (SO) coupling and Zeeman field.Sato2009 (); Sato2010 (); Sau2010 (); Alicea2010 (); Alicea2012 () Sato, Takahashi, and Fujimoto have proposed the tightbinding modelSato2009 (); Sato2010 () that can describe an wave superfluid of ultracold fermionic atoms in an optical lattice, where an effective SO interaction can be generated by spatially varying laser fields, or 2D wave TSC in a solid device. Such TSC can be realised, as proposed by Sau et al. in terms of the continuum model,Sau2010 () in a semiconductor heterostructure, where a semiconductor with Rashba SO coupling is sandwiched between a conventional wave superconductor and a ferromagnetic insulator. While wave superconductivity is induced by the proximity effect, the ferromagnetic insulator can generate Zeeman coupling via exchange interactions, thus affecting only the spin degree of freedom in the semiconductor. A vortex in the model hosts a zeroenergy Majorana bound stateSau2010 (); Tewari2010 () and hence vortices in a 2D wave topological superconductor obey nonAbelian exchange statistics, like those in chiral wave superconductors.Kopnin1991 (); Read2000 (); Ivanov2001 (); Stern2004 (); Stone2006 ()
Inherent in the 2D wave TSC model are the two chiralities, in the tightbinding modelFujimoto2008 (); Sato2009pwave () or in the continuum modelZhang2008 (); Alicea2010 (); Shitade2015 () (also for in the tightbinding model), where is the wavevector. The noninteracting Fermi surface is split sideways by the Rashba SO interaction, which causes winding of spin as one goes around each Fermi surface,Frigeri2006 () and the Zeeman field perpendicular to the 2D system favours one spin component along its direction over the other, resulting in two separate bands. With large enough Zeeman splitting and the chemical potential in the gap between the two bands, the system has a single Fermi surface on which spin is fixed for each and thus becomes effectively spinless – condition required for TSC that can support Majorana zero modes.Alicea2012 () Although the Fermi surface has a certain chirality (spin winding), in principle the two intrinsic chiralities are always present in a 2D wave TSC state.Shitade2015 ()
It is known that in a spintriplet chiral wave superconductor, vortices in the and states, which are degenerate at zero field, are not equivalent.Takigawa2001 (); Ichioka2002 () The free energy is lower and the upper critical field is higher when the vorticity, i.e., the angular momentum carried by the supercurrent, is antiparallel to the chirality, i.e., the angular momentum of Cooper pairs in the condensate.Ichioka2002 () Moreover, bound states in such an antiparallel vortex are more robust against nonmagnetic impurities due to cancellation of angular momenta between the supercurrent and the condensate; which makes the vortex core region wavelike and allows the Anderson theoremAnderson1959 () to take effect.Kato2000 (); Kato2002 () Most notably, in each of the chirality domains, the order parameter of the other chirality is induced around the vortex centre, more prominently for lower applied field.Takigawa2001 (); Ichioka2002 ()
The effects of nonmagnetic impurities on vortex bound states in a 2D wave topological superconductor, depending on the major chirality with respect to the vorticity, have been studied recently for a single vortex in the continuum modelSau2010 () by means of a Greenfunction technique for calculating the impurity selfenergy.Masaki2014 (); Masaki2015 () It has been found that when the major chirality is opposite to the vorticity, vortex bound states are less influenced by nonmagnetic impurities for relatively weak SO coupling, compared to the case where the major chirality is in the same direction as the vorticity. As to spinless superconducting states per se, Ivanov has stated that though with slightly different structure of the quasiparticle eigenfunctions, the two types of vortices have the same lowenergy spectra and braiding statistics.Ivanov2001 ()
Volovik has studied the effects of a single nonmagnetic impurity in a spinful chiral wave superconductor, where Cooper pairs have a definite angular momentum, using quasiclassical theoryVolovik1999 () and has found that the Majorana zero mode in the vortex core is not affected by a nonmagnetic impurity. The robustness of the Majorana fermion at zero energy has also been shown for a coreless vortex in an “antidot” on the surface of a disordered threedimensional (3D) topological insulator,Ioselevich2012 () where superconductivity is induced by proximity to a superconducting film with a circular hole with radius larger than the coherence length and the mean free path. Such a 3D TSC state is in the symplectic class of AII, where the topological invariant belongs to , and the system has no particlehole symmetry and is odd under time reversal.Schnyder2008 () It is intriguing to study the effects of a nonmagnetic impurity in a 2D wave topological superconductor in the presence of vortices, which belongs to symmetry class D and has particlehole symmetry, but no timereversal symmetry.Schnyder2008 () In this system there exist two underlying chiralities ( in the continuum model) and the effective (spinless) wave nature of the system varies depending on the material parameters.
The purpose of the present work is to examine the effects of the two chiralities inherently present in 2D wave TSC states on vortex structure, by solving the Bogoliubovde Gennes (BdG) equationsdeGennes () on the tightbinding model of Sato, Takahashi, and FujimotoSato2009 (); Sato2010 () selfconsistently for the superconducting order parameter in the vortex lattice. The tightbinding model is versatile and useful for modelling real systems in that band structure and the filling factor can easily be incorporated in terms of hopping amplitudes and chemical potential. We assume that the 2D system has Rashba SO coupling and a pairing interaction that drives wave superconductivity or superfluidity, and is under Zeeman field, e.g., generated by proximity to a ferromagnetic insulator in a heterostructure. Our model is applicable to systems such as wave superfluids of fermionic atoms created by wave Feshbach resonance in an optical latticeSato2009 (), oneatomlayer TIPb on Si(111)Matetskiy2015 (), and ionicliquid based electronic doublelayer transistors.Li2016 (); Nagai2016 ()
Solving the BdG equations for TSC has high numerical demand as the dimension of the BdG Hamiltonian matrix is four times the total number of lattice sites, and also the system area needs to be large enough to sustain two vortices that are well separated to have a pair of isolated Majorana fermions as vortex bound states. Thus, the conventional way of solving the BdG equations by direct diagonalization is not feasible. We utilise the Chebyshev polynomial expansion methodCovaci2010 (); Nagai2012 () for solving for the order parameter selfconsistently, as well as calculating the local density of states (LDOS) after selfconsistency has been achieved. Furthermore, we use the numerically efficient algorithm developed by Sakurai and Sugiura Sakurai2003 (); Nagai2013 () to obtain quasiparticle spectra within an energy window of one’s choice.
Ii Model
We use the tightbinding model for a 2D wave topological superconductor:Sato2009 (); Sato2010 (); Nagai2015 ()
(1)  
where we consider hopping among nearestneighbour lattice sites only with the hopping amplitude , is the chemical potential, is the singleparticle potential due to a nonmagnetic impurity at site , is the Zeeman field, is the Rashba SO coupling strength, is the wave superconducting order parameter at site , and H.c. stands for the Hermitian conjugate. We set the lattice constant to be unity, and and are the unit vectors in the and directions. and creates and annihilates, respectively, the electron at site with spin (). We solve the BdG equations with the Hamiltonian (1) and solve for the order parameter selfconsistently for a given coupling constant for the pairing interaction, :
(2) 
When the system has translational symmetry, the realspace Hamiltonian in Eq. (1) can be Fouriertransformed to momentum space and written asSato2010 ()
(3) 
where and
(4) 
Here and are the creation and annihilation operators of the electron with momentum and spin , , , and and are the Pauli matrices. The above Hamiltonian in momentum space can be diagonalized to obtain the quasiparticle spectrum as
(5) 
where and we have assumed an isotropic wave order parameter, . Depending on the values of , , and , the system can be in a trivial or nontrivial (Abelian or nonAbelian) topological phase according to the first Chern number or the ThoulessKohmotoNightingaleNijs (TKNN) number,TKNN1982 () which we denote as , as classified for four different band regions in Ref. Sato2010, . The topological invariant Schnyder2008 () can be calculated byIshikawa1987 (); Volovik2009 (); Gurarie2011 ()
(6)  
where . The spectral gap is the minimum value of in Eq. (5) and e.g., as is varied for a given set of , , and , the system transitions from one topological (trivial or nontrivial) phase to another every time vanishes. The system is in Abelian and nonAbelian phase when is even and odd, respectively ( and in this model), and in trivial phase when . Achieving nonAbelian phase with in most of the band region would require relatively large Zeeman fieldSato2010 (); Nagai2014oval () and accordingly large values of and/or that may be unrealistic for actual materials in order to overcome the Pauli depairing effect. Thus, in this work we focus on nonAbelian states with in the band regions and , whereSato2010 ()
(7)  
(8) 
The vortex lattice is formed by uniform magnetic field applied in the direction, . The electron wavefunction acquires the Peierls phase factor while traversing from site to (coordinate from to ) due to the associated vector potential so that the hopping amplitude is modified as
(9) 
where the vector potential in the symmetric gauge. To ensure obtaining Majorana fermions that come in pairs as a solution to the BdG equations,Stone2004 (); Sato2009pwave () we place two vortices in our system that is the vortex unit cell, i.e., two flux quanta within the system area. The field strength is controlled by the system size as with , where and are the number of lattice sites in the and directions. For the results presented below, we impose the periodic boundary condition for the vortex latticeTakigawa2000 (); Kawakami2016 () such that there is one vortex in the centre of a square lattice and a quarter vortex at each of the four corners centred right outside the corner site, e.g., at . For all the results shown, an odd number of lattice sites has been used so that the vortex inside the system is centred at the centre site of the square lattice. We have used the convention that the electron charge is negative () and hence the quasiparticle bound states carry the angular momentum of about the vortex centre. The phase winding of the order parameter is referred to as the vorticity ( in our calculation) hereafter.
To study the effects of nonmagnetic impurities, we place a single nonmagnetic impurity with potential at the centre site of the system. This is where one of the vortices is centred at, which is a reasonable assumption as a vortex tends to be pinned by an impurity or a defect in real materials, and can be thought of as a pinning potential. We stop selfconsistent iterations for the order parameter at the th iteration step, when the order parameter as a complex vector of length satisfies
(10) 
where we set . We have found that results do not change with tighter convergence criteria with as small as .
(a)

(b)

The energy spectrum (5) depends only on the magnitude of the Zeeman field and so the phase transition between different topological phases is independent of the sign of . In a given topological phase, however, which one of the two underlying chiralities, ,Fujimoto2008 (); Sato2009pwave () is more manifest is determined by the sign of Wu2012 () as well as the sign of .Sato2010 () This can be seen by expressing in the “chirality basis” that diagonalizes the normalstate Hamiltonian :Sato2010 ()
(11) 
where and
(12) 
In the normal state, the eigenspectrum in Eq. (5) reduces to , where in is necessary for obtaining the correct eigenvalues and eigenfunctions as can be seen by taking the limit . We note that this factor is missing in the discussion of the chirality basis in Ref. Sato2010, where was supposed to be positive. It can be seen in Eq. (12) that the intraband pairing in the band () has the chirality of (), while the interband pairing is purely wave. This is analogous to the two chiralities present in the nonAbelian phase of the continuum model.Alicea2010 (); Shitade2015 (). Although orbital angular momentum is not a good quantum number in the presence of SO coupling, it has been found in Ref. Shitade2015, that when is small enough, the average angular momentum carried by a Cooper pair can be close to in nonAbelian states dominated by , as in chiral wave superconductors.
We illustrate in Fig. 1 for and the two cases, (a) and (b) , where the normalstate Fermi surface is formed by and , respectively. This value of is small enough so that the respective chirality, and , associated with the Fermi surface can be dominant over the other, as discussed in Sec. III.1 for (where the dominant chirality is and , respectively, for and ). Due to the interband pairing, however, both chiralities are present in a TSC state in general. We will demonstrate in Secs. III.1 and III.3 that for relatively weak SO coupling, vortex structure and effects of a nonmagnetic impurity can be influenced strongly by one of the chiralities and , which we call the chirality of and , respectively, in the remainder of the paper.
The major chirality in a nonAbelian TSC state stems from the chirality (spin winding) of the single Fermi surface in the normal state.Sato2010 (); Alicea2010 (); Shitade2015 (). The latter is for , and , , and for , and , : The first and second combinations correspond to the convention used in the continuum model in Refs. Masaki2014, , Masaki2015, and Ref. Shitade2015, , respectively. In the case of a trivial TSC state, the chirality of the dominant one of the two Fermi surfaces can manifest itself. The TKNN number reverses its sign under , whereas it does not under . Whether trivial or nontrivial, the underlying chirality in the TSC state – rather than the TKNN number itself – governs how the system responds to vorticity and nonmagnetic impurities in the presence of a vortex.
Iii Results
iii.1 Chirality vs Vorticity
We have studied the vortex lattice states for a wide variety of parameter sets (, , , and the resulting bulk order parameter and spectral gap ) within our tightbinding model (1), in particular exploring the effects of the strength of SO coupling and the sign of the Zeeman field . We find that for the sign of hardly makes any difference in a pure vortex lattice state (i.e., with no impurity). The converged order parameter, the excitation spectrum, and the groundstate energy are practically the same for positive and negative in the absence of impurity. We will discuss the influence of chirality on the vortex structure that can be apparent even in a pure system for relatively small in Sec. III.3.
In Fig. 2 we show the quasiparticle spectra for in the vortex lattice in a pure 5151site system for , , and (, , and ), where the abscissa is an index numbering the discrete eigenvalues. The major change caused by reversal of the Zeeman field while keeping all other parameters the same is the average numbers of spinup and spindown electrons in the system being interchanged. Accordingly, the overall magnitude of the spinup and spindown probability amplitudes of the Majorana bound state (and lowenergy quasiparticle excitations) in each vortex core changes when the sign of is flipped: e.g., if the spinup component of a Majorana is dominant over the spindown component for a given , then vice versa for . Other than this, the overall structure of each of the spinup and spindown Majorana wave functions is little affected by the direction of the Zeeman field. It is just barely discernible in Fig. 2 that the first excited state has a slightly lower energy for than for . We find that the opposite is true for positive , though the difference tends to be very small. The Majorana boundstate energy is and , respectively, for and in this example. These energy levels are not exactly zero due to nonzero overlap of the Majorana wavefunctions bound to nearestneighbour vortices on the finite lattice. To demonstrate that they are indeed Majorana bound states, the minimum eigenvalue in units of the bulk spectral gap is plotted in Fig. 3 as a function of the distance between nearestneighbour vortices for , 41, 51, 55, and 61; for the two systems shown in Fig. 2 and for , , , and (, , and ). It can be seen that the minimum eigenvalue approaches zero as the system size increases.
Presented in Fig. 4 is the magnitude of the converged order parameter of the vortex lattice in a 4141site system for , , , ( and ) as a function of and coordinates; with (a) no impurity and (b) a single nonmagnetic impurity with potential at the centre of the lattice, where one of the vortices is centred at. The contour projection of the order parameter magnitude onto the plane is in steps of . For this value of , the vortex structure in a pure system barely depends on the sign of , and for is very similar to the one shown in Fig. 4(a) for . In the presence of a nonmagnetic impurity, however, one of the two chiralities can manifest itself in the vortex core structure. For , while the vortex structure of the order parameter does not change much from the pure case for (not shown), it can be seen clearly in Fig. 4(b) that for , the vortex core shrinks in comparison to the noimpurity case shown in Fig. 4(a).
The LDOS at the vortex centre (where the nonmagnetic impurity is placed at) is plotted as a function of quasiparticle energy in Fig. 6 for the system presented in Fig. 4 and in Fig. 5 for the same system but for , where the (a) spinup and (b) spindown components of the LDOS are compared between and . The Lorentz kernelWeisse2006 (); Covaci2010 () has been used for calculation of the LDOS with the corresponding Lorentzian smoothing width of . The number of terms summed over in the Chebyshev polynomial expansion of the Green function components for the LDOS is 3000, while the cutoff of 2000 terms was used for the order parameter. The bulk spectral gap in this system (irrespective of the sign of ) is . We first note that the Majorana bound state is dominated by its spindown component at the vortex centre for both [Figs. 5(b) and 6(b)]. It turns out to be always the case for that spin of the Majorana zero mode is dominantly down at the vortex centre, regardless of the sign of the Zeeman field or the chemical potential. As mentioned above in regard to the Zeeman field, one spin component can be dominant over the other in a Majorana for a given and , and the dominance is reversed under the change , or . Whether up spin or down spin is dominant in the wave function overall, spin of the Majorana bound state is mostly down at the vortex centre for .
This spin polarization at the vortex centre is consistent with the findings of Ref. Nagai2014top, , where Majorana bound states in a threedimensional topological superconductor with spinorbit coupling have been found to be spinpolarised in a vortex core. Assuming circular symmetry around the vortex centre, the spinup and spindown wavefunctions are given by the Bessel function of the first kind, and , respectively, where and are the orbital angular momentum carried by the spinup and spindown components of a quasiparticle, for the chirality of .Nagai2014top () The Majorana bound state has either or , whichever yields the lowest (zero) energy, depending on the vorticity. In Ref. Nagai2014top, the vorticity of was used that resulted in the Majorana wavefunction dominated by its spinup component in the vortex core.
It can be seen in Fig. 5 that both spinup and spindown (Majorana) boundstate energies are hardly affected by the impurity for . This can be understood in terms of the major chirality being : the angular momentum carried by Cooper pairs in the condensate is mostly cancelled by the angular momentum carried by the supercurrent, making the system effectively wavelike in the vortex core region. In contrast, the spinup (nonMajorana) boundstate energies are shifted substantially by the presence of the impurity for [Fig. 6(a)], where the chirality has the same sign as the vorticity. Also in this case, however, the Majorana bound state is robust in that its energy is unchanged, albeit with reduced amplitude at the impurity site [Fig. 6(b)].
Figures 7(a) and 7(b) show the quasiparticle spectra for the systems presented in Figs. 5 and 6, respectively. It can be seen that even for , the first excited state is pushed up in energy though slightly, resulting in an increased minigap, by the nonmagnetic impurity at the vortex centre. Note that each energy level is doubly degenerate as there are two vortices in the lattice, though numerically not exactly as the two vortices are not equivalent: Only one of them has its centre in the system, where we place a single nonmagnetic impurity, hence shifting the energy of only one of the two states (in each level). The only perceptible change caused by the impurity for is the small shift of the first excited level, and this is reflected in the LDOS in Fig. 5(a) as the slight shift in the peak position. For , on the other hand, the first few excited states are affected by the impurity, and this results in the substantial shift of the boundstate peak in the LDOS at the vortex centre, as can be seen in Fig. 6(a).
For both signs of , the energy of the Majorana bound state remains practically the same ( on the 4141 lattice) with or without the impurity. This turns out to be the case also for stronger impurity potential, or stronger spinorbit coupling as discussed in Sec. III.2, where the vortex bound states are significantly affected by the impurity for as well as for . We also find that the minigap is increased by the presence of a nonmagnetic impurity at the vortex centre, regardless of the (nonzero) magnitude, or the sign of the potential – in case of a positive potential, as long as the magnitude is large enough so that the order parameter is suppressed at the impurity site in the absence of vortices. The order parameter can be enhanced at the impurity site for a relatively weak, positive potential,Hu2013 (); Goertzen2016 () and the BdG equations tend to have difficulty converging when such an impurity is placed at the vortex centre, especially for presumably due to competition between the dominant chirality and the vorticity. By including a Gaussian impurity potential in the core of a single vortex in the continuum model,Sau2010 () the authors of Ref. Mao2010, showed that the minigap simply increased continuously as the magnitude of the positive Gaussian potential was increased from zero. This result may have been possible because these authors did not solve the BdG equations selfconsistently for the vortex state within the topological superconductor.
iii.2 Strong SO Coupling
For the systems studied in the previous subsection, the chemical potential is fairly close to the top of the band and also the SO coupling constant is relatively small so that the Fermi surface in the normal state is nearly perfectly circular. This means that the system has rotational symmetry to a good approximation and the difference between the two intrinsic chiralities can be apparent in vortex states. In this subsection, we examine a case where the chemical potential is far away enough from the top or bottom of the band and/or SO coupling is relatively strong so that the Fermi surface is nowhere close to having circular symmetry. For this purpose we use and . The normalstate Fermi surface is shown in Fig. 8, which is centred at like the Fermi surface () for and (not shown). One can see that a large portion of this Fermi surface is almost flat and perpendicular to the directions .
The converged order parameter for , , and (, , and ) in the 4141 vortex lattice with at the centre of the system is shown in Fig. 9 for (a) and (b) . The contour projection onto the plane is again in steps of . It can be seen that the vortex core size is not much different for the positive and negative in this case: It also turns out that it does not change much from the pure case () for either sign of . Furthermore, in contrast to the results for , , and in Fig. 4, the contour projection shows that the order parameter has more squarelike symmetry rather than circular about the vortex centre. We note that it is more so for negative than positive , and even in the result for shown in Fig. 4 the order parameter has square symmetry ( rotated with respect to the axes) very close to the vortex centre.
Figures 10 and 11 present the (a) spinup and (b) spindown components of the LDOS at the vortex centre as a function of quasiparticle energy for the systems shown in Fig. 9, in comparison with the pure case. The influence of the impurity is visible for both , with some energy shifts for spinup quasiparticles. In fact, vortex bound states are more affected by the impurity for positive than for negative : the shift of the spinup LDOS peak for is a reflection of the first excited state moved up in energy from to , thus increasing the minigap by about . Interestingly, unlike the results for , , and in Figs. 5 and 6, the Majorana bound state for and some of the spinup bound states for become more bound to the vortex centre by the presence of the impurity. Yet, the energy of the Majorana bound state ( in this example) is barely changed by the impurity for both signs of .
We show in Figs. 12 and 13 the Majorana wave functions for the systems presented in Figs. 10 and 11, respectively. In each of Figs. 12 and 13, the spinup (a) particle and (b) hole probability amplitudes and the spindown (c) particle and (d) hole probability amplitudes for are plotted as a function of and coordinates for the entire 4141site system. Respective probability amplitudes are plotted for in (e) and (f), and (g) and (h). We first note that the particle and hole probability amplitudes in a given spin component are practically identical for all the cases shown, as expected for a Majorana fermion. Fourfoldsymmetric extension of the wave functions can be discerned especially for the spinup components (and the spindown components for and ), reflecting in the Fermi wave vector. The boundstate energy for this finitesize system being slightly higher than the counterpart for , , and mentioned at the end of Sec. III.1 can be attributed to the nonzero overlap of the Majorana wave functions among the nearestneighbour vortices.
Masaki and Kato have foundMasaki2014 (); Masaki2015 () in terms of the continuum modelSau2010 () that for weak SO coupling, vortex bound states are more robust against nonmagnetic impurities when the chirality is opposite to the vorticity, compared to the case where the chirality and vorticity are in the same direction, and that this difference between the two chiralities diminishes as SO coupling is made stronger. They have also found that the scattering rates of zeroenergy bound states are very small regardless of the major chirality and the strength of SO coupling.Masaki2015 () Our results, though only with a single nonmagnetic impurity at the vortex centre, are consistent with their findings.
iii.3 Weak SO Coupling
In this subsection, we illustrate that when SO coupling is relatively weak, vortex structure can be markedly different for the two opposite directions of the Zeeman field even in a pure vortex state. It has been found in Ref. Shitade2015, that the weaker the SO coupling, the closer the average angular momentum per Cooper pair to in the ()dominated states. Thus, in terms of the direction of the Zeeman field (or alternatively the external field applied to create the vortex lattice) one can make one of the two intrinsic chiralities strongly manifest in the vortex core structure.
Shown in Fig. 14 is the order parameter in the pure 5151 vortex lattice for , , and (, , and ); for (a) and (b) . Compared with the order parameter for and shown in Fig. 4(a) (for , but it is very similar to the order parameter for ), it can be seen in Fig. 14(a) that for , the coherence length in the sense of a recovery length of the order parameter from the vortex centre to the bulk value is much larger for than . Note that the system size, namely, the size of the unit cell of the vortex lattice is 4141 and 5151, respectively, in Figs. 4(a) and 14(a). In fact, we have tried the coupling constant of for that yields a similar bulk order parameter as for ; however, the coherence (recovery) length in this case is so large that there is substantial overlap of vortex cores and it is not useful for comparison of the vortex structure with the case.
Striking is the difference between the and cases for [Figs. 14(b) vs 4(b)] and also between positive and negative in Figs. 14(a) and 14(b). As can be seen in Fig. 14(b), for as the order parameter tends toward zero at the vortex centre, it oscillates and is enhanced slightly around the vortex centre. At first glance, this is reminiscent of the order parameter induced in the vortex core with vorticity in the domain of a chiral wave superconductor.Ichioka2002 () In our calculation, however, the vorticity is and thus parallel to the major chirality for . Moreover, we only have the superconducting order parameter stemming from spinsinglet wave pairing, and the enlargement of the vortex core for and the enhancement around the vortex centre for are happening within the same wave order parameter. Both underlying chiralities are always present and mixed except for ,Shitade2015 () and it is not clear as to whether there is a way to define a unique order parameter for each of the two chiralities separately. It is apparent nonetheless that some extra order is induced in the system inside the vortex core, increasing the superconducting order somewhat.
The suppression and enhancement of the order parameter around the vortex centre for and , respectively, seen in Figs. 14 (a) and 14(b) are reflected in the quasiparticle spectra presented in Fig. 15. One can see that the first few excited levels of the vortex bound states are higher for than for , with the minigap of and , respectively ( for for both ). By comparison, the Majorana boundstate energy is not much different between () and (). The groundstate energy of the system and equivalently the average energy gain per electron are also hardly different for and . Björnson and BlackSchaffer have solved the BdG equations on the Hamiltonian (1) selfconsistently (without the vector potential) for a vortex in a square lattice with open boundaries and have also found asymmetry in lowenergy spectra between positive and negative for a given vorticity for .Bjornson2013 ()
The Majorana wavefunctions are shown in Fig. 16 for the above two systems; where the spinup (a) particle and (b) hole probability amplitudes and the spindown (c) particle and (d) hole probability amplitudes for , and respective probability amplitudes for in (e) and (f), and (g) and (h), are plotted as a function of and coordinates for the entire 5151 lattice. Once again, the particle (left column) and hole (right column) probability amplitudes are virtually identical in each spin component. For , the boundstate peak of the spindown component is strongly localized at the vortex centre [(c) and (d)], while the spinup component, though much smaller in magnitude, extends substantially outside the core region [(a) and (b)].
Remarkably, the spindown probability amplitudes are not peaked at the vortex centre for , as clearly visible in Figs. 16(g) and 16(h). They are peaked at surrounding sites and though small in magnitude, have substantial extension outside the vortex core. The spinup amplitudes are a little more confined: They are also suppressed right at the vortex centre, although it is not discernible in Figs. 16(e) and 16(f). Furthermore, the spinup amplitudes are one order of magnitude larger than the spindown amplitudes, contrary to the case as well as the systems in Figs. 12 and 13 with strong SO coupling. The LDOS at the vortex centre as a function of quasiparticle energy shown in Fig. 17(b) confirms the absence of the spindown Majorana component at the vortex centre for , in stark contrast to the examples shown in Figs. 5, 6, 10 and 11, as well as the LDOS for in Fig. 17(a). The peak in the spinup LDOS for seen in Fig. 17(b) corresponds to the first excited state (the minigap at ).
Placing a nonmagnetic impurity at the vortex centre can enhance the unusual structure of the vortex core for . This is illustrated in Fig. 18, where the order parameter is shown for the systems analogous to those presented in Fig. 14, but with at the centre of the 4141 vortex lattice, for (a) and (b) . For , compared to the pure system (with 4141 lattice sites), the size of the vortex core is little changed by the nonmagnetic impurity, while its shape is modified slightly. For , on the other hand, the enhancement of the order parameter around the vortex centre is exaggerated by the presence of the impurity so much that the order parameter at its peaks is significantly larger than the bulk value.
Finally, we demonstrate the manifestation of chirality in the vortex lattice in the trivial phase. Shitade and NagaiShitade2015 () have found in the ()dominated states that when , the average angular momentum per Cooper pair can be close to already in the trivial phase as () approaches the critical value for the phase transition from the trivial to nonAbelian phase. We present in Fig. 19 the order parameter in the pure 4141 vortex lattice for , , and , for (a) and (b) . These systems are in the trivial phase with . The coupling constant has been chosen so as to make the bulk order parameter () similar to that in the nonAbelian systems for , , and discussed above. It can be seen in Figs. 19(a) and 19(b) that the suppression and enhancement of the order parameter around the vortex centre for and , respectively, are more substantial than in the nonAbelian systems shown in Figs.14(a) and 14(b). This signifies the fact that the angular momentum carried by Cooper pairs governs vortex core structure, even though the current TSC model can be mapped onto a spinless wave superconductor with chiralities only in the nontrivial phase.Sato2010 ()
iii.4 Minigap
We show in Fig. 20 the quasiparticle spectra for and in the pure 5151 vortex lattice; for (red plus), (blue cross), and (magenta star). The coupling constant has been chosen to be and , respectively, for and so that the bulk order parameter as for (). The bulk spectral gap (independent of the sign of ) is , , and for , , and , respectively. As mentioned in relation to Fig. 2 in Sec. III.1, in these systems and the minigap is larger for than for by a very small amount (of the order of ). The minigap is , , and , respectively, for , , and .
For these and other systems we have examined, we have not found any direct correlation between the minigap and the bulk order parameter nor the spectral gap. In this particular example (also for ), simply the smaller the , the larger the minigap. The formula for the minigap for wave or chiral wave superconductors, , where is the magnitude of the Fermi wave vector and is the coherence length,deGennes (); Tinkham () has been used in the literature for the wave TSC model with Rashba SO coupling and Zeeman field. For a given and , is determined uniquely from Eq. (5) once is fixed. However, also depends on and . Thus, it is unclear as to whether this common formula applies to the minigap in wave TSC vortices with the two inherent chiralities.
As discussed at the end of Sec. III.1, the minigap is increased by a nonmagnetic impurity which the vortex is pinned by, and the increase in the minigap can be substantial as illustrated in Figs. 6 and 10, depending on the major chirality and the strength of SO coupling. The larger the minigap, the shorter the minimum time required for braiding operation of vortices to avoid nonadiabatic transitions of the Majorana zero mode to excited states.Masaki2015 () Therefore, further exploration of the effects of and on the minigap, especially in the presence of a nonmagnetic impurity in the vortex core, is of great interest from the application point of view for topological quantum computation.
Iv Conclusions
In summary, we have performed a selfconsistent study of the vortex lattice in a twodimensional topological superconductor with an wave pairing interaction, Rashba spinorbit coupling, and Zeeman field. When a vortex is pinned at a nonmagnetic impurity, one of the two intrinsic chiralities in the nonAbelian topological superconducting state can manifest itself in the vortex structure and affect the lowenergy excitation spectrum for relatively weak spinorbit coupling. One can make one chirality dominant over the other by changing the direction of the Zeeman field, or alternatively, the sign of the chemical potential. In such states, (nonMajorana) vortex bound states are less influenced by the impurity if the dominant chirality is opposite to the vorticity, compared to the case where it is in the same direction as the vorticity. For stronger spinorbit coupling, where orbital angular momentum is less of a “good quantum number”, lowenergy spectra tend to be more affected by the impurity regardless of the direction of the Zeeman field. The Majorana zero mode effectively remains a zeroenergy bound state, and its spin is polarised according to the vorticity unless spinorbit coupling is rather weak, in which case the spin polarization depends on the major chirality as well.
We have shown in Sec. III.3 that when spinorbit coupling is weak, the vortex core structure can be strikingly different for the two directions of the Zeeman field: The order parameter is suppressed and enhanced around the vortex centre when the major chirality is antiparallel and parallel, respectively, to the vorticity. We have demonstrated this phenomenon not only in nonAbelian topological phase, but also in trivial phase where the TKNN number is zero. The major chirality in such a trivial TSC state is the chirality of the dominant (or single) Fermi surface. The enhancement of the order parameter implies some extra superconducting order being induced. This occurs, however, in the vortex core region where the chirality and the vorticity are in the same direction, and appears to be counterintuitive in view of the vortex physics in spintriplet superconductors. Moreover, it is not immediately obvious how to map out chiralitybased order parametersBjornson2015 (); Shitade2015 () from the sole order parameter in the model that is of wave pairing symmetry. The chiral nature of 2D wave TSC states and the role that the angular momentum carried by Cooper pairs plays in determining various properties of the system are to be explored in further detail in a future publication.
Finally, the suppression (the enlargement of the vortex core) and enhancement of the selfconsistent order parameter around the vortex centre will be reflected in the supercurrent and hence the field distribution in the vortex core area. Thus, particularly in a 2D wave topological superconductor where spinorbit coupling is relatively weak, manifestation of different chiralities can be detected by probing the field distribution in the vortex lattice, e.g., by NMR and by switching the direction of the applied magnetic field.
V Acknowledgements
We thank S. L. Goertzen and T. Kawakami for helpful discussions. K.T. is grateful to CCSE at Japan Atomic Energy Agency for hospitality, where part of the research was performed. This work was enabled in part by support provided by WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca). Part of the calculation was also performed on the supercomputing system SGI ICE X at the Japan Atomic Energy Agency. The research was supported by the Natural Sciences and Engineering Research Council of Canada, and partially by JSPS KAKENHI Grant No. 26800197 and the “Topological Materials Science” (No. 16H00995) KAKENHI on Innovative Areas from JSPS of Japan.
References
 (1) M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. Lett. 103, 020401 (2009).
 (2) M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. B 82, 134521 (2010).
 (3) J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Phys. Rev. Lett. 104, 040502 (2010).
 (4) J. Alicea, Phys. Rev. B 81, 125318 (2010).
 (5) J. Alicea, Rep. Prog. Phys. 75, 076501 (2012).
 (6) S. Tewari, J. D. Sau, S. Das Sarma, Ann. Phys. 325, 219 (2010).
 (7) N. B. Kopnin and M. M. Salomaa, Phys. Rev. B 44, 9667 (1991).
 (8) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
 (9) D. A. Ivanov, Phys. Rev. Lett. 86, 268 (2001).
 (10) A. Stern, F. von Oppen, and E. Mariani, Phys. Rev. B 70, 205338 (2004).
 (11) M. Stone and S.B. Chung, Phys. Rev. B 73, 014505 (2006).
 (12) S. Fujimoto, Phys. Rev. B 77, 220501(R) (2008).
 (13) M. Sato and S. Fujimoto, Phys. Rev. B 79, 094504 (2009).
 (14) C. Zhang, S. Tewari, R. M. Lutchyn, and S. Das Sarma, Phys. Rev. Lett. 101, 160401 (2008).
 (15) A. Shitade and Y. Nagai, Phys. Rev. B 92, 024502 (2015).
 (16) P. A. Frigeri, D. F. Agterberg, I. Milat, and M. Sigrist, Eur. Phys. J. B 54, 435 (2006).
 (17) M. Takigawa, M. Ichioka, K. Machida, M. Sigrist, Phys. Rev. B 65, 014508 (2001); see also references therein.
 (18) M. Ichioka and K. Machida, Phys. Rev. B 65, 224517 (2002).
 (19) P. W. Anderson, J. Phys. Chem. Solids 11, 26 (1959).
 (20) Y. Kato, J. Phys. Soc. Jpn. 69, 3378 (2000).
 (21) Y. Kato and N. Hayashi, J. Phys. Soc. Jpn. 71, 1721 (2002).
 (22) Y. Masaki, Y. Kato, J. Phys.: Conf. Ser. 568, 022028 (2014).
 (23) Y. Masaki and Y. Kato, J. Phys. Soc. Jpn. 84, 094701 (2015).
 (24) G. E. Volovik, JETP Lett. 70, 609 (1999).
 (25) P. A. Ioselevich, P. M . Ostrovsky, and M. V. Feigelâman, Phys. Rev. B 86, 035441 (2012).
 (26) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
 (27) P. G. de Gennes, Superconductivity of Metals and Alloys (Westview Press, Boulder, 1999).
 (28) A. V. Matetskiy, S. Ichinokura, L. V. Bondarenko, A. Y. Tupchaya, D. V. Gruznev, A. V. Zotov, A. A. Saranin, R. Hobara, A. Takayama, and S. Hasegawa, Phys. Rev. Lett. 115, 147003 (2015).
 (29) L. J. Li, E. C. T. O’Farrell, K. P. Loh, G. Eda, B. Özyilmaz, and A. H. Castro Neto, Nature (London) 529, 185 (2016).
 (30) Y. Nagai, S. Hoshino, and Y. Ota, Phys. Rev. B 93, 220505 (2016).
 (31) L. Covaci, F. M. Peeters, and M. Berciu, Phys. Rev. Lett. 105, 167006 (2010).
 (32) Y. Nagai, Y. Ota, and M. Machida, J. Phys. Soc. Jpn. 81, 024710 (2012).
 (33) T. Sakurai and H. Sugiura, J. Comput. Appl. Math. 159, 119 (2003).
 (34) Y. Nagai, Y. Shinohara, Y. Futamura, Y. Ota, and T. Sakurai, J. Phys. Soc. Jpn. 82, 094701 (2013).
 (35) Y. Nagai, Y. Ota, and M. Machida, J. Phys. Soc. Jpn. 84, 034711 (2015).
 (36) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
 (37) K. Ishikawa and T. Matsuyama, Nucl. Phys. B 280, 523 (1987).
 (38) G. E. Volovik, The Universe in a Helium Droplet (Oxford University Press, Oxford, 2009).
 (39) V. Gurarie, Phys. Rev. B 83, 085426 (2011).
 (40) Y. Nagai, Y. Ota, and M. Machida, J. Phys. Soc. Jpn. 83, 094722 (2014).
 (41) M. Stone and R. Roy, Phys. Rev. B 69, 184511 (2004).
 (42) M. Takigawa, M. Ichioka and K. Machida, J. Phys. Soc. Jpn. 69, 3943 (2000).
 (43) T. Kawakami, private communication.
 (44) L.H. Wu, Q.F. Liang, Z. Wang and X. Hu, J. Phys.: Conf. Ser. 393, 012018 (2012).
 (45) A. Weiße, G. Wellein, A. Alvermann and H. Fehske, Rev. Mod. Phys. 78, 275 (2006).
 (46) Y. Nagai, H. Nakamura, and M. Machida, J. Phys. Soc. Jpn. 83, 064703 (2014).
 (47) H. Hu, L. Jiang, H. Pu, Y. Chen, and X.J. Liu, Phys. Rev. Lett. 110, 020401 (2013).
 (48) S. L. Goertzen, Y. Nagai, and K. Tanaka (unpublished).
 (49) L. Mao and C. Zhang, Phys. Rev. B 82, 174506 (2010).
 (50) K. Björnson and A. M. BlackSchaffer, Phys. Rev. B 88, 024501 (2013).
 (51) M. Tinkham, Introduction to Superconductivity (McGrawHill, New York, 1996).
 (52) K. Björnson and A. M. BlackSchaffer, Phys. Rev. B 91, 214514 (2015).