Loop currents and anomalous Hall effect from time-reversal symmetry-breaking superconductivity on the honeycomb lattice

We study a tight-binding model on the honeycomb lattice of chiral $d$-wave superconductivity that breaks time-reversal symmetry. Due to its nontrivial sublattice structure, we show that it is possible to construct a gauge-invariant time-reversal-odd bilinear of the pairing potential. The existence of this bilinear reflects the sublattice polarization of the pairing state. We show that it generates persistent loop current correlations around each lattice site and opens a topological mass gap at the Dirac points, resembling Haldane's model of the anomalous quantum Hall effect. In addition to the usual chiral $d$-wave edge states, there also exist electron-like edge resonances due to the topological mass gap. We show that the presence of loop-current correlations directly leads to a nonzero intrinsic ac Hall conductivity, which produces the polar Kerr effect without an external magnetic field. Similar results also hold for the nearest-neighbor chiral $p$-wave pairing. We briefly discuss the relevance of our results to superconductivity in twisted bilayer graphene.


I. INTRODUCTION
Chiral superconductors, which possess order parameters that break time-reversal symmetry, are currently the subject of much attention due to their nontrivial topological properties. 1,2 The best known example of a chiral pairing state is the A phase of superfluid 3 He. 3 Here Cooper pairs have the orbital angular momentum quantum numbers L = 1 and L z = ±1, and the pairing potential has (p x ± ip y )-wave symmetry. A direct solid-state analogue of this phase is believed to be realized in the triplet superconductor Sr 2 RuO 4 . 4 Chiral superconductivity can also be obtained for pairing with higher-order orbital angular momentum. For example, the low-temperature superconducting phase of UPt 3 may realize a chiral f -wave state, 5,6 while chiral d-wave superconducting states have been proposed for URu 2 Si 2 , 7 SrPtAs, 8 and twisted bilayer graphene. 9 Many other materials have been predicted to show chiral superconductivity, such as water-intercalated sodium cobaltate Na x CoO 2 · yH 2 O, 10 the half-Heusler compound YPtBi, 11 and transition metal dichalcogenides. [12][13][14] The breaking of time-reversal symmetry in a chiral superconductor can be revealed by a number of experimental techniques, e.g. muon spin rotation or Josephson interferometry. 2 In the last dozen years, measurements of the polar Kerr effect have emerged as a key experimental probe. 15 It gives evidence for an anomalous ac Hall conductivity at zero external magnetic field, which is a signature of broken time-reversal symmetry. A number of superconductors have been shown to display a nonzero Kerr signal below their critical temperatures, specifically Sr 2 RuO 4 , 16 UPt 3 , 17 URu 2 Si 2 , 18 PrOs 4 Sb 12 , 19 and Bi/Ni bilayers. 20 Although these observations give clear evidence for broken time-reversal symmetry, there is ongo-ing debate over the mechanism underlying the polar Kerr effect in chiral superconductors. An extrinsic Kerr effect may originate from impurity scattering, [21][22][23] whereas an intrinsic Kerr effect is possible for clean multiband superconductors. [24][25][26][27][28][29][30][31] The latter mechanism requires that the pairing potential depends on electronic degrees of freedom beyond the usual spin index, e.g. orbital or sublattice indices. However, it remains unclear what general model-independent conditions these additional electronic degrees of freedom have to satisfy in order to produce a Kerr effect. Here we develop a general condition for this and then apply it to a minimal model of a chiral d-wave superconductor in order to clarify the underlying physics.
Such a minimal theoretical model of a chiral superconductor is provided by the extended Hubbard model on the honeycomb lattice. 32,33 Various theoretical techniques [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] applied to this system generally agree on the existence of a spin-singlet chiral d-wave state at a doping level close to the van Hove singularity. Closer to halffilling, however, different methods have yielded singlet and triplet pairings, 37,44,[47][48][49][50][51] pair-density-wave Kekule order, 52,53 or an unconventional coexistence with antiferromagnetism. [54][55][56] The purpose of our paper is not to further interrogate the phase diagram, but rather to examine the properties of the chiral d-wave state in the case where the nearest-neighbor pairing dominates. Such inter-sublattice pairing would satisfy the multiband requirement 25 for the anomalous Hall conductivity. Thus, chiral d-wave pairing on the honeycomb lattice provides a minimal model of the intrinsic Kerr effect, in contrast to the more complicated multiband models of Sr 2 RuO 4 [24][25][26][27][28] and UPt 3 . [29][30][31] The recent discovery of superconductivity in twisted bilayer graphene, 9 which has been proposed to realize a chiral d-wave state, [57][58][59][60][61][62] makes this study timely. We discuss the relationship between our model and these arXiv:1802.02280v3 [cond-mat.supr-con] 3 Apr 2019 proposals in more detail near the end of the paper. Using this minimal model as an example, we show how to construct a gauge-invariant time-reversal-odd term by taking the product of the pairing potential and its Hermitian conjugate. The existence of such a bilinear is a prerequisite for the experimental detection of time-reversal symmetry-breaking superconductivity in a clean and homogeneous system. In the honeycomb model, the bilinear arises from the varying participation of the two sublattices in the pairing across the Brillouin zone and describes spontaneous breaking of the discrete Z 2 timereversal symmetry. The presence of this term results in the opening of a topological mass gap at the Dirac points and the emergence of persistent loop current correlations, in a striking analogy to Haldane's model of the anomalous Hall insulator. 63 Furthermore, we show that the loop current correlations imply a nonzero anomalous Hall conductivity, hence connecting the polar Kerr effect in superconductors with the time-reversal-odd bilinear product of the pairing potentials.
The paper is organized as follows. We start in Sec. II by introducing the model of spin-singlet chiral d-wave pairing on the honeycomb lattice. In Sec. III we define a gauge-invariant bilinear product of the superconducting pairing potentials that breaks time-reversal symmetry. As a consequence of the existence of this bilinear, we demonstrate the opening of the mass gap at the Dirac point in Sec. IV and the existence of loop currents in Sec. V. The anomalous ac Hall conductivity is calculated in Sec. VI. A phenomenological description of the loop currents is outlined in Sec. VII. The relationship of our work to proposals of chiral d-wave superconductivity in twisted bilayer graphene is discussed in Sec. VIII. We conclude in Sec. IX with a brief discussion of the broader implications of our work. In Appendix A we present similar results for a spin-triplet chiral p-wave state on the honeycomb lattice. In Appendix B we show how the bilinear discussed in Sec. II applies to a broader class of Hamiltonians. More general expressions for the loopcurrent order and the Hall conductivity in the case of inequivalent sublattices are given in Appendix C. The high-frequency small-gap limit of the ac Hall conductivity is derived in Appendix D.

