Defect Chemistry and Na-Ion Diffusion in Na3Fe2(PO4)3 Cathode Material

In this work, we employ computational modeling techniques to study the defect chemistry, Na ion diffusion paths, and dopant properties in sodium iron phosphate [Na3Fe2(PO4)3] cathode material. The lowest intrinsic defect energy process (0.45 eV/defect) is calculated to be the Na Frenkel, which ensures the formation of Na vacancies required for the vacancy-assisted Na ion diffusion. A small percentage of Na-Fe anti-site defects would be expected in Na3Fe2(PO4)3 at high temperatures. Long-range diffusion of Na is found to be low and its activation energy is calculated to be 0.45 eV. Isovalent dopants Sc, La, Gd, and Y on the Fe site are exoergic, meaning that they can be substituted experimentally and should be examined further. The formation of Na vacancies and Na interstitials in this material can be facilitated by doping with Zr on the Fe site and Si on the P site, respectively.


Introduction
Rechargeable sodium ion batteries (SIBs) have gained considerable attention for the development of large-scale energy storage applications due to the high abundance, low cost, and non-toxicity of sodium [1][2][3][4][5]. In practice, there are only a few electrode materials that have been reported. This is due to the larger ionic radius of Na compared to that of Li. Designing a new Na-based electrode material consisting of an appropriate transition metal with a high electrochemical performance could make this material promising.
A number of iron-based phosphate cathode materials [6][7][8][9][10] were proposed for Li-ion batteries due to the structural stability and high redox potential provided by the PO 4 3− matrix. Although similar iron-based cathode materials can be prepared for NIBs in theory, only a few of them, including NaFePO 4 [11], Na 2 FeP 2 O 7 [12], Na 3 V 2 (PO 4 ) 3 [13] and Na 4 Fe 3 (PO 4 ) 2 P 2 O 7 [14], have been reported in the literature. Recently, monoclinic phase Na 3 Fe 2 (PO 4 ) 3 was synthesized using the solid state method and examined as a cathode material for SIBs [15]. This material showed a very high cyclic stability and a reversible discharge capacity of 40 mA·h·g −1 with a flat plateau at about 2.5 V. Rajagopalan et al. [16] studied the electrochemical performance and reversible capacity of Na 3 Fe 2 (PO 4 ) 3 . Their study shows that the discharge specific capacity and cycling stability can be improved by wrapping Na 3 Fe 2 (PO 4 ) 3 with conducting carbon. There are no further experimental or theoretical studies available on Na 3 Fe 2 (PO 4 ) 3 for the use of this material in rechargeable SIBs.
A fundamental understanding of Na 3 Fe 2 (PO 4 ) 3 gained through computational simulation techniques based on the classical pair-potentials can be used to optimize its performance as these techniques have been useful in experimental characterization, the prediction of ion pathways, and the the determination of promising dopants. This technique has been used in a variety of oxide materials, including electrode materials for lithium and sodium ion batteries . In this work, we have used classical pair-potential simulation to examine possible defects that can be observed in Na3Fe2(PO4)3, Na ion migration pathways and the solution of trivalent dopants (Al 3+ , Ga 3+ , Sc 3+ , Y 3+ , Gd 3+ and La 3+ ) on the Fe site and tetravalent dopants (Si 4+ , Ge 4+ , Ti 4+ , Sn 4+ , Zr 4+ and Ce 4+ ) on the Fe and P sites.

Computational Methods
The calculations employed in this study are based on the classical pair wise potentials. The General Utility Lattice Program (GULP) [39] was used. This code uses ion-ion interactions in the form of long-range (i.e., Coulombic) attraction and short-range repulsion (i.e., Pauli electron-electron) and attraction (dispersion). We used the Buckingham potentials (refer to Table S1 in the supplementary information) to model short-range interaction.
Bulk Na3Fe2(PO4)3 and defect configurations were optimized using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [40], as implemented in the GULP code. The forces on the atoms were below 0.001 eV/Å in all cases. The Mott-Littleton method [41] was employed to model point defects. This methodology has been well-explained in previous studies. The current simulation is expected to overestimate the defect enthalpies. This is due to the spherical shape of the ions with a low concentration. Nevertheless, trends in relative energies are expected to be consistent.
The present atomistic simulations use isobaric parameters in the calculations of formation and migration energies. The detailed thermodynamic relations associated with isobaric parameters are discussed in previous theoretical work [42][43][44][45][46][47].

