Van Vleck analysis of angularly distorted octahedra using VanVleckCalculator

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 first-order Jahn–Teller distortions, but determination of the Van Vleck modes of an octahedron is complicated by the presence of angular distortion of the octahedron. This problem is most commonly resolved by calculating the bond distortion modes (Q₂, Q₃) along the bond axes of the octahedron, disregarding the angular distortion and losing information on the octahedral shear modes (Q₄, Q₅ and Q₆) 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 package for doing this, VanVleckCalculator, is introduced and some examples of its use are given. Finally, it is shown that octahedral shear and angular distortion are often, but not always, correlated, and a parameter η is proposed as the shear fraction. It is demonstrated that η can be used to predict whether the values will be correlated when varying a tuning parameter such as temperature or pressure.

1. Introduction

The Van Vleck distortion modes (Van Vleck, 1939) describe all possible displacements of octahedrally coordinated ligands about a core atom. They are particularly useful in the context of the Jahn–Teller (JT) effect (Jahn & Teller, 1937), which in general occurs when a high-symmetry coordination is de­stabilized with respect to a deviation to lower symmetry as a consequence of electronic degeneracy. The JT effect distorts the crystal structure via the JT distortion. While the JT distortion is not unique to octahedra in bulk crystalline materials, it was in octahedra that it was first observed experimentally (Bleaney & Bowers, 1952), and it was in materials with JT-distorted octahedra that colossal magneto­resistance (Millis et al., 1996) and high-temperature superconductivity (Fil et al., 1992; Keller et al., 2008) were discovered.

A transition metal cation in an octahedral configuration will have its d orbitals split into three t_2g orbitals at lower energy and two e_g orbitals at higher energy. [Here 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.] It will have a number n of electrons in these d orbitals (hereafter described as dⁿ). For certain values of n and, where applicable, certain low- or high-spin characters,¹ there will exist multiple orbitals that could be occupied by an electron or an electron hole with equal energy. This degeneracy is destabilizing, 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 minimized energy. This is shown for a low-spin d⁷ transition metal cation (such as Ni³⁺ or Co²⁺) in Fig. 1, though such distortions may occur for any value of n in dⁿ where there is a degenerate occupancy. The stabilization energy due to the JT effect is larger for e_g degeneracy than t_2g degeneracy, and so the effect will remain prominent up to higher temperatures, and is hence more widely studied, in JT-active materials with e_g degeneracy (Castillo-Martínez et al., 2011).

Figure 1 The orbital rearrangement due to a tetragonal elongation for an octahedrally coordinated low-spin d^7 transition metal ion, which typically occurs as a result of the first-order JT effect.

In the literature, various techniques for parameterizing the JT distortion are used. An often-used example (Kimber, 2012; Lawler et al., 2021; Nagle-Cocco et al., 2022; Genreith-Schriever et al., 2023) is the bond-length distortion index, defined by Baur (1974) as

where n is the number of ligands (6 for an octahedron), l_i is the distance between the core ion and the ith coordinated ion, and l_av is the average of all the distances between the core ion and coordinated ions.

A similar parameter (Shirako et al., 2012; Sarkar et al., 2018; Nagle-Cocco et al., 2022) is the effective coordination number ECoN, which for an octahedron deviates from 6 only when there is bond-length distortion. ECoN is defined by Hoppe (1979) as

where is a modified average distance defined as

Here, l_min is the smallest core–ligand bond length in the octahedron.

Finally, a third parameter used to quantify the JT distortion (Schofield et al., 1997; Kyono et al., 2015; Mikheykin et al., 2015) is the quadratic elongation 〈λ〉, defined by Robinson et al. (1971) as

where l₀ 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 (Cumby & Attfield, 2017), involving fitting an ellipsoid to the positions of the ligands around a coordination polyhedron, calculating the three principal axes of the ellipsoid R₁, R₂ and R₃, where R₁ ≤ R₂ ≤ R₃, and using the variance of these three radii as a metric for the distortion. This has been applied to the first-order JT distortion by Pughe et al. (2023).

These parameterizations each have merits. However, they are not sensitive to the symmetry of the octahedral distortion. The Van Vleck modes are conceptually different from each of these for quantifying the JT distortion because they can be used to quantify distortion with the precise symmetry of the transition metal e_g orbitals. This is important because JT 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 (Child & Roach, 1965; Bacci et al., 1975; Holland et al., 2002; Halcrow, 2009; Teyssier et al., 2016; Schmitt et al., 2020; Streltsov et al., 2022), although there is relatively little unambiguous experimental evidence for a JT-induced shear compared with more typical E_g distortion.