II. MICROSCOPIC MODEL
The Bogoliubov-de Gennes (BdG) Hamiltonian of superconducting pairing on the honeycomb lattice is where Ψ k = (a k,↑ , b k,↑ , a † −k,↓ , b † −k,↓ ) T , and the operator a k,σ (b k,σ ) annihilates an electron with momentum k = (k x , k y ) and spin σ on the A (B) sublattices. In Eq. (1), H 0 (k) and ∆(k) are 2×2 matrices in the sublattice space, and the absence of spin-orbit coupling allows the spin variables to be factored out.
Using the Pauli matrices s λ to encode the sublattice degree of freedom, we write the normal-state Hamiltonian as Here µ is the chemical potential, t is the nearest-neighbor hopping amplitude, and the R j are the vectors of length a connecting an A site to its neighboring B sites, see Fig. 1.
For generality, we also include the Semenoff term 64 as a staggered potential δ s . This makes the A and B sites inequivalent, hence breaking inversion symmetry and lowering the point group from D 6h to D 3h . We consider chiral spin-singlet superconducting pairing on the nearest-neighbor bonds along the directions R j shown in Fig. 1. This gives the pairing term The magnitude ∆ of the pairing potential is the same for each bond j, but the phase is φ j = (j − 1)2π/3, and the two choices of sign in Eq. (4) define degenerate pairing potentials with opposite chiralities. A similar chiral spintriplet pairing is discussed in Appendix A. The two pairing potentials in Eq. (4) can be written in terms of the Pauli matrices s λ where and a is the distance between neighboring sites. When projected onto the states near the Fermi surface, the basis states ∆ x 2 −y 2 (k) and ∆ xy (k) have the forms of d x 2 −y 2 and d xy waves, so ∆ ± (k) can be regarded as a chiral (d x 2 −y 2 ± id xy )-wave pairing state. The matrices s x and s y are multiplied by the functions that are even and odd with respect to k → −k. This ensures that the pairing potentials are even under inversion, e.g.

