Combined study of phase transitions in the P2-type Na X Ni 1/3 Mn 2/3 O 2 cathode material: experimental, ab-initio and multiphase-ﬁ eld results

The research of new electrode materials such as sodium intercalation compounds is key to meet the challenges of future demands of sustainable energy storage. For these batteries, the intercalation behavior on the micro-scale is governed by a complex interplay of chemical, electrical and mechanical forces strongly in ﬂ uencing the overall cell performance. The multiphase-ﬁ eld method is a suitable tool to study these multi-physics and bridge the scale from ab-initio methods to the cell level. In this work, we follow a combined approach of experiments, density functional theory (DFT) calculations and multiphase-ﬁ eld simulations to predict thermodynamic and kinetic properties for the P2-type Na X Ni 1/ 3 Mn 2/3 O 2 sodium-ion cathode material. Experimentally, we obtain the thermodynamic potential and diffusion coef ﬁ cients at various sodium contents using electrochemical techniques and discuss limitations of the experimentally applied methods. DFT is used to identify stable phases by calculating an energy hull curve. Then, the in ﬂ uence of long-range dispersion interactions and the exchange-correlation functional on the voltage curve is investigated by comparison with experimental results. Finally, multiphase-ﬁ eld simulations are performed based on inputs from experiments and DFT. The ﬁ tting of phase-speci ﬁ c chemical free energies from DFT calculations and experimental data is discussed. Our results highlight the thermodynamic

electronic or structural phase transitions 2 and pronounced lattice distortion upon insertion and extraction of ions.Depending on the applied voltage, the occurring phase transitions can be partially irreversible 3,4 , which crucially limits cell stability and performance and thus needs to be investigated in more detail.
Layered oxides based on the structure Na X TMO 2 are promising candidates for sodium-ion intercalation batteries and have received considerable attention throughout the past couple of years.However, first works date back to the emergence of rechargeable batteries 2 as both lithium and sodium intercalation compounds were initially investigated.O3 and P2-type layered oxides are the most common structural polymorphs 5 where O denotes the octahedral and P the prismatic coordination of sodium-ions, respectively.The number corresponds to the amount of transition metal layers per unit cell 6 .In O3-type layered oxides, the sodium-ion octahedra share edges with the transition metal (TM) octahedra.P2-type layered oxides on the other hand exhibit two different prismatic sodium-ion environments, either sharing faces (denoted as Na f ) or sharing edges (Na e ) with the transition metal octahedra.Transition metals can be (Cu, Cr, Fe, Co, Ni, Mn, Ti, V) or combinations thereof [2][3][4]7,8 . Som of the most promising cathode materials for sodium-ion batteries are based on a combination of Ni and Mn due to their comparatively high capacity, attractive working potential, high ionic Na + mobility and low toxicity 4,9,10 .
Lu and Dahn 7 first investigated the P2-Na X Ni 1/3 Mn 2/3 O 2 layered oxide and found high reversibility which matches the previous findings for Na X CoO 2
Little effort has been addressed to modeling and simulation of this material, which is surprising given that it is one of the most often experimentally investigated cathode materials for future sodium-ion batteries.In the field of lithium-ion material research the amount of simulative works is continuously increasing, covering all time-and length-scales from the atomic to the cell level.For investigation of material properties on the atomic level, ab-initio methods such as density functional theory (DFT) have been proven to be important in finding new promising structures 27 and identifying the electro-chemical as well as mechanical properties of known compounds [28][29][30][31] .First-principles calculations are well-established to compute phase diagrams 32 , ordering effects 15,33 , diffusion barriers 15 , and open circuit voltages (OCV) 34 .However, the quality of computed material properties depends strongly on the level of exchange and correlation functionals E xc used.The most frequently employed E xc functional 35 is the Perdew-Burke-Ernzerhof functional 36 (PBE) based on the generalized gradient approximation 37 (GGA), which fails in the description of transition metal oxides with localized electrons.It has been shown that the Strongly Constrained and Appropriately Normed semilocal (SCAN) 38 functional outperforms PBE for the case of the lithium layered oxides such as LiNiO 2 , and LiMnO 2 39,40   .Furthermore, the introduction of a Hubbard U correction provides further improvement and is an appropriate method for predicting the electronic properties of oxides containing strongly correlated 3d electrons for both PBE and SCAN functionals 41,42 .While the exact U parameter depends on chemical factors such as the oxidation state of the transition metal 43 and adapts to the charge carrier concentrations.The convex hull for the entire range of charge carrier concentrations is determined using a datadriven approach, taking into account a single fixed U parameter across the entire range 44,45 .In addition, in layered compounds it is often crucial to take dispersion corrections into account as the van der Waals (vdW) interaction significantly contributes to the bonding between the layers 46 .However, dispersion-corrected PBE functionals typically overestimate this interaction.In contrast, the SCAN functional shows an improved description of intermediate-range VdW interactions and, in combination with the revised Vydrov-van Voorhis functional (SCAN+rVV10) 47,48 , additionally accounts for the long-range dispersion interactions.
DFT is limited to a certain amount of atoms and mostly concerned with bulk properties of crystalline structures.Free surfaces and interfaces are only accessible with high computational effort and under very strong assumptions while phase transformations and the rich interplay of surface reactions, diffusive transport and mechanical deformations are completely out of range 28 .To the best of our knowledge, little effort has been made to bridge the gap between the atomistic and the continuum level 29 , which is a crucial factor for the computational screening of cathode materials.Some fundamental aspects are well established, e.g., the relation of energies above hull with electro-chemical potentials and the resulting theoretical voltage profile 28 .The phase-field method is well suited to study intercalation dynamics in single crystals and agglomerates as has been demonstrated in studies on the commercialized LiFePO 4 (LFP) cathode material for lithiumion batteries 49,50 .This framework is especially useful for the simulation of materials undergoing phase transformations as moving surfaces are inherently difficult to track otherwise.The phase-field method has a strong thermodynamic foundation as the evolution of phases, concentration and stresses is based on minimization of the Gibbs free energy.Multi-physics problems coupling chemical, electrical and mechanical driving forces can be investigated 51,52 .The transfer of material data from DFT to Cahn-Hilliardtype phase-field methods has been discussed by Hörmann et al. 29 .
This work is dedicated to the investigation of non-equilibrium processes occurring during charge and discharge of Na X Ni 1/3 Mn 2/3 O 2 , which undergoes multiple phase transitions and charge orderings.For this purpose we experimentally derive the thermodynamically relevant potential and solid diffusion coefficients for Na X Ni 1/3 Mn 2/3 O 2 at various sodium contents in the "Electrochemical testing" section.In the "DFT" section, we use DFT to calculate an energy hull curve for stable phases and investigate the influence of long-range dispersion interactions and the exchange-correlation functional on the voltage curve with respect to the experimental data.Results from experiments and DFT are used to perform multiphase-field simulations in the "Gibbs free energy fitting for multiphase-field simulations" and "Kinetics of phase transformations" sections.First, the fitting of phasespecific chemical free energies from DFT calculations is discussed with respect to experimental results.Following, the inclusion of a Butler-Volmertype boundary condition allows to study the kinetics of the system as a competition of surface reaction, bulk diffusion and elastic deformation upon phase transformations.To the best of our knowledge, this is the first report which successfully simulates multiple phase transitions for an intercalation cathode material for sodium-ion batteries.The combination of experimental data and ab-initio results provides solid ground for highly valuable phase-field simulations results.

