The role of disorder in NaO$_2$ and its implications for Na-O$_2$ batteries

TL;DR

This study uses DFT with Hubbard U correction to explore structural and magnetic disorder in NaO₂, revealing complex magnetic phases and their implications for Na-O₂ battery performance.

Contribution

It introduces a 3-body Ising Model combined with Monte Carlo simulations to accurately describe disorder in NaO₂, advancing understanding of its properties in batteries.

Findings

NaO₂ exhibits significant rotational and magnetic disorder.

Ferromagnetic phases near room temperature form alternating orthogonal O₂ dimers.

Bulk NaO₂ structures are insulating with some ferromagnetic structures showing a small electronic gap.

Abstract

We present a DFT study utilizing the Hubbard U correction to probe structural and magnetic disorder in $\mathrm{NaO_{2}}$, primary discharge product of Na-O$_2$ batteries. We show that $\mathrm{NaO_{2}}$ exhibits a large degree of rotational and magnetic disorder; a 3-body Ising Model is necessary to capture the subtle interplay of this disorder. MC simulations demonstrate that energetically favorable, FM phases near room temperature consist of alternating bands of orthogonally-oriented $\mathrm{O_{2}}$ dimers. We find that bulk structures are insulating, with a subset of FM structures showing a moderate gap ($<2$ eV) in one spin channel.

Tables

Table 1. Table 1: Calculaed bandgap values (eV) for various ferromagnetic and antiferromagnetic configurations of NaO 2 , for U=3 eV using HSE06 functional.

Structure	Spin 1(+)	Spin 2(-)
AAAA(PPPP)	6.65	1.26
AAAB(PPPP)	7.43	3.69
AABB(PPPP)	7.03	3.96
AABC(PPPP)	7.3	3.92
ABCD(PPPP)	5.91	1.32
AAAA(PPNN)	3.82	3.82
AAAB(PPNN)	3.91	3.91
AABB(PPNN)	4.13	4.13
AABB(PNPN)	4.21	4.21
AABC(PPNN)	4.39	4.13
ABCD(PPNN)	3.94	3.94

Equations

$$\mathrm{N}{\mathrm{a}}^{+}+{\mathrm{O}}_2+{\mathrm{e}}^{-}\rightleftharpoons {\mathrm{Na}{\mathrm{O}}_2}_{(\mathrm{s})}$$

$$\mathrm{\Delta }{\mathrm{H}}_{{\mathrm{F}}_{\mathrm{Na}{\mathrm{O}}_2}}={\mathrm{E}}_{\mathrm{Na}{\mathrm{O}}_2}^{\mathrm{DFT}}-{\mathrm{E}}_{\mathrm{Na}}^{\mathrm{Ref}}-{\mathrm{E}}_{{\mathrm{O}}_2}^{\mathrm{Ref}}+\mathrm{\Delta pV}$$

$${\mathrm{E}}_{{\mathrm{O}}_2}^{\mathrm{Ref}}=2{\mathrm{E}}_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{DFT}}-2{\mathrm{E}}_{{\mathrm{H}}_2}^{\mathrm{DFT}}-\mathrm{\Delta }{\mathrm{H}}_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{Exp}}$$

$${\mathrm{E}}_{\mathrm{N}}{\mathrm{a}}^{\mathrm{Ref}}={\mathrm{E}}_{\mathrm{N}}\mathrm{aC}{\mathrm{l}}^{\mathrm{DFT}}-\frac{1}{2}{\mathrm{E}}_{{\mathrm{Cl}}_2}^{\mathrm{DFT}}-\mathrm{\Delta }{\mathrm{H}}_{\mathrm{N}}\mathrm{aC}{\mathrm{l}}^{\mathrm{Exp}}$$