In this paper, we present a Python (Van Rossum & Drake, 2009) 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 η, 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.

2. Theory

Within an octahedron, we can split the six 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 multiple of 90°), there would exist a basis where each of the three axes exist directly along the x, y and z axes, 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 Fig. 2, where each pair of ions is represented by a different colour.

Figure 2 The six Van Vleck modes exhibited for an octahedron, with sites labelled using the notation in the Theory section. For the octahedra exhibiting Q1, Q2 and Q3 distortions, there is no angular distortion; for the octahedra exhibiting Q4, Q5 and Q6 distortions, there is no bond-length distortion. For the octahedral shear (Q4, Q5 and Q6) 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.

For each of the six ligands, we define a set of coordinates , and , where α is a, b or c denoting the pair in which the ligand is, and β is − or + denoting the 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 (capitalized to distinguish from true coordinates) which give the displacement of the ion within an axis away from its ideal position. For instance, for the ion with α = a and β = −, = , = and = . See Fig. 2 for clarification of the ion notation.

Using these coordinates, the first six Van Vleck modes (Q_j, j = 1–6) are defined as follows (Van Vleck, 1939):

We only discuss these first six Van Vleck modes, which are shown in Fig. 2. Q₁ to Q₃ describe bond-length distortions, whereas Q₄ to Q₆ describe octahedral shear distortions. Q₁ is a simple expansion/contraction mode which does not affect symmetry and will not be discussed further.

Q₂ and Q₃ are a planar rhombic distortion and a tetragonal distortion, respectively; they are considered degenerate due to the Hamiltonian, which is discussed for instance by Kanamori (1960). These two modes form a basis for distortions describing different octahedral configurations with the symmetry of the transition metal e_g orbitals (Goodenough, 1998; Khomskii & Streltsov, 2021). These modes are of most relevance for first-order JT distortions occurring due to degenerate e_g orbitals. A phase space of possible octahedral configurations can be constructed using these two parameters (Kanamori, 1960), as shown in Fig. 3. Here the magnitude of the distortion ρ₀ can be calculated as

and the angle ϕ of this distortion, being of purely Q₃ character, can be calculated by

(Note that this angle does not represent a physical angle within the octahedron.)

Figure 3 The Q2–Q3 phase space for elongated octahedra, with a representation of the values ρ0 and ϕ. Based on a figure from an article by Goodwin (2017).

All possible combinations of the Q₂ and Q₃ modes correspond to a particular angle ϕ and hence a particular configuration as shown in Fig. 3. The structural effect of a rotation of ϕ within a range of 120° can be quite significant, as shown in Fig. 3; such changes can manifest as a JT-elongated (compressed) octahedron with four short (long) and two long (short) bonds [such as NiO₆ in NaNiO₂ (Nagle-Cocco et al., 2022)] or two short, two medium and two long bonds [such as LaMnO₃ (Rodriguez-Carvajal et al., 1998)].

A characteristic of JT 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 JT distortion will have some planar rhombic (Q₂) or tetragonal (Q₃) 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 (Child & Roach, 1965; Bacci et al., 1975; Holland et al., 2002; Halcrow, 2009; Teyssier et al., 2016; Schmitt et al., 2020; Streltsov et al., 2022). For the more commonly studied case of a degeneracy in the e_g orbitals, the effect of a rotation of ϕ similarly changes the symmetry of the d orbitals. Fig. 1 shows the splitting of the d orbitals in an octahedrally coordinated d⁷ transition metal due to an elongation-type first-order JT distortion, where the tetragonal elongation occurs along the z axis. Note that the unpaired e_g electron occupies the d_z² orbital. In the opposite case of a compression-type first-order JT distortion along the z axis, the lower-energy, and hence singly occupied, orbital would be the ; this would correspond to a rotation in ϕ of 180°. More generally, as a function of ϕ, there exist a set of special angles separated by a 60° rotation corresponding to a particular e_g orbital being singly occupied by a d electron. These are tabulated in Table 1. An octahedron for which ϕ does not correspond to one of these special angles exhibits orbital mixing (Rodriguez-Carvajal et al., 1998; Zhou & Goodenough, 2008b).

Table 1
The special angles in the Q₂–Q₃ phase space (Fig. 3) as a function of = arctan(Q₂/Q₃), with the associated singly occupied e_g orbital, for d⁴ and low-spin d⁷ 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.

