London Penetration Depth as a Test of Order Parameter Symmetry in Sodium Cobaltate Superconductors

Temperature dependence of the magnetic field penetration depth λ was calculated for water intercalated sodium cobaltate superconductor Na x CoO 2 · y H 2 O. Assuming that the system is in the chiral d+id–wave superconducting state, it was shown that the shifting of the excitation spectrum nodal points off the normal phase Fermi surface due to variation of the sodium content x changes the functional form of the temperature dependence of λ − 2 from exponential to linear in the low temperatures region. It is argued that this change in the functional form of T–dependence of the λ − 2 can serve as a proof for the chiral symmetry of the superconducting order parameter in the sodium cobaltate.


Introduction
Unusual properties of the superconducting state of the electron subsystem in sodium cobaltate Na x CoO 2 are due to the quasi-two-dimensional nature of the current carrier dynamics, as well as the triangular symmetry of the conducting layers of the crystal lattice formed by cobalt ions [1,2]. Cooper instability in these compounds occurs in the sodium content region 1 4 < x < 1 3 under water intercallation [3][4][5], and the maximum critical temperature T c does not exceed 5 K.
In the paper [18] it was shown that the mentioned experimental features can be described on the basis of a microscopic approach, assuming that the system is in the chiral d x 2 −y 2 +id xy -wave superconducting state, but the pairing occurs between the holes which are not on the nearest, but on the next nearest cobalt ions. Otherwise, the nodal (or rather Dirac) points are located only at the edges and in the center of the Brillouin zone and thus could not appear on the Fermi surface of the normal phase.
Complete suppression of Cooper pairing on the nearest sites (as it was supposed in [18]) is possible due to strong Coulomb repulsion V 1 of holes from the first coordination sphere. Since, however, in sodium cobaltates V 1 is not large, the study of superconductivity in these systems cannot be limited by taking into account the pairing between holes only from the second coordination sphere, it is also necessary to include pairings of the holes located on the nearest ions. In this case, as it was shown in [19], it is also possible for nodal points, at a certain value x, to be on the Fermi surface of normal phase, but for this it is necessary that the pairing amplitude for holes on the next nearest cobalt ions ∆ 2 exceeded the pairing amplitude for holes on the nearest ions ∆ 1 .
Vanishing only at certain points in the Brillouin zone is an important feature of the complex order parameter with d x 2 −y 2 +id xy -type of symmetry. For example, for d x 2 −y 2 -and d xy -pairings (separately forbidden by the triangular lattice symmetry) the whole lines of zeros are formed in the Brillouin zone and, therefore, the nodal points in the excitation spectrum of Bogolyubov quasiparticles exist on the Fermi surface at any current carriers concentration. For d+id pairing, nodal points on the Fermi surface can appear only at certain values of x = x c .
Given the above, we can assume that if the chiral order parameter is actually realized in sodium cobaltates, and the nodal points of this order parameter are on the Fermi surface at some value of x (from the interval [1 4, 1 3]), as stated in [18], then changing x should obviously shift these nodal points off the Fermi surface. This should necessarily lead for the gap to appear in the Bogolyubov quasi-particles excitation spectrum on the entire Fermi surface and, therefore, to the changes in the temperature dependence of various thermodynamic quantities in the superconducting state.
To verify the possibility of experimental observation of this effect in the sodium cobaltate Na x CoO 2 ⋅ yH 2 O, and, therefore, to prove that the order parameter in these materials has a d+id-wave symmetry, we study theoretically modifications of the temperature dependence of the magnetic field penetration depth when x changes in the range [1 4, 1 3] in which Cooper instability occurs.
The paper is organized as follows. In the second section, a microscopic model of sodium cobaltates is formulated. The third section discusses in detail the formula for calculating the temperature dependence of the London penetration depth. The fourth section presents the results of numerical calculations of λ −2 as a function of temperature. The final fifth section presents the main conclusion of the study. For convenience, cumbersome calculations were moved in Appendixes A and B.

