Superconductivity Emerging from Excitonic Mott insulator
— Theory of Alkaline Doped Fullerene
A three-orbital model derived from the two-dimensional projection of the ab initio Hamiltonian for alkaline doped fullerene AC with A=Cs,Rb,K is studied by a variational Monte Carlo method. We correctly reproduce the experimental isotropic -wave superconductivity around the ab initio parameters. With narrowing the bandwidth, the transition to an insulator is also reproduced, where orbital symmetry is found to be spontaneously broken with emergence of an excitonic Mott insulator for two orbitals and an antiferromagnetic insulator nearly degenerate with a spin liquid for the third orbital. The superconductivity is a consequence of exciton melting.
pacs:74.20.Pq, 74.70.Wz, 75.10.Kt
Introduction.–High critical temperature (high-) superconductivity induced by strong electronic correlations is one of the central topics in condensed matter physics and it has been intensively studied since its discovery in a family of cuprate compounds Bednorz and Müller (1986); Imada et al. (1998); Lee et al. (2006), where crucial role of strong electron correlations are prominent in various aspects. Recent discovery of the high- superconductivity in the iron-based superconductors Kamihara et al. (2008), where the five iron 3 orbitals contribute to the low-energy degrees of freedom, renewed the interest of the multi-orbital effects on the high- superconductivity. Multi-orbital physics such as Hund’s physics as well as the orbital differentiation was proposed to play essential roles for understanding the normal state properties Misawa et al. (2012); Georges et al. (2013); Haule et al. (2008) and superconductivity Misawa and Imada (2014a).
In the alkali-doped fullerides C (A=K, Rb, Cs), low-energy physics is also described by degenerate orbitals — each alkaline atom donates one electron to the three-fold degenerate lowest unoccupied molecular orbital with the symmetry at the fullerene atom. Thus the three orbitals become half filled on average.
The unique feature of the solid C in comparison to the above cuprates and the iron-based superconductors is a small band width ( eV) ascribed to small overlap of molecular orbitals each at the neighboring fullerene atoms, while the largest Jahn-Teller phonon frequency coupled to the orbitals is comparable and as large as 0.2 eV. It makes the role of phonons substantially large. The discovered superconductivity with the highest K is believed to have the isotropic -wave symmetry supporting the crucial role of phonons Gunnarsson (1997); Ganin et al. (2008); Takabayashi et al. (2009); Ganin et al. (2010); Zadik et al. (2015); Mitrano et al. (2016).
The superconductivity is, however, found next to the Mott insulator, suggesting the role of electron correlation as well. In addition, the superconductivity is found when the three orbitals are half-filled, in sharp contrast to the cuprates and iron-based superconductors, where typically the carrier doping into half filled bands is required for the superconductivity.
It was proposed that the relatively strong Jahn-Teller electron-phonon coupling in the alkali-doped fullerides leads to an effectively negative (inverted) Hund’s rule coupling (IHRC) Capone et al. (2002, 2009). By analyzing a degenerate three-orbital model with the IHRC with the dynamical mean-field theory (DMFT), Capone et al. showed that the IHRC actually induces the isotropic -wave high- superconductivity Capone et al. (2002). Recently, quantitative evaluation of the interaction parameters including the IHRC from ab initio calculations has been done Nomura et al. (2015). By solving the obtained ab initio Hamiltonian with the extended DMFT (E-DMFT), calculated critical temperatures were shown to be consistent with the experimental results.
In general, superconductors with relatively high in the ratio to the energy scale of electron band width is found exclusively in strongly correlated systems and their small coherence lengths require serious account of spatial correlation effects, while DMFT does not consider them. Furthermore, severe competitions with magnetic and charge orders known in the cuprates and the iron-based superconductors urge us to seriously examine the spatial correlation and fluctuation effects.
In this Letter, we employ the many-variable variational Monte Carlo (mVMC) method Tahara and Imada (2008); Misawa and Imada (2014b); mVM () to take account of both spatial and quantum fluctuations. We first analyze the two-dimensionally projected ab initio Hamiltonian for the face-centered-cubic alkali-doped fullerides (fcc-CsC). Then, to capture the essence of the superconductivity, we also analyze simplified three- and two-orbital models on triangular and square lattices, respectively, in which each orbital has transfer in only one direction and forms one-dimensional chains in each triangular bond direction and they interact with each other only at the same site.
We reveal the conditions for enhancing superconductivity: First a realistic unique three-orbital structure, where the electronic transfer at each nearest-neighbor atomic bond is governed by only one orbital-diagonal transfer, depending on the bond direction, is important to reproduce the -wave superconductivity in the realistic parameter region. The superconductivity is replaced by the excitonic Mott insulator with antiferromagnetic order in the strong coupling region near the ab initio parameters and replaced by excitonic insulator without antiferromagnetism in the two-orbital model. The superconductivity is universally stabilized through the melting of the excitonic insulators realized by the carrier doping. The excitonic insulator is the mother state of the superconductivity in the multi-orbital models with the IHRC and offers the clue for obtaining higher- superconductivity in fullerides and/or other multi-orbital systems with the IHRC.
Model and Method.–We first study three-orbital low-energy Hamiltonian for triangular-lattice cross section of fcc-fullerides, which is defined by
where () is the creation (annihilation) operator on the th orbital with spin at th site and and are the number operator.
The transfer integrals in the fullerides are given by 33 matrices Nomura et al. (2012). Although the Hamiltonian is defined on the fcc lattice, to reduce the numerical cost, we project fcc lattice into the triangular lattice [see Fig. 1(a)], i.e., we only consider [1,0,1], [-1,1,0], [0,1,1] vectors in the fcc lattice and regard them as [1,0], [0,1], [1,1] vectors in the triangular lattice. (The triangular lattice is represented by mapping to the square lattice with diagonal bonds in one direction.) In this study, we take only nearest-neighbor transfer integrals (in the original triangular lattice representation) because the amplitudes of the next-nearest-neighbor transfer integrals are small. The matrices of the hopping integrals are given as
In this study, we use the ab initio parameters for fcc-CsC, which shows maximum , denoted by fcc-Cs() in the literature Nomura et al. (2012). In fcc-Cs(), is 0.0372 eV and we take this value as energy unit. Other hopping parameters are given as , , .
According to the calculations Nomura et al. (2015), interaction parameters that is renormalized by the electron phonon interactions are estimated as , , and with . They nearly satisfies the rotational symmetry, i.e., . In our calculations, we systematically monitor the interaction parameters with the constraint around the above realistic parameter region. We take sites with periodic-periodic boundary conditions, where () denotes number of sites (the number of orbitals). We define the electron density as .
In the mVMC method Tahara and Imada (2008); Misawa and Imada (2014b); mVM (), the variational wave function is defined as , where and are the Gutzwiller factors Gutzwiller (1963) defined as and the Jastrow factors Jastrow (1955); Capello et al. (2005), defined as , respectively. The pair-product part is the generalized pairing wave function defined as , where denotes the variational parameters. For details of wave functions, see Refs. Gros, 1989; Bajdich et al., 2008; Tahara and Imada, 2008; mVM, . In this Letter, we have independent variational parameters for pair-product part. All the variational parameters are simultaneously optimized by using the stochastic reconfiguration method Sorella (2001); Tahara and Imada (2008). The variational function can flexibly describe several phases such as correlated paramagnetic metals, antiferromagnetic/charge-ordered phases and superconducting phases as well as their coexistence.
Results.– First, we examine the stability of the superconducting phase at half filling in the two-dimensionally projected Hamiltonian (we call this model full model). To detect the superconductivity, we calculate the isotropic -wave superconducting correlations, which is defined as , where . We also calculate the average value of at long distance (), which is denoted by . As shown in Fig. 1(b), we find that the superconducting correlations approach a nonzero constant value at long distance for . The size dependence of the saturated value is small and this indicates that the superconducting phase is stable for . As we show later, insulating phase competes with superconducting phase and superconductor-insulator transition occurs around .
Next, we systematically study the interaction dependence. As shown in Fig. 1(c), the superconducting correlations are largely enhanced by increasing in the moderately strong coupling region (). This result indicates that the superconductivity is induced by the Suhl-Kondo mechanism Suhl et al. (1959); Kondo (1963). The enhancement is mainly induced by (see SM (), S.1).
To analyze the nature in the strong coupling region, the doping dependence of the chemical potential defined by , which is shown in Fig. 1(d). Here, is the total energy at filling and . The insulating phase indicated by the gap in appears in the strong coupling region.
In Fig. 2(a), we show dependence of energies for the superconducting and insulating phases. The energy crossing indicates the first-order phase transition at . Because the ab initio interaction parameter is estimated as Nomura et al. (2015), the insulating phase looks overestimated on the triangular lattice. However, this is naturally understood from the bandwidth () of the triangular lattice, which is a half of the fcc lattice. Since the ratio primarily determines the metal (superconductor)-insulator transition, the transition found at here corresponds to for the fcc lattice. Therefore the ab initio value is close to the transition but in the superconducting phase in agreement with the experimental result. On the other hand, we note that IHRC is required to keep a value comparable to the ab initio value to stabilize the superconductivity even for the present model on the two-dimensional triangular lattice.
To further examine the nature of the insulating phase, we show the orbital occupancies and orbital-dependent doublon densities as a function of in Fig. 2(b), which are defined as and . shows orbital differentiation and thus the orbital symmetry breaking in the insulating phase.
In this orbital-differentiated Mott insulator, two orbitals are disordered and the antiferromagnetic order exists only in one orbital with the small doublon density [see the inset in Fig. 2(a) and S. 3 in SM ()]. Because orbital has larger transfer in direction () compared with orbital, short range antiferromagnetic correlation is induced by the proximity effect of the antiferromagnetic order in orbital [see S. 4 in SM ()]. This short-range antiferromagnetic correlation induces the small but finite differences in between the disordered orbitals ( and ).
We further note that a genuine orbital-differentiated Mott insulator without any spin symmetry breaking exists as a low-energy excited state and it’s energy relative to the antiferromagnetic ground state is very small () (see SM (), S.5). This near degeneracy of the magnetically ordered and disordered states is consistent with the experimental results where the antiferromagnetic transition occurs at a very low temperature for fcc-CsC (K) Ganin et al. (2010); its ordered moment is considerably small and it was suggested that antiferromagnetic phase coexists with the magnetically disordered phase Kasahara et al. (2014). We will discuss the nature of the orbital-differentiated Mott insulator later.
Next, to capture the essence of the superconductivity and the insulating phase found in the full model, we study a simplified model. One characteristic feature of the transfer matrices in the full model is that the directions of the largest hopping depend on the orbitals. If we consider only the largest hopping , the three orbitals form three chains as shown in Fig. 1(a). In this model, if one maps the triangular lattice to the square lattice with the diagonal hopping in one direction , each orbital has hopping only in either , , or  directions, depending on the orbitals. Electrons on different chains interact only at the same site through , , and terms. We call this simplified model the model.
To reveal the multi-orbital effects on the superconductivity, we also consider a two-orbital model on the square lattice analogously to the model. In the two-orbital model, the orbital 0 (1) has hopping only in the  () direction. The and two orbital models have some resemblance to the Kitaev model Kitaev (2006) in the sense that they have orbital-dependent spatial anisotropy.
In the case of the two-orbital model, insulating charge gap opens for any positive as shown in an example in Fig. 3(a). The nature of the insulator is understood as a doublon and a holon locally bound and resonating in the form of Frenkel excitons at each site with the degenerate two orbitals, while they do not break any translational symmetry, where charge/spin structure factors do not have appreciable peaks (see SM (), S.6). This local configuration indeed does not have energy loss of the interorbital Coulomb interaction . Therefore, the ground state of the two-orbital model is interpreted as an excitonic insulator phase without any spatial symmetry breaking and stabilized by electron (originally including electron-phonon) correlations where the orbitals of doubly occupied (doublon) and empty (holon) orbitals do not spatially order, but is represented by the linear combination at each site, where for instance expresses the doublon (holon) at the orbital 1 (2).
By doping carriers into the excitonic (bosonic) Mott insulator, the superconducting phase immediately appears as shown in Fig. 3(b), where doublons (local Cooper pair) condense. Although the local Cooper pair is already formed at half filling, it is frozen.
Discussion and Summary.– In the three-orbital models (full and models), because the three electrons occupy a site at half filling at large and , an exciton (one doublon and one holon) formation leaves one electron unpaired at the third orbital. This introduces much larger charge fluctuations than the two-orbital model. We find that this charge fluctuation easily causes charge melting of the excitonic Mott insulator phase found in the two-orbital model and the superconductivity appears even at half filling.
If the local repulsions and become sufficiently strong, in the full model, one unpaired orbital finally becomes the antiferromagnetic Mott insulator and other orbitals with Frenkel exciton forms an excitonic insulator similarly to the two orbital model. This phase may be called an orbital-differentiated Mott insulator. A similar phase in a three-orbital model was discussed before Hoshino and Werner (2017).
In the three-orbital model, a disordered state exists as a low-energy excited state in the strong coupling region, where spin order is absent. Since the averaged electron filling per unit cell is odd (three), this spin-disordered state is a genuine Mott insulator, which is not adiabatically connected to the band insulator. In other words, this phase is a promising candidate of quantum spin liquid if its energy can be lowered below the antiferromagnetic state and further detailed study of tuning is an intriguing issue but left for future studies.
AC forms a fcc or A15 structure instead of the 2D lattice considered here. We expect that the essential physics is the same because each bond has basically only one dominant diagonal hopping and the same excitonic Mott insulator is likely located nearby. Detailed study on the real lattice structure is an important future issue.
In summary, for the two-dimensional version of the Hamiltonian (full model) for the alkaline doped fullerene, we have found that the isotropic -wave superconducting phase becomes the ground state in the realistic parameter region. In stronger coupling region, a superconductor-insulator phase transition occurs, where the insulating phase is interpreted as an orbital differentiated Mott insulator consisting of a two-orbital excitonic Mott insulator and a one-orbital antiferromagnetic insulator. The antiferromagnetic insulator in the third orbital is nearly degenerate with the spin-liquid phase. To extract the essence, we have also analyzed simple three- and two-orbital models, consisting of network of chains with orbital-dependent anisotropic transfer. The comparison of the two- and three-orbital models suggests that the melting of the excitonic insulator by the carrier doping or by introducing the third frustrating orbital stabilize isotropic -wave superconductivity. Although the origin of the superconductivity is ascribed to the strong attractive force arising from the IHRC, the attractive force generally induces the insulating phase at half filling. To melt the insulator at half filling, the odd number of orbitals is helpful. Our study on the multi-orbital systems with the IHRC shows a new route to obtain stable superconductors as well as a quantum spin liquid in the systems with the IHRCs.
This work was financially supported by a Grant-in-Aid for Scientific Research (No. 22104010, No. 16H06345 and No. 16K17746) from Ministry of Education, Culture, Sports, Science and Technology, Japan. This work was also supported in part by MEXT as a social and scientific priority issue (Creation of new functional devices and high-performance materials to support next-generation industries GCDMSI) to be tackled by using post-K computer. The authors thank the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo for the facilities. We thank the computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (hp150173, hp150211, hp160201,hp170263) supported by Ministry of Education, Culture, Sports, Science, and Technology, Japan. TM is supported by Building of Consortia for the Development of Human Resources in Science and Technology from the MEXT of Japan.
- Bednorz and Müller (1986) J. G. Bednorz and K. A. Müller, Z. Phys. 64, 189 (1986).
- Imada et al. (1998) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
- Lee et al. (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
- Kamihara et al. (2008) Y. Kamihara, T. Watanabe, M. Hirano, and H. Hosono, J. Am. Chem. Soc. 130, 3296 (2008).
- Misawa et al. (2012) T. Misawa, K. Nakamura, and M. Imada, Phys. Rev. Lett. 108, 177007 (2012).
- Georges et al. (2013) A. Georges, L. d. Medici, and J. Mravlje, Annu. Rev. Condens. Matter Phys. 4, 137 (2013).
- Haule et al. (2008) K. Haule, J. H. Shim, and G. Kotliar, Phys. Rev. Lett. 100, 226402 (2008).
- Misawa and Imada (2014a) T. Misawa and M. Imada, Nat. Commun. 5, 5738 (2014a).
- Gunnarsson (1997) O. Gunnarsson, Rev. Mod. Phys. 69, 575 (1997).
- Ganin et al. (2008) A. Y. Ganin, Y. Takabayashi, Y. Z. Khimyak, S. Margadonna, A. Tamai, M. J. Rosseinsky, and K. Prassides, Nat. Mat. 7, 367 (2008).
- Takabayashi et al. (2009) Y. Takabayashi, A. Y. Ganin, P. Jeglič, D. Arčon, T. Takano, Y. Iwasa, Y. Ohishi, M. Takata, N. Takeshita, K. Prassides, et al., Science 323, 1585 (2009).
- Ganin et al. (2010) A. Y. Ganin, Y. Takabayashi, P. Jeglič, D. Arčon, A. Potočnik, P. J. Baker, Y. Ohishi, M. T. McDonald, M. D. Tzirakis, A. McLennan, et al., Nature 466, 221 (2010).
- Zadik et al. (2015) R. H. Zadik, Y. Takabayashi, G. Klupp, R. H. Colman, A. Y. Ganin, A. Potočnik, P. Jeglič, D. Arčon, P. Matus, K. Kamarás, et al., Sci. Adv. 1, e1500059 (2015).
- Mitrano et al. (2016) M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M. Riccò, S. R. Clark, et al., Nature 530, 461 (2016).
- Capone et al. (2002) M. Capone, M. Fabrizio, C. Castellani, and E. Tosatti, Science 296, 2364 (2002).
- Capone et al. (2009) M. Capone, M. Fabrizio, C. Castellani, and E. Tosatti, Rev. Mod. Phys. 81, 943 (2009).
- Nomura et al. (2015) Y. Nomura, S. Sakai, M. Capone, and R. Arita, Sci. Adv. 1, e1500568 (2015).
- Tahara and Imada (2008) D. Tahara and M. Imada, J. Phys. Soc. Jpn. 77, 114701 (2008).
- Misawa and Imada (2014b) T. Misawa and M. Imada, Phys. Rev. B 90, 115137 (2014b).
- (20) http://ma.cms-initiative.jp/en/application-list/mvmc.
- (21) See Supplemental Materials for details of and dependence of the superconducting correlations, doping dependence of , size dependence of the spin structure factors, short-range antiferromagnetic correlations of disordered orbital, low-energy excited state, and spin and charge structure factors in the two-orbital model.
- (22) T. Ohgoe and M. Imada, arXiv:1703.08899.
- Nomura et al. (2012) Y. Nomura, K. Nakamura, and R. Arita, Phys. Rev. B 85, 155452 (2012).
- Gutzwiller (1963) M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
- Jastrow (1955) R. Jastrow, Phys. Rev. 98, 1479 (1955).
- Capello et al. (2005) M. Capello, F. Becca, M. Fabrizio, S. Sorella, and E. Tosatti, Phys. Rev. Lett. 94, 026406 (2005).
- Gros (1989) C. Gros, Ann. Phys. 189, 53 (1989).
- Bajdich et al. (2008) M. Bajdich, L. Mitas, L. K. Wagner, and K. E. Schmidt, Phys. Rev. B 77, 115112 (2008).
- Sorella (2001) S. Sorella, Phys. Rev. B 64, 024512 (2001).
- Suhl et al. (1959) H. Suhl, B. T. Matthias, and L. R. Walker, Phys. Rev. Lett. 3, 552 (1959).
- Kondo (1963) J. Kondo, Prog. Theor. Phys. 29, 1 (1963).
- Kasahara et al. (2014) Y. Kasahara, Y. Takeuchi, T. Itou, R. H. Zadik, Y. Takabayashi, A. Y. Ganin, D. Arčon, M. J. Rosseinsky, K. Prassides, and Y. Iwasa, Phys. Rev. B 90, 014413 (2014).
- Kitaev (2006) A. Kitaev, Ann. Phys. 321, 2 (2006).
- Hoshino and Werner (2017) S. Hoshino and P. Werner, Phys. Rev. Lett. 118, 177002 (2017).
- Koga and Werner (2015) A. Koga and P. Werner, Physical Review B 91, 085108 (2015).
Appendix A S.1 Pair hopping and exchange interaction dependence of the superconducting correlations
In the main text, we show that the superconducting correlations are mainly governed by the IHRC. The IHRC consists of two parts, i.e., the exchange Hund’s rule coupling , and the pair hopping term in Eq. (1). As shown in Fig.S 1, we show that superconducting correlations are largely suppressed by decreasing the amplitude of the pair hopping terms and it becomes nearly zero for . In contrast to that, does not alter the superconducting correlations. This result clearly shows that the pair hopping is the crucial interaction for the superconductivity. We note that the previous studies with DMFT for similar multi-orbital systems claims that superconductivity appears without pair hoppings Koga and Werner (2015); Hoshino and Werner (2017). In the light of our result, these results should be reexamined whether they are artifacts of the DMFT because the DMFT generally overestimate the superconductivity due to its mean-field nature.
Appendix B S.2 Doping dependence of the orbital-dependent doublon densities
We show the doping dependence of the orbital-dependent doublon densities in Fig.S 2. Around half filling, the spontaneous symmetry breaking of orbital still exists and it vanishes at large doping region, where the superconducting phase becomes stable. A first-order phase transition between the excitonic phase with orbital-differentiated doublons and the superconducting phase occurs at finite doping. The doping dependence of the chemical potential has kink structure as we mentioned in the main text.
Appendix C S.3 Size dependence of the spin structure factors in the full model
In Fig.S 3, we show the spin structure factors in the strong coupling region for several different systems sizes. For small system sizes (), we find that the spin structure factors have peak at while they have peak at in large system sizes (). Although it is hard to perform calculations for larger system sizes, it is plausible that stripe magnetic order () becomes stable in the thermodynamic limit.
The size extrapolation of the spin structure factors is shown in Fig.S 4. Although the positions of the Bragg peaks depend on the system sizes, the size dependence of the peak value is smooth. We estimate the thermodynamic value of the spin structure factors as .
Appendix D S.4 Short-range antiferromagnetic correlations of disorders orbitals in the full model
In Fig.S 5, we show how the antiferromagnetic order in orbital affects the spin structure factors on the other disordered orbitals and . Because orbital has relatively larger hopping in direction compared to orbital, orbital is largely affected by the antiferromagnetic order and shows short-range antiferromagnetic order as shown in Fig.S 5(a).
Appendix E S.5 Low-energy excited state in the full model
Here, we mention the low-energy excited state in the insulating phase of the full model. For , the low-energy excited state without spin/charge order. The spin and charge structure factors are shown in Fig.S 6. The charge structure factors is defined as
We do not find any significant peaks in the spin and charge structure factors. Its energy per site is (energy unit is ) while the ground state energy is . The energy difference is given by . This nearly degenerate state may make the Neel temperature low and induce the peculiar coexistence of the magnetic order and the disordered state observed in experiment Kasahara et al. (2014).