Superconductivity in ionic-Hubbard model on honeycomb lattice

Superconductivity in ionic-Hubbard model on honeycomb lattice

Physica C 484 (2013) 56–58 Contents lists available at SciVerse ScienceDirect Physica C journal homepage: www.elsevier.com/locate/physc Superconduc...

300KB Sizes 0 Downloads 37 Views

Physica C 484 (2013) 56–58

Contents lists available at SciVerse ScienceDirect

Physica C journal homepage: www.elsevier.com/locate/physc

Superconductivity in ionic-Hubbard model on honeycomb lattice T. Watanabe a,⇑, S. Ishihara b a b

Department of Natural Science, Chiba Institute of Technology, Shibazono 2-1-1, Narashino 275-0023, Japan Department of Physics, Tohoku University, Aramaki-aza-aoba 6-3, Aoba-ku, Sendai 980-8578, Japan

a r t i c l e

i n f o

Article history: Accepted 17 February 2012 Available online 27 February 2012 Keywords: Layered nitride superconductor Honeycomb lattice Ionic-Hubbard model Variational Monte Carlo method

a b s t r a c t To consider the origin of superconductivity (SC) in layered nitride superconductors, we study an ionicHubbard model on a honeycomb lattice, which consists of two sublattices with a level offset, using an optimization variational Monte Carlo method. The parameter values for the Coulomb interaction and the level offset, to realize a band insulator, is evaluated in the non-doped system. In carrier-doped band insulating systems, a spin singlet d-wave SC state is stabilized. The spin fluctuation remains relatively strong even for large level offsets. We find that a renormalization effect of the energy band primarily contributes to the spin fluctuation. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Experimental research on layered nitride superconductors, LixMNCl (M = Hf, Zr), which is composed of alternate stacking of MN and Cl bilayers, has revealed many interesting features in the superconducting (SC) state [1]. They are quasi two-dimensional superconductors, in which electron conduction occurs in the MN bilayer with a honeycomb lattice structure consisting of M and N [2–4]. The mother compound b-MNCl is a band insulator and shows SC by doping of electrons into the MN bilayer, namely by the Na and Li intercalation. They exhibit relatively high Tc, up to 15 K for M = Zr and 25 K for M = Hf, regardless of the very small specific heat coefficient 1 mJ/mol K2 [5]. Furthermore, the isotope effect has been found to be weak [6,7] and the SC gap ratio is 2D0/ kBTc = 4.6–5.2 [5], which is unexpectedly high in comparison with the BCS theory. These results show that the SC in LixMNCl is not of a conventional BCS type. NMR experiments suggest that the SC gap is of an unconventional spin singlet pairing [6,8]. Tunneling spectroscopy [9] and specific heat measurements [5,10] suggest that it is a fully opened but an anisotropic gap. If the pairing symmetry is of the d-wave type similar to the cuprate, the spin fluctuation is considered to be a strong candidate for an origin of the SC. However, the present SC state is not realized near a magnetic Mott insulator, but a nonmagnetic band insulator. In order to clarify the origin, the stability for SC close to the metal to band insulator transition should be investigated in an appropriate theory. To this end, we study an ionic-Hubbard model on a singlelayered honeycomb lattice, which consists of two sublattices with ⇑ Corresponding author. Tel.: +81 47 454 9619; fax: +81 47 454 9689. E-mail address: [email protected] (T. Watanabe). 0921-4534/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.physc.2012.02.021

a level offset. We adopt an optimization variational Monte Carlo (VMC) method, which enables us to accurately evaluate a metal to band-insulator transition as well as a Mott transition, despite the coupling strength. First, we determine the parameter values for the Coulomb interaction and the level offset to realize the band insulator, at a half-filling. We study a stability of SC and a spin correlation function for different doping into the band insulator. 2. Formulation As an effective model for LixMNCl, we consider a Hubbard model on a single layer honeycomb lattice consisting of alternating ‘‘M’’ and ‘‘N’’ orbitals with a level offset. The Hamiltonian is given as follows,

