Relation of Superconducting Pairing Symmetry and Non-Magnetic Impurity Effects in Vortex States

Yasuaki Sera1, Takahiro Ueda 1, Hiroto Adachi 1,2 and Masanori Ichioka 1,2,* 1 Department of Physics, Okayama University, Okayama 700-8530, Japan; yasuaki-sera@s.okayama-u.ac.jp (Y.S.); takahiro-ueda@s.okayama-u.ac.jp (T.U.); hiroto.adachi@okayama-u.ac.jp (H.A.) 2 Research Institute for Interdisciplinary Science, Okayama University, Okayama 700-8530, Japan * Correspondence: ichioka@okayama-u.ac.jp


Introduction
Superconductivity is caused by the long-range order of the pair potential, which is wave function of Cooper pair. In type II superconductors [1][2][3], magnetic fields penetrate into the superconductor as quantized flux lines, and form a flux line lattice under applied fields between upper and lower critical fields. Around a flux line with a vortex of supercurrent, the complex value pair potential has the phase winding 2π around the vortex center, and the amplitude of the pair potential vanishes at the vortex center. At the vortex core where the pair potential is suppressed, low-energy bound states called Caroli-de Gennes Matricon states appear [4][5][6]. The local electronic structure of the bound state at the vortex core was experimentally observed by the scanning tunneling microscope (STM) observation by Hess et al. [7][8][9]. The STM observations of the vortex core image have been performed to study exotic properties of superconductivity in many superconductors [10][11][12][13], including high-T c cuprate superconductors [14][15][16][17][18][19], iron-based superconductors [20][21][22][23][24], and topological superconductors [24,25].
In unconventional superconductors with pairing mechanism other than the conventional electron-phonon interaction, anisotropic pairing symmetries such as p-, d-, and f -wave symmetries are realized rather than the isotropic s-wave. Therefore, to understand the pairing mechanism of unconventional superconductivity, we have to identify the pairing symmetry including the node positions of the pairing function. For the purpose, we can obtain important feature of the anisotropic The purpose of this work is to theoretically study the non-magnetic impurity scattering effects on the spatial structure of vortex states. Especially, we focus the contribution of the sign change of the pairing function on the Fermi surface, comparing the results of p x -wave paring symmetry ϕ p x (θ) with sign change, and the anisotropic s-wave ϕ |p x | (θ) without sign change. We evaluate the spatial structure of the pair potential and the LDOS in the vortex state by quantitative calculation of selfconsistent Eilenberger theory in the vortex lattice state [28,29,[41][42][43][44] and study the dependence on the scattering rate 1/τ 0 of non-magnetic impurities both in the Born and the unitary limits. Around the twofold symmetric vortex core, the anisotropy ratio is easily evaluated comparing the core radius of the vortex core image along the x and the y directions. From the results of the vortex structure, we discuss contributions from the sign change of ϕ (θ). The study about the anisotropic s-wave pairing can be applied to the twofold symmetric vortex core such as in NbNi 3 [13] and FeSe [20][21][22], which have deep minimum in the pairing function.
This paper is organized as follows. After this introduction, we explain our formulation in Section 2. In Section 3, we see the impurity effects in the uniform states. In Section 4 we study impurity effects on the pair potential and the LDOS in the vortex states for the cases of ϕ p x (θ) and ϕ |p x | (θ). The last section is devoted to discussion and conclusion.

Eilenberger Theory with Non-Magnetic Impurity Scatterings
Our numerical calculations of the vortex state are based on the Eilenberger theory [28,29,[41][42][43][44]. Eilenberger equation is derived from Gor'kov equation in the superconducting state under the assumption ∆ 0 /E F 1, which is satisfied in many superconductors, for the superconducting gap ∆ 0 and Fermi energy E F . The Eilenberger theory is suitable to treat inhomogeneous superconducting states, such as surface, interface, and vortex states. There, pair potential of superconductivity shows spatial variation in the length scale of the coherence length, which is longer than the atomic scale. To include contributions of non-magnetic impurities, we consider self-energies of the impurity scattering, as explained below. In addition to the Born limit [41,42], we cover also the unitary limit of the impurity scatterings by the t-matrix approximation [45][46][47][48][49][50].
Eilenberger equation is given by iv · ∇ĝ + iω nτ3 −∆ −Σ,ĝ =0 (1) in Nambu matrix form [45,50].ω n = ω n + iv · A with Matsubara frequency ω n , and the vector potential A. v = v F /v F with Fermi velocity v F and v 2 F = v 2 F k . · · · k indicates the Fermi surface average. The quasi-classical Green's functionsĝ and the pair potential∆ are, respectively, In our formulation using the Eilenberger unit, length, temperature, and magnetic field are, respectively, measured in unit of ξ 0 , T c , and B 0 . Here, ξ 0 =hv F /2πk B T c and B 0 = φ 0 /2πξ 2 0 with the flux quantum φ 0 . T c is superconducting transition temperature in the clean limit at a zero magnetic field. The energy E, pair potential ∆ and ω n are in unit of πk B T c .

Methods of Numerical Calculations in the Vortex Lattice
We calculate the spatial structure of vortices in the vortex lattice state by selfconsistent Eilenberger theory [28,29,[41][42][43][44]51,52]. We determine the spatial structure of vortices to be selfconsistent with the quasi-classical Green's functions of electronic states, solving the Riccati differential equations. We introduce contributions of the impurity scatterings by self-energies from non-magnetic s-wave impurity scatterings [43][44][45][46][47][48][49][50]. This method of Eilenberger theory can be applied to all temperature and magnetic field range in the vortex state, and appropriately capture contributions of vortex core and inter-vortex interaction. In many previous researches, local electronic states in the vortex states were studied in the clean limit. The researches of the impurity effect in the vortex states were done in pairing symmetries with isotropic gap, such as isotropic s-wave and chiral p-wave. Although there were some studies of impurity effect in the vortex state in the d x 2 −y 2 -wave pairing symmetry [44,47], we need systematic study of the impurity effect about the impurity scattering rate dependence and of the identification of contributions from the sign-change in the pairing function of anisotropic superconductors with line nodes.
In our calculation, magnetic fields B are applied z axis direction, so that B = (0, 0, B). The vector potential is given by A(r) = 1 2 B × r + a(r) in the symmetric gauge. For simplicity, we consider a two-dimensional cylindrical Fermi surface with a Fermi velocity v F = (v x , v y , 0) ∝ (cos θ, sin θ, 0) at a Fermi wave number k = (k x , k y , k z ) = (k F cos θ, k F sin θ, k z ). We consider two cases of the pairing function, ϕ(θ) = ϕ p x (θ) and ϕ |p x | (θ).
As the first step in our selfconsistent calculations of the vortex states, we calculate the quasiclassical Green's functions in the Eilenberger theory, under given ∆, A, G, F, and F † . To obtained the solution a(ω n , k, r), we numerically integrate Equation (11) along the trajectory from r − vs 0 to r [26,51,52]. As demonstrated in Figure 8 of [52], choosing enough long length s 0 of these trajectories we obtain unique solution independent of s 0 and initial values at r − vs 0 . The solution b(ω n , k, r) is obtained in the similar way by integrating Equation (12) along the trajectory from r + vs 0 to r. From the solutions a and b, we obtain quasiclassical Green's functions g, f , and f † in Equation (10) at r.
As an initial state, we start our calculations from the Ablikosov lattice solution, and a(r) = 0 of the vector potential A(r) = 1 2 B × r + a(r) in the symmetric gauge. Unit vectors of the vortex lattice are given by r 1 = (a x , 0) and r 1 = (a x , ζa y ) with ζ = 1 2 and Ba x a y = φ 0 . We set a x /a y = √ 3/2 in the triangular vortex lattice. r 0 = (x 0 , y 0 ) = − 1 2 (r 1 + r 2 ), so that one of vortex centers is located at 1 2 (r 1 + r 2 ) + r 0 = 0. A unit cell in our calculation is set to be u 1 (r 1 − r 2 ) + u 2 r 2 with −0.5 ≤ u i < 0.5 (i = 1, 2).
Using the quasiclassical Green's functions g, f , and f † obtained from Equations (10)-(12), we estimate new values of ∆(r) and A(r) by the gap equation (14) and the current equation We set in Equation (14), so that the value of the cut-off ω cut does not seriously affect the numerical results.
In the present calculation, ω cut = 20k B T c . We set the Ginzburg-Landau (GL) parameter κ = 30 assuming typical type-II superconductors. For the large κ, since variation of internal fields due to a(r) is very small compared with the applied field as is shown in Ref. [44], contributions of small a(r) do not seriously affect on results of the pair potential and the LDOS around a vortex. The self-energies G, F and F † are obtained from Equations (5) and (8) using the quasiclassical Green's functions. We calculate ∆(r), a(r), and self-energies at mesh points r in a unit cell of the vortex lattice, and estimate them at arbitrary points by interpolation of the values on the mesh points. Once we obtain ∆(r) and a(r) in a unit cell, we can know ∆(r) and a(r) in other unit cells by the lattice translation R = mr 1 + nr 2 (m and n are integers) as ∆(r + R) = ∆(r)e iχ(r,R) , a(r + R) = a(r), (17) where in the symmetric gauge. First, we solve Equations (9)-(15) at T = 0.5T c , and obtain selfconsistent solutions of ∆(r), A(r), and quasi-classical Green's functions. We perform the calculations for some values of 1/τ 0 in the Born and the unitary limits, in addition to the clean limit 1/τ 0 = 0. Using the renewed pair potential, vector potential, and self-energies, we solve the Riccati Equations (11) and (12) again to obtain new quasi-classical Green's functions. We iterate these calculation steps, until a sufficiently self-consistent solution is obtained.
Next, we study local electronic states. Using the selfconsistently obtained ∆(r) and A(r) in the calculation of ω n , we solve Riccati Equations (11) and (12) with iω n → E + iη to obtain g(iω n → E + iη, k, r), f (iω n → E + iη, k, r), and f † (iω n → E + iη, k, r) as a function of real energy E. η is an infinitesimal constant. Then, we calculate self-energies in Equation (8). These calculations of Riccati equations and the self-energies are iterated, until a selfconsistent solution is obtained. Using the selfconsistently obtained g, the LDOS is given by N(E, r) = N(E, k, r) k with k-resolved LDOS

Impurity Effects in Uniform States
Before studying the vortex states, we evaluate the strength of the superconductivity to the non-magnetic impurity scattering in uniform states without vortices at a zero magnetic field. Figure 2a shows the 1/τ 0 -dependence of uniform ∆ at T/T c = 0.5 to see how the impurity scattering suppresses p x -wave and anisotropic s-wave superconductivities. For the p x -wave pairing symmetry, ∆ is slightly smaller in the unitary limit compared with the Born limit. In both limits ∆ vanishes at 1/τ 0 ∼ 0.19. On the other hand, by the Anderson theorem, ∆ is not changed by the non-magnetic impurity scattering in isotropic s-wave symmetry. Even in the anisotropic s-wave pairing symmetry ϕ |p x | (θ), ∆ shows very weak suppression as is seen in Figure 2a. Therefore, s-wave superconductors are strong to the non-magnetic impurity scattering. Figure 2b shows ∆ in wider range until higher 1/τ 0 . There we see that ∆ for ϕ |p x | (θ) survives until 1/τ 0~2 2. The drop near 1/τ 0~2 2 is related to ω cut . We see that ∆ is smaller for larger ω cut when 1/τ 0 is large. Our study is performed at 1/τ 0 ≤ 1, where the ω cut -dependence of ∆ is negligible. The Born and the unitary limits have almost the same 1/τ 0 dependence for ϕ |p x | (θ). The difference between two limits comes from L ≡ g 2 k + f k f † k in Equation (5). For example, in isotropic s-wave pairing, since L = 1 from the relation g 2 + f f † = 1 of the quasiclassical Green's functions in Equation (10), 1/τ = 1/τ 0 in both limits.
Impurity effect on the DOS spectrum N(E) in uniform states for the p x wave pairing ϕ p x (θ) is presented in Figure 3a for the Born limit and in Figure 3b for the unitary limit. N(E) for the corresponding anisotropic s-wave ϕ |p x | (θ) is shown in Figure 3c for the Born limit, and in Figure 3d for the unitary limit. The horizontal axis is E/∆ with ∆ for each 1/τ 0 . Both ϕ p x (θ) and ϕ |p x | (θ) have the same N(E) in the clean limit (not shown in the figure), as the DOS is determined only by the amplitude of the pairing function, |ϕ(θ)|. We see a peak of N(E) at the gap edge E/∆ = √ 2, as the maximum of |ϕ p x (θ)| is √ 2. The peak is smeared with increasing the impurity scattering rate 1/τ 0 both in the Born and unitary limits. The smearing is seen in all cases of Figure 3. In the presence of line nodes, it is know that N(E) ∝ E near E = 0 in the clean limit 1/τ 0 = 0. For the p x -wave symmetry in Figure 3a in the Born limit, N(E) = 0 even when 1/τ 0 becomes finite. With increasing 1/τ 0 , N(E) at low E ( = 0) is enhanced. In the unitary limit in Figure 3b, finite DOS appears at E = 0 [53,54]. Flat N(E) near E = 0 becomes larger with increasing 1/τ 0 . On the other hand, for anisotropic s-wave symmetry in Figure 3c,d, the Born and the unitary limits show the similar behaviors. Even in the unitary limit of the impurity scattering, the zero-energy states do not appear, as N(E = 0) = 0. Rather, at finite 1/τ 0 gap opens, since N(E) = 0 at finite E near E = 0. With increasing 1/τ 0 , gap amplitude becomes larger.
∆ in the unitary limit shows almost the same 1/τ 0 -dependence.  The reason for the finite gap is because the gap structure is determined by ∆ϕ + F as is seen in Eilenberger Equations (9). In uniform states, the quasiclassical anomalous Green's function f has the same symmetry as the pairing function ϕ(θ). By the Fermi surface average, the anomalous self-energy F = 0 for ϕ(θ) = ϕ p x (θ). However, for ϕ(θ) = ϕ |p x | (θ), F = 0 and finite gap are induced due to the finite F, which has isotropic s-wave symmetry on the Fermi surface. However, as the anomalous f and F show different behaviors at the vortex core from those in the uniform states, we need to study the impurity effect around the vortex core by detailed calculations of the spatial structure in the vortex state.

Impurity Effects in Vortex
States for p x -Wave Pairing ϕ p x (θ) and Anisotropic s-Wave ϕ |p x | (θ) Figure 4a presents calculated spatial structure of the pair potential's amplitude |∆(r)| within a unit cell of the vortex lattice in xy plane for the p x -wave pairing symmetry ϕ p x (θ) at a low field B = 0.01B 0 . The amplitude |∆(r)| decreases with approaching the vortex center. The vortex core shape shows twofold symmetric shape, reflecting the twofold symmetric pairing function ϕ p x (θ) on the Fermi surface. This vortex core shape indicates the relation ξ x > ξ y , where ξ x (ξ y ) is the coherence length relating to the vortex core radius in x (y) direction. This anisotropy of the coherence length is consistent to the estimate ξ is larger on the Fermi surface. Even when the impurity effect is stronger, vortex core shape keeps similar twofold symmetric anisotropy both in the Born and the unitary limits. The suppression by the impurity scattering is stronger in the unitary limit compared with the Born limit also around the vortex core. These behaviors are also seen from the comparison of the profile between x and y directions shown in the lower panels of Figure 4a. In each panel, |∆(r)| and the slope around the vortex core are smaller along the x axis, compared to those along the y axis. When we compare the Born (solid lines) and the unitary (dashed lines) limits, the slope of |∆(r)| at the vortex center is almost the same in the both limits. Approaching outside region of the vortex core, deviation of |∆(r)| between two limits becomes larger. There, |∆(r)| in the Born limit is larger compared with the unitary limit. Figure 4b presents |∆(r)| for anisotropic s-wave pairing ϕ |p x | (θ). When the impurity scattering is weak with small 1/τ 0 , the vortex core shape shows similar twofold symmetric shape as in Figure 4a. However, approaching the dirty case with larger 1/τ 0 , the vortex core becomes circular shape. This is because contributions of the self-energy F with isotropic s-wave symmetry becomes dominant in the effective pair potential ∆ϕ + F. Differences between the Born and the unitary limits are small in this anisotropic s-wave pairing ϕ |p x | (θ). These behaviors are also seen from the comparison of the profile between x and y directions shown in the lower panels of Figure 4b. In each panel, as lines for the Born limit overlap with those for the unitary limit, the x and y dependences of ∆(r) are same in both limits. Although |∆(r)| is smaller along the x-direction, the deviation between x and y directions becomes small with increasing 1/τ 0 to be circular vortex core shape.
From the fitting of the slope of ∆(r) at the vortex center, we define the core radius ξ x along the x-direction as ∆(x, y = 0) = ∆ max x/ξ x near x = 0. ∆ max is maximum of |∆(r)| within a unit cell of the vortex lattice. Similarly we define the core radius ξ y along the y direction as ∆(x = 0, y) = ∆ max y/ξ y near y = 0. We present 1/τ 0 dependence of ξ x and ξ y in Figure 5a. Although each ξ increases as a function of 1/τ 0 , and ξ x is longer than ξ y . For p x -wave pairing symmetry, ξ x and ξ y rapidly increase toward the critical point 1/τ 0 ∼ 0.19. For anisotropic s-wave pairing symmetry, ξ x and ξ y gradually increase in this range of 1/τ 0 , as the critical point 1/τ 0 ∼ 22 is far. The anisotropy ratio ξ y /ξ x is presented in Figure 5b as a function of 1/ξ 0 . For the p x -wave, since the ratio ξ y /ξ x ∼ 0.65 in the range 0 ≤ 1/τ 0 < 0.19 of the superconductivity, twofold symmetric vortex core shape is unchanged by the impurity scattering, as seen Figure 4a. On the other hand, for anisotropic s-wave ϕ |p x | (θ), the ratio ξ y /ξ x gradually increases towards 1 changing to circular vortex core shape seen in Figure 4b.
The zero-energy LDOS N(E = 0, r) around a vortex core is presented in Figure 6a for the p x -wave pairing symmetry. N(E = 0, r) has a peak at the vortex center, and extends toward the y-axis, which is the node direction of the pairing function ϕ p x (θ) on the Fermi surface. In the vortex lattice, N(E = 0, r) is slightly suppressed on the line connecting nearest neighbor vortices [28,29]. This suppression is seen in N(E = 0, r) on the y axis at 1/τ 0 = 0.01.    Figure 6a shows how the local electronic states around a vortex are changed by the non-magnetic impurity scattering effects in the Born and the unitary limits, respectively. The sharp structure of N(E = 0, r) around the vortex core is smeared with increasing 1/τ 0 , and the smearing effect is greater in the unitary limit than in the Born limit. With increasing the impurity effect, orientation of the twofold symmetric vortex core shape changes. At larger 1/τ 0 , N(E = 0, r) extends toward the x direction. This is similar anisotropy to that of |∆(r)| in Figure 4a. In the unitary limit, we see enhancement of N(E = 0, r) outside of vortex core by the impurity effect, which occurs even in uniform states as shown in Figure 3b. These behaviors are also seen in lower panels of Figure 4b. There, we see that N(E = 0, r) at larger 1/τ 0 is larger on the x-axis, compared with the y-axis direction, because twofold symmetric vortex core image extends toward the x-direction. In the unitary limit, peak height at the vortex center is lower, and N(E = 0, r) is higher outside of the vortex core, compared with the Born limit.
On the other hand, impurity effect on the zero-energy LDOS N(E = 0, r) is presented in Figure 6b for anisotropic s-wave pairing ϕ |p x | (θ). When the impurity effect is weak with small 1/τ 0 , N(E = 0, r) shows similar anisotropic core image to that in Figure 6a, as the gap anisotropy is same as |ϕ p x (θ)| = |ϕ |p x | (θ)|. With increasing 1/τ 0 of the impurity effect, peak structure of N(E = 0, r) around the vortex core is smeared, and the vortex core image of N(E = 0, r) transforms to circular shape. In this process, the extension of N(E = 0, r) toward the node direction along the y axis is weakened, because contributions of the isotropic self-energy F become dominant. This changes to circular core shape is consistent to that of the pair potential in Figure 4b. The impurity effect in the unitary limit is qualitatively same as in the Born limit for the anisotropic s-wave pairing symmetry. These behaviors are also seen in lower panels of Figure 6b. Although peak height at the vortex center is higher in the unitary limit, the difference between the Born and the unitary limits becomes smaller with increasing 1/τ 0 . Thus, N(E = 0, r) has almost the same spatial structure in the two limits at 1/τ 0 = 1. To see the 1/τ 0 dependence quantitatively, we present the peak height N(E = 0, r = 0) at the vortex center as a function of 1/τ 0 in Figure 7a. For the p x wave pairing ϕ p x (θ) in the range 0 ≤ 1/τ 0 ≤ 0.15, with increasing 1/τ 0 , the peak height is suppressed toward normal state value 1. It is slightly smaller in the unitary limit than in the Born limit. On the other hand, for the anisotropic s-wave pairing ϕ |p x | (θ) in the range 0 ≤ 1/τ 0 ≤ 1, the suppression effect of the peak height by the impurity scattering is greater in the Born limit than in the unitary limit. Although the peak height is rapidly suppressed as in the p x -wave pairing in the Born limit, it decreases slowly in the unitary limit. These behavior in the anisotropic s-wave is consistent to results of isotropic s-wave pairing [47].
To quantitatively estimate the twofold symmetric vortex core shape in the zero-energy LDOS, the vortex core radii ξ x and ξ y are, respectively, determined as a radius of half width by N(E = 0, x = ξ x , y) = 0.5N(E = 0, x = y = 0) on the x axis and N(E = 0, x = 0, y = ξ y ) = 0.5N(E = 0, x = y = 0) on the y axis. The 1/τ 0 dependences of ξ x and ξ y are presented in Figure 7b. There we see clear differences between the Born and the unitary limits, as the peak height is different between two limits as shown in Figure 7a. For p x -wave pairing symmetry, ξ x and ξ y rapidly increase toward the critical point 1/τ 0 ∼ 0.19. The increase is greater in the unitary limit because N(E = 0, r) largely appear outside of the vortex. For anisotropic s-wave pairing symmetry, ξ x and ξ y gradually increase in this range of 1/τ 0 . The increase is smaller in the unitary limit than in the Born limit.  Figure 7. (a) Peak height N(E = 0, r = 0) at the vortex center as a function of impurity scattering rate 1/τ 0 for ϕ p x (θ) in the range 0 ≤ 1/τ 0 ≤ 0.15 and ϕ |p x | (θ) in the range 0 ≤ 1/τ 0 ≤ 1. The solid (dashed) lines are for the Born (unitary) limit. (b) 1/τ 0 dependence of vortex core radius ξ x and ξ y , which are determined from the zero energy LDOS N(E = 0, r) around a vortex core as described in text, for ϕ p x (θ) and ϕ |p x | (θ). Solid (dashed) lines are for the Born (unitary) limit. (c) The anisotropy ratio ξ y /ξ x as a function of 1/τ 0 , from ξ x and ξ y in panel (b).
The anisotropy ratio ξ y /ξ x from N(E = 0, r) is presented in Figure 7c as a function of 1/ξ 0 . In the clean limit, ξ y /ξ x ∼ 1.7, which is different anisotropy from the order parameter. For the p x -wave ϕ p x (θ), the ratio ξ y /ξ x rapidly decreases as a function of 1/τ 0 to ξ y /ξ x ∼ 0.5, which is the same orientation of twofold symmetry as in ∆(r). On the other hand, for anisotropic s-wave ϕ |p x | (θ), the ratio decrease until ξ y /ξ x ∼ 1.0 to be circular vortex core image. The decrease of ξ y /ξ x is slower in the unitary limit, compared with the Born limit.
Last, we study the spatial structure of the LDOS N(E, r) at finite energy E = 0.01, 0.1, 0.2, and 0.5. It is presented in Figure 8 for the p x -wave pairing ϕ p x (θ), and in Figure 9 for the anisotropic s-wave pairing ϕ |p x | (θ). These figures focus the vortex core region. At low energy E = 0.01, spatial structure of the LDOS is almost the same as the zero-energy LDOS N(E = 0, r) at 1/τ 0 = 0.01 and 0.1 presented in Figure 6. At 1/τ 0 = 0.05 in the Born limit, N(E, r) shows similar spatial structure both for the p x -wave ϕ p x (θ) in Figure 8 and the anisotropic s-wave ϕ |p x | (θ) in Figure 9. There, twofold symmetric vortex core image of the LDOS at E = 0.01 extends along the y axis, which is node-direction of the pairing function ϕ p x (θ). At E = 0.1 and higher E, the vortex core shape is seen to be rotated 90 • , and it extends to the x-direction. At E = 0.1 and 0.2, N(E, r) has large intensity near the vortex core, but it is suppressed at the vortex center. At higher energy E = 0.5, N(E, r) is rather smaller in the vortex core region, and has peak on the y axis far from the vortex center. These behaviors are also seen in the upper panels of Figures 8b and 9b. With increasing E, the peak position is shifted from the vortex center r = 0 at E = 0.01 to finite r. The distance between the peak position and the vortex center is longer at higher E. At 1/τ 0 = 0.01, N(E, r) in the unitary limit (dashed lines) has almost the same spatial structure as in the Born limit (solid lines) in Figures 8 and 9.
The spatial pattern of the LDOS and the E-dependence at 1/τ 0 = 0.01 are qualitatively explained by picture of the quasiparticle trajectory [27,31], which is schematically presented in Figure 10 for the p x -wave pairing symmetry ϕ p x (θ). This is the same also for ϕ |p x | (θ). As schematically presented in Figure 10a, the k-resolved LDOS N(E, k, r) for k with the velocity v = (cos θ, sin θ, 0) is mainly distributed on the line with distance r ⊥ (E, θ) from the vortex center. Roughly speaking, r ⊥ is given by the relation E = |∆(r)ϕ p x (θ)|. If we assume |∆(r)| = ∆ max r/ξ for simplicity, we obtain r ⊥ = ξ(E/∆ max )/|ϕ p x (θ)|. When a quasiparticle is propagating toward an antinode direction θ = 0 or π, r ⊥ takes the minimum value ξE/∆ max . When propagating toward a node direction θ = ±π/2, r ⊥ → ∞. On a line with angle α from the x-axis as shown in Figure 10a, distance s between the main distribution of N(E, k, r) and the vortex center is given by The LDOS N(E, r) is given by the sum of N(E, k, r) distribution in the range 0 ≤ θ < 2π. The peak of the LDOS appears at the inner edge of the LDOS distribution, where distance s on a line of angle α takes the minimum s min (E, α) as a function of θ. For the p-wave pairing symmetry ϕ p x (θ), Figure 10b,c, s min on a line of angle α is plotted for 0 ≤ α < 2π, as quasiparticle trajectories. On the trajectories, first, in-coming quasiparticles is pointing to the direction θ = −π/2 or π/2, where ϕ p x (θ) = 0, far from the vortex center. On going around the vortex core, when the direction θ of the quasiparticle is increasing, the distance between the trajectories and the vortex center is narrower since |ϕ p x (θ)| becomes larger. When the direction of the trajectory is θ = 0 or π, the distance is the minimum at y = ±ξE/∆ max . Then, trajectories go away from the vortex center towards the direction of θ = π/2 or 3π/2, with increasing θ. For small E in Figure 10b, tails of the trajectories extend toward ±y directions. This pattern of the trajectories corresponds to twofold symmetric vortex core image extending to the y-direction in the spatial structure of N(E, r) at E = 0.01 in Figure 8. On the other hand, for larger E in Figure 10c, shape of trajectories within the crossing points of two trajectories gives twofold symmetric vortex core image extending to the x-direction in N(E, r) at larger E in Figure 8, as tails of the trajectories far from the vortex are smeared.    E=0.01 E=0.1 E=0.2 E=0.5 that the peak height at the quasiparticle trajectories becomes lower at 1/τ 0 = 0.1, compared with that of 1/τ 0 = 0.01. The large differences between the Born and the unitary limits are only seen near the vortex core.

Discussion and Conclusions
To clarify contributions from the sign-change of the pairing function ϕ(θ) in the non-magnetic impurity scattering effects on the vortex state, we studied the twofold symmetric vortex core structure of the pair potential |∆(r)| and the LDOS N(E, r) in anisotropic superconductors with p x -wave pairing symmetry ϕ p x (θ) and anisotropic s-wave pairing ϕ |p x | (θ). The LDOS is accessible by the STM observation. In the STM experiments, twofold symmetric vortex core was observed in NiBi 3 [13] and FeSe [20][21][22], which may be anisotropic s-wave superconductors. If the impurity effects are experimentally studied in these superconductors, the theoretical results in this work may be examined by the STM observation. If a chiral p-wave superconductor exists, the p x -wave pairing superconductivity may be realized under uniaxial pressure, as tried in Sr 2 RuO 4 [55]. We note that twofold symmetric vortex core shape may also come from the contributions of anisotropic Fermi surface, in addition to anisotropic pairing symmetry. The study about the contribution of anisotropic Fermi surface belongs to future study.
As anisotropic superconductors with higher symmetry of the pairing, d x 2 −y 2 -wave pairing is realized in high T c cuprate superconductors [56,57], and some of heavy fermion superconductors [58,59]. As d x 2 −y 2 -wave pairing has four vertical line node on the Fermi surface, vortex core becomes four-fold symmetric shape [26,27]. As the corresponding anisotropic s-wave superconductors, four-fold symmetric vortex core image was observed in YNi 2 B 2 C by STM [10,11,32]. Therefore, detailed studies about the impurity effect of the vortex states are expected also in these superconductors. The twofold or four-fold symmetric vortex core structure may be related to various physical quantities in the anisotropic superconductors under magnetic fields, such as vortex lattice configuration [60] and angle-resolved specific heat measurement [58,59]. Therefore, studies about the impurity effects on these physical quantities are also interesting.
In summary, non-magnetic impurity scattering effect on the vortex core states are studied by the Eilenberger theory in superconductors with p x -wave pairing symmetry, as well as the corresponding anisotropic s-wave symmetry. From the comparison of the pair potential and the local electronic states around a vortex, we examine the differences between anisotropic superconductors with and without sign-change of the pairing function, and find how twofold symmetric vortex core changes with increasing the impurity scattering rate both in the Born and the unitary limits. For example, we found that twofold symmetric vortex core image of zero-energy LDOS changes the orientation of the twofold symmetry with increasing the scattering rate when the sign change occurs in the pairing function. Without the sign change, the vortex core shape reduces to circular one with approaching dirty cases. These results clarify the contributions from the sign-change of the pairing function in anisotropic superconductors, which are helpful to identify the pairing symmetry in unconventional superconductors, in relation to the STM vortex core observation.