Structural properties of hypothetical CeBa₂Cu₃O₇ compound from LSDA+DMFT calculations

(Rare-earth)Ba₂Cu₃O₇ (R123) was the first series of high-T_C superconducting materials for which it was found out that a diamagnetic yttrium ion Y could be replaced by a magnetic rare-earth element without influencing the superconductivity [1]. It was also noted that there were three exceptions, namely, Ce, Pr and Tb. Interestingly, ab initio electronic structure calculations showed that electronic bands of PrBa₂Cu₃O₇ [2–5] or the hypothetical stoichiometric CeBa₂Cu₃O₇ [6] and TbBa₂Cu₃O₇ [7] systems have metallic character like in the parent superconducting YBa₂Cu₃O₇. Hence, one could suspect that these compounds, if synthesized as a single “123” phase, should exhibit superconducting properties, too. The discovery of superconductivity in Pr123 [8–12] supports this reasoning. From the technological point of view, however, the situation is much more complicated since there is a strong experimental evidence that some fraction of Ba sites in the R123 lattice is occupied by the rare-earth ions (R_Ba defects) [13]. In this picture, Ce+³ and Pr+³ are unique because of their large size which allows them to occupy Ba sites. It was demonstrated that the partial substitution of Pr or Ce has a detrimental influence on sample conductivity and superconductivity, whereas Tb was found to substitute for Y without significant effects on either the normal state conductivity or T_C [14]. Recently, it has been reported that (Y_1-xPr_x)123 and (Y_1-xCe_x)123 polycrys-talline bulk samples exhibit two close and sharp genuine superconducting transitions as a consequence of the phase separation [15, 16]. Ce takes a special place among the other rare-earths substituting for Y in the R123 family of compounds for one more reason. Cerium metal with a single electron in the f shell is a classic example of an electronic transition induced by temperature and/or pressure. At temperatures less than 600 K and pressures less than 2 GPa, it undergoes a transition between two isostructural phases, a high-pressure α phase and a low-pressure γ phase. Noteworthy, the high-temperature γ phase has approximately 15 % larger volume and displays a Curie-Weisslike temperature dependence of the magnetic susceptibility suggesting the existence of local magnetic moments while the α phase shows a Pauli-like temperature independent paramagnetism [17, 18]. It was noted that in the γ phase 4f electrons are localized whereas upon pressure they are more and more delocalized [19]. According to our best knowledge, it is not known whether similar situation could also take place in case of cerium in the Ce123 system.

In the present paper we extend the latter approach [6] to the structural properties of the hypothetical stoichiometric CeBa₂Cu₃O₇ system and present the results of our ab initio investigations with the particular emphasis put on the structural properties at high pressures and temperatures. We have used a LSDA+DMFT method [20–26] which merges the local spin density approximation (LSDA) with the dynamical mean-field theory (DMFT) [27–29] to account for the electronic correlations. This approach has been developed from the band structure and many-body communities have joined two of the most successful approaches of their respective fields. So far, LSDA+DMFT has been successfully employed to properly describe not only pure cerium [22, 24, 30] but also some other Ce-based compounds, such as, e.g., Ce₂O₃ oxide [23, 26], CeP, CeAs, CeSb, CeBi [31] or CeIrIn₅, CeCoIn₅ and CeRhIn₅ [32]. It is commonly known that strongly correlated electron systems could exhibit many unusual properties. They are extremely sensitive to small changes in their control parameters resulting e.g. in tendency to phase separation. This makes the study of Ce123 system still challenging.