ϕ (°)	Ψ(ϕ)
0	dz^2
60
120	dy^2
180
240	dx^2
300

The Q₄ to Q₆ modes describe shear of the octahedra, i.e. the effect whereby paired ligands on opposite sides of a central ion are displaced in opposite directions and have trigonal T_2g character. The shear modes may be used to quantify the JT distortion in octahedra where the degeneracy occurs within t_2g orbitals (Child & Roach, 1965; Teyssier et al., 2016). The magnitude of the calculated shear is typically correlated with angular distortion, which is commonly quantified using the metric called the bond-angle variance (BAV) (Robinson et al., 1971), defined here as

m is the number of bond angles (i.e. 12 for an octa­hedron), ζ_i is the ith bond angle and ζ₀ is the ideal bond angle for a regular polyhedron (i.e. 90° for an octahedron). However, for direct comparison with the shear modes, it is more appropriate to use the standard deviation σ_ζ.

For an octahedron with non-zero T_2g(Q₄, Q₅, Q₆) 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 η, demonstrated in Fig. 4 and defined below.

Figure 4 Three possible octahedral shear/anti-shear distortions, with the associated value of the shear fraction η as defined in equation (22). In the case where η = 0, the only distortion is anti-shear within a single plane. In the case where η = 0.5, there are two planes in which there is distortion, a shear and an anti-shear distortion equal in magnitude. In the case where η = 1, there is a plane with a purely shear distortion.

First, we must define a set of shear and `anti-shear' angular indices, which are modifications of equations (8) to (10) in terms of angles rather than displacements. The indices are represented with Δ and a subscript corresponding to the plane in which rotation occurs: the ab plane corresponds to the Q₄ mode, the ac plane to the Q₅ mode and the bc plane to the Q₆ mode. The absence or presence of a prime symbol ′ designates whether the index represents shear or anti-shear, respectively. Finally, the δ 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₄ to Q₆) would occur. These are defined thus (see the supporting information, Fig. S7):

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

From here, we define the shear fraction η as

This η parameter will be important in interpreting the relation between the angular distortion σ_ζ and the Van Vleck shear modes Q₄ to Q₆.

3. Implementation

In this section, the algorithm used to calculate the Van Vleck distortion modes is discussed. It is written using Python 3 (Van Rossum & Drake, 2009) as a package called VanVleckCalculator, with the full code available on GitHub (Nagle-Cocco, 2023) and also presented with annotations in the supporting information. Data handling and some calculations make use of NumPy (Harris et al., 2020) and crystal structures are handled using PyMatGen (Ong et al., 2013).

A flow chart showing the octahedral rotation algorithm can be found in Fig. 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 the supporting information.

3.1. Selecting an origin

Selection of the origin is a key step in calculating Van Vleck modes. The most common approach, for an MX₆ 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 (Rietveld, 1969) of Bragg diffraction data, where the M ion is likely to occur on a high-symmetry Wyckoff site. A third, similar, option would be to choose the average position of the six ligands as the origin in space. An example of when this may be a desirable choice would be for systems exhibiting a pseudo-JT effect (also called the second-order JT 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 `big-box' pair distribution function (PDF) analysis (Tucker et al., 2007) and molecular dynamics (MD) (Bocharov et al., 2020) 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 JT distortion. In a crystallographic unit cell, thermal motion of atoms is typically represented by variable atomic displacement parameters (ADPs) (Peterse & Palm, 1966). In contrast, a supercell should reflect local phenomena, for instance exhibiting local JT distortions in a system with a non-cooperative JT 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 therefore not be appropriate to set the core ion as the centre of the octahedron in a supercell, 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 as a result of random motion. The alternative option would simply be to 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 six ligands or a specified set of coordinates.

3.2. 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 that are not integer multiples of 90°). 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; see, for example, previous work on the perovskite LaMnO₃ (Goodenough et al., 1961; Rodriguez-Carvajal et al., 1998; Capone et al., 2000; Chatterji et al., 2003; Zhou & Goodenough, 2008b; Zhou et al., 2011; Snamina & Oleś, 2016; Fedorova et al., 2018; Lindner et al., 2022), other perovskites (Alonso et al., 2000; Wang et al., 2002a; Tachibana et al., 2007; Zhou & Goodenough, 2008a; Castillo-Martínez et al., 2011; Franchini et al., 2011; Chiang et al., 2011; Dong et al., 2012; Fedorova et al., 2015; Ji et al., 2019; Xu et al., 2020; Ren et al., 2021) or non-perovskite materials (Moron et al., 1993; Cussen et al., 2001; Wang et al., 2002b) (we note that some reports 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₂ and Q₃ are defined according to the following equations which were first expressed by Kanamori (1960), 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 so that they are mathematically equivalent to equations (6) and (7)]:

This relies on the implicit assumption that the bond lengths are orthogonal. This is clearly a reasonable approximation in many cases, particularly when the angular distortion is very small. For instance, in LaMnO₃, the corner-sharing octahedral connectivity enables mismatched polyhedra to tessellate via octahedral tilting [Fig. 5(e) in Section 4.1] 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.

3.3. 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₄, Q₅ and Q₆) to be zero.

First, three orthogonal axes are taken as the x, y and z axes.² 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° about the x axis, the user would be recommended to give as axes [1, 0, 0], and . 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 assigned to one of these three axes. This assignment is performed according to the magnitude of the projection of a vector between the two atoms in a pair along a given axis, with the z axis assigned first, then the y axis from amongst the two pairs not assigned to the z axis, and finally 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 the first along the assigned vector and the ligand with the positive distance occurs second.

Second, the octahedron is rotated repeatedly about the x, y and z directions of the basis 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 × 10⁻⁴ rad in Van­VleckCalculator); if the total rotation exceeds the tolerance, the step is repeated.³ 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 minimize 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.

3.4. 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₂ and LaMnO₃ with both a method that ignores angular distortion and calculates modes along bond directions [consistent with the Q₂ and Q₃ equations defined by Kanamori (1960)] and a method that uses Cartesian coordinates in order to take angular distortion into account. Table 2 shows this for these two materials. Firstly, for the Van Vleck modes calculated without ignoring angular distortion, we can see that the octahedral shear modes (Q₄, Q₅, Q₆) 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₄, Q₅ and Q₆ modes, it makes negligible difference for the calculation of Q₂ and Q₃ modes and the associated ρ₀ and ϕ parameters. It is therefore likely to be a reasonable approximation to take, particularly for calculation of ϕ as is common in the literature, even for octahedra which exhibit higher angular distortion. However, there is a definite loss of information in assuming that the shear modes Q₄ to Q₆ are zero. The impact of this is assessed in the case studies.

Table 2
A comparison between calculated ϕ, ρ₀, Q₂, Q₃, Q₄, Q₅, Q₆ and η values for NaNiO₂ and LaMnO₃ 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 were originated by Kanamori (1960)]

The centre of the octahedron is taken as the central ion position. Values were calculated using crystal structures reported in the ICSD. To demonstrate the difference in angular distortion, the bond-angle variance BAV [defined in equation (13)] is also tabulated. The BAV is rounded to the third significant figure; Q modes and related parameters are rounded to the fourth decimal place.

	NaNiO2		LaMnO3
	Kanamori	Cartesian	Kanamori	Cartesian
ICSD code	415072		50334
Reference	Sofin & Jansen (2005)		Rodriguez-Carvajal et al. (1998)
Octahedron	NiO6		MnO6
JT-active	Yes		Yes
Connectivity	Edge		Corner
BAV (°^2)	35.2		0.45
Q2 (Å)	0.0000	0.0000	0.2745	0.2745
Q3 (Å)	0.2834	0.2833	−0.0860	−0.0860
Q4 (Å)	0	0.2078	0	0.0130
Q5 (Å)	0	0.2001	0	0.0114
Q6 (Å)	0	0.2001	0	0.0361
ϕ (°)	0.0000†	0.0000	107.3929	107.4034
ρ0 (Å)	0.2834	0.2833	0.2877	0.2876
Δshear (Å)	N/A	0.3534	N/A	0.0389
Δanti-shear (Å)	N/A	0	N/A	0
η	N/A	1.0	N/A	1.0

†Note that ϕ = 0° is equivalent to 120 or 240°.

4. Case studies

4.1. Temperature dependence of octahedral shear in LaAlO₃

Perovskite and perovskite-like crystal structures are amongst the most important and most widely studied crystalline material classes in materials science today. Perovskite crystal structures have ABX₃ chemical formulae, with A and B being ions at the centres of dodecagons and octahedra, respectively, and the X anions constituting the vertices of these polyhedra. The BX₆ octahedra interact via corner-sharing interactions. There are also perovskite-like crystal structures such as the double perovskites, A₂BB′X₆ (King & Woodward, 2010; Koskelo et al., 2023), for which many of the same principles apply.