Electrochemical testing
The chemical composition of the prepared cathode active material was analyzed with ICP-OES, from which a chemical formula according to Na X Ni Y Mn 1−Y O 2 is calculated.The resulting molecular formula is Na 0.69 Ni 0.34 Mn 0.66 O 2 , which is in good agreement with our targeted stoichiometry of Na 2/3 Ni 1/3 Mn 2/3 O 2 .Powder X-ray diffraction (XRD) was performed to characterize the crystal structure of the obtained cathode active material.The obtained diffraction pattern is presented in Fig. 1.The diffraction pattern can be fully described with space group P6 3 /mmc and matches well with ICDD 00-054-0894.The two additional reflections at 27.2°and 28.3°are associated with "large zigzag" Na-vacancy ordering 15 .No further reflections are observed, indicating the phase-pure nature of the material.A close inspection of the diffraction pattern reveals a slight broadening of the (10l) reflections, hinting toward stacking faults 53,54 .A Rietveld refinement accounting for stacking faults was performed following the model presented in our previous publication 54 .A satisfying fit is obtained, the refined lattice parameters are a = 2.890 Å and c = 11.147Å and the stacking fault probability was found to be as low as 4.2%.The lattice parameters are in good accordance with literature 11,54,55 .The refinement of the crystallite size from diffraction data is hampered due to the presence of these stacking faults.However, no influence of these stacking faults on the electrochemical performance is expected 53,54 .
SEM was performed to characterize the morphology and particle architecture of the cathode active material powder as presented in Fig. 1.The powder consists of spherical secondary particles in the range of 5-20 μm in diameter.The crystals offer distinct facets mimicking the hexagonal layered crystal structure of the material.Based on SEM images the crystal size can be roughly estimated to ~1 μm in (100) and (010) direction and 0.5 μm in (001) direction.
To electrochemically obtain the potential of P2-Na X Ni 1/3 Mn 2/3 O 2 at various sodium stoichiometries X, galvanostatic intermittent titration technique (GITT) is performed in coin cells with sodium metal serving as counter electrode.In GITT measurements, short defined current pulses are applied followed by a relaxation period to obtain the potential of the material in equilibrium state (Fig. 2a).Based on the assumption, that every measured electron corresponds to the (de)intercalation of one sodium-ion (no side reactions), the sodium stoichiometry can be calculated for each step based on the measured charge and the mass of the active material in the electrode.By this means, the equilibrium potential of the active material can be experimentally explored over a broad range of composition in one single experiment.The transition metal oxide host structure is assumed to be stable (no breaking of TM-O bonds at room temperature), leading to metastable compositions such as O2-Na 0 Ni 1/3 Mn 2/3 O 2 , which are not accessible by high-temperature synthesis 16 .
A GITT measurement is represented in Fig. 2b for P2-Na X Ni 1/3 Mn 2/3 O 2 using a two-electrode coin cell with sodium metal as counter electrode.Low loadings and specific currents (17.3 mA g −1 ) are applied to minimize any contribution of the counter electrode or diffusion through the porous working electrode composite.Therefore, the measured cell voltage closely reflects the potential of the P2-Na X Ni 1/3 Mn 2/3 O 2 cathode active material.Both, the potential after relaxation (OCV) and during the current pulse (C/10) are presented in Fig. 2c.For the OCV curve (blue), multiple distinct voltage steps at sodium contents of X ≈ 2/3, X ≈ 1/2 and X ≈ 1/3 are apparent.In literature, these voltage steps are associated with Na-vacancy ordering 15 .At X = 2/3 the Na f sites form a "large zigzag" pattern, at X = 1/2 the ordering can be described as alternating rows of Na f and two rows of Na e 15 .At X = 1/3, single layer row ordering of Na e is reported 15 .At sodium contents X < 1/3 the voltage plateau indicates a two-phase reaction in accordance with diffraction experiments 7,17 .Small voltage plateaus are observed before and after the above mentioned Na-vacancy orderings as indicated by sharp peaks in the dQ/dV plot (Inset Fig. 2c).These voltage plateaus are located at 3.15 V, 3.31 V, 3.62 V, 3.66 V and 4.15 V corresponding to sodium stoichiometries of 0.66, 0.52, 0.48, 0.36, 0.32, respectively.
In Fig. 2c, a difference in cell voltage between the C/10 curve (red) and the OCV curve (blue) is observed due to overpotentials (mainly) on the working electrode, caused by phenomena such as charge transfer resistance or solid diffusion.At medium sodium contents of 0.35 ≤ X ≤ 0.6, the difference between the equilibrated voltage and the measured voltage under applied current is nearly negligible, indicating fast kinetics.At high sodium contents of X > 0.7, the difference between the measured voltages under current and in equilibrium is significant (>100 mV), indicating slow kinetics of the active material.After the first charge, the measured relaxed potentials during discharge differ from the ones obtained in the preceding charge, which is most likely caused by side reactions at high voltage causing a mismatch of the calculated sodium stoichiometries.
Furthermore, based on GITT measurements solid diffusion coefficients can be calculated as demonstrated by Weppner and Huggins 22 : where D app is the calculated apparent diffusion coefficient, τ is the duration of the current pulse, m B is the mass of the active material, M B is the atomic weight of P2-Na X Ni 1/3 Mn 2/3 O 2 , V M is the molar volume of P2-Na X Ni 1/3 Mn 2/3 O 2 , S is the electrochemically active area, ΔE S is the potential change due to the change of the stoichiometry (relaxed state), ΔE t is the voltage change during the current pulse neglecting the iR drop and L is the diffusion length.Note, that several assumptions are made by Weppner and Huggins during the derivation of this equation to calculate solid diffusion coefficients from GITT measurements based on the Nernst equation and Fick's second law.Using Nernst equation to couple the concentration of active ions at the electro-chemically active surface with the measured potential indirectly assumes continuous change of the potential with the surface concentration of the active ion (solid solution like behavior of the active material).As boundary conditions, semi-infinite solid diffusion (diffusion length L) is assumed.The derivative of equilibrated potential with the surface concentration is approximated to be linear, which only holds true for (infinitesimal) small changes in stoichiometry.The same approximation is made for the derivative of the transient potential with time, which limits the duration of the current pulse to (infinitesimal) small time frames.Additionally, the molar volume is considered to be constant.
To calculate apparent diffusion coefficients based on the equation by Weppner and Huggins, we assume a crystal morphology in the form of cuboids with 1 μm in dimension in a-b plane and 0.5 μm in c-direction.The selected crystal shape and size is in accordance with SEM images.To estimate the electro-chemically active surface area, Kr-Physisorption measurements have been performed to derive the specific surface area of the active material using the BET method.The obtained surface area is 0.51 m 2 g −1 .In a layered structure, only the prismatic surface is available for sodium-ion (de)intercalation 56 .Based on the chosen crystal morphology, we would expect half of the crystal surface to be prismatic.Therefore, the electrochemical active surface area is estimated to be half of the obtained Kr-BET area.Please note that the selected crystal morphology would result in higher surface area than measured via Kr-BET, if solely isolated crystals would be present within the sample.However, as the crystals form secondary particles due to agglomeration (compare SEM), the surface area of the single crystals is partially reduced.The molar volume is calculated based on our refinement of XRD pattern and kept constant for all compositions.The iR drop of the cell is calculated based on the high frequency resistance measured with PEIS and the applied current.Based on these assumptions, we have derived apparent diffusion coefficients from two cells with different τ (150 s and 30 s, respectively) as presented in Fig. 2d.Within a medium sodium content (0.35 ≤ X ≤ 0.6), the apparent diffusion coefficient is as high as 2 × 10 −10 to 10 −13 cm 2 s −1 .At distinct sodium contents of X ≈ 1/3, X ≈ 1/2 and X ≈ 2/3, high apparent diffusion coefficients are observed.This observation matches very well with reports on Na-vacancy orderings at these stoichiometries 9,15 .For narrow sodium compositions close to X ≈ 1/3, X ≈ 1/2 and X ≈ 2/3, the material follows a solid solution like behavior.The corresponding voltage plateaus adjacent to this distinct sodium stoichiometries (light blue background in Fig. 2d) each represent a order-disorder phase transition between the respective Na-vacancy ordered phase and the corresponding solid solution phase.Within the course of these order-disorder phase transitions significant lower diffusion coefficients are calculated.At high (X > 0.7) and low (X < 0.3) sodium content, the apparent diffusion coefficient is significantly lower (10 −15 to 10 −12 cm 2 s −1 ).These observations match very well with the observed polarization during the current pulses, which is low at medium sodium content and higher at low and high sodium content.As described above, the electrochemical (de)sodiation at low sodium content (1/3 < X < 0, light red background in Fig. 2d) follows a two-phase reaction 7 .Kinetics in first-order phase transformations might be either controlled by the migration of the phase boundary or the diffusion across this phase boundary 57 .Additionally, for open systems such as battery intercalation materials kinetics of two-phase reactions might also be controlled by solid diffusion or by reaction kinetics 58,59 .Since the derivation of the equation assumes a solid solution regime and Fick's law of diffusion, the calculated values in this region most likely do not represent the real diffusion coefficients of the compound 56,60,61 .Furthermore, in phase transformation materials the assumption that all particles do simultaneously participate when the electrode is being (dis)charged does not necessarily hold true 56 .Therefore, the active surface area, active mass and eventually also the diffusion length might differ over the course of a two-phase reaction, further hindering the evaluation of diffusion coefficients using the equation by Weppner and Huggins.In theory, solid diffusion coefficients for both end members of the two-phase reaction could be derived from measurements in their respective single-phase boundary solubility region 60 .Simulations by Han et al. demonstrate, that GITT and PITT can deliver reasonably accurate diffusion coefficients in the single-phase region of phase change materials 61 .
By applying the obtained diffusion coefficients to τ ≪ L 2 /D app , one must check, if the assumed boundary conditions (semi-infinite diffusion) is fulfilled 22 .In our case (τ = 30 s, L = 0.5 μm), this assumption holds true for all diffusion coefficients D app < 8.33 × 10 −11 cm 2 s −1 , which covers nearly all of our measurement points.For the highest obtained diffusion coefficients the boundary condition of semi-infinite diffusion cannot be guaranteed.Therefore, the measured values might even underestimate the real diffusion in the active material.Even shorter current pulses are necessary to fulfill the assumed boundary conditions.Extensions of the GITT method by Weppner and Huggins have been published for the calculation of diffusion coefficients of porous electrodes 62 and phase-transformation electrodes 60 .Unfortunately, both extensions require either additional assumptions or measured values, eventually resulting in high uncertainty for the calculated diffusion coefficients due to error propagation.
Overall, for P2-Na X Ni 1/3 Mn 2/3 O 2 we find high diffusion coefficients in the range of 2 × 10 −10 to 10 −12 cm 2 s −1 for 1/3 < x < 2/3 and smaller diffusion coefficients in the range of 10 −13 to 10 −14 cm 2 s −1 for x > 2/3, which is in accordance with reported trends in the literature 9,12,25 .For fast ionconductors such as P2-Na X Ni 1/3 Mn 2/3 O 2 , care must be taken to respect the assumption of semi-infinite diffusion when calculating diffusion coefficients.Within two-phase reactions (light blue and light red background in Fig. 2d) the assumption of Fick's law does not hold true and diffusion coefficients are most likely underestimated 60 .By pointing out these restrictions and limitations, we hope to arise awareness in the scientific community for possible limitations of experimentally derived diffusion coefficients in the literature.