In this work zero-temperature ground-state spin polarized electronic structure calculations were carried out using LMTART program [33] in the framework of the density functional theory (DFT) [34, 35]. Our calculations were performed within the atomic sphere approximation (ASA), version of the linear muffin-tin orbitals (LMTO) method [36] coupled to the technique of so-called empty spheres, which enables the treatment of the open structures on an equal footing with the packed structures. ASA uses overlapping atomic spheres (muffin-tin spheres, MT), where the potential is expanded in spherical harmonics inside the spheres but any contribution from the interstitial region is neglected. In this approximation, however, the computational results depend on the values of the muffin-tin radii. Nevertheless, it is regarded as a quite good approach, especially in case of d and f shells, which are typically strongly localized inside the atomic spheres. The MT radii are automatically found by the LMTART program at the beginning of the self-consistent calculations so that the potential is approximately spherical inside the MT sphere. Such method is fast and provides reasonably good energy bands [33]. In case of Ce123, empty spheres were located in the anti-chain O(5) positions. The exchange correlation energy of electrons was described in the local spin density approximation (LSDA) with the parameterization of Perdew and Wang [37]. The spin-orbit coupling matrix elements were also taken into account using a variational procedure proposed by Andersen [38]. The charge density and the effective potential expanded in spherical harmonics up to l = 6 (4, 4, 1) divisions of the Brillouin zone along the three directions for the tetrahedron integration were used in all the calculations. Self-consistency was achieved by demanding the convergence of the total energy to be smaller than 10^-5 Ry.

The relative positions of the constituent atoms in the unit cell, the muffin-tin sphere radii (r_MT), as well as the valence and semicore orbitals included in the initial basis set are listed in Table 1. As a starting point, we assumed the lattice parameters of Ce123 to be comparable with the ones obtained for the orthorhombic Pr123 from the systematic neutron diffraction study of R123 superconductors [39]. Thus, we used b/a = 1.00957 and c/a = 3.03535 as an input and those ratios were kept unchanged during our calculations of the system for different cell volumes V. There was no relaxation of the positions of atoms in the unit cell, as well.

In our ground-state calculations (T = 0 K) we used LSDA+U approach [40–42] in which the underestimation of the interband Coulomb 4f interactions is corrected by the Hubbard U parameter and the self-consistent procedure corresponds to solving the multiband Hubbard model in the Hartree-Fock (HF) approximation. The choice of U is, however, not unambiguous and it is not trivial to determine its value a priori. Hence, U is often fitted to reproduce a certain set of experimental data, for example band gaps and structural properties. Unfortunately, in case of Ce123 there is no appropriate data available. Therefore we decided to apply the value which is comparable with the one reported in the previous investigations of Pr123 system (U p 0.4 Ry) [43, 44]. Four input parameters are required when the on-site Coulomb and exchange interactions between 4f electrons are considered in our calculations. They are Slater integrals F₀, F₂, F₄ and F₆. These integrals are directly related to on-site Coulomb interaction (U) and exchange interaction (J) by relations U = F₀, J = (286F₂+ 195F₄+ 250F₆)/6435 [45]. For rare-earth elements F₄/F₂ = 0.624 and F6/F2 = 0.448 [46]. Finally, the following Slater integrals were used for the Ce123 system: F₀ = 0.44118, F₂ = 0.78812, F₄ = 0.49004, F₆ = 0.35131. Additionally, the Hartree-Fock approximation to the self-energy in the self-consistent calculations was used. The HF self-energy is obtained for f state occupancy, which can be different from the atomic one, therefore, the actual self-consistent f-band occupancy was used in the expression for the double-counting correction. In our finite-temperature investigations (T = 500 K and T = 700 K) the LSDA+DMFT calculations were performed using the so-called Hubbard I (HI) impurity solver [22, 29] and 500 Matsubara frequencies were used, what is approximately twice as much as in the analogous calculations of the compressed lanthanides Ce, Pr, and Nd [30].

In this study the equilibrium cell volume V₀ and the bulk modulus B0 were determined using the most common Birch-Murnaghan (B-M) isothermal equation of state (EOS) [47]. Next, the corresponding lattice parameter a0 was evaluated. The additional advantage of the B-M EOS is that it can be used for very high pressures where pressure derivative B’ varies significantly with pressure.