H ¼ t

X y X X X ðci;r ciþs;r þ h:c:Þ þ U ni" ni#  D ni þ D ni : i;s;r

i

i2A

ð1Þ

i2B

Here, s runs over all nearest-neighbor sites of the i site and ni = ni" + ni; is the number operator with ni;r ¼ cyi;r ci;r . In this model, the level offset of two orbitals is 2D, where D is the staggered potential between A and B sublattices. This model almost reproduces the conduction and valence bands closest to the Fermi level obtained by a first principles band calculation [11]. To this model, we apply an optimization VMC method [12]. As a variational wave function, we use a Jastrow type: W = PU, where U is an onebody (Hartree–Fock) part, and P is a many-body correlation factor. As Q for P, in addition to the onsite (Gutzwiller) factor, PG(g) = i[1  (1  g)ni"ni;], we introduce a crucial intersite factor PQ (P = PGPQ), which binds a doublon (doubly occupied site) to a holon (empty site) within nearest neighbor sites [13]. Regarding U, we consider two trial states: (i) a state in which antiferromagnetic (AF) and charge orders coexist, UAF+CO(DAF,

57

T. Watanabe, S. Ishihara / Physica C 484 (2013) 56–58

e Þ. UAF+CO is obtained by diagoDCO), and (ii) a BCS state UBCS ðDd ; D nalization of the Hartree–Fock Hamiltonian with mean fields for two long-range orders, DAF and DCO. Here, DAF controls the staggered spin density of an AF order, and DCO the charge density wave consisting of charge-rich A and charge-poor B sites [14]. The BCS P state is written as UBCS ¼ k;a¼ ðuak eyka" eyka# ÞNe =2 , in which Ne is a

the results about the d- and d + id waves because only the BCS states are stabilized relative to the normal state. The d + id-wave pairing is suggested by a theory for a honeycomb-lattice Hubbard model with D = 0 [15] and a triangular-lattice t–J model [16]. In Fig. 2a, a phase diagram for the most stable state of UBCS is depicted in the plane of D and d, where d is the doping concentration per site. The d- and d + id-wave states are realized near half-filling and the d-wave state is significantly stabilized up to the band insulating region, namely even for D > Dc. As for the d-wave state, in Fig. 2b, we depict a long-range average of a d-wave nearest-neighbor pair correlation function,

fixed electron number, and eykar is a quasi-particle operator obtained by diagonalization of the kinetic energy part of Eq. (1). qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e2  l uak ¼ aDðkÞ=½Eak þ ðEak Þ2 þ jDðkÞj2  with Eak ¼ a jAk j2 þ D 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 ik s ~ 2 , where Ak ¼ t and DðkÞ ¼ ReAk DðkÞ= jAk j þ D and l0 is se

Pd ðrÞ ¼

chemical potential for the non-interacting system. The variation e is optimized as a renormalization parameter of enparameter D

3. Results In Fig. 1a, we show a phase diagram for the most stable of UAF+CO at half-filling. In all insulating states, we have confirmed the quadratic behavior in q near the C point of the charge structure factor N(q), which indicates the existence of a charge gap (not shown), except for the Dirac point (D, U)=(0, 0). The boundary between the AF Mott and band insulating regions is determined by the optimized parameters DAF and DCO: an AF Mott insulator for DAF > 0 and arbitrary DCO, while a band insulator for DAF = 0 and DCO > 0. In Fig. 1b, order parameters for the long-range AF and charge orders, P P MSDW ¼ 2=N S i ð1Þi Szi and M CDW ¼ 1=N S i ð1Þi ðni  1Þ, respectively, are shown as a function of D at U/t = 7. MSDW has a sufficiently large value for small D and vanishes at Dc (=1.9t), which is a boundary between Mott and band insulating regions. On the other hand, MCDW gradually increases with D and becomes large magnitude in the band insulating region. Now we show the results for SC in the doped band insulators. Because the present model has a particle–hole symmetry, the results are independent of the kind of carriers. Based on the trigonal symmetry of the system, we checked the stabilities for different pairing symmetries: s, d, d + id, p and f waves. Here, we present

