Pressure dependency of localization degree in heavy fermion CeIn3: A density functional theory analysis

Two dramatic discrepancies between previous reliable experimental and ab initio DFT results are identified to occur at two different pressures in CeIn3, as discussed through the paper. We physically discuss sources of the phenomena and indicate how to select an appropriate functional for a given pressure. We show that these discrepancies are due to the inaccuracy of the DFT + U scheme with arbitrary Ueff and that hybrid functionals can provide better agreement with experimental data at zero pressure. The hybrid B3PW91 approach provides much better agreement with experimental data than the GGA + U. The DFT + U scheme proves to be rather unreliable since it yields completely unpredictable oscillations for the bulk modulus with increasing values of Ueff. Our B3PW91 results show that the best lattice parameter (bulk modulus) is obtained using a larger value of α parameter, 0.4 (0.3 or 0.2), than that of usually considered for the AFM phase. We find that for hybrid functionals, the amount of non-local exchange must first be calibrated before conclusions are drawn. Therefore, we first systematically optimize the α parameter and using it investigate the magnetic and electronic properties of the system. We present a theoretical interpretation of the experimental results and reproduce them satisfactorily.


Computational Details
All the calculations presented here are performed in the framework of the density functional theory 20,21 by employing the full potential linearized augmented plane wave (LAPW) method 22,23 , as embodied in the WIEN2k code 24 , in the presence of the spin-orbit coupling (SOC). Since we intend to investigate the compound under question in its AFM phase, we construct a 2 × 2 × 2 supercell which contains two nonequivalent Ce atoms, one with spin up (Ce ↑ ) and the other with spin down (Ce ↓ ). We consider a mesh of 182 special k points in the irreducible wedge of the first Brillouin zone that corresponds to the grids of 12 × 12 × 12 in the Monkhorst-Pack scheme 25 for the supercell. We choose the muffin-tin radii (R MT 's) to be 2.2 a.u. for In and 2.8 a.u. for Ce. The cutoff parameters K max = 7/R MT and l max = 10 are used for the expansions of the wave functions inside the muffin-tin spheres and the interstitial region in terms of the lattice harmonics and the plane waves, respectively. The periodic charge density and potential are Fourier expanded up to G max = 16 Ry. In addition to the LDA and GGA approximations, we use the LDA + U approach with several U eff and B3PW91 hybrid with several α as well as B3LYP functionals. In the hybrid functionals generally some portion of the exchange term in the semilocal functional is replaced by the Hartree-Fock (HF) exchange so that general form of these functionals can be written as: where ∆E x B88 is the Becke's 1988 gradient correction to the exchange functional 26 , E c LYP is the correlation functional of Lee, Yang and Parr 27 , E c VWN is the VWN local-density approximation to the correlation functional 28 , and ∆E c PW91 is the Perdew-Wang gradient correction for the correlation functional 29 . Becke optimized the coefficients to α = 0.2, β = 0.72 and γ = 0.81 18 . In this work we allow to vary the α parameter in the range of 0 to 0.4 to seek better optimization for our case. One would notice that the LDA + U and hybrid can be considered as semiempirical schemes since by varying the adjustable U parameter in LDA + U and α parameter in hybrid schemes one can try to tune the methods to predict more closely some experimental data while deviating more for other data. Thus, in this case it is important to optimize these parameters by taking different physical quantities into account and making sure that the optimized parameters can optimally reproduce the experimental data corresponding to the different physical quantities. In structural properties section, we will follow this strategy and optimize the α parameter in B3PW91 to obtain both lattice parameter and bulk modulus in closer agreement with the experimental data. In magnetic properties section, we will show that our optimized α parameter can predict also the magnetic moment in better agreement with the experiment compared to the other values of the α parameter.
Structural properties. Here, we intend to evaluate the accuracy of the B3PW91 hybrid functional 18,19 as our selected successful method in the rest of this work in predicting the lattice constant and bulk modulus of the CeIn 3 compound in the FM and AFM phases. To this end, the results, as tabulated in Table 1, are compared with the experimental data and theoretical results reported by the others. However, one should keep in mind that the correctness of a theoretical description can be only partially validated by comparison with the experimental data or, eventually, by comparison to theoretical results obtained with significantly more sophisticated methods. The structural properties are obtained by calculating the total energies of the AFM primitive unit cell 30  computational details section) as a function of its volume and fitting the data with the Birch-Murnaghan isothermal equation of state 31 within a variety of exchange-correlation functionals. The results together with the previous experimental 32, 33 and theoretical 11 data are presented for comparison in Table 1. The calculated lattice constants show that their predicted values are not very sensitive to the magnetic ordering as imposed in FM and AFM phases, see Table 1. Hence, we may only discuss the AFM results. The results show that the predicted lattice parameter with the PBE-GGA is larger than those of LDA and WC-GGA. Although our calculated lattice parameter by PBE-GGA is surprisingly identical to the experimental value, PBE-GGA is not successful enough in prediction of bulk modulus compared to that of WC-GGA. In order to simultaneously reproduce the experimental bulk modulus as well as lattice parameter, we use the so called LDA + U and hybrid approaches as the other alternatives. In these methods the degree of localization can be tuned by their adjustable parameters, i.e., U eff in LDA + U and α in hybrid schemes. M. Caffarel and coworkers recently showed that the spatial distribution of the DFT spin density distribution was very sensitive to and depended on the fraction of HF exchange used in the functional employed for CuCl 2 molecule having a single unpaired electron 34 . They showed that for providing a meaningful description for this physical quantity in this molecule using the B3LYP-DFT, the percentage of HF exchange used had to be increased up to about 40%. Therefore, we perform a systematic exploration of the role of the arbitrary value of U eff in GGA + U and dependence of α in B3PW91 functional schemes for the AFM phase of the system. In order to provide a better understanding of the results, lattice parameter and bulk modulus are also presented in Fig. 1(a,b) showing how these basic properties change when U eff in GGA + U and α in B3PW91 vary. The lattice constant is more increased by increasing the U eff of the GGA + U in worse agreement with the experimental data than PBE-GGA, see Table 1 and Fig. 1(a). The results also reveal that the bulk moduli predicted by GGA + U with all of the values considered for its U eff parameter are farther away from the experimental result compared to that of PBE-GGA, see Table 1 and Fig. 1(b). This indicates that the effects of GGA + U on the lattice constant and bulk modulus are in opposite directions of the results improvement over the PBE-GGA results. If we change the base of the exchange-correlation energy from PBE-GGA to LDA and use LDA + U instead of GGA + U, however, we can anticipate to overcome the problem. Our results verify this anticipation, because as can be seen from Table 1 35 the larger (smaller) lattice parameter of the fcc γ-Ce (α-Ce) in agreement with experiment. The B3LYP with α = 0.2 lattice parameter is larger and bulk modulus is smaller than those of GGA + U. Therefore, the onsite hybrid B3LYP functional only aggravates the unsatisfactory GGA + U results and thereby yields the worst presented results in Table 1. In the hybrid B3LYP functional, LYP and VWN5 correlation functionals are used, while in the hybrid B3PW91 functional PW91-LDA and PW92-LDA correlation functionals are used. As discussed above, the overall effect of the included correlations in LDA + U is more successful than that of GGA + U for the system under study. From this we can anticipate that the overall effect of the included correlations in B3PW91 due to its included PW91-LDA and PW92-LDA terms can be also more efficient than that of B3LYP functional for the compound under study. This anticipation is confirmed by the results presented in Table 1 and Fig. 1. According to results presented in Table 1 and Fig. 1(a), the lattice parameter increases as α parameter increases with the B3PW91 XCF. However, this is not the case for the bulk modulus, see Fig. 1(b). The bulk modulus of the system, as shown in Fig. 1(b), fluctuates as α parameter increases from zero to 0.4 by step 0.05. The B3PW91 XCF with α = 0.4 predicts the lattice parameter in the best agreement with the experimental result compared to other values of α, however, it fails to satisfactorily predict the bulk modulus. Figure 1 shows that the best agreement between theory and experiment may not be necessarily achieved only by increasing the α parameter and both of the lattice parameter and bulk modulus would be taken into account for optimizing the α parameter. We also calculate the structural properties using the LDA-EECE method. However, this method, as shown in Table 1, is not very successful in prediction of CeIn 3 lattice constant and bulk modulus. The phase transition pressure and coexistence pressure range were predicted 14,16 for ε and ζ solid phases of O 2 by B3PW91 with 20% of HF exchange to be 75 GPa and 75-145 GPa in an acceptable agreement with the experimental data 36 of 95 GPa and 95-110 GPa, respectively. Structural parameters of the ε and ζ phases were also predicted 14 by B3PW91 in excellent accordance with experiments. These reports in agreement with our results show that the performance of hybrid functional methods strongly depends on the amount of non-local HF exchange. Now, by taking the above discussed points into consideration, the time seems apt to conclude this section by attaining without proof that the overall effects of the hybrid B3PW91 functional with α = 0.2 on the structural properties of the CeIn 3 compound are more impressive and better than the other functionals.

