Spin–orbital superexchange physics emerging from interacting oxygen molecules in KO2

We propose that the spin–orbital–lattice coupled phenomena, widely known for the transition-metal oxides, can be realized in molecular solids, comprising of orbitally degenerate magnetic O2− ions. KO2 is one such system. Using the first-principles electronic structure calculations, we set-up an effective spin–orbital superexchange model for the low-energy molecular bands and argue that many anomalous properties of KO2 indeed replicate the status of its orbital system in different temperature regimes.


Introduction
The prospect of making new magnetic materials with new functionalities has attracted much attention to nontraditional solutions without conventional d-or f-elements [1]- [4]. Among them, the idea of molecular crystals, where a molecule participates as both structural and magnetic unit, takes a special place [5]- [9]. Unlike the atom, the molecule is an essentially nonspherical object. This nonsphericity can be that additional degree of freedom, which can be potentially used in order to control and design the properties of the new magnetic materials [10,11].
The simplest example is provided by diatomic molecules. The majority of such molecules form a nonmagnetic spin-singlet ground state. Nevertheless, there are two known exceptions: NO and O 2 [12]. The O 2 molecule has two unpaired electrons in the doubly degenerate π g level resulting in the spin-triplet ground state, as is required by the first Hund rule. If these molecules form a crystal (by either cooling or pressurizing), it may become magnetic. Such a situation is indeed realized in solid oxygen [5,6].
Since O 2 is a good oxidizer and can easily take an additional electron when it is brought into contact with alkali elements, there is another way of making crystalline O 2 , in the form of ionic crystals, which are pretty much similar to the notorious NaCl, but distorted due to the additional elastic interactions involving the oxygen molecules. Such materials do exist. One of the typical examples is KO 2 , which has been used as a chemical oxygen generator in rebreathers, spacecraft and life support systems. The magnetic and structural properties of KO 2 (and other alkali hyperoxides) were intensively studied in the 1970s [7]- [11], [13], and have become practically forgotten today.
Early investigations on KO 2 revealed a rich variety of structural and magnetic properties, which were nearly overshadowed by many unresolved problems 1 . KO 2 has five structural phase transitions [7,8]. However, the details are known only for the room-temperature bct phase (figure 1). The lowering of temperature (T ) causes some canting of the oxygen molecules, which is characterized by the angle θ formed by the oxygen molecules and the c-axis of the crystal. However, many details of this canting are also unknown. In order to interpret the temperature dependence of the inverse paramagnetic susceptibility (χ −1 m ), it is typically divided into a series of straight-line segments, corresponding to each crystallographic phase [8,13]. Then, each segment is fitted in terms of the Curie-Weiss law. The corresponding Curie-Weiss temperature (T CW ) changes from −44 K, in the high-temperature region, to 2.3 K, in the lower-temperature region, indicating some change of intermolecular magnetic interactions. However, no profound change of the magnetic susceptibility has been observed around the phase transitions, indicating that perhaps there is a more fundamental reason for the change of T CW . A sizable orbital moment persists down to 12 K, and is not quenched by the crystal distortion [8,13]. A long-range antiferromagnetic order has been detected below 7 K [8,15]. However, its type is also unknown. It was further suggested that the paramagnetic-to-antiferromagnetic transition is accompanied by a sharp (from about 20 • to 30 • ) change of the angle θ formed by the oxygen molecules and the c-axis of the crystal [8]- [11]. It was also reported that below 12 K, the magnetic field of about 70 kOe stabilizes the ferromagnetic order [13]. Left inset shows results of LSDA calculations. Right inset shows the dispersion of the π g bands. Fermi level is located at zero energy. The high-symmetry points of the Brillouin zone are denoted as = (0, 0, 0), X = (π/a, π/a, 0), N = (π/a, 0, π/c), P = (π/a, π/a, π/c) and Z = (0, 0, 2π/c) [14].
There are a number of good reasons to reconsider the situation around the alkali hyperoxides today.
1. The theoretical environment has changed dramatically over the past three decades.
Particularly, the purpose of this work will be to show that many properties of KO 2 can be successfully reinterpreted in terms of novel spin-orbital-lattice coupled phenomena, which have been overlooked in earlier studies on KO 2 and are intensively discussed in the context of transition-metal oxides today [16]- [19]. Our analysis will be based entirely on the first-principles electronic structure calculations. Namely, we will establish the lowenergy model for KO 2 , derive parameters of this model from the first-principles electronic structure calculations [20], and then apply it to the analysis of the magnetic properties of KO 2 . A similar strategy has been successfully applied recently to the series of distorted perovskite oxides [21]. 2. The experimental potentials have also risen, and could probably help us to re-examine the behavior of KO 2 in many new details today. We hope that our work will revive an interest in this activity. 3. Although KO 2 itself may look somewhat exotic today, the idea of spin-orbital superexchange (SE) interactions between molecular complexes itself is certainly more generic and can find its application in other molecular compounds. 4. The new aspect of spin-orbital SE physics in molecular systems is definitely related to the rotational degrees of freedom of the molecular complexes. Even in the single-orbital case, the rotations of the oxygen molecules may lead to a new type of magnetoelastic phenomena, known as magnetogyration [10,11]. In the multi-orbital case, the situation is even more intriguing because these rotations will also lift the degeneracy of the molecular levels and couple not only to spin, but also to the orbital degrees of freedom. An analog of this phenomenon in the conventional strongly correlated systems is the Jahn-Teller effect, which is typically caused by nonequivalent displacement of the oxygen atoms around transition-metal sites [16]. Nevertheless, such similarities with the Jahn-Teller 4 effect are only formal, because the molecular rotations may have a different origin, driving mechanisms and consequences for the properties of KO 2 and other molecular compounds.

