Quantifying the stability of the anion ordering in SrVO₂H

As a computational model, we considered possible anion configurations in Sr₈V₈O₁₆H₈, by replacing one third of oxygens with hydrogens in the $2\times 2\times 2$ supercell of SrVO₃, and optimized their structures in the way described in Sec. II.2. Because of a formidable number of possible anion configurations even for our relatively small-sized simulation cell, we restricted our analysis to the configurations satisfying that each vanadium atom has two nearest-neighboring hydrogen atoms. This constraint is energetically natural and often applied for analysis of anion (dis)ordering in materials with $M$O${}_{4}X_{2}$ octahedra ($X$: anion other than oxygen) corr_disorder . By using this assumption, the number of possible (symmetrically-inequivalent) crystal structures becomes around 250.

We focused on the following two quantities to characterize anion configurations. One is the number of trans configurations of hydrogens around vanadiums, $n_{trans}$, which is up to eight since our supercell includes eight vanadium atoms. The other one is the number of a pair of two hydrogens with specific relative positions as shown in Fig. 1(c) with a solid circle, $n_{neighbor}$. Note that the atoms used for defining $n_{neighbor}$ are two vanadiums and two hydrogens surrounded with the solid circle in Fig. 1(c), and the other atoms shown in this figure are irrelevant to the definition. In other words, $n_{neighbor}$ is the number of the square-shaped relative position of these four atoms: V-V-H-H. Here, we avoid double-counting pairs of the equivalent hydrogens due to the periodic boundary condition, and then, e,g., $n_{neighbor}$ is eight for the $2\times 2\times 2$ supercell constructed from the experimental structure of SrVO₂H shown in Fig. 1(a).

These two quantities were chosen by the following reason. First, we can represent the hydrogen configuration in Sr₈V₈O₁₆H₈ by a vector ${\bm{x}}=\{x_{1},x_{2},\dots,x_{24}\}$ where $x_{i}$ denotes whether the $i$-th oxygen site in Sr₈V₈O₂₄ (i.e., a $2\times 2\times 2$ supercell of SrVO₃) is replaced with hydrogen ($x_{i}=1$) or not ($x_{i}=0$). One can also uniquely specify a hydrogen configuration by a set of the values of all the $k$-th monomial $x_{i_{1}}x_{i_{2}}\dots x_{i_{k}}$ up to 8th order, where $k$-th monomial corresponds to the $k$-body correlation among hydrogen (e.g., $x_{5}x_{7}x_{9}=1$ means that hydrogens occupy the 5th, 7th, and 9th anion sites). Thus, we expect the relation,

$$E({\bm{x}})=E_{0}+\sum_{i}c_{i}x_{i}+\sum_{i<j}c_{ij}x_{i}x_{j}+\dots,$$

(1)

where $E({\bm{x}})$ is the total energy with the hydrogen configuration ${\bm{x}}$ and $E_{0},c_{i},c_{ij},\dots$ are unknown parameters (coefficients). This is the way used in the popular cluster-expansion method cluster , where atomic species in multicomponent systems such as alloy are denoted by a discrete variable and the energy of the system is represented with polynomials of those variables. There is no need to consider terms including $x_{i}^{n}$ ($n\geq 2$) since $x_{i}^{n}=x_{i}$. Because a small supercell was used in this study under the periodic boundary condition, we should concentrate on short-range correlations represented with a low-degree monomial. There is no first order monomial that can be used to characterize anion configurations. This is because the equivalency among all the oxygen sites in cubic SrVO₃ requires that not $x_{i}$ but $\sum_{i=1}^{24}x_{i}$ should be used to represent the total energy of this system (to say, all the coefficients $c_{i}$ in $E({\bm{x}})$ should be the same by symmetry), but $\sum_{i=1}^{24}x_{i}$ is constant (8) in our simulation. By a similar discussion for the second order monomials, the number of equivalent hydrogen pairs $\sum_{\langle ij\rangle}x_{i}x_{j}$ can be used to represent the energy of this sytem, where the $i$-th and $j$-th hydrogens having specific relative positions that are equivalent by the crystal symmetry, are summed over. Then, only the two quantities, $n_{trans}$ and $n_{neighbor}$, are allowed if one considers monomials representing a pair occupation of hydrogens with their distance less than or equal to $a$, a lattice constant of cubic primitive cell. Note that the inter-hydrogen length for the initial cubic structure is used in explanation just for simplicity. We also note that the number of hydrogen pairs with their distance being $a/\sqrt{2}$ as shown with a dotted circle in Fig. 1(c), is nothing but the number of cis configurations, and thus is uniquely determined by $n_{trans}$ because we restrict the anion configurations to those where all the vanadiums have two adjacent hydrogens. For simplicity, we omit the higher-order monomials in the following analysis. It is an interesting future task to perform more elaborated statistical analysis such as the multiple regression analysis using many factors omitted in this study, while it might require a larger supercell without the restriction introduced in our simulation. In the latter analysis, we shall see that, even by our simplification, the total energies well correlate with the two quantities focused in this study with a clear physical interpretation.

