d-wave superconductivity in the frustrated two-dimensional periodic Anderson model

Superconductivity in heavy-fermion materials can sometimes appear in the incoherent regime and in proximity to an antiferromagnetic quantum critical point. Here we study these phenomena using large scale determinant quantum Monte Carlo simulations and the dynamical cluster approximation with various impurity solvers for the periodic Anderson model with frustrated hybridization. We obtain solid evidence for a $d_{x^2-y^2}$ superconducting phase arising from an incoherent normal state in the vicinity of an antiferromagnetic quantum critical point. There is a coexistence region and the width of the superconducting dome increases with frustration. Through a study of the pairing dynamics we find that the retarded spin-fluctuations give the main contribution to the pairing glue. These results are relevant for unconventional superconductivity in the Ce-$115$ family of heavy-fermions.


I. INTRODUCTION
d-wave superconductivity in proximity to a quantum critical point has been found in many compounds, such as layered and quasi one-dimensional organic superconductors, the iron pnictides, high-temperature superconductors and in the heavy fermion compounds. There are many examples of quantum critical points in heavy fermion materials so they are an important testing ground for theories of quantum criticality and its relation to superconductivity. [1] This is the main problem that we consider.
Heavy fermion behavior arises when highly localized felectron bands are hybridized with conduction-electron bands. This hybridization leads to the RKKY interaction between f electron local moments but also to the eventual screening of the local moments by the Kondo effect. These competing tendencies are summarized by Doniach's phase diagram [1] where an antiferromagnetic and a Fermi-liquid phase both appear. The zerotemperature phase transition between these phases is generally believed to be a quantum critical point (QCP).
Remarkably, superconductivity appears in the vicinity of this QCP, [2] notwithstanding the exact nature of this QCP as determined by the Kondo deconstruction energy scale E * loc . [3,4] [5-7] Amongst the many heavy fermion superconductors, [8] there are several compounds, for example UBe 13 , PuCoGa 5 and CeCoIn 5 , that show the peculiarity [1] that d x 2 −y 2 pairing symmetry develops out of an incoherent metallic state. [9][10][11][12] The latter quasi-2D Ce-115 materials have an easily accessible superconducting transition temperature (∼ 2K) where superconducting and magnetic properties can be precisely measured [13] and they are especially interesting, because the itinerant-tolocalized transition of 4f electrons can be readily obtained by applying ambient pressure or changing the chemical composition [14,15]. Moreover, the observation of an evolving superconducting dome in the vicin-ity of the magnetic QCP strongly suggests that it is a candidate for AFM spin-fluctuation mediated superconductivity [16] [17,18]. The proximity between antiferromagnetism and d-wave superconductivity is also observed in many strongly correlated systems, such as hightemperature superconductors and layered organic superconductors.
In some unusual cases [19][20][21][22], superconducting phases in heavy fermions have been found in the absence of an obvious nearby QCP. Then the composite fermion pairing theory [23][24][25] might be more appropriate than the spin-fluctuation hypothesis to describe the underlying physics. Nevertheless, there are cases where, despite the absence of an obvious QCP, spin-fluctuation mediated pairing might be important, as exemplified by the superconducting phase of CeIrIn 5 . [26,27] This is quite different from the "SC2" phase of CeCu 2 Si 2 [19,28] or the Y b heavily doped CeCoIn 5 [21].
On the theoretical side, the two-dimensional periodic Anderson model (PAM), the Kondo-lattice model (KLM) or the degenerate Coqblin-Schrieffer model [29] are expected to capture the essential physics of spin-fluctuation mediated Ce-115 superconductivity. Previous analytical studies include large-N approches [1], mean-field theory [30] and phenomenological models of fermions coupled to fluctuating bose modes [31]. The PAM/KLM models have also been treated with the variational method [32], exact diagonalization (ED) and DMRG calculations on small clusters [33]. In these studies, an Heisenberg exchange is usually artificially added to simulate the RKKY interaction.
Without artificially adding an Heisenberg exchange term, here we show, using large scale determinant quantum Monte Carlo (DQMC) simulations [34] as well as the dynamical cluster approximation (DCA) [35] with a continuous time quantum Monte Carlo (CTQMC) [36] impurity solver, [37] that d x 2 −y 2 superconductivity can arise out of an incoherent metallic phase in the frustrated PAM. We demonstrate that the width of the supercon-ducting dome surrounding the QCP can be increased by strong frustration. Based on the analysis of the magnetic susceptibility and of the anomalous self-energy, we also show that the driving force for pairing in this model comes primarily from retarded antiferromagnetic spin fluctuations, which reinforces the hypothesis that this is the case for Ce-115.

