Rethinking the magnetic properties of lepidocrocite: A density functional theory and cluster expansion study

Iron oxides are a common class of minerals whose applications span an array of disciplines from biotechnology to environmental science to electronics. Lepidocrocite, $γ$-Fe(III)OOH is a naturally occurring iron oxy-hydroxide that is most common in rocks, soils, and rusts.¹ While it has wide general relevance, its industrial application derives from its high temperature thermal decomposition to maghemite,² whose thin films are ubiquitous in magnetic storage, transistors, and other modern electronics applications.³ 

The magnetic properties of lepidocrocite have been studied extensively,^4–8 but there remain discrepancies between the experimentally measured behavior and what might be expected from this class of mineral. The magnetochemistry of iron oxides is governed primarily by the superexchange interactions.⁹ Superexchange provides an effective antiferromagnetic coupling across the Fe sublattice. Hence, if superexchange is the dominant effect, an antiferromagnetic ordering is expected for the ground state. The magnitude of exchange defines the magnetic coupling strength and is heavily dependent on the Fe–O–Fe bond angles and distances. Iron oxides often exhibit Néel/Curie temperatures that are commensurate with their respective exchange parameters as determined by these Fe–O–Fe bonds.¹⁰ Due to the strength of superexchange in these systems, magnetic disorder typically manifests at temperatures above 800K in the oxides and around 300K for the oxy-hydroxides.¹ The Néel temperature for lepidocrocite, however, is estimated to be between 50 and 77 K.^4,6 This is a surprisingly low value given its crystal structure, in which double chains of edge-sharing FeO$6$ octahedra along the [100] direction are connected by corner-sharing along [001] to comprise layers in its orthorhombic unit cell.¹¹ These layers are stacked along [010] and connected by hydrogen bonding. In contrast to $α$-FeOOH (goethite) and $β$-FeOOH (akaganeite), whose structures consist of double chains of iron octahedra connected by edge-sharing,¹ the bulk structure of lepidocrocite is comprised of two-dimensionally periodic edge-sharing –(Fe–O–Fe–O)– bonding interactions along both its a and c axes. Given this, it is reasonable to predict magnetic order to be maintained at higher temperatures than its polymorphic counterparts.

A number of theories have been presented to account for the experimentally measured low temperature transition of lepidocrocite. These include amorphous crystallinity,⁴ the presence of water impurities,⁵ or small concentrations of other magnetic ion inclusions.¹² Superparamagnetism arising from the nanoscale domain sizes in lepidocrocite platelets is a likely culprit, but the effect is absent in AC susceptibility measurements⁶ and size dependent Mössbauer studies.¹³ 

In this work, we seek to better define the magnetic interactions between Fe(III) in lepidocrocite in order to resolve the differences between the experimental results and the expected magnetic behavior. We accomplish this using density functional theory in combination with cluster expansion methods to compute the ground state magneto-crystalline unit cell. Using these methods, we determined the different strengths of magnetic coupling pairs as detailed in Sec. III A. As a result, we propose a new ground state magneto-crystalline unit cell symmetry in Sec. III B, the one that matches experimental data well and improves upon the magnetic descriptions used in the previous computational literature.^14–17 Using this structure, we then recalculate a host of various properties of lepidocrocite and show improved agreement with respect to the experiment. In Sec. III G, we postulate from the results of the relative calculated coupling strengths and consideration of magnetic dimensionality that the finite temperature properties of lepidocrocite can be modeled as a Heisenberg quasi-1D spin 1/2 AFM chain, leading to a good agreement with experimental results. Our findings are summarized in Sec. IV.

A physically complete theoretical description of magnetism in the condensed matter is enormously challenging. By itself, first-principles determination of the ground state magnetic unit cell, which can be larger than that of the atomic structure, is not immediately straightforward. Finite temperature properties require surmounting even larger hurdles. In this work, we approach magnetic properties through the approximation of the real system by a collinear disordered local moment (DLM) model.¹⁸ Although noncollinear, magnon, and spin wave effects are ignored, thermodynamic properties have successfully been calculated for a number of systems using this simplified model.^15,19–23

