Impact of static disorder and dynamic disorder on the thermal conductivity of sodium superoxide (NaO₂)

Crystal structures that contain molecules possess rotational degrees of freedom (DOF) that influence their vibrational and thermal properties.^1,2 While such crystals have a well-defined lattice and basis, disorder can emerge as temperature increases and thermal fluctuations overcome rotational energy barriers. At low temperatures, the molecules vibrate and librate (i.e., exhibit small rotations) around their equilibrium positions, maintaining an orientationally ordered configuration. We refer to this state as the librating regime. At high temperatures, the molecules may freely rotate, in what we refer to as the dynamic disorder regime. Between the librating and dynamic disorder regimes, some crystals take on a structure where the molecules transition between different orientations on a time scale larger than the librational period. We refer to this behavior as the static disorder regime.

As an example, consider the C $60$ molecular crystal. At temperatures below 260 K, experimental measurements indicate that the C $60$ molecules are locked into an orientationally ordered configuration where they vibrate and librate.^3,4 In this regime, thermal conductivity decreases with increasing temperature, resembling the behavior of a typical atomic crystal. At temperatures above 260 K, the molecules freely rotate³ and thermal conductivity is measured to be temperature independent, as in an amorphous material, with a measured value of 0.4 W/m K.⁴ These two thermal conductivity trends have also been observed in molecular dynamics (MD) simulations of the C $60$ molecular crystal.^5–8 The simulations of Liang et al. suggest the existence of an intermediate static disorder regime, where the thermal conductivity is also temperature independent.⁶ 

Similar thermal conductivity trends have been measured in superatomic crystals built from C $60$ and chalcogen-based superatoms.^9,10 Using MD simulations, Liang et al. found that a crystal-like thermal conductivity only exists in the librating regime for superatomic crystals where the C $60$ molecules are orientationally ordered.⁶ An amorphous-like thermal conductivity is found when there is an orientational disorder in C $60$ molecules, be it in the librating, static disorder, or dynamic disorder regimes.

Another material whose thermal properties are impacted by rotational DOF is methylammonium lead halide (MAPbI $3$⁠), a metal halide perovskite of interest in photovoltaic and thermoelectric applications.^11–14 MD simulations have been conducted to calculate the room temperature thermal conductivity of the tetragonal phase, where the MA $+$ ions display dynamic disorder. Hata et al.’s treatment of the MA $+$ ions as being fully flexible led to a thermal conductivity of 0.18 W/m K.¹⁵ When they treated the MA $+$ ions as diatomic molecules, point masses, or removed them altogether, larger thermal conductivities of 0.21, 0.30, and 0.33 W/m K were obtained. The largest increase coincided with the removal of rotation. Conversely, Caddeo et al. found that thermal conductivity only increased notably when the internal degrees of freedom of the MA $+$ ions were removed (i.e., rotation was not the defining factor).¹⁶ Giri et al. found that the temperature-independent thermal conductivity at higher temperatures can be attributed to rotations of the MA $+$ ions that disappear under compression.¹⁷ The impact of rotational DOFs on thermal transport in MAPbI $3$⁠, thus, remains an open question.

Static disorder and dynamic disorder are also found in alkali metal superoxides (i.e., XO $2$ where X = Na, K, Rb, Cs) and alkali metal cyanides (e.g., NaCN). Minimal attention has been given to studying the vibrational and thermal transport properties of such crystals.

We, herein, investigate sodium superoxide (NaO $2$⁠), a crystal built from Na $+$ ions and bonded O $2 −$ ions that displays the regimes of libration, static disorder, and dynamic disorder. NaO $2$ recently emerged as a material of interest after it was observed as a discharge product of sodium-air batteries.^18–20 As NaO $2$ can be recharged at lower over-potentials than the expected peroxide (Na $2$O $2$⁠), this finding suggests the possibility of highly reversible batteries. Understanding the vibrational structure and thermal transport in NaO $2$ will enable the specification of working conditions that favor its production.