Electronic structure in the local-density approximation (LDA)
The oxygen molecule appears to be the building block of not only crystals, but also of the electronic structure of KO 2 in the LDA [22]. The strong hybridization within the molecule leads to the formation of the molecular levels. The interaction between the molecules is considerably weaker, so that the molecular levels form a group of narrow nonoverlapping bands (figure 1) 2 . Thus, there is a clear analogy with the atomic limit in the theory of strongly-correlated systems [16], except that now, the localized electrons (or holes) reside on the molecular orbitals, which are distributed between several atomic sites. Similar to solid oxygen [5], the doublydegenerate π g band located near the Fermi level is formed by antibonding molecular p x and p y orbitals. However, in comparison with the solid oxygen, the potassium atom donates an extra electron into the π g band. Therefore, the band is 3 4 filled. Within the local-spin-density approximation (LSDA), the system is ferromagnetic and half-metallic. Therefore, the Curie temperature can be estimated in the frameworks of the double exchange (DE) model [23,24], which yields T C ∼ 320 K. 3 Thus, KO 2 is expected to be a good half-metallic ferromagnet. The obtained Curie temperature is extremely high, in direct contrast with the experimental data. Moreover, KO 2 is known to be an insulator (a yellow solid). Nevertheless, such discrepancies between LSDA and the experimental data are not surprising because LSDA is known to largely oversimplify the physics of Coulomb correlations in the narrow-band compounds. Moreover, due to the peculiar 3 4 filling of the molecular π g band, not only spin but also orbital degrees of freedom are expected to be active and contribute to spin and lattice dynamics of KO 2 [16]. These effects are totally missing in LSDA [25].