In practice, our model system entails a combination of local spin-up and spin-down configurations, creating some stoichiometric relation of $Fe(↑)xFe(↓)1−x$OOH in the bulk. By assuming these collinear states are eigenstates, the total energy of the system can then be represented as

where $σi=±1$ denotes the sign of the projection onto an arbitrary quantization axis of the electronic spin at site $i$ (the projection is assumed to have magnitude $5ℏ/2$ since we expect five unpaired electrons of the same spin at each Fe due to large intra-atomic exchange effects).

Using cluster expansion methods,²⁴ as implemented in the Alloy Theoretical Automated Toolkit (ATAT), $Jij$ of Eq. (1) can be represented as effective cluster interactions (ECIs) of an “alloy” of spin-up and -down sites. Ideally, these ECIs represent the interaction energies for all possible cluster sizes and configurations, and are exact for the complete basis. They are determined by solving the set of linear equations relating the total energy of a configuration E to the cluster interaction energy J through the cluster correlation term $∏$⁠, where $∏$ is related to the probability of finding some cluster pattern $α$ in some configuration $σ$ and E is calculated by ab initio methods. This series of interactions is truncated to some limit of clusters that give a reasonable cross-validation score, (⁠$<0.025$⁠), the result of predicted energies being within reasonable agreement of calculated energies for a series of configurations.

With a given set of ECIs, the predicted energy of any configuration can be rapidly calculated. This allows an efficient search of all possible ground state configurations up to a reasonable supercell size. This is valuable for determining the magnetic order of anti-ferromagnetic or ferrimagnetic systems in anisotropic crystal structures, where all possible configurations within even small supercell sizes would be unfeasible to calculate directly. These predicted energies can then be verified with ab initio calculations. An additional benefit is the ability to include a large number of exchange terms into spin wave calculations allowing the determination of magnetic excitation energies and thermodynamic properties. Prior literature has primarily included only a few of these exchange terms as surmised from the experiment or from fitting to only a dozen or so configurations.^25–27 Cluster expansion methods allow the calculation of a large number of exchange terms including three-body and higher interactions.

Electronic structure calculations were performed using the Vienna Ab initio Simulation Package (VASP)²⁸ within the framework of plane wave density functional theory in periodic boundary conditions, using the projector augmented wave (PAW)²⁹ method to represent the core electrons. Due to the large number of calculations required for the determination of ECIs, the exchange correlation is represented by the computationally inexpensive generalized gradient approximation (GGA+U) using the Perdew, Burke, and Ernzerhof (PBE)³⁰ variant. A U–J Hubbard correction value of 4 eV was used, consistent with prior computational work of Fe(III) oxides.^31–33 The ATAT VASP wrapper was modified to incorporate these Hubbard corrections into the CE calculations. The crystal structure is described using the $Cmc21$ space group. Unit cell parameters and fractional atomic coordinates were simultaneously varied in the energy relaxation. Calculations using the PBE0 functional³⁴ were performed to verify 0K properties independent of the Hubbard U parameter.

The plane wave basis energy cutoff for all structures was converged to energy differences of less than 1 meV per 50 eV increase. The maximum required energy cutoff, 900 eV, was used across all systems to ensure accurate energy comparisons. To avoid Pulay stress error, the energy cutoff was increased by 30% for all geometry optimizations. A static calculation at the original cutoff was then performed to obtain the total energy. Brillouin zone sampling was performed with gamma-centered k-points using the Monkhorst–Pack scheme³⁵ and again converged to energy differences within 1 meV for all unit cell crystal structures with a k-point mesh of $6×2×6$⁠. All supercell structures were sampled inversely proportional to their size in each cell dimension. Lattice constants were determined from the cell geometry optimization from VASP using an stopping criterion of $10−4eV$ between ionic steps.