Early research efforts on NaO $2$ focused on its structure and phase transitions.^21,22 Below a temperature of 196 K, the crystal displays a marcasite structure. At temperatures above 196 K, NaO $2$ adopts a pyrite structure with a 12-atom unit cell, characterized by a simple cubic lattice formed by the centers of mass of the ions. In between temperatures of 196 and 223 K, each O $2 −$ ion aligns along one of the four cube diagonals. Due to thermal fluctuations, the O $2 −$ ions transition between the four orientations. This structure, thus, displays static disorder, as depicted in Fig. 1(b). Above a temperature of 223 K, the O $2 −$ ions freely rotate, leading to a dynamically disordered structure, as shown in Fig. 1(c).

Fracassi and Klein conducted MD simulations of the disordered pyrite phase and lattice dynamics calculations of the ordered pyrite phase using an interatomic potential with rigid O $2 −$ ions.²³ The calculated acoustic phonon frequencies were in fair agreement with the previous experimental data.²⁴ Measured and calculated dispersion curves for NaO $2$^21,23,25 display an acoustic branch with an upward curvature near the Gamma point, which arises due to translational–rotational coupling of the O $2 −$ ions. Kang et al. reported density functional theory-based zero-temperature dispersion curves for the pyrite structure.²⁶ They found a large number of imaginary frequencies because they did not include finite temperature effects in force constant calculations. We are not aware of any work that has measured or calculated the thermal conductivity of NaO $2$⁠.

The objective of this study is to calculate the thermal conductivity of pyrite NaO $2$ in its librating, static disorder, and dynamic disorder regimes, as shown in Figs. 1(a)–1(c). The transitions between the regimes will be induced by increasing (i) the temperature and (ii) the lattice constant, so as to investigate the impacts of thermal fluctuations and the rotational energy barriers. The rotational dynamics that accompany the static disorder and the dynamic disorder challenge the standard use of lattice dynamics calculations and thermal transport described using a phonon or Wigner formulation,^27,28 which require that the constituent particles vibrate and/or librate around equilibrium positions. We, therefore, use MD simulations, which make no assumptions about the nature of the thermal transport, to calculate the thermal conductivity. MD simulations are classical and, thus, do not include quantum effects on the vibrational mode populations or heat capacity, leading to a temperature-independent heat capacity.²⁹ Due to the large system sizes and simulation times required, the MD simulations will be driven by empirical potentials. We do not believe that the use of MD simulations or empirical potentials will impact our ability to resolve the fundamental nature of thermal transport in NaO $2$⁠, although they may introduce some errors when compared to future experimental measurements.

The remainder of this paper is organized as follows. The simulation details are presented in Sec. II A, including how the rigid bodies are modeled. The methodologies for quantifying the structural dynamics, predicting thermal conductivity from MD simulations, and performing lattice dynamics calculations are provided in Secs. II B–II D. The MD results are presented in Secs. III A and III B, where we explore the impact of rotations on the thermal conductivity for a range of temperatures and lattice constants. We further analyze thermal transport in the dynamic disorder regime using instantaneous normal modes and Allen–Feldman theory in Sec. III C. The findings are summarized in Sec. IV.

The MD simulations are performed using large-scale atomic/molecular massively parallel simulator (LAMMPS)³⁰ at temperatures between 50 and 400 K with periodic boundary conditions and a time step of 1 fs. To allow for a consistent interpretation, the pyrite structure is considered at all temperatures. The O $2 −$ ion is treated as a rigid body with a bond length of 1.31 Å. The rigid assumption is reasonable because the ion’s intramolecular vibrational frequency (34.6 THz)²³ is much higher than the thermal vibrational frequency $k BT/h=8.33$ THz at a temperature of 400 K, where $k B$ is the Boltzmann constant, $T$ is the temperature, and $h$ is the Planck constant, indicating that the intramolecular mode will not be activated at the temperatures of interest.

A challenge in modeling thermal transport in NaO $2$ using LAMMPS is that it contains point masses and linear rigid bodies. The SHAKE and RATTLE algorithms, which are commonly used for small rigid molecules, cannot constrain a bond angle to be 180  $°$⁠. Furthermore, the RIGID algorithm, which must then be used, cannot calculate the heat flux required for the Green–Kubo thermal conductivity calculations. To address this issue, we made two modifications to the RIGID algorithm. First, as described in Ref. 5, the heat flux calculation was augmented to include the angular velocities and torques associated with rigid bodies. Second, the minimum number of atoms per rigid molecule was decreased to one. This change treats both the non-bonded Na $+$ ions and bonded O $2 −$ ions as rigid bodies. We validated the second modification by (i) calculating the thermal conductivity of crystalline Lennard-Jones argon and comparing to a standard calculation (Sec. S3 in the supplementary material) and (ii) confirming that the Na $+$ ions had zero rotational kinetic energies (KEs) and torques during the simulations. The modified scripts are publicly available.³² 