malized down to around half of D, even in the region of D > Dc. Thus, because the renormalization effect remarkably reduces the original band gap, the spin fluctuation and the d-wave SC correlation survive in the carrier-doped band insulating state. 4. Summary and discussions To elucidate the origin of SC in LixMNCl, we evaluated the stability of the SC state in an ionic-Hubbard model on a single-layered honeycomb lattice, using an optimization VMC method. We investigated the D–U parameter region in which the band insulator is realized in the non-doped system. The values of U and D for MNCl are estimated by comparing the band calculation, as D  0.83t and U  6t [10]. In Fig. 1a, these values are found to be fairly close to the boundary between Mott and band insulating regions, although it is located in the Mott insulating area. Therefore, we suppose that MNCl is a band insulator close to the phase boundary. When

Δc

(b) 1

(a) 10 8

ð2Þ

pffiffiffi where Dys ðRi Þ ¼ ðcyi" cyiþs# þ cyiþs" cyi# Þ= 2 and ^r1 and ^r2 are primitive translation vectors in a real space. The pair correlation in the region of D/t < 2.5, in which the d-wave state appears, is distinctly enhanced against that for D > 2.5t. This indicates that the d-wave SC state is stabilized in the carrier-doped band insulator. In Fig. 2c,we plot the maximum of the spin structure factor P SðqÞ ¼ 1=NS i;j hSzi Szj ieiqðRi Rj Þ , when the peak is located at qmax = (2p/3,2p/3). Although S(qmax) is rapidly reduced as D approaches Dc, it remains significant magnitude even in the region of D > Dc, namely up to D  3.5t. This is why the d-wave SC state is stabilized in the carrier-doped band insulator. The predominant spin fluctuation originates from the renormalization effect of bands due to the electron correlation. The renormalized band gap parame is depicted as a function of D in Fig. 3. The curve for the case eter D ~ is renorof U = 0 is given by a broken line. By the injection of U, D

ergy bands due to the electron correlation, independently of D given in Eq. (1). As for the gap function D(k) = DBCSzk, we study different pairing symmetries, e.g. zk = cos k1  cos k2 for a dx2y2 (d) wave and zk ¼ cos k1 þ ei2p=3 cos k2 þ ei4p=3 cos k3 for a d + id wave. DBCS is a parameter to be optimized and k1, k2 and k3 are primitive translation vectors in a wave space. e , are optimized by an All parameters, g, l, DAF, DCO, Dd and D optimization VMC scheme with 2–5  105 samples. The systems of Ns = 96–196 sites with periodic–antiperiodic boundary conditions are used.

0.8

AF Mott insulator

6

MSDW

M

U/ t

0.6

4

M CDW

0.4

band insulator

2 0

0 1 X X ð1Þ1dðs;s Þ hDys ðRi ÞDs0 ðRi þ rÞi; 4NS i s;s0 ¼^r ;^r 1 2

0.2

0

1

2

3

Δ /t

4

5

0

0

1

2

3

Δ /t

Fig. 1. (a) Phase diagram spanned by D and U for UAF+CO, and (b) MSDW and MCDW as a function of D for U/t = 7, at half-filling. The vertical arrow indicates Dc. For both (a) and (b), the system of Ns = 196 is used.

58

T. Watanabe, S. Ishihara / Physica C 484 (2013) 56–58

Δc

(a) 0.15

(c) 5

0.05

4

0

3

Δc

[

S (qmax)

δ

0.1

-4

10 ]

(b)

4

2

Pdave

3 1

2 1 0

0 0 0

1

2

1

3

2

3

4

Δ /t

Δ /t

