Density-functional Theory for f electron Systems: the {\alpha}-{\gamma} Phase Transition in Cerium

The isostructural {\alpha}-{\gamma} phase transition in cerium is analyzed using density-functional theory with different exchange-correlation functionals, in particular the PBE0 hybrid functional and the exact- exchange plus correlation in the random-phase approximation [(EX+cRPA)@PBE0] approach. We show that the Hartree-Fock exchange part of the hybrid functional actuates two distinct solutions at zero temperature that can be associated with the {\alpha} and {\gamma} phases of cerium. However, despite the relatively good structural and magnetic properties, PBE0 predicts the {\gamma} phase to be the stable phase at ambient pressure and zero temperature, in contradiction with low temperature experiments. EX+cRPA reverses the energetic ordering, which emphasizes the importance of correlation for rare- earth systems.

The first-principles description of f electron materials is a considerable challenge and a highly debated topic in condensed-matter physics. The simultaneous presence of itinerant spd states and localized partially occupied f states and their mutual interaction in rare-earth materials give rise to a rich variety of physical phenomena that continue to be a testing ground for electronic-structure theories. Cerium is one of the most prominent representatives in this regard and, even more intriguingly, undergoes an isostructural (fcc) α-γ phase transition accompanied by a volume collapse of 15% at room temperature and ambient pressure [1,2]. The α phase is characterized by enhanced Pauli paramagnetism and has a smaller volume, while the larger-volume γ phase follows a Curie-Weiss behavior for the magnetic susceptibility.
At zero temperature, first-principles calculations have so far been unable to produce a double minimum in the total energy versus volume curve, that would be a direct indication of the phase transition, within a single theoretical framework. In local or semilocal (LDA or GGA) functionals of density-functional theory (DFT) the f electrons are always delocalized, due to the strong self-interaction error of the functionals, and only the α phase is described with some confidence [3,4]. The selfinteraction corrected local spin density approximation [3,5,6] and Hubbard U augmented local or semilocal DFT calculations (LDA/GGA+U) [4,7] enforce localization of the f electrons. They subsequently yield a phase, whose volume and magnetic moment are consistent with the γ phase, but the description of the α phase requires a different treatment, namely LDA. Dynamical mean field theory in combination with LDA (LDA+DMFT) has been applied to the study of the phase transition at finite temperatures [8][9][10][11], but could so far not be extended to the zero temperature limit. Whether a double minimum exists or would only emerge in the free energy curve at finite temperature, due to entropic effects as suggested by Amadon et al. [10], is therefore still a matter of debate.
In this Letter we show that hybrid density functionals [12][13][14], that incorporate a fraction of exact exchange yield a double minimum within a single theoretical and computational framework. The results are further improved quantitatively by employing exact exchange plus correlation in the random-phase approximation (EX+cRPA) (see Ren et al. [15] and references therein) [16]. In our approach all electrons are treated on the same quantum mechanical level, in contrast to LDA/GGA+U or LDA+DMFT studies. We obtain two distinctly different solutions, whose structural, electronic and magnetic properties are consistent with experimental results for the α and γ phase, respectively.
All calculations in this work were performed with the all-electron code FHI-AIMS (Fritz-Haber-Institut ab initio molecular simulation) [17,18], that is based on numeric atom-centered orbitals. Relativistic effects are treated at the level of the scaled zero-order regular approximation [19]. Here we present results obtained using the PBE0 hybrid functional [13] for both cluster and periodic systems [20] and show that the HSE hybrid functional [14] yields a similar description. For comparison we also applied the local-density approximation in the parameterization of Perdew and Zunger [21] and the Perdew-Burke-Ernzerhof generalized gradient approximation (PBE) [22]. Periodic calculations were performed with a 6×6×6 k mesh. This gives energies that are converged to within 5 meV, which is sufficient for the energy scale of interest here [cf. Fig. 1 (c)]. The hybrid functional calculations were carried out with a tier 1 numeric atom-centered orbitals basis [17], whereas for (EX+cRPA)@PBE0 it proved necessary to go up to tier 3 [23]. In both cases ferromagnetic ordering is assumed in our spin-polarized calculations.
In cerium the 4f , 5d, and 6s states all lie in the vicinity of the Fermi energy, giving rise to different electronic configurations that are very close in energy. In cases like this, approaches that are based on the density matrix rather than the density, such as DFT+U , hybrid functionals or Hartree-Fock are more susceptible to local minima in  the potential-energy landscape of the electrons [24][25][26][27][28][29][30]. This fact can be exploited to search for distinct, stable electronic configurations. For a given lattice constant a 0 our PBE0 calculations are initialized with the density matrix from a preceding PBE calculation and with a high electronic temperature (T ∼1000 K). The electronic temperature yields a broadening of the one-electron energy levels. In subsequent calculations the temperature is gradually reduced until a particular solution is stabilized at T = 0. When scanning a range of lattice constants (a 0 ) for cerium metal, the binding energies from different values of a 0 fall into two different smooth curves. In practical calculations, the solution at one particular lattice constant can be used to initialize the calculations at a neighboring lattice constant. In this way, the two binding energy curves can be stabilized with relative ease. It is one of the core results of this work that the above technique provides two different stable solutions. In Fig.  1 the cohesive energies (E coh ) obtained from LDA, PBE, PBE0, and HSE06 are presented as a function of the lattice constant. Our LDA and PBE results are in agreement with previous calculations [3,4] and exhibit only one minimum. The associated volume is consistent with the α phase, although the actual value is underestimated. Constraining the magnetic moment does not introduce a second minimum. In contrast, in PBE0 and HSE06 two stable solutions are found. One solution has a minimum approximately coinciding with the LDA or PBE minimum, while the second assumes its minimum at a much larger lattice constant, consistent with the one of the γ phase. The magnitude of the cohesive energy systematically reduces from LDA to PBE, and to PBE0.
The two PBE0 solutions differ in their electronic structure as e.g. reflected in the density of states (see Ref. [31]) and the magnetic moment m (cf Fig. 3 (c)). m of the low volume phase lies around 0.2 µ 0 , while in the high volume phase m is close to one. A rapid change of the local magnetic moment across the α-γ phase transition was also observed in LDA+DMFT [8]. Also the number of f electrons is approximately one in both phases, as suggested by positron annihilation experiments [32]. Figure 2 shows the difference of the PBE0 electron density between the α-and the γ-like solutions, n α (r) − n γ (r), projected onto a volume slice at the same lattice constant of 4.6Å. The α-like phase has a higher density in the interstitial region, a clear indication of electron delocalization, whereas in the γ-phase solution the electron is more localized around the Ce atom. More importantly, however, the density difference has the shape of an f orbital of xyz, z(x 2 − y 2 ) symmetry, as evidenced by its three-dimensional projection (not shown). This provides a strong indication that the balance between localization and delocalization of the f electrons plays a key role in the emergence of the double minimum in the cohesive energy curve.
Further inspection of Fig. 1 reveals that the energetic ordering of the two PBE0 solutions is not consistent with the experimental phase diagram [33], which shows that the α phase is lower in energy than the γ phase at low temperatures. The opposite is true for the PBE0 results. To overcome this discrepancy we resort to a more accurate treatment of exchange and correlation. As the name indicates, within EX+cRPA the exchange term is treated exactly (i.e. not reduced by a factor as in hybrid functionals), and correlation is treated at the level of the random-phase approximation. The mixing factor in the hybrid functionals that controls the fraction of ex-FIG. 2. Volume slice through the difference between the bulk Ce electron densities of the α and γ phases at the same lattice constant of 4.6Å, at which both phases have the same energy (see Fig. 1). The α phase has a larger contribution in the interstitial region, whereas the γ phase is more localized and exhibits a clear f -orbital shape. act exchange may be regarded as a simplified screening function, which is replaced by a physical and systemdependent screening in EX+cRPA. EX+cRPA is so far only implemented for finite systems in the current version of FHI-AIMS. We therefore adopt the strategy of modeling bulk Ce with systematically increasing cerium clusters. The clusters are cut from the fcc crystal structure, with one atom in the center surrounded by shells of first, second, and third nearest neighbors (i.e. thirteen, nineteen, and forty-three atoms in total). To reduce edge effects, we use the formula for evaluating the effective cohesive energy of clusters (see, e.g., [34,35]) where E is the total energy, N c the number of atoms in the cluster with c nearest neighbors, and E atom c is the atomic total energy for a c-fold-coordinated atom. Basis-set superposition errors are corrected when evaluating E atom c . For each cluster two distinct solutions were always found. Figure 3 demonstrates how the two sets of PBE0 cohesive energies, equilibrium lattice constants, and magnetic moments (on the central atom) of the clusters converge towards the corresponding bulk values as the cluster size grows. It can be seen from Fig. 3 that the essential physics is already captured with the smallest cluster size. In the following, the (EX+cRPA)@PBE0 results for the Ce 19 cluster will be shown, which suffices for the discussion below.
The EX+cRPA calculations are performed non-selfconsistently with the orbitals and eigenenergies of the two PBE0 solutions as input. Figure 4 shows the (EX+cRPA)@PBE0 cohesive energy for the 19-atom cluster. Since the density of states near the Fermi level is higher in the low-volume phase, and consequently screening and the RPA correlation energy are larger, the lowvolume phase moves down in energy relative to the highvolume one and the correct energetic order is restored. According to the extrapolation of the experimental data [10,33], the difference in internal energy (∆U ) between the two phases should lie between 20 and 30 meV. The (EX+cRPA)@PBE0 value for the 19-atom cluster amounts to ∆U ≃ 45 meV. Although the difference is larger than the experimental estimation of the maximum energy difference between the two phases, it is comparable to the entropy contribution T ∆S at ambient conditions [33]. Therefore, our results on the electronic contribution to the phase transition do not rule out that entropy might play a noticeable role in the phase transition at finite temperature as proposed in Ref. 10. The common tangent construction to the (EX+cRPA)@PBE0 cohesive energy curves leads to a transition pressure of P t ≃ −0.74 GPa at zero temperature, in good agreement with the extrapolated experimental P t ≃ −0.8 GPa. The calculated lattice constants for the α-and γ-like phases are 4.45 and 5.03Å, respectively. This corresponds to a volume collapse of ≃ 30% at zero temperature, to be con- trasted with the 15% observed experimentally at ambient conditions [1]. Over the years many experimental and theoretical studies have addressed the origin of the transition but a conclusive answer is still lacking. The two prevalent propositions are: a Mott transition for the f electrons [37,38] and the Kondo volume collapse [39,40]. In the Mott picture the hybridization between f orbitals is believed to change across the transition, leading to one phase with delocalized f electrons (α phase) and the other with localized f states (γ phase). In the Kondo volume collapse model, the f electrons are assumed to be localized in both phases, and the change in spin screening of the localized moments by the conduction spd electrons is responsible for a change in the Kondo temperature across two orders of magnitude, with an associated change in system properties. Although the two models are different in nature, recent theoretical works suggest that the resulting scenarios are quite similar at finite temperature [41] and both are consistent with available experimental data [2,42]. The common belief that the phase transition should be driven by changes in the electronic structure alone has also been questioned, but the contribution of e.g. lattice vibrations has been estimated to be lower than 30% of the total energy change [43][44][45]. Since our calculations are performed at zero temperature they can only describe a pressure induced phase transition at zero temperature but not a temperature induced transition at ambient pressure as observed in the experimental phase diagram. However, the occurrence of a double minimum in the cohesive energy curve and the strong signature of f states in the accompanying density difference shown in Fig. 2 suggest that the localization of f electrons is a strong contributing factor to the phase transition. Density-functional theory (being a ground state theory) does not need to give the right electronic excitation spectrum to describe the structural phase transition at zero temperature, even if the spectrum is essential for the distinction between the Mott transition and the Kondo volume collapse at finite temperature. Nevertheless, our zero temperature results are consistent with the scenario of a Mott transition and show that advanced DFT exchange-correlation functionals can indeed capture such Mott physics.
In conclusion, PBE0 hybrid functional calculations combined with EX+cRPA produce multiple distinct solutions of the electronic structure of bulk cerium at zero temperature. These can be discriminated by their magnetic moment and the degree of f electron localization and have been tentatively associated with the α and γ phases of cerium. At the PBE0 level, the energetic ordering of the two solutions is reversed compared to the zero temperature extrapolation of the experimental phase diagram. EX+cRPA recovers the right ordering, which highlights the role of correlation in rare-earth systems.
An interesting aspect that emerged from the cluster extrapolation approach is the presence of a volume collapse at the nanoscale down to the dimer. This is a very unusual feature for a first-order phase transition in a metal and opens the question of whether other lanthanides would exhibit the same behavior. This prediction suggests further investigation both theoretically and experimentally.