II. PERIODIC ANDERSON MODEL WITH FRUSTRATED HYBRIDIZATION
In a model description of heavy fermion systems, frustration, Kondo coupling strength and the f -orbital degeneracy determine the quantum phase of the ground state [6]. We thus consider the frustrated periodic Anderson model on a two-dimensional square lattice. [38] The Hamiltonian is where we take k, σ, i as the momentum, spin and lattice site indices respectively. n σ f denotes the occupation number operator for f -orbitals. The conduction band dispersion relation is chosen as k = −2t[cos(k x ) + cos(k y )], with the nearest-neigbour hopping integral t taken as the energy unit throughout this paper. We assume no dispersion of f -orbitals and we neglect the f -orbital degeneracy. The f energy level f is set to zero so that the filling is n f ∼ 1, n c ∼ 0.9. The strength of the Kondo coupling and of the RKKY interaction is determined by the combined effects of the f − c hybridization V k = V + 2V [cos(k x ) + cos(k y )] and of U , the Coulomb repulsion between f electrons. In an antiferromagnetic configuration of the conduction electrons, the on-site hybridization V and the hybridization with nearest-neighbors V lead to conflicting effective interactions with the f electron. The frustration is maximal when V and V are of the same order of magnitude, as in Ce-115. [10] Further discussion of the model appears in Appendix A.

III. RESULTS
After we present evidence for d x 2 −y 2 pairing, we will discuss the question of coherence. We then present the phase diagram and conclude with the origin of pairing. To reveal the many-body correlation effect on superconductivity in an unbiased way, we show the result of DQMC calculations for the pairing susceptibility P , de- fined by where the greek indices represent either conduction c or f electron operators, g(r) is the form factor in real space and G is the normalization factor, G = r |g(r)| 2 . Let P 0 be the bubble contribution, dropping the vertex contribution. For a given pairing channel, the sign of P − P 0 basically reflects whether the pairing is favoured (positive sign) or not (minus sign) [39]. Fig. 1 displays our DQMC results for d x 2 −y 2 , s and extended s wave pairing susceptibilities as a function of temperature T for V = 0.3, V = 0.75 and U = 4. We learn that among those various pairing channels, the d x 2 −y 2 wave dominates [10] since the effective pairing interaction P − P 0 increases rapidly as T is lowered (inset of Fig. 1), suggesting that a divergent susceptibility would occur at the Berezinsky -Kosterlitz -Thouless (BKT) transition temperature T BKT . [40]. This contrasts with local s-wave pairing where the effective pairing interaction P − P 0 decreases, becoming more negative as T is lowered. This rules out s-wave pairing.
[41] By varying V and V , the itinerant character of felectrons can be adjusted. This is done in the inset of Fig. 1 where DQMC results suggest that the pairing strength has a maximum for intermediate V and V , namely for V = 0.3 and V =0.75, in the vicinity of an antiferromagnetic phase transition, as we will see.
Further insight into the nature of pairing and on the phase diagram is provided by DCA [35,42], which allows us to reach much lower T than DQMC. Following Ref. 43, the pairing susceptibilities are obtained from small pinning fields in the linear response region. The results are shown in Fig. 2 for temperatures as low as T ∼ 0.015, where the inset of Fig. 2a shows that T c can be extrapolated from the diverging susceptibilities. Since in Ce−115 materials the Coulomb repulsion U of 4f orbitals is expected to be significantly larger than the LDA bandwidth, (U ∼ 5eV ) vs (W ∼ 400meV ), [44] we study the evolution of pairing upon approaching the extended Kondo limit of the PAM model. More specifically, in Fig. 2a we display the d x 2 −y 2 pairing susceptibility in the (f f, f f ) channel, keeping V 2 /U constant, V 2 /U ≈ 0.0225, and increasing U [45]. The behavior as a function of T clearly suggests that d x 2 −y 2 pairing is stable in the extended Kondo limit. In fact, the estimated T c grows with increasing U , despite the fact that at large T the susceptibility is suppressed with increasing U . At large T , where vertex corrections are unimportant, the on-site U reduces the low-energy density of states (DOS), decreasing the pairing susceptibility. By contrast, at low T , increasing U drastically diminishes the higher order frustrated exchange terms, therefore enhancing the magnetic interaction vertices that mediate Cooper pairing. Since a heavy fermion is composed of both f and c electrons, a † kσ = u k c † kσ + v k f † kσ , the heavy fermion Cooper pair a † k,↑ a † −k,↓ has four different d x 2 −y 2 susceptibilities that should diverge simultaneously at T c . This is shown in Fig. 2b. Since the RKKY coupling between the local moments of the f -orbitals is emergent in the PAM, the magnetic "pairing glue" originates from the f -orbitals. This is consistent with Fig. 2b. Note that at large T , when Kondo screening is weak, the (f c, f f ) and (cc, f f ) channels are strongly suppressed and the effective pairings P −P 0 in the (cc, cc) and (f c, f c) channels, shown in the inset, are small compared with the (f f, f f ) channel in the inset of Fig. 1. This differs from the prediction of the composite pairing theory [23] for a two-channel Kondo model. That model has a different source of pairing force, leading to a f c pairing channel that is the strongest.

