Vortex Core Structure in Multilayered Rashba Superconductors

We numerically study the electronic structure of a single vortex in two dimensional superconducting bilayer systems within the range of the mean-field theory. The lack of local inversion symmetry in the system is taken into account through the layer dependent Rashba spin-orbit coupling. The spatial profiles of the pair potential and the local quasiparticle density of states are calculated in the clean spin-singlet superconductor on the basis of the quasiclassical theory. In particular, we discuss the characteristic core structure in the pair-density wave state, which is spatially modulated exotic superconducting phase in a high magnetic field.


Introduction
Intensive studies on locally noncentrosymmetric systems (LNCS) has been conducted [1,2]. The LNCS are realized in various materials such as multilayered cuprates [3], the doped topological insulator Cu x Bi 2 Se 3 [4], the pnictide superconductor SrPtAs [5] and so on. Especially, the heavy-fermion CeCoIn 5 /YbCoIn 5 superlattice is fabricated [6] and the artificial control of the inversion symmetry breaking was realized recently [7]. Therefore, one can seek experimentally the theoretically predicted exotic superconducting phase in a magnetic field [8,9]. The inhomogeneity due to a magnetic field is not taken into account in the existing theoretical researches on the superconductivity in LNCS [8,9]. However, the exotic superconducting phase stabilizes in a magnetic field. So understanding the vortex states in LNCS is necessary. The observation of quasiparticle states in the vicinity of a vortex gives us a lot of information on the exotic superconducting phase. We focus on the exotic superconducting phase called the pair-density wave (PDW) phase, which is stabilized in a high magnetic field perpendicular to the layers [8]. In the PDW phase, the phase of the superconducting order parameter changes by π between the two layers, whereas no phase difference of the order parameter in the BCS phase in a low magnetic field. In this paper, we investigate the vortex core structure of the PDW and the BCS phases, and try to obtain the characteristics of the PDW phase.

Formulation
We consider the multilayered system characterized by the layer-dependent spin-orbit coupling (SOC) strength α m with the layer index m. In this paper, the number of the layer is fixed to N = 2. The layer dependence of the SOC (α 1 , α 2 ) = (α, −α) reflects the local noncentrosymmetricity of the system [10]. We assume the simple form of the Rashba type SOC described by g(k) = (−k y , k x , 0)/k F , which satisfies the normalization condition g(k) k = 1 on the two dimensional Fermi surface (FS). k F is the Fermi wave number and · · · k is the average over the FS.
We consider spin-singlet s-wave pairing states within the layer. The pairing potential in the spin space is expressed as ∆ s 1 s 2 (r) = ∆(r)(iσ y ) s 1 s 2 . We assume that the spatial variation of the pair potential ∆(r) is the same in each layer. Then, the pair potential is expressed aŝ ∆(r) = ∆(r)iσ y ⊗ D, where D is the N × N diagonal matrix in the space composed of the layer (band) degree of freedom. In the case of N = 2, D = diag(1, s) with s = 1 for the BCS state and s = −1 for the PDW one.· denotes the 2N × 2N matrix in the spin and band space. A vortex line perpendicular to the two dimensional superconducting layer has the form, ∆(r) = |∆(r)|exp(iφ r ) in each layer.
In this paper, we consider the superconductor in which the paramagnetic pair breaking effect is dominant. In addition, we assume the system is in the type-II limit and neglect the vector potential. Then, we incorporate only the paramagnetic pair breaking effect into the quasiclassical theory through the Zeeman term. There is three factors which split the 4-fold degenerated FS. The spin degeneracy is lifted by both the SOC and the Zeeman field µ B H, and the band degeneracy is lifted by the interlayer hopping t ⊥ . We assume that the strength of these three factors is sufficiently weak and incorporate them into the quasiclassical theory as perturbations. In this situation, the split of the FS is small and the difference of the Fermi velocity between the weakly split FS is negligibly small. So we can put v F = v Fk . Then, we can expand the quasiclassical approximation described in Ref. [11] into the multilayered system.
The quasiclassical Green's function depends on the directionk = (cos φ k , sin φ k ) with the azimuthal angle φ k and is written as a 4N × 4N matrix in the particle-hole space: where ω n is the Matsubara frequency. The quasiclassical Green's function obeys the following Eilenberger equation with the spin quantization axis parallel to the z axis [12].
Here, I N ×N is the unit matrix, T ⊥ = offdiag(1, 1) and S d = diag(1, −1). offdiag(·, ·) indicates the 2 × 2 matrix which has the offdiagonal component only.Ĥ * SO (−k) = −g(k) · σ tr ⊗ S d , which is obtained from the relation g(−k) = −g(k). µ B is the Bohr magneton. We use the unit in which = k B = 1. For the numerical calculation, we employ the Riccati formalism. Regardinǧ K(k) as the self energy, one can use the same manner described in Ref. [13].
Starting from the pairing interaction in Ref. [14], the gap equation for the spin-singlet pair potential ∆(r) is given by [15] · · · k denotes the average on the circular arc with its radius 1. The coupling constant is determined from where T c0 is the superconducting transition temperature for α = t ⊥ = µ B H = 0 and n c (T ) = (ω c /πT − 1)/2. We fix the cut off frequency to ω c = 7πT c0 . The local density of states per spin and layer is given by where N F is the density of states per spin and layer at the Fermi level in the normal state.