We performed the structural optimization using the density functional theory with the Perdew-Burke-Ernzerhof parametrization of the generalized gradient approximation (PBE-GGA) PBE and the projector augmented wave (PAW) method paw as implemented in the Vienna ab initio simulation package (VASP) vasp1 ; vasp2 ; vasp3 ; vasp4 . Both cell deformation and optimization of the atomic coordinates were allowed unless noted. The following orbitals for each element are treated as core electrons in the PAW potentials: [Ar]$3d^{10}$ for Sr, [Ne] for V, and [He] for O. We started our calculation from the G-type antiferromagnetic configurations of the V-$d$ orbitals. Here, we need to include the spin polarization to describe the $d^{2}$ occupation into the crystal field shown in Fig. 1(b), which was pointed out to be a key to realize the trans preference of this system in the previous studies as mentioned in the introduction. We included the $+U$ correction DFTU1 ; DFTU2 with $U-J=3$ eV for the V-$d$ orbitals, which is a typical value for the $3d$ orbitals and is also consistent with our first-principles estimation by constrained random-phase approximation cRPAochi . Here, we used the simplified rotationally invariant approach to the DFT$+U$ method introduced by Dudarev et al. DFTU2 , where the correction term only depends on $U-J$ rather than the values of $U$ and $J$ themselves. Note that SrVO₂H becomes gapless in simple PBE-GGA calculation without the $+U$ correction even when the spin polarization is taken into account, owing to the well-known underestimation of the band gap in PBE-GGA, while the system is insulating in experiments SVOH_Angew . Further discussion on our choice of the $U-J$ parameter is presented in Sec. III.1. While our simple DFT$+U$ calculation cannot describe the electron correlation effects completely, we consider that our simulation captures important factors determining the anion ordering as we shall see: a large structural change by introduction of hydrogen and the different crystal fields between trans and cis configurations. This is to some extent guaranteed by the spin-polarized insulating state observed in experimental studies at low temperature, where such a static approximation of the electron-electron interaction is somewhat validated. Electron correlation effects beyond our treatment, e.g., renormalization of the size of the crystal-field splitting, is an interesting issue but requires high computational cost and then out of scope of this study. The initial crystal structure was generated in the way described before with a small random displacement, around $0.01$ in the crystal coordinate of the supercell for each direction, added to break the symmetry before the structural optimization. The plane-wave cutoff energy of 550 eV and a $6\times 6\times 6$ $\bm{k}$-mesh were used. We applied Gaussian smearing with a smearing width of 0.1 eV. Structural optimization was continued until the Hellmann-Feynman force on each atom becomes less than 0.01 eVÅ^-1 for each direction. The spin-orbit coupling was not taken into account throughout this paper.

Figure 2(a) presents the total energies of computed structures plotted against $n_{trans}$ colored using a value of $n_{neighbor}$. Here, the two crystal structures with all-trans configurations called structures A and A^′ in this paper are also shown. We define $\Delta$ as the total energy of the most stable structure except A and A^′, relative to that for structure A. The corresponding crystal structure with $n_{trans}=4$, called structure B in this paper, is shown in Fig. 2(b). Note that structure B is not used for calculating $\Delta$ under different conditions where another structure with $n_{trans}\neq 8$ is more stable than structure B, as shown in the following sections. Because of the periodic boundary condition applied to our small supercell, $n_{trans}$ always takes an even value in our simulation. We find that structure A has the lowest energy in our simulation, which is consistent with experiments. The total energies of the calculated configurations tend to be lowered by increasing $n_{trans}$. However, the all-trans structure A^′ has a sizably higher energy than that for structure A. This means that the explanation of the anion ordering in the previous studies based on the stable $d^{2}$ crystal field for the trans configurations, is partially correct from the viewpoint of the negative correlation between the total energy and $n_{trans}$ as mentioned above but insufficient to describe the stability of anion ordering in this system. We shall come back to detailed analysis on the $n_{trans}$ dependence later in this paper, but first we focus on the difference between structures A and A^′, which should be a key to understand what is missing in the description by $n_{trans}$.