B. Coherence
We stress that complete screening of the local moments is not a necessary condition for diverging pairing susceptibilities. In fact, we find that complete screening of f moments does not occur even when the QCP is approached. This is confirmed by the fact that the magnetic susceptibility increases with decreasing T (not shown) meaning that in our case the Fermi surface is small at the antiferromagnetic to paramagnetic transition. [46] This is further illustrated by DOS obtained from maximum entropy analytical continuation [47] of the local Matsubara Green functions in Fig. 3. The noninteracting result is repeated for reference on all panels. As shown on the left-top panel of Fig. 3, at large T the heavy fermion quasiparticle is absent, and because of the intense scattering by f local moments, the effective hybridization between the conduction band and f -band is reduced, leading to a strong suppression of the hybridization gap and to a larger DOS for the conduction band (red line) at the Fermi level. This can be understood in the scenario of f -orbital selective Mott transition (OMST) [48,49]. The d -wave superconducting gap that appears on the emergent quasiparticle peak at low T is displayed on the lower-right panel. On the leftbottom panel the superconducting order parameter is suppressed. One then sees that in the underlying normal state, the DOS of the conduction band develops a peak at the Fermi level, reflecting the increased damping of the low-energy quasiparticles on the corresponding f -orbital that appears when superconductivity is suppressed.

C. Phase diagram, QCP and frustration
To mimic the Doniach phase diagram, where the ratio of Kondo and RKKY couplings is the control parameter, we plot in Fig. 4 the phase boundaries in the T − V plane for two fixed values of the frustration, V /V = 2 and V /V = 5. In the limit of small hybridization, the Kondo screening and antiferromagnetic RKKY correlations in PAM/KLM both vanish, leaving only local moment fluctuations in the orbital-selective Mott insulator with large entropy at non-zero temperatures [48]. As V increases, the Kondo energy scale and RKKY interaction both increase. First, RKKY dominates over Kondo screening and the antiferromagnetic ground state appears in coexistence with superconductivity. Increasing V again drives the system across the QCP. Then, long-range RKKY correlations are gradually quenched by Kondo screening and, eventually, superconductivity is destroyed when f electrons become too itinerant.
It is striking that the superconducting dome follows the change in the position of QCP with changing magnetic frustration (V /V ratio). This result vividly illustrates the intrinsic connection between the QCP and superconductivity in the PAM model. We also point out that the superconducting dome widens in parameter space when the system is more frustrated (V /V =2). This phenomenon might explain the appearance of a superconducting phase in CeIrIn 5 despite the absence of a nearby magnetic QCP. We also note that the right-end of the superconducting dome does not move toward much larger V even when the QCP does. This can be understood as indicating that in this region the f electrons become more itinerant so that short-ranged antiferromagnetic fluctuations are suppressed, even when there is a nearby QCP.