Model
Microscopic three-band model for sodium cobaltates, accounting for all t 2g orbitals of cobalt ions and 2p orbitals of oxygen ions, was elaborated in [20] (LDA), [21] (DMFT), [22] (FRG). It was found that only one band, formed predominantly by 3d-orbitals (d xy , d yz , d zx ), intersects Fermi level, so that Fermi surface is represented by only one hole pocket around Γ-point in the Brillouin zone. At the same time, another LDA-calculations at x = 0.33 [23,24] besides the large Fermi surface around Γ-point of the Brillouin zone give also six small e ′ g -pockets in the vicinity of K-points. However, according to ARPES data [25][26][27] in the wide doping range x, there is only one wide Fermi surface in the vicinity of Γ-point. Thus we assume that electronic structure of the sodium cobaltates in the doping range of interest can be well described within the effective one-band or a 1g -band model [28]. The dispersion law of the hole-like quasiparticles of this band in the tight-binding approximation can be written as where t -is the tunneling integral between nearest-neighbor cabalt ions, and the wavenumbers k x and k y are given in the units of triangular lattice parameter a. The tunneling integral t ′ between next-nearest-neighbors was considered to be small here and hence was discarded. Though the energy spectrum (1) resembles the spectrum of non-interacting quasiparticles [29], the value of t is assumed to take into account renormalizations due to all interactions in the system. We chose t to be equal to 0.123 eV [28]. Dispersion of t k in the direction perpendicular to the triangular planes is not taken into account because spacing between the planes in water intercalated sodium cobaltates is fairly large, about 20 Angstroms [3].
Since, in the absence of sodium, the cobaltates are paramagnetic metal, and the compound NaCoO 2 is a band dielectric [1], it is reasonable to assume that the value of x = 0 corresponds to exactly one hole per cobalt ion in the a 1g -band (n h = 1), and at x = 1 the number of holes n h is zero. Therefore, there is a relationship between x and n h : x = 1 − n h . Since the interval for x variation we are interested in does not include the vicinity of the point x = 1 2, where the unusual dielectric phase arises, we can assume that the model of a 1g -band (1) describes correctly the spectral properties of the sodium cobaltate in the region of the phase diagram where superconductivity is observed.
Without going into the detail of the nature of Cooper instability, we just postulate the existence in the ensemble of holes of a superconducting pairing in the d+id-channel, which gives rise to one superconducting gap in the quasiparticle spectrum. According to this the amplitudes of the superconducting order parameter ∆ n (0) (n = 1, 2), at T = 0, will be considered below as parameters of the model. The amplitude ∆ 2 (0) will be defined by the critical temperature T c using the well-known, in the BCS theory [30], relation: 2∆ 2 (0) T c = 3.52 (see also [31,32]). The ratio of amplitudes ∆ 1 (0) ∆ 2 (0) will be determined so that the nodal points of the chiral order parameter were right on the Fermi surface at the optimal value of x (i.e. a value of x corresponding to the maximum critical temperature T c ). Since the Cooper instability in the sodium cobaltite occurs in the interval 1 4 < x < 1 3, and the maximum T c equal to 5 K corresponds to the middle of this interval, it seems natural to approximate the concentration dependence of the critical temperature in the given interval by parabola: The change in the amplitude of the order parameter with temperature will be modelled by the formula also following from BCS theory. When calculating the value of λ, we will also need the distance d between the CoO 2 -layers.
According to experimental data d = 19.6207