We built a cubic simulation box with $4×4×4$ unit cells, corresponding to 256 Na $+$ ions and 256 O $2 −$ ions. This box size was sufficient to remove the size effects in thermal conductivity calculations (Sec. S4 in the supplementary material). A lattice constant of 5.46 Å was used for the temperature series.²³ For the density series, lattice constants between 5.10 and 5.70 Å were used at a temperature of 150 K. The system is brought to the desired temperature by running the simulation in the NVT ensemble (i.e., at constant mass, volume, and thermodynamic temperature) with a Nosé–Hoover thermostat for 10 $5$ time steps followed by 10 $5$ time steps in the NVE ensemble (i.e., at constant mass, volume, and total energy). The production run consists of 10 $6$ further time steps in the NVE ensemble.

The root mean square displacements (RMSDs) of the Na $+$ ions and the individual oxygen atoms were calculated every ten time steps in the production run using the average atomic positions over this time as the references. The reported RMSDs are the average over all particles of a given species and over the production run.

As will be shown in Sec. III B, the NaO $2$ thermal conductivity in the dynamic disorder regime is temperature independent. This behavior is typical of disordered solids at high temperatures.³⁷ To investigate this behavior, we will apply harmonic lattice dynamics-based calculations that are typically performed for disordered solids.

A solid’s vibrational spectrum is determined from its equilibrium atomic positions, where the net force on each atom is zero. A dynamical matrix is built from the second-order (i.e., harmonic) force constants. Its eigenvalues are the squares of the mode frequencies and its eigenvectors provide the relative motions of the atoms in each mode. For a crystal, the analysis is performed at the unit-cell level. The modes are phonons, which are propagating and non-localized. For a disordered system, the analysis is performed on a supercell. The eigenvectors can be used to broadly classify the modes as propagons (phonon-like), diffusons (non-localized, non-propagating), or locons (localized, non-propagating).^38–40 From the frequencies, one can build the density of states (DOS) for any solid, and the phonon dispersion for a crystal.

In the dynamic disorder regime of NaO $2$⁠, the O $2 −$ ions freely rotate and do not have equilibrium positions. The vibrational spectrum can be calculated, however, by building a dynamical matrix based on the positions at a moment in time (i.e., a snapshot). This procedure is called instantaneous normal mode (INM) analysis.^41–43 We calculate the INMs using GULP⁴⁴ for snapshots taken from the MD simulations and average over five independent snapshots (sampled every 100 fs) to get good statistics in the subsequent calculations.

The oxygen atom and Na $+$ ion RMSDs for a lattice constant of 5.46 Å are plotted vs temperature in Fig. 2. Three distinct regimes are evident for the oxygen atoms, which we interpret with additional information obtained by visualizing the simulations. Below a temperature of 110 K, the RMSD scales as T $0.5$, as expected in the harmonic crystal approximation for vibrating and/or librating atoms.⁴⁶ Between temperatures of 110 and 200 K, the RMSD increases rapidly with increasing temperature. This is the static disorder regime, where the O $2 −$ ions occasionally rotate to a different equilibrium position. At and above a temperature of 220 K, the O $2 −$ ions freely rotate in the dynamic disorder regime. This temperature is close to the experimentally observed transition temperature of 223 K.²³ In this regime, the RMSD is slightly larger than half of the O–O bond length and continue to increase slowly due to the vibrations, following a $T 0.045$ dependence. This transition in the RMSD behavior from the librating regime to the dynamic disorder regime is similar to that observed for the C $60$ molecular crystal and superatomic crystals.⁶ 