Fig. 1, Fig. 2, Fig. 3 and Fig. 4 show the total energy versus volume Etot(V) curves obtained during the LSDA, LSDA+U and LSDA+DMFT at T = 500 K and T = 700 K calculations, respectively. All single-point energies were computed by keeping the c/a ratio constant for different cell volumes V. As the starting point, the cell volume of V_ini = 1188.227 a.u.³ was always taken, which is very close to the experimental value obtained for Pr123 system. Noteworthy, we have found out that during the LSDA calculations the default basis set, consisting of Ce 6s, 5p, 5d, 4f valence states and 5s semicore states, leads to the artificial energy minimum which appears at V₀ ≈ =3 1143 a.u.³ < V_ini (not shown). According to experimental data [39], the lattice parameters in the R123 family of compounds increase with the rare-earth ionic radius and thus one could rather expect that the cell volume of Ce123 should be greater than Vini. When Ce 5p semicore states as well as Ce 6p valence states are included (Table 1), the minimum is shifted towards higher values giving much more realistic equilibrium cell volume. Similar modifications of the LMTO basis set were done in the previously reported calculations of Ce metal and Ce₂O₃ oxide [23]. It was noted by the authors of the cited work that without the Ce 5p semicore there was no minimum of the total energy versus volume curve in Ce₂O₃. Also the underestimation of the lattice parameter of pure Ce metal occurred in their calculations when Ce 6p states were neglected.

Fig. 1

As the energy versus volume curve is considered in the LSDA+U and the LSDA+DMFT, we can split the energy into two contributions. The first one includes the interaction U and the second includes all the other terms. If we separate these two contributions as a function of volume, the first term in the DFT+DMFT calculations is practically negligible. It is because the number of f electrons found in the Hubbard I solver is about one in this approximation for each volume (more precisely, it slightly decreases with the increased volume from 1.08 e to 1.05 e). On the other hand, in the LSDA+U the energy decomposition shows a slightly different behavior. It comes directly from the fact that the number of correlated electrons in this approach decreases from 1.39 e to 1.31 e as the volume increases. Generally speaking, the analysis of the variation of the energy as a function of volume leads to the conclusion that it is mainly correlated to the variation of the number of f electrons. Moreover, the orbital anisotropies are different in the LSDA+U and the LSDA+DMFT calculations. In contrast to the LSDA+U, the f electrons in the LSDA+DMFT approach are more homogeneously distributed because of the possible fluctuations between different localized atomic states, hence the diagonal terms in the occupation matrix of the local orbital basis are similar. Recently, these effects have been analyzed in detail for the analogous calculations of Ce and Ce₂O₃ systems [26].

Fig. 2

Fig. 3

To analyze structural properties of the Ce123 system the equilibrium cell volume V₀ and bulk modulus B₀ were determined within the methodology described before. All values of V₀ and B₀ obtained from the Birch-Murnaghan fits as well as the evaluated lattice parameters a₀, are collected in Table 2. As can be easily seen from Fig. 1 and Fig. 2, the distinct global minimum, corresponding to the zero-temperature static stability of the Ce123 system, exists in both the LSDA or the LSDA+U approaches, and the Birch-Murnaghan equation of state gives an excellent description of the system in the full range of the examined cell volumes. Unfortunately, there is no available structural data for the Ce123 compound to be compared directly with our results. However, we can use the interpolation of the existing experimental data for R123 homologues to mimic the general tendency and to enable some comparison with the theoretical values obtained in this work. It is known from the literature that the lattice parameters of the nearest Ce123's “neighbors” in the R123 series (Pr123 and La123), evaluated from the low-temperature (10 K) neutron-diffraction measurements [39], have the following values: a(Pr123) = 7.2925 a.u. and a(La123) = 7.3756 a.u. One can notice that the LSDA overestimates the lattice parameter (a₀ w a(La123)), whereas the LSDA+U seems to give the lattice parameter with the reasonable intermediate value.

Fig. 4

Fig. 5