Equation for the London Penetration Depth
Evaluation of the magnetic field penetration depth λ is based on the London equation [33] describing in the local approximation the linear relationship of the superconducting current density j with a magnetic field vector-potential A. A comprehensive description of the methods for calculating the London penetration depth can be found in the textbooks [31,32] and original papers [34][35][36][37]. We shall follow the paper [37] where a simple and physically clear derivation of the formula for λ is given for arbitrary dispersion laws and gap functions of any symmetry. The formula reads as follows where When writing Expressions (4) and (5), the following notations were used: e-electron charge, c-speed of light, ̵ h-Planck's constant, T-temperature in energy units. The Bogoliubov quasiparticle spectrum E k is expressed in the usual way via the dispersion of holes t k and the gap function ∆ k where ξ k = t k − µ, µ -chemical potential.
Note that summation in Formula (4) runs over three components of the wave vector k: k x , k y and k z . Since the integrand function does not depend on k z , the triple summation can be reduced to a double one: where S is the area of the CoO 2 -layer, and d is the, defined above, distance between layers. Derivation of Formulas (4) and (5) could be found in Appendix A. The right part of Expression (4) can be conveniently written as the sum of two terms The first term gives the value of λ −2 at T = 0, and the second term describes temperature variation of λ −2 . In Formula (9): In the case of spin-singlet pairing two types of symmetry of the superconducting order parameter ∆ k are usually considered: s-and d+id-type. For triangular lattice we write down the order parameter of interest, with d+id-wave symmetry, in the complex form [18,38]: where ∆ n are the above-introduced pairing amplitudes between n-th neighbors, and the functions β ′ nk and β ′′ nk are defined by expressions Functions β ′ nk are related to d x 2 −y 2 -symmetry of the order parameter, and functions β ′′ nkd xy -symmetry. It should be remembered here that the second harmonics (n = 2) in (10) are necessary to be accounted for to shift the nodal points of the gap function inside the Brillouin zone on the Fermi contour. In this regard, it is worth clarifying the seemingly contradicting conclusion of the paper [22] on the absence of higher harmonic contributions to the d+id-wave pairing. The fact is that the paper [22] does not take into account both the implied by us intersite interactions, leading to the d+id-pairing, and the intersite Coulomb repulsion, which leads to the suppression of the first harmonic amplitude ∆ 1 as compared to the second harmonic amplitude ∆ 2 (see formula (10)). As a consequence, the authors of [22] did not obtain the nodal points of the gap function on the Fermi contour for d+id-pairing.
In the paper [37], Formula (4) was used to analyze the temperature dependence of the magnetic field penetration depth in cuprate HTSCs. The current carriers spectrum t k was constructed as a superposition of square lattice invariants and the tunneling integrals were chosen so as to fit the dispersion curves to the result of ARPES-experiments. The order parameter ∆ k for the cuprates was written as a difference of cosines: cos k x − cos k y . In our case, related to the triangular lattice, quasimomentum dependence of the function t k is defined by Formula (1) and the order-parameter-by Equation (10).

Temperature Dependence of the London Penetration Depth for the Chiral Order Parameter
Given the results of the paper [19] we will choose the ratio between the amplitudes of the ∆ 1 and ∆ 2 so that the Dirac points of the spectrum E k were strictly on the Fermi contour at the value x = 0.29, corresponding to the middle of the interval 1 4 < x < 1 3 and, respectively, to the maximum value of T c = 5 K. For the quasiparticle spectrum (1) and the gap function (10), this condition is satisfied at ∆ 1 ∆ 2 = 0.23 (see Figure 1). It can be shown (see Appendix B) that in this case the low-temperature corrections to the function λ −2 (T) are linear in T, i.e., δλ −2 (T) ∼ −T at T → 0. When the value of x changes, the nodal points shift and, as a result, the superconducting gap becomes different from zero on the entire Fermi contour. In this case the low-temperature part of the function λ −2 (T) should become exponential.
This means that if the order parameter in the sodium cobaltates has a chiral symmetry of d+id-type, as it was assumed in the paper [18], then with x variation we should also expect changes in the functional form of temperature dependence of λ −2 at T → 0: for x = 0.29 the function λ −2 (T) should be linear, but for values x other than 0.29 -exponential.
Note that for pure d x 2 −y 2 (or d xy ) pairing, as well as for s-pairing, the function λ −2 (T) should not change its analytical form in the low-temperature limit. In the first case, the nodal points of the gap function on the Fermi surface exist at any value of x, which means the δλ −2 (T) for all x is linear in T, and in the second case, there is always a gap on the entire Fermi surface and, therefore, for any x the T-dependence of δλ −2 must be exponential.
An important practical question that arises in this regard is the possibility to observe the modification of the functional form of δλ −2 (T) with changing the value of x in the narrow range [1 4, 1 3], in which sodium cobaltate manifests Cooper instability. Figure 2 shows the temperature dependence of λ −2 , calculated at the concentration x = 0.29 using the formula (4). The value of λ at zero temperature turns out to be 1.75 ⋅ 10 −5 cm. The analysis of the presented curves (see also the insert in Figure 2) clearly indicates the change in the analytical form of the temperature dependence of λ −2 (T) in the low-temperature region. It is seen that the solid line has a finite slope relative to the horizontal axis at T = 0, which indicates a linear dependence of the function λ −2 (T) at low temperatures and thus reveal the nodal points in the excitation spectrum on the Fermi surface. The dashed and dash-dotted lines corresponding to the x values lying on the boundaries of the superconducting region on the phase diagram are almost parallel to the horizontal axis in the T → 0 limit. This means that for these values of x in the Bogolyubov excitation spectrum, there is a gap on the entire Fermi surface.

