Revisiting electromagnetic response of superconductors in mean-field approximation

In the standard mean-field treatment of superconductors, the electron-electron interactions are assumed to be written in terms of local density operators. However, more general interactions, such as pair-hopping interactions, may exist or may be generated in a low-energy effective Hamiltonian. In this work, we study the effect of correlated hopping interactions toward the electromagnetic response of superconductors. When only the Hamiltonian after the mean-field approximation is provided, one cannot unambiguously determine its electromagnetic response whenever such interactions are allowed. This work demonstrates that such interactions induce additional terms in the current operator, leading to modifications in the Meissner weight and optical conductivities that deviate from conventional expectations. These results underscore the need for caution when incorporating gauge fields into the BdG Hamiltonian.

Introduction.-One of the most remarkable features of superconductors is the Meissner effect, which is the expulsion of an applied magnetic field from the bulk of the sample [1].The Bardeen-Cooper-Schrieffer (BCS) theory successfully explained the mechanism based on the mean-field approximation [1].Although this treatment apparently breaks the U (1) symmetry, the vertex correction restores the gauge invariance of the response kernel [2].
In the study of superconductors at the mean-field level, one often starts with the Bogoliubov-de Gennes (BdG) Hamiltonian without specifying the Hamiltonian before the mean-field approximation.The BdG Hamiltonian is not invariant under U (1) phase rotation and thus its coupling to the gauge field is not uniquely determined.This results in ambiguities in the electromagnetic response of superconductors described by the BdG Hamiltonian.
To see the problems in details, let us revisit the meanfield approximation from microscopic perspective.The mean field approximation is a common technique to handle many-body interactions.However, when employing the mean field approximation, certain ambiguities arise as shown in Fig. 1.One source of ambiguity pertains to how we define the order parameters.Depending on specific choice of order parameters, different mean-field Hamiltonians can be derived, even when applied to the same microscopic model [Fig.1(a)].In this paper, our primary focus is on the superconducting order parameter, so this particular issue is not the central concern of our study.Another form of ambiguity arises when we only know a mean-field Hamiltonian because many microscopic models can end up with the same mean-field Hamiltonian [Fig.1(b)].This, in turn, leads to ambiguity in the physical observables such as electromagnetic responses.Furthermore, in the case of superconductive case, there is an additional ambiguity introduced when dealing with the incorporation of gauge fields [Fig.1(c)].Here, H k describes band dispersions for the normal phase, ∆ k is the gap function that satisfies ∆ −k = −β∆ T k , and C is a constant.The parameters α, β and the Nambu spinor Ψk are given by α = 1, β = −1, Ψk := (ĉ k↑ , ĉ † −k↓ ) T for spinful electrons, and α = 1/2, β = 1, Ψk := (ĉ k , ĉ † −k ) T for spinless electrons.The particle-hole symmetry P satisfies P 2 = β1.Note that the Hamiltonian is not invariant under U (1) gauge transformation(ĉ kσ → e iθ ĉkσ ), because the offdiagonal terms change under this transformation.