DFT
In the P2 phase of Na X Ni 1/3 Mn 2/3 O 2 , the transition metals form a honeycomb structure with the Ni atoms in the center of the honeycomb surrounded by the Mn atoms as shown in Fig. 3a.As described above, the Na ions occupy prismatic sites that are either edge-sharing (Na e ) or face-sharing (Na f ) with the transition metal octahedra.The honeycomb arrangement of the transition metals leads to the two distinct Na f sites, one sharing a face with two Mn octahedra (Na f MnÀMn ), the other with one Mn octahedron and one Ni octahedron (Na f NiÀMn ), while all Na e are equivalent in structure (see Fig. 3b).
The presence of different Na sites in the structures gives rise to distinct Na arrangements for the Na intercalation in which the relative superstructures were previously identified for Na contents X = 1/3, 1/2, and 2/3 15 .Furthermore, at low Na content, the P2 structure undergoes a phase transformation by sliding the transition metal layers against each other, forming the O2 phase 7,9,15 .To derive the convex hull for the Na intercalation, it is necessary to investigate the possible Na arrangements as well as to consider both structure types at low Na contents.The convex hull was generated based on both the conventional PBE + U + D3 level of theory and the more accurate SCAN and SCAN+rVV10 functionals.In this way, 20 Na configurations were generated based on the previously identified superstructures 15 for both Na contents X = 1/2 and 2/3.For the Na range from X = 0 to X = 1/3, we considered 10 Na configurations for both the P2 and the O2 structure type at Na contents X ∈ [0, 0.04, 0.08, 0.16, 0.25, 1/3] applying the same procedure as for the recently studied related compound P2-Na X Mn 3/4 Ni 1/4 O 2