III. TIME-REVERSAL-ODD BILINEAR
A central goal of our work is to understand how broken time-reversal symmetry in the particle-particle superconducting channel can lead to observable effects in the particle-hole channel, e.g., the anomalous Hall conductivity and the polar Kerr effect. For such effects, it is not sufficient to consider the pairing potential ∆(k) alone, since it is not gauge-invariant. Rather, these observable must depend on a time-reversal symmetry-breaking bilinear combination of ∆(k) and ∆ † (k).
In order to define the time-reversal operation, let us label the second-quantized electron operators ψ c,σ (k) by the sublattice index c and the spin index σ = ±. The time-reversal operation involves the substitution ψ c,σ (k) → σψ c,−σ (−k) and complex conjugation of the matrix elements in the BdG Hamiltonian. 66 Then the off-diagonal term in Eq. (1) transforms as follows where summation over repeated indices is implied. Note that, to obtain the second line, we anticommuted the fermion operators and then swapped the sublattice indices. Thus we obtain a BdG Hamiltonian of the same form with upon time reversal. The simplest bilinear product of the pairing potential with its Hermitian conjugate is ∆(k)∆ † (k). 67 The timereversal-odd part of this bilinear product, which we abbreviate as TROB, obtained as the difference between ∆(k)∆ † (k) and its time-reversed counterpart, is a commutator (12) Due to its gauge invariance and odd time-reversal behavior, a non-zero TROB permits broken time-reversal symmetry in the particle-particle channel to manifest in the particle-hole channel. In Appendix B, we show that the expression for the TROB in Eq. (12) applies to more general Hamiltonians, which may include spin-orbit coupling and more electronic degrees of freedom, or break inversion symmetry. In the second-quantized formalism, the TROB matrix from Eq. (12) appears in the timereversal-odd part of the commutator of the pairing terms where Θ is the time-reversal operation, and We immediately see that the TROB in Eq. (12) always vanishes for a single-band spin-singlet superconductor where ∆(k) is just a complex function. This implies that any probe of time-reversal symmetry breaking, e.g., the Hall conductivity or polar Kerr effect, must vanish if such a system is clean, so that the momentum k is a good quantum number. Hence, the experimental detection of time-reversal symmetry breaking in single-band superconductors must rely upon inhomogeneities not conserving k, e.g., scattering off impurities. 22 However, for a clean multiband system, where the pairing potential can be expressed as a matrix in the band indices, it is possible for the commutator in Eq. (12) to take on nonzero values. This is the case for the honeycomb lattice model of Sec. II, for which we obtain where the wedge product [a ∧ b] = a x b y − a y b x is used for the two-component vector ∆ k = (∆ x k , ∆ y k ) from Eqs. (7) and (8). In the second line, the sum is taken over the pairs of nearest-neighbor bonds in Fig. 1. The nonzero TROB in Eq. (15) implies the existence of a time-reversal symmetry-breaking sublattice polarization of the pairing state, which we define as The sublattice polarization Ξ ± (k) is crucially important for the physical effects discussed in the rest of the paper. It quantifies the relative participation in the pairing of electrons on the A and B sites. Pairing at the M point of the Brillouin zone involves both sublattices equally, and so Ξ ± (k M ) = 0. In contrast, pairing at the K (K ) point involves exclusively the B (A) sublattice for the ∆ + potential, and so Ξ + (k K ) = −Ξ + (k K ) = 9|∆| 2 ; the sublattice polarization is reversed for ∆ − . 56 This can be considered as a generalization to non-spin internal degrees of freedom of the spin polarization of a single-band nonunitary triplet state. 68 It has recently been pointed out that such a polarization generically arises in multiband time-reversal symmetry-breaking superconductors, where it can have dramatic effects on the low-energy nodal structure. 69 Although the effect of the polarization on the electronic structure is confined to high energies in our fully-gapped pairing state, we shall see below that it plays a key role in generating the Hall conductivity. Further insight into the implications of a nonzero TROB is provided by the concept of the superconducting fitness, which has recently emerged as a way to characterize the pairing state in multiband materials. 70,71 For our system, where the normal-state Hamiltonian H 0 (k) is time-reversal invariant, a superconducting state is said to have perfect fitness when i.e., the normal-state Hamiltonian H 0 (k) commutes with the pairing potential ∆(k). Then these two matrices can be simultaneously diagonalized in the normal-state band basis, and so there is no interband pairing in the case of perfect fitness. In this basis, a multiband BdG Hamiltonian with even-parity spin-singlet pairing reduces to a collection of decoupled single-band terms, so the TROB must therefore vanish. We conclude that the lack of perfect fitness, i.e., a violation of Eq. (17) and the presence of interband pairing, is a necessary (but not sufficient) condition for a nonvanishing TROB. The presence of interband pairing has been previously noted as crucial for the existence of the polar Kerr effect in clean chiral superconductors. 24,25 The chiral d-wave pairing potential in our model does violate the superconducting fitness condition: where k = ( x k , y k ), complex conjugation in the righthand side applies only for the negative chirality, and we set the Semenoff term to zero for simplicity (i.e., δ s = 0). The violation of the fitness condition is due to the nontrivial phases φ j along the nearest-neighbor bonds, making the complex vector ∆ k in Eqs. (7) and (8) not parallel to the real vector k in Eq. (3). Although the presence of both intraband and interband pairing in the chiral dwave state is energetically disadvantageous due to mismatch of the energies of different bands, 72 it can emerge in a mean-field BCS theory due to the short range of realspace interaction between electrons. Indeed, the pairing potential Eq. (4) naturally arises from a mean-field decoupling of the nearest-neighbor exchange interactions in a t-J model. 46 It is instructive to compare our results to a chiral dwave state with purely intra-sublattice (i.e., next-nearestneighbor) pairing, as proposed in Ref. 39. For this state, the pairing potential is proportional to the unit matrix in sublattice space. As such, this potential commutes with the normal-state Hamiltonian and so has perfect fitness. Thus, despite the fact that∆ ± (k) breaks time-reversal symmetry and has a nonzero phase winding around the Fermi surface, this state does not display an intrinsic polar Kerr effect because TROB = 0. The pairing potential ∆ ± (k) and the TROB describe spontaneous breaking of the continuous U (1) gauge symmetry and the discrete Z 2 time-reversal symmetry, respectively. In the mean-field BCS theory, both symmetries are broken simultaneously. In a more general framework, however, these two symmetries may be broken at separate phase transitions taking place at different temperatures. For example, the TROB may acquire a nonzero expectation value at a higher temperature by selecting positive or negative chirality (which can be detected experimentally by observing the polar Kerr effect), while the expectation value of the pairing potential ∆ ± (k) is still zero due to phase fluctuations. The superconducting properties, such as supercurrent and Meissner effect, would emerge at a lower temperature, where ∆ ± (k) acquires a non-zero expectation value. This scenario is discussed in more detail in Sec. VII.
The above considerations are not limited to spinsinglet even-parity superconductivity. In Appendix A, we show that an odd-parity spin-triplet chiral p-wave pairing has a similar TROB and sublattice polarization.