$$\mathrm{E}=-\sum _{⟨\mathrm{ik}⟩}{\mathrm{j}}_{2_{\mathrm{i},\mathrm{k}}}{\sigma }_{\mathrm{i}}{\sigma }_{\mathrm{k}}-\sum _{⟨\mathrm{ikl}⟩}{\mathrm{j}}_{3_{\mathrm{i},\mathrm{k},\mathrm{l}}}{\sigma }_{\mathrm{i}}{\sigma }_{\mathrm{k}}{\sigma }_{\mathrm{l}}$$

$$\Omega =\prod _{j=1}^M\left(\genfrac{}{}{0ex}{}{{\sum }_{i=j}^M{N}_i}{N_j}\right)$$

Full text

Electrification of road transport recently led to increased interest in high-energy-density rechargeable batteries, in particular, metal-air batteries. Luntz and McCloskey (2014); Sapunkov et al. (2015) Li-O2 battery, one of the most promising, owing to its high energy density, suffers from electrolyteBurke et al. (2015); Abraham and Jiang (1996); Abraham and Alamgir (1990); Bryantsev et al. (2011, 2013); Bryantsev and Faglioni (2012); Walker et al. (2013); Elia et al. (2014); Khetan et al. (2014); Burke et al. (2015) and electrodeOttakam Thotiyl et al. (2012); Thotiyl et al. (2013); McCloskey et al. (2012a, b) instability, as well as limited discharge capacity.Viswanathan et al. (2011); Aetukuri et al. (2015); Chen et al. (2013); Bruce et al. (2012); Johnson et al. (2014) Hartmann et al. recently demonstrated a rechargeable Na-O2 battery, with sodium superoxide (NaO2) as the primary observed discharge product, which showed superior cycle life to state-of-the-art Li-O2 batteries.Hartmann et al. (2013a, b) Subsequently, there have been numerous theoretical and experimental efforts to understand the fundamental mechanisms that govern electrochemistry in Na-O2 batteries.Kang et al. (2014); Yang and Siegel (2015); Hartmann et al. (2013b); Das et al. (2014); Hartmann et al. (2013a); Ellis and Nazar (2012); Kundu et al. (2015); Kwak et al. (2015); Ortiz-Vitoriano et al. (2015)

Accurate description of the Na-O2 battery reaction mechanismsKim et al. (2013); Kang et al. (2014) requires detailed understanding of electronic structure throughout the phase space of sodium oxides. The roles of nucleation,Krishnamurthy et al. (2016); Ortiz-Vitoriano et al. (2015); Hartmann et al. (2015) nanoscale stabilization,Kang et al. (2014) and surface energetics of various sodium-oxygen compounds have been examined through a combination of density functional theory (DFT) calculations and electrochemical measurements. This has led to an improved understanding of Na-O2 battery reactions, which, as proposed, constitute a combination of surface and solution mechanisms of the NaO2 discharge product.Adelhelm et al. (2015); Sapunkov et al. (2015) In certain cases, it should be noted, Na2O2 has been reported as the discharge product,Jian et al. (2014); Araujo et al. (2015) and selectivity between NaO2 and Na2O2 is not yet fully characterized.

Alkali superoxides are known to be highly disordered materials, both in geometric structure and magnetic ordering.Smith et al. (1966); Ren and Wu (2013) Presently, there is very limited understanding of geometric and magnetic disorder in NaO2, due to challenges in experimentally synthesizing NaO2. It is crucial to map out disorder in room-temperature NaO2, considering its importance in determining electronic structure, surface energetics, and growth properties, relevant for Na-O2 batteries. In this work, we aim to fill the gaps in the current understanding of the electronic structure of NaO2. We perform DFT calculations incorporating the Hubbard U correction and use that to understand the geometric and magnetic disorder in NaO2. We employ an Ising model accounting for magnetic and geometric degrees of freedom and find that both 2-body and 3-body nearest-neighbor interactions, which take into account both the geometric and magnetic disorder of the structure, are critical to accurately describe bulk NaO2. The Ising model is used within a metropolis monte carlo framework to characterize the complete phase space as a function of temperature. In order to determine the effect of disorder on the electronic properties of NaO2, we employ hybrid density functional theory calculations and find that the calculated bandgap is strongly affected by disorder. All investigated structures exhibit a direct bandgap around 4 eV, with the exception of a small subset of ferromagnetic structures, whose direct bandgap is below 2 eV in one spin channel.