Conclusions
In conclusion, we note that demonstrated possibility to observe the change in the analytical form of the temperature dependence of the inversed square of the London penetration depth in sodium cobaltate Na x CoO 2 ⋅ yH 2 O with sodium doping x change may provide an evidence for the chiral order parameter with d x 2 −y 2 +id xy -wave symmetry to be realized in this compound. For comparison we remind that in cuprate HTSC the Cooper instability occurs in the d x 2 −y 2 -channel. In this case, there are always nodal points in the spectrum of Bogolyubov excitations on the Fermi surface, and therefore the linear dependence of λ −2 (T) in the low temperature region is preserved at any doping. In paper [39], the linearity of the function λ −2 (T) was proposed to be considered as a proof of d x 2 −y 2 -wave pairing in cuprates. Similarly, in s-wave superconductors the temperature dependence of λ −2 is always described by an exponential function.
The insufficient amount of experimental data on the London penetration depth in sodium cobaltates is due to weak motivation of such experiments, we believe. We are aware of just a few of the experiments [40,41] with not good enough accuracy of measurements. The goal of the paper is to suggest to experimentalists a practical way to identify the order parameter symmetry in sodium cobaltate superconductors by analysing the changes in temperature dependence of the inverse square of the London penetration depth with doping. Derivation of the formula for the magnetic field penetration depth in superconductors on the triangular lattice is carried out following the paper [37]. As is known [42], the account for the magnetic field in the conduction electron Hamiltonian in the tight binding approximation leads to the renormalization of the tunneling amplitude t f m : where R f and R m are the radius-vectors of the sites f and m between which the electron tunnels, A is the magnetic field vector-potential. In this case the kinetic energy operator for electrons is where the operator c + f σ (c f σ ) creates (annihilates) an electron on the site f with spin projection σ = ±1 2. In the regime of weak magnetic fields (for which only the Peierls substitution (A1) is valid), the vector-potential A changes little on the scale of interatomic distances. This fact allows us to expand the phase factor in Expression (A1) in powers of A where R f m = R f − R m , and the value of the vector A is taken at the point (R f + R m ) 2. In accordance with the expansion (A3), we represent the Hamiltonian (A2) as Corrections of the first order in A to the kinetic energy Hamiltonian in the quasimomentum representation can be written in the form [37] H (1) In this formula, N is the number of crystal lattice sites, and Fourier transforms of the vector-potential and the second quantization operators are determined by expressions: A(r) = 1 N ∑ q e iqr A q and c f σ = 1 √ N ∑ q e iqR f c qσ respectively. The quasiparticle spectrum t k for the triangular lattice under consideration is given by Formula (1). Comparing Expression (A5) for H (1) kin with the formula for the energy of charged particles in the vector-potential field: we find an expression for the paramagnetic part of the electric current density operator where V is the crystal volume. For obtaining the average value of the superconducting current density j p q using the Equation (A7), the thermodynamic average ⟨c + kσ c k+q,σ ⟩ should be calculated at least in the linear approximation in A, since in the zero approximation j p q = 0. To calculate the average ⟨c + kσ c k+q,σ ⟩ in the superconducting phase with the required accuracy, we add to the Hamiltonian of a superconductor the operator H (1) kin [31,37]: where E k = ξ 2 k + ∆ k 2 is the Bogoliubov quasiparticle spectrum, ξ k = t k − µ, γ + kσ (γ kσ ) are the creation (annihilation) operators of Bogoliubov quasiparticles with quasimomentum k and quantum number σ.
Next, we express the operator H kin via γ-operators, and notice that to find the penetration depth λ, it is sufficient to consider the long-wavelength limit for the vector-potential Fourier transform A q . Assuming q = 0 we find Taking into account Expression (A9), the formula for the Hamiltonian of a superconductor (A8) can be represented in the form N ∇ k t k is the quasiparticle spectrum subject to a weak magnetic field. If now in Formula (A7) forĵ p q we also transform to the Bogolyubov quasi-particle operators, then for the average value of the paramagnetic part of the current density in the long-wavelength limit we obtain Since the thermodynamic averaging in (A11) involves a density matrix with the Hamiltonian (A10), then ⟨γ + kσ γ kσ ⟩ = f (Ẽ k ), where f is the Fermi-Dirac distribution function. Finally given the smallness of A, we obtain the expression for the paramagnetic part of the current density: To calculate the diamagnetic part j d of the current density j = j p + j d , we write down the expression for H (2) kin in the long-wavelength limit using Formulas (A2) and (A3): Let us rewrite Expression (A13) in terms of γ-operators and average it over thermodynamic ensemble. Comparing the result with the current density definition (A6), after a few simple transformations, we find: where u k and v k are coefficients of the Bogolyubov u − vtransformation: Combining Expressions (A12) and (A14), we obtain the formula for the net superconducting current density It is seen that in the normal phase (at ∆ k = 0), the superconducting current density j q=0 , as it should be, turns into zero.
Formally, from Expression (A15) it follows that the superconducting current density j q=0 is proportional not to the vector-potential A q=0 , but to the superposition of the dispersion gradient ∇ k t k and the gap function gradient ∇ k ∆ k 2 . This fact gives reason to suspect the nondiagonality of the tensor connecting the current density j q=0 with the vector-potential A q=0 . However, if we consider each projection j α q=0 (α = x, y) of the supercurrent density in Expression (A15) separately and take into account the symmetry of the integrand, then each of these projections turns out to be proportional only to the corresponding projection of the vector-potential A α q=0 : Numerical calculations show that, for both s -and d + idwave order parameter symmetry, the coefficient between j α q=0 and A α q=0 in Formula (A16) is actually independent of the value of α. Choosing, for definiteness, the direction of the vector A q=0 along x-axis, and comparing the expression for the superconducting current density (A16) to the London Equation (3) we obtain Formula (4) for the magnetic field penetration depth λ.