IV. TOPOLOGICAL MASS GAP
Let us set the Semenoff term in Eq. (2) to zero first: δ s = 0. In this case, there is no gap at the Dirac points K and K in the normal state. However, the energy spectrum of the BdG Hamiltonian in Eq. (1) shows an unexpected gap opening at the Dirac points near E = ±µ = ±0.5t in Fig. 2(a), far away from the usual superconducting gap at the Fermi level E = 0. Note that the momentum q = k − k K is measured relative to the K point in Fig. 2.
To gain insight into the nature of this unexpected gap, we derive an effective Hamiltonian for the states near the Dirac points, perturbatively including the superconducting pairing in the limit ∆ |µ|. Our starting point is the formal expression for the electron-like component of the Green's function, To find the energy spectrum in the vicinity of the Dirac points, we replace ω → −µ in the last term of Eq. (20) and obtain an effective Hamiltonian Near the K point, we can expand the first term to linear order in the relative momentum with the correction due to superconductivity Near the K point, the expansion of the unperturbed Hamiltonian is identical except for the reversed sign in front of q y , and the correction is (24) Note that Eqs. (23) and (24) can be obtained from the last term in Eq. (20) only in the vicinity of the Dirac points, where H 0 in Eq. (22) is proportional to the unit matrix in the limit of vanishing q. Equation (21) can be interpreted as an effective normal-state Hamiltonian with the second-order perturbative correction due to superconducting pairing.
The perturbative correction given by Eqs. (23) and (24) is proportional to the matrix product ∆(k)∆ † (k). Its time-reversal-even part, proportional to the unit matrix s 0 , shifts the energy of the Dirac point. In contrast, the time-reversal-odd part (i.e., the TROB), proportional to s z , opens a mass gap. This demonstrates the appearance of time-reversal symmetry breaking in the particlehole channel due to the nonzero TROB. The gapped energy dispersion derived via this perturbative argument, shown by the dashed red curve in Fig. 2(a), is in excellent agreement with the dispersion of the full model near to the Dirac point.
The mass gaps at the K and K points introduced by the superconductivity [Eqs. (23) and (24)] have opposite signs. This suggests a topologically nontrivial state, as in Haldane's model of the quantum anomalous Hall state on the honeycomb lattice. 63 The topological nature of the mass gap is confirmed by calculation of the Chern numbers for the different bands and observation of chiral edge states within the energy gaps via the bulk-boundary correspondence. With the opening of the mass gap, the four eigenstates of the BdG Hamiltonian in Eq. (1) are everywhere nondegenerate, so a Chern number ν α can be defined for each band α = 1, 2, 3, 4, as labeled in Fig. 2(a). As shown in Fig. 3, each band has a nonzero Chern number for δ s = 0, i.e. in the absence of the Semenoff term. The sum of the Chern numbers of the occupied bands 3 and 4 below the chemical potential is −2, consistent with the chiral d-wave superconductivity. Correspondingly, the two topologically-protected chiral edge states within the superconducting gap are clearly visible near E = 0 in the energy spectrum weighted by the integrated probability density of the electron-like wave function components near the surface, as shown in Fig. 2(b) for the armchair edge. The nonzero Chern numbers of the outer bands 1 and 4, which are separated by the mass gap from the inner bands 2 and 3, imply that the mass gap is topological. Thus, we would expect to find a single chiral edge state within each mass gap. However, due to the spectrum doubling in the superconducting state, the hole-like states generally overlap with the energy range of the mass gap and can hybridize with the edge state. Nevertheless, the edge state persists as a predominantly electron-like edge resonance inside the mass gap between bands 1 and 2 in Fig. 2(b).
A combination of a nonzero Semenoff term δ s = 0 in Eq. (2) and the superconducting corrections in Eqs. (23) and (24) produces different magnitudes of the mass gaps at the two Dirac points K and K . At a critical value δ c = 1 2 ( 9|∆| 2 + µ 2 − |µ|), the gap at one of the Dirac points passes through zero and changes sign. Correspondingly, as shown in  of the Chern numbers of the occupied bands 3 and 4 remains −2. This is consistent with the mass gaps at K and K having the same sign, which is topologically trivial. Accordingly, we do not observe any edge resonance within the gap, as shown in Fig. 2(c).
Repeating the calculations for a zigzag edge, we also find evidence for Haldane states. However, they are mixed with the standard flat-band edge states that exist at the zigzag edges of a hexagonal lattice, making their interpretation more complicated.

V. LOOP CURRENTS
It was argued in the previous Section that the energy gaps observed at the Dirac points are similar to the energy gaps in Haldane's model of the quantum anomalous Hall insulator. 63 They arise in Haldane's model due to the presence of a time-reversal symmetry-breaking nextnearest-neighbor hopping term, resulting in loop currents around each lattice site shown by the arrows in Fig. 4. In second quantization, the time-reversal-odd part of this hopping term is proportional to the dimensionless operator χ lc , which is defined as Here c 1 and c 2 are the primitive lattice vectors (see Fig. 1), the operator a r,σ (b r,σ ) destroys a spin-σ electron on the A (B) site of the unit cell corresponding to the lattice vector r, and τ 0 is the unit matrix in Nambu space. The sign convention in Eq. (25) matches the convention for the link directions in Fig. 4. In Haldane's model, the operator in Eq. (25) has a nonzero expectation value χ lc = 0, resulting in the loop currents shown in Fig. 4. In our model, the operator χ lc appears in the commutator (13) with the TROB given by Eq. (15). The commutator of the pairing terms on the adjacent nearest-neighboring links generates electron transfer between the next-nearest neighboring sites with a complex amplitude carrying the phase difference of pairing potentials shown in Fig. 1. The analogy between our system and Haldane's model suggests that χ lc also has a non-zero expectation value in the chiral d-wave state, which is readily verified to be Here E k,α=1,2 > 0 are the quasiparticle dispersions corresponding to the upper two bands shown in Fig. 2 Since Ξ ± (k) has the same momentum dependence as the term in front of the fraction in Eq. (26), the summand has the same sign everywhere in the Brillouin zone, and thus the expectation value χ lc is nonzero. The essential importance of the TROB in ensuring χ lc = 0 is consistent with the role of the TROB in generating the energy gaps at the Dirac points. As such, the inclusion of the Semenoff term does not alter the conditions for a nonzero expectation value of χ lc , but Eq. (26) is replaced by a more complicated expression given in Appendix C. We note that χ lc = 0 was calculated in Ref. 51 for the closely-related nearest-neighbor chiral p-wave state introduced in Appendix A. In contrast, χ lc is zero for the intra-sublattice chiral d-wave state described by Eq. (19), where the TROB vanishes. However, unlike in Haldane's model, a nonzero expectation value of χ lc in our system only implies the presence of loop current correlations. Since the normal-state Hamiltonian (2) does not contain next-nearest-neighbor hopping, there are no current operators between nextnearest sites in our model. We can remedy this by introducing a next-nearest-neighbor hopping term with a small real amplitude t : The corresponding current operators between the nextnearest-neighbor sites m and n belonging to the sublattice c = A, B are hence obtained from Eq. (27) as where ψ c,σ (r) is the annihilation operator for spin-σ electrons on sublattice c of unit cell r. Adding the current operators in Eq. (28) with the signs corresponding to Fig. 4, we introduce the total current operator I tot as Then the expectation value I of the microscopic current on one link is obtained as where we divide by 6 because there are six currents of equal magnitude in a unit cell. The current I is very small, because it is proportional to the small hopping amplitude t and χ lc /N ∼ µ∆ 2 /t 3 1 from Eq. (26). Another physical consequence of χ lc = 0 is the existence of a nonzero anomalous Hall conductivity in the absence of an external magnetic field, which is calculated in the next Section. Unlike the current I in Eq. (30), the Hall conductivity does not require t = 0, so we set t = 0 in the rest of the paper to simplify calculations.

