Van Vleck Analysis of Angularly Distorted Octahedra using VanVleckCalculator

Liam. A. V.Nagle-Cocco[Email<EMAIL_ADDRESS>ORCID: 0000-0001-9265-1588.]
Siân E.Dutton[Email<EMAIL_ADDRESS>ORCID: 0000-0003-0984-5504.]

Cavendish Laboratory, University of Cambridge, JJ Thomson Avenue, Cambridge, CB3 0HE, United Kingdom

A method and associated Python script, VanVleckCalculator, is described for parametrising octahedral shear and first-order Jahn-Teller distortions in crystal structures.

Van Vleck modes describe all possible displacements of octahedrally-coordinated ligands about a core atom. They are a useful analytical tool for analysing the distortion of octahedra, particularly for the first-order Jahn-Teller distortion.
Determination of Van Vleck modes of an octahedron is complicated by the presence of angular distortion of octahedra however. This problem is most commonly resolved by calculating the bond distortion modes ($Q_2$, $Q_3$) along the bond axes of the octahedron, disregarding the angular distortion and losing information on the octahedral shear modes ($Q_4$, $Q_5$, and $Q_6$) in the process. In this paper, the validity of assuming bond lengths to be orthogonal in order to calculate the van Vleck modes is discussed, and a method is described for calculating Van Vleck modes without disregarding the angular distortion. A Python code for doing this, VanVleckCalculator, is introduced, and some examples of its use are given. Finally, we show that octahedral shear and angular distortion are often, but not always, correlated, and propose a parameter as the shear fraction, $\eta$. We demonstrate that $\eta$ can be used to predict whether the values will be correlated when varying a tuning parameter such as temperature or pressure.

§ INTRODUCTION

The van Vleck distortion modes <cit.> modes describe all possible displacements of octahedrally-coordinated ligands about a core atom. They are particularly useful in the context of the Jahn-Teller effect <cit.>, which in general occurs when a high-symmetry coordination is destabilised with respect to a deviation to lower symmetry as a consequence of electronic degeneracy. The Jahn-Teller effect distorts the crystal structure via the Jahn-Teller distortion. While the Jahn-Teller distortion is not unique to octahedra in bulk crystalline materials, it is in octahedra that it was first observed experimentally <cit.>, and it is in materials with Jahn-Teller-distorted octahedra that colossal magnetoresistance <cit.> and high-temperature superconductivity <cit.> were discovered.

A transition metal (TM) cation in an octahedral configuration will have its $d$ orbitals split into three $t_{2g}$ orbitals[In this paper, we use the notation that lower case symmetry descriptors (such as $e_g$ or $t_{2g}$) refer to orbitals with this symmetry, and upper case descriptors (such as $E_g$ or $T_{2g}$) refer to the symmetry more generally.] at lower energy and two $e_g$ orbitals at higher energy. It will have a number, $n$, of electrons in these $d$ orbitals (hereafter described as $d^n$). For certain values of $n$ and, where applicable, certain low- or high-spin characters[In the low-spin case, $t_{2g}$ orbitals fill fully before $e_g$ orbitals gain electrons; in the high-spin case, once the $t_{2g}$ orbitals are singly-occupied, the next two electrons will populate the $e_g$ orbitals.], there will exist multiple orbitals that could be occupied by an electron or an electron hole with equal energy. This degeneracy is destabilising, resulting in the most stable configuration of atomic sites being one in which the ligands distort from their high-symmetry positions in order to rearrange the orbitals into a non-degenerate system with minimised energy. This is shown for a low-spin $d^7$ TM cation (such as Ni$^{3+}$ or Co$^{2+}$) in Figure <ref>, though such distortions may occur for any value of $n$ in $d^n$ where there is a degenerate occupancy. The stabilisation energy due to the Jahn-Teller effect is larger for $e_g$ degeneracy than $t_{2g}$ degeneracy, and so the effect is prominent to higher temperatures, and hence more widely-studied, in JT-active materials with $e_g$ degeneracy <cit.>.

The orbital rearrangement due to a tetragonal elongation for an octahedrally-coordinated low-spin $d^7$ transition metal ion, which typically occurs due to the first-order Jahn-Teller effect.

In the literature, various techniques for parameterising the Jahn-Teller distortion are used. An often-used example <cit.> is the bond length distortion index, defined by baur1974geometry , as:

\begin{equation}\label{BLDI_equation}
D=\frac{1}{n} \sum_{i=1}^n \frac{|l_i - l_\mathrm{av}|}{l_\mathrm{av}}
\end{equation}

where $l_i$ is the distance between the core ion and the $i$th coordinated ion, and $l_\mathrm{av}$ is the average of all the distances between the core ion and coordinated ions.

A similar parameter <cit.> is the effective coordination number, which for an octahedron deviates from 6 only when there is bond length distortion, defined by hoppe1979effective  as:

\begin{equation}\label{ECoN_equation}
\mathrm{ECoN} = \sum_{i=1}^n \exp \left[ 1 - \left(\frac{l_i}{l'_\mathrm{av}} \right)^6 \right]
\end{equation}

where $l'_\mathrm{av}$ is a modified average distance defined as:

\begin{equation}
l'_\mathrm{av} =
\frac{  \sum_{i=1}^n l_i \exp \left[ 1 - \left(\frac{l_i}{l_\mathrm{min}} \right)^6  \right]}
{ \sum_{i=1}^n   \exp \left[ 1 - \left(\frac{l_i}{l_\mathrm{min}} \right)^6 \right]}
\end{equation}

Finally, a third parameter used to quantify the Jahn-Teller distortion <cit.> is the quadratic elongation, $<\lambda >$, defined by robinson1971quadratic  as:

\begin{equation}\label{QuadElon_equation}
<\lambda > = \frac{1}{n} \sum_{i=1}^n \left(\frac{l_i}{l_0}\right)^2
\end{equation}

where $l_0$ is the centre-to-vertex distance of a regular polyhedron of the same volume.

More recently, an alternative approach to modelling polyhedral distortion has been described <cit.>, involving fitting an ellipsoid to the positions of the ligands around a coordination polyhedron, calculating the three principal axes of the ellipsoid, $R_1$, $R_2$, and $R_3$, where $R_1\le R_2 \le R_3$, and using the variance of these three radii as a metric for the distortion. This has been applied to the first-order Jahn-Teller distortion in pughe2023partitioning .

These parameterisations each have merits. However, they are not sensitive to the symmetry of the octahedral distortion. The van Vleck modes are conceptually different to each of these for quantifying the Jahn-Teller distortion because they can be used to quantify distortion with the precise symmetry of the transition metal $e_g$ orbitals. This is important because Jahn-Teller distortions typically follow a particular symmetry. When the distortion is due to degeneracy in the $e_g$ orbitals it will be of $E_g$ symmetry; when it is due to degeneracy in the $t_{2g}$-degenerate orbitals it may be either $E_g$ or $T_{2g}$ symmetry <cit.>, although there is relatively little unambiguous experimental evidence for a Jahn-Teller-induced shear as compared with more typical $E_g$ distortion..

In this paper, we present a Python <cit.> package, VanVleckCalculator, for calculating the van Vleck distortion modes. We show that the approach to calculating the modes which is commonly used in the literature is a reasonable approximation for octahedra with negligible angular distortion, but results in the loss of information in other cases. We propose a new metric, the shear fraction $\eta$, for understanding the correlation between octahedral shear and angular distortion. Finally, we re-analyse some previously-published data in terms of the van Vleck modes to show that these can be an effective way of understanding octahedral behaviour.

The 6 van Vleck modes exhibited for an octahedron, with sites labelled using the notation in the Theory section. For the octahedra exhibiting $Q_1$, $Q_2$, and $Q_3$ distortions, there is no angular distortion; for the octahedra exhibiting $Q_4$, $Q_5$, and $Q_6$ distortions, there is no bond length distortion. For the octahedral shear ($Q_4$, $Q_5$, and $Q_6$) modes, axes are drawn to show where the bond directions would be if undistorted. An octahedron can exhibit several, or all, of these distortions simultaneously.

§ THEORY

Within an octahedron, we can split the 6 ligand ions into three pairs, where the two ions within the pair are opposite one another. In the absence of angular distortion (i.e., assuming all ligand-core-ligand angles are an integer number of 90$^\circ$), there would exist a basis where each of the three axes exist directly along the $x$-, $y$-, and $z$-axis, and where the origin in space is defined as the centre of the octahedron.

Each pair within an octahedron can therefore be assigned to an axis and labelled as the $a$, $b$, or $c$ pair respectively. Within a pair, ions can be labelled as $-$ or $+$ depending on whether they occur at a negative or positive displacement from the origin, along the axis, respectively. This notation is demonstrated in Figure <ref>, where each pair of ions is represented by a different colour.

For each of the 6 ligands, we define a set of coordinates: $x^\alpha_\beta$, $y^\alpha_\beta$, and $z^\alpha_\beta$, where $\alpha$ is $a$, $b$, or $c$ denoting the pair in which the ligand is, and $\beta$ is $-$ or $+$ denoting which ion within the pair.

The ideal positions of the six ligands are: ($R$,0,0), ($-R$,0,0), (0,$R$,0), (0,$-R$,0), (0,0,$R$), and (0,0,$-R$), where $R$ is defined as the distance between the centre of the octahedron and the ligand in an ideal octahedron (in practice, this is taken as the average of the core-ligand bond distances). This results in 18 independent variables. Using these, we further define a set of van Vleck coordinates (capitalised to distinguish from true coordinates) which is the displacement of the ion within an axis away from its ideal position. For instance, for the ion with $\alpha=a$ and $\beta=-$: $X^a_-$ = $x^a_- + R$, $Y^a_-$ = $y^a_-$, and $Z^a_-$ = $z^a_-$. See Figure <ref> for clarification of the ion notation.

Using these coordinates, the first six van Vleck modes ($Q_j$; $j=1-6$) are defined as follows <cit.>:

\begin{equation}\label{Q1_equation_definition}
Q_1 = X^a_+ - X^a_- + Y^b_+ - Y^b_- + Z^c_+ - Z^c_-
\end{equation}

\begin{equation}\label{Q2_equation_definition}
Q_2 = \frac{1}{2} \left[ X^a_+ - X^a_- - Y^b_+ + Y^b_- \right]
\end{equation}

\begin{equation} \label{Q3_equation_definition}
Q_3 = \frac{1}{\sqrt{3}} \left[ \frac{1}{2} \left( X^a_+ - X^a_- + Y^b_+ - Y^b_- \right) - Z^c_+ + Z^c_- \right]
\end{equation}

\begin{equation}\label{Q4_equation_definition}
Q_4 = \frac{1}{2} \left[ X^b_+ - X^b_- + Y^a_+ - Y^a_- \right]
\end{equation}

\begin{equation}\label{Q5_equation_definition}
Q_5 = \frac{1}{2} \left[ Z^a_+ - Z^a_- + X^c_+ - X^c_- \right]
\end{equation}

\begin{equation}\label{Q6_equation_definition}
Q_6 = \frac{1}{2} \left[ Y^c_+ - Y^c_- + Z^b_+ - Z^b_- \right]
\end{equation}

The $Q_2$-$Q_3$ phase space for elongated octahedra, with a representation of the values $\rho_0$ and $\phi$. Based on a figure from an article by Goodwin_two_temps .

We only discuss these first six van Vleck modes, which are shown in Figure <ref>. $Q_1$ to $Q_3$ describe bond length distortions, whereas $Q_4$ to $Q_6$ describe octahedral shear distortions. $Q_1$ is a simple expansion/contraction mode which does not affect symmetry and will not be discussed further.

$Q_2$ and $Q_3$ are a planar rhombic distortion and a tetragonal distortion respectively; they are considered degenerate due to the Hamiltonian, which is discussed for instance in kanamori1960crystal . These two modes form a basis for distortions describing different octahedral configurations with the symmetry of the transition metal $e_g$ orbitals <cit.>. These modes are of most relevance for first-order Jahn-Teller distortions occurring due to degenerate $e_g$ orbitals. A phase space of possible octahedral configurations can be constructed using these two parameters <cit.>, as shown in Figure <ref>. Here the magnitude of the distortion $\rho_0$ can be calculated as follows:

\begin{equation} \label{rho_0}
\rho_0 = \sqrt{{Q_2}^2+{Q_3}^2} % call it \rho_0 to be consistent with Goodenough and Zhou
\end{equation}

and the angle[Note that this angle does not represent a physical angle within the octahedron.] $\phi$ of this distortion from being of purely $Q_3$ character can be calculated by:

\begin{equation}\label{phi_equation}
\phi = \arctan{\left(\frac{Q_2}{Q_3}\right)}
\end{equation}

All possible combinations of the $Q_2$ and $Q_3$ modes correspond to a particular angle $\phi$, and hence a particular configuration as shown in Figure <ref>. The structural effect of a rotation of $\phi$ within a range of 120$^\circ$ can be quite significant, as shown in Figure <ref>; such changes can manifest as a Jahn-Teller-elongated{compressed} octahedron with 4 short{long} and 2 long{short} bonds (such as NiO$_6$ in NaNiO$_2$ <cit.>) or 2 short, 2 medium, and 2 long bonds (such as LaMnO$_3$ <cit.>).

The special angles in the $Q_2$-$Q_3$ phase space [Figure <ref>], as a function of $\phi=\arctan\left(Q_2/Q_3\right)$, with the associated singly-occupied $e_g$ orbital, for $d^4$ and low-spin $d^7$ octahedral complexes. Note that for angles which are not special angles, there will be mixing of the orbital states of the nearest two special angles.

$\phi$ ($^\circ$)     $\Psi(\phi)$

0      $d_{z^2}$

60     $d_{y^2-z^2}$

120     $d_{y^2}$

180     $d_{x^2-y^2}$

240     $d_{x^2}$

300     $d_{z^2-x^2}$

A characteristic of the Jahn-Teller distortion is that, in the absence of external distortive forces, the symmetry of the structure matches the symmetry of the orbitals involved. Typically, any $d$-orbital Jahn-Teller distortion will have some planar rhombic ($Q_2$) or tetragonal ($Q_3$) character. However, sometimes when the degeneracy occurs in the $t_{2g}$ orbital, there may instead be a trigonal component to the symmetry of the distortion, which manifests as an angular distortion instead <cit.>. For the more commonly-studied case of a degeneracy in the $e_g$ orbitals, the effect of a rotation of $\phi$ similarly changes the symmetry of the $d$ orbitals. Figure <ref> shows the splitting of the $d$ orbitals in an octahedrally-coordinated $d^7$ transition metal due to an elongation-type first-order Jahn-Teller distortion, where the tetragonal elongation occurs along the $z$-axis. Note that the unpaired $e_g$ electron occupies the $d_{z^2}$ orbital. In the opposite case of a compression-type first-order Jahn-Teller distortion along the $z$ axis, the lower-energy, and hence singly-occupied, orbital would be the $d_{x^2-y^2}$; this would correspond to a rotation in $\phi$ of 180$^\circ$. More generally, as a function of $\phi$, there exist a set of special angles separated by a 60$^\circ$ rotation corresponding to a particular $e_g$ orbital being singly-occupied by a $d$ electron. These are tabulated in Table <ref>. An octahedron for which $\phi$ does not correspond to one of these special angles exhibits orbital mixing <cit.>.

The $Q_4$ to $Q_6$ modes describe shear of the octahedra, i.e. the effect whereby paired ligands at opposite sides of a central ions are displaced in opposite directions
, and have trigonal $T_{2g}$ character. The shear modes may be used to quantify the Jahn-Teller distortion in octahedra where the degeneracy occurs within $t_{2g}$ orbitals <cit.>. The magnitude of the calculated shear is typically correlated with angular distortion, which is commonly quantified using the $\sigma_\zeta^2$ metric called the Bond Angle Variance <cit.> (BAV), defined here as:

\begin{equation}\label{BAV_equation}
\sigma^2_\zeta = \frac{1}{m-1} \sum_{i=1}^m (\zeta_i - \zeta_0)^2
\end{equation}

where $m$ is the number of bond angles (i.e. 12 for octahedra), $\zeta_i$ is the $i$th bond angle, and $\zeta_0$ is the ideal bond angle for a regular polyhedron (i.e. 90$^\circ$ for an octahedron). However, for direct comparison to the shear modes, it is more appropriate to use the standard deviation $\sigma_\zeta$.

Three possible octahedral shear/anti-shear distortions, with the associated value of the shear fraction $\eta$ as defined in Equation <ref>. In the case for $\eta=0$, the only distortion is anti-shear within a single plane. In the case for $\eta=0.5$, there are two planes in which there is distortion, a shear and anti-shear distortion equal in magnitude. In the case for $\eta=1$, there is a plane with a purely shear distortion.

For an octahedron with non-zero $T_{2g}(Q_4,Q_5,Q_6)$ modes, increasing their magnitude will increase the angular distortion, but an octahedron may have angular distortion without exhibiting octahedral shear. To analyse the extent to which angular distortion in an octahedron is due to shear, we propose a shear fraction parameter $\eta$, demonstrated in Figure <ref> and defined below.

First, we must define a set of shear and “anti-shear" angular indices, which are modifications of Equations <ref> to <ref> in terms of angles rather than displacements. The indices are represented with $\Delta$ and a subscript corresponding to the plane in which rotation occurs: the $ab$-plane corresponds to the $Q_4$ mode, the $ac$-plane to the $Q_5$ mode, and the $bc$-plane to the $Q_6$ mode. The absence or presence of a prime symbol, $\prime$, designates whether the index represents shear or anti-shear respectively. Finally, the $\delta$ angle is the rotation of the ligand from its ideal van Vleck coordinate in a clockwise direction, within the plane in which the corresponding van Vleck shear ($Q_4$ to $Q_6$) would occur. These are defined thus (see SI, Figure S7):

\begin{equation}\label{Delta_ab_eq}
\Delta_{ab} = \frac{1}{2} \left[ \delta^b_+ - \delta^b_- + \delta^a_+ - \delta^a_- \right]
\end{equation}

\begin{equation}\label{DeltaPrime_ab_eq}
\Delta^\prime_{ab} = \frac{1}{2} \left[ \delta^b_+ + \delta^b_- - \delta^a_+ - \delta^a_- \right]
\end{equation}

\begin{equation}\label{Delta_ac_eq}
\Delta_{ac}  = \frac{1}{2} \left[ \delta^a_+ - \delta^a_- + \delta^c_+ - \delta^c_- \right]
\end{equation}

\begin{equation}\label{DeltaPrime_ac_eq}
\Delta^\prime_{ac}  = \frac{1}{2} \left[ \delta^a_+ + \delta^a_- - \delta^c_+ - \delta^c_- \right]
\end{equation}

\begin{equation}\label{Delta_bc_eq}
\Delta_{bc}  = \frac{1}{2} \left[ \delta^c_+ - \delta^c_- + \delta^b_+ - \delta^b_- \right]
\end{equation}

\begin{equation}\label{Delta_prime_bc}
\Delta^\prime_{bc}  = \frac{1}{2} \left[ \delta^c_+ + \delta^c_- - \delta^b_+ - \delta^b_- \right]
\end{equation}

We then quantify the shear and “anti-shear" distortions using the following equations:

\begin{equation}\label{Deltamathrmshear2}
\Delta_\mathrm{shear}^2 = \Delta_{ab}^2 + \Delta_{ac}^2 + \Delta_{bc}^2
\end{equation}

\begin{equation}\label{Deltamathrmantishear2}
\Delta_\mathrm{anti-shear}^2 = \Delta_{ab}^{\prime2} + \Delta_{ac}^{\prime2} + \Delta_{bc}^{\prime2}
\end{equation}

From here, we define the shear fraction $\eta$ as follows:

\begin{equation}\label{eta_equation}
\eta  = \frac{\Delta_\mathrm{shear}^2}{\Delta_\mathrm{shear}^2+\Delta_\mathrm{anti-shear}^2}
\end{equation}

This $\eta$ parameter will be important in interpreting the relation between the angular distortion, $\sigma_\zeta$, and the van Vleck shear modes $Q_4$ to $Q_6$.

§ IMPLEMENTATION

In this section, the algorithm used to calculate Van Vleck distortion modes is discussed. It is written using Python 3 <cit.> as a
package called VanVleckCalculator, with the full code available on GitHub <cit.>, and also presented with annotations in the Supplementary Information.
Data handling and some calculations make use of NumPy <cit.>, and crystal structures are handled using PyMatGen <cit.>.

A flow chart showing the octahedral rotation algorithm can be found in Supplementary Information, Figure S1.

Besides calculating the van Vleck modes and the angular shear modes described in this paper, VanVleckCalculator can also calculate various other parameters as described in Supplementary Information.

§.§ Selecting an origin

Selection of the origin is a key step in calculating van Vleck modes. The most common approach, for an $M$X$_6$ octahedron, is to take the $M$ ion as the origin. This is a reasonable approach, given that $M$ ions are typically positioned at, or very close to, the centre of an octahedron. This is particularly appropriate for unit cells derived from Rietveld refinement <cit.> of Bragg diffraction data, where the $M$ ion is likely to occur at a high-symmetry Wyckoff site. A third, similar, option would be to choose the average position of the 6 ligands as the origin in space. An example of when this may be a desirable choice would be for systems exhibiting a pseudo Jahn-Teller effect (also called the second-order Jahn-Teller effect), where the central cation is offset from the centre of the octahedron.

In some instances, a crystal structure may be simulated using a supercell. Examples include so-called “big box" Pair Distribution Function (PDF) analysis <cit.> and Molecular Dynamics (MD) <cit.> simulations. Such a supercell typically retains the periodicity which is an axiom of a typical crystallographic unit cell, but will exhibit local variations. For instance, a unit cell obtained by analysis of Bragg diffraction data is typically regarded as an “average" structure, insensitive to local phenomena such as thermally-driven atomic motion or disordered atomic displacements such as a non-cooperative Jahn-Teller distortion. In a crystallographic unit cell, thermal motion of atoms is typically represented by variable Atomic Displacement Parameters (ADPs) <cit.>. In contrast, a supercell should reflect local phenomena, for instance exhibiting local Jahn-Teller distortions in a system with a non-cooperative Jahn-Teller distortion, and representing thermal effects not with ADPs but rather by distributing equivalent atoms in adjacent repeating units in slightly different positions. In this regard, a supercell can be considered a “snapshot" of a crystal system at a point in time. It may not be appropriate to set the core ion as the centre of the octahedron in a supercell, therefore, as the positioning of both core and ligand ions is in part due to thermal effects, and so the “centre" of the octahedron will be displaced due to random motion. The alternative option would be to simply use the crystallographic site of the central ion and fix this as independent of the precise motion of the central ion locally.

In VanVleckCalculator, the user has the option to take as the centre of the octahedron either the central ion, the average position of the 6 ligands, or a specified set of coordinates.

§.§ Calculating van Vleck modes along bond directions

The calculation of the van Vleck modes, as described in the Theory section, requires that the basis in space be the octahedral axes (i.e. the three orthogonal axes entering the octahedron via one vertex, passing through the central ion, and exiting via the opposite vertex). For a given crystal structure, this may require that an octahedron be rotated about each of the three axes making up the basis, until the octahedral axes perfectly align with the basis. This becomes more complicated when the octahedron exhibits angular distortion (i.e. exhibits ligand-core-ligand angles are not integer multiples of 90$^\circ$). In this case, it is impossible to define octahedral axes according to the strict criteria previously defined.

In the literature, this problem is generally evaded by simply calculating the Van Vleck modes on the basis of bond directions rather than Cartesian coordinates; for example, previous work on the perovskite LaMnO$_3$ <cit.>, other perovskites <cit.>, or non-perovskite materials <cit.>.[We note that some works use a different variation which still uses Kanamori's approximation. Papers cited here include those which use the approximation, even if the precise definitions differ.] In this case, $Q_2$ and $Q_3$ are defined according to the following equations which were first expressed by kanamori1960crystal , where $l$, $m$, and $s$ are the short, medium, and long bond lengths respectively[The equations presented here differ from Kanamori's as they have been multiplied by a factor of $\frac{\sqrt{2}}{2}$, so that they are mathematically equivalent to Equations <ref> and <ref>.]:

\begin{equation}\label{schmitt_Q2}
\textcolor{black}  {   Q_2 = l-s   }
\end{equation}

\begin{equation}\label{schmitt_Q3}
\textcolor{black}  {   Q_3 = \frac{(2m-l-s)}{\sqrt{3}}   }
\end{equation}

This relies on the implicit assumption that bond lengths are orthogonal. This is clearly a reasonable approximation in many cases, particularly when angular distortion is very small. For instance, in LaMnO$_3$, the corner-sharing octahedral connectivity enables mismatched polyhedra to tessellate via octahedral tilting [Figure <ref>(e)] rather than intra-octahedral angular distortion. However, for systems with greater angular distortion, for instance those with edge- or face-sharing interactions, it is not so clear that this approximation is valid.

§.§ Calculating van Vleck modes within Cartesian coordinates

In VanVleckCalculator, we have written an algorithm for rotating an octahedron about three Cartesian axes with a defined origin within the octahedron, such that the ligands are as close as possible to the axes (within the constraint that there is angular distortion). This allows for calculation of van Vleck modes in a way that does not artificially constrain the octahedral shear modes ($Q_4$, $Q_5$, and $Q_6$) to be zero.

First, three orthogonal axes are taken as the $x$-, $y$-, and $z$- axes[We note that, for a set of three orthogonal vectors chosen as the axes, the choice to assign each to $x$, $y$, or $z$ will not affect the value of $\rho_0$, but will affect the value of $\phi=\arctan{\left(Q_2/Q_3\right)}$ by an integer multiple of 120$^\circ$, plus a reflection about the nearest special angle (see Table <ref>) if there is $Q_2$-$Q_3$ mixing.]. By default, these are the [1,0,0], [0,1,0], and [0,0,1] axes respectively, but alternative sets of orthogonal vectors can be given by the user; for instance, for regular octahedra rotated 45$^\circ$ about the x axis, the user would be recommended to give as axes [1,0,0], [0,$\sqrt{2}$,$-\sqrt{2}$], and [0,$\sqrt{2}$,$\sqrt{2}$]. This vector is given as a Python list with shape (3,3). For consistency, the cross product of the first two axes should always be parallel with the third given vector; if anti-parallel, the algorithm will automatically multiply all elements in the third vector by -1. The three pairs of the octahedron (as defined in the Theory section) are each then assigned to one of these three axes on the basis of which pair has the largest projection of its displacement (the vector between two on a particular axis, with the $z$-axis assigned first, then the $y$-axis from amongst the two pairs not assigned to the $z$-axis, then the $x$ axis is automatically assigned to the remaining pair). Within each pair, the ligands are then ordered such that the ligand with the negative distance is along the assigned vector first, then the ligand with positive distance occurs second.

Second, the octahedron is rotated about the $x$-, $y$-, and $z$- directions of the basis repeatedly until the orthogonal axes supplied in the previous step match the basis precisely. This is performed in a while loop structure, with the rotation angles about the three axes summed in quadrature and compared with a defined tolerance (by default, $3\times10^{-4}$ radians in VanVleckCalculator), and if the total rotation exceeds the tolerance, the step is repeated[This is because rotation operations do not commute, and so a single rotation about each axis is unlikely to result in the defined axes being superimposed over the basis vectors.]. This step is usually unnecessary, and can be skipped by leaving the default set of orthogonal axes, which are [1,0,0], [0,1,0], and [0,0,1] (meaning no rotation will occur).

Third, an automatic rotation algorithm will further minimise the effect of angular distortion. For each of the three axes, the four ligands not intended to align with that axis are selected. The angle to rotate these four ligands about the origin such that each is aligned with its intended axis within the plane perpendicular to the axis of rotation is calculated. The octahedron is then rotated about this axis by the average of these four angles. This occurs iteratively until, for a given iteration, the sum (in quadrature) of the three rotation angles is less than the already-mentioned defined tolerance.

At this point, the octahedron is optimally aligned with the basis (given the limitation that there may be angular distortion) and the van Vleck modes can be calculated.

§.§ Ignoring or including angular distortion: a comparison

To evaluate the utility of calculating the van Vleck modes without disregarding the angular distortion, we perform a comparison between the two approaches. We have calculated the van Vleck distortion modes and associated parameters for octahedra in NaNiO$_2$ and LaMnO$_3$ with both a method that ignores angular distortion and calculates modes along bond directions (consistent with the $Q_2$ and $Q_3$ equations defined by kanamori1960crystal ), and a method that used Cartesian coordinates in order to take angular distortion into account. Table <ref> shows this for these two materials. Firstly, for the van Vleck modes calculated without ignoring angular distortion, we can see the octahedral shear modes ($Q_4$, $Q_5$, $Q_6$) are larger for the material with higher angular distortion (as quantified using bond angle variance). While the effect of ignoring angular distortion is significant for the $Q_4$, $Q_5$, and $Q_6$ modes, it makes negligible difference for the calculation of $Q_2$ and $Q_3$ modes, and the associated $\rho_0$ and $\phi$ parameters. It is therefore likely a reasonable approximation to take, particularly for calculation of $\phi$ as is common in literature, even for octahedra which exhibit higher angular distortion. However, there is a definite loss of information in assuming the shear modes $Q_4$ to $Q_6$ are zero. The impact of this is assessed in the case studies.

A comparison between calculated $\phi$, $\rho_0$, $Q_2$, $Q_3$, $Q_4$, $Q_5$, $Q_6$, and $\eta$ values for NaNiO$_2$ and LaMnO$_3$ at room-temperature, calculated using orthogonal axes as described in this report (“Cartesian" method), and alternatively by ignoring angular distortion and calculating van Vleck modes along bond lengths (“Kanamori" method[So named because the equations originate in kanamori1960crystal ]). The centre of the octahedron is taken as the central ion position. Values were calculated using crystal structures reported on the Inorganic Crystal Structure Database (ICSD). To demonstrate the difference in angular distortion, the Bond Angle Variance (defined in Equation <ref>) is also tabulated. BAV is rounded to the third significant figure; $Q$ modes and related parameters are rounded to the 4th decimal place.

2cNaNiO$_2$     2cLaMnO$_3$

ICSD code      2c415072     2c50334

Ref.      2csofin2005new      2crodriguez1998neutron

Octahedron      2cNiO$_6$     2cMnO$_6$

JT-active      2cYes     2cYes

Connectivity     2cedge      2ccorner

BAV ($^{\circ2}$)     2c35.2      2c0.45

Method     Kanamori      Cartesian     Kanamori     Cartesian

$Q_2$ (Å)     0.0000     0.0000     0.2745     0.2745

$Q_3$ (Å)     0.2834     0.2833     -0.0860     -0.0860

$Q_4$ (Å)     0     0.2078     0     0.0130

$Q_5$ (Å)     0     0.2001     0     0.0114

$Q_6$ (Å)     0     0.2001     0     0.0361

$\phi$ ($^\circ$)     0.0000[Note that $\phi=0^\circ$ is equivalent to 120$^\circ$ or 240$^\circ$.]     0.0000     107.3929     107.4034

$\rho_0$ (Å)     0.2834     0.2833     0.2877     0.2876

$\Delta_\mathrm{shear}$ (Å)     N/A     0.3534     N/A     0.0389

$\Delta_\mathrm{anti-shear}$ (Å)     N/A     0     N/A     0

$\eta$     N/A      1.0     N/A     1.0

§ CASE STUDIES

§.§ Temperature-dependence of octahedral shear in LaAlO$_3$

Perovskite and perovskite-like crystal structures are amongst the most important and widely-studied crystalline material classes in materials science today. Perovskite crystal structures have ABX$_3$ chemical formulae, with A and B being ions at the centres of dodecagons and octahedra, respectively, with the X anion constituting the vertices of these polyhedra. The BX$_6$ octahedra interact via corner-sharing interactions. There are also perovskite-like crystal structures such as the double perovskites, $A_2BB^\prime X_6$ <cit.>, for which many of the same principles apply.

The ideal perovskite system would be cubic, with space group $Pm\bar{3}m$, but many related structures with lower symmetry are known. This typically occurs in three situations <cit.>:

* when there is a mismatch between the ionic radii of the octahedrally-coordinated BX$_6$ cation and the dodecagonally-coordinated AX$_{12}$ cation, resulting in tilting of the octahedra; see Figure <ref>(e).

* when there is displacement of the central cation from the centre of the octahedron, typically due to the pseudo Jahn-Teller effect.

* when the ligands of the octahedron are distorted by electronic phenomena such as the first-order Jahn-Teller effect.

In this case study, we focus on the first case, where a size mismatch results in octahedral tilting. Octahedra are often modelled as rigid bodies, but in practice they are not rigid in all systems, and the octahedral tilting will often induce strain resulting in angular distortion. This is typically far smaller than that seen in edge-sharing materials such as NaNiO$_2$, but it is large enough that it cannot be disregarded when attempting to fully understand the structure of the material. As was noted by Darlington <cit.>, this angular distortion commonly manifests as shear.

The results of our analysis on LaAlO$_3$ as a function of temperature. (a) octahedral tilting angle as reported by hayward2005transformation  and extracted using DataThief III <cit.>. (b) radii of the minimum bounding ellipsoid fitted to the O anions of the AlO$_6$ octahedra using PIEFACE <cit.>. (c) octahedral shear parameter $Q_5$ of the AlO$_6$ octahedra, where $Q_5=-Q_4=-Q_6$, calculated using VanVleckCalculator, compared with $\sigma_\zeta$ the bond angle standard deviation (orange). (d) shear fraction $\eta$, defined in Equation <ref>. (e) the transition between low-symmetry (tilting) and high-symmetry (tilt-free) perovskite structures, adapted with permission from angel2005general, APS Physical Review Letters, 95, 025503, copyright 2005 American Physical Society. (f) the perovskite crystal structure of LaAlO$_3$ at 4.2 K from hayward2005transformation .

LaAlO$_3$ is a perovskite-like $ABX_3$ material which is cubic (space group $Pm\bar{3}m$) above around $\sim$830 K, but which exhibits a rhombohedral distortion below this temperature (with space group $R\bar{3}c$) due to octahedral tilting <cit.>, see see Figure <ref>(e) and (f). Throughout both temperature regimes, there is the absence of bond length distortion; a calculation of the bond length distortion index would yield a value of zero at all temperatures. In the low-temperature regime, the magnitude of the distortion continuously decreases with increasing temperature, reaching zero at the transition temperature. Most commonly in the literature, the tilting angle between the octahedral axis and the $c$-axis (0$^\circ$ in the cubic phase) is used to quantify this distortion; for LaAlO$_3$, this is shown in Figure <ref>(a). The strain induced by this distortion results in intra-octahedral angular distortion. hayward2005transformation  model this in terms of strain tensors, finding a linear temperature dependence below the transition temperature, which differs from the temperature-dependence of the tilting angle (which resembles an exponential decline), implying the two are related but distinct phenomena. cumby2017ellipsoidal  instead model the octahedral distortion for this same dataset using the radii of a minimum-bounding ellipsoid, and also find approximately linear temperature dependence of the long and short radii as they approach convergence (see Figure <ref>(b)).

Here, we calculate the van Vleck shear modes. Due to the symmetry of the octahedral tilting, there is only one independent shear mode, and $Q_5=-Q_4=-Q_6$. We compare this with the bond angle standard deviation given in Equation <ref>, see Figure <ref>. We see that despite being distinct parameters, the temperature dependence of both is entirely identical. We attribute this to the shear fraction, $\eta$, being precisely 1 for all temperatures where there is angular distortion, meaning that shear is completely correlated with angular distortion.

§.§ Big box analysis of Pair Distribution Function data on LaMnO$_3$

The Jahn-Teller distortion in LaMnO$_3$, a perovskite-like $ABX_3$ material which has the crystal structure shown in Figure <ref>(a), occurs as a consequence of degeneracy in the $e_g$ orbitals on the high-spin $d^4$ Mn$^{3+}$ ion. At ambient temperatures, it is a prime example of a cooperative Jahn-Teller distortion, exhibiting long-range orbital order where the elongation of the Jahn-Teller axis alternates between the $a$ and $b$ directions for neighbouring MnO$_6$ octahedra, never occurring along the $c$ direction <cit.> [Figure <ref>(b)]. With heating through $\sim$750 K, the Jahn-Teller distortion can no longer be observed in the average structure obtained from Bragg diffraction <cit.>. However, the Jahn-Teller distortion persists locally as has been shown by pair distribution function <cit.> and EXAFS <cit.> measurements. This transition is one of the most widely-studied orbital order-disorder transitions for the first-order Jahn-Teller distortion. The high-temperature orbital regime has been described theoretically in terms of a three-state Potts model <cit.>, a view supported by big box analysis of combined neutron and x-ray pair distribution function data <cit.>, as performed using RMCProfile <cit.>.

The perovskite-like structure of LaMnO$_3$, as obtained from ICSD structure 50334, is shown in (a). (b) the orbital ordering at room-temperature in LaMnO$_3$, and is reprinted with permission from khomskii2020orbital, Chem. Rev. (2021), 121, 5, 2992–3030, copyright 2021 American Chemical Society. (c) and (d) show polar plots with a point representing the calculated $\phi$ and $\rho_0$ values for each MnO$_6$ octahedron in a $10\times10\times8$ supercell of LaMnO$_3$ at room-temperature, as obtained from reverse Monte Carlo analysis of neutron Pair Distribution Function data in thygesen2017local . In (c), orthogonal axes were used (i.e. angular distortion was included in the calculation, using the method described in this manuscript), whereas in (d) the Mn-O bond directions were taken as the axes regardless of orthogonality. (e) the Mn$^{3+}$ octahedra which exhibit a mixed $Q_2$-$Q_3$ type distortion due to the first-order Jahn-Teller effect, manifesting as three different bond lengths, labelled in ascending order of length as $s$ (orange), $m$ (grey), and $l$ (green). (f) a histogram of the smallest to largest Mn-O bond length within each octahedron in the $10\times10\times8$ supercell, with the blue vertical lines indicating the bond lengths in the average structure.

In this case study, we take a $10\times10\times8$ supercell of LaMnO$_3$, obtained using RMCProfile against total scattering data obtained at room-temperature, and previously published in the aforementioned work <cit.>. Results are shown in Figure <ref>. We repeat the analysis of this supercell from the perspective of the $E_g(Q_2,Q_3)$ van Vleck distortion modes, using two different approaches: (1) the algorithm for automatically determining a set of orthogonal axes is applied to each octahedron individually

, and (2) following the van Vleck equations <ref> and <ref> proposed by kanamori1960crystal  where angular distortion is disregarded. In each of these cases the crystallographic site of the supercell is taken as the origin, and so thermally-driven variations in the Mn position will not affect the result.

As can be seen in Figures <ref>(d)-(f), there are two clusters of octahedra within the polar plot, occurring at $\phi\approx\pm107^\circ$. This corresponds to occupation of the $d_{y^2}$ orbitals ($+$) and of the $d_{x^2}$ orbitals ($-$). In both cases, the superposition of perpendicular $Q_3$ compression and elongation modes results in an octahedron with mixed $Q_2$-$Q_3$ character. This finding is consistent with previous works which placed MnO$_6$ octahedra from LaMnO$_3$ into the framework of an $E_g(Q_2,Q_3)$ polar plot <cit.>.

Figure <ref>(e) shows the MnO$_6$ octahedron in the average structure of LaMnO$_3$ at room temperature, with the three different bond lengths plotted in Figure <ref>(f) along with a histogram of all the bond lengths in the supercell. This shows how the combination of the $Q_2$ and $Q_3$ distortion modes manifests in the octahedral distortion.

The $Q_2$ contribution to the distortion, as seen from the three different Mn-O bond lengths in LaMnO$_3$, is also present in Jahn-Teller-distorted $A$CuF$_3$ ($A$=Na,K,Rb) <cit.> and even in some Jahn-Teller-undistorted perovskites <cit.>, indicating it is related to the structure. It is not intrinsic to Jahn-Teller-distorted manganates, as it is absent in high-spin $d^4$ Mn$^{3+}$ with edge-sharing octahedral interactions and colinear orbital ordering such as $\alpha$-NaMnO$_2$ and LiMnO$_2$ (checked using ICSD references 15769 and 82993 <cit.> respectively). The $Q_2$ component to the octahedral distortion is therefore likely intrinsic to the crystal structure <cit.>, which occurs as a result of octahedral tilting reducing the symmetry from cubic $Pm\bar{3}m$ to $Pnma$. In LaMnO$_3$, the combination of the $Q_2$ component to the distortion and the orbital ordering [Figure <ref>(b)] are a possible distortion of the $Pnma$ space group. In this way, the orbital ordering may be coupled to the octahedral tilting, a link previously made by lufaso2004jahn .

Finally, we also calculate the $Q_4$ to $Q_6$ octahedral shear modes for all octahedra in the supercell, presented as a histogram in Figure S2
in Supplementary Information. We present the average and standard deviation, as calculated assuming orthogonal axes and with the automated octahedral rotation: $Q_4=-0.02\pm0.13$ Å, $Q_5=0.02\pm0.10$ Å, and $Q_6=-0.00\pm0.11$ Å. In each case, the magnitude of the distortion is zero within standard deviation, and also contains the value from the average structure presented in Table <ref> within the range of error. This low level of shear generally supports the validity of calculating the $E_g(Q_2,Q_3)$ van Vleck modes along bond directions rather than a Cartesian coordinate system for a system like LaMnO$_3$.
It is interesting to note that the standard deviation is higher for $Q_4$, which quantifies the shear within the plane in which there is orbital ordering.

§.§ Effect of pressure on the JT distortion in NaNiO$_2$

In recent years, there have been several studies looking at the effect of applied pressure on the Jahn-Teller distortion in crystalline materials <cit.>. Most of these have shown that, as a general rule, pressure reduces the magnitude of the Jahn-Teller distortion as a consequence of the elongated bond being more compressible than the shorter bonds.

zhou2011jahn  use van Vleck modes to quantify the effect of pressure on the Jahn-Teller distortion in the corner-sharing perovskite-like compounds LaMnO$_3$ and KCuF$_3$. While application of pressure reduces the magnitude of the distortion, as quantified using $\rho_0$ (Equation <ref>), they argue that it does not change the orbital mixing $\phi$ (Equation <ref>)
. KCuF$_3$ has similar orbital ordering to LaMnO$_3$, except the degeneracy is due to the $d^9$ hole rather than an electron. The variable-pressure crystal structures for KCuF$_3$ are available on ICSD (catalog codes 182849-182857), and are utilised here.

The crystal structures of Jahn-Teller-active $C2/m$
NaNiO$_2$ and inactive $R\bar{3}c$ Fe$_2$O$_3$ are shown in (a) and (b) respectively. (c) and (d) show a comparison of various metrics for quantifying the degree of Jahn-Teller distortion as a function of pressure, for NiO$_6$ octahedra in NaNiO$_2$ and FeO$_6$ octahedra in Fe$_2$O$_3$ respectively. The parameters subject to comparison are the magnitude $\rho_0$, bond length distortion index, effective coordination number, and quadratic elongation. Dashed lines indicate a linear fit to the data, whereas solid lines connect data points.

An $E_g(Q_2,Q_3)$ radial plot comparing the pressure-dependence of the $M$O$_6$ ($M$=Ni,Cu) octahedra for KCuF$_3$ and NaNiO$_2$ between 0 and 5 GPa, where $\rho_0$ is normalised to the value at the lowest measured pressure and the dashed lines represent the average $\phi$ for each material within this pressure range.

The pressure-dependence of the shear and angular distortion in (a) Jahn-Teller-distorted NiO$_6$ octahedra in NaNiO$_2$ and (b) Jahn-Teller-undistorted FeO$_6$ octahedra in Fe$_2$O$_3$. Shear distortion is represented with the $Q_4$, $Q_5$, and $Q_6$ modes for the octahedra, and angular distortion is represented by bond angle variance. Dashed lines indicate a fitted straight line to the data, whereas solid lines are plotted from point to point. $\eta$ is the angular shear fraction defined in Equation <ref>. Note that for the NiO$_6$ octahedra, $Q_5=Q_6$, whereas for FeO$_6$ octahedra, $Q_4=Q_6=-Q_5$. For Fe$_2$O$_3$, the average position of the O ligands were taken as the centre of the octahedron.

We previously studied the effect of pressure on the Jahn-Teller distortion in NaNiO$_2$ <cit.>, by performing Rietveld refinement <cit.> of neutron diffraction data from the PEARL instrument <cit.> at the ISIS Neutron and Muon Source. However, we did not utilise the van Vleck distortion modes, instead quantifying the Jahn-Teller distortion using the bond length distortion index <cit.> and the effective coordination number <cit.>. In that study, we found no deviation from the ambient-pressure space group $C2/m$ <cit.>, shown in Figure <ref>(a), for all pressure points at room-temperature up to $\sim$4.5 GPa. This space group permits only four short{long} and two long{short} bonds or 6 equal bond lengths, depending on the angle $\beta$, and so throughout the measured pressure range there exists no $Q_2$ character to the Jahn-Teller distortion, consistent with the principle that hydrostatic pressure does not change orbital mixing <cit.>.

Here, we perform a fresh analysis of the variable-pressure octahedral behaviour as a function of pressure at room temperature in NaNiO$_2$ in terms of the $E_g(Q_2,Q_3)$ van Vleck distortion modes. For a reference we sought a material which does not exhibit a first-order Jahn-Teller distortion but does exhibit bond length distortion; for this purpose, we selected Fe$_2$O$_3$, the pressure dependence of which was previously studied in finger1980crystal , and which exhibits bond length distortion due to its face- and edge-sharing octahedral connectivity. Fe$_2$O$_3$ contains high-spin $d^5$ Fe$^{3+}$ cations within octahedra which interact via both face- and edge-sharing interactions. It should be noted that Fe$_2$O$_3$ likely exhibits some very subtle pseudo Jahn-Teller distortion (related to, but distinct from the first Jahn-Teller effect discussed here) on account of the Fe$^{3+}$ ions <cit.>, but this does not impact the discussion in any meaningful way.

In Figure <ref>(c) we compare (for NaNiO$_2$) $\rho_0$ with three other parameters (bond length distortion index, quadratic elongation, and effective coordination number) which are often used to parametrise the magnitude of the Jahn-Teller distortion. The trend for each is near identical, although the magnitudes differ greatly, indicating that each is a reasonable parameter for quantifying the magnitude of the Jahn-Teller distortion. This can be compared to Figure <ref>(e) which shows the same parameters for the Jahn-Teller-undistorted FeO$_6$ octahedra in Fe$_2$O$_3$, where it can be seen that $\rho_0$ remains approximately at zero throughout the measured pressure range, despite a high level of bond length distortion as represented by the bond length distortion index, effective coordination number, and quadratic elongation (a similar plot for KCuF$_3$ can be seen in SI, Figure S3
). This means that, while these parameters are valid for quantifying the magnitude of Jahn-Teller distortion, they are also sensitive to other kinds of distortion. $\rho_0$ is calculated using $Q_2$ and $Q_3$ which have $E_g$ symmetry, and so $\rho_0$ will only be non-zero for a distortion with $E_g$ symmetry. Thus, it is arguably the ideal choice for parameterising the magnitude of this type of Jahn-Teller distortion. However, while $\rho_0$ is more reliable than the other parameters shown in Figures <ref>(c,d) for demonstrating the presence of a Jahn-Teller distortion, it is not always strictly zero for a Jahn-Teller-inactive octahedron, as it will have a non-zero value if the octahedron is distorted with an $e_g$ symmetry. For example, the NaO$_6$ octahedron in $C2/m$ NaNiO$_2$ has the same symmetry as the NiO$_6$ octahedron, and so exhibits a value of $\rho_0$ between 0.065 and 0.05 within the studied pressure range [Figure S4
in Supplementary Information], and Jahn-Teller-inactive FeO$_6$ octahedra in $R$FeO$_3$ perovskites have non-zero $\rho_0$ due to the $E_g$ symmetry of the distorted octahedra, as shown in zhou2008intrinsic .

Figure <ref> shows a polar plot for the behaviour of NaNiO$_2$ and KCuF$_3$ in the range 0 to 5 GPa (the measured range for NaNiO$_2$). It can be seen that within this pressure range, the magnitude of the Jahn-Teller distortion decreases far more for KCuF$_3$ than NaNiO$_2$; this reflects the fact that KCuF$_3$ is more compressible, with a bulk modulus 57(1) GPa <cit.> compared with 121(2) GPa for NaNiO$_2$ <cit.>, as obtained by a fit to the third-order Birch-Murnaghan equation-of-state <cit.>. Within this pressure range we see that $\phi$ does not change with pressure for either material, and that this property is true regardless of whether $\phi$ is or is not a special angle (as in Table <ref>), consistent with the interpretation of zhou2011jahn .

Finally, in the previous study <cit.>, we showed using specific O-Ni-O bond angles that pressure reduces the angular distortion for NaNiO$_2$.
Here, we show that pressure also reduces the related shear distortion in NaNiO$_2$. This is demonstrated in Figure <ref> where we plot the octahedral shear $Q_4$, $Q_5$, and $Q_6$ modes
for NaNiO$_2$ and Fe$_2$O$_3$ against the bond angle standard deviation, $\sigma_\zeta$, defined in Equation <ref>. Unlike the AlO$_6$ octahedra in LaAlO$_3$ [Figure <ref>], for NiO$_6$ octahedra in NaNiO$_2$ there is no perfect correlation between the shear modes and angular distortion despite $\eta\approx1$, because there is more than one independent shear mode, but we can see that shear distortion and angular distortion are still highly correlated. However, for Fe$_2$O$_3$ the shear fraction $\eta << 1$ and there is no correlation between the shear distortion modes and angular distortion. This difference in behaviour likely arises because the main driver of the change is a continuous decrease in the Jahn-Teller distortion in NaNiO$_2$, as compared to Fe$_2$O$_3$ where positions of the oxygen anions are determined by the reduced degrees of freedom arising from trying to satisfy multiple face- and edge-sharing interactions. This result could only be achieved by calculating the van Vleck modes in a Cartesian coordinate system as outlined in this paper, as opposed to calculating the distortion modes along bond directions, indicating the relevance of calculating the van Vleck modes in this way, and of the shear fraction $\eta$ we propose in this work.

§ CONCLUSION

We present VanVleckCalculator, a code package written in Python 3 for the calculation of octahedral van Vleck distortion modes. These modes are particularly important for understanding the behaviour of the Jahn-Teller distortion, and we have shown that the parameter $\rho_0$ (which is based on the van Vleck $Q_2$ and $Q_3$ modes) is a more reliable way of quantifying the Jahn-Teller distortion than other oft-used parameters such as the bond length distortion index.

We show the importance of using a Cartesian set of coordinates for this calculation, instead of calculating the modes along bond directions, as is often done in the literature. This is because calculating the van Vleck distortion modes along bond directions relies on the assumption that there is no angular distortion or octahedral shear, which is often a false assumption and artificially constrains the $Q_4$, $Q_5$, and $Q_6$ modes to be zero. We show that there is value in calculating these later modes, for instance in understanding the effect of octahedral tiling on octahedra in perovskite-like materials. These shear modes will also be useful for parameterising the Jahn-Teller effect when the degeneracy occurs in the $t_{2g}$ orbitals and results in a trigonal distortion, because their symmetry matches the distortion.

We also show that octahedral shear correlates with angular distortion for materials under the influence of tuning parameters such as pressure or temperature where there is a continuously-varying distortion, such as octahedral tilting (as in LaAlO$_3$) or first-order Jahn-Teller distortion (as in NaNiO$_2$). However, there is no correlation when the distortion is due to competing interactions due to face- or edge-sharing octahedra (as in Fe$_2$O$_3$). We propose a new parameter, the shear fraction $\eta$ (defined in Equation <ref>), which can be used to predict whether there will be correlation between octahedral shear modes and angular distortion.

The authors thank Andrew L. Goodwin (University of Oxford) for useful discussions and for sharing the LaMnO$_3$ supercell from thygesen2017local .
The authors thank James Cumby (University of Edinburgh) for sharing the crystal structures for LaAlO$_3$ originally published in hayward2005transformation  and subsequently analysed in cumby2017ellipsoidal .
The authors acknowledge comments on earlier drafts of this manuscript from James M. A. Steele, Fiamma Berardi, and Venkateswarlu Daramalla, all at the University of Cambridge.
LNC acknowledges Annalena R. Genreith-Schriever (University of Cambridge) and Ben Tragheim (University of Warwick) for useful discussions.

LNC acknowledges a scholarship EP/R513180/1 to pursue doctoral research from the UK Engineering and Physical Sciences Research Council (EPSRC).

Graphical software

Graphs and radial plots were prepared using Matplotlib <cit.>. Crystal structures figures were made using Vesta-III <cit.>.