The Na $+$ ion RMSD also follows the T $0.5$ scaling in the librating regime. In the static disorder regime, the Na $+$ ion RMSD continues to scale as T $0.5$. Because the time between O $2 −$ ion transitions [ $O$(10 ps), Sec. S10 in the supplementary material] is much longer than the atomic vibration time scale [ $O$(1 fs)], the dynamics of the Na $+$ ions are not affected. It is only when the O $2 −$ ions approach free rotation at a temperature of 220 K that the Na $+$ ions experience a change in their potential energy landscape, providing them with more space to move and a sudden increase in their RMSD. The Na $+$ ion RMSD continues to increase slowly in the dynamic disorder regime, following a $T 0.23$ dependence.

The ion-averaged oxygen RACFs at a lattice constant of 5.46 Å in the librating and dynamic disorder regimes are plotted in Figs. 3(a) and 3(b) for a range of temperatures. RACFs for individual ions in these two regimes are plotted in Figs. S3(a) and S3(c) in the supplementary material and show similar behavior to the averaged RACFs.

In the librating regime, the RACF converges to a value just below unity that decreases with increasing temperature. This result is consistent with the ions exhibiting small amplitude oscillations and maintaining a crystal-like, orientationally ordered structure. The Fourier transform of the RACF at a temperature of 50 K in this regime is plotted in Fig. S4 in the supplementary material. The spectrum shows broadband activity below frequencies of 6 THz. This range aligns with the acoustic and low-level optical phonon branches in the dispersion curves, which are plotted in Fig. S10 in the supplementary material.

In the static disorder regime, between temperatures of 110 and 200 K, the RACF curves for individual ions do not show consistent behavior. As shown in Fig. S3(b) in the supplementary material, the individual RACF curves transition between +1 and $−$1, indicating that the ions switch between states that are parallel and anti-parallel to the initial state (i.e., parallel or anti-parallel to the same body diagonal). Because these transitions happen at different times for different ions, averaging the individual ion RACFs does not lead to any meaningful insight.

The calculated thermal conductivities for a lattice constant of 5.46 Å are plotted in Fig. 4 vs temperature. In the librating regime, where NaO $2$ is fully crystalline, thermal conductivity decreases from 4 to 1 W/m K as the temperature increases from 50 to 110 K. This behavior is typical of electrically insulating crystals where the dominant mode of thermal transport is through phonons⁴⁷ and was also observed in the librating regimes of superatomic crystals built from C $60$ and chalcogen-based superatoms.^5–10 As temperature increases, phonon mode population also increases, resulting in greater phonon–phonon scattering and a reduction in thermal conductivity.

Above a temperature of 200 K, in the dynamic disorder regime, thermal conductivity takes on a uniform value of 0.7 W/m K. This temperature independence is typical of amorphous materials at high temperatures,^39,48,49 where it is a result of their structural disorder limiting vibrational lifetimes. We hypothesize that the free rotations and the accompanying disorder lead to a uniform thermal conductivity. This hypothesis is investigated in Sec. III C.

In the static disorder regime, which begins at a temperature of 110 K, the thermal conductivity initially decreases to 0.7 W/m K and then plateaus above a temperature of 180 K. We attribute the initial thermal conductivity decrease to the presence of a small amount of orientational disorder in the librating, crystalline system. The RACF calculations at a temperature of 160 K, shown in Fig. S3(b) in the supplementary material, show that only a few of the O $2 −$ rotate to a different orientation, supporting this interpretation. Between temperatures of 50 and 180 K, the thermal conductivity shows a $T − 1.38$ scaling, which is shown in Fig. 4. A stronger scaling than the $T − 1$ behavior associated with three-phonon scattering in a classical system suggests the presence of higher-order scattering events.

To further investigate the effect of changing the lattice constant, we conducted MD simulations where it was varied between 5.10 and 5.70 Å at a temperature of 150 K. We used RMSDs (Fig. 5), RACFs (Fig. S5 in the supplementary material), and visualizations of the simulations to distinguish between the librating, static disorder, and dynamic disorder regimes.

The calculated thermal conductivities are plotted in Fig. 6. Thermal conductivity monotonically decreases with increasing lattice constant, as would be expected based on the weakening interactions between the ions. In the librating regime, thermal conductivity decreases from 1.5 to 1.4 W/m K as the lattice constant increases from 5.10 to 5.25 Å.