Bulk Na3Fe2(PO4)3 Structure
Bulk Na3Fe2 (PO4)3 belongs to the monoclinic crystal system (space group C2/c). The crystal structure of Na3Fe2(PO4)3 is shown in Figure 1. Its lattice parameters (a = 15.070 Å, b = 8.740 Å, c = 8.724 Å, α = γ = 90.0° and β = 125.1°) were determined by Fanjet et al. [48] in their powder neutron diffraction study. The crystal structure consists of FeO6 octahedra and PO4 tetrahedra units sharing their corners. In order to validate the interatomic potentials used in this study, a bulk Na3Fe2(PO4)3 structure was optimized under constant pressure. The difference between the calculated equilibrium lattice constants and its corresponding experimental values is less than 1.5%, showing the good In order to validate the interatomic potentials used in this study, a bulk Na 3 Fe 2 (PO 4 ) 3 structure was optimized under constant pressure. The difference between the calculated equilibrium lattice constants and its corresponding experimental values is less than 1.5%, showing the good reproduction and the suitability of these potential parameters for defect modeling. Table 1 lists the calculated and experimental values together with the error percentages.

Intrinsic Defect Process
Here, we consider the formation of point defects (vacancies and interstitials) to calculate the Schottky and Frenkel defect processes. The Na-Fe anti-site defect process is also considered. The intrinsic point defects are important as they stimulate the ions to diffuse in the lattice. The point defects were combined to construct reaction energy processes (Schottky, Frenkel and anti-site) using Kröger-Vink notation [49]. The reaction equations are as follows: Na/Fe antisite (cluster) : Na X Na + Fe X Fe → Na Fe : The energetics for the defect processes are shown in Figure 2. Our calculations show that the Na Frenkel is the lowest defect energy process (0.45 eV/defect), indicating that the formation of Na vacancies is facilitated by this process. Thus, this process would accelerate the vacancy-aided Na diffusion. The Na-Fe anti-site defect is found to be the second most favorable defect energy process (1.12 eV/defect). In this defect, a small population of Na on the Fe site and Fe on the Na site would be observed. This defect has been found in a variety of battery materials both experimentally and theoretically [17][18][19][50][51][52][53]. The formation of Na 2 O is calculated to be a 2.46 eV/defect, suggesting that there is a possibility of Na 2 O loss in this material at high temperatures. Other defect processes are relatively high energy, meaning that they cannot be observed under standard battery operating conditions.

Sodium Ion Diffusion
The diffusion of Na ions with a low migration barrier is a key requirement for a promising cathode material. There is no experimental report on the Na ion migration pathways in Na3Fe2(PO4)3. Using the current methodology, it is possible to calculate the Na ion diffusion pathways and activation energies. Three different local Na vacancy migration hops (refer to Figure 3) were calculated. Table 2 lists the Na-Na separations together with the activation energies. Energy profile diagrams together with activation energies for local Na hops are shown in Figure 4. The lowest activation energy (0.44 eV) is calculated for hop A in which Na migrates along the bc plane. Hop B exhibits an activation energy of 0.45 eV, which is very close to the value calculated for hop B. The highest activation energy (2.37 eV) is calculated for hop C. Both hops B and C are also along the bc plane.

Sodium Ion Diffusion
The diffusion of Na ions with a low migration barrier is a key requirement for a promising cathode material. There is no experimental report on the Na ion migration pathways in Na 3 Fe 2 (PO 4 ) 3 . Using the current methodology, it is possible to calculate the Na ion diffusion pathways and activation energies. Three different local Na vacancy migration hops (refer to Figure 3) were calculated. Table 2 lists the Na-Na separations together with the activation energies. Energy profile diagrams together with activation energies for local Na hops are shown in Figure 4. The lowest activation energy (0.44 eV) is calculated for hop A in which Na migrates along the bc plane. Hop B exhibits an activation energy of 0.45 eV, which is very close to the value calculated for hop B. The highest activation energy (2.37 eV) is calculated for hop C. Both hops B and C are also along the bc plane.