Vibrational energy contributions were determined from frozen phonon mode calculations. VASP was used to calculate the dynamical matrix of the AF system at a series of volumes about the equilibrium volume. Thermodynamic properties were then calculated using Phonopy.³⁶ The constant pressure vibrational heat capacity was calculated with the quasi-harmonic approximation using the calculated values of the bulk modulus and thermal expansion coefficient.

The $J$s of Eq. (1) are given by the Effective Cluster Interactions (ECIs) between the Fe(III) spins. The results of the ATAT calculation are shown in Fig. 1. The estimated total energies of configurations using the ECIs were compared with ab initio calculations resulting in an excellent cross-validation score of $3.8×10−4$⁠, confirming the model’s predictive power. While multi-site interactions of Fe(III) clusters of size 3 and greater were calculated, only pairwise interactions were found to have contributions to the total energy greater than 1 meV. Table I details the properties of all pair interactions with energies greater than 1 meV and are visually represented in Fig. 2. Pairs 2 and 3 are both AF edge-sharing irons with weak interactions likely due to the effect of near $90°$ Fe–O–Fe bond angles on superexchange. The strongest interaction is AF pair 4 oriented along the c-axis. The Fe–O–Fe bond angle of $157°$⁠, the closest bond angle to $180°$ in the system, greatly contributes to the superexchange interaction. These findings are consistent with experimental work relating AF strength to Fe–O–Fe bond angle.¹⁰ Pairs 5 and 6 are the weakest interactions, occurring over large Fe–Fe distances within the crystal. Although weak, the multiplicity of the interactions would still contribute significantly to the total energy.

To determine the ground state magnetic order, configurations of increasing supercell size were exhaustively searched by ATAT, and the ECIs were used to quickly predict the energies of these systems. Discovered lower energy configurations were then confirmed with ab initio calculations.

The lowest energy structure found was a net zero magnetic moment AF system. Its structure departs from the prior computational work in that the crystallographic symmetry is broken by the inclusion of magnetic order, resulting in the smallest repeating unit of a $2×1×2$ supercell of the conventional cell. In the conventional cell, Fig. 3(a), the system is restricted to FM ordering along the a and c axes. Expansion to the supercell allows for the AF ordering along these same axes [Fig. 3(b)]. Of special significance is the allowed AF ordering along the c axis, where the Fe–O–Fe bonds are near $180°$⁠. Here, pair interactions are strongest, leading to a large decrease in the total energy. Exclusion of this strong interaction would likely have an impact on calculations of this system. Also of note is the predicted presence of AF layering, as described in the experimental literature.¹⁰ This layering occurs diagonally to the primary crystal axes in strips oriented between the a and c axes and perpendicular to the (010) surface as shown in Fig. 3(c).

Given the new predicted magneto-crystalline structure, it is of interest to calculate how the crystal properties may differ from those based upon the conventional cell. In Secs. III C–III F, we compare calculated values for the newly found lowest energy structure, referred to as AF, with those of the lowest energy structure found having the conventional cell symmetry, AF$′$⁠.

The density of states (DOS) is sensitive to the magnetic order of a material resulting in possible changes to peak shapes and bandgap sizes. Bandgaps previously calculated for lepidocrocite, on the basis of the conventional unit cell, have been used to determine the U–J Hubbard term used in GGA+U calculations by fitting to measured bandgaps of, for example, 2.06 eV.¹ However, such experimental measurements are typically performed at room temperature, and it is unclear what the magnetic disorder is at these temperatures and how it may influence bandgap measurements. Shown in Fig. 4, there is a significant difference between the calculated DOS of AF and AF$′$⁠. AF displays much sharper peaks and a significantly larger bandgap of 0.5 eV.