An O $2 −$ ion’s ability to rotate depends on the associated energy barrier and the thermal energy in the system. To aid in interpreting the results in the static disorder regime, we estimated the O $2 −$ ion rotational energy barriers, $E a$⁠. The details of the calculation are provided in Sec. S8 in the supplementary material and the results are plotted in Fig. 6. The energy barrier decreases from $0.47±0.01$ to $0.35±0.02$ eV as the lattice constant increases from 5.25 to 5.60 Å. In contrast, the energy barrier in the static disorder regime in the temperature series at 140, 160, and 180 K is constant (⁠ $0.41±0.01$⁠, $0.42±0.02$⁠, and $0.42±0.02$ eV, respectively).

In the static disorder regime for the lattice constant series, the thermal conductivity decreases from 1.4 to 0.6 W/m K. This decrease is larger than the decrease observed in the static disorder regime in the temperature series. Examination of the oxygen atom RMSD shows a similar behavior in the static disorder regime in both the temperature and the lattice constant series, with the value increasing quickly from 0.20 to 0.72 Å, even though the origin of the increase is different (i.e., higher temperature vs more free space). In the temperature series, the interaction strengths remain the same and the vibrational mode populations increase as the rotations emerge. In the lattice constant series, the interaction strengths decrease and the vibrational mode populations remain the same as the rotations emerge. These results suggest that the impact on the thermal conductivity of the interaction strengths is stronger than that of mode populations.

Above a lattice constant of 5.60 Å, in the dynamic disorder regime, the thermal conductivity continues to decrease as the strength of the interactions decreases.

It is notable that the Na $+$ RMSD, plotted in Fig. 5, does not change as the lattice constant increases across the three regimes. This observation suggests that even though there is more space for the Na $+$ ions to occupy, their vibrations are limited by their thermal energy.

Analyses in Secs. III A and III B led to the conclusion that NaO $2$ exhibits characteristics typical of an amorphous material in its dynamic disorder regime. To further explore thermal transport in this regime, we now turn to the INM analysis outlined in Sec. II D. The calculations were performed at a temperature of 300 K for a lattice constant of 5.46 Å. The results are minimally changed when considering other temperatures in the dynamic disorder regime.

The dynamic disorder DOS is plotted in Fig. 7(a). Also plotted is the DOS for the pyrite crystal at zero temperature. The dynamic disorder DOS scales as the square of the frequency up to 2.4 THz, which is within the acoustic phonon part of the crystal phonon dispersion (Fig. S10 in the supplementary material). This scaling is predicted under the Debye approximation for a crystal, which assumes linear dispersion curves, suggesting that these low-frequency modes are propagating (i.e., phonon-like). As the frequency increases, the dynamic disorder DOS contains one broad peak, which is in contrast to the sharp features in the crystal DOS.

The participation ratios (PRs) [Eq. (6)] at a temperature of 300 K are plotted in Fig. 7(b). 87% of the PRs are greater than 0.1, suggesting non-localized modes. Smaller values, which correspond to localized modes, occur at the highest frequencies, consistent with the calculations on amorphous materials (e.g., amorphous calcium silicate hydrates,⁴⁸ amorphous silica,⁴⁹ amorphous silicon,³⁹ and amorphous carbon⁴⁹). The DOS and PR results in Figs. 7(a) and 7(b) suggest that there is a significant spectrum of non-localized and non-propagating modes (i.e., diffusons) in the dynamic disorder regime of NaO $2$ that will contribute to its thermal conductivity.

In previous investigations of disordered systems,^48,50 a transition between propagons and diffusons was identified by examining the frequency scaling of the mode lifetimes. The lifetimes in these studies were obtained from normal mode analysis in MD simulations. This analysis is not possible here as the atoms do not have well-defined equilibrium positions. It is, thus, challenging to establish a propagon–diffuson transition point.

The diffusivities of all vibrational modes at a temperature of 300 K are plotted in Fig. 7(c). The diffusivity decreases as the frequency increases up to 2.4 THz, corresponding to the range where we believe the modes are propagating based on the DOS. Between frequencies of 2.4 and 11 THz, the diffusivities vary between 10 $− 9$ and 10 $− 8$ m $2$/s and then show a pronounced drop. This drop is consistent with the localization of the high frequency modes shown in Fig. 7(b). We categorize these intermediate frequencies as diffusons, whose contribution to thermal conductivity can be calculated from Eq. (8). The diffusivity trends are comparable to other disordered materials (e.g., amorphous silicon and amorphous silica,⁵⁰ and amorphous carbon⁵¹).