Results and discussions
In the numerical results, we fix the temperature, the SOC and the interlayer hopping to T /T c0 = 0.1, α/T c0 = 2 and t ⊥ /T c0 = 1. In Fig. 1, we show the spatial variations of the pair potential amplitude around a vortex for (a) the BCS and (b) the PDW states. The horizontal axis is normalized by the coherence length ξ 0 = v F /T c0 . Each plot corresponds to the different magnetic field. The amplitude of the pair potential gets smaller with increasing the Zeeman field due to the paramagnetic pair breaking effect in both pairing states. One can notice that the pair potential amplitude of the PDW state is smaller than that of the BCS state. It is because the PDW state stabilizes in a high magnetic field under the influence of the paramagnetic depairing. Though the amplitude of the pair potential in the PDW state is smaller than that in the BCS state, the core size in the PDW state does not spread [compare the profiles for µ B H/T c0 = 1.5]. The core size in the PDW state for µ B H/T c0 = 1.5 and 2 is smaller than that in the BCS state for µ B H/T c0 = 1.5. Thus, the sudden decrease in the core radius with increasing the magnetic field is a signature of the BCS-PDW phase transition. The BCS state is completely destroyed for µ B H/T c0 2, whereas in the PDW state, the nonzero self-consistent solutions of the pair potential exist even for the high magnetic field µ B H/T c0 2. The profiles of the pair potential is not altered for α/T c0 = 50. In the uniform system without vortices, the stable pairing state prefers the larger bulk amplitude of the pair potential. In this paper, though there is a vortex in the system, we consider that the pair potential with the larger bulk amplitude is the more stable solution of the gap equation for a given magnetic field. The T -H phase diagram of bilayer system is displayed in Ref. [8], which is not altered qualitatively in the vortex state when the system belongs to type-II limit and has the large Maki parameter [8,16]. According to the field dependence of the pair potential amplitude in the bulk obtained by the quasiclassical theory, the non-zero self-consistent solution of the gap equation exists only in the PDW state for µ B H/T c0 1.8. This is consistent with the results calculated from the Bogoliubov-de Gennes equation [8].
In Fig. 2, we display the energy and spatial dependence of the local density of states (LDOS) around a vortex core. On the basis of the above consideration, the calculation of the LDOS is conducted for (a) µ B H/T c0 = 1.5 in the BCS state and (b) µ B H/T c0 = 2 in the PDW one. As shown in Fig. 2(a), the zero energy peak (ZEP) splits off due to the Zeeman field. The Zeeman split of the ZEP gets suppressed with increasing the SOC strength α for a fixed magnetic field, since the spin quantization axis is locked within the x-y plane by the Rashba SOC.
Contrastingly, though the PDW state emerges in a high magnetic field, the zero energy peak exists, which is quite different from the LDOS structure in the BCS state. If the zero energy quasiparticle states within a core are experimentally observed by scanning tunneling microscopy/spectroscopy (STM/STS), the STM/STS image of a vortex in the PDW state can be drastically different from that in the BCS state. To obtain the clue to the LDOS structure  in the PDW state, we start from the case with α/T c0 = 0 for simplicity. In this case, the four folding ZEP splits into four peaks. First, for µ B H/T c0 = 0, the band degeneracy is lifted by the interlayer hopping and the ZEP splits into the two peaks. This split of the ZEP is regarded as the level repulsion of the zero energy bound states. In the BCS state, the band degeneracy is not lifted only by the interlayer hopping. Then, the π phase difference between the layers is essential for the repulsion of the zero energy bound states. Second, each peak is split off due to the Zeeman field. Hence the four peaks appear due to the interlayer hopping and the Zeeman field. Further increasing the magnetic field, the two peaks in the low energy side shift to the lower energy side and those in the high energy side shift to the higer energy one. Eventually the two peaks in the low energy side get combined and a single zero energy peak appears.
In the real materials, for example, the CeCoIn 5 /YbCoIn 5 superlattice, the superconducting layer is CeCoIn 5 , which is experimentally identified as a d-wave superconductor [17]. The LDOS structure obtained in this study changes considerably, since the quasiparticle states around a vortex are sensitive to the pairing symmetry. The more realistic calculation assuming the d-wave pairing symmetry is the future work.

Summary
We have numerically investigated the electronic structure of a vortex core in bilayer Rashba superconductors by means of the self-consistent quasiclassical calculation. We found that the local density of states (LDOS) structure in the pair-density wave (PDW) state is quite different from that in the BCS state. The zero energy vortex bound state exsists in the PDW state, whereas it is absent in the BCS state due to the Zeeman effect. This characteristic LDOS structure can be observed by scanning tunneling spectroscopy at low temperature. Another feature of the PDW state is that the core size is small compared with that in the BCS state in the vicinity of the BCS-PDW phase transition. These features in a magnetic field can be a key to identify the exotic superconducting phase.