26
. This procedure is based on the fact that at low Na contents, the Na arrangements depend primarily on both the Na site energies and the in-plane electrostatic repulsion between Na ions.Therefore, the Na site energies for low Na contents corresponding to one or two Na atoms per supercell were first determined and found to increase in energy in the order Na e < Na f NiÀMn < Na f MnÀMn .Then the structure was successively filled to account for the different Na contents by considering the in-plane electrostatic Na-Na repulsion by maximizing the Na-Na distance.In the next step, the five lowest energy configurations at each Na content X = 0, 1/3, 1/2, and 2/3 were chosen as starting geometry for the calculations using the SCAN and SCAN+rVV10 functionals in combination with a data-driven Hubbard-type U correction 45 .The convex hulls obtained with the different E xc functionals, shown in Fig. 4, predict stable Na orderings at the Na contents X = 1/3, 1/2, and 2/3.Notably, the compositional formation energies for the orderings at X = 1/3 and 1/2 are significantly lower in the case of the SCAN+U calculations than for the PBE+U+D3 and SCAN+rVV10+U calculations.At a Na content of X = 1/3, all considered functionals predict the O2 phase and the P2 phase to be energetically degenerate within the accuracy of DFT, while experimentally the P2 phase is observed at this Na content 7,9,15 .
The OCVs for the Na concentration ranges X = 0-1/3, X = 1/3-1/2, and X = 1/2-2/3 are obtained from the respective lowest energy configurations on the convex hull and are shown in

GPa: ð3Þ
The full tensor is included in the multiphase-field simulations in the following section.

Gibbs free energy fitting for multiphase-field simulations
The primary particles of Na 2/3 Ni 1/3 Mn 2/3 O 2 typically exhibit a hexagonal shape and platelet-like morphology with the c-axis oriented along the thin dimension.In this study, we consider a single crystalline platelet with a simplified cuboid shape with dimensions according to our experimental results (1 μm × 1 μm × 0.5 μm).The fluxes resulting from the intercalation reaction are applied on all surfaces but bulk diffusion is limited to the layers, i.e., D 001 = 0. Hence, sodium can adsorb on the large crystal facets in the (001)-direction but the bulk of the crystal is only accessible through insertion at (100)-and (010)-facets followed by bulk diffusion.
Table 1 | OCV calculated based on the PBE+U+D3, SCAN+U, and SCAN+rVV10+U levels of theory for Na content ranges   As described above, Na 2/3 Ni 1/3 Mn 2/3 O 2 crystallizes in the hexagonal P63/mmc space group 15 .The primitive hexagonal unit cell contains two formula units (f.u.) as it spans two layers of sodium-ions/transition metals and its volume is V ¼ ffiffi ffi 3 p =2 a 2 c ¼ 80:4Å 3 .Consequently, we compute the reference concentration as c ref = 2/(N A V P2-2/3 ) = 41,319 mol m −3 which is used to convert from eV to J m −3 or, alternatively, de-dimensionalize simulation input parameters.The lattice constants of the intercalation host are both phase-and concentration-dependent.Normal eigenstrains along the primary crystal axes have been calculated with reference to the P2-Na 2/3 Ni 1/3 Mn 2/3 O 2 (P2-2/3) phase, e.g., ϵ 0 a ¼ ða O2 À a P2À2=3 Þ=a P2À2=3 as listed in Supplementary Table 1.Strong anisotropy becomes obvious from the calculated values.The stiffness tensors computed from DFT are directly included as an input for multiphase-field simulations.The following simulation studies are based on the assumption that all deformations are elastic and occurring interfaces are coherent with an interfacial energy of γ αβ ≈ 0.1 J m −2 63 .The standard exchange current density j 00 ¼ Fk 0 c max which is associated with the energetic barrier for intercalation, is assumed to be ≈0.1 A m −2 64 .
The chemical energetic landscape is approximated by phase-wise fitting functions where we include the O2-phase and three ordered states of P2 (1/3: single-row, 1/2: double-row, 2/3: large-zigzag).In this work, two types of ansatz functions are used to fit the phase-specific chemical energy, namely a quadratic and a logarithmic expression: The diffusion potential is defined as the derivative with respect to the site filling fraction x α , which yields: for the two cases given in Eqs. ( 4) and ( 5).We employ a fit purely based on quadratic functions (4) with parameters given in the Supplementary Table 2 and shown in Fig. 5a.These parameters have been determined from the hull energy results of DFT calculations based on the SCAN+U +rVV10 functional.To account for entropic effects at room temperature (T = 300 K), logarithmic terms have been added to the formation energy data points such that: The parameter x max in Eq. ( 5) is chosen to be 0.68 to limit the composition in the relevant voltage range to x ∈ [0, 2/3].The equilibrium in two-phase regions is given by the tangent construction depicted by dashed lines, which results in the voltage plateaus according to Eq. ( 13).If we replace the energies of the first (O2) and last phase (P2-2/3) with logarithmic functions (Eq.( 5)) and the parameters in Supplementary Table 2, the corresponding voltage curve exhibits more realistic slopes for X → 0 and X → 2/3 as shown in Fig. 5b.
The resulting voltage in thermodynamic equilibrium matches the experimental OCV (relaxed points from GITT measurement) remarkably well.The largest deviation is observed at X → 2/3 where the predicted plateau is higher compared with experiments.The two function fits differ for X → 0 and X → 2/3 where the logarithmic fits approach −∞ and ∞ while the derivation of quadratic energies yields a linear relation with respect to x α .As a consequence, the logarithmic functions are naturally bounded to the interval x α ∈ [0, 2/3] while for quadratic functions the concentration values can be negative.
Although the chemical energies f α chem in Fig. 5 are individual convex functions, the overall multiphase-field energy landscape (compare Eq. ( 15)) is continuous in the multi-dimensional space of volume fractions ϕ α and site filling fraction x.The obstacle potential f ob represents an energetic barrier between phases which penalizes phase coexistence and, thus, leads to a miscibility gap.In total, the homogenous free energy landscape is given by f chem + f ob .This is an important difference compared to the Cahn-Hilliard model where the homogenous energy is formulated as one continuous function of the composition x.The separation into chemical and obstacle energy density allows for an effective re-scaling of the diffuse interfacial width and an accurate description of interfacial energy for every αβ-phase pair 65 .

