Inhomogeneous Electronic Distribution in High-T Cuprates
We theoretically investigate the doping evolution of the electronic state of high-T cuprate on both sides of the half-filling on the basis of the three-dimensional three-band Hubbard model with a layered structure using the Hartree-Fock approximation. Once a small amount of holes or electrons are doped into the half-filled state, our model exhibits the charge-transfer insulator-to-metal transition along with a chemical potential jump. At the same time, the doped holes or electrons are inhomogeneously distributed, and they tend to form clusters in the vicinity of the half-filling. This suggests the possibility of microscopic phase separation with the separation between the metallic and the insulating regions.
As shown by the photoemission spectroscopy, the undoped high-T superconducting cuprate (HTSC) is a charge-transfer insulator, where the Cu 3d electrons are almost localized by the strong electron correlation.  When the rare-earth element is substituted out of the two-dimensional (2D) CuO layers and the number of holes or electrons doped into the CuO layers is increased, Cu 3d electrons hybridized with O 2p electrons achieve itinerancy and display superconductivity. Moreover, in slightly hole-doped LaSrCuO with , the chemical potential shift suppression is observed by photoemission spectroscopy (PES). [2, 3] The electronic phase separation between the antiferromagnetic insulating phase and the superconducting phase is considered to be one of the reason for the suppression of the shifting of the chemical potential.  The phase separation assumes that the doped carriers are inhomogeneously distributed due to the strong electron correlation.
Many experimental findings have suggested that the electrons under such circumstances favor some types of spontaneous ordering in certain doped regions. For instance, in order to elucidate the anomalous suppression of the superconducting transition temperature in LaBaCuO,  the spin and charge correlations in LaBaCuO or (La,Nd)SrCuO with have been intensively studied, and the stripe order has been observed by neutron scattering. [6, 7, 8, 9, 10, 11, 12] The stripe order has also been observed by x-ray scattering. [13, 14, 12, 15, 16, 17] Furthermore, an electron paramagnetic resonance study showed that microscopic electronic phase separation occurs in LaSrCuMnO with ,  and a Cu nuclear magnetic resonance (NMR) study suggested that a large charge droplet (’blob’) is formed in the electron-doped NdCeCuO. 
Much theoretical works has also been performed to study the behavior of the doped carriers in HTSC. Pioneering works adopting the Hartree-Fock approximation (HFA) have studied the stripe order in LaBaCuO on the basis of the 2D one-band Hubbard model [20, 21] or the 2D two-band Hubbard model.  Furthermore, dynamical mean field theory (DMFT) has been exploited to study the stripe phase on the basis of the 2D Hubbard model with non-equivalent sites, where .  The DMFT approach has also been adopted to analyze the three-band Hubbard model. [24, 25, 26, 27] In some of these works, [25, 27] the DMFT approach was combined with the local density approximation (LDA). Another study  considered the possibility of two-sublattice antiferromagnetism. All of these works have successfully reproduced the Zhang-Rice singlet band. [24, 25, 26, 27] This shows that the three-band Hubbard model is an appropriate model for HTSC near the half-filling and that the DMFT is a powerful tool for analyzing its electronic state. However, when we investigate the inhomogeneous electronic distribution near the half-filling, we need to adopt the model for a large number of non-equivalent sites. In general, the DMFT costs much more than the HFA to analyze the model with a large number of non-equivalent sites.
In this paper, we analyze the normal ground state of the 3D three-band Hubbard model with a single-layered perovskite structure in order to study the evolution of the electronic state when holes or electrons are doped into the undoped HTSC. We consider 256 non-equivalent copper sites for each rectangular parallelepiped super cell, and adopt the HFA for these conditions. We performed the calculation without any assumptions about the electronic distribution, and we obtained fully self-consistent solutions except near the 1/8-filling. These solutions showed the chemical potential jump at half-filling, which means that the electron suddenly becomes itinerant when a small number of holes or electrons are doped. Moreover, the doped holes or electrons tend to form clusters in the vicinity of the half-filling. These clusters are considered to form a metallic region, and are surrounded by the insulating region. This suggests the possibility of microscopic electronic phase separation in HTSC near half-fillig.
Our 3D three-band Hubbard model Hamiltonian, , is composed of d-electrons at each Cu site and p-electrons at each O site. To consider the spatial inhomogeneity, we introduce the rectangular parallelepiped super cell containing Cu and O sites as a unit cell. Thus, is defined as follows:
Here we use the abbreviations and , where and are the annihilation (creation) operators for the -orbital and -orbital electron on the -th site, as specified by the momentum and spin , respectively. , , and are the on-site Coulomb repulsion between the -orbitals, the number of -space lattice points in the first Brillouin zone (FBZ), and the chemical potential, respectively. The FBZ is defined in the reciprocal space to the lattice whose unit cell contains Cu and O sites. The unit cell is schematically shown in Fig. 1. The two non-equivalent CuO layers, indicated with and , are alternatively stacked along -axis. On each CuO layer, the size of the unit cell along -axis and -axis is and Cu sites, respectively. Thus, when we set and , . We take the primitive translation vectors for the unit cell to be , , and , where , , and are the unit vectors for the -axis, -axis, and -axis, respectively. in Eq. (1) is defined as follows:
where is the hybridization gap energy between the - and -orbitals. , where , is the coordinate of the -orbital electron on the -th site. is Kronecker’s delta, i.e., it is for and for . In the following, we take both and to be a unit of length and set . Then, we can represent
where is the transfer energy between a p-orbital and a p-orbital, is that between a d-orbital and a p-orbital, and is that between d-orbitals, respectively. In this study, is the unit of energy.
We adopt the HFA with respect to every Cu site and two spin states, and we only consider collinear spin states. Thus, we define
and we approximate Eq. (1) as follows:
for a given total number of doped holes .
3 Results and Discussion
In the numerical calculations, we divide the FBZ into equally-spaced rectangular parallelepiped. The parameter sets are selected as , , , , and . In order to obtain self-consistent solutions, we need to set initial values for fields , , and and carry out iterative calculation until all of these fields have sufficient accuracy. Here, we chose the initial values for the fields from the uniform-distributed random numbers. In this manner, we obtained the solutions for the doping region near half-filling; and , where . For these solutions, all of the obtained fields , , and have three digits of accuracy. Once they get three digits of accuracy, all of these fields rapidly converge. Thus, we can consider the solutions with these fields as fully self-consistent. We also tried to obtain the solutions for other doping region e.g. but failed. For these doping regions, the ground state accompanied by certain types of long-period superlattice structure would be stable. The fields corresponding to such ground state could be hardly obtained by our iterative calculation. We henceforth concentrate our discussion on the doping region near half-filling.
The doping dependence of the chemical potential for our fully self-consistent solutions is shown in Fig. 2. The chemical potential changes rapidly at half-filling, which is consistent with experimental results obtained by comprehensive PES studies on HTSC.  The magnitude of this change is about , and the value is in accordance with , which is equal to in our calculations. This fact can be explained by the doping dependence of the density of states (DOS), which is shown in Fig. 3. In this figure, we show the mean DOS over the solutions with similar doping and chemical potential as DOS for each doping state, and we indicate each doping state with the mean doping over these solutions. All these labels in Fig. 3 are summarized in Table 1 For instance, DOS for in Fig. 3 is the mean DOS over the seven solutions for which we found and as a mean with standard error of the mean of doping and chemical potential, respectively.
In the slightly electron-doped case, , the Fermi level crosses the upper Hubbard d-band, while in the slightly hole-doped case, , the Fermi level crosses the p-band. The upper d-level is about larger than the p-level when , due to the Hubbard splitting. In the energy range between the p-band and the upper Hubbard d-band, the density of states is almost zero. In the undoped case, , the Fermi level locates in this energy gap, and it causes the chemical potential to jump by at half-filling. The quantum Monte-Carlo calculation of the two-band Hubbard model in the limit of infinite dimensions gave the same result for the chemical potential jump as our calculation.  Thus, the chemical potential jump at half-filling is the characteristic behavior of the multi-band Hubbard model composed of both Cu 3d electrons and O 2p electrons independent of the dimensionality of the model.
The inhomogeneous distribution of every obtained field can be observed in our solutions. The doping dependence of the distribution of the obtained fields can be shown by histogram, as in Fig. 4. In this figure, we show the accumulated histograms over the solutions with similar doping and chemical potential as the histogram for each doping state, and we indicate each doping state with the mean doping over these solutions as well as in Fig. 3. At half-filling, where the ground state is insulating, almost all have the same value near and almost all and have the same value near . In the electron-doped case, with the increase of the doped electrons, only the peak for broadens and its center shifts higher. In contrast, in the hole-doped case, with increasing hole density, not only the peak for but also the peaks for and broaden and their centers shift lower. Except for their variances, the doping dependence of the average of is similar to that of obtained by the LDA+DMFT calculation, and the doping dependence of the average of or is similar to that of . 
The spacial distribution of the magnetic moment also becomes inhomogeneous by doping. In the same manner as in Fig. 4, the doping dependence of the distribution of the magnetic momenta magnitude can be shown by histogram in Fig. 5. At half-filling, one half of is exactly at and the other half is exactly at . This indicates that the d-electron spins are fully polarized. Both in the electron-doped case and in the hole-doped case, with an increase of the doped carrier, the two peaks for become broader and their centers shift toward zero. This means that the number of the fully-polarized d-electron spins decreases and that the number of the doubly-occupied Cu sites increases. We should note that the tails of in the electron-doped case extend more than those in the hole-doped case.
The difference between the electron-doped case and the hole-doped case is caused by the doped carriers being differently distributed in the unit cell. In order to show this distribution, we define the distribution function of the doped carriers as follows:
The doping dependences of for the electron-doped case and for the hole-doped cases are shown in Figs. 6 and 7, respectively. In these figures, the red spots indicate the location of doped holes and the blue ones indicate the location of doped electrons. In the electron-doped case, as shown in Fig. 6, the doped electrons form blobs even in the slightly electron-doped case, with . These blobs can be identified as the ones in NdCeCuO, suggested by the Cu NMR study.  On the other hand, as shown in Fig. 7, the doped holes stay within a single CuO cluster in the slightly hole-doped case, with . Hence it can be recognized that the cluster formed by the doped electrons is larger than the one formed by the doped holes. We can quantitatively clarify this tendency with the doping dependence of the number of the sites having more than certain value of excess carriers in Fig. 8 since such a number should increase with the size of the cluster formed by the doped carriers. The difference between the electron-doped case and the hole-doped case is explained in the following; In both cases, the doped carriers basically tend to form an extended cluster, since the strong on-site Coulomb repulsion hinders the double occupancy on each site, and instead align them next to each other to gain more kinetic energy. The study on Coulomb gas ordering in a 3D layered system by the Brownian dynamics approach supports our explanation.  In our lattice model, the size of such a cluster depends on the number of adjacent orbitals where the doped carriers are allowed to occupy. In the electron-doped case, the doped electrons have room to sit only on the Cu sites, because all O sites are fully filled by electrons. On the other hand, in the hole-doped case, the doped holes have room to sit on both Cu sites and O sites. Therefore, when the number of doped electrons is as many as the number of doped holes, the cluster formed by the doped electrons is larger than the one formed by the doped holes.
In this paper, we conducted HFA calculations for the 3D three-band Hubbard model with a single-layered perovskite structure considering a large number of non-equivalent sites. We obtained the fully self-consistent solutions both for the electron-doped and hole-doped cases at or near half-filling. Our solutions show the chemical potential jump at half-filling. The jump can be explained by the DOS dependence on doping, which is characteristic to the multi-band Hubbard model composed of both Cu 3d electrons and O 2p electrons independent of the dimensionality of the model. Inhomogeneous electronic distributions near half-filling are observed in our solutions. There is a remarkable difference in the inhomogeneous electronic distributions between the electron-doped and hole-doped cases. That is, the clusters formed by doped carriers extend more in the electron-doped cases than in the hole-doped cases. The difference between the electron-doped and hole-doped case is caused by the difference in the species of orbitals the electron and hole are allowed to occupy, which should be explained only on the basis of the multi-band Hubbard model. Thus, the theoretical approach on the basis of the multi-band Hubbard model can explain both the chemical potential jump at half-filling and inhomogeneous electronic distributions near half-filling in a comprehensive way.
The authors are grateful to Dr. N. Nakai and Prof. Y. Aiura for their stimulating discussions. We are also grateful to an anonymous reviewer of the first manuscript for providing insightful comments and directions which have resulted in the revised manuscript.
-  A. Fujimori, E. Takayama-Muromachi, Y. Uchida, and B. Okai, Phys. Rev. B 35, 8814 (1987).
-  A. Ino, T. Mizokawa, A. Fujimori, K. Tamasaku, S. Uchida, T. Kimura, T. Sasagawa, and K. Kishio, Phys. Rev. Lett. 79, 2101 (1997).
-  A. Fujimori, A. Ino, T. Mizokawa, C. Kim, Z.-X. Shen, T. Sasagawa, T. Kimura, K. Kishio, M. Takaba, K. Tamasaku, H. Eisaki, and S. Uchida, J. Phys. Chem. Solids 59, 1892 (1998).
-  A. Fujimori, A. Ino, T. Yoshida, T. Mizokawa, Z.-X. Shen, C. Kim, T. Kakeshita, H. Eisaki, and S. Uchida, in Open Problems in Strongly Correlated Electron Systems, eds. J. Bonča, P. Prelovšek, A. Ramšak, and S. Sarkar (Kluwer Academic Publishers, Dordrecht, 2001) Chap. 1.
-  A. R. Moodenbaugh, Y. Xu, M. Suenaga, Y. J. Folkerts, and R. N. Shelton, Phys. Rev. B 38, 4596 (1988).
-  J. M. Tranquada, B. J. Sternlieb, J. D. Axe, Y. Nakamura, and S. Uchida, Nature 375, 561 (1995).
-  J. M. Tranquada, N. Ichikawa, and S. Uchida, Phys. Rev. B 59, 14712 (1999).
-  M. Matsuda, M. Fujita, K. Yamda, R.-J. Birgeneau, Y. Endoh, and G. Shirane, Phys. Rev. B 65, 134515 (2002).
-  N. B. Christensen, H. M. Rønnow, J. Mesot, R. A. Ewings, N. Momono, M. Qda, M. Ido, M. Enderle, D. F. McMorrow, and A. T. Boothroyd, Phys. Rev. Lett. 98, 197003 (2007).
-  S. R. Dunsiger, Y. Zhao, B. D. Gaulin, Y. Qiu, P. Bourges, Y. Sidis, J. R. D. Copley, A. Kallin, E. M. Mazurek, and H. A. Dabkowska, Phys. Rev. B 78, 092507 (2008).
-  M. Kofu, S.-H. Lee, M. Fujita, H.-J. Kang, H. Eisaki, and K. Yamada, Phys. Rev. Lett. 102, 047001 (2009).
-  M. Hücker, M. v. Zimmermann, G. D. Gu, Z. J. Xu, J. S. Wen, G. Xu, H. J. Kang, A. Zheludev, and J. M. Tranquada, Phys. Rev. B 83, 104506 (2011).
-  P. Abbamonte, A. Rusydi, S. Smadici, G. D. .Gu, G. A. Sawarzky, and D. L. Feng, Nature Phys. 1, 155 (2005).
-  M. Hücker, M. v. Zimmermann, M. Debessai, J. S. Schilling, J. M. Tranquada, and G. D. Gu, Phys. Rev. Lett. 104, 057004 (2010).
-  S. B. Wilkins, M. P. M. Dean, J. Fink, M. Hücker, J. Geck, V. Soltwisch, E. Schierle, E. Weschke, G. Gu, S. Uchida, N. Ichikawa, J. M. Tranquada, and J. P. Hill, Phys. Rev. B 84, 195101 (2011).
-  M. P. M. Dean, G. Dellea, M. Minola, S. B. Wilkins, R. M. Konik, G. D. Gu, M. Le Tacon, N. B. Brookes, F. Yakhou-Harris, K. Kummer, J. P. Hill, L. Braicovich, and G. Ghiringhelli, Phys. Rev. B 88, 020403(R) (2013).
-  G. Fabbris, M. Hücker, G. D. Gu, J. M. Tranquada, and D. Haskel, Phys. Rev. B 88, 060507 (2013).
-  A. Shengelaya, M. Bruun, B. I. Kochelaev, A. Safina, K. Conder, and K. A. Müller, Phys. Rev. Lett. 93, 017001 (2004).
-  O. N. Bakharev, I. M. Abu-Shiekah, H. B. Brom, A. A. Nugroho, I. P. McCulloch, and J. Zaanen, Phys. Rev. Lett. 93, 037002 (2004).
-  D. Poilblanc and T. M. Rice, Phys. Rev. B 39, 9749(R) (1989).
-  M. Kato, K. Machida, H. Nakanishi, and M. Fujita, J. Phys. Soc. Jpn. 59, 1047 (1990).
-  J. Zaanen and O. Gunnarsson, Phys. Rev. B 40, 7391 (1989).
-  M. Fleck, A. I. Lichtenstein, and A. M. Oleś, Phys. Rev. B 64, 134528 (2001).
-  M. B. Zölfl, T. Maier, T. Pruschke, and J. Keller, Eur. Phys. J. B 13, 47 (2000).
-  C. Weber, K. Haule, and G. Kotliar, Phys. Rev. B 78, 134519 (2008).
-  L. de’ Medici, X. Wang, M. Capone, and A. J. Millis, Phys. Rev. B 80, 054501 (2009).
-  C. Weber, K. Haule, and G. Kotliar, Nature Phys. 6, 574 (2010).
-  A. Fujimori, A. Ino, J. Matsuo, T. Yoshida, K. Tanaka, and T. Mizokawa, J. Electron Spectrosc. Relat. Phenom. 124, 127 (2002).
-  A. Georges, G. Kotliar, and W. Krauth, Z. Phys. B 92, 313 (1993).
-  Yu. G. Pashkevich and A. E. Filippov, Phys. Rev. B 63, 113106 (2001).