Magnetic properties.
It is generally believed that magnetic ordering can be changed by pressure 10,[37][38][39] . This is the case for the compound under question too, since experimental results 10 show that the magnetic ordering Néel temperature, T N , of the CeIn 3 compound decreases by applying pressure. The previous theoretical calculations 11 using PBE-GGA functional also confirm this fact that the magnetic moments of CeIn 3 decrease by imposing pressure. Our PBE-GGA calculations also validate the latter theoretical results, as shown in Table 2 and Fig. 2.
So far everything has been qualitatively fine as expected. But, two sizable discrepancies can be observed when the theoretical predictions are quantitatively compared with the experimental data. The first discrepancy is concerned with the T N of the compound. The experimental results show that the T N of CeIn 3 compound decreases with increasing pressure and may be completely suppressed at the critical pressure P c = 2.5 GPa, viz.  Table 2 and Fig. 2). A comparison elucidates that the theoretically predicted P s = 16 GPa by PBE-GGA is much larger than the corresponding experimental P c = 2.5 GPa. It was tried to improve this discrepancy by using the WC-GGA functional 11 . Although the theoretical WC-GGA P s = 9 GPa is closer to the corresponding experimental P c = 2.5 GPa, the predicted P s using the WC-GGA functionals is still very far from P c . The second discrepancy is concerned with the TOT of the system. Reliable neutron-diffraction results show that the TOT of the CeIn 3 is 0.65 ± 0.1 μ B at T = 5 K and zero pressure 13 . However, previous theoretical study 11 at zero pressure showed that the TOT was 0.170 μ B within PBE-GGA and 0.112 μ B within the WC-GGA. Our PBE-GGA calculations lead to the TOT = 0.24 μ B at zero pressure which is fairly comparable with the previous PBE-GGA result. In addition to spin magnetic moment per atom in the muffin-tin shpere, we have also considered and included the spin magnetic moment per atom in the interstitial region and this is why our calculations a little bit differ from the previous theoretical work. Regardless of the latter small difference between theoretical results, a comparison elucidates that the theoretically predicted TOT's by PBE-GGA and WC-GGA are much smaller than the corresponding experimental TOT = 0.65 μ B . Therefore, these functionals are not enough successful in predicting the TOT at zero pressure.
As discussed above there are two dramatic inconsistencies between theory and experiment: 1) theoretical P s is far from the experimental P c , and 2) there are a large discrepancy between the theoretical and experimental TOT's at zero pressure. We show that the discrepancies between experiment and theory originate from high pressure dependency of the Ce 4f-electrons localization. On the other hand, the validity of the selected exchange-correlation functional for simulating CeIn 3 are highly sensitive to the pressure regime at which the DFT calculations are performed. Therefore, these discrepancies cannot be well resolved simultaneously only by increasing the accuracy of the calculations within a single specific currently available functional as they occur in two different pressure regimes.
The main idea to overcome these inconsistencies is inspired from the following experimental observations. The dHvA experiments revealed that the character of 4f-electrons was changed from localized to itinerant above the P c 40,41 . These results elucidate that the degree of localization can play important roles in the properties of this compound near the P c . This has motivated us to investigate the effects of the localization degree of Ce 4f-electrons and pressure on the magnetic properties of CeIn 3 . For this purpose, we calculate magnetic moments per Ce atom of this compound versus pressure within several exchange-correlation functionals (XCFs). These XCFs contain PBE-GGA + U with varying U eff from 0.5 to 5.5 eV by step 0.5 eV, B3PW91 hybrid functional with varying α from 0 to 0.4 by step 0.05, B3LYP hybrid functional with α = 0.2, LDA, LDA-EECE, and PBE-GGA. The degree    of localization is adjusted by U eff and α parameters of the GGA + U and hybrid methods, respectively. Our results are presented in Table 2 together with experimental 13 and other theoretical 11,12 results for comparison. The GGA + U results presented in Table 2 show that TOT decreases as pressure increases. The same result can be seen for B3PW91 and B3LYP. Even though it is possible (at least in principle) to find the P s using the GGA + U, B3LYP and B3PW91, we anticipate that the P s using these functionals must be much larger than experimental P c . This anticipation can be verified by the results given in Table 2, since TOT is not zero at pressures very higher than P c for all of the considered U eff values using the GGA + U and all of the α parameters with the B3PW91. Consequently, like GGA + U, the hybrid functionals cannot also predict the P s in acceptable agreement with experiment. These results indicate that increasing the 4f-correlation using GGA + U and hybrid functionals even make worse the predicted P s than that of the PBE-GGA with respect to the P c . Therefore, the 4f-correlation must be decreased to improve the P s compared to the experimental P c . As another alternative, we then use the LDA functional which has the least correlation. The LDA calculations predict that P s = 0, because TOT is zero above the zero pressure, see Table 2. The correlation of LDA is smaller than that of GGA and this enables LDA to substantially improve the P s prediction in much better agreement with the experimental value P c compared to those of PBE-GGA and even WC-GGA as well as the other considered functionals. Hence, LDA could solve the first inconsistency problem between theoretical P s and experimental P c by decreasing the degree of 4f localization.
We now concentrate on the second problem, i.e., the disagreement between theoretical and experimental TOT at zero pressure. The LDA functional leads to zero TOT at zero pressure. Thus, LDA deteriorates the TOT predictions of PBE-GGA and WC-GGA at zero pressure compared to the experiment 13 , while conversely GGA + U and hybrid approaches improve them, see Table 2. The results, as tabulated in this table, show that the B3PW91 with α = 0.2 is the most successful functional for the TOT prediction at zero pressure among the other considered functionals, because in this case TOT becomes 0.67 μ B , which is in more agreement with the experimental result 12 of 0.65± 0.1 μ B . As can be seen from the results presented in Table 2, the TOT at zero pressure is not improved by changing the α parameter to values more or less than 0.2. Therefore, α = 0.2 would be considered as an optimized parameter for B3PW91 in this case in agreement with the discussion presented in structural properties section. Thus, B3PW91 with α = 0.2 can cause to overcome the second inconsistency by increasing the degree of localization.
As discussed above, the LDA functional can only predict the P s in good agreement with P c , but it is unsuccessful in reproducing the TOT at zero pressure. Conversely, the B3PW91 with α = 0.2 can only predict the TOT at zero pressure, but it fails to predict the P s compared to the experimental P c . All these indicate that increasing the degree of localization increases the TOT and makes the theoretical prediction consistent with the experimental result at zero pressure. This is in agreement with the recent theoretical work 42 . Conversely, reducing the degree of 4f-localization decreases the predicted P s and thereby decreases the inconsistency between the theory and experiment. Hence, the degree of 4f-localization can influence the effects of pressure on the magnetic properties. Therefore, it is not possible to describe magnetic properties of CeIn 3 for every pressure only by one fixed degree of 4f-localization. Thus, the degree of localization depends on the pressure and should be tuned in terms of the pressure. High pressure is equivalent to the low 4f localization, and vice versa. This is in agreement with the mentioned dHvA observations which show the changes in the degree of 4f-localization above the P c 40,41 . These changes have been also observed in the CeSn x In 3−x using the Gd 3+ electron resonance experiments at the critical x c ≈ 0.65 where T N → 0 43 . Since atomic radius of Sn is bigger than the atomic radius of In, increasing the Sn concentration is equivalent to increasing the pressure, which is known as the chemical pressure. Furthermore, the dependency of XCF to the pressure regimes in agreement with our results was also previously observed in ε and ζ phases of solid O 2 14,16 . Electronic structure. As discussed in magnetic properties section, the pressure and degree of 4f-localization can affect the magnetic properties. The pressure effects in turn depend on the degree of 4f-localization. In order to elucidate the effects of pressure and 4f-localization as well as their dependencies on each other, here we turn our attention to the electronic structure of CeIn 3 . To this end, we investigate the electronic structure of CeIn 3 at two different pressures within three different functionals corresponding to three different degrees of 4f-localization. The selected two pressures are P = 0 corresponding to a = 4.68 Å and P = 14 GPa corresponding to a = 4.44 Å. These pressures are calculated by B3PW91 with α = 0.2. The latter 14 GPa pressure is equivalent to the predicted LDA pressure 3.5 GPa, which is close to the experimental P c = 2.5 GPa. Since we intend to investigate the effects of 4f-localization degree on the electronic properties, three functionals are selected. These functionals are B3PW91 with α = 0.2 having high degree of 4f localization, B3PW91 with α = 0.1 having intermediate degree of 4f localization, and LDA having low degree of 4f localization. The total and 4f partial DOSs (Fig. 3), band structures (Fig. 4) and Fermi surfaces (Fig. 5) of CeIn 3 are calculated using these mentioned three methods and two pressures.
Total and 4f partial DOSs. Total DOSs are shown in semicore [see Fig. 3 (a ij )] and valence [see Fig. 3 regions separately and 4f partial DOSs are also shown in Fig. 3 (c ij ), where i = 1, 2 and j = 1 to 3. The DOSs calculated by B3PW91 with α = 0.2 at zero pressure, as shown in Fig. 3(a 11 ), (b 11 ), and (c 11 ), are in good agreement with those of previously calculated 30 by GGA + U with U eff = 5.5 eV. There are two peaks in the semicore region near − 15 eV and − 14 eV energies for the considered three functionals and two pressures, see Fig. 3 (a ij ) (i = 1, 2 & j = 1 to 3). Pressure causes these peaks to be broadened and decreased in magnitude, as can be seen by comparing Fig. 3(a 11 ) with (a 21 ), and Fig. 3(a 12 ) with (a 22 ), as well as Fig. 3(a 13 ) with (a 23 ). However, pressure has no remarkable effect on the locations of these peaks. The results obtained from the considered three functionals show that the semicore states look precisely alike at P = 0, see Fig. 3 (a 1j ) (j = 1 to 3). This observation at zero pressure remains correct at P = 14 GPa as well, see Fig. 3 (a 2j ) (j = 1 to 3). These show that the reduction of the 4f localization cannot significantly affect the semicore states. All these show that the effects of pressure on the semicore states are more considerable than those of the 4f-localization degree. The semicore states are mostly composed of the d-In states while the valence states are mostly made up of 4f-Ce states. Thus, it is expected that the effects of the 4f localization are more considerable on the valence states than the semicore states. The reduction of the 4f localization at zero pressure causes the total DOSs in the valence region to be narrower and pile up with taller height around the Fermi level, see Fig. 3 (b 1j ) (j = 1 to 3). The same result can be clearly seen at P = 14 GPa by comparing the narrow and piled up DOSs around the Fermi level in Fig. 3(b 23 ) with the broad total DOSs shown in Fig. 3(b 21 ,b 22 ). Although the effects of pressure on the valence states is almost negligible using the B3PW91 with α = 0.2 [compare Fig. 3(b 11 ) with (b 21 )], the valence DOSs can be also changed by pressure. For example, the LDA results show that the total valence LDA DOS is substantially affected by the pressure, compare Fig. 3(b 13 ) with (b 23 ). Thus, the effects of pressure on the valence states increase by reduction of the 4f localization. This implies that the effects of pressure on the valence DOS depend on the degree of 4f localization. This dependency is what we have emphasized on it in the magnetic properties of CeIn 3 , see magnetic properties section. The 4f partial DOSs more clearly reveal the pressure dependence of the valence states on the degree of 4f localization. There are two nonequivalent Ce atoms in the considered unit cell. These Ce atoms are labeled as Ce1 and Ce2. The results show that up and down 4f DOSs of Ce1 are asymmetric with respect to each other at zero pressure using all considered functionals, see Fig. 3 (c 1j ) (j = 1 to 3). This is the case for Ce2 at zero pressure as well, see Fig. 3 (c 1j ) (j = 1 to 3). These imply that the up 4f DOS of each Ce atom cannot be canceled out by its own down 4f DOS and thereby yields nonzero total magnetic moment per Ce atom. But, the up 4f DOS of the Ce1 cancels the down 4f DOS of the Ce2 out and vice versa which yields zero for the total magnetic moment per unit cell as expected for the AFM phases. The same results are also observed at P = 14 GPa using the B3PW91 with α = 0.2 and 0.1, see Fig. 3 (c 2j ) (j = 1 to 3). However, an interesting different behavior is observed at P = 14 GPa using LDA, see Fig. 3(c 23 ). In contrast to the LDA 4f DOSs at zero pressure, the LDA 4f DOSs at P = 14 GPa show that up and down 4f DOSs of Ce1 are symmetric with respect to each other. Thus, the up 4f DOS of each Ce atom can be canceled out by its own down 4f DOS and thereby yields zero total magnetic moment per Ce atom. Therefore, LDA predicts that CeIn 3 would behave as a paramagnetic (PM) system at the considered pressure (14 GPa). But, one should not here conclude that the phase transition form AFM to PM occurs at 14 GPa, since LDA results in Table 2 clearly show that the total magnetic moment per Ce atom and its components are all zero for every pressure P ≥ 0. Thus, the LDA predicts that the AFM to PM phase transition occurs at zero pressure. All these indicate that both the pressure and the degree of 4f localization must be simultaneously taken into consideration for a theoretical investigation of the AFM to PM phase transition in CeIn 3 which evidently exhibits the dependence of pressure effects on the 4f localization and vice versa.
Band structures. The band structures per spin are calculated using the considered three functionals and two pressures. The results are presented in Table 3 and Fig. 4. The 4f states are also characterized in Fig. 4; the larger thickness of the bands is corresponded to the more contribution of the 4f states. It would be also noted that the E min and E max of the bands in Table 3 are measured with respect to the Fermi energy, i.e., the negative (positive) energies are located below (above) the Fermi level. For convenience, we below start the band structure discussion by considering the results obtained from the most considered localized functional, i.e., B3PW91 with α = 0.2 (in this work). Then, we report on the results obtained by the intermediately localized functional, i.e., B3PW91 with α = 0.1, and finally concentrate on the least localized functional, i.e., LDA. Thus, the degree of 4f localization is gradually reduced step by step through the discussion presented in this section.
As the first step for describing the band structures using the B3PW91 with α = 0.2 at P = 0 and P = 14 GPa, we concentrate only on Fig. 4(a,b) as well as first to third rows of Table 3. The results show that three bands cross the Fermi level at zero pressure which are labeled as γ 1 , γ 2 and γ 3 . The minimum energy of these bands are negative (E min < 0) while their maximum energies are positive (E max > 0). Their occupancy numbers are also smaller than unity. Applying the pressure causes the bandwidths, occupancy numbers and the shape of the bands to be changed. For example, pressure increases the bandwidth and occupancy number of γ 1 . Furthermore, pressure causes the E max of γ 1 and γ 2 bands to decrease which means that pressure pushes these bands downward below the Fermi level. The E min and E max of γ 3 band are also increased by pressure, which means that pressure pushes the γ 3 upward above the Fermi level. But, the pressure cannot cause these bands to be located completely below or   above the Fermi level. Thus γ 1 , γ 2 and γ 3 cross the Fermi level again at P = 14 GPa. In addition, pressure increases the 4f states around the Fermi level. The thickness of the bands around the Fermi level is larger at P = 14 GPa than P = 0. However, the 4f states at the Fermi level are negligible at the two considered pressures and they are mostly dispersed away from the Fermi level.
We now go to the second step and describe the results with the B3PW91 functional with α = 0.1 at P = 0 and P = 14 GPa. Therefore, we concentrate only on Fig. 4(c,d) besides fourth to sixth rows of Table 3. The results show that at zero pressure the E max of γ 1 , γ 2 and γ 3 are positive while their E min are negative. Their occupancy numbers are also smaller than unity. These results imply that these three bands cross the Fermi level at zero pressure. This is the same as the situation of the first step. The results reveal an interesting event at P = 14 GPa in this step. The E max of the γ 1 band is negative and its occupancy number is unity at P = 14 GPa. This means that pressure causes the γ 1 to be located below the Fermi level, which is quite different from that of the first step. Pressure also changes the E min and E max of γ 2 and γ 3 , but it cannot prevent these bands to cross the Fermi level. Therefore, two bands γ 2 and γ 3 cross the Fermi level at P = 14 GPa compared to the P = 0 situation where three bands (γ 1 , γ 2 and γ 3 ) do this. A manifold of the 4f bands is displayed in Fig. 4(c) which is dispersed near the Fermi level up to 1.5 eV. The same result can be seen at P = 14 GPa. A comparison of Fig. 4(c,d) reveals that the thickness of the bands at Fermi level at P = 14 GPa is larger than those at zero pressure. This implies that the pressure increases the 4f states at the Fermi level.
Let us now consider the last step to describe the LDA results. Therefore, we concentrate on Fig. 4(e,f) together with the seventh to ninth rows of Table 3. The results show that the E max of γ 1 is negative and its occupancy number is unity at zero pressure. This implies that the γ 1 is located below the Fermi level and does not cross it. However, the E max s of γ 2 and γ 3 are positive and their E min s are negative. The occupancy numbers of γ 2 and γ 3 are also smaller than unity. This implies that γ 2 and γ 3 cross the Fermi level. The same results are also obtained at P = 14 GPa. This means that pressure cannot change the number of the bands which cross the Fermi level using the LDA functional. However, pressure changes the position, bandwidths, occupancy numbers and the shape of the bands. For example, pressure increases the bandwidths of the γ 1 , γ 2 and γ 3 bands. The results also show a manifold of the 4f bands which straddles the Fermi level at zero pressure. This manifold is dispersed up to 1 eV above   Table 3. Maximum energy (E min ) and minimum energy (E max ), bandwidth and occupancy number (Occup) of three γ 1 , γ 2 and γ 3 bands using the B3PW91 with α = 0.2 and α = 0.1 as well as LDA at P = 0 and P = 14 GPa. E min and E max are measured with respect to the Fermi level so that positive energies are located above the Fermi level while negative energies are located below the Fermi level.
the Fermi level. The same result is also seen at P = 14 GPa. This implies that pressure cannot change significantly the position and dispersion of this manifold using the LDA. In essence, we have shown that the effects of pressure on the band structure can strongly depend on the used functional. This implies that the effects of pressure can be influenced by the degree of 4f localization. Every functional has its own specified degree of 4f localization which will be assumed for and applied on the case under the theoretical study. For example, comparison of Fig. 4(a,c) with (e) shows that the reduction of 4f localization causes the shape of the bands to change and the 4f states to pile up near the Fermi level. The reduction of 4f localization also decreases the bandwidths of γ 1 , γ 2 and γ 3 at zero pressure. We conclude the section by the fact that it is crucial to first specify the pressure, and then for each specified pressure use a suitable functional with an appropriate degree of 4f-localization.
Fermi surfaces. Electronic properties of compounds can be more clearly analyzed by shedding light into their Fermi surface (FS) topologies which contain important experimental information. Several dHvA experiments 40,41,[44][45][46] show that the FS of CeIn 3 compound is changed as pressure increases from P < P c to P ≥ P c , which can be explained by changing the character of 4f electrons from localized for P < P c to itinerant for P ≥ P c . A variety of approaches have been proposed to describe the FS topology of cerium based compounds. The FS measurements can be well described by the DFT + LDA/GGA band structure calculations in some compounds for which their f-electron shells are completely either filled or empty so that contributions of their f-electrons to the corresponding FSs are negligible. Even if the f-electron shells are neither completely filled nor completely empty, they can be still well described by the DFT + LDA/GGA, provided that the f-electrons of the compounds under question are itinerant and behave as band-like electrons. It was observed that the Ce based compounds with f-itinerant electrons had larger FS than those with f-localized electrons 47,48 . These observations revealed that the f-electrons and pressure would play key roles in the FS formation. These have motivated us to investigate the effects of the 4f localization and pressure on the FS of CeIn 3 . Then, we have calculated the FS using the B3PW91 with α = 0.2, B3PW91 with α = 0.1 and LDA functionals at zero pressure and P = 14 GPa. The FS of a compound includes some branches which are produced by the partially filled bands. The branches of CeIn 3 Fermi surface are calculated using the three considered functionals for the two pressures and displayed as S 1 , S 2 and S 3 in Fig. 5. The S j branch is produced by the γ j band provided that the γ j band crosses the Fermi level, as expected from the definition of FS. We have also tabulated contributions of the partial DOSs and interstitial (Int) DOS to the DOS(E F ) in Table 4 using these functionals and pressures.
There are three FS branches with the B3PW91 with α = 0.2 at P = 0, as depicted in Fig. 5 (a 1j ) for j = 1 to 3. These three branches also exist at P = 14 GPa using this functional, see Fig. 5 (a 2j ) for j = 1 to 3. The shapes of these branches at P = 14 GPa are approximately the same as those at P = 0, but the sizes of these branches are changed by pressure. For example, a comparison of Fig. 5(a 11 ) and (a 21 ) shows that pressure only reduces the sizes of the small packets of S 1 branch (located in the corners of the first Brillouin zone) and does not drastically change their shapes. Note that the packets shapes remain unchanged, as pressure does not also remarkably change the shapes of the γ 1 , γ 2 and γ 3 bands at the Fermi level and the points where these bands cross the Fermi level using this functional, see Fig. 4(a,b). The different sizes of the branches originate from the character change of the bands by pressure. As shown in Table. 3, pressure increases the occupation numbers of γ 1 and γ 3 bands but decreases the occupation number of γ 2 band using the B3PW91 functional with α = 0.2. Moreover, the results in Table 4 show that pressure increases the DOS f (E F ) while it decreases the DOS p (E F ), DOS d (E F ) and DOS Int (E F ) using this functional. The DOS s (E F ) is also increased by pressure which is negligible with respect to the DOS f (E F ).
There are also three S 1 ,S 2 and S 3 branches at zero pressure using the B3PW91 with α = 0.1, as shown in Fig. 5 (b 1j ) for j = 1 to 3. A comparison of Fig. 4(a,c) reveals that the shapes of the γ 1 and γ 3 bands at the Fermi level are not significantly changed by decreasing the α parameter at P = 0. This causes the shapes of the S 1 and S 3 branches to remain approximately unchanged at P = 0 by reducing α parameter from 0.2 to 0.1, compare Fig. 5(a 11 ) with (b 11 ) and Fig. 5(a 13 ) with (b 13 ). But, the shape of the γ 2 band is changed at the Fermi level by decreasing the α parameter at P = 0. For example, the γ 2 band is touching the Fermi level between W and L points in Fig. 4(c) while it is not in Fig. 4(a). This causes the shape of the S 2 branch using the B3PW91 with α = 0.2 to be different from that of B3PW91 with α = 0.1, compare Fig. 5(a 12 ) with (b 12 ). This change is due to the degree change of the 4f localization. The results in Table 4 show that the DOS f (E F ) is significantly increased by decreasing the α parameter using the B3PW91 functional while change in the magnitude of the other partial DOSs is negligible at the Fermi level. This causes the character of the bands at the Fermi level and thereby size of the branches to be changed by changing the α parameter. These results are in agreement with previous experimental and theoretical measurements 49 of CeIn 3 FS at AFM phase. In contrast to the case of P = 0, there are only two S 2 and S 3 branches for the case of P = 14 GPa using the B3PW91 with α = 0.1, i.e., S 1 branch is not reproduced at P = 14 GPa after reducing α, see Fig. 5 (b 2j ) j = 1 to 3. Pressure causes the γ 1 band to be located below the Fermi level at P = 14 GPa using the B3PW91 with α = 0.1, as discussed in band structures section. This is why the S 1 branch is disappeared at P = 14 GPa when reducing α, see Fig. 5(b 21 ). A comparison of Fig. 5(b 12 ,b 22 ) elucidates that pressure changes the S 2 branch. This change originates from the effects of pressure on the shape of the γ 2 band at the Fermi level. The γ 2 band is tangent to the Fermi level between W and L points in Fig. 4(c) while it is not in Fig. 4(d). Pressure also fills the holes of the S 3 branch at P = 0; the existed empty spaces in Fig. 5(b 13 ) are filled in Fig. 5(b 23 ) by increasing pressure. The results in Table 4 show that the DOS f (E F ) is remarkably increased at P = 14 GPa while the DOS s (E F ), DOS d (E F ) and DOS Int (E F ) are decreased with respect to those obtained at P = 0 using the B3PW91 with α = 0.1.
Using the LDA functional, as illustrated in Fig. 5 (c ij ) for i = 1, 2 and j = 1 to 3, there are only two S 2 and S 3 branches at P = 0 or similarly at P = 14 GPa, i.e., the S 1 branch is absent. The γ 1 band is located below the Fermi level within this functional and as a result the S 1 branch is absent in the FS. The results as presented in Table 4 show that the magnitude of DOS f (E F ) is much larger than those of the other states. In addition, Fig. 4(e,f) show that the 4f states migrate toward the Fermi level using this functional. These figures also show that the 4f states disperse more uniformly at the Fermi level using the LDA compared to the other functionals. These results indicate that contributions of the 4f states in the two γ 1 and γ 2 bands are high and approximately equal. The shapes of these bands using the LDA are very different from those of the other functionals. Thus, the cross-sectional areas of the FS are filled by LDA much more than the other functionals.
We conclude this section by the following points. There are only two branches for the CeIn 3 FS with LDA at both P = 0 and P = 14 GPa (in agreement with the experimental results 44 in the PM phase), while there are three branches in FS obtained by B3PW91 with α = 0.2 at both P = 0 and P = 14 GPa as well as the FS obtained by B3PW91 with α = 0.1 at P = 0. Although the FSs obtained by B3PW91 with α = 0.1 at P = 14 GPa appears to be also in agreement with the experiment by giving only two branches, we notice that 14 GPa is much less than the predicted P s by B3PW91 with α = 0.1; B3PW91 failed (even much worse than GGA) to predict the experimental P c (see magnetic properties section). This means that P = 14 GPa is greater than P s using the LDA but less than P s using the B3PW91. This implies that for P = 14 GPa the system lies in the PM phase using the LDA, but it still lies in the AFM phase using the B3PW91. These results reveal that the FSs using those functionals that treat f-electrons as itinerant electrons such as LDA are in better agreement with experiment for P > P c than those functionals that treat f-electrons as localized electrons such as B3PW91. Conversely, FSs using those functionals that treat f-electrons as localized electrons such as B3PW91 are in better agreement with experiment for P~0 than those functionals that treat f-electrons as itinerant electrons such as LDA.
Electric field gradient. In this section we take the main point of this work into account and show that the electric field gradient (EFG) as an extremely sensitive quantity can be calculated as a function of pressure and show how the EFG can be affected before and after P c . The electric field gradient (EFG) is a tensor of rank 2 which has only two independent components in the principal axes system. These components are the main component of EFG, V zz , and the asymmetry parameter . For the case under study, η is zero. Thus, we only focus on V zz . The EFG contribution of the electrons inside (at the surface and outside) of the muffin-tin sphere is called the valence (lattice) EFG which is denoted by V zz val (V zz lat ). The EFG quantity is extremely sensitive to the anisotropic charge distributions of the core electrons 50 as well as to the aspherical electron density distribution of valance electrons 51 , and as a result to the valance electronic structure 30 . The EFG, thereby, can serve as a powerful gauge for measuring such a degree of localization 30 .
Here, we calculate the EFGs of CeIn 3 compound using the B3PW91 with α = 0.2, B3PW91 with α = 0.1 and LDA functionals. These calculations are performed by taking the equilibrium lattice constant a = 4.68 Å into account which is predicted by B3PW91 with α = 0.2. The calculated V zz , and their components V zz val and V zz lat are presented in Table 5 along with the other theoretical and experimental results for comparison. The results show that high 4f localized functionals such as LDA + U and hybrid functionals are more successful in EFG prediction at zero pressure than low 4f localized functionals such as LDA or GGA. This result confirms the high localization degree of CeIn 3 4f electrons at zero pressure. Among the high localized functionals, B3PW91 with α = 0.2 is in better agreement with experiment than previous theoretical calculations even GGA + U + SO calculations. The success of this functional can be explained by the relation between EFG and total DOS(E F ). The previous calculations have revealed that EFG ∝ DOS(E F ) 30 . The total DOS(E F ) calculated by B3PW91 (α = 0.2) is .
19 38 [States/(Spin Ry)] which is smaller than that of other functionals 11,30 even GGA + U + SO functional 30 . Therefore, EFG using this functional is also smaller and thereby closer to the experiment than those obtained using other functionals. A very small EFGs at Ce site are also reported in Table 5 which is produced by the SOC interactions. Based on the relation EFG ∝ DOS(E F ) we also anticipate that EFG must increase as the 4f localization decreases, as the DOS(E F ) is increased by decreasing the 4f localization (see Fig. 3). Our calculations confirm this anticipation. The results reveal that EFG at In site increases as α parameter reduces from 0.2 to 0.1 using the B3PW91 functional. Furthermore, the LDA functional with lowest degree of 4f localization yields maximum value of EFG among all considered functionals. Moreover, the previous calculations show that EFG at In site using the PBE-GGA is larger than that of PBE-GGA + U. These results imply that the reduction of 4f localization causes the EFG at In site to increase. To clarify this point, we calculate the EFG at In site with respect to the α parameter using the B3PW91 functional, see Fig. 6(a). This figure shows that the EFG increases as the α parameter decreases. This   Table 5. The main component of EFG, V zz , besides its valence, V zz val , and lattice, V zz lat , components in 10 21 Vm −2 unit at Ce and In sites using several exchange-correlation functionals (XCFs). a Reference [11]. b Reference [30]. c Reference [53]. Calculations in the present work are done using the equilibrium lattice constant a = 4.68 Å as evaluated by B3PW91with α = 0.2.
figure also shows that the EFG is dominated by the p-electrons. The anisotropy functions of Δ p(E F ) and Δ d(E F ) are also calculated using the three mentioned functionals at P = 0 and P = 14 GPa. The results are tabulated in Table 6. As evidently shown in this table, the magnitude of Δ p(E F ) is one order of magnitude greater than that of Δ d(E F ) at In site. This result is consistent with the previous calculations 30 .
Finally, we have calculated the EFG at In site using the B3PW91 with α = 0.2 as a function of pressure, as shown in Fig. 6(b). The results show that EFG increases versus pressure in a bowing shape. This bowing behavior of EFG was also previously seen using the PBE-GGA and WC-GGA 11 . However, as discussed in the previous sections, pressure can change the degree of localization. Hence, it is not appropriate to plot a sensitive quantity such as EFG versus pressure without paying attention to the changes made (by pressure) in the degree of localization. To this end, we have also plotted the EFG curve versus pressure by selecting more appropriate functionals point by point as pressure changes. Thus, EFG is calculated using the B3PW91 with α = 0.2 at zero pressure, and using the B3PW91 with α = 0.1 at P = 2 GPa, as well as using the LDA at P = 5 and 10 GPa. As it can be clearly seen from Fig. 6(b), the behavior of the EFG curve obtained using a variety of functionals is qualitatively and quantitatively different from those obtained using only a single functional. We close this section by summarizing that a specific functional can be suitable only at a special range of pressure for the strongly correlated system. Summary and concluding remarks. We have used the CeIn 3 compound as an appropriate sample to show how the degree of localization can be strongly affected by pressure in a highly correlated system within the density functional theory. To this end, we have systematically calculated the effects of pressure and 4f localization degree on the structural, magnetic, and electronic properties of this compound. The calculations were performed employing a variety of localized and itinerant functionals versus pressure. Our PBE-GGA results in agreement with the previous results reported by the others predict that the total magnetic moment per cerium atom (TOT) is suppressed for this compound at a pressure of about 16 GPa. However, this prediction of PBE-GGA is very far from the experimental value of 2.5 GPa. This inconsistency between theory and experiment can be improved using LDA functional. The LDA functional yields better results for the suppressed pressure in closer agreement with experiment than those of the other considered functionals. The Fermi surface calculations using the LDA functional are also found to be consistent with the experimental measurements. In addition to the aforementioned inconsistency, the PBE-GGA functional also fails to satisfactorily predict the TOT at zero pressure. The PBE-GGA predicts that the TOT of the compound is 0.24 μ B at zero pressure, which is inconsistent with the experimental value of 0.65 μ B . In this case, the LDA functional not only dose not improve the predicted TOT at zero pressure, but decreases the accuracy compared to that of PBE-GGA. Our results show that high localized B3PW91 functional with α = 0.2 can well predict the TOT of the system at zero pressure (0.67 μ B ) in better agreement with the experimental value (0.65 μ B ). Lattice parameter and bulk modulus of the system are systematically studied using the GGA + U and B3PW91 functionals by varying their corresponding U eff and α parameters, respectively. Our GGA + U calculations show that the bulk modulus of the system unreliably fluctuates as U eff parameter increases. This behavior which does not allow to select an appropriate value for the U eff indicates that the DFT + U scheme is rather unreliable. The hybrid B3PW91 approach provides much better agreement with the experimental data than the GGA + U scheme. , shows that the best lattice parameter is obtained using a larger value of α (0.4) than that of usually considered for the AFM phase, see Table 1. Our results show that the GGA + U fails to produce acceptable results by increasing the U eff parameter since not only its predicted bulk modulus unpredictably fluctuates by increasing U eff but also the fluctuations are performed around a value which are far from the experimental bulk modulus. However, it is worth mentioning that the hybrid B3PW91 gives reliable results for the AFM phase of the system using the larger value of α = 0.3. The latter result together with the excellent prediction of the lattice parameter by even a larger value of α (0.4) in agreement with the previous study of CuCl 2 compound 34 confirms that the percentage of HF would be increased for obtaining more reliable results. Therefore, it is crucial to calibrate the amount of HF exchange for predicting reliable results. By following this strategy after a systematic study we found that bulk modulus of the system in its AFM phase is well predicted by α = 0.2 and 0.3. By considering both of the lattice parameter and bulk modulus results and calibrating them to the corresponding experimental data we found that the B3PW91 with α = 0.2 can predict the structural properties more consistent with the experimental data. The agreement with the experimental results is acceptable with the standard B3PW91 approach for the FM phase. It is worth mentioning that these consistencies can only partially validate the correctness of the reported results. Furthermore, the Fermi surface calculations using this functional are in agreement with experimental measurements for the AFM phase. In addition, the calculated electric field gradients at In sites by this hybrid functional are found to be very closer to the experimental data at zero pressure compared to the other considered functionals. However, this hybrid functional makes the predicted suppressed pressure worse than that of PBE-GGA compared to the experimental result. Thus, the high localized B3PW91 functional with α = 0.2 is suitable functional in predicting the CeIn 3 properties at zero pressure while the LDA is a suitable functional for high pressures. This means that the character of 4f electrons is changed by applying the pressure in this compound in agreement with the dHvA experiments. Therefore, the degree of 4f localization must be tuned in terms of the pressure. High localization is suitable for low pressures and vice versa; viz., low localization is suitable for high pressures. Based on these results we conclude that the GGA failure originates from the degree of 4f localization which is considered in GGA functional. The GGA is not low localized enough to become suitable for studying the CeIn 3 at high pressures. The 4f-electrons are more localized using the GGA than those of LDA and thereby GGA predicts the critical pressure much larger than those of experiment and LDA. On the other hand, GGA is not highly localized enough to become suitable for studying the CeIn 3 at very low pressures around the P = 0. Thus, the TOT is predicted by the PBE-GGA functional at zero pressure to be smaller than those of the experiment and B3PW91 functional with α = 0.2. Our electronic structure, including DOS, band structure, and Fermi surface analyses reveal that if a fixed exchange-correlation functional is used for every pressure, then the effects of pressure on the electronic structure will not be well considered for the highly correlated systems. This is in contrast to what happens in nature, because if the degree of 4f-electron localization decreases, the electronic structure can be changed. Our results show that the GGA functional can be considered as an appropriate functional for studying the CeIn 3 compound at intermediate pressures but not for every pressure regime. In essence, this state-of-the-art DFT study demonstrates that the variation of a physical quantity with respect to pressure can sensitively depend on the used exchange-correlation functional in such a system. This in turn implies that the behavior of a physical quantity as a function of pressure can be valid within DFT only over a specific pressure range. Therefore, a physical quantity-versus-pressure curve cannot always be obtained reliably by only one currently available exchange-correlation functional over an arbitrary range of pressure. Instead, it would be more reliable to first divide the arbitrary pressure interval into some specific pressure intervals. Then, the degree of localization is determined for each pressure interval. Eventually, based on the latter determination the best exchange-correlation functional is selected and/or tuned for every pressure interval. Otherwise, consistency between experimental and theoretical DFT results using a specifically fixed exchange-correlation functional cannot be guaranteed for any pressure interval so that even one may expect to encounter an extraordinary large discrepancy between theory and experiment. By this strategy, we have explained why there are two dramatic discrepancies between previous reliable experimental and accurate ab initio DFT results in CeIn 3 and discussed how these discrepancies can be resolved.