Kinetics of phase transformations
The chemical energy fits can be used within the framework of the multiphase-field method to study the dynamic behavior of Na X Ni 1/3 Mn 2/3 O 2 during charge and discharge.To study the interplay of surface reaction and bulk diffusion, a simulation study with varying C-rates is conducted.Diffusivities are assumed to be phase-wise constant and are approximated from GITT results in the "Electrochemical testing" section as D O2 = 8 × 10 −14 cm 2 s −1 and D P2 = 2 × 10 −10 cm 2 s −1 , respectively.All simulations are based on the function fit including logarithmic terms due to better agreement with the experimental results in the equilibrium case.
For comparison, extreme current rate tests without formation cycles have been performed at C/10, C/5, 1C and 5C (Fig. 6a-d).At slow rates (C/10 and C/5) the potential profiles run as expected with large initial charge capacities of 172 mAh g −1 for C/10 and 156 mAh g −1 for C/5, respectively.Under faster current rates (1C and 5C), the performance drastically decreases to capacities of 113 mAh g −1 and 86 mAh g −1 .The electrochemical curves in Fig. 6c, d show a drastic increase in overpotential and hysteresis.The poor rate performance, especially at higher potentials can result from the low ionic diffusion coefficient and poor electronic conductivity 66 .
The simulation results in Fig. 6e exhibit some general features caused by the interplay of surface reaction and bulk diffusion.With increasing Crate, the deviation between the applied voltage and equilibrium voltage increases which can be rationalized by two factors.A higher charging rate corresponds to a larger flux of sodium from the active material to the electrolyte.As the overpotential is the driving force for intercalation, the applied voltage must be higher to extract more sodium (see Eq. ( 18)).Secondly, the higher flux of ions at the surface leads to a higher concentration gradient from the surface toward the bulk of the particle.Therefore, the difference between the surface concentration x surf and the average composition X in Na X Ni 1/3 Mn 2/3 O 2 becomes larger.The applied voltage in Fig. 7e depends on x surf and is plotted over X.As a result, the voltage hysteresis increases, plateaus are sloped and steps in the voltage curve are less pronounced in the diffusion-controlled regime 67 .
The phase transitions are exemplarily shown for 1C in the inset pictures in Fig. 6e.The transition between P2 ordered states proceeds through a shrinking core mechanism at 1C and 5C.At low C-rates, the shrinking core is replaced by a wave-like front moving through the crystal in accordance with simulations for LCO 68 .The first order phase transition between the O2 and P2-phase preferentially proceeds perpendicular to the layers which is a result of the strong lattice contraction in the c-direction (see Fig. 7).At high C-rates, sodium extraction becomes diffusion limited which causes a transition from the layer-by-layer mechanism to a shrinking core as shown in Fig. 7.At 5C, the O2-phase forms at all active facets which hinders the extraction of the remaining sodium due to the comparatively slow diffusivity D O2 .This corresponds to the large capacity loss at 5C in Fig. 6 and, additionally, causes severe stresses (bending of particle in Fig. 7).
In all cases, phase transitions are initiated at corners and edges of the crystal where the curvature is high.In simulations where the elastic deformation is neglected, the continuous removal of sodium ions always leads to a ring of the O2 phase forming at the active facets of the crystal (i.e., shrinking core mechanism).This hinders the extraction of ions from the bulk of the particle and results in a blocking effect caused by the slow diffusivity in the O2-phase.As a result, the accessible capacity at all C-rates is significantly underestimated if elastic deformation is not accounted for (see Supplementary Fig. 2).The experimentally observed capacities can only be achieved through layer-by-layer filling, which is preferential for the minimization of stored elastic energy.
The lattice mismatch between phases leads to stress hotspots at locations with high concentration gradients.In brittle materials, the maximum tensile stress, which is the largest eigenvalue of the local stress tensor, can be used as an estimate for possible fracture.In our simulations, we find high tensile stresses close to the top and bottom facets during charge and in the center of the particle during discharge.The temporal evolution of stresses is shown in the Supplementary Videos.The observed stress distributions are consistent with the observation that cracks can form inside the particle and reach the surface through continued cycling 69 .The simulated voltage curves feature strong resemblance with the experimentally obtained ones from which we conclude that the model covers the dominant effects and can be used for predictive material simulation.The differences between the curves in Fig. 6a-e give some hints for further model refinement in future works.Very low overpotentials are observed in the range of X ∈ [0.33, 0.5], which suggests that the charge transfer coefficient k 0 might be phase-specific.The deviation at high filling fractions X > 0.6 results from the energy fit in the "Gibbs free energy fitting for multiphase-field simulations" section.Furthermore, there is an asymmetry between charge and discharge in the high voltage plateau ≈4.16 V, which becomes apparent through the high overpotential in the discharge even at low rates (Fig. 6a, b).To the best of our knowledge the reason is not yet fully understood and could be clarified by future detailed simulation studies.