D. Retardation and origin of pairing
The origin of pairing may be deduced from the dynamical processes entering the zero-frequency real part of the anomalous self-energy I Σ (ω) [50] and from the cumulative order parameter I G (ω) [51], defined respectively by where Σ a (ω) is the imaginary part of the anomalous self-energy Σ i,i+r a (ω) and F (ω) is the imaginary part of the retarded lattice Gorkov function F (ω) = Im β 0 dτ c i+r (τ )c i (0) e i(ω+iη)τ , with i and i + r nearest neighbours.
These functions are plotted in Fig.5 along with the imaginary part of the antiferromagnetic spin susceptibility χ (q = (π, π), ω). As in the case of cuprates [50,51], we find the that the Cooper pairs initially form over an energy range comparable to those where antiferromagnetic fluctuations develop. The dependence of pairing on the RKKY interaction strength can clearly be seen in Fig. 5, where increasing U , hence decreasing the RKKY interaction, shifts the peaks of both I Σ (ω) and χ (q = (π, π), ω) towards the low-energy side. Frustrating magnetism in the conduction band by adding a nextnearest neighbour hopping t also shows the same correlation between the two quantities: the characteristic frequency of the spin fluctuations decreases along with characteristic frequencies in both I Σ (ω) and I G (ω). Our results thus confirm that the retarded spin-fluctuations mediate d -wave superconductivity in heavy fermion superconductors.
In addition to the contribution of low-frequency (retarded) spin-fluctuations, Fig.5b shows that there can be a significant gain in pairing for an energy scale set by the upper Hubbard band of the f electrons. This highfrequency (more instantaneous) contribution to pairing is much larger than what is found in the case of cuprates. [51]. It is probably because the upper Hubbard band seen from the f electron point of view is still in the conduction band. In other words, it is enhanced by the large RKKY interaction [50] that results from intermediate U and moderate conduction band frustration. In realistic Ce−115 materials the RKKY interaction is believed to be small [15], hence the high-frequency contribution to pairing should be less important.

IV. CONCLUSION
In summary, we investigated the 2D PAM model with frustration induced by the real-space structure of the hybridization V, V between f and c electrons. The unbiased DQMC large-cluster simulations suggest that in this model, which has a Fermi surface resembling the α sheet of Ce−115 materials [52], d x 2 −y 2 pair correlations increase rapidly at low T , and are enhanced when the anti-ferromagnetic QCP is approached. The DCA-CTQMC calculations further confirm the existence of a d x 2 −y 2 superconducting phase that appears out of an incoherent metallic phase. In the T − V plane, the superconducting phase has a dome shape surrounding the QCP of the antiferromagnetic phase. Finally, through an analysis of the frequency dependence of pairing, we show that the d-wave superconductivity in this model is mediated by retarded spin fluctuations.

ACKNOWLEDGMENTS
We would like to thank A. Georges, M. Aichhorn, S. Burdin, P. Coleman, C. Lacroix, J. Schmalian, C. Pépin Green functions. The Lanczos solver at zero temperature on a 2 × 2 cluster with 7-8 bath sites is used in the section where the real-frequency functions I Σ (ω) and I G (ω) are computed. Cellular DMFT (CDMFT) on 2 × 2 cluster has also been carried out to make comparisons to our DCA result. Qualitatively consistency is obtained, though the superconducting transition T c in CDMFT is lower than that of DCA. This may be because DCA preserves the periodic boundary conditions while CDMFT breaks it. In order to calculate the pairing susceptibilities, we have used the pinning field approach, [43] i.e., we observe the response of the system as small pairing fields are applied. To make sure that the response resides in the linear region, we used pinning fields for three different values, 5 × 10 −5 , 2 × 10 −4 and 1 × 10 −3 to monitor the changes of the pair response.