Calculated lattice parameters for the AF and AF$′$ structures are shown in Table II and are compared to experimental values. Lattice parameters were taken to be the cell dimensions of the crystallographic unit cells comprising the $2×1×2$ supercell. Notably, all cell parameters in the ground state structure are within 1% of experiment, an improvement relative to previous studies.^14,16,17 The errors become larger for the AF$′$ structure, indicating the importance of choosing the correct ground state magneto-crystalline cell. In particular, a large reduction in error occurs along the c axis entailing the near $180°$ Fe–O–Fe magnetic interaction. There is also an especially good match along the b axis where lepidocrocite is layered by hydrogen bonded sheets. This could possibly be due to two reasons. Deformation along the b axis may have an especially weak impact on the total energy due to relatively soft hydrogen bonding interactions connecting the layers. This would allow the interlayer magnetic interactions, though weak, to become impactful on equilibrium distances. Moreover, the c axis spacing may also be correlated with the interlayer distance. Although not practical to implement during the full ATAT analysis, the lattice parameters for the AF system were recalculated using the PBE0 functional and were found to match closely to those found using GGA+U.

The bulk modulus was calculated by fitting the Birch–Murnaghan equation of state to the energy–volume dependence and elastic moduli were calculated using the stress–strain method as implemented in VASP. To the best of our knowledge, experimentally measured elastic constants have not been reported. They are also not commonly reported in the computational literature. In Table III, we compare the calculated lattice constants for the AF system to these previous computational studies. For the bulk modulus, the Birch–Murnaghan EOS fit to the energy vs volume curve is shown in Fig. 5. The calculated bulk modulus is in good agreement with the result by Otte et al.³³ Although this work did not definitively state the magnetic order of their ground state, the agreement suggests the use of the correct ground state structure. In contrast, the bulk modulus is more than twice that reported by Guo and Barnard,¹⁶ where only configurations of crystallographic unit cell AF systems were considered. This is reflective of a broadening of the energy vs volume curve for the AF system when including the proper pair interactions. The elastic moduli were calculated directly from VASP and also show large differences from those computed for the crystallographic cell. A large difference between the c$33$ values is present as might be expected when comparing cells consisting of AF pairs vs FM pairs along this axis. More surprisingly, the largest difference is between c$22$ moduli, indicating more compressibility of the hydrogen bonded stacking layers in the AF system.

The cohesive enthalpies determined for both AF and AF$′$ are significantly different by 0.1 eV/(f.u.). However, although generally good agreement with experiment, the lower energy AF structure has a larger deviation. It should be noted that the experimental value is calculated from standard enthalpies of formation, as a $Δ$H for lepidocrocite extrapolated to 0 K is unavailable. Therefore, it is inconclusive from the energetics calculations which system is more representative of the real system.

The calculated ECIs correspond to Nearest Neighbor (NN) and Next Nearest Neighbor (NNN) exchange constants. The relative magnitudes of these calculated values are consistent with those from Goodenough–Kanamari rules for superexchange but are inconsistent with the low experimental Néel temperature. This contradiction may be explained by a closer examination of the unique symmetries present in lepidocrocite, which is not found in other iron oxides.

First, note that both the J3, J5, J6 couplings contain frustrated ferromagnetic bonds balanced by equal numbers of anti-ferromagnetic bonds leading to net zero contributions to the total energy. The remaining exchange terms, J3 and J4 each couple in ways that contribute to the total magnetic dimensionality of the system through the a and c lattice vectors respectively. Further J4 is nearly 10 times the magnitude of J3. In terms of $J/kb$⁠, they are 11 K and 146 K respectively, giving $J4≫J3$⁠. The consequence of this is that lepidocrocite likely behaves as a quasi-2D magnetic system at low temperature regimes but in the region approaching the experimental critical temperature, it is expected that lepidocrocite would behave as a quasi-1D system. Numerous studies have been performed on compounds possessing similar features^37–39 and the reduction of Néel temperature is specifically noted as a property of these systems,⁴⁰ but lepidocrocite does not appear to have been identified as such in prior literature. Hence, we present an analysis in which we use a 1D magnetic model and neglect all exchange terms except J4.

The thermodynamics of a Heisenberg quasi-1D spin 1/2 AFM chain can be modeled as excitations of spinons. The dispersion relation is given by Eq. (2),⁴¹ where J is the intra-chain coupling, k is the spinon wave vector, and a the lattice length along the chain,