Discussion
In this work, we combine experimental investigations with simulations on two different length scales to shed more light on the phase transition mechanisms involved in the cycling of Na X Ni 1/3 Mn 2/3 O 2 .The results highlight the consistency of thermodynamic modeling with experimental results.Hull energies for stable phases obtained through DFT simulations can be used to fit Gibbs free energies that are used as an input for multiphase field simulations.The accuracy of the OCV was shown to depend on the exchange-correlation functional, with SCAN+U and SCAN +rVV10+U providing voltages in good agreement with experimentally measured ones.However, the relative phase stability, which is of particular importance in a multi-scale approach, depends strongly on the inclusion of long-range dispersion interactions.The resulting equilibrium potential exhibits a close match with the OCV obtained from GITT measurements.
Furthermore, the dynamic micro-battery model reflects the overall electrode behavior at various C-rates surprisingly well.The presented multiphase-field model covers several phase transitions and successfully predicts the expected overpotential at technically relevant charging rates.The simulation results highlight that the inclusion of elastic deformation is crucial to capture the phase transformation mechanism of the O2-P2 transition.Due to the large lattice mismatch in the c-direction, the phase transformation proceeds perpendicularly to the diffusion direction at low to moderate charging rates.Simulations, which neglect the elastic energy contribution falsely predict a strong blocking effect upon charge resulting from the slow diffusion in the O2-phase.At high charging rates, the capacity resulting from the O2-P2 phase transformation cannot be utilized which is a material-intrinsic limitation.The anisotropy of layered materials restricts diffusion to the ab-plane while deformation mainly occurs along the c-direction.Therefore, the phase transformation preferentially proceeds perpendicular to the diffusion direction.In the diffusion limited regime, the O2-phase forms at all active facets and reduces the mobility of sodium due to the reduced interlayer distance.These effects might be mitigated if an OP4-phase is formed instead 26 .In future investigations, the established simulation framework could be used to study the multi-particle behavior during phase transformations to identify the rate limiting step and improve the understanding of GITT measurements.

Equilibrium thermodynamics
Let us first consider a sodium-ion battery under open-circuit conditions.It is important to note that the driving force for the operation of a rocking-chair battery such as a sodium-ion battery is the energy gain upon the transfer of a Na atom from the anode to the cathode 70,71 .The open-circuit voltage V OC can directly be derived from this energy gain ΔG via: where F is the Faraday constant and z is the charge of the ion in the electrolyte, which is z = 1 for alkali metals such as Li, Na, of K.If ΔG is expressed in eV, then this energy directly corresponds to the open circuit potential.Assuming that the enthalpic contributions dominate the free energy difference ΔG, V OC can be easily derived from total energy calculations based on density function theory (DFT) 72 by taking the difference in the binding energies of Na in the anode and the cathode.Note furthermore, that this energy difference and hence also the open-circuit voltage is a function of the concentration of Na in cathode and anode and hence of the state of charge of the battery.Equivalently, V OC can also be expressed as the difference between the chemical potentials μ A and μ C of Na in the anode and the cathode, respectively 70 : with e being the elemental charge.The chemical potential of Na can be decomposed into the electrochemical potentials of the Na + ion and the electron according to μ Na ¼ μNa þ þ μe À .In a battery under open-circuit conditions, anode and cathode are connected by an ion conducting and electronically isolating electrolyte.Hence in equilibrium the electrochemical potentials of the Na + ion has to be constant throughout the whole cell: which will be achieved through the formation of electric double layers (EDL) at both anode and cathode 73 .
In this work, we consider the intercalation of sodium into a layered oxide Na X TMO 2 cycled against a sodium metal electrode.At constant temperature and pressure, the chemical potential of metallic sodium is constant (μ A Na ¼ μ É;A Na ) while the chemical potential in the cathode depends on the site filling fraction of vacancies x: In the continuum approximation, the diffusion potential μ(x) is related to the partial derivative of the Gibbs free energy density with respect to the site filling fraction μ(x) = ∂f/∂x 74,75 , i.e., μ(x) is a continuous function of the local average of filled sites in the phase-field method.Under these assumptions, Eq. ( 10) can be expressed as: where the reference cell voltage V É cell ¼ ðμ É;A Na À μ É;C Na Þ=e is independent of x and thus constant.In the non-equilibrium case, a gradient in chemical potential can arise in the cathode from the surface into the bulk of active particles due to diffusion limitation and ∇μ is the respective driving force for diffusion.This thermodynamic background facilitates the comparison of experimental results that are close to equilibrium (i.e., low C-rate) with total energies calculated by DFT in the "DFT" section.Furthermore, both can be used as an input to formulate the Gibbs free energy density f for phase-field simulations which is discussed for Na X Ni 1/3 Mn 2/3 O 2 in the "Gibbs free energy fitting for multiphase-field simulations" section.In the following, we denote the local average of filled sites with x while the average composition in the cathode is given by X = (∫xdV)/V.