The discharge process of the Na-O2 battery involves the $\mathrm{Na^{+}}$ ion coupled electron transfer reaction with O2 at the cathode, which produces NaO2, the primary reported discharge product.Hartmann et al. (2013a, b)

[TABLE]

A few studies do report the formation of Na2O2,Jian et al. (2014) which is the thermodynamically stable structure according to the Na-O phase diagram at standard temperature and pressure.Kang et al. (2014) The preferential formation of NaO2 as a discharge product in room-temperature Na-O2 batteries, over the thermodynamically stabler Na2O2, still remains a puzzle. The stable phase of bulk NaO2 itself is also temperature-dependent.Wriedt (1987); Carter and Templeton (1953); Templeton and Dauben (1950); Kang et al. (2014) Below 196 K, NaO2 takes the $\mathrm{Pnnm}$ space group. Between 196 and 223 K, NaO2 takes the Pa3̄ space group. Above 223 K, NaO2 takes the Fm3̄m space group, the primary relevant structure for room-temperature Na-O2 batteries.

To explore the reaction mechanism of the Na-O2 battery, we need to understand the bulk composition of NaO2 and accurately describe its electronic structure and dependent properties. The oxygen dimer in NaO2 behaves similar to the superoxide anion, O${}_{2}^{-}$.Partridge et al. (1992) Within the molecular orbital picture of ground-state (triplet) O2, the 2 highest occupied orbitals are the $\mathrm{\pi_{2p}^{*}}$ antibonding orbitals, with one electron in each orbital, following Hund’s Rule of Maximum Multiplicity.Boyd (1984) The highest occupied molecular orbital of O${}_{2}^{-}$ remains $\mathrm{\pi_{2p}^{*}}$, now occupied by 3 electrons, as illustrated in the appendix. This means that the O${}_{2}^{-}$ dimer has $1$ unpaired electron, making O${}_{2}^{-}$ magentic. The magnetic ordering of these spins is one of the key contributions to the configurational disorder of bulk NaO2.

We considered bulk NaO2 in the Fm3̄m space group structure. In this structure, sodium occupies the body centers and corners of the cubic cell, and oxygen dimers occupy face centers and edge centers of the cell, as illustrated in the appendix. We introduce the following naming scheme for these structures to describe both their geometric and magnetic arrangement. The first 4 letters, which can be A, B, C, or D, refer to the 4 possible geometric orientations of oxygen dimers. The next 4 letters, which can be P or N, refer to the net positive or negative magnetic moment of the corresponding dimers. Thus, a structure designated AAAA-PPPP has all dimers mutually parallel to each other (AAAA), and all dimers of identical, positive magnetization (PPPP), providing a ferromagnetic structure. Alternatively, a structure designated ABCD has all dimers mutually orthogonal to each other, while a structure designated PNPN or PPNN is antiferromagnetic. Formation enthalpies were calculated for 19 possible structures to explore the phase space of bulk NaO2.

Self-consistent DFT calculations were performed using the Projector Augmented Wave Method as implemented in GPAW,Enkovaara et al. (2010) with the Revised PBE (RPBE) exchange correlation functional.Hammer et al. (1999) To correct electron localization in NaO2, we incorporated the Hubbard U, applied on oxygen $2p$ states in NaO2.Hubbard (1963); Anisimov et al. (1991); Himmetoglu et al. (2014); García-Mota et al. (2012); Meredig et al. (2010) Calculations were run with a real-space grid of 0.18 Å, and 6$\times$6$\times$6 k-point sampling, following the Monkhorst-Pack scheme.Monkhorst and Pack (1976) Fermi-Dirac smearing of 0.01 eV was used to facilitate convergence and broyden-type mixing of electron densities was used in the calculation.Johnson (1988); Srivastava (1984)