One noticeable difference between structures A and A^′ is whether hydrogens are placed onto the same ($ab$) plane or not, which can be captured by $n_{neighbor}$. The colors of data points in Fig. 2(a) show that not only these two structures but also the other data points exhibit a tendency that a large $n_{neighbor}$ stabilizes the energy. In Fig. 2(c), the total energy is plotted against $n_{neighbor}$ instead of $n_{trans}$. For each set of data points with the same $n_{trans}$, we can verify the negative correlation between the total energy and $n_{neighbor}$. A possible interpretation is as follows. It is known that the size of the H^- and O^2- ions are different and thus the lattice constants of SrVO₂H are quite anisotropic: $a=3.9290(6)$Å and $c=3.6569(5)$Å  in experiment SVOH_Angew . In other words, introduction of hydrogen inclines the crystal to shrink along the V-H-V direction. In structure A, the hydrogen occupation along the same direction allows the crystal to coherently shrink along the same direction ($c$ direction), while such coherent shrinkage is not allowed in structure A^′. The coherent shrinkage by concomitant occupation of the neighboring hydrogen sites takes place when $n_{neighbor}$ is large. We note that an existence of a hydrogen pair captured by $n_{neighbor}$ in the present small supercell means that hydrogens occupy anion sites in a row with an infinite length. Thus, here we see that, not only a two-dimensional occupation of hydrogen in structure A, i.e., all oxygens in some $ab$ planes are fully replaced with hydrogens, but also a one-dimensional occupation of hydrogen as captured by $n_{neighbor}$ effectively lowers the total energy. While the small supercell used in this study is only capable of representing limited patterns of hydrogen configurations, we can partially see how the ordering of hydrogen occupation stabilizes the crystal structure.

Before moving on to further analysis, to check the robustness of our conclusion against the $U-J$ value, we calculated the total energy difference between structures A, A^′, and B, using $U-J=5$ eV. For this purpose, we did not optimize the crystal structure but just used the optimized structure using $U-J=3$ eV. As a result, we find that the energy difference between structures A and B is $1.10$ eV for $U-J=5$ eV while it is $1.13$ eV for $U-J=3$ eV, and that between structures A and A^′ is $0.97$ eV for $U-J=5$ eV while it is $0.99$ eV for $U-J=3$ eV. Therefore, we believe that our conclusion is less affected by the choice of the $U-J$ value in simulation.

To verify our view, we also calculated the total energies using a fixed cubic lattice, with $a=b=c=$ 3.864 Å and 3.999 Å for Figs. 3(a) and (b), respectively, while atomic coordinates were relaxed. Although the overall correlation between the total energies and $n_{trans}$, $n_{neighbor}$ are similar to Fig. 1(a), we find that the energy difference between structures A and A^′ are drastically decreased: from $0.99$ eV in Fig. 2(a) to $0.29$ and $0.16$ eV in Figs. 3(a)–(b), respectively. This is because the tetragonal deformation, i.e., $a=b\neq c$, is preferred for structure A but is not allowed in the cubic lattice. It is interesting that first-principles calculation of unsynthesized KTiO₂H, although with restricted number of sampled structures, also pointed out that the different size between O^2- and H${-}$ ions and a resulting shrinkage of the lattice along Ti-H-Ti are important in stabilizing structure A Tsuneyuki , where the crystal-field effect cannot take place for Ti-$d^{0}$ occupation therein.

Note that, structure A with the cubic lattice constants exhibits an octahedral rotation with the V-O-V angle of 168.5 and 170.1 degrees for Figs. 3(a) and (b), respectively, unlike the that for the relaxed tetragonal lattice shown in Fig. 2(a). It is interesting that the lowest energy among the all-cis configurations with $n_{trans}=0$ relative to the total energy of structure A becomes much lower in Fig. 3(a) ($0.20$ eV) than in Fig. 3(b) ($0.36$ eV). A definitive conclusion cannot be reached within our limited calculation data, but it seems that a smaller ionic radius of H^- alleviates a cramped environment for atoms, and then stabilizes the crystal structure with the cis configurations, while it is not so active for structure A with VO₂ planes without hydrogen. In the experimental study on SrCrO₂H, it was pointed out that the introduction of hydrogen into SrCrO₃ makes its tolerance factor closer to unity and then raises the Néel temperature CrOH . It is an interesting future issue to investigate the effect of the octahedral rotation and tilt onto the stability of anion configuration, requiring some quantities to characterize such kinds of local distortion. Here we just point out that this issue has some relevance to the coherent shrinkage captured by $n_{neighbor}$ dependence of the total energy and by the energy difference between structures A and A^′: both of them originate from the different anion size and the resultant structural change.