Let us consider superconductors described by the BdG Hamiltonian
To examine the electromagnetic response of the superconductor, it is customary to introduce the gauge field A arXiv:2304.07432v2 [cond-mat.supr-con]26 Oct 2023 by replacing H BdG k with (see, for example, Refs.[3][4][5]) and define the paramagnetic current operator and the kinetic energy operator by derivatives with respect to A: However, the A-dependence of H BdG k (A) in Eq. ( 2) is puzzling in two ways: A does not entirely appear in the form of k + A, which normally follows by the minimum coupling, and ∆ k does not depend on A even though it may have a nontrivial k-dependence.
The purpose of this work is to revisit these points and clarify the subtlety behind the BdG Hamiltonian.We argue that the gap function ∆ k and the constant C may depend on A and contribute to the paramagnetic current operator and the kinetic operator whenever the electronelectron interactions are not solely given in terms of the electron density nxσ := ĉ † xσ ĉxσ .As a consequence of these contributions, we find that the Meissner weight and the optical conductivity get modified from the standard results.
Furthermore, the BdG Hamiltonian for A = 0 is not sufficient to determine the A dependence of ∆ k and C unless the Hamiltonian before the mean-field approximation is provided.This means that the electromagnetic response of superconductor described by a BdG Hamiltonian alone is not completely well-defined.We clarify these points through the discussion of two illuminating examples at zero temperature T = 0.
s-wave superconductor with pair-hopping.-Asthe simplest example, let us consider a BCS type superconductor for single-band spinful electrons on d-dimensional hypercubic lattice: where ĉkσ 's are annihilation operators of electrons satis- is the nearestneighbor hopping, µ (|µ| < 2td) is the chemical potential that might be seen as the ν = 0 component of the gauge field A ν , U 0 is the onsite density-density interaction, and J is nearest-neighbor pair-hopping interaction.This model for the U 0 = A i = 0 case is called the Penson-Kolb model [6,7] and its electromagnetic response is studied in Refs.[8,9].We assume the periodic boundary condition with the length L i in i-th direction.The continuum version of this model is included in Appendix A.
The Hamiltonian has the U(1) symmetry associated with the electron density nx := σ=↑,↓ nxσ , which is necessary to fix the A dependence.The local current operator for the link between x and x + e i is given by and satisfies the continuity equation i nx , The second term in the current operator originates from the pair-hopping interaction.The model also possesses the spin rotation symmetry and, when A i = 0, the time-reversal symmetry and the inversion symmetry.
The superconducting order can be characterized by In this work, ϕ is assumed to be position-independent at least when A = 0.However, ϕ may depend on x when A i ̸ = 0.In fact, the large gauge transformation As we are interested in the large L i limit, for the consistency of our assumption we will keep |A i | much smaller than 2π/L i and will set A i = 0 at the end of the calculation.
The Hamiltonian Ĥ(A) in Eq. ( 5) can be converted to the BdG form in Eq. ( 2) by the mean-field approximation followed by the Fourier transformation ĉxσ := k ĉkσ e ik•x .We find that H BdG k (A) is given by the band dispersion When A = 0, the self-consistent equation for ϕ at T = 0 reads as where is the excitation energy of Bogoliubov quasi-particle, and ∆ := −U tot ϕ is the gap function for A = 0.The self-consistent equation contains only U tot , not U 0 nor J separately.Hence, if only ĤBdG is given without Ĥ, one cannot judge if U tot comes from the density-density interaction U 0 or the pair-hopping interaction J, as shown in Fig. 1(a).
p + ip topological superconductor in two dimension.-Theabove discussions are not restricted to s-wave superconductors.As a more nontrivial example, let us discuss a single-band model of spinless electron where U 0 describes the density-density interaction and J is a correlated hopping term that favors the p + ip order.The model has both U (1) symmetry and the four-fold rotation symmetry.The band dispersion is still given by ξ k above.The current operator and its continuity equation are summarized in Appendix B. Suppose that ϕ i := ⟨ĉ x+ei ĉx ⟩ is nonzero and position independent.The constant C for this model is given by When J = 0, the relative phase between ϕ 1 and ϕ 2 are arbitrary, while J > 0 favors (ϕ 1 , ϕ 2 ) = (iϕ, ϕ) with nonzero ϕ.This form of the gap function corresponds to the topological superconductor with the half-quantized thermal Hall conductance [10].Hence, the correlated hopping term (or some equivalent interaction) is necessary for the desired nontrivial topology.Assuming this type of order, we find Eq.( 2) with We stress that A does not enter these quantities through the replacement of k with k + A, even though it is introduced by the minimal coupling in the U(1) symmetric Hamiltonian.
Electromagnetic response.-Aswe have seen, when the starting Hamiltonian with U(1) symmetry contains interactions not solely written in terms of density operators, the superconducting gap ∆ k and the constant C in the BdG Hamiltonian generally depend on the gauge field A. Such dependence gives rise to additional terms in the current operator and the kinetic energy operator.For example, the current operator with a finite momentum q is given by with γ i,k+q,k := v i,k+ q 2 σ 0 − 2Jϕ sin qi 2 iσ 2 , which contains the contribution from the pair-hopping interaction in the off-diagonal part, in addition to the standard band velocity term v i,k := ∂ ki ξ k in the diagonal part.On the other hand, the charge density operator with γ 0,k+q,k := σ 3 is not affected by J.In the reminder of this work we discuss the physical consequences of these additional terms using the spinful electron model in Eq. ( 5).Without loss of generality, we set ∆ to be real using the U(1) symmetry.Meissner Weight.-Letus consider the linear response kernel of the current operator toward the gauge field with a frequency ω = q 0 and a momentum q: Here and hereafter we write q = (q 0 , q) for short.According to the linear response theory, the response kernel is given by where is the retarded current correlation function in the meanfield approximation and is the diamagnetic contribution giving the Meissner weight.We find and M BdG µν = 0 if µ or ν is 0. This result also applies to the spinless model and hence generalizes the result in Refs.[8,9].In addition to the usual term associated with the electron density ⟨n k ⟩ = α(1 − ξ k /E k ) and the band curvature ∂ ki ∂ kj ξ k , the second term arises from the A dependence of ∆ k and C. If the standard form in Eq. ( 2) were assumed instead, this term would be missed.This Here, we used t = 0.6, µ = 0, ∆ = 1.We expand Re[σ11(ω, q1)] in the Taylor series of q 2 1 and here we show the coefficient of the q 2 1 term.Red, green, and blue markers represent J = Utot, J = 0.5Utot, and J = 0. (b) The ratio of the difference of conductivity between J = Utot and J = 0 as a function of t at ω = 2.1.
effect might be measured through the penetration depth λ of an external magnetic field.
Response kernel with vertex correction.-As is wellknown, the kernel K BdG µν (q) in Eq. ( 15) in the mean-field treatment does not respect the gauge invariance.That is, the induced current d ν=0 K BdG µν (q)A ν (q) is not invariant under the gauge transformation A 0 (q) → A 0 (q) − ωχ(q) and A i (q) → A i (q) + 2 sin qi 2 χ(q).To obtain the gaugeinvariant optical conductivity σ ij (q) towards the electric field E j (q) = −iωA j (q), let us take into account the vertex correction following the steps in Refs.[1,2].
First, we define the vertex function Γ µ by and Here, is the time-ordered Green function.The continuity equation for the current operator [1,2] implies the generalized Ward identity (GWI): For reader's convenience, we review the derivation in Refs.[1,2] in Appendix A. The vertex function can be obtained by solving the Bethe-Salpeter equation i,k+q,k := v i,k+ q 2 σ 0 and v i,k := ∂ ki ξ k .
Once Γ µ is obtained, the time-ordered current correlation function P µν (q) can be expressed as Then the retarded correlation function R µν (q) is given by ReR µν (q) = ReP µν (q) and ImR µν (q) = sgn ω ImP µν (q) (see Appendix C).The optical conductivity σ ij (q) including the vertex correction is then given by for t = 1.5, µ = 0, ∆ = 1 (Utot is fixed by the self-consistent equation) in two dimensional system.We expand Re[σ11(ω, q1)] in the Taylor series of q 2 1 and here we show the coefficient of the q 2 1 term.Red squares, green triangle, and blue circles are for J = Utot, J = 0.5Utot, and J = 0, respectively.The insets in panels (a)-(b) are the log-log plots of the absolute value |Re[σ11(ω, q1)/q 2 1 ]|.Gray lines are obtained by fitting to determine the power of the decay Re[σ11(ω, q1)/q 2 1 ] ∝ ω −n .
Using the GWI ( 22) and the self-consistent equation ( 8), we have implying that the gauge invariance is restored: 2 and 3 show our numerical results for the optical conductivity σ ij (q) in Eq. (26) for q 2 = • • • = q d = 0.In single-band models with inversion and the time-reversal symmetry, the optical conductivity vanishes when q = 0 [3].Hence, here we assume spatially modulating field with q 1 ̸ = 0 and focus on the coefficient of q 2 1 in the Taylor series of Re[σ 11 (ω, q 1 )] with respect to q 1 .Here, we fix ∆ = 1, and U tot is fixed by the selfconsistent equation.
We found that nonzero J have a significant impact on the conductivity when the normal velocity v i,k is small, indicating that the band width of the normal state or t is small compared to J.In Fig. 2(a) we present a comparison of the results for three cases: J = 0, J = 0.5U tot , and J = U tot , while keeping U tot , t, and ∆ constant.Additionally, we examine the t-dependence of the ratio of the conductivity difference between J = U tot and J = 0 at ω = 2.1 in Fig. 2(b).The results show that when the normal state band is relatively flat in comparison to ∆, the nonzero J have a pronounced impact on conductivity.Recently, there has been a growing interest in models featuring flat or quasi-flat band superconductors [11][12][13].This finding underscores the need for caution when applying standard practices to address electromagnetic re-sponse in such models.
Figure 3 shows the results for the J = 0.5U tot case (red squares) and the J = 0 case (blue circles) with the large band width of the normal state compared to J in two and three dimensions.We found that nonzero J significantly affect the bare conductivities [(a),(d)] and the vertex corrections [(b),(e)] separately, although such differences are mostly suppressed in the sum σ ij (q) := σ In particular, the contribution from J usually decay with smaller power of ω as shown in the insets.However, the cancellation is not perfect and there are still finite differences originating from nonzero J.
Conclusion.-In this work, we revisited the electromagnetic response of superconductors.In general, the correspondence of the microscopic models and the BdG Hamiltonians after the mean-field approximation is "many to one": in the case of the spinful electron model in Eq. ( 5), as far as the renormalized parameter U tot = U 0 + Jd is fixed, models with different choices of on-site Coulomb interaction U 0 and the pair-hopping interaction J lead to the same BdG Hamiltonian.However, we found that the Meissner weight and optical conductivities are sensitive to the specific value of J, as shown in Eq. (18), Figs. 2 and 3.This means that the response toward U(1) gauge field has ambiguities unless the microscopic model with U(1) symmetry is provided.Our results call for caution in the standard practice in introducing the gauge-field A to the BdG Hamiltonian as in Eq. (2).