Based on the dynamic disorder DOS and diffusivities, we chose three transition frequencies (2.4 THz, 2.7 THz, and 3.0 THz) for calculating the AF thermal conductivity. These frequencies are at the high end of the acoustic phonon modes in the perfect crystal (Fig. S10 in the supplementary material), where the propagon–diffuson transition typically occurs.^50,51 The cumulative AF thermal conductivity for these cutoffs normalized by the Green–Kubo value in the dynamic disorder regime (0.7 W/m K) is plotted in Fig. 7(d). The total AF thermal conductivities at the three cutoffs are $0.63±0.01$ W/m K (2.4 THz), $0.59±0.01$ W/m K (2.7 THz), and $0.56±0.01$ (3.0 THz), values that are all less than the Green–Kubo value. From these calculations, we conclude that 75%–85% of the total thermal conductivity in the dynamic disorder regime can be attributed to diffusons. This result aligns with the findings on diffuson contribution in amorphous materials (e.g., 67% in amorphous silicon,⁵⁰ 94% in amorphous silica,⁵⁰ 70% in calcium silicate hydrates,⁴⁸ and up to 60% in low-density amorphous carbon⁵¹).

We conducted equilibrium MD simulations and lattice dynamics calculations to explore the impact of temperature and lattice constant on the rotational dynamics and thermal transport in the pyrite phase of NaO $2$⁠. Three distinct regimes were observed for the dynamics of the O $2 −$ ions: librating, static disorder, and dynamic disorder (i.e., freely rotating) (Fig. 1).

In the temperature series, as shown in Fig. 4, a crystal-like thermal conductivity (i.e., decreasing with increasing temperature) manifests in the presence of orientational order when the O $2 −$ ions are librating. An amorphous-like thermal conductivity (i.e., temperature-independent) arises in the presence of orientational disorder, which is found in the static disorder and dynamic disorder regimes. Decomposition of the thermal conductivity calculated from the Green–Kubo method indicates that the rotations do not significantly contribute to thermal transport at any temperature. In the dynamic disorder order regime, INM analysis suggests that 80% of the thermal conductivity is contributed by diffusons [Fig. 7(d)], modes that are non-propagating and non-localized.

As the lattice constant is increased, the O $2 −$ ions pass through the same three structural regimes. As shown in Fig. 6, the thermal conductivity decreases monotonically in this series, due to weakening interaction strengths in all three regimes and the emergence of rotations due to a decreasing rotational energy barrier in the static disorder regime.

The MD simulations and INM analysis approaches presented here can be applied to other crystals that contain rotational degrees of freedom (e.g., alkali metal cyanides, fullerenes, and hybrid organic-inorganic perovskites). Our results also suggest the potential for thermal conductivity tuning and switching by controlling molecular rotational dynamics through external stimuli such as light, stress, and electric fields.

See the supplementary material for information about the interatomic potential and O 2 − ion structure, validation of the LAMMPS modifications, thermal conductivity size effect analysis, Debye temperature calculation, additional RACFs, temperature dependence of thermal conductivity at different densities, Green–Kubo thermal conductivity decomposition, the rotational energy barrier calculation, phonon dispersion curves, uncertainty quantification in Allen–Feldmann calculations due to the broadening factor, and the calculation of the transition time in the static disorder regime.

H.R. and A.J.H.M. were supported by the National Science Foundation (NSF) (Award No. DMR-2025013). C.S. was supported by the National Natural Science Foundation of China (NNSFC) (Grant No. 52306102), the Excellent Young Scientists Fund (Overseas) of Shandong Province (2023HWYQ-107), and the Taishan Scholars Program of Shandong Province.

The authors have no conflicts to disclose.

Hariharan Ramasubramanian: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Visualization (lead); Writing – original draft (lead). Cheng Shao: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Software (equal); Validation (equal); Writing – review & editing (equal). Alan J. H. McGaughey: Conceptualization (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – review & editing (lead).

Modified fix RIGID scripts developed for this work and a sample input file for LAMMPS are publically available at https://github.com/rfhari/ThermalTransportSodiumSuperxoide, Ref. 32.