Electrochemical testing
A spherical phase pure P2-type Na X Ni 1/3 Mn 2/3 O 2 material was synthesized following a scalable two-step process.First, a Ni 1/3 Mn 2/3 (OH) 2 precursor was prepared by coprecipitation in a continuously stirred tank reactor (CSTR, 1L).Appropriate ratios of metal nitrate solution (Mn(NO 3 ) 2 × 4 H 2 O, Ni(NO 3 ) 2 × 6 H 2 O, both Carl Roth), sodium hydroxide solution (NaOH, Carl Roth) and ammonia solution (NH 4 OH, Carl Roth) are mixed under vigorous stirring.The obtained Ni 1/3 Mn 2/3 (OH) 2 precipitate was thoroughly washed, filtered and dried.In the second step, this precursor was mixed with appropriate amounts of NaOH (Carl Roth) solution (wet impregnation).The obtained mixture was dried, homogenized and calcined for 10 h at 900 °C in pure oxygen in a box furnace (Carbolite Gero).The resulting powder was directly transferred into a Büchi glass oven, where the powder was kept under dynamic vacuum at 200 °C overnight.Subsequently, the dried powder was transferred into an Ar filled glovebox (MBraun, O 2 < 0.1 ppm, H 2 O < 0.1 ppm) to avoid any surface impurities resulting from reaction with ambient moisture and CO 2 76,77   .Electrode preparation and cell assembly were performed inside the same Ar filled glovebox to avoid any reaction with the ambient and therefore minimize side reactions during electrochemical characterization.
Chemical and physical characterization of the Na X Ni 1/3 Mn 2/3 O 2 powder was performed outside the glovebox using fresh samples of the same batch for every characterization technique.The chemical composition was analyzed with ICP-OES (Spectro Arcos SOP) in a diluted aqua regia solution.Crystal structure was characterized using powder XRD in Bragg-Brentano geometry with a Cu Kα X-ray tube on a Bruker D8Advance.Rietveld refinement of the obtained diffraction pattern was performed using FAULTS software 78 .Throughout this manuscript, Vesta Software 79 is used to prepare depictions of crystal structures.Particle architecture and morphology were investigated using SEM at 5 kV acceleration voltage on a Zeiss Leo 1530 VP equipped with an Everhart-Thornley-SE detector.The specific surface area of the Na X Ni 1/3 Mn 2/3 O 2 powder was derived from Kr-Physisorption at 77 K on a 3Flex 3500 from Micromeritics based on the BET method.
Electrode sheets were prepared by mixing polyvinylidene difluoride binder (PVDF, Solvay Solef P5130), conductive carbon (Timcal Super-P-Li) and active material in the weight ratio 8:8:84 dispersed in appropriate amount of N-methyl-2-pyrrolidone (NMP, anhydrous, Sigma-Aldrich).The resulting slurry was thoroughly mixed for ~3 h and then cast on Al-foil using the doctor blade technique.Electrodes with 12 mm diameter are punched and dried overnight in a Büchi glass oven at 130 °C under dynamic vacuum.The active material loading of these electrodes is ~3-4 mg cm −2 .Coin cells (CR2032, Hohsen) are prepared using these electrodes, two layers of glass fiber separator (GFA, Whatman) soaked with 150 μl 1M NaPF 6 in PC + 5% FEC and a 16 mm diameter sodium metal counter electrode (Acros Organics).
For GITT measurements, a VMP3 potentiostat (Biologic) is used.Current pulses at C/10 rate (1C = 173 mA g −1 ) are applied for 150 s each, where 4.3 V and 1.5 V were used as upper and lower cut-off, respectively, to avoid pronounced side reactions with the electrolyte.GITT measurements with shorter current pulses of 30 s at the same rate are performed using the same equipment to evaluate diffusion coefficients for sodium stoichiometries 0.3 < X < 0.69.Complete relaxation was assumed, when the voltage change was below 0.1 mV h −1 .The ohmic drop was calculated based on impedance measurements at high frequency (200 kHz) and the applied current.Typically, the ohmic resistance of our cells was in the order of 4 Ω.Electrochemical symmetric rate tests have been performed for Na X Ni 1/3 Mn 2/3 O 2 at selected rates C/10, C/5, 1C and 5C.T-cells are used with a 12 mm sized working electrode, two glass fiber (GFA, Whatman) separators, Na-metal as counter and reference electrode and 300 μl 1M NaPF6 in PC + 5% FEC as electrolyte.

DFT calculations
Periodic first-principles calculations based on density functional theory (DFT) 80,81 are performed on O2 and P2 phases of Na x Mn 2/3 Ni 1/3 O 2 containing various Na arrangements.The PBE 36 , SCAN 38 , and SCAN+rVV10 47,48 exchange and correlation functionals are used with the projector augmented-wave (PAW) 82 method to compare the results of optimized geometries, calculated open circuit voltages (OCV), as well as band gaps.In addition, the influence of dispersion effects was captured using the Grimme D3 correction 46 in combination with the PBE functional as implemented in Vienna Ab initio Simulation Package (VASP) [83][84][85] .The strongly correlated 3d-electrons of the Ni and Mn atoms have been treated by the Hubbard-type U corrections 86 .Since the U values depend on the choice of the E xc we used the following data-driven Hubbard-type U corrections of U Mn = 3.9 eV and U Ni = 6.2 eV 44,87 for the PBE functional and U Mn = 1.99 eV and U Ni = 0.41 eV 45 for the SCAN functional.A plane wave cutoff of 520 eV and an electronic convergence criterion of 1 × 10 −6 were chosen, while the Brillouin zone was sampled using a 3 × 3 × 3 k-point mesh.A large supercell containing 24 formula units and ferromagnetic ordering was used, and the structures were relaxed without any restriction until all forces converged below 0.01 eV Å −1 .The OCV of the cathode material is derived from the respective convex hull for the Na intercalation range constructed from the compositional formation energies E f x .The decomposition reaction into the phases O2-Na 0 Ni 1/3 Mn 2/3 O 2 and P2-Na 2/3 Ni 1/3 Mn 2/3 O 2 is chosen as reference and E f x is calculated from the energies of the reference phases E 0 and E 2/3 and the energies E x at the respective Na content x by: Once the convex hull is constructed by calculating the formation energy per formula unit, the OCV for each stable configuration is obtained following Eq.( 9).The chemical potential may be approximately expressed by the DFT total energy by neglecting entropic contributions.In addition, the elastic stiffness tensors were calculated for O2-Na 0 Ni 1/3 Mn 2/3 O 2 and P2-Na 1/3 Ni 1/3 Mn 2/3 O 2 by applying an energy convergence criterion of 1 × 10 −7 eV and a force convergence criterion of 1 × 10 −3 eV Å −1 .The width of the displacement of each ion was set to 0.015 to minimize the influence of higher-order terms on the elastic constants.