To get insight into the $n_{trans}$ dependency, we performed the same simulation as Fig. 2(a) but with different numbers of electrons. This is because the energy stabilization by the trans crystal field should be affected by a change of the number of electrons. Figures 4(a) and (b) present the total energies calculated using the number of electrons increased or decreased by two in the supercell, respectively. In other words, the electron occupation in the V-$d$ bands is $d^{1.75}$ and $d^{2.25}$ in average, for Figs. 4(a) and (b), respectively. For these calculation, the same amount of the background positive charge was introduced. We find that $\Delta$ is decreased from $1.13$ eV in Fig. 2(a) to $0.61$ and $0.81$ eV in Figs. 4(a)–(b), respectively. The enhanced stability of the trans configuration in the $d^{2}$ occupation pointed out and tested for some limited anion configurations in previous studies, is verified by our simulation for more intensive set of anion configurations. In fact, experimental studies showed that hydrogens are randomly distributed in $A$TiO_3-xH_x ($A=$ Ca, Sr, Ba, Eu) TiOH1 ; TiOH2 ; TiOH3 ; TiOH4 , SrV_1-xTi_xO_1.5H_1.5 SVTiOH , and SrCrO₂H CrOH , where the number of $d$ electrons are increased or decreased from SrVO₂H. Note that $\Delta$ can be decreased if one uses a larger supercell that can represent anion configurations not representable in our simulation. In this sense, $\Delta$ estimated in our study is its upper bound, and so the stability of the trans configuration against the change of the electron number can be lost more easily than our simulation.

Here we point out that a fewer electron occupation than $d^{2}$ can induce a Jahn-Teller instability in the trans crystal field that makes $d_{xz}$ and $d_{yz}$ inequivalent, and in fact it takes place in our simulation: V-O length becomes different by up to 0.09 Å along the $a$ and $b$ directions in the optimized crystal structure. It is interesting that a small reduction of electrons can provide an interesting platform of Jahn-Teller physics, while reducing a large number of electron will destabilize the crystal structures with the trans configuration. We note that a Jahn-Teller instability can take place also for the $d^{n}$ ($1<n<3$) occupation in the cis configuration, but the crystal structure with cis configurations has a low symmetry in general even without the Jahn-Teller instability, and so it is difficult to quantify such an instability for the cis configurations.

To verify our view on the crystal-field effect, we here make an alternative estimation of the energy stabilization by the crystal-field splitting. In our previous study cRPAochi , we evaluated the size of the crystal-field splitting between the $d_{xy}$ and $d_{xz/yz}$ orbitals, $\Delta_{cf}=0.45$ eV, by constructing the Wannier functions Wannier1 ; Wannier2 . Suppose for simplicity that the energy levels of the $t_{2g}$ orbitals in the VO₄H₂ octahedron are simply determined by the number of hydrogens on the plane where the orbitals are extended (e.g., the $xy$ plane for the $d_{xy}$ orbital), the energy diagrams for the trans and cis configurations are obtained as shown in Fig. 4(c). Therefore, the energy difference between the trans and cis configurations with $d^{2}$ occupation is around $\Delta_{cf}/2=0.23$ eV per vanadium. This estimate corresponds to $-0.23$ $n_{trans}$ (eV) dependence of the total energy in the supercell, which amounts to the energy stabilization of $1.9$ eV for $n_{trans}=8$ compared with $n_{trans}=0$. This is roughly consistent with the energy separation among the data set of $n_{trans}=0$ and $n_{trans}=8$ at $n_{neighbor}=4$ and $n_{neighbor}=8$ in Fig. 2(c). On the other hand, at $n_{neighbor}=4$ in Fig. 2(c), the energy separation between the data set of $n_{trans}=2$ and $n_{trans}=4$, and that between $n_{trans}=4$ and $n_{trans}=8$, are around 0.2 eV, which is a few times smaller than our rough estimate. One possible reason is that, as we have seen, the introduction of hydrogen can induce a large structural change, which might alter the size of the crystal-field splitting from our rough estimate. In addition, we can do a similar estimate of how $\Delta$ changes by adding or removing two electrons from or into the supercell. By considering the cis and trans crystal fields shown in Figs. 4(c)–(d), the resultant change of $\Delta$ is $\Delta_{cf}=0.45$ eV for the both cases (i.e., making $d^{1}$ or $d^{3}$ filling for two vanadiums as we verified in the crystal structures used for calculating $\Delta$). This rough estimate is again consistent with our simulation for the change of $\Delta$, $0.52$ eV for Fig. 4(a) and $0.32$ eV for Fig. 4(b). These observations, to some extent, validate our interpretation regarding the crystal field effect. Because the introduction of hydrogen can induce a large structural change, further investigation using more extensive data sets and characteristic quantities in a larger supercell would be an important issue to distinguish the crystal field effect and others more clearly, while the trans-preference is clearly verified in our simulation.