Two schemes were used to calculate formation enthalpy of bulk NaO2. The formation enthalpy of NaO2 is given by:

[TABLE]

The pressure-volume work term, $\mathrm{\Delta pV}$, can be disregarded, as it is typically about 5 orders of magnitude smaller than internal energy contributions in formation enthalpy calculations.Aydinol et al. (1997); Obrovac and Chevrier (2014) Internal energies of Na, O2, and NaO2 were calculated using DFT simulations.

It is well known that molecular oxygen is poorly described in DFT.Droghetti et al. (2008) The energy of O2 can be more accurately calculated using water as the reference molecule.Rossmeisl et al. (2005) In the first scheme used, reference energy of O2 was computed using the DFT-calculated internal energies of H2O and gaseous H2, as well as the experimental formation enthalpy of H2O:

[TABLE]

The Hubbard U correction was applied on the oxygen atoms in both H2O and NaO2.

The second scheme was used to calculate reference energy of bulk Na. It was demonstrated in prior work that formation enthalpies of alkali oxides, peroxides and superoxides are best described when the oxidation state of the metal in the reference compound is matched to the oxidation state of the metal in the compound under investigation.Christensen et al. (2015) Following this scheme, the Na reference energy was calculated using the calculated internal energies of NaCl and gaseous Cl2, as well as the experimental formation enthalpy of NaCl:

[TABLE]

Currently, there is no universally accepted method for selecting the optimal Hubbard U to fit calculations to experimental data. In this study, we follow the scheme that compares the calculated formation enthalpy of a substance under investigation to the experimentally measured formation enthalpy.Wang et al. (2006) The experimentally measured formation enthalpy of NaO2 is $\mathrm{{\Delta}_{f}H^{0}_{solid}=-2.7eV}$.Chase (1998)

We find that the formation enthalpies of all considered bulk configurations of NaO2 lie within 0.3-0.6 eV/NaO2, at all investigated values of the Hubbard U, as shown in Fig. 1(a). The fully-organized, ferromagnetic (AAAA-PPPP) and fully disorganized, antiferromagnetic (ABCD-PPNN) structures were found to be least energetically stable, while the moderately disorganized structures with intermediate magnetic order, such as AABB-PPPN, were most stable. Many of these moderately disorganized structures were energetically degenerate at all examined values of the Hubbard U. We observe that U = 3 eV is matches the experimental formation energy and will be used for further analysis. It is also worth highlighting that the conclusions regarding disorder in bulk NaO2 remain consistent for all examined values of the Hubbard U.

Given the high degeneracy observed among examined structures of NaO2, it becomes crucial to examine the role of thermal disorder and its role in determining electronic properties of NaO2. The computation of disorder in larger supercells directly through DFT is computationally intractable. We employ the approach of determining a lattice HamiltonianHamer and Barber (1981) that can be used to describe the energetics of bulk NaO2.

In order to map out energetic interactions between the magnetic and geometric degrees of freedom, we utilize a modified Ising Model for the lattice Hamiltonian. As implemented, the model consists of a lattice of $N$ sites $i$, whose filling is described by occupation terms, $\sigma_{i}$. System energy contributions due to particle interactions in neighboring sites are captured by the 2-body and 3-body interaction terms, $j_{2_{i,k}}$ and $j_{3_{i,k,l}}$. The total energy of the $N$-site lattice is calculated as:

[TABLE]