Sodium Ion Diffusion
The diffusion of Na ions with a low migration barrier is a key requirement for a promising cathode material. There is no experimental report on the Na ion migration pathways in Na3Fe2(PO4)3. Using the current methodology, it is possible to calculate the Na ion diffusion pathways and activation energies. Three different local Na vacancy migration hops (refer to Figure 3) were calculated. Table 2 lists the Na-Na separations together with the activation energies. Energy profile diagrams together with activation energies for local Na hops are shown in Figure 4. The lowest activation energy (0.44 eV) is calculated for hop A in which Na migrates along the bc plane. Hop B exhibits an activation energy of 0.45 eV, which is very close to the value calculated for hop B. The highest activation energy (2.37 eV) is calculated for hop C. Both hops B and C are also along the bc plane.    Table 2. Na-Na separations and their activation energies for the Na ion migration, as shown in Figure  3.

Migration Path
Na-Na Separation (Å) Activation Energy ( We examined possible long-range Na ion pathways linked by local Na hops. Three possible long-range Na ion paths and their corresponding overall activation energies are listed in Table 3 Clark et al. [54] used classical atomistic simulation to calculate the activation energy of Na ions in Na2FeP2O7. In their study, three-dimensional long-range Na ion migration paths were observed in different directions with activation energies in the range of 0.33-0.49 eV. Sodium ion migration paths together with activation energies were calculated in "olivine" NaFePO4 material [55]. The lowest energy path was observed along the [010] direction with the activation energy of 0.30 eV. The current activation energy value of 0.45 eV calculated for long-range Na ion migration in Na3Fe2(PO4)3 suggests that this material is also a promising cathode material competitive with the other Fe-based polyanions materials for Na-ion batteries. Table 3. Long-range Na ion diffusion paths and their overall activation energies.

Long-Range Path Direction Overall Activation Energy (eV)
A→A→B→B bc plane 0.45 A→C→C→ A bc plane 2.37 C→C→C→C bc plane 2.37

Isovalent Doping
Solutions of isovalent dopants (Al, Ga, Sc, Y, Gd, and La) were considered on the Fe site. The selection of trivalent dopants with a wide range of ionic radii is based on previous studies that considered these dopants in different materials, including phosphate-based battery materials [23,[56][57][58]. Though they exhibit a high atomic weight, low abundance, or high cost, low-level doping of Na3Fe2(PO4)3 can exhibit improvement in electronic conductivity. The solution enthalpy was calculated using the following reaction: Figure 5 reports the solution enthalpies of M 3+ ions on the Fe site. The lowest solution energy is calculated for Sc. Solutions of La and Gd are also promising as they exhibit negative values. Yttrium shows a small but negative solution enthalpy. The solution enthalpy of Al is highly endoergic, suggesting that it cannot be doped under normal conditions. Dopants exhibiting endoergic solution  Table 2. Na-Na separations and their activation energies for the Na ion migration, as shown in Figure 3.

Migration Path
Na-Na Separation (Å) Activation Energy (eV) We examined possible long-range Na ion pathways linked by local Na hops. Three possible long-range Na ion paths and their corresponding overall activation energies are listed in Table 3. Clark et al. [54] used classical atomistic simulation to calculate the activation energy of Na ions in Na 2 FeP 2 O 7 . In their study, three-dimensional long-range Na ion migration paths were observed in different directions with activation energies in the range of 0.33-0.49 eV. Sodium ion migration paths together with activation energies were calculated in "olivine" NaFePO 4 material [55]. The lowest energy path was observed along the [010] direction with the activation energy of 0.30 eV. The current activation energy value of 0.45 eV calculated for long-range Na ion migration in Na 3 Fe 2 (PO 4 ) 3 suggests that this material is also a promising cathode material competitive with the other Fe-based polyanions materials for Na-ion batteries.

Isovalent Doping
Solutions of isovalent dopants (Al, Ga, Sc, Y, Gd, and La) were considered on the Fe site. The selection of trivalent dopants with a wide range of ionic radii is based on previous studies that considered these dopants in different materials, including phosphate-based battery materials [23,[56][57][58]. Though they exhibit a high atomic weight, low abundance, or high cost, low-level doping of Na 3 Fe 2 (PO 4 ) 3 can exhibit improvement in electronic conductivity. The solution enthalpy was calculated using the following reaction: Figure 5 reports the solution enthalpies of M 3+ ions on the Fe site. The lowest solution energy is calculated for Sc. Solutions of La and Gd are also promising as they exhibit negative values.
Yttrium shows a small but negative solution enthalpy. The solution enthalpy of Al is highly endoergic, suggesting that it cannot be doped under normal conditions. Dopants exhibiting endoergic solution enthalpies can be doped experimentally to prepare [Na 3 (Fe x M 1−x ) 2 (PO 4 ) 3 ; M = Sc, La, Gd and Y; x = 0-1] composites. Such composites may have the different chemical, electronic, and mechanical properties required for different purposes. Gaining knowledge on the exact composition requires experiments.

Tetravalent Doping
Both the Fe and P sites were considered for tetravalent dopants (Si, Ge, Ti, Sn, Zr, and C). It is possible to form Na vacancies by doping with M 4+ ions on the Fe site. Such vacancies facilitate Na ion migration in this material. The following reaction equation was constructed: Solution enthalpies are shown in Figure 6a. In all cases, high solution enthalpies (>3 eV) are observed, suggesting that the formation of Na vacancies are unlikely at normal temperatures. The lowest solution enthalpy is calculated for Zr. Solution enthalpies for Ge and Sn are very close to the value calculated for Zr. The Si exhibits the highest solution enthalpy. Doping with M 4+ ions on the P site can lead to the formation of Na interstitials, as shown in Equation (11). This engineering strategy in turn would enhance the capacity of Na3Fe2(PO4)3.

Tetravalent Doping
Both the Fe and P sites were considered for tetravalent dopants (Si, Ge, Ti, Sn, Zr, and C). It is possible to form Na vacancies by doping with M 4+ ions on the Fe site. Such vacancies facilitate Na ion migration in this material. The following reaction equation was constructed: Solution enthalpies are shown in Figure 6a. In all cases, high solution enthalpies (>3 eV) are observed, suggesting that the formation of Na vacancies are unlikely at normal temperatures. The lowest solution enthalpy is calculated for Zr. Solution enthalpies for Ge and Sn are very close to the value calculated for Zr. The Si exhibits the highest solution enthalpy.

Tetravalent Doping
Both the Fe and P sites were considered for tetravalent dopants (Si, Ge, Ti, Sn, Zr, and C). It is possible to form Na vacancies by doping with M 4+ ions on the Fe site. Such vacancies facilitate Na ion migration in this material. The following reaction equation was constructed: Solution enthalpies are shown in Figure 6a. In all cases, high solution enthalpies (>3 eV) are observed, suggesting that the formation of Na vacancies are unlikely at normal temperatures. The lowest solution enthalpy is calculated for Zr. Solution enthalpies for Ge and Sn are very close to the value calculated for Zr. The Si exhibits the highest solution enthalpy. Doping with M 4+ ions on the P site can lead to the formation of Na interstitials, as shown in Equation (11). This engineering strategy in turn would enhance the capacity of Na3Fe2(PO4)3.
2RO + 2P + Na O → 2R + 2Na • + P O (11) Doping with M 4+ ions on the P site can lead to the formation of Na interstitials, as shown in Equation (11). This engineering strategy in turn would enhance the capacity of Na 3 Fe 2 (PO 4 ) 3 .
Solution enthalpies calculated for this process are shown in Figure 6b. The promising candidate is Si as its solution enthalpy is 0.41 eV. The Ge exhibits a slightly higher (1.23 eV) solution enthalpy. Other dopants have solution enthalpies greater than 2.50 eV, meaning that they should be doped at high temperatures. The highest solution enthalpy (5.98 eV) is calculated for Ce.

Conclusions
Atomistic simulation techniques were employed to examine the defects, Na ion migration pathways, and a variety of isovalent and isovalent dopants on the Fe site in Na 3 Fe 2 (PO 4 ) 3 . The lowest energy defect process is the Na Frenkel, suggesting that the Na diffusion in this material would be assisted by Na vacancies. The Na-Fe anti-site defect is calculated to be the second lowest energy process, meaning that a small population of Na and Fe ions exchange their positions. The diffusion of Na is calculated to be low in this material and Na ion migrates via the bc plane with the activation energy of 0.45 eV. The favorable isovalent dopants on the Fe site are Sc, La, Gd, and Y, meaning that the synthesis of [Na 3 (Fe x M 1−x ) 2 (PO 4 ) 3 M = Sc, La, Gd and Y] is worth investigating experimentally. Doping with Zr on the Fe site can increase the concentration of Na vacancies needed for the Na ion diffusion, while doping with Si on the P site can facilitate the formation of Na interstitials required for the improvement in the capacity of Na 3 Fe 2 (PO 4 ) 3 .
Funding: This research was financially supported by the European Union's H2020 Programme under Grant Agreement no. 824072-HARVESTORE.