Although both our simulation and experiment show that the anion-ordered structure A is the most stable, random anion distribution can be to some extent stabilized by entropy. Here, we provide a rough estimate on the effects of possible two kinds of entropy as follows.

First, the configuration entropy for fully random anion configuration is given by $S_{\mathrm{conf}}=-k_{B}\sum_{i}x_{i}\ln x_{i}$ per single anion site, where $x_{i}$ is the fraction of the anion $i$, i.e., $x_{i}=2/3$ for oxygen and $1/3$ for hydrogen in our case. Thus, an energetic gain $-TS_{\mathrm{conf}}$ for fully random distribution compared with structure A with perfect anion ordering is $-390$ meV for the $2\times 2\times 2$ supercell at $T=300$ K, which is much smaller than $\Delta$ in Fig. 2(a) (1.13 eV). Note that the structures with random anion configurations have, in average, a larger energy by around 2 eV than structure A as shown in Fig. 2(a).

Second, the $d^{2}$ occupation in the cis configuration realized for random anion configuration can yield the electronic entropy. To evaluate the upper bound of the electronic entropy effect, we assume that all the vanadiums have the cis configuration of hydrogens without a sufficiently large Jahn-Teller instability that breaks the degeneracy of the energy levels. We also assume that the G-type antiferromagnetic ordering with a high-spin state $S=1$ is kept, where $d^{2}$ occupation is two-fold degenerate in the $t_{2g}$ manifold of the cis configuration (see Fig. 4(c)). Then, the electronic entropy is given by $S_{el}=k_{B}\ln 2$ per vanadium, resulting in $-143$ meV stabilization for the $2\times 2\times 2$ supercell at $T=300$ K. This is in the same order of magnitude but smaller than the stabilization by the configuration entropy described above. While we assume the high-spin state ($S=1$) here, note that the low-spin state ($S=0$) rather stabilizes the all-trans configurations in terms of the electronic entropy effect. Additional entropy released for the paramagnetic state is not considered here because this is also acquired for the all-trans configurations above the transition temperature.

By considering the two types of entropies here, we can conclude that the entropy effects can be regarded as a relatively small correction to the total energy in SrVO₂H, although this effect can be important when $\Delta$ is small, e.g., in the case of different transition metal species and/or at high temperatures.

It is often the case that a relatively large zero-point vibrational energy due to a light mass of hydrogen can affect the stability of the crystal structure including hydrogen. To check how the zero-point vibrational energy of hydrogens affects our conclusion, we performed phonon calculation for structures A and B because these two structures determine the value of $\Delta$ for Figs. 2(a) and Figs. 4(a)–(b), while the lowest-energy structure with $n_{trans}\neq 8$ is different from structure B in Fig. 3. For simplicity, we calculated the $\Gamma$-point phonon frequencies to give a rough estimate of the zero-point vibrational energy: $\sum_{i}\hbar\omega_{i}/2$, where $i$ and $\omega_{i}$ denote the index of the phonon mode at the $\Gamma$ point and its frequency, respectively. For this purpose, we employed the finite displacement method as implemented in the Phonopy phonopy software in combination with VASP. Other computational conditions are the same as those for the structural optimization presented in Sec. II.2.

Figure 5 presents the calculated phonon frequencies at the $\Gamma$ point for structures A and B. While the phonon frequencies are similar as a whole between the two structures, we find that the eight highest frequencies exhibit a sizable change of around 20 meV. As a result, the zero-point vibrational energy, $\sum_{i}\hbar\omega_{i}/2$, for these eight phonon modes in structure B is smaller by 79.6 meV than that for structure A. Because these high-frequency phonon modes represent the hydrogen moving toward the V-H-V direction, such an increase of the phonon frequency should take place when the V-H-V length becomes shorter. Therefore, the energy stabilization by the shrinkage of the V-H-V length that can be captured by $n_{neighbor}$, is partially canceled by the increase of the zero-point vibrational energy. However, the effect of the zero-point vibrational energy seems to be smaller than the energy profits discussed so far, hence does not affect our conclusion.