Formation enthalpies of $\mathrm{NaO_{2}}$ calculated by DFT simulations were used to fit Ising Model coefficients for our system. Formation enthalpy data at each value of the Hubbard U was used to derive a set of corresponding interaction coefficients, though a least-squares regression fit. The sole differences among the structures lay in the relative geometric and magnetic arrangement of the oxygen dimers, and these differences are reflected in the nearest-neighbor interaction terms $j$. A total of 4 different types of $j_{2_{i,k}}$ terms are identified, along with 6 different types of $j_{3_{i,k,l}}$ terms. Our naming scheme relies on comparison of dimers, e.g. $j_{2AAPP}$, corresponds to the interaction between 2 dimers of parallel geometry and identical magnetization, while $j_{2ABPN}$, corresponds to the interaction between 2 dimers of orthogonal geometry and opposite magnetization. Likewise, $j_{3AAAPPP}$, corresponds to the interaction among 3 dimers of mutually parallel orientation and mutually identical magnetic moment, while $j_{3ABCPPN}$ corresponds to the interaction among 3 oxygen dimers of mutually orthogonal orientation, where 2 dimers share an identical magnetization that differs from that of the third dimer. The plot of fitted interaction coefficients is shown in Fig. 1(b). The plot demonstrates that among the 2-body interaction terms, $j_{2AAPN}$ was most energetically stabilizing, while among the 3-body interaction terms, $j_{3AABPPP}$ was most energetically stabilizing.

Our DFT calculations suggested that the most thermodynamically stable configurations correspond to bulk structures of intermediate disorder, both in geometric and magnetic arrangement (e.g. AABB-PPPN). Subsequent Ising Model calculations supported this observation, with $j_{2AAPP}$ the most-stabilizing 2-body interaction term, and $j_{3AABPPP}$ the most-stabilizing 3-body interaction term. It should be noted that using a purely 2-body Ising Model of NaO2 finds that the most-stabilizing nearest-neighbor interaction term was $j_{2ABPN}$, implying that the fully-disorganized antiferromagnetic structure (e.g. ABCD-PNPN) would be the predominant energetically stable structure. This is not consistent with our DFT calculations, and therefore shows the need to incorporate 3-body interaction terms into the model.

Interaction coefficients $j$ derived from the Ising Model were used in Metropolis-Hastings Monte Carlo (MHMC) simulations, to characterize larger bulk cells with higher degree of disorder, at finite temperature. Supercell structures were constructed as N$\times$N$\times$N arrays of cubic unit cells of $\mathrm{Na_{4}O_{8}}$, the structures studied in DFT simulations. Simulations were run for supercells of size 4 to 7, to study convergence of bulk properties as the overall simulated system size increases. Periodic boundary conditions in all directions were implemented and the systems were studied using the Markov Chain Monte Carlo Method,Gilks (2005); Hastings (1970) modified with the Metropolis-Hastings Algorithm.Chib and Greenberg (1995) Interaction coefficients corresponding to a Hubbard U correction of 3 eV were used with the MHMC simulations, to most closely match the experimental formation enthalpy of NaO2. The results would remain qualitatively consistent for other values of U, and there would be a minor change in the observed phase transition temperatures.

At the initialization of each MHMC simulation, a fully-organized supercell bulk structure was set up, with all dimers identical (e.g. AAAA-PPPP). Each individual trial step consisted of switching a randomly selected $\mathrm{O_{2}}$ dimer in the bulk structure to one of the other available configurations; orientational or magnetic change. The effect of this switch on system energy was calculated. The standard Metropolis Algorithm was used to decide whether to accept the trial step or not. At the end of every step, system formation enthalpy and entropy were recorded for further analysis.

Configurational entropy was calculated as the logarithm of the number of configurations with the same proportions of oxygen dimers of different types. If there are $M$ distinct dimer types available, the total number of possible configurations, $\Omega$, of the full structure can be calculated using the complete multinomial coefficient:

[TABLE]

where $i$, $j$ refer to the available configurations of dimers and $N_{i}$, $N_{j}$ refer to the number of dimers of a particular configuration present in the supercell. The configurational entropy is then simply calculated by Boltzmann’s Entropy Formula: $S_{Conf}=k_{B}\log{\Omega}$.

To explore high-temperature, high-energy phases of bulk $\mathrm{NaO_{2}}$, each system simulation was initially raised to a temperature of 1252 K. Once the energy of the structure was stabilized at the initial high temperature, the system was annealedVanderbilt and Louie (1984) in temperature steps of 0.25 K, down to 2 K. At each 50 K increment, formation enthalpy and configurational entropy data was collected for analysis. This scheme worked robustly and the thermodynamic quantities converged for all supercell sizes investigated.