The ideal perovskite system would be cubic, with space group , but many related structures with lower symmetry are known. This typically occurs in three situations (Woodward, 1997):

(i) when there is a mismatch between the ionic radii of the octahedrally coordinated B cation and the dodecagonally coordinated A cation, resulting in tilting of the octahedra [see Fig. 5(e)];

Figure 5 The results of our analysis of LaAlO3 as a function of temperature. (a) The octahedral tilting angle as reported by Hayward et al. (2005) and extracted using DataThief III (Tummers, 2006). (b) The radii of the minimum bounding ellipsoid fitted to the O anions of the AlO6 octahedra using PIEFACE (Cumby & Attfield, 2017). (c) The octahedral shear parameter Q5 of the AlO6 octahedra, where Q5 = −Q4 = −Q6, calculated using VanVleckCalculator, compared with the bond-angle standard deviation σζ (orange). (d) The shear fraction η, defined in equation (22). (e) The transition between low-symmetry (tilting) and high-symmetry (tilt-free) perovskite structures, adapted with permission from Angel et al. (2005), copyright (2005) American Physical Society. (f) The perovskite crystal structure of LaAlO3 at 4.2 K, drawn using the structure reported by Hayward et al. (2005).

(ii) when there is displacement of the central cation from the centre of the octahedron, typically due to the pseudo-JT effect;

(iii) when the ligands of the octahedron are distorted by electronic phenomena such as the first-order JT 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₂, but it is large enough that it cannot be disregarded when attempting a full understanding of the structure of the material. As was noted by Darlington (1996), this angular distortion commonly manifests as shear.