Appendix B
The general idea of calculating the low-temperature corrections to the thermodynamic properties of superconductors is explained in the book [43]. To analyze the low-temperature behavior of the inverse London penetration depth squared in quasi-two-dimensional superconductors on a triangular lattice, consider the situation when the nodal points of the Bogolyubov quasiparticle spectrum are located right on the Fermi surface (see Figure 1). It is obvious that in the limit T → 0, the main contribution to the integral in Expression (9) for the function δλ −2 (T) arises from a small neighborhoods of six nodal points, marked in Figure 1 with bold circles. To be definite, consider the contribution to the integral (9) from a small neighborhood U C of the nodal point C located in the positive domain of the k x -axis and having coordinates (k C , 0). It is seen that in this point the gradients of the spectrum ξ k and the real part of the gap function ∆ k are directed along the k x -axis, and the gradient of the imaginary part of ∆ k -along the k y -axis. This allows one to expand these functions in the neighborhood of the point C as: ξ k = ξ C (k x − k C ), Re∆ k = ∆ ′ C (k x − k C ), Im∆ k = ∆ ′′ C k y , where ξ C , ∆ ′ C and ∆ ′′ C are the values of the spectrum, real and imaginary parts of the gap function in the point C respectively. Substituting these expansions in the formula (9) and introducing new variables: q x = k x − k C , q y = k y , we get: where E q = (ξ 2 C + (∆ ′ C ) 2 )q 2 x + (∆ ′′ C ) 2 q 2 y , and the letters U C next to the symbol of integral mean that the integration is over small neighborhood U C around the point C.
Further we note that, in the limit T → 0, the main contribution to the integral over q y in Expression (A17) originates on the interval q y < where x = Dq x T. Since in the integral over x, the main contribution is due to the neighborhood of zero, the limits of integration over x can be extended to infinity. Given that in this case ∫ dx (e x + 1) = 2 ln 2, we finally find: Obviously integration over neighborhoods of other nodal points in formula (9) will also lead to linear low-temperature corrections to the value of λ −2 .