Fig. 2. (a) Phase diagram in the D and d parameter space for UBCS, where a solid circle exhibits that only the d-wave state is stable against the normal state, and an empty triangle indicates that both d- and d + id-wave states are stable. (b) Long-range average of a d-wave pair correlation function P ave and (c) the maximum of the spin structure d factor S(qmax), as a function of D for U/t = 7 and d = 0.04. For all figures, the vertical arrow points Dc and the system of Ns = 96 is used.

d > 0.10. To pursue this phenomenon, other factors, such as intersite and inter-layer interactions, might be necessary.

Δc

Acknowledgment

4

Δ=Δ

The author appreciate Prof. Y. Iwasa, and Y. Taguchi for their valuable discussion. This work is supported by KAKENHI (23104703 and 21224008) from MEXT, and Grand Challenges in Next-Generation Integrated Nanoscience. A part of numerical calculations has been done in supercomputing systems in ISSP, the University of Tokyo, and Kyoto University.

~

Δ /t

3

2

References

1

0 0

1

2

3

4

Δ /t e as a function of D for U/t = 7 and Fig. 3. Renormalized band gap parameter D d = 0.04.A broken line is for the case of U = 0. The vertical arrow indicates Dc and the system of Ns = 96 is used.

carriers are doped into such a band insulator, the spin fluctuation survives with sufficient magnitude and the d-wave SC is stabilized near the half-filling. The primary origin is in a renormalization effect of the energy band. We think that the band gap estimated by a first principle calculation for MNCl is significantly reduced due to the electron correlation. The renormalized ban gap may favor SC mediated by spin fluctuation. However, within the present model, it seems not easy to understand the experimental observation of the doping concentration d-dependence of Tc: Tc is enhanced near the insulating phase, d  0.025, and is nearly d independent for

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]

[15] [16]

S. Yamanaka, K. Hotehama, H. Kawaji, Nature 392 (1992) 580. I. Hase, Y. Nishihara, Phys. Rev. B 60 (1999) 1573. C. Felser, R. Seshadri, J. Mater. Chem. 9 (1999) 459. R. Weht, A. Filipetti, W.E. Pickett, Europhys. Lett. 48 (1999) 320. Y. Taguchi, M. Hisakabe, Y. Iwasa, Phys. Rev. Lett. 94 (2005) 217002. H. Tou, Y. Maniwa, S. Yamanaka, Phys. Rev. B 67 (2003) 100509. R. Y. Taguchi, T. Kawabata, T. Takano, A. Kitora, K. Kato, M. Tanaka, Y. Iwasa, Phys. Rev. B 76 (2007) 064508. K. Hotehama, T. Koiwasaki, K. Uemoto, S. Yamanaka, H. Tou, J. Phys. Soc. Jpn. 79 (2010) 014707. (a) T. Ekino, T. Takasaki, H. Fujii, S. Yamanaka, Physica C 388–389 (2003) 573; (b) T. Ekino, T. Takasaki, H. Fujii, S. Yamanaka, Physica B 328 (2003) 23. Y. Kasahara, T. Kishiume, T. Takano, K. Kobayashi, E. Matsuoka, H. Onodera, K. Kuroki, Y. Taguchi, Y. Iwasa, Phys. Rev. Lett. 103 (2009) 077004. K. Kuroki, Phys. Rev. B 81 (2010) 104502. C.J. Umrigar, K.G. Wilson, J.W. Wilkins, Phys. Rev. Lett. 60 (1988) 1719. (a) H. Yokoyama, H. Shiba, J. Phys. Soc. Jpn. 56 (1987) 1490; (b) H. Yokoyama, Prog. Theor. Phys. 108 (2002) 59. In the Hartree-Fock approximation, the mean field densities are defined as follows: hnir i ¼ n þ rDAF þ DCO for i 2 A and hnir i ¼ n  rDAF  DCO for i 2 B, where n ¼ 1  d. S. Pathak, V.B. Shenoy, G. Baskaran, Phys. Rev. B 81 (2010) 085431. T. Watanabe, H. Yokoyama, Y. Tanaka, J. Inoue, M. Ogata, J. Phys. Soc. Jpn. 73 (2004) 3404.