Effective low-energy models
In order to go beyond the conventional LSDA and incorporate the physics of Coulomb correlations, we first set-up an effective lattice fermion model for the low-energy π g bands: and derive all parameters of this model from the first-principles electronic structure calculations using the method proposed in [20]. The basic idea is to relate each lattice-point (i or j) to a single oxygen molecule and formulate the problem in a Wannier basis constructed from antibonding molecular p x (m = 1) and p y (m = 2) orbitals. This is the main difference from the 2 In all the calculations, we have used the experimental lattice parameters a = 4.030 Å, c = 6.697 Å and d OO = 1.306 Å, where d OO is the O-O distance in the oxygen molecule [7]. 3 The magnetic properties of half-metallic compounds can be described in terms of DE and SE interactions, which are related with first and second moments of the occupied density of states (m (1) and m (2) , respectively) [24]. The contribution of the DE and SE interactions to the Curie temperature can be easily estimated within the frameworks of the mean-field theory [23] as T DE C = −4m (1) /15k B ∼ 880 K and T SE C = −m (2) /3k B ex ∼ −560 K, where ex ∼ 1 eV is the spin splitting in LSDA, and m (1) and m (2) are calculated using the density of states of the partly occupied minority-spin π g band (figure 1). formulation which is typically adopted for conventional strongly correlated systems, where each lattice point is typically associated with a single transition-metal site. In what follows, c † iα (c iα ) in (1) creates (annihilates) an electron in the Wannier state α, specified by the combination of indices (m, σ ) representing the antibonding molecular orbital m = 1 or 2 and the spin σ =↑ or ↓. Then, the one-electron partĥ i j = h αβ i j of the model Hamiltonian (1) can be derived by using the generalized downfolding method [20,26] and starting from the electronic structure in the LDA. As usual, the intersite matrix elements ofĥ i j stand for the transfer integrals (or intersite hoppings)t i j = t mm i j , which do not depend on the spin indices. The form of these transfer integrals is explained in figure 2. Basically, all transfer interactions are restricted by the four nearest neighbors. The longer-range interactions are considerably smaller. The on-site part ofĥ i j stands for the relativistic spin-orbit interaction, formulated on the basis of |p x ↑ , |p y ↑ , |p x ↓ and |p y ↓ orbitals, and the 'crystal-field splitting', caused by the rotations of the oxygen molecules, if the latter are applied. The spinorbit coupling constant ξ , defined through the gradient of the LDA potential on the basis of the Wannier functions [21], is of the order of 17 meV. The screened Coulomb interactions in the π g band have been computed in two steps [20]. First, we derive the interaction parameters between oxygen 2p-orbitals in the atomic limit (i.e. by neglecting all the hybridization effects involving the oxygen 2p-orbitals) by using the constraint-LDA method. This yields the intraatomic Coulomb interaction u ≈ 11.37 eV, the interatomic intramolecular Coulomb interaction v ≈ 1.26 eV, and the intraatomic exchange interaction j ≈ 2.30 eV. Then, we switch on the hybridization and consider the screening of the atomic interactions in the π g band by other bands, which are constructed from the same oxygen 6 2p-orbitals 4 . This defines the full matrixÛ = U αβγ δ of the screened Coulomb interactions between antibonding molecular p x and p y orbitals. This part is done in the random-phase approximation (RPA), starting from the interaction parameters obtained in the constraint-LDA method [20] 5 . The obtained matrixÛ = U αβγ δ is fully specified by two Kanamori parameters of the intraorbital Coulomb interaction U ≈ 3.55 eV and the exchange interaction J ≈ 0.62 eV [28]. The interorbital Coulomb interaction U is related to U and J by the identity U = U − 2J .
Since the Coulomb interactions U and U are clearly the largest parameters in the problem, the Hubbard model (1) can be further converted into a spin-orbital SE model, by starting from the limit of isolated molecular levels and treating all transfer integrals as a perturbation [16,18]. Due to the 3 4 filling, it is more convenient to use the hole representation and ascribe to each molecular center a single hole spin-orbital where |a| 2 + |ã| 2 = |b ↑ | 2 + |b ↑ | 2 = |b ↓ | 2 + |b ↓ | 2 = 1. Then, the energy gain caused by the virtual hoppings in the bond i j can be computed in the second order of the perturbation theory as [18,21]: where G(α i , α j ) is the Slater determinant constructed from the hole orbitals {α i }, 6 E j M and | j M stand for eigenvalues and eigenfunctions of the excited two-hole configurations at the molecular center j, andP j is a projector operator, which enforces the Pauli principle and prevents any hoppings into the already occupied hole orbital α j [21]. Generally, the two-hole state | j M is a linear combination of several Slater determinants, and E j M takes into account the correct multiplet structure of the excited configurations. Therefore, the energy gain (3) takes into account some correlation effects beyond the mean-field approximation, for example underlying the commonly used LDA + U method [25]. These correlation effects are especially important for stabilizing the antiferromagnetic interactions [21]. Our next goal is the elimination of the orbital degrees of freedom, corresponding to the b-variables in (2), and the construction of an effective spin-only model separately for each temperature regime. This can be generally done by averaging (3) [16,18]. In principle, D(b ↑ ,b ↑ , b ↓ ,b ↓ , T ) can be derived from thermodynamic principles. However, general formulation along this line is rather complex. Nevertheless, it is always possible to consider two limiting cases, which will be investigated 4 Namely, the π u band, that is a bonding combination of the p x -and p y -orbitals, as well as the σ * g and σ * u bands, constructed from the p z -orbitals (figure 1). 5 Meanwhile, in the process of RPA calculations, we switch off all contributions to the screening associated with the transitions between π g bands in the polarization function [27]. This needs to be done in order to get rid of the parasitic metallic screening, which is always present in RPA if one starts with the metallic LDA band structure. The advantage of RPA is that in the frameworks of this approach one can easily manipulate with different channels of screening. Thus, even by starting from the LDA band structure, which is obviously far from being perfect in the case of strongly correlated systems, there is a chance to enforce the correct type of screening for the Coulomb interactions by switching off all unphysical channels. This procedure is sometimes called the constraint-RPA. 6 Note that in the limit of isolated oxygen molecules, there is only one hole residing at each molecular center. This is essentially a noninteracting system, which can be described by the single Slater determinant. 7 in detail in the next section. One is the orbital order in the limit T → 0, and another is the complete orbital disorder in the limit T → ∞.