BdG Hamiltonian
The BdG Hamiltonian can be diagonalized as where γk↑ γ † −k↓ It follows that γ † kσ (t) = γ † kσ e iE k t in the Heisenberg picture.The annihilation operators of electrons can be expressed in terms of operators for Bogoloiubov quasi-particles: From these expressions, we find In the derivation, we used )

FIG. 1 .
FIG.1.Schematic illustrations of (a) one-to-many relationship between microscopic model and mean-field Hamiltonian, (b) many-to-one relationship between microscopic models and BdG Hamiltonian, and (c) electromagnetic response of superconductors.Given the U (1) symmetric Hamiltonian Ĥ, one can unambiguously introduce the gauge field A and derive the electromagnetic response with or without the mean-field approximation.In contrast, given the BdG Hamiltonian ĤBdG alone, one cannot uniquely determine ĤBdG (A) and cannot discuss its electromagnetic response without ambiguities.
FIG. 2. (a)The optical conductivities Re[σ11(ω, q1)] towards nonuniform field with the vertex correction for various J.Here, we used t = 0.6, µ = 0, ∆ = 1.We expand Re[σ11(ω, q1)] in the Taylor series of q2  1 and here we show the coefficient of the q 2 1 term.Red, green, and blue markers represent J = Utot, J = 0.5Utot, and J = 0. (b) The ratio of the difference of conductivity between J = Utot and J = 0 as a function of t at ω = 2.1.