VI. HALL CONDUCTIVITY
The existence of loop-current correlations in Eq. (26) for the chiral d-wave state naturally suggests the presence of an intrinsic Hall conductivity. Indeed, the nontrivial sublattice structure of the BdG Hamiltonian (1) is consistent with the conditions outlined in Ref. 25 for the existence of an intrinsic Hall effect.
As shown by the Feynman diagrams in Fig. 5, the Hall conductivity can be obtained as the difference of the current-current correlation functions where ω n is a bosonic Matsubara frequency and S is the total area of the crystal. Here J a is the a-component of the current operator where V k is the velocity vertex in Nambu notation and the velocity components are obtained from Eq. (3).
A straightforward evaluation of the Feynman diagrams in Fig. 5 (for the vanishing Semenoff term δ s = 0) yields the Hall conductivity The sign of σ H (ω) correlates with the sign of the chemical potential µ, and the Hall conductivity vanishes at µ = 0 (at the Dirac point) due to particle-hole symmetry. The real and imaginary parts of the Hall conductivity calculated from Eq. (36) are shown in Fig. 6. This expression is consistent with Eq. (24) of Ref. 29 for the Hall conductivity in UPt 3 in the limiting case where spin-orbit coupling and intrasublattice hopping terms are neglected. As the point groups of UPt 3 and the honeycomb lattice are both D 6h , such terms are also allowed in our model, but we neglect them for simplicity. Equation (36) shows some similarity to a general formula 73 for the intrinsic ac Hall conductivity in terms of the Berry curvature for a nonsuperconducting two-band system, which includes Haldane's model. However, that formula is not directly applicable to our superconducting case, because the effective two-band model derived in Sec. IV is only suitable near to the Dirac points.
From the numerator of Eq. (36), it is clear that the anomalous Hall conductivity σ H (ω) is nonzero only when the sublattice polarization Ξ ± (k) in Eq. (16) has common irreducible representations with the product of velocities where ϕ j are the geometric angles between the vectors R j and the x axis, which, in our model, are the same as the phases φ j in Eq. (4). Thus, Eq. (37) has the same momentum dependence as the TROB in Eq. (15) and Ξ ± (k) in Eq. (16).
The full result Eq. (36) is rather complicated, but it simplifies in the high-frequency limit ω t. This regime is experimentally relevant, as the polar Kerr effect measurements detailed in Ref. 15 are performed at infrared frequency ω = 0.8eV which is typically large compared to the hopping integrals in a strongly-correlated material. As shown in Ref. 74, the Hall conductivity in this limit is given by Taking into account Eq. (37), we find that the commutator of the x-and y-components of the current operator appearing in this expression is directly proportional to the loop-current operator in Eq. (25) [ We hence find that the high-frequency Hall conductivity is proportional to the expectation value of the loopcurrent operator: . Equations (38)-(40) establish a direct connection between the Hall conductivity and the loop currents discussed in Sec. IV. As shown in Fig. 6, the agreement between Eqs. (36) and (40)  The inset compares the approximate and exact results in more detail. We use the same parameters as in Fig. 2(a) and set the temperature kBT = 0.05t.
limit of small ∆, where the Green's functions appearing in Fig. 5 can be expanded to the second order in the pairing potential. This approach yields Eq. (D3) similar to Eq. (40), but with the BdG energies E k replaced by the normal-state energies k . The high-frequency Hall conductivity in Eq. (40) is real, but the polar Kerr effect is primarily sensitive to the imaginary part of the Hall conductivity when the refraction index is predominantly real. 15 Although it is not possible to directly associate the imaginary part of the Hall conductivity at a given frequency to the loop currents in the superconductor, an indirect connection is provided by the sum rule 75 Again, the right-hand side of this equation is proportional to the expectation value of the loop-current operator, and we hence conclude that the existence of the loop-current correlations results in a nonzero imaginary Hall conductivity.
It should be noted that, in contrast to nonsuperconductors, the dc Hall conductivity in superconductors is not directly related to the Chern number, as discussed in Appendices A and B of Ref. 22. Thus, the topological phase diagram shown in Fig. 3 in terms of the Semenoff term δ s is not particularly relevant for the calculation of the Hall conductivity. A generalization of Eq. (36) to a nonzero Semenoff term in Appendix C shows that the ac Hall conductivity σ H (ω) is nonzero for any value of δ s . Moreover, in the high-frequency limit, Eqs. (38) and (39) are still valid for δ s = 0, so the Hall conductivity remains proportional to the expectation value of the loop-current operator, which is mainly sensitive to the pairing potential and only weakly dependent upon the Semenoff term.

VII. PHENOMENOLOGICAL TREATMENT
In Eq. (26) we obtained a nonzero expectation value of the loop-current operator from a microscopic theory of the chiral d-wave state at the level of the singleparticle Green's functions. The appropriate interactions would, however, lead to true long-range loop-current order. The interplay of this order with the superconductivity could then be understood within the framework of a phenomenological Landau expansion of the free energy density where f 0 is the normal-state free energy density. The first two lines describe the superconductivity, where η 1 and η 2 are the order parameters corresponding to the two states in the E 2g irreducible representation. The term with β 2 > 0 stabilizes the time-reversal symmetrybreaking configuration (η 1 , η 2 ) ∝ (1, ±i) studied in this paper. The coupling to the loop-current order parameter δ lc ∝ χ lc is given by γ, and κ > 0 implies that this order is subdominant. Minimization of f with respect to δ lc shows that the loop-current order becomes induced in the time-reversal-breaking superconducting state. As already mentioned in Sec. III, a more intriguing possibility could be that the loop-current order preempts the superconductivity. We speculate that fluctuating superconducting order may cause the discrete Z 2 time-reversal symmetry to be broken with δ lc = 0 at a higher temperature than the continuous U (1) gauge symmetry, which is rigorously permitted only at zero temperature in two dimensions. Similar scenarios were discussed for multiband superconductors in Ref. 76 and for pair-density wave order in the underdoped cuprates in Ref. 77.

VIII. RELEVANCE TO SUPERCONDUCTIVITY IN TWISTED BILAYER GRAPHENE
It has been proposed theoretically 57-62 that the superconducting state observed in twisted bilayer graphene 9 (TBLG) realizes chiral d-wave pairing. Given that our analysis concerns hypothetical superconductivity in monolayer graphene, it is worthwhile to survey theories of TBLG briefly and explore possible links to our work. Although most proposals include more electronic degrees of freedom than our model, in many cases they show a qualitative resemblance, implying that the physics discussed in our paper may be applicable.
Some of the earliest proposals, such as Refs. 57 and 58, assume SU (4) symmetry of a single-particle Hamiltonian, for which the physics we discuss does not apply. However, SU (4) symmetry-breaking terms may change this conclusion and are currently under consideration.
Phenomenological models with orbital or sublattice degrees of freedom have been considered in Refs. 78 and 59, respectively. Due to the presence of these additional electronic degrees of freedom, the pairing potential may have a nonzero TROB, thus resulting in similar physics to that discussed here. The model of Ref. 59, based upon a three-site sublattice to simulate the AA and AB regions of the Moiré pattern of TBLG, resembles our model most closely and, indeed, reduces to it in the limit t = ∆ j = 0, where the triangular lattice of the AA regions is neglected.
Several papers 79-82 proposed a low-energy description of the normal-state electronic structure based on an emergent honeycomb lattice with two additional electronic orbital degrees of freedom at each lattice site. Reference 60 considers such a model with p x and p y orbitals on a honeycomb lattice with local electronic interactions and finds that a d-wave chiral pairing state emerges. The paper argues that the mechanism is similar to that for chiral d-wave pairing in single-layer graphene at quarter doping, for which our model applies.
Adopting an alternative approach, Refs. 61 and 62 numerically analyze how a nearest-neighbor chiral d-wave state in each of the two graphene layers is modified by the Moiré structure in TBLG. Interestingly, these papers find intra-unit-cell supercurrent loops. Since the intralayer pairing state is identical to ours, we suggest that these currents may be related to the loop current correlations found in our work.
Although a theoretical description of the superconducting state in TBLG remains unsettled, the presence of multiple electronic degrees of freedom indicates that a chiral d-wave state is likely to have a nonvanishing TROB. Consequently, much of the physics discussed in our paper may be applicable. An experimental measurement of the polar Kerr effect in TBLG would be particularly useful to verify whether the superconducting pairing breaks time-reversal symmetry.

IX. CONCLUSIONS
In this paper we have examined the appearance of the polar Kerr effect in a minimal model of time-reversal symmetry-breaking chiral d-wave superconductivity on the honeycomb lattice. We have demonstrated that the existence of a gauge-invariant time-reversal-odd bilinear (TROB) constructed from the pairing potential is an essential requirement for the polar Kerr effect. In the context of the honeycomb lattice, the TROB reflects the sublattice polarization of the pairing. The key physical manifestation of the TROB is the appearance of an emergent nonsuperconducting order in conjunction with the superconductivity, which we identify as loop currents similar to those in Haldane's model of a quantum anomalous Hall insulator. 63 This is directly evidenced in the energy spectrum, where we observe the opening of a topological gap with opposite signs at the Dirac points K and K . The Kubo formula calculation of the intrinsic ac Hall conductivity in the absence of an external magnetic field shows that it is directly proportional to the expectation value of the loop current operator. Thus we establish an explicit relation connecting the emergent loop-current correlations and both the real and imaginary parts of the Hall conductivity.
The model considered here is another example of a time-reversal symmetry-breaking superconducting state with an intrinsic Hall conductivity, and generalizes these analyses to an even-parity pairing state. The first example is Sr 2 RuO 4 , where different pairing in the Ru d xz and d yz orbitals implies a polarization in the d xz -d yz orbital space. [24][25][26][27][28] More recently, a theoretical treatment of the Kerr effect in UPt 3 has identified the time-reversal-odd sublattice dependence of the pairing potential permitted by the nonsymmorphic symmetry as an essential ingredient. 29,65 The similarity of these models to the simpler case considered in our paper suggests the intriguing possibility that the Hall conductivities in these systems can be also understood in terms of loop-current correlations induced by a TROB. Although we have considered only two internal degrees of freedom here, loop currents can also arise in materials with more complicated unit cells, 83 including twisted bilayer graphene. The observation of the polar Kerr effect in many unconventional superconductors therefore suggests that pairing states supporting nonzero TROBs may be realized in a broad range of materials.

ACKNOWLEDGMENTS
The authors acknowledge useful discussions with Fengcheng Wu, Carsten Timm and Henri Menke. We thank the hospitality of the summer 2016 program "Multi-Component and Strongly-Correlated Superconductors" at Nordita, Stockholm, where this work was initiated. DSLA acknowledges support of ERC project DM-321031. DFA acknowledges support from the NSF via DMREF-1335215. Other chiral pairing states on the honeycomb lattice also display a polarization in their sublattice degrees of freedom due to a non-zero TROB. For example, chiral p-wave triplet superconductivity with nearest-neighbor pairing has been considered by several authors. 47,49,51,56 Here we examine the case where the vector d of the triplet pairing is oriented along the z axis, which is also the spin quantization axis. In this representation, the triplet pairing takes place between the opposite spins and is described by the BdG Hamiltonian in Eq. (1) with the pairing potential Note that one of the off-diagonal terms has opposite sign compared with Eq. (4) for the pairing potential of the chiral d-wave state, but the phases along each bond φ j = (j − 1)2π/3 are the same. As shown in Fig. 7, the phase of the pairing on each bond winds by 2π as one moves around a hexagonal plaquette, in contrast to the chiral d-wave state where the phase winds by 4π. The pairing potential in Eq. (A1) can be decomposed into basis states of the irreducible representation E 1u Projected onto the states at the Fermi surface, the basis functions ∆ p x (k) and ∆ p y (k) appear as p x -wave and p y -wave triplet states, respectively. Like the basis functions for E 2g discussed in Sec, II, these states contain the matrices s x and s y , but here with odd-and even-parity coefficients, respectively. This ensures that the pairing potentials are odd under inversion, i.e. I † ∆ p x (k)I = −∆ p x (−k). The E 1u basis functions can be obtained from the basis functions for E 2g by multiplying them with s z . This follows from the direct product rules for the point group D 6h , since s z belongs to the irreducible representation B 1u and E 1u = B 1u × E 2g . Thus, the TROB of the chiral p-wave state is the same as for the chiral d-wave case in Eq. (15): The physics arising from the existence of the TROB in the chiral p-wave state is thus essentially the same as for the chiral d-wave state discussed in the main part of the paper.
We finally note that, in the presence of a Semenoff term, the reduced symmetry of the lattice due to lack of inversion implies that both the p-wave and d-wave pairing potentials are basis states of the same irreducible representation E of the point group D 3h . As such, a chiral state can generally involve a mixture of the two. 54,56 The time-reversal rule in Eq. (B4) is similar to that for a non-superconducting Hamiltonian, 66 because U is real.
According to Eq. (B4), the pairing potential transforms aŝ where we used the fermion exchange relation. It is convenient to define ∆(k) without the "hat" as which corresponds more closely to the pairing potential introduced in Sec. II. Combining Eq. (B5) with Eq. (B6) and using U 2 = −1, we reproduce Eq. (11) Finally, we obtain the TROB as the difference between ∆(k)∆ † (k) and its time-reversed counterpart: where we used Eqs. (B5) and (B6). Equation (B8) reproduces Eq. (12) for the TROB, which is the result we wanted to show in general form.
As discussed in Sec. III, the superconducting fitness restricts opportunities for a nonzero TROB. The general condition for superconducting fitness 70,71 is Here ∆ Clearly, TROB = 0 only when the vector d * k is not parallel to d k . The pairing with d k × d * k = 0 is known in the literature 68 as nonunitary pairing, because the product ∆(k)∆ † (k) is not proportional to the unit matrix. Obviously, TROB in Eq. (B8) vanishes for a unitary pairing, so nonunitarity is a necessary, but generally not sufficient condition for TROB = 0.
A similar construction can be obtained for spin-singlet pairing in a two-band model, where the pairing potential is expanded in the Pauli matrices s λ for sublattice space: The honeycomb lattice model described in Sec. II is a special case where ∆ (0) = 0 and the vector ∆ k has only two components. An evaluation of Eq. (B8) for Eq. (B13) gives a formula similar to Eq. (B12) For TROB = 0, the pairing vector ∆ k must be not parallel to its complex conjugate.
To evaluate the fitness condition for this model, we take the normal-state Hamiltonian to be spinindependent and write it as where h (0) k and h k are necessarily real, because H 0 is Hermitian. Then the fitness condition (B10) is Perfect fitness is achieved when the pairing vector ∆ k is parallel to the real vector h k , which makes TROB vanish in Eq. (B14), so perfect fitness is incompatible with TROB = 0.

Appendix C: Effect of Semenoff term
For simplicity, the main text gives Eq. (26) for the loop-current operator expectation value and Eq. (36) for the Hall conductivity only in the absence of the Semenoff term, i.e. for δ s = 0. Here we present the general expressions for δ s = 0, which may be useful for applications to transition metal dichalcogenides, [12][13][14] where the A and B sites are strongly inequivalent.
In the presence of the Semenoff term, the expectation value of the loop-current operator is given by where the coefficients of the quartic polynomial in the fermionic frequency ν m in the denominator are The numerator in Eq. (C1) is no longer directly proportional to the sublattice polarization Tr{∆ † (k)s z ∆(k)}, but now also contains a term proportional to the Semenoff term δ s . Nevertheless, the contribution from this additional term, which is also proportional to fermionic frequency ν m , is only nonzero if the coefficient c 1 of the linear term in the denominator is also nonzero. As this is only the case if the pairing potential has a time-reversal symmetry-breaking sublattice polarization, the key role of the non-zero TROB in producing the loop current correlations is robust to the presence of the Semenoff term.
The Hall conductivity in the presence of the Semenoff term is given by .
Similarly to Eq. (C1), a nonzero Semenoff term again results in a new term proportional to δ s in the numerator. The coefficient of δ s in the numerator of Eq. (C2) has the full symmetry of the lattice, whereas the prefactor [v k ∧ v * k ] belongs to the irreducible representation A 2 of the point group D 3h . The contribution from this new term will thus be vanishing, unless the denominator also contains a term in the irreducible representation A 2 . Such a term is only present if the linear coefficient c 1 of the polynomials in the denominator is nonzero, which requires a sublattice polarization of the pairing. Thus, the nonzero Hall conductivity remains a signature of a finite TROB in the presence of the Semenoff term.
Appendix D: High-frequency small-∆ limit of the Hall conductivity The high-frequency limit of the Hall conductivity was derived in Ref. 74 from the general form of the currentcurrent correlation function. Here we present an alternative derivation based upon approximation of the Green's functions in the Feynman diagrams shown in Fig. 5. Specifically, in the high-frequency limit |ω| |∆|, the Hall conductivity Eq. (36) should only weakly depend upon the modification of energy spectrum in the superconducting state. We thus expect that a perturbative expansion in the pairing Hamiltonian will quickly converge. To achieve this, we first note that the full Green's function G is related to the Green's function of the normal system G 0 by the Dyson's equation where H ∆ is the pairing part of the BdG Hamiltonian Eq. (1). Expanding this to second order in H ∆ , we approximate G by Note that the normal part of the Green's function in Eq. (D2) reproduces Eq. (20). Using the approximate Eq. (D2) to replace the full Green's function in the current-current correlator π xy (iω n ), we obtain the expansion shown in Fig. 8. Performing the analytic continuation iω n → ω + i0 + , the first diagram on the right hand side is ∼ 1/ω in the high-frequency limit, as the external frequency passes through a single normal-state Green's function. The next diagram is also ∼ 1/ω, since a redefinition of the internal frequency (see second line) also allows the external frequency iω n to pass through a single normal-state Green's function. In contrast, the external frequency in the third diagram must necessarily pass through two Green's functions, and this diagram can be shown to be at least ∼ 1/ω 2 in the high-frequency limit.
Keeping only the first two diagrams, therefore, we approximate π xy (iω n ) as shown in the second line of Fig. 8. We observe that, in performing the Matsubara summation over the internal frequency, the residue of the poles of the Green's function containing the external frequency will be at least ∼ 1/ω 3 , whereas the residue of the poles of the other Green's functions will be ∼ 1/ω. Since the ∼ 1/ω contribution only arises from the unit matrix (i.e. τ 0 ⊗ s 0 ) component of the Green's function containing the external frequency, we make the approximation G 0 (k, iν m ± iω n ) ≈ (±iω n ) −1 τ 0 ⊗ s 0 and hence factor the external frequency out of the Matsubara sum. This yields the diagram involving the commutator of the velocity vertices and the second-order Green's function correction G 0 H ∆ G 0 H ∆ G 0 . This product is proportional to the expectation value of the dimensionless loop-current operator Eq. (26) expanded to lowest order in the pairing potential. Evaluating this diagram, we obtain the Hall conductivity π xy (iω n ) ≈ e 2 iω n 1 Sβ k,νm , where k,1(2) = +(−) ( x k ) 2 + ( y k ) 2 − µ are the dispersions in the normal state. A similar analysis yields 8. Diagrammatic derivation of the high-frequency small-∆ limit. The current-current correlator πxy(iωn) is expanded in powers of the pairing Hamiltonian H∆, which is treated as a perturbation. The leading-order contribution in the high-frequency limit comes from the first two terms on the right hand side, as shown in the second line. Note the redefinition of the internal frequency iνm = iνm + iωn in the second term. After the Green's functions containing the external frequency are factored out, the result is expressed in terms of the expectation value of the loop-current operator χ lc . The double line represents the full Green's function G, and the single line denotes the Green's function G0 of the normal system. π yx (iω n ) = −π xy (iω n ). We hence obtain the Hall con-ductivity .
We recognize this as the lowest-order term in the expansion of Eq. (40) in powers of the gap magnitude.