Implications to the properties of KO 2
First, let us consider the low-temperature limit k B T ξ , 7 is fully controlled by the relativistic spin-orbit interaction. Therefore, at each lattice point, The basic idea behind this treatment is that the relativistic spin-orbit interaction is the only interaction, which directionally splits the π g levels in the limit of isolated oxygen molecules. The Coulomb interaction will additionally split the levels. However, for the one-hole state, these processes are fully isotropic and, by themselves, do not determine the direction of the splitting. Therefore, without intersite hoppings, the relativistic spin-orbit interaction fully controls the symmetry of the occupied hole orbitals. The hoppings can also contribute to the splitting of the molecular π g levels. However, in order to do so, the hole should hop to the neighboring molecule and then come back. These hoppings are penalized by the large Coulomb repulsion between the holes when they are located at the same molecular center. Therefore, the proper energy scale for these (SE) processes ist 2 /Ū , wheret 2 is the square of an effective transfer integral, constructed from the individual matrix elements of t mm i j , andŪ is of the order of U − 3J or U , depending on whether the interacting holes have the same or opposite spin. Then, using the obtained values of U, J , and the transfer integrals shown in figure 2, one can easily verify that for the individual paths, starting from some particular molecular orbital and ending with the same orbital after the hoppings on to one of the neighboring molecules,t 2 /Ū does not exceed 6 meV and is considerably smaller than the relativistic spin-orbit coupling ξ . In principle, there could be several SE paths, which contribute at once to the on-site energy. However, in order to split the π g levels and compete with the relativistic spin-orbit interaction, the symmetry of SE interactions should be broken by either spin or orbital ordering. This reduces the number of paths which can actual contribute to the π g -levels splitting. Thus, it is reasonable to start the analysis by considering the limit of large relativistic spin-orbit interaction. The main advantage of such a treatment is that it greatly simplifies the problem (with only a small loss of the accuracy) and brings all the discussions basically to the level of textbook results for the renormalized spin-wave theory. Below we will certainly verify this picture by considering more rigorous arguments without any assumptions about the superiority of the relativistic spinorbit interaction over other interactions in the system. We would also like to note that a similar strategy has been applied to the analysis of SE interactions in the limit of large crystal-field splitting [29]. 8 Since each hole orbital α is confined to the two-dimensional subspace spanned by |p + ↑ and |p − ↓ , the energies (3) can be further mapped on to the anisotropic Heisenberg model with the pseudospin 1/2: are the Pauli matrices in the basis of |p + ↑ and |p − ↓ orbitals 8 . The parameters of SE interactions can be derived by computing the energies (3) for several different orientations of the pseudospins. It is important to note that although the subspace of the occupied orbitals {α} is restricted by |p + ↑ and |p − ↓ , all four orbitals were considered in the process of virtual hoppings to the excited state in (3). For example, by denoting as (τ + kx , τ − kx ), (τ + ky , τ − ky ) and (τ + kz , τ − kz ) the pseudospin states of the kth oxygen molecule corresponding to the (positive, negative) directions parallel to the x-, y-and z-axes (see figure 2), 9 the parameters in (4) can be found as . They are summarized in figure 3. The direct exchange interactions are considerably weaker and can be neglected [30]. One can clearly see that the largest interaction is the secondneighbor coupling J 2 , which stabilizes the easy-axis ferromagnetic state. Other interactions are antiferromagnetic and frustrated on the bct lattice.
Then, the Curie temperature can be easily estimated in the frameworks of renormalized spin-wave theory as (k B T C ) −1 = q 2/ω 0 q , where ω 0 q = 2(J 0 − J ⊥ q ) is the unrenormalized magnon frequency, J 0 = k J k , and J ⊥ q is the Fourier image of {J ⊥ k } [31] 10 . At the Brillouin zone boundaries, ω 0 q is largely reduced by antiferromagnetic interactions J ⊥ k ( figure 3). Moreover, due to the frustration, there are many such excited configurations with comparable energies. This explains the relatively small value of T C ∼ 70 K. It does not seem to be fully consistent with the available experimental data [7,8,13]. However, there are also a number of factors, beyond conventional renormalized spin-wave theory, which will further reduce the theoretical T C . 10 The idea is to replace in the expression for the magnon frequency the value of spin (or pseudospin) moment by its thermal average τ z : ω q = 2 τ z (J 0 − J ⊥ q ). On the other hand, τ z is expressed through the averaged number of magnons, yielding the equation τ z q coth ω q /2k B T = 1, which is solved self-consistently for each T . 1. The hole orbitals {α i } can further relax to simultaneously minimize the energies of relativistic spin-orbit and SE interactions: Thus, contrary to the previous treatment, here we do not make any assumptions about the superiority of the relativistic spin-orbit interaction over other interactions in the system. Since ξ is finite, the SE interactions can cause some deformation of the relativistic spinorbitals. The role of this deformation can be estimated in the supercell calculations. The results of these calculations are shown by symbols in figure 3. First, we would like to note that the large-ξ limit already reproduces the main details of the spectrum of magnetic excitations. The deformation of the relativistic spin-orbitals leads only to some quantitative corrections. Nevertheless, it does reduce the gap, separating excited antiferromagnetic configurations from the ferromagnetic ground state, and T C by about 17%. 11 The largest change occurs in the X-point of the Brillouin zone, where the bct-symmetry is broken by the antiferromagnetic order. 2. Below 231 K the actual symmetry of KO 2 is lower than bct. Although the details of the low-temperature structures are largely unknown, a number of experimental facts speak in favor of uniform rotations (or canting) of the oxygen molecules [7]- [9]. We follow these experimental suggestions and consider the rigid rotations of the oxygen molecules around the crystallographic axes [0, 1, 0] and [−1, 1, 0] (figure 4) 12 . Then, we can readily explain the role of these rotations and even estimate their magnitude in different temperature  figure 1, except for X * = (π/a, −π/a, 0), N * = (0, π/a, π/c) and P * = (π/a, −π/a, π/c).
regimes. Indeed, the analysis of the magnetic susceptibility suggests that the orbital moment remains largely unquenched practically in the whole paramagnetic region [8]. This imposes a severe constraint on the upper bound of the molecular rotations, which should be smaller than 20 • . Then (and only then), the crystal-field splitting of the molecular π g levels caused by the rotations is much weaker than the relativistic spin-orbit coupling, and the structure of the magnetic interactions in KO 2 is only weakly affected by the rotations. In this regime, the magnetic properties of KO 2 can still be described by the renormalized spinwave theory, formulated on the basis of two relativistic spin-orbitals |p + ↑ and |p − ↓ , but in the rotated coordinate frame. The rotations gradually decrease ω 0 q ( figure 4) and, therefore, T C . Larger rotations make the ferromagnetic phase unstable. Thus, θ < 20 • is the distortion, which respects both relativistic spin-orbit coupling and the ferromagnetic character of interatomic magnetic interactions. Therefore, we believe that this picture of moderate rotations takes place in the wide paramagnetic region, 12 K < T < 231 K.
Nevertheless, the situation for very low temperatures, T < 12 K, must be fundamentally different. In this region, the crystal distortion appears to be large, as it is clearly manifested by almost complete quenching of the orbital magnetic moments [8]. According to the phase diagram in figure 4, such a situation should correspond to the molecular rotations θ 30 • , where the molecular levels splitting is fully controlled by the crystal-field 12

Conclusions
Using the first-principles electronic structure calculations, we have set-up an effective spinorbital model for the low-energy molecular bands of hyperoxide KO 2 . On the basis of this model, we have argued that the magnetic properties of KO 2 provide the first example of the spinorbital SE physics realized in a molecular solid. The magnetic properties of KO 2 largely depend on the orbital state of the O − 2 ions and evolve, with increasing temperature, from a quenched antiferromagnetic phase to the paramagnetic region, where the character of intermolecular interactions gradually changes from mainly ferromagnetic, and driven by the relativistic spinorbit interaction, to antiferromagnetic, and corresponding to the picture of an independent spin and orbital disorder. The geometry of molecular orbitals adds new functionalities to the classical problem of SE interactions. There are also a number of open questions. Particularly, any quantitative theory describing the reorientation of the oxygen molecules is missing. Why do the oxygen molecules tend to rotate at low temperature? Why are their orientations different in the case of KO 2 and NaO 2 [8]? These problems prompt further investigations.