Multiphase-field model
In the general picture, the multiphase-field method 88,89 covers transitions between various physical states, e.g., solid-liquid, order-disorder or different crystal structures.These phase transitions are described based on the evolution of their respective volume fractions captured by the time-and spacedependent components ϕ α of the phase-field vector ϕ ¼ fϕ 1 ; :::; ϕ α ; :::; ϕ N g T .Based on our previous works 52,63 , we formulate an energy functional: consisting of interfacial f int ¼ f grad þ f ob and bulk energy density contributions f bulk = f chem + f elast .A detailed model description is given in the Supplementary Methods.Phase transitions in battery materials (especially lithium iron phosphate) are often modeled by the Cahn-Hilliard (CH) model 49,90 .Practical challenges are, however, that for a battery material with multiple phase transitions such as Na X Ni 1/3 Mn 2/3 O 2 and other sodium intercalation compounds, it can be quite challenging to find an appropriate function fit for the chemical free energy.Secondly, the CH model is a microscopic model 91 in the sense that the diffuse interface represents the natural scale of diffuseness between two crystalline phases (i.e., interfacial width of few nanometers).This small length scale has to be resolved for the interfacial energy to be accurate which leads to high computational cost.The multiphase-field approach based on Allen-Cahn equations 88,89 on the other hand allows for effective re-scaling of the interfacial width while preserving the physical interfacial energy γ αβ which enables the simulation of microstructures in the micrometer range.A comparison of the CH model and the multiphase-field approach for lithium iron phosphate can be found in ref. 63.
The relevant kinetic processes of intercalation include ionic and electronic transport, interfacial ion transfer and, possibly, phase transformations 67 .Furthermore, elastic stresses due to a change in lattice parameters strongly couple with diffusion and phase transformation kinetics 49 .The evolution of phases is governed by Allen-Cahn equations 51,74 which are coupled to the evolution of the diffusion potential for mass conservation 92 .As mechanical relaxation proceeds on a much faster timescale, the static momentum balance is solved in every timestep.Purely elastic deformations are considered while explicit modeling of plasticity, fracture or dislocations is beyond the scope of this work.Detailed expressions of the evolution equations and respective boundary conditions are given in the Supplementary Methods.
Particularly important in the context of this work, is the formulation of sodium intercalation dynamics.While section "Equilibrium thermodynamics" holds true for both DFT and the phase-field method, the consideration of non-equilibrium dynamics is beyond the scope of DFT.The insertion reaction: is driven by the overpotential η, which is given by the difference of electrochemical potentials of Na + ions at the electrolyte-electrode interface η ¼ ðμ Furthermore, assuming good electronic wiring in the cathode and fast electron conduction compared to the time scales of ion intercalation and diffusion, μC e À is constant across the cathode.The voltage difference between electrodes measured by a potentiometer corresponds to the difference of electro-chemical potentials of the electrons V ¼ μA e À À μC e À =e 93 .Employing Eqs. ( 12) and ( 13), we get: where V is the applied voltage.μ surf is the diffusion potential of sodium at the surface of the active material and will be spatially resolved in simulations.
The previously discussed model assumptions lead to the notion of the micro-battery model sketched in Fig. 8.The reaction is modeled by the Butler-Volmer equation 93,94 assuming symmetry of the forward and backward reaction (α = 0.5; more details in Supplementary Methods): The assumption of galvanostatic (dis-)charge is fulfilled by restricting the global current: into the active material.n denotes the inward surface normal and j ⋅ n = j N .
For any given C-rate, we first solve Eq. ( 19) to compute the globally applied voltage V.This voltage is then inserted into Eq.( 18) together with the local surface site filling fraction x surf and the diffusion potential μ surf , to compute the surface reaction of sodium R surf = j N /F.

Fig. 2 |
Fig. 2 | GITT measurement results.a Schematic diagram of the voltage change as a function of time during a galvanostatic pulse and of the applied current.b Typical voltage response (discharge-charge-discharge) obtained during a GITT experiment of P2-Na X Ni 1/3 Mn 2/3 O 2 cathode.c Comparison of the difference in cell voltage between C/10 (red) and OCV (blue).The corresponding dQ/dV plot obtained from OCV curve is presented as inset.d Apparent diffusion coefficients derived from GITT measurements of two cells with τ = 150 s (red, black) and τ = 30 s (orange).

Fig. 3 |
Fig. 3 | DFT crystal structure of Na 2/3 Ni 1/3 Mn 2/3 O 2 .a The honeycomb configuration of the transition metals in the P2 phase of Na 2/3 Ni 1/3 Mn 2/3 O 2 .b The DFT-optimized supercell with two different Na sites: Na in a face with the Mn and Ni octahedra Na f NiÀMn , the top layer, and Na in an edge-sharing site Na e , the lower layer.The Ni atoms are colored in gray, Mn atoms in magenta, Na atoms in yellow, and O atoms in red.

Fig. 4 |
Fig. 4 | Formation energies (F.E.) per formula unit with respect to O2-Na 0 Ni 1/3 Mn 2/3 O 2 (blue markers) and P2-Na 2/3 Ni 1/3 Mn 2/3 O 2 (red markers).a PBE+U+D3, b SCAN+U and SCAN+rVV10+U levels of theory.The convex hull resulting from the PBE+U+D3 calculations is shown as a gray dashed line, also in (b) for comparison.The SCAN+U convex hull is shown as a dotted and the SCAN +rVV10+U hull as a solid black line.

Fig. 5 |
Fig. 5 | Gibbs free energy fittings (top row) and resulting equilibrium voltages (bottom row).a Parabolic fitting and b logarithmic fitting.Black dashed lines depict common tangents indicating two-phase coexistence.Blue shaded areas mark the remaining single phase regions.The DFT data in the top row corresponds to the SCAN+U+rVV10 results shown in Fig. 4b.In the bottom row, the measured V OC from GITT is shown for comparison (green markers).

Fig. 6 |
Fig. 6 | C-rate tests performed for C/10, C/5, 1C and 5C.a-d Experimental and (e) simulated results.The black solid line in (a)-(d) are the relaxed OCV points from the GITT measurement in Fig. 2. The black reference in subfigure (e) is the equilibrium potential from Fig. 5b.Inset pictures show the local filling fraction of sodium ions for charging with 1C.

Fig. 7 |
Fig. 7 | Schematic of the phase transformation mechanism.Bottom row shows the simulated phase boundary (white dashed lines) and local composition x at various C-rates but equal average composition X = 0.18.A transition from the layer-by-layer mechanism to the shrinking core is observed with increasing charging rate.

Fig. 8 |
Fig. 8 | Schematic of the micro-battery model.a Single Na X TMO 2 cathode particle cycled against Na metal anode assuming fast transport in the electrolyte and fast reaction kinetics at the anode side.b Intercalation fluxes can be highly inhomogenous depending on phase transformations and the site filling fraction x surf at the surface.
Table. 1.The PBE+U+D3 calculations overestimate the OCV by about 30-40 mV while the calculations based on the SCAN and SCAN+rVV10 functionals underestimate the experimental OCV by about 15 mV.In Table 2, the c lattice constants and band gaps for the lowest energy structures of O2-Na 0 Ni 1/3 Mn 2/3 O 2 , P2-Na 1/3 Ni 1/3 Mn 2/3 O 2 and P2-Na 2/3 Ni 1/3 Mn 2/3 O 2 obtained by the different E xc functionals are

Table 2
Experimental c-lattice parameters are included for comparison.