LaAlO₃ is a perovskite-like ABX₃ material which is cubic (space group ) above around ∼830 K but exhibits a rhombohedral distortion below this temperature (with space group ) due to octahedral tilting (Hayward et al., 2005) [Figs. 5(e) and 5(f)]. There is no bond-length distortion in either temperature regime; 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° in the cubic phase) is used to quantify this distortion; for LaAlO₃ this is shown in Fig. 5(a). The strain induced by this distortion results in intra-octahedral angular distortion. Hayward et al. (2005) modelled this in terms of strain tensors, finding a linear temperature dependence below the transition temperature. This differs from the temperature dependence of the tilting angle (which resembles an exponential decline), implying the two are related but distinct phenomena. Cumby & Attfield (2017) instead modelled the octahedral distortion for this same data set using the radii of a minimum-bounding ellipsoid and also found an approximately linear temperature dependence of the long and short radii as they approach convergence [see Fig. 5(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₅ = −Q₄ = −Q₆. We compare this with the bond-angle standard deviation given in equation (13) [Fig. 5(c)]. We see that, despite being distinct parameters, they have identical temperature dependences. We attribute this to the shear fraction η being precisely 1 for all temperatures where there is angular distortion, meaning that shear is completely correlated with angular distortion.

4.2. Big-box analysis of PDF data on LaMnO₃

The JT distortion in LaMnO₃, a perovskite-like ABX₃ material which has the crystal structure shown in Fig. 6(a), occurs as a consequence of degeneracy in the e_g orbitals on the high-spin d⁴ Mn³⁺ ion. At ambient temperatures it is a prime example of a cooperative JT distortion, exhibiting long-range orbital order where the elongation of the JT axis alternates between the a and b directions for neighbouring MnO₆ octahedra, never occurring along the c direction (Khomskii & Streltsov, 2021) [Fig. 6(b)]. With heating to ∼750 K, the JT distortion can no longer be observed in the average structure obtained from Bragg diffraction (Rodriguez-Carvajal et al., 1998). However, the JT distortion persists locally, as has been shown by PDF (Qiu et al., 2005) and EXAFS (García et al., 2005; Souza et al., 2005) measurements. This transition is one of the most widely studied orbital order–disorder transitions for first-order JT distortion. The high-temperature orbital regime has been described theoretically in terms of a three-state Potts model (Ahmed & Gehring, 2006, 2009), a view supported by big-box analysis of combined neutron and X-ray PDF data (Thygesen et al., 2017), as performed using RMCProfile (Tucker et al., 2007).

Figure 6 (a) The perovskite-like structure of LaMnO3 as obtained from ICSD structure 50334. (b) The orbital ordering at room temperature in LaMnO3. Reprinted with permission from Khomskii & Streltsov (2021), copyright (2021) American Chemical Society. (c) and (d) Polar plots with each point representing the calculated ϕ and ρ0 values for each MnO6 octahedron in a 10 × 10 × 8 supercell of LaMnO3 at room temperature, as obtained from reverse Monte Carlo analysis of neutron PDF data collected by Thygesen et al. (2017). In panel (c), orthogonal axes were used (i.e. angular distortion was included in the calculation, using the method described in this manuscript), whereas in panel (d) the Mn—O bond directions were taken as the axes, regardless of orthogonality. (e) The Mn^3+ octahedron which exhibits a mixed Q2–Q3-type distortion due to the first-order JT 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 lengths within each octahedron in the 10 × 10 × 8 supercell, with the blue vertical lines indicating the bond lengths in the average structure.

In this case study, we take a 10 × 10 × 8 supercell of LaMnO₃, obtained using RMCProfile against total scattering data obtained at room temperature and previously published in the aforementioned work (Thygesen et al., 2017). The results are shown in Fig. 6. We repeat the analysis of this supercell from the perspective of the E_g(Q₂, Q₃) Van Vleck distortion modes, using two different approaches: (i) applying the algorithm for automatically determining a set of orthogonal axes to each octahedron individually, and (ii) following the Van Vleck equations (23) and (24) proposed by Kanamori (1960) 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 Figs. 6(c) and 6(d), there are two clusters of octahedra within the polar plot, occurring at ϕ ≃ ±107°. This corresponds to occupation of the d_y² orbitals (+) and of the d_x² orbitals (−). In both cases, the superposition of perpendicular Q₃ compression and elongation modes results in an octa­hedron with mixed Q₂–Q₃ character. This finding is consistent with previous work which placed MnO₆ octahedra from LaMnO₃ onto the framework of an E_g(Q₂, Q₃) polar plot (Zhou & Goodenough, 2008a; Zhou et al., 2011).

Fig. 6(e) shows the MnO₆ octahedron in the average structure of LaMnO₃ at room temperature, with the three different bond lengths plotted in Fig. 6(f) along with a histogram of all the bond lengths in the supercell. This shows how the combination of the Q₂ and Q₃ distortion modes manifests in the octahedral distortion.

The Q₂ contribution to the distortion, as seen from the three different Mn—O bond lengths in LaMnO₃, is also present in JT-distorted ACuF₃ (A = Na, K, Rb) (Lufaso & Woodward, 2004; Marshall et al., 2013; Khomskii & Streltsov, 2021) and even in some JT-undistorted perovskites (Zhou & Goodenough, 2008a), indicating it is related to the structure. It is not intrinsic to JT-distorted manganates, as it is absent in high-spin d⁴ Mn³⁺ with edge-sharing octahedral interactions and col­linear orbital ordering such as α-NaMnO₂ and LiMnO₂ [checked using ICSD (Inorganic Crystal Structure Database, FIZ-Karlsruhe, Germany; https://icsd.fiz-karlsruhe.de/index.xhtml) references 15769 (Jansen & Hoppe, 1973) and 82993 (Armstrong & Bruce, 1996) respectively]. The Q₂ component to the octahedral distortion is therefore probably intrinsic to the crystal structure (Zhou & Goodenough, 2006, 2008a), which occurs as a result of octahedral tilting reducing the symmetry from cubic to Pnma. In LaMnO₃, the combination of the Q₂ component to the distortion and the orbital ordering [Fig. 6(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 Lufaso & Woodward (2004).

Finally, we also calculate the Q₄ to Q₆ octahedral shear modes for all octahedra in the supercell, presented as a histogram in Fig. S2. We present the average and standard deviation, calculated assuming orthogonal axes and with the automated octahedral rotation: Q₄ = −0.02 ± 0.13 Å, Q₅ = 0.02 ± 0.10 Å and Q₆ = −0.00 ± 0.11 Å. In each case, the magnitude of the distortion is zero within the standard deviation, and the value from the average structure presented in Table 2 also falls within the range of error. This low level of shear generally supports the validity of calculating the E_g(Q₂, Q₃) Van Vleck modes along bond directions rather than a Cartesian coordinate system for a system like LaMnO₃. It is interesting that the standard deviation is higher for Q₄, which quantifies the shear within the plane in which there is orbital ordering.

4.3. Effect of pressure on the JT distortion in NaNiO₂

In recent years, there have been several studies looking at the effect of applied pressure on the JT distortion in crystalline materials (Åsbrink et al., 1999; Loa et al., 2001; Choi et al., 2006; Zhou et al., 2008, 2011; Aguado et al., 2012; Mota et al., 2014; Caslin et al., 2016; Zhao et al., 2016; Collings et al., 2018; Bhadram et al., 2021; Lawler et al., 2021; Scatena et al., 2021; Ovsyannikov et al., 2021; Nagle-Cocco et al., 2022). Most of these have shown that, as a general rule, pressure reduces the magnitude of the JT distortion as a consequence of the elongated bond being more compressible than the shorter bonds.

Zhou et al. (2011) used Van Vleck modes to quantify the effect of pressure on the JT distortion in the corner-sharing perovskite-like compounds LaMnO₃ and KCuF₃. While the application of pressure reduces the magnitude of the distortion, as quantified using ρ₀ [equation (11)], they argue that it does not change the orbital mixing ϕ [equation (12)]. KCuF₃ has similar orbital ordering to LaMnO₃, except the degeneracy is due to the d⁹ hole rather than an electron. The variable-pressure crystal structures for KCuF₃ are available from the ICSD (catalogue codes 182849–182857) and are utilized here.

We previously studied the effect of pressure on the JT distortion in NaNiO₂ (Nagle-Cocco et al., 2022) by performing Rietveld refinement of neutron diffraction data from the PEARL instrument (Bull et al., 2016) at the ISIS Neutron and Muon Source (Oxfordshire, UK). However, we did not utilize the Van Vleck distortion modes, instead quantifying the JT distortion using the bond-length distortion index (Baur, 1974) and the effective coordination number (Hoppe, 1979). In that study, we found no deviation from the ambient-pressure space group C2/m (Dick et al., 1997; Sofin & Jansen, 2005), shown in Fig. 7(a), for all pressure points at room temperature up to ∼4.5 GPa. This space group permits only four short (long) and two long (short) bonds or six equal bond lengths, depending on the angle β, and so throughout the measured pressure range there is no Q₂ character to the JT distortion, consistent with the principle that hydrostatic pressure does not change orbital mixing (Zhou et al., 2011).

Figure 7 The crystal structures of (a) JT-active C2/m NaNiO2 and (b) inactive Fe2O3. (c) and (d) Comparisons of various metrics for quantifying the degree of JT distortion as a function of pressure, for NiO6 octahedra in NaNiO2 and FeO6 octahedra in Fe2O3, respectively. The parameters subject to comparison are the magnitude ρ0, the bond-length distortion index, the effective coordination number and quadratic elongation. Dashed lines indicate a linear fit to the data, whereas solid lines connect data points.

Here, we perform a fresh analysis of the variable-pressure octahedral behaviour as a function of pressure at room temperature in NaNiO₂ in terms of the E_g(Q₂, Q₃) Van Vleck distortion modes. For a reference we sought a material that does not exhibit a first-order JT distortion but does exhibit bond-length distortion; for this purpose, we selected Fe₂O₃, the pressure dependence of which was previously studied by Finger & Hazen (1980) and which exhibits bond-length distortion due to its face- and edge-sharing octahedral connectivity. Fe₂O₃ contains high-spin d⁵ Fe³⁺ cations within octahedra which interact via both face- and edge-sharing interactions. Note that Fe₂O₃ probably exhibits some very subtle pseudo-JT distortion (related to, but distinct from, the first JT effect discussed here) on account of the Fe³⁺ ions (Cumby & Attfield, 2017; Bersuker & Polinger, 2020), but this does not impact the discussion in any meaningful way.

In Fig. 7(c) we compare (for NaNiO₂) ρ₀ with three other parameters (bond-length distortion index, quadratic elongation and effective coordination number) which are often used to parameterize the magnitude of the JT 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 JT distortion. This can be compared with Fig. 7(d), which shows the same parameters for the JT-undistorted FeO₆ octahedra in Fe₂O₃; it can be seen that ρ₀ 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₃ can be seen in Fig. S3). This means that, while these parameters are valid for quantifying the magnitude of the JT distortion, they are also sensitive to other kinds of distortion. ρ₀ is calculated using Q₂ and Q₃ which have E_g symmetry, so ρ₀ will only be non-zero for a distortion with E_g symmetry. Thus, it is arguably the ideal choice for parameterizing the magnitude of this type of JT distortion. However, while ρ₀ is more reliable than the other parameters shown in Figs. 7(c) and 7(d) for demonstrating the presence of JT distortion, it is not always strictly zero for a JT-inactive octahedron as it will have a non-zero value if the octahedron is distorted with an e_g symmetry. For example, the NaO₆ octahedron in C2/m NaNiO₂ has the same symmetry as the NiO₆ octahedron and so exhibits a value of ρ₀ between 0.065 and 0.05 within the studied pressure range (Fig. S4), and JT-inactive FeO₆ octahedra in RFeO₃ perovskites have non-zero ρ₀ due to the E_g symmetry of the distorted octahedra, as shown by Zhou & Goodenough (2008a).

Fig. 8 shows a polar plot for the behaviour of NaNiO₂ and KCuF₃ in the range 0–5 GPa (the measured range for NaNiO₂). It can be seen that within this pressure range, the magnitude of the JT distortion decreases far more for KCuF₃ than for NaNiO₂; this reflects the fact that KCuF₃ is more compressible, with a bulk modulus of 57 (1) GPa (Zhou et al., 2011) in contrast to 121 (2) GPa for NaNiO₂ (Nagle-Cocco et al., 2022), as obtained by a fit to the third-order Birch–Murnaghan equation of state (Birch, 1947). Within this pressure range we see that ϕ does not change with pressure for either material and that this property is true regardless of whether ϕ is or is not a special angle (as in Table 1), consistent with the interpretation of Zhou et al. (2011).

Figure 8 An Eg(Q2, Q3) radial plot comparing the pressure dependence of the MO6 (M = Ni, Cu) octahedra for KCuF3 and NaNiO2 between 0 and 5 GPa, where ρ0 is normalized to the value at the lowest measured pressure and the dashed lines represent the average ϕ for each material within this pressure range.

Finally, in the previous study (Nagle-Cocco et al., 2022) we showed, using specific O—Ni—O bond angles, that pressure reduces the angular distortion for NaNiO₂. Here, we show that pressure also reduces the related shear distortion in NaNiO₂. This is demonstrated in Fig. 9 where we plot the octahedral shear Q₄, Q₅ and Q₆ modes for NaNiO₂ and Fe₂O₃ against the bond-angle standard deviation σ_ζ, defined in equation (13). Unlike the AlO₆ octahedra in LaAlO₃ (Fig. 5), for NiO₆ octahedra in NaNiO₂ there is no perfect correlation between the shear modes and angular distortion despite η ≃ 1, 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₂O₃ the shear fraction η ≪ 1 and there is no correlation between the shear distortion modes and angular distortion. This difference in behaviour probably arises because the main driver of the change is a continuous decrease in the JT distortion in NaNiO₂, while in Fe₂O₃ the 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 η we propose in this work.

Figure 9 The pressure dependence of the shear and angular distortion in (a) JT-distorted NiO6 octahedra in NaNiO2 and (b) JT-undistorted FeO6 octahedra in Fe2O3. Shear distortion is represented with the Q4, Q5 and Q6 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. η is the angular shear fraction defined in equation (22). Note that for the NiO6 octahedra, Q5 = Q6, whereas for FeO6 octahedra, Q4 = Q6 = −Q5. For Fe2O3, the average position of the O ligands was taken as the centre of the octahedron.

5. Conclusion

We have presented 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 JT distortion, and we have shown that the parameter ρ₀ (which is based on the Van Vleck Q₂ and Q₃ modes) is a more reliable way of quantifying the JT distortion than other oft-used parameters such as the bond-length distortion index.

We have shown 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₄, Q₅ and Q₆ modes to be zero. We have shown that there is value in calculating these latter modes, for instance in understanding the effect of octahedral tilting on octahedra in perovskite-like materials. These shear modes will also be useful for parameterizing the JT effect when the degeneracy occurs in the t_2g orbitals and results in a trigonal distortion, because their symmetry matches the distortion.

We have also shown 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₃) or first-order JT distortion (as in NaNiO₂). However, there is no correlation when the distortion is caused by competing interactions due to face- or edge-sharing octahedra (as in Fe₂O₃). We propose a new parameter, the shear fraction η [defined in equation (22)], which can be used to predict whether there will be correlation between octahedral shear modes and angular distortion.

6. Related literature

For further literature related to the supporting information, see Halasyamani (2004), Koçer et al. (2019) and Swanson & Peterson (1980).

Footnotes

¹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.

²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 ρ₀, but will affect the value of = arctan(Q₂/Q₃) by an integer multiple of 120°, plus a reflection about the nearest special angle (see Table 1) if there is Q₂–Q₃ mixing.

³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.