AAxion electrodynamics in p + is superconductors Chao Xu and Wang Yang ∗ Department of Physics, University of California San Diego,9500 Gilman Drive, La Jolla, California 92093, USA Department of Physics and Astronomy and Stewart Blusson Quantum Matter Institute,University of British Columbia, Vancouver, B.C., Canada, V6T 1Z1
We perform a systematic study of axion electrodynamics in p + is superconductors. Unlike thesuperconducting Dirac/Weyl systems, the induced electric field does not enter into the axion action.Furthermore, in addition to the usual axion angle which is defined as the phase difference betweenthe superconducting phases on the two Fermi surfaces of different helicities, the axion field containsan additional sinusoidal term. Our work reveals the differences for axion electrodynamics betweenthe relativistic cases and the p + is superconductors. I. INTRODUCTION
Spin triplet superconductivity and paired superfluid-ity have a complex spin-orbit entangled structure in theCooper pair wavefunctions , leading to exotic behav-iors like topological properties . Another interest-ing type of Cooper pairing is the one with spontaneoustime reversal symmetry breaking, which arises when twoor more channels of pairing instabilities compete andcoexist . A mixture of triplet and singlet Cooperpairings which breaks time reversal symmetry has beentheoretically studied in different contexts, showing ex-otic properties including nontrivial bulk electromagneticand gravitational responses , quantized surface ther-mal Hall effects , chiral Majorana fermions propa-gating along the magnetic domain wall on the surface ,and high order topology . Recently, there has beenexperimental evidence of triplet pairing gap functionswith spontaneous time reversal symmetry breaking inreal materials .Axion as an elementary particle was proposed in thehigh energy context more than four decades ago ,which has been considered as a candidate for dark mat-ter and dark energy, thought its existence still remainsinconclusive. On the other hand, the dynamical axionfield has been proposed to exist in topological systemsas a condensed matter realization of axions . Thecoupling between the axion angle θ and the electromag-netic field is of the form (cid:82) d xθ (cid:126)E · (cid:126)B (up to an overallnumerical constant factor), resulting in various magne-toelectric effects, where (cid:126)E and (cid:126)B are the electric andmagnetic fields, respectively. In particular, the three-dimensional (3D) p + is superconductor has been con-sidered as a superconducting platform which hosts axionfield , where the triplet pairing component is invari-ant under spin-orbit coupled rotations analogous to thepairing of the He-B superfluid.In this work, we perform a systematic derivation ofthe coupling between the axion angle and the electro-magnetic field in p + is superconductors, including thecontributions from both the orbital and spin channels.The axion electrodynamics has been already studied insuperconducting 3D Dirac and Weyl systems with mixed- parity pairing which breaks time reversal symmetry .On the other hand, we find that there are several dif-ferences between the p + is superconductors and the su-perconducting Dirac/Weyl case. Firstly, the electromag-netic part of the action is not of the (cid:126)E · (cid:126)B form, but is ∇ ( φ − (cid:126) ∂ t Φ /e ) · (cid:126)B , where φ and Φ are the electric poten-tial and the superfluid phase, respectively. In particular,the induced electric field ∂ t (cid:126)A/c does not appear in theaction. The second difference is that the axion angle Θ ax in the p + is superconductors is not just θ ax , defined asthe phase difference between the superconducting phaseson the two Fermi surfaces of different helicities. In addi-tion to θ ax , Θ ax also acquires a sinusoidal term sin( θ ax ).These differences arise from a lack of Lorentz symmetryin p + is superconductors. II. MODEL HAMILTONIAN
We consider a superconducting 3D centrosymmetricelectronic system which exhibits a mixture of singlet andtriplet pairing symmetries. The band dispersion is ξ α ( (cid:126)k ) = (cid:126) m k − (cid:15) F , (1)in which (cid:15) F = (cid:126) m k f is the Fermi energy where k f isthe Fermi wavevector, and α = ↑ , ↓ is the spin index. Thepairing Hamiltonians ˆ P s ( (cid:126)k ) and ˆ P p ( (cid:126)k ) for the s -wave and He-B like p -wave pairing gap functions are defined asˆ P † s ( (cid:126)k ) = ( iσ ) αβ c † α ( (cid:126)k ) c † β ( − (cid:126)k ) , ˆ P † p ( (cid:126)k ) = 1 k f k j ( iσ j σ ) αβ c † α ( (cid:126)k ) c † β ( − (cid:126)k ) , (2)respectively, in which: α, β are spin indices; c † α ( (cid:126)k ) is theelectron creation operator with momentum (cid:126)k and spin α ; σ j ’s ( j = 1 , ,
3) are the three Pauli matrices in spinspace; and repeated indices imply summations. We notethat ˆ P p ( (cid:126)k ) is invariant under spin-orbit coupled SO(3)rotations, which has the same form as the pairing in the He-B superfluid . a r X i v : . [ c ond - m a t . s up r- c on ] S e p The pattern of the mixed-parity pairing gap functioncan be determined by a Ginzburg-Landau free energyanalysis. Keeping up to quartic terms and neglectingterms involving temporal and spatial derivatives, themost general form of the free energy invariant under bothtime reversal ( T ) and inversion ( P ) symmetries is givenby F = − α s ∆ ∗ s ∆ s − α p ∆ ∗ p ∆ p + β s | ∆ s | + β p | ∆ p | + γ | ∆ p | | ∆ s | + γ (∆ ∗ p ∆ ∗ p ∆ s ∆ s + c.c ) , (3)in which ∆ λ ( λ = s, p ) are the pairing gap functions inthe λ -channel. When the instabilities in the s - and p -wave channels coexist, both α s and α p are negative. Attree level, the coefficients β λ , γ j ( λ = s, p , j = 1 ,
2) areall determined by the electronic band structure, indepen-dent of the interactions. In particular, close to the super-conducting transition point, γ = ζ (3)8 π T c N F is genericallypositive , where ζ ( z ) ( z ∈ C ) is the Riemann zeta func-tion; T c is the superconducting transition temperature(where for simplicity, a degenerate transition tempera-ture for both s - and p -wave channels is assumed); and N F = mk f π (cid:126) is the density of states at the Fermi levelfor a single spin component. An important implicationof a positive γ is that a relative π/ s and ∆ p is energetically favorable as can bereadily seen from Eq. (3), leading to a superconductingpairing of the p ± is form. Notice that the p ± is pairingspontaneously breaks both time reversal and inversionsymmetries, but remains invariant under the combined PT -operation up to an overall gauge transformation .We note that in addition to the intrinsic p ± is supercon-ductors, the p ± is pairing symmetry can also be realizedvia proximity effects, as discussed in Ref. 24.In the remaining parts of this article, | ∆ λ | is denotedas ∆ λ for short, and the π/ p + is pairing Hamiltonian ∆ p ˆ P † p ( (cid:126)k ) + i ∆ s ˆ P † s . In theBogoliubov-de-Gennes (BdG) formalism, the mean fieldHamiltonian acquires the form H BdG = 12 (cid:88) (cid:126)k Ψ † ( (cid:126)k ) H ( (cid:126)k )Ψ( (cid:126)k ) , (4)in whichΨ † ( (cid:126)k ) = ( c †↑ ( (cid:126)k ) c †↓ ( (cid:126)k ) c ↑ ( − (cid:126)k ) c ↓ ( − (cid:126)k )) , (5)and the matrix kernel H ( (cid:126)k ) is H ( (cid:126)k ) = ( (cid:126) m k − (cid:15) F ) γ + ∆ p k f (cid:126)k · (cid:126)γ + ∆ s γ , (6)in which (cid:126)γ = ( γ , γ , γ ), and the matrices γ µ ( µ =0 , , , ,
4) are defined as γ = σ τ , γ = − σ τ , γ = − σ τ , γ = σ τ , γ = σ τ , (7) where τ i ( i = 1 , ,
3) are the Pauli matrices in the Nambuspace, and σ and τ denote the 2 × γ µ ( µ = 0 , , , ,
4) satisfy the anticommutation relations { γ µ , γ ν } = 2 δ µν . (8)Later we will also consider spatially and temporallyvarying pairing gap functions ∆ λ ( (cid:126)r, t ) ( λ = s, p ). Includ-ing the minimal coupling to electromagnetic potentials { A µ } µ =0 , , , , the Hamiltonian in real space becomes H BdG = 12 (cid:90) d r (cid:88) (cid:126)r Ψ † ( (cid:126)r ) ˆ T ρ γ ρ Ψ( (cid:126)r ) , (9)in which the summation is over ρ = 0 , , , ,
4; Ψ † ( (cid:126)r ) isthe Fourier transform of Ψ † ( (cid:126)k ); andˆ T = (cid:126) m ( − i ∇ + ec (cid:126)A ) − µ ( (cid:126)r ) , ˆ T j = 12 k f { ∆ p ( (cid:126)r, t ) , − i∂ j } γ j , ( j = 1 , , , ˆ T = ∆ s ( (cid:126)r, t ) . (10)In addition to Eq. (9), there is also a Zeeman term inthe presence of a magnetic field, i.e., H ZM = 12 gµ B (cid:90) d r (cid:88) (cid:126)r Ψ † ( (cid:126)r )[ − i (cid:15) ijk B i γ j γ k ]Ψ( (cid:126)r ) , (11)in which g is the Land´e factor, µ B is the Bohr magneton, (cid:126)B = c ∇ × (cid:126)A is the magnetic field, and the 1 / / III. TRANSVERSE SUPERCURRENT IN THEPRESENCE OF SPATIAL INHOMOGENEITIES
In this section, we perform a quick real space calcu-lation of the transverse supercurrent induced by staticelectric field and spatial inhomogeneity of ∆ s . The cor-responding term in the axion action S ax can then bedetermined from the relation j i = − c δS ax δA i ( i = 1 , , c is the light velocity. Later in Sec. IV, we willcarry out a systematic path integral calculation. In ad-dition to the supercurrent discussed in this section, thereis also bound current originating from the spin magneticmoment which is not related to electron transport, as willbe discussed in detail in Sec. VI.In BdG form, the matrix kernel of the operator of theelectric current density j i ( (cid:126)x ) ( i = 1 , ,
3) at position (cid:126)x isˆ j i ( (cid:126)x ) = − e (cid:126) m (cid:2) δ (ˆ (cid:126)r − (cid:126)x )( − i ∇ (cid:126)x ) + ( − i ∇ (cid:126)x ) δ (ˆ (cid:126)r − (cid:126)x ) (cid:3) σ τ , (12)in which ˆ (cid:126)r is the coordinate operator. The expectationvalue of ˆ j i ( (cid:126)x ) is (cid:68) ˆ j i ( (cid:126)x ) (cid:69) = Tr(ˆ j i ( (cid:126)x ) ˆ G ) , (13)in which the Green’s function ˆ G isˆ G = 1 − ∂ τ − H , (14)where τ is the imaginary time, and H is the matrix ker-nel of the Hamiltonian H BdG in the presence of a spa-tially varying electric potential φ ( (cid:126)r ) and s -wave pairinggap function ∆ s ( (cid:126)r ). We emphasize that the symbol “Tr”denotes the trace operation of an operator, which, in ad-dition to the trace of the matrix structure in the spinand Nambu spaces, also involves the integral over spatialcoordinates. In what follows, we use “tr” to indicate thetrace which is only taken over the 4 × G can be rewritten asˆ G = ( − ) 1 − ∂ τ + H ( − ∂ τ + H ) . (15)The Hamiltonian squared H can be separated as H = H + Q, (16)in which H = ˆ T ρ ˆ T ρ ,Q = (cid:88) ≤ ρ<σ ≤ [ ˆ T ρ , ˆ T σ ] γ ρ γ σ , (17)where ˆ T ρ is defined in Eq. (10) and the anticommutationrelations Eq. (8) is used. Straightforward calculationsshow thatˆ T ρ ˆ T ρ = [ − (cid:126) m ∇ − µ ( (cid:126)r )] + ∆ p k f ( −∇ ) + [∆ s ( (cid:126)r )] , [ ˆ T , ˆ T i ] = − i ∆ p k f ∂ i µ, [ ˆ T , ˆ T ] = − (cid:126) m ( ∇ ∆ s + 2 ∇ ∆ s · ∇ ) , [ ˆ T i , ˆ T ] = − i ∆ p k f ∂ i ∆ s , (18)in which µ ( (cid:126)r ) = (cid:15) F + eφ ( (cid:126)r ), where φ ( (cid:126)r ) is the electricpotential.Expanding (cid:68) ˆ j i ( (cid:126)x ) (cid:69) in powers of Q , we obtain (cid:68) ˆ j i ( (cid:126)x ) (cid:69) = ∞ (cid:88) n =0 ( − ) n +1 Tr (cid:2) ˆ j i ( (cid:126)x )( 1 − ∂ τ + H Q ) n × − ∂ τ + H ( − ∂ τ + H ) (cid:3) . (19)In what follows, we will only keep terms up to n = 2.The n = 0 term in Eq. (19) vanishes: The odd-in- ∂ τ term vanishes after Matsubara frequency summation;and the other terms contain a trace of a single γ -matrix,hence are also zero.The n = 1 term also vanishes. Removing the odd-in- ∂ τ term, the n = 1 term contains one ∆ H and one H under the trace operation. However, ∆ H and H are productsof two and one γ -matrices, respectively. Hence the resultvanishes since the trace of a product of three γ -matricesis zero.Finally we consider the n = 2 term (cid:68) ˆ j (2) i ( (cid:126)x ) (cid:69) = − Tr (cid:2) ˆ j i ( (cid:126)x ) 1 − ∂ τ + H Q − ∂ τ + H Q − ∂ τ + H H (cid:3) . (20)To lowest order in the gradient expansion , µ ( (cid:126)r ) and∆ s ( (cid:126)r ) can be taken as constants in H . The only wayto have a nonzero trace is a multiplication of all the five γ -matrices. Since the electric current operator ˆ j i ( (cid:126)x ) con-tains a − i∂ i , the H term must contribute ∆ p k f ( − i∂ i ) γ i sothat the trace is nonzero. Therefore, (cid:68) ˆ j (2) i ( (cid:126)x ) (cid:69) = 2( ∆ p k f ) Tr (cid:2) ˆ j i ( (cid:126)x ) − ∂ τ + H ∂ j µ − ∂ τ + H ∂ k ∆ s × − ∂ τ + H ( − i∂ i ) (cid:3) tr( γ γ j γ k γ γ i ) , (21)in which the trace of the product of five γ -matrices canbe easily evaluated usingtr( γ γ i γ j γ k γ ) = − (cid:15) ijk . (22)Using Eq. (12), (cid:68) ˆ j (2) i ( (cid:126)x ) (cid:69) can be evaluated as (cid:68) ˆ j (2) i ( (cid:126)x ) (cid:69) = − (cid:15) ijk e (cid:126) m ( ∆ p k f ) (cid:90) d(cid:126)yd(cid:126)z (cid:8) (cid:104) (cid:126)x | ( − i∂ i ) 1 − ∂ τ + H | (cid:126)y (cid:105) ∂ j µ ( (cid:126)y ) (cid:104) (cid:126)y | ( − i∂ i ) 1 − ∂ τ + H | (cid:126)z (cid:105)× ∂ k ∆ s ( (cid:126)z ) (cid:104) (cid:126)z | − ∂ τ + H ( − i∂ i ) | (cid:126)x (cid:105) + (cid:104) (cid:126)x | − ∂ τ + H | (cid:126)y (cid:105) ∂ j µ ( (cid:126)y ) (cid:104) (cid:126)y | − ∂ τ + H | (cid:126)z (cid:105) ∂ k ∆ s ( (cid:126)z ) × (cid:104) (cid:126)z | − ∂ τ + H ( − i∂ i ) | (cid:126)x (cid:105) (cid:9) . (23)In momentum space, Eq. (23) becomes (cid:68) ˆ j (2) i ( (cid:126)x ) (cid:69) = − (cid:15) ijk e (cid:126) m ( ∆ p k f ) × β (cid:88) ω n (cid:90) d(cid:126)yd(cid:126)z∂ j µ ( y ) ∂ k ∆ s ( z ) (cid:2) Π α =1 (cid:90) d(cid:126)k α (2 π ) (cid:3) × ( k i k i + k i ) (cid:2) Π α =1 ω n + [ H ( (cid:126)k α )] e ik α · ( (cid:126)r α − (cid:126)r α +1 ) (cid:3) , (24)in which β is the inverse temperature; ω n = (2 n + 1) π/β is the fermionic Matsubara frequency; (cid:126)r , , = (cid:126)x, (cid:126)y, (cid:126)z ;and (cid:126)r = (cid:126)r .To lowest order in the gradient expansion, we can set (cid:126)y = (cid:126)z in ∂ j µ ( (cid:126)y ) ∂ k ∆ s ( (cid:126)z ) within Eq. (24), then integrat-ing over (cid:126)z gives a momentum delta function δ (3) ( (cid:126)k − (cid:126)k ).Furthermore, if the higher order terms in the gradi-ent expansion are neglected, then (cid:126)y can be set as (cid:126)x in ∂ j µ ( (cid:126)y ) ∂ k ∆ s ( (cid:126)y ), and the integration over (cid:126)y gives δ (3) ( (cid:126)k − (cid:126)k ). As a result, we obtain j i ( (cid:126)x ) = (cid:80) n =0 (cid:68) ˆ j ( n ) i ( (cid:126)x ) (cid:69) as j i ( (cid:126)x ) = e D (∆ s , ∆ p ) (cid:15) ijk ∂ j φ ( (cid:126)x ) ∂ k ∆ s ( (cid:126)x ) , (25)in which ∂ j µ = − e∂ j φ is used, and the coefficient D is D (∆ s , ∆ p ) = 4 (cid:126) m ( ∆ p k f ) × β (cid:90) d k (2 π ) k / ω n + ξ ( (cid:126)k ) + ( ∆ p k f ) k + ∆ s ] , (26)where the replacement k i → k / β (cid:80) n can be replaced by (cid:82) dω π . In theweak pairing limit ∆ s , ∆ p (cid:28) (cid:15) F , k can be simply taken as k f in ( ∆ p k f k ) in the denominator, and (cid:82) dk = (cid:82) dξ/ ( (cid:126) v f )where v f = (cid:126) k f /m is the Fermi velocity. Eq. (26) canthen be evaluated, yielding D (∆ s , ∆ p ) = 16 π (cid:126) ∆ p (∆ p + ∆ s ) . (27)The expression of j i can be obtained from j i = − cδS ax /δA i , where S ax = − e c (cid:90) d xD (∆ s , ∆ p ) (cid:15) ijk A i ∂ j φ∂ k ∆ s . (28)Writing D (∆ s , ∆ p ) as a derivative of Θ ax , i.e., D (∆ s , ∆ p ) = − π (cid:126) ∂ Θ ax ∂ ∆ s , (29)and performing an integration by part, we arrive at S ax = − α π (cid:90) d x Θ ax (cid:15) ijk ∂ i φ∂ j A k , (30)where α = e (cid:126) c is the fine structure constant, and (cid:82) d x = (cid:82) d (cid:126)rdt . Imposing the boundary condition Θ ax = 0 fora pure s -wave superconductor (i.e., ∆ s (cid:29) ∆ p ), Θ ax isdetermined to beΘ ax (∆ s , ∆ p ) = 24 π (cid:90) ∞ ∆ s dxD ( x, ∆ p )= π − s ∆ p ∆ s + ∆ p − s ∆ p ) . (31)Therefore, defining θ ax = π − s ∆ p ) , (32)Θ ax can be written asΘ ax = θ ax − sin( θ ax ) . (33) Here we have two comments regarding the axion an-gle in Eq. (33) and the action in Eq. (30). First notethat the pairings on the two degenerate Fermi surfaces ofdifferent helicities (i.e., (cid:126)k · (cid:126)σ = ±
1) are | ∆ | e ± iθ ax where | ∆ | = (cid:113) ∆ s + ∆ p . However, the axion angle Θ ax is notjust θ ax ,which is the difference between the supercon-ducting phases on the two Fermi surfaces, but acquiresan additional sin( θ ax ) term. Secondly, the conventionalaxion term in a relativistic system is of the form ∼ (cid:126)E · (cid:126)B where (cid:126)E = −∇ φ − ∂ t (cid:126)A . On the other hand, Eq. (30)only contains the ∇ φ · (cid:126)B term, which is reasonable since inthis section we did not include a time-dependent vectorpotential. However, we will see in Sec. IV that in fact,the ∂ t (cid:126)A · (cid:126)B term is missing even in the general situation.Therefore, replacing φ by φ − (cid:126) e ∂ t Φ to keep gauge invari-ance, Eq. (30) is already the full axion action, where Φ isthe U(1)-breaking superconducting phase mode. As willbe discussed in Sec. IV, the difference between the p + is and the relativistic cases is the result of a lack of Lorentzinvariance in p + is superconductors. IV. PATH INTEGRAL FORMULATION
In this section, we formulate the systematic path in-tegral approach to the axion action in p + is supercon-ductors. We assume the four-fermion interaction to be ofthe form H int = − g s L (cid:88) (cid:126)k P † s ( (cid:126)k ) P s ( (cid:126)k ) − g p L (cid:88) (cid:126)k P † p ( (cid:126)k ) P p ( (cid:126)k ) , (34)in which L is the linear size of the system in space and g λ > λ = s, p ) are coupling constants, and P † λ ( (cid:126)k )’s( λ = s, p ) are defined in Eq. (2). After performinga Hubbard-Stratonovich transformation, the partitionfunction in the imaginary time formalism can be writ-ten as Z = (cid:90) D [ c † , c ] D [∆ ∗ s , ∆ s ] D [∆ ∗ p , ∆ p ] e −S , (35)in which S = S ∆ + S f , (36)where S ∆ = (cid:88) λ = s,p g λ (cid:90) dτ d r | ∆ λ | , (37)and the fermionic part S f is S f = 12 (cid:90) dτ d r Ψ † ( τ, (cid:126)r )( ∂ τ + H f )Ψ( τ, (cid:126)r ) , (38)in which Ψ † ( τ, (cid:126)r ) is a set of Grassmann numbers definedthrough the Fourier transform of Ψ † ( (cid:126)k ) in Eq. (5), and H f is given by H f = 12 m ( − i (cid:126) ∇ + ec (cid:126)A ) − (cid:15) F − ieφ + gµ B (cid:126)B · (cid:126)σ,H f = [ 12 k f e − i Φ { ∆ p , − i ∇} e − i Φ · (cid:126)σ + i ∆ s e − i (2Φ+Φ l ) ] iσ ,H f = H † f (1 , ,H f = − [ 12 m ( − i (cid:126) ∇ − ec (cid:126)A ) − (cid:15) F − ieφ + gµ B (cid:126)B · (cid:126)σ T ] , (39)in which H ijf ( i, j = 1 ,
2) is the ( i, j )-block of H f ; iφ is theelectric potential in imaginary time; { A, B } = AB + BA denotes the anticommutator of the operators A and B ;and Φ and Φ l are the superconducting phase mode andthe Leggett mode, respectively.The phase mode Φ can be absorbed into electromag-netic potentials by performing a gauge transformationthrough the following replacements: − ieφ → φ (cid:48) = − ieφ + ∂ τ Φ ,ec (cid:126)A → (cid:126)A (cid:48) = ec (cid:126)A + (cid:126) ∇ Φ . (40)Assuming a background ∆ λ ( λ = s, p ) and includingsmall fluctuations of the different modes, H f becomes H f = H f + ∆ H f , in which H f is simply Eq. 6, and∆ H f = ∆ H (1) A + ∆ H (2) A + ∆ H φ + ∆ H Z +∆ H p + ∆ H s + ∆ H l , (41)where ∆ H (1) A = (cid:126) m { (cid:126)A (cid:48) ( τ, (cid:126)r ) , − i ∇} σ τ , ∆ H (2) A = 12 m [ (cid:126)A (cid:48) ( τ, (cid:126)r )] γ , ∆ H φ = φ (cid:48) ( τ, (cid:126)r ) γ , ∆ H Z = − i gµ B (cid:15) ijk B i γ j γ k , ∆ H p = 12 k f { δ ∆ p ( τ, (cid:126)r ) , − i ∇} · (cid:126)γ, ∆ H s = δ ∆ s ( τ, (cid:126)r ) γ , ∆ H l = ∆ s δ Φ l ( τ, (cid:126)r ) σ τ . (42)We note that in momentum space, ∆ H (1) A ( (cid:126)q ) = (cid:126) m (2 (cid:126)k + (cid:126)q ) · (cid:126)A (cid:48) ( (cid:126)q ) and ∆ H p ( (cid:126)q ) = k f ( (cid:126)k + (cid:126)q/ · δ ∆ p ( (cid:126)q ).Since the fermionic quasiparticles are fully gapped, theaction for the collective bosonic degrees of freedom canbe obtained by integrating over the fermions, resulting in Trln( − ∂ τ − H f ). In what follows, we will only considerthe axion terms. They arise in the third order terms inthe Trln-expansion, i.e., S (3) f = −
16 Tr[( G ∆ H f ) ] , (43) in which G = ( ∂ τ + H ) − . Here we note that as dis-cussed in Sec. V and Sec. VI, in addition to the axionterms, there are other nonvanishing terms in S (3) f whichinvolve two spacetime derivatives, as a consequence ofa lack of Lorentz symmetry. We do not explicitly calcu-late these additional terms since the calculations are verycumbersome. A list of such terms based on a symmetryanalysis is included in Appendix A. V. ORBITAL CONTRIBUTION TO AXIONELECTRODYNAMICS
In this section, we calculate the orbital contributionto the axion action based on the path integral approach.There are three diagrams which contain two A µ ’s ( µ =0 , , , A is φ ) and one δ ∆ λ ( λ = s, p ), as shownin Fig. 1 (I, II, III). Diagram III – though not zero –does not contribute to the axion action, hence we neglect.We will only calculate the terms involving two spacetimederivatives in diagrams I, II in Fig. 1. A. Diagram I
This diagram potentially can contribute to (cid:82) d x∂ t (cid:126)A (cid:48) · (cid:126)B in the axion action. However, we demonstrate that infact, this contribution vanishes. λ = s One term contributing to the λ = s case is − Tr[ G ∆ H (1) A G ∆ H (1) A G ∆ H s ]. Including the combina-toric factor of 3, we obtain D I s = 12 1 β (cid:88) ω n (cid:90) d (cid:126)k (2 π ) tr (cid:2) G ( ω n , (cid:126)k ) (cid:126) m (2 (cid:126)k + (cid:126)q ) · (cid:126)A (cid:48) (Ω , (cid:126)q ) G ( ω n + Ω , (cid:126)k + (cid:126)q ) × (cid:126) m (2 (cid:126)k − (cid:126)q ) · (cid:126)A (cid:48) ( − Ω − Ω , − (cid:126)q − (cid:126)q ) × G ( ω n − Ω , (cid:126)k − (cid:126)q ) δ ∆ s (Ω , (cid:126)q ) γ (cid:3) , (44)in which the minus sign coming from the fermion loopcancels with the sign in Eq. (43), and G ( ω n , (cid:126)k ) = a χ ( ω n , (cid:126)k ) γ χ a χ ( ω n , (cid:126)k ) a χ ( ω n , (cid:126)k ) , (45)where the summation of χ is over χ = 0 , , , , , γ = σ τ ; and a = ξ (cid:126)k , a i = ∆ p k f k i ( i = 1 , , , a = ∆ s , a = iω n . (46) FIG. 1: Diagrams potentially contributing to axion electrodynamics where λ = s, p in ∆ λ . Plugging Eq. (45) into Eq. (44), the trace of the numer-ators of G givestr (cid:2) a χ ( ω n , (cid:126)k ) γ χ a χ ( ω n + Ω , (cid:126)k + (cid:126)q ) γ χ a χ ( ω n − Ω , (cid:126)k − (cid:126)q ) γ χ γ (cid:3) (47)We note that the O (Ω j ) ( j = 1 ,
2) term in Eq. (47) is8 ω n ∆ s (Ω − Ω ).To generate (cid:82) d x∂ t (cid:126)A · ∇ × (cid:126)A , we must considerterms in Eq. (44) which involve one Ω j and one (cid:126)q j ( j = 1 , O (Ω j )term in 1 / (cid:80) i =0 [ a i ( ω n + Ω , (cid:126)k + (cid:126)q )] (Ω = − Ω , Ω and (cid:126)q = − (cid:126)q , (cid:126)q ) is proportional to ω n Ω, similar as the O (Ω j )term in the trace in Eq. (47). Since the Matsubara sum-mation of terms involving odd powers of ω n vanishes, weconclude that there is no contribution to (cid:82) d x∂ t (cid:126)A ·∇× (cid:126)A from Diagram I for λ = s . λ = p The analysis for λ = p in Diagram I is exactly similarlyand the contribution again vanishes. B. Diagram II
This diagram potentially can contribute to (cid:82) d x ∇ φ (cid:48) · (cid:126)B in the axion action. We show that it gives exactly theaxion action derived in Sec. III. λ = s One term contributing to the λ = s case is − Tr[ G ∆ H (1) A G ∆ H φ G ∆ H s ]. Including the combina-toric factor of 6 and using Eq. (45), we obtain D II s = 1 β (cid:88) ω n (cid:90) d (cid:126)k (2 π ) (cid:126) m (2 (cid:126)k + (cid:126)q ) · (cid:126)A (cid:48) (Ω , (cid:126)q ) × φ (cid:48) ( − Ω − Ω , − (cid:126)q − (cid:126)q ) δ ∆ s (Ω , (cid:126)q ) × tr (cid:2) a χ ( ω n , (cid:126)k ) γ χ a χ ( ω n + Ω , (cid:126)k + (cid:126)q ) γ χ γ · a χ ( ω n − Ω , (cid:126)k − (cid:126)q ) γ χ γ (cid:3) , × Π ν =1 a χ ( ω n + Ω (cid:48) ν , (cid:126)k + (cid:126)q (cid:48) ν ) a χ ( ω n + Ω (cid:48) ν , (cid:126)k + (cid:126)q (cid:48) ν ) , (48)in which Ω (cid:48) = 0, Ω (cid:48) = Ω , Ω (cid:48) = − Ω , and (cid:126)q (cid:48) = 0, (cid:126)q (cid:48) = (cid:126)q , (cid:126)q (cid:48) = − (cid:126)q . Recall that we want a term ∼ (cid:15) ijk ∂ i φ∂ j A k in the action. Such term can be generated from the traceof a multiplication of the product of five γ -matrices asshown in Eq. (22). Then the (cid:15) ijk term within the tracein Eq. (48) can be evaluated astr[ ... ] = −
4( ∆ p k f ) (cid:15) ijk k i ( k j + q j )( k k − q k )= 4( ∆ p k f ) (cid:15) ijk k i q j q k . (49)Since Eq. (49) already contains a product of twowavevector q ’s, we can set (cid:126)q j , Ω j to be zero in the re-maining parts of Eq. (48). This gives D II s = D (∆ s , ∆ p ) (cid:15) ijk q j q k × A (cid:48) i (Ω , (cid:126)q ) δ ∆ s ( ω , (cid:126)q ) φ (cid:48) ( − Ω − Ω , (cid:126)q − (cid:126)q ) , (50)in which D (∆ s , ∆ p ) is given exactly by Eq. (26). Inte-gration by part and transforming to the real space, thecorresponding term in the action is (cid:90) d xD (∆ s , ∆ p ) (cid:126)A (cid:48) · ( ∇ φ (cid:48) × ∇ δ ∆ s ) . (51) λ = p When λ = p , D II p can be obtained from D II s by replac-ing δ ∆ s in Eq. (48) with δ ∆ p , and changing the traceto tr (cid:2) a χ ( ω n , (cid:126)k ) γ χ a χ ( ω n + Ω , (cid:126)k + (cid:126)q ) γ χ γ · a χ ( ω n − Ω , (cid:126)k − (cid:126)q ) γ χ k f ( (cid:126)k − (cid:126)q · (cid:126)γ (cid:3) . (52)The (cid:15) ijk term in Eq. (52) can be straightforwardly eval-uated as − s (∆ p ) k f (cid:15) ijk q i q j k k . (53)As a result, the corresponding term in the axion actionin real space is − D (cid:48) (∆ s , ∆ p ) (cid:126)A (cid:48) · ( ∇ φ (cid:48) × ∇ δ ∆ p ) , (54)in which D (cid:48) (∆ s , ∆ p ) = ∆ s ∆ p D (∆ s , ∆ p ) . (55) C. Diagram III
This diagram does not contribute to the axion action asexplained at the beginning of this section, though it doescontribute to non-axion terms as discussed in AppendixA.
D. Orbital contribution to the axion action
Combining Eq. (51,54) together and using ∇ δ ∆ λ = ∇ ∆ λ , we obtain S o ax = (cid:90) d x (cid:126)A (cid:48) · [ ∇ φ (cid:48) × ( D ∇ ∆ s − D (cid:48) ∇ ∆ p )] , (56)in which D ( D (cid:48) ) is D (∆ s , ∆ p ) ( D (∆ s , ∆ p )) for short.Plugging in the expression of D given in Eq. (27), wehave D ∇ ∆ s − D (cid:48) ∇ ∆ p = 16 π (cid:126) ∆ s ∆ p ) ] ∇ ( ∆ s ∆ p ) . (57) Using the integral (cid:90) ∞ ∆ s ∆ p dx (1 + x ) = 14 [ π − s ∆ p ) − s / ∆ p s / ∆ p ) ] , (58) S o ax becomes S o ax = − α π (cid:90) d x (cid:126)A (cid:48) · [ ∇ φ (cid:48) × ∇ Θ o ax ] , (59)where Θ o ax coincides exactly with the expression in Eq.(33). Integrating by parts, employing Eq. (40), andtransforming to the real time (i.e., iφ → φ ), we obtain S o ax = − α π (cid:90) d x Θ o ax ∇ ( φ − (cid:126) e ∂ t Φ) · (cid:126)B, (60)where ∇ × ∇ Φ = 0 is used. This is exactly what we haveobtained in Eq. (30).The differences between Eq. (60) and the axion actionin the relativistic case have been discussed by the end ofSec. III. The essential reason is a lack of Lorentz symme-try in the p + is case. Here we note that on a technicallevel, the band dispersion in the normal metal phase of p + is superconductor does not contain negative energystates. Therefore, unlike the relativistic case, there is nointer-band transition (from negative to positive energybands) in the p + is superconducting case, which leads tothe different behaviors between the two situations. VI. ZEEMAN CONTRIBUTION TO AXIONELECTRODYNAMICS
In this section, we calculate the Zeeman contributionto the axion action based on the path integral approach.
A. Diagram IV
This diagram potentially can contribute to (cid:82) d x ∇ φ (cid:48) · (cid:126)B in the axion action. We will derive its explicit expression. λ = s One term contributing to the λ = s case is − Tr[ G ∆ H Z G ∆ H φ G ∆ H s ]. Including the combina-toric factor of 6, we obtain D IV s = 1 β (cid:88) ω n (cid:90) d (cid:126)k (2 π ) tr (cid:2) G ( ω n , (cid:126)k ) φ (cid:48) (Ω , (cid:126)q ) γ G ( ω n + Ω , (cid:126)k + (cid:126)q ) × ( − ) 14 igµ B (cid:15) ijk γ j γ k B i ( − Ω − Ω , − (cid:126)q − (cid:126)q ) × G ( ω n − Ω , (cid:126)k − (cid:126)q ) δ ∆ s (Ω , (cid:126)q ) γ (cid:3) . (61)Only considering the axion term (cid:82) d x ∇ φ (cid:48) · (cid:126)B , we canset Ω = Ω = 0, and it is enough to expand D IV s upto linear order in (cid:126)q α ( α = 1 , α, β = 1 , i, j, k = x, y, z ; k i (cid:54) = k j (cid:54) = k k ) D IV s = E si iq i B i + E si iq i B i + O ( q αj q βk ) , (62)in which E αsi = − gµ B ∆ p k f β (cid:88) ω n (cid:90) d (cid:126)k (2 π ) M αsi , (63)where M si = − ω n + ξ (cid:126)k + ( k i − k j − k k ) ∆ p k f − ∆ s [ ω n + ξ (cid:126)k + ( ∆ p k f ) k + ∆ s ] ,M si = ω n − (cid:126) m ( k f − k + 4 k i ) ξ (cid:126)k + ∆ p k f ( k − k i ) − ∆ s [ ω n + ξ (cid:126)k + ( ∆ p k f ) k + ∆ s ] . (64)We note that: the q i B j ( i (cid:54) = j ) terms vanish in D IV s afterthe integration over the solid angle of (cid:126)k ; and k i can bereplaced by k / M αsi ≡ M αs , E αsi ≡ E αs , (65)where i = x, y, z .We will only keep the leading order terms in an expan-sion over ∆ λ /(cid:15) f ( λ = s, p ). Then in E s , all k can bereplaced by k f in the integrand M s , and (cid:82) d (cid:126)k/ (2 π ) canbe set as N f (cid:82) dξ , where N f = mk f π (cid:126) . On the other hand, M si contains a term − (cid:126) m ( k f + k ) ξ (cid:126)k , which is one orderless than the other terms in the numerator. Therefore, inthis case, (cid:82) d (cid:126)k/ (2 π ) should be set as N f (cid:82) (1 + ξ (cid:15) f ) dξ ,and M s → − ξ (cid:126)k (cid:15) f [ ω n + ξ (cid:126)k + ( ∆ p k f ) k + ∆ s ] + ω n − ξ (cid:126)k + ∆ p − ∆ s [ ω n + ξ (cid:126)k + ( ∆ p k f ) k + ∆ s ] + 4∆ p ξ (cid:126)k [ ω n + ξ (cid:126)k + ( ∆ p k f ) k + ∆ s ] , (66)where the first term can combine with N f (cid:82) ξ (cid:15) f dξ within (cid:82) d (cid:126)k/ (2 π ) giving a nonzero contribution, and the thirdterm in Eq. (66) comes from expanding the denominatorusing ∆ p k /k f = ∆ p (1 + ξ (cid:126)k /(cid:15) f ). Again at zero tem-perature, β (cid:80) ω n = (cid:82) dω π . Performing the integrations (cid:82) dω (cid:82) dξ , we obtain E s (∆ s , ∆ p ) = gµ B m π (cid:126) ∆ p (3∆ s + ∆ p )(∆ p + ∆ s ) ,E s (∆ s , ∆ p ) = gµ B m π (cid:126) ∆ p (∆ s − ∆ p )(∆ p + ∆ s ) . (67) Transforming back to the real space, the action be-comes S Z ax = (cid:90) d x (cid:2) E s ∇ δ ∆ s · (cid:126)Bφ (cid:48) + E s δ ∆ s (cid:126)B · ∇ φ (cid:48) (cid:3) . (68)Integrating by parts, dropping total derivative terms, andusing ∇ · (cid:126)B = 0, we obtain the axion action (cid:90) d x ( E s − E s ) φ (cid:48) ∇ δ ∆ s · (cid:126)B. (69) λ = p Similar as the λ = s case, the expression with a spa-tially varying δ ∆ p is D IV p = 1 β (cid:88) ω n (cid:90) d (cid:126)k (2 π ) tr (cid:2) G ( ω n , (cid:126)k ) φ (cid:48) (Ω , (cid:126)q ) γ G ( ω n + Ω , (cid:126)k + (cid:126)q ) × ( − ) 14 igµ B (cid:15) ijk γ j γ k B i ( − Ω − Ω , − (cid:126)q − (cid:126)q ) × G ( ω n − Ω , (cid:126)k − (cid:126)q ) 1 k f δ ∆ p (Ω , (cid:126)q )( (cid:126)k − (cid:126)q · (cid:126)γ (cid:3) . (70)Again setting Ω = Ω = 0, and keeping the linear in (cid:126)q α ( α = 1 , D IV p = E p iq i B i + E p iq i B i + O ( q αj q βk ) , (71)in which E αp = − gµ B ∆ s k f β (cid:88) ω n (cid:90) d (cid:126)k (2 π ) M αp , (72)where ( k i (cid:54) = k j (cid:54) = k k ) M p = ω n + ξ (cid:126)k + ( ∆ p k f ) ( k − k i ) + ∆ s [ ω n + ξ (cid:126)k + ( ∆ p k f ) k + ∆ s ] ,M p = 4[ − ( ∆ p k f ) ( k j + k k ) + (cid:126) m k i ξ (cid:126)k ][ ω n + ξ (cid:126)k + ( ∆ p k f ) k + ∆ s ] . (73)We note that again, the values of Eq. (73) do not dependon i = x, y, z . Calculations show that (up to leadingorder in ∆ λ /(cid:15) f , where λ = s, p ) E αp = − ∆ s ∆ p E αs . (74)Correspondingly, the contribution to the axion action is (cid:90) d x ( E p − E p ) φ (cid:48) ∇ δ ∆ p · (cid:126)B. (75) B. Diagram V
This diagram potentially can contribute to (cid:82) d x∂ t (cid:126)A (cid:48) · (cid:126)B in the axion action. However, in a way similar withthe discussions in Sec. V A, it can be seen that the axioncontribution from this diagram vanishes, since the inte-grand is odd with respect to the fermionic Matsubarafrequency. C. Zeeman contribution to the axion action
Combining Eqs. (69,75) and using Eq. (74), we obtain S Z ax = (cid:90) d xφ (cid:48) E Z (∆ s , ∆ p ) ∇ ( ∆ s ∆ p ) · (cid:126)B, (76)in which E Z = ∆ p ( E s − E s )= gµ B m π (cid:126) s / ∆ p ) ] . (77)Using (cid:90) ∆ s ∆ p dx (1 + x ) = 12 (cid:2) arctan(∆ s / ∆ p ) + ∆ s / ∆ p s / ∆ p ) (cid:3) , (78)and plugging in µ B = e (cid:126) mc , S Z ax becomes S Z ax = ge π (cid:126) c (cid:90) d xφ (cid:48) ∇ Θ Z ax · (cid:126)B, (79)in which α is again the fine structure constant, andΘ Z ax = 2 arctan( ∆ s ∆ p ) + 2∆ s ∆ p ∆ s + ∆ p . (80)Further performing an integration by part and transform-ing to real time, we obtain S Z ax = gα π (cid:90) d x Θ Z ax ∇ ( φ − (cid:126) e ∂ t Φ) · (cid:126)B, (81)in which g = 2 in vacuum, but can be somewhat arbitraryin solid state materials. The action S Z ax in Eq. (81) canalso contribute to transverse supercurrent, similar to itsorbital counterpart as discussed in Sec. III. However,we note that unlike the orbital case, the current arisingfrom Eq. (81) are bound currents, which are not relatedto wavepacket transport of electrons.Here we note that Θ Z ax = π − Θ o ax , where the expres-sion of Θ o ax is given in Eq. (33). This difference comesfrom the integral (cid:82) ∆ s / ∆ p dx in Eq. (78), unlike the onein Eq. (58) which is (cid:82) ∞ ∆ s / ∆ p dx . The choice of the lowerand upper bounds of the integration in Eq. (78) is con-sistent with Ref. 24, where the free energy describing themagnetoelectric effect in the spin channel is calculated in the vicinity of T c . It is shown in Ref. 24 that for a pure p -wave superconductor, the spin magnetoelectric effectvanishes. Therefore, in the current situation, Θ Z ax shouldbe zero for a pure p -wave superconductor in the absenceof ∆ s . VII. CONCLUSION
In conclusion, we have studied the coupling betweenthe axion angle and the electromagnetic field in p + is superconductors. Including both the orbital and Zeemancontributions, the axion action is derived as S ax = − α π (cid:90) d x Θ ax ∇ ( φ − (cid:126) e ∂ t Φ) · (cid:126)B, (82)in which α is the fine structure constant, and the axionangle is Θ ax = Θ o ax − Θ Z ax = − gπ + (1 + g )( θ ax − sin θ ax ) , (83)where g is the Land´e factor and θ ax = π − s / ∆ p ). The axion action in Eq. (82) has twodifferences compared with the relativistic case. The in-duced electric field does not appear, therefore S ax doesnot have Lorentz symmetry. In addition, the axion angleΘ ax also contains a sinusoidal term sin( θ ax ), not just θ ax .Our work reveals the crucial difference for axion electro-dynamics between the p + is superconductors and thesuperconducting Dirac/Weyl systems. Acknowledgments
CX is supported by The Office ofNaval Research under grant N00014-18-1-2722. WYis supported by the postdoctoral fellowship at StewartBlusson Quantum Matter Institute, University of BritishColumbia.
Appendix A: Symmetry allowed non-axion terms
In this appendix we examine all the symmetry al-lowed terms which contain one ∆ λ ( λ = s, p ), two A µ ’s( µ = t, x, y, z ), and two spacetime derivatives. Up to anoverall factor F (∆ s , ∆ s / ∆ p )∆ λ (where F is a functionwhich can be determined by calculating the correspond-ing diagram), the terms invariant under 3D rotations and PT -operation are ∂ t (cid:126)A (cid:48) · ∂ t (cid:126)A (cid:48) δ ∆ λ ∂ t (cid:126)A (cid:48) · (cid:126)A (cid:48) δ ∆ λ ∂ t φ (cid:48) ∂ t φ (cid:48) δ ∆ λ ∂ t φ (cid:48) φ (cid:48) δ ∆ λ (A1)( ∇ · (cid:126)A (cid:48) )( ∇ · (cid:126)A (cid:48) ) δ ∆ λ ∇ (cid:126)A (cid:48) · (cid:126)A (cid:48) δ ∆ λ ∇ φ (cid:48) · ∇ φ (cid:48) δ ∆ λ ∇ φ (cid:48) φ (cid:48) δ ∆ λ (A2)0 ∂ t (cid:126)A (cid:48) · ( ∇ × (cid:126)A (cid:48) ) δ ∆ λ (cid:126)A (cid:48) · ( ∇ × ∂ t (cid:126)A (cid:48) ) δ ∆ λ (A3)( ∇ · (cid:126)A (cid:48) ) ∂ t φ (cid:48) δ ∆ λ (A4)( ∇ × (cid:126)A (cid:48) ) · ( ∇ φ (cid:48) ) δ ∆ λ , (A5)in which λ = s, p . The terms in Eq. (A3) vanish asdiscussed in Sec. V A, and Eq. (A5) is the axion termwhich has been calculated and discussed in the main text.Here we make a comment on the order of the coeffi-cients of the non-axion terms. In superconductors, theleading term in the action for the phase mode is (cid:90) d x (cid:2) N f ( (cid:126) ∂ t Φ − eφ ) + n s m ( (cid:126) ∇ Φ + e (cid:126)Ac ) (cid:3) , (A6)in which n s is the superfluid density. Notice that inthe long wavelength limit, additional spacetime gradient terms are suppressed by factor of ( qξ c ) n ∼ ( (cid:126) v f q/ ∆) n ,where ξ c is the coherence length, and q can be either | (cid:126)q | or | Ω | /v f . For simplicity, consider the term ∂ t (cid:126)A · ∇ × (cid:126)A (although this term vanishes as discussed in Eq. (A3),it works as an illustration for the other non-axion terms,regarding the order of the coefficients). The coefficient C of this term should be on order of ∼ ( v f ∆ ) n s m . Us-ing n s ∼ k f , it is straightforward to obtain C ∼ ( (cid:15) f / ∆) .Other other hand, recall that the coefficients of the axionterms are of O [(∆ /(cid:15) f ) ]. Therefore, generically, the non-axion terms can be much larger, i.e., enhanced by a factorof ( (cid:15) f / ∆) compared with the axion terms. However, wenote that none of the non-axion terms can contribute tothe effects like transverse supercurrent as discussed inSec. III, and in fact, they do not exhibit magnetoelectriceffects. ∗ Electronic address: [email protected] A. J. Leggett, Rev. Mod. Phys. , 331 (1975). D. Vollhardt and P. W¨olfle,
The Superfluid Phases of He-lium 3 (Taylor and Francis, London, 1990). G. E. Volovik,
The Universe in a Helium Droplet (OxfordUniversity Press, New York, 2003). A. P. Mackenzie and Y. Maeno, Rev. Mod. Phys. 75, 657(2003). A. J. Leggett,
Quantum Liquids (Oxford University, NewYork, 2006). M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045(2010). X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057(2011). Y. Ando and L. Fu, Annu. Rev. Condens. Matter Phys. ,361 (2015). M. Sato and S. Fujimoto, J. Phys. Soc. Jpn. , 072001(2016). C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev.Mod. Phys. , 035005 (2016). M. Sato and Y. Ando, Rep. Prog. Phys. , 076501 (2017). G. Volovik, Phys. Lett. A , 277 (1988). G. Volovik and V. Yakovenko, J. Phys. Condens. Matter , 5263 (1989). D. A. Ivanov, Phys. Rev. Lett. , 268 (2001). N. B. Kopnin and M. M. Salomaa, Phys. Rev. B , 9667(1991). M. Stone and S.-B. Chung, Phys. Rev. B , 014505(2006). S. Tewari, S. Das Sarma, C. Nayak, C. Zhang, and P.Zoller, Phys. Rev. Lett. , 010506 (2007). C. Zhang, S. Tewari, R. M. Lutchyn, and S. Das Sarma,Phys. Rev. Lett. , 160401 (2008). L. Fu and C. L. Kane, Phys. Rev. Lett. , 096407 (2008). M. Cheng, K. Sun, V. Galitski, and S. Das Sarma, Phys.Rev. B , 024504 (2010). X.-L. Qi, T. L. Hughes, S. Raghu, and S.-C. Zhang, Phys.Rev. Lett. , 187001 (2009). Y. Wang and A. Chubukov, Phys. Rev. B , 035149(2014). Y. Wang and L. Fu, Phys. Rev. Lett. , 187003 (2017). W. Yang, C. Xu, and C. Wu, arXiv:1711.05241 (2017). W.-C. Lee, S.-C. Zhang, and C. Wu, Phys. Rev. Lett. ,217002 (2009). Y. Li, D. Wang, C. Wu, New Journal of Physics 15 (8),085002 (2013). R. Thomale, C. Platt, W. Hanke, and B. A. Bernevig,Phys. Rev. Lett. , 187003 (2011). C. Platt, R. Thomale, C. Honerkamp, S.-C. Zhang, andW. Hanke, Phys. Rev. B , 180502 (2012). M. Khodas and A. V. Chubukov, Phys. Rev. Lett. ,247003 (2012). R. M. Fernandes and A. J. Millis, Phys. Rev. Lett. ,127001 (2013). A. Hinojosa, R. M. Fernandes, and A. V. Chubukov, Phys.Rev. Lett. , 167001 (2014). V. Stanev and Z. Teˇsanov´ıc, Phys. Rev. B , 134522(2010). S.-Z. Lin and X. Hu, Phys. Rev. Lett. , 177005 (2012). M. Marciani, L. Fanfarillo, C. Castellani, and L. Benfatto,Phys. Rev. B , 214508 (2013). S. Maiti and A. V. Chubukov, Phys. Rev. B , 144511(2013). F. Ahn, I. Eremin, J. Knolle, V. B. Zabolotnyy, S. V.Borisenko, B. B¨uchner, and A. V. Chubukov, Phys. Rev.B , 144513 (2014). J. Garaud and E. Babaev, Phys. Rev. Lett. , 017003(2014). S. Maiti, M. Sigrist, and A. Chubukov, Phys. Rev. B ,161102 (2015). R. B. Laughlin, Phys. Rev. Lett. , 5188 (1998). T. Senthil, J. B. Marston, and M. P. A. Fisher, Phys. Rev. B , 4245 (1999). B. Horovitz and A. Golub, Phys. Rev. B , 214503 (2003). Y. Jiang, D.-X. Yao, E. W. Carlson, H.-D. Chen, and J.Hu, Phys. Rev. B , 235420 (2008). M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. B ,134521 (2010). A. M. Black-Schaffer, Phys. Rev. Lett. , 197001 (2012). R. Nandkishore, L. S. Levitov, and A. V. Chubukov, Nat.Phys. 8, (2012). W.-S. Wang, Y.-Y. Xiang, Q.-H. Wang, F. Wang, F. Yang,and D.-H. Lee, Phys. Rev. B , 035414 (2012). M. L. Kiesel, C. Platt, W. Hanke, and R. Thomale, Phys.Rev. Lett. , 097001 (2013). F. Liu, C.-C. Liu, K. Wu, F. Yang, and Y. Yao, Phys. Rev.Lett. , 066804 (2013). A. M. Black-Schaffer and C. Honerkamp, J. Phys.: Con-dens. Matter. C. Wu and J. Hirsch, Phys. Rev. B , 020508 (2010). A. J. Niemi and G. W. Semenoff, Phys. Rep. , 99(1986). S. Ryu, J. E. Moore, and A. W. W. Ludwig, Phys. Rev. B , 045104 (2012). X.-L. Qi, E. Witten, and S.-C. Zhang, Phys. Rev. B ,134519 (2013). Z. Wang, X. L. Qi, and S. C. Zhang, Phys. Rev. B , 014527 (2011). N. Read and D. Green, Phys. Rev. B , 10267 (2000). M. Stone, Phys. Rev. B , 184503 (2012). B. Roy, Phys. Rev. B , 220506 (2020). T. Shang, M. Smidman, S. K. Ghosh, C. Baines, L. J.Chang, D. J. Gawryluk, J. A. T. Barker, R. P. Singh, D.Mck. Paul, G. Balakrishnan, E. Pomjakushina, M. Shi,M. Medarde, A. D. Hillier, H. Q. Yuan, J. Quintanilla,J. Mesot, and T. Shiroka, Phys. Rev. Lett. , 257002(2018). S. Sundar, S Salem-Sugui Jr , M K Chattopadhyay, S BRoy, L S Sharath Chandra, L F Cohen, and L Ghivelder,Supercond. Sci. Technol. R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. , 1440(1977). S. Weinberg, Phys. Rev. Lett. , 223 (1978). F. Wilczek, Phys. Rev. Lett. , 279 (1978). R. Li, J. Wang, X.-L. Qi, and S.-C. Zhang, Nat. Phys. ,284 (2010). P. Goswami, B. Roy, Phys. Rev. B , 041301(R) (2014). M. Stone, P. L. S. Lopes, Phys. Rev. B , 174501 (2016). A. Atland and B. Simons,