MHMC simulations were used to determine effect of temperature on the structural disorder and in turn the free energy of bulk $\mathrm{NaO_{2}}$, as shown in Fig. 1(c) and Fig. 1(d). When the temperature is high, most of the structures could have become accessible, as the entropic term dominates. Thus, the structure at high temperature was predominantly fully-disordered, both in geometry and in magnetic moment as shown in Fig. 1(e). Around 650 K, the structure underwent a phase change, wherein the bulk $\mathrm{NaO_{2}}$ formed large, parallel planes of alternating geometry of $\mathrm{O_{2}}$ dimers, as depicted in Fig. 1(f). In this phase, the magnetic moment was still largely disordered, within and across the monolithic-geometry planes. Around 350 K, the bulk structure underwent a second phase change, where the magnetic moment of all dimers was locked in phase together and the geometry maintained the alternating planar structure. The degree of disorder in the system could strongly affect the electronic properties of NaO2, so we explored this further using hybrid density functional theory calculations.

Non-self-consistent hybrid density functional theory calculating employing the HSE06Heyd et al. (2003) functional was used to determine direct bandgap for all phases of NaO2, as an indicator of electrical conductivity. HSE06 has been shown to accurately capture bandgaps for various semiconducting and insulating materials with mean absolute error (MAE) varying from 0.26$-$0.4 eV for different classes.Heyd et al. (2005); Tran and Blaha (2009); Henderson et al. (2011); Chan and Ceder (2010) Bandgap values were calculated across all k-points and the minimum was reported for each phase of NaO2 for Hubbard U = 3-7 eV. We find that there is an increase in bandgap with increasing U as shown in Table II in S.I. The widening of the bandgap with increasing U is due to enhanced localization of electrons in the $\pi^{*}_{2p}$ oxygen bands in NaO2. We also observe that the bandgap in both spin channels is identical for antiferromagnetic structures, and markedly different for ferromagnetic structures. This behavior is observed as the additional electrons in O${}_{2}^{-}$ ions only occupy one spin channel. In the AFM structures, electrons are equally distributed in the spin channels, and the gap is between pairs of $\pi^{*}_{2p}$ states. In the FM structures, the bandgap in the occupied spin channel is between $\pi^{*}_{2p}$ and $\sigma^{*}_{2p}$ state and is thus higher than AFM; while for the other channel, it is between deeper $\pi^{*}_{2p}$ states and is thus lower than that in AFM. For mixed cases (e.g. PPPN), we expect the bandgap to be bounded by the AFM and FM cases, due to electron occupation of the $\pi^{*}_{2p}$ bands. Finally, we observe that the most symmetric structures, AAAA-PPPP and ABCD-PPPP.

In conclusion, we demonstrate the need and effectiveness of a 2-body and 3-body nearest-neighbor Ising Model in describing energetics of bulk NaO2 across a large temperature range. Our study demonstrates that at low temperatures, close to room temperature and below, the predominant phase of NaO2 is ferromagnetic, with alternating planes of oxygen dimers in consistent geometric orientations. The system is expected to exhibit some degree of magnetic disorder even at room temperature. Our study on the electronic properties shows that NaO2 is a wide-bandgap insulator, with a bandgap around 4 eV and we expect it to have poor electrical conductivity at room-temperature. In the context of Na-O2 batteries, this implies that growth of the discharge product is most likely occurring due to the solution mechanism pathway involving a chemical dissolution of NaO2 into Na+ and O${}_{2}^{-}$, similar to what is seen in Li-O2 batteries. However, our analysis shows that preferential nucleation of certain magnetic phases could be possible through appropriate electrode choice. Further, our approach is expected to be important for other metal-oxygen batteries, such as K-O2 batteries.

Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund and National Science Foundation CAREER award CBET-1554273 for partial support of this research.