The heat capacity can then be found using statistical mechanics. For the Heisenberg quasi-1D spin 1/2 AFM chain it has been derived by explicit ansatz (famously due to Bethe⁴²) and using conformal field theory.⁴³ For spinons, the heat capacity can be approximated by the more usable Eq. (3),⁴⁴ which closely reproduces the results from the Bethe ansatz solution.⁴⁵ Here, the spinons are considered to have fermionic excitations, therefore, $f$ is given as the Fermi–Dirac distribution with the zero chemical potential and $ε$ is given by Eq. (2),

Equation (3) was evaluated numerically using the trapezoidal rule. The maximum of the resulting curve (plotted in Fig. 6) is the predicted Néel temperature. From the heat capacity curve calculated for J4, a critical temperature of 69 K was found in excellent agreement with experimental value of 68 K from the heat capacity studies by Majzlan, where a slight anomalous bump was noted in the region of the Néel temperature.⁴⁶ The total heat capacity calculations are compared to experimental heat capacity curves in Fig. 6. The calculated contributions to the total heat capacity are comprised of the vibrational and magnetic portions. Electronic contributions were omitted due to the large bandgap, precluding occupancy of conduction bands at low temperatures. The addition of the calculated magnetic heat capacity to the total curve matches well with the experimental data in the short temperature region leading up to the critical point. The difference between the magnetic heat capacities is emphasized in Fig. 7, where the calculated vibrational heat capacity has been subtracted out from the experimental data. At low temperatures, the quasi-1D model greatly overestimates the heat capacity and this region is likely better fit by a proportional $T2$ model. Heat capacity contributions have previously been fitted to the low temperature regions,⁴⁷ but relied upon inclusion of a linear scaled electronic contribution attributed to oxygen vacancies as well as a 3D AFM $T3$ magnetic term in order to fit the experimental data. A complex magnetic landscape at these low temperatures however may make a physically accurate fitting difficult in practice. At and above the critical temperature, the heat capacity begins to diverge greatly from the calculated results. This is likely due to the inability to properly model excitations in both the magnetic and the vibrational terms.

We have used density functional theory coupled with cluster expansion methods to determine parameters for all significant cluster interactions within a lepidocrocite crystal having high spin $Fe↑xFe↓1−x$OOH stoichiometry aligned parallel to the c axis. These parameters were then used to determine a ground state AF magnetic order, which was found to be consistent with experiment, but requiring a larger supercell of the conventional cell, a finding that therefore supersedes prior computational work. Lattice properties of this new ground state structure were calculated and found to yield improvements in the lattice constants, and significant differences in the elastic moduli, cohesive enthalpy, and electronic DOS.

Magnetic thermodynamics properties were determined through the direct calculation of the magnetic heat capacity assuming a Heisenberg quasi-1D spin 1/2 AFM chain as the dominant model near disordering temperatures. A Néel temperature of 69 K was determined, consistent with experimental results. The magnetic dimensionality view of this system provides a reasonable explanation resolving the conflict between the expected strong superexchange and the low temperature onset of magnetic disorder. This view would also explain the large disparity with respect to the magnetic critical temperatures of other iron oxides, of which none contain $180°$ Fe–O–Fe linear chain morphologies.

See the supplementary material for modifications made to the ATAT code, complete listing of elastic constants and BM fit data, tabulated experimental total heat capacity data, calculated vibrational heat capacities, and calculated spin wave dispersion data.

This work was supported by IDREAM (Interfacial Dynamics in Radioactive Environments and Materials), an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences (BES). Computational work was performed, in part, at the Center for Institutional Research Computing at Washington State University and using supercomputing resources at the Environmental Molecular Sciences Laboratory (EMSL), a national scientific user facility at Pacific Northwest National Laboratory (PNNL). EMSL is sponsored by the Department of Energy’s Office of Biological and Environmental Research. PNNL is a multi-program national laboratory operated for DOE by Battelle Memorial Institute under Contract No. DE-AC05-76RL0-1830.

The data that support the findings of this study are available in the supplementary material.