On the other hand, it is clear from Fig. 3 and Fig. 4 that in the LSDA+DMFT approach at the finite temperatures of 500 K and 700 K the situation is more complicated and some discrepancy between the B-M EOS and the calculated E-V curves is visible. It is why the B-M fits cannot be used for all the computational data to receive the acceptable accuracy. Instead, we performed the fits only for the single-energy points in the nearest vicinity of the minimum. At 500 K one can see the global minimum indicating the existence of the stable phase with the equilibrium cell volume higher than in the LSDA+U approach. Noteworthy, in this case the corresponding lattice parameter is higher than a(La123). What is more, the evident “softness” in the total energy is visible for V < 1128.8 a.u.³ (i.e. at the high pressure), which could be interpreted as a trace of the metastable phase. In Fig. 5 pressure versus cell volume P(V) curve obtained from the Birch-Murnaghan equation of state at T = 500 K is shown. As could be seen, the hypothetical low-volume metastable state appears at the high pressure of P > 30 GPa. One can suspect that this “volume collapse” of Ce123 should have similar nature as the effect that was observed in the pure cerium (i.e. y-a transition). At first glance, the E-V curve at 700 K, presented in Fig. 4, is much more distorted than the one at 500 K. It is also worth to note that at 700 K the B-M fit is relatively poor. Although the distinct global minimum exists in the E-V curve too, it is difficult to determine its position and curvature with the same precision as it was done at T = 500 K. Nevertheless, one can notice that the lowest single-point energy corresponds to V_min(700 K) = 1188.227 a.u.³, which is 5 % smaller when compared with the analogous value at T = 500 K (V_min(500 K) = 1247.639 a.u.³), as could be seen in Fig. 3. The evaluated lattice parameter, corresponding to the global minimum at T = 700 K, is practically identical as the one obtained from the zero-temperature LSDA+U computational data (Table 2). This result seems to be important because it could indicate that this material displays the unusual property of contracting in volume on heating (or strong expanding on cooling) at high temperatures, that is, the occurrence of large negative volumetric thermal expansion coefficient is possible (very rough estimation leads to a_V w -290 x 10^-6 K^-1). Moreover, the distinct signature of the structural transition appears for V w 1116 a.u.³ in the E-V curve at 700 K in our LSDA+DMFT calculations (Fig. 4). On the other hand, rapid “softening” of the system is visible for the large volumes of V > 1307.050 a.u.³, as well. From the technological point of view, this could mean that the investigated Ce123 crystal structure undergoes very dramatic transitions during the cooling to the room temperature after the high-temperature synthesis (~1200 K), namely, ~7 % cell volume expansion at approximately 500 K followed by the immediate contraction at lower temperatures. It is worth to mention that the LSDA+DMFT calculations at T = 100 K and T = 300 K (not presented here) enable one to evaluate the lattice parameters to be a₀ = 7.2974 a.u. and a₀ = 7.3020 a.u., respectively, which are very close to the value obtained from the LSDA+U fit (Table 2). This unusual behavior could be considered as a potential source of problems with the “123” phase stability. On the other hand, the hypothetical low-volume metastable Ce123 phase should be less defected, because the smaller rare-earth ions have difficulty occupying Ba site in R123 structure being too small to efficiently bond to the nearest neighbors [13]. Therefore it would be interesting to perform the high-pressure synthesis of this intriguing compound in the future.

Our LSDA calculations give the bulk modulus B₀ (Ce123) = 140.7 GPa, as it is shown in Table 2, which is slightly smaller than B₀(Pr123) = 151 GPa reported for Pr123 from the similar LDA computations [48]. One can also see that when the Hubbard correction U is applied for Ce 4f states, the obtained bulk modulus is approximately 20 GPa higher than the one in the LSDA approximation. Although the zero-pressure bulk modulus at 500 K, estimated from the B-M fit in the LSDA-DMFT approach, reaches a little bit higher value it is still comparable with the previous results. Somewhat surprisingly, B₀ estimated from the B-M fit at T = 700 K is enormous, comparable with the bulk modulus of diamond B₀(diamond) = 442 GPa [49]. One can suspect that this overestimation stems mainly from the limitation of the B-M fit to the LSDA+DMFT results obtained at T = 700 K. It is important to notice at this place that the experimental value B₀(Pr123) = 49.9 GPa was reported from ultrasonic velocity measurements of polycrystalline Pr123 [50] and the significant discrepancy between the theory and the experiment is noticeable. It was argued [51] that the results of zero-temperature ground-state electronic structure calculations are not directly comparable with experimental measurements that include zero-point phonon effects and are often taken at room temperature. The temperature and phonon effects can modify the calculated bulk modulus and could be a source of the frequently reported overes-timation of bulk moduli obtained in the framework of the DFT. One should also keep in mind that applied ASA approximation could not describe correctly the bulk modulus in LSDA [23]. In case of the DMFT approach some limitation of the impurity solver should be taken into account, as well, because the Hubbard I approximation is valid for the well localized systems and the low hybridization (i.e. high volumes) [26]. So, the results concerning the energy versus volume curves should be taken only as trends.

And last but not least, in Fig. 6 density of states (DOS) for Ce 4f “spin-up” and “spin-down” states against a background of total DOS of the Ce123 system, obtained from the LSDA+DMFT calculations at T = 700 K for different cell volumes, is presented. As could be seen from the total DOS, valence band broadens by ~2 eV under 20 % compression from V = 1247.639 a.u.³ (V/V₀ = 1.05) to V = 1009.993 a.u.³ (V/V₀ = 0.85). At the largest volume, one can suspect the appearance of pure j = 5/2 lower Hubbard bands (LHB) near —2 eV, and the mixed j = 5/2 and j = 7/2 upper Hubbard band (UHB) above the Fermi level (E_F), with the splitting consistent with U = 0.4 Ry (~5.5 eV ). The UHB is essentially a set of three atomic-like “multiplet” states.

Fig. 6

Similar band structure was reported for Ce₂O₃ from the fully self-consistent LDA+DMFT calculations with Hubbard I approximation [23]. It is visible in Fig. 6 that under the high compression both Hubbard bands lose their spectral weights. At the smallest cell volume the LHB practically disappears at the expense of the new growing peak at E_F and the first UHB peak is shifted towards the Fermi level. It is known from the literature [30] that rapid growth in 4f spectral weight at the Fermi level in Ce-based systems is attributed to the so-called Kondo resonance. This effect appears at the expense of the Hubbard side bands, followed by the associated “softness” in the total energy. The most dramatic changes in these signatures for case of Ce coincide with the two-phase region in the γ-α transition. Thus, the appearance of the 4f peak at E_F in our LSDA+DMFT calculations for the high compression could be interpreted, to some extent, as a sign of the isostructural transition in the examined Ce123 system. What is more, this effect is consistent with the previously discussed E-V results concerning the structural properties of this system at the low-volume region. The spin-polarized DMFT calculations can also indicate the following “magnetic scenario” for the possible phase transition in Ce123. For large unit cells the f-electron spectrum is split by the Hubbard repulsion (V/V₀ = 1.05 in Fig. 6), resulting in the appearance of local spin moments at cerium sites (μ = 0.0023 μ_B). With decreasing volume (V/V₀ = 0.85) a quasi-particle is formed at the Fermi level, what leads to the decrease in the spin moment (μ = 0.0019 μ_B). Noteworthy, this effect cannot be reproduced in the LSDA and LSDA+U calculations. In the both cases spin moments of Ce increase during the volume compression from approximately 0.0013 μ_B to 0.0015 μ_B in the LSDA and from 0.0016 μ_B to 0.0052 μ_B in the LSDA+U approach.

In this work we have investigated the structural properties of the hypothetical stoichiometric CeBa₂Cu₃O₇ (Ce123) compound, which has not been synthesized as a single phase yet. We have presented the results of our ab initio calculations carried out by the use of the density functional theory (DFT). We utilized a method which merges the local spin density approximation (LSDA) with the dynamical mean-field theory (DMFT) to account for the electronic correlations using the Hubbard I solver. The LSDA+DMFT calculations were performed for T = 500 K and T = 700 K. The particular emphasis was put on the pressure-induced changes in the electronic band structure related to strongly correlated 4f states. As a reference, we applied the zero-temperature ground-state conventional LSDA scheme and the LSDA+U approach with the Hubbard U correction, where the self-consistent procedure corresponds to solving multiband Hubbard model in the Hartree-Fock approximation. For the simplicity, we have not relaxed the positions of atoms and the cell volume was the only variable of the system. Next, we performed the analysis of the structural properties using the Birch-Murnaghan (B-M) isothermal equation of state (EOS).

The most important results could be listed as follows:

The LSDA+DMFT calculations indicate the occurrence of a large negative volumetric thermal expansion coefficient near T = 500 K. As a consequence, the Ce123 material could display the unusual property of strong expanding on cooling, what may be responsible for the structural instability during the synthesis.

To conclude, we can expect that some improvement of the accuracy in the structural parameters determination should be achieved by the use of the more sophisticated and computationally costly Quantum Monte Carlo (QMC) solver [22, 29] in the LSDA+DMFT approach. It is also worth to point out at this place that some phase transition under pressure into a structure with different symmetry cannot be excluded and such a scenario should be analyzed in the future, as well. So far, our results concerning the energy versus volume should be treated only as trends.