An ab initio investigation on the electronic structure, defect energetics, and magnesium kinetics in Mg₃Bi₂

Abstract

We present a comprehensive ab initio investigation on Mg₃Bi₂, a promising Mg-ion battery anode material with high rate capacity. Through combined DFT (PBE, HSE06) and G₀W₀ electronic structure calculations, we find that Mg₃Bi₂ is likely to be a small band gap semiconductor. DFT-based defect formation energies indicate that Mg vacancies are likely to form in this material, with relativistic spin–orbit coupling significantly lowering the defect formation energies. We show that a transition state searching methodology based on the hybrid eigenvector-following approach can be used effectively to search for the transition states in cases where full spin–orbit coupling is included. Mg migration barriers found through this hybrid eigenvector-following approach indicate that spin–orbit coupling also lowers the migration barrier, decreasing it to a value of 0.34 eV with spin–orbit coupling. Finally, recent experimental results on Mg diffusion are compared to the DFT results and show good agreement. This work demonstrates that vacancy defects and the inclusion of relativistic spin–orbit coupling in the calculations have a profound effect in Mg diffusion in this material. It also sheds light on the importance of relativistic spin–orbit coupling in studying similar battery systems where heavy elements play a crucial role.

---

1 Introduction

Extensive research has been performed on lithium-ion battery (LIB) systems in the last few decades, resulting in their integration into electric vehicles and handheld devices.¹ However, supply restrictions, high cost, and safety limitations of LIBs have kindled research into battery systems using alternative metal ions. Recently, batteries based on multivalent ion chemistries such as Mg²⁺, Al³⁺, and Ca²⁺ have gained attention due to their low cost and high volumetric densities, features important for large-scale applications such as electric vehicles and grid-based energy storage systems.² Of these, Mg-ion batteries (MIBs) are an attractive system due to the high abundance of Mg in the Earth's crust and the large volumetric capacity of Mg metal. However, currently at least three factors hinder the application of working MIB systems: (1) absence of electrolytes compatible with high-voltage cathodes, (2) slow diffusion of divalent Mg²⁺ ions in electrode materials, and (3) passivation of Mg metal surfaces due to electrolyte decomposition reactions and the concomitant difficulties in stripping and plating Mg.³

While possible high-voltage MIB cathode candidates have been identified,² it remains extremely challenging to identify an electrolyte that is suitable for use with both high voltage cathodes and Mg metal. This incompatibility of electrolytes that are suitable for high voltage cathodes with Mg metal and vice versa has been a problem in developing a working MIB system. Hence, a number of studies have focused on using an alternative anode material other than Mg to bypass this problem.³ In this respect, bismuth metal is a promising anode material with a low discharge voltage of 0.2 V versus Mg metal.^4–6 Bismuth alloys with magnesium to form an intermetallic Mg₃Bi₂ phase, resulting in a theoretical capacity of 385 mA h g⁻¹. However, the most interesting feature of a bismuth anode is that it displays fast Mg ion insertion and de-insertion, a feature not present in most other Mg-ion electrodes.^4,7 Recently, our group has used ²⁵Mg Nuclear Magnetic Resonance (NMR) spectroscopy to study Mg₃Bi₂ and reported fast exchange between the two crystallographically distinct Mg sites in this material; a dependence of the exchange rate and activation energy on the preparation condition (ball-milling and electrochemical synthesis) was also demonstrated.⁸

As ion dynamics are properties intimately related to the defect chemistry (which in turn depends heavily on the preparation conditions) and ultimately the electronic structure, we present an in-depth ab initio investigation on the Mg₃Bi₂ system, the end product of the discharge process. The results are presented in four sections: in the first section, we discuss the structural characteristics of Mg₃Bi₂, identifying the possible Mg diffusion pathways. Next, electronic structures and defect energetics are described. Due to the well known problem of band gap underestimation in semilocal DFT, we use three different methods (PBE, HSE06, and G₀W₀) to calculate the electronic structure. In particular, the Greens function based G₀W₀ method is used to give the most accurate value for the minimum band gap, overcoming the band gap problem in semilocal DFT.

Bismuth is a heavy element where relativistic spin–orbit coupling is expected to have a significant influence on the properties relevant to battery chemistry, as shown for lead in a recent ab initio study of the lead–acid battery.⁹ Thus, calculations were performed with and without the explicit spin–orbit coupling Hamiltonian. Inclusion of spin–orbit coupling was found to reduce the band gap, defect formation energies, and also the Mg migration barrier which was calculated using the hybrid eigenvector-following approach. Using these results we are able to reconcile and rationalize previous reports on the Mg mobility in Mg₃Bi₂, in which very different Mg mobilities were assumed to occur due to the different concentrations of Mg vacancies in the samples. We also demonstrate that the hybrid eigenvector-following method can be a very efficient approach for locating transition states in systems where spin–orbit coupling is likely to play an important role.

2 Methodology

2.1 Computational details

All Density Functional Theory (DFT) calculations were performed with the VASP code^10,11 employing the projector-augmented wave (PAW) method.¹² Spin-polarized Perdew–Burke–Ernzerhof (PBE) and Heyd–Scuseria–Ernzerhof (HSE06) exchange–correlation functionals were adopted.^13,14 For the energy and force calculations, PAW pseudopotentials treating the Mg 3s and Bi 5d as valence states were used, with a plane-wave basis cutoff of 350 eV. All lattice relaxations were performed with 1.3 times the ENMAX value as defined in the pseudopotential file. An additional support grid was used for the evaluation of augmented charges (ADDGRID =.TRUE.). Self-consistent field (SCF) cycles were converged with an energy tolerance of 10⁻⁴ eV. Monkhorst–Pack k-point sampling of <0.05 Å⁻¹ was used in the Brillouin zone. For the density of states and band structure calculations, Mg 2p states were also treated as valence states in the PAW potential and an increased cutoff of 550 eV was used. SCF cycles were converged to a 10⁻⁶ eV limit. A Γ-centered k-point sampling of <0.03 Å⁻¹ was used.

Single-shot G₀W₀ calculations were performed with 250 frequency grid points and a 360 eV plane-wave cutoff for the response function calculations. HSE06 wavefunctions were used as starting wavefunctions for the G₀W₀ calculations. Band gap convergence with respect to the number of frequency grid points, number of empty bands, and plane-wave cutoffs were checked. Quasiparticle energy iterations were converged to a 10⁻⁸ eV limit.

Relativistic corrections to the electronic structure were taken into account through two levels of theory: a ‘scalar relativistic’ correction which only includes the mass-velocity and Darwin terms of the relativistic Hamiltonian, as implemented in VASP;¹⁵ and a ‘full relativistic’ correction that explicitly includes the spin–orbit coupling term alongside the scalar relativistic correction.

The initial structure of Mg₃Bi₂ was fully relaxed until the energy differences between the subsequent steps are converged to 10⁻⁵ eV per cell and the forces are <0.05 eV Å⁻¹. For the defect calculations, two supercells each containing 40 atoms (16 Mg, 24 Bi) and 120 atoms (48 Mg, 72 Bi) were created and again relaxed to the above criteria. Defect notations follow that of Kroger and Vink, with charge superscripts omitted for clarity (e.g. Bi_Mg(oct) denotes a neutral Bi sitting on an octahedral Mg site).¹⁶ All defect energies were referenced to the respective bulk metals (Mg, Bi) at the same level of theory as the defect calculations. Only the neutral defects are considered since the system has a vanishing band gap under the PBE level of theory, with the metallic behaviour enhanced under the inclusion of spin–orbit coupling. Whereas the formation energies under the HSE06 level of theory (where a finite band gap was observed) could give a chemical potential dependence, this was not possible due to the computational resources available.

2.2 Transition state searching

Transition state and steepest-descent pathway were found using the hybrid eigenvector-following approach as implemented in the OPTIM code.¹⁷ Details on the algorithm are available elsewhere.^18,19 In brief, this method finds the smallest eigenvalue of the Hessian matrix without calculating the full matrix, by minimizing the Rayleigh–Ritz ratio. Here we use a low-memory Broyden–Fletcher–Goldfarb–Shannon (LBFGS) scheme to minimize this ratio on a given point to find the uphill path. The minimization was deemed to have occurred when the root-mean-square gradient is smaller than 0.025 eV Å⁻¹. Then up to five LBFGS minimization steps are performed in the tangent space until the root-mean-square gradient is smaller than 10⁻³ eV Å⁻¹. Initial guesses of the transition state structure were produced by removing a nearby Mg atom and putting the migrating Mg atom in the middle of the proposed diffusion path. The steepest descent pathway from the transition state is found by displacing the atoms by 0.1 Å from the transition state in the parallel and antiparallel directions to the eigenvector. Local minima are then found by the LBFGS algorithm with energy convergence of 10⁻³ eV.

We note that the transition state geometry and eigenvectors obtained from the significantly cheaper scalar relativistic calculations could be re-used as an initial guess for the fully relativistic calculations. We have observed that only <5 force evaluations are typically needed for the relativistic calculations when this approach is taken, which significantly reduces the computational requirement. In addition, unlike the classical nudged elastic band (NEB) method where only a few points along the pathway are sampled, the HEF method can sample the path at an arbitrary step size. This allows us to capture the finer details of the energy variation along the pathway.

3 Results and discussion

3.1 Crystal structure

Mg₃Bi₂ is the last in a series of magnesium pnictide compounds Mg₃X₂ (X = P, As, Sb, Bi) and is a Zintl-type compound with closed shell ions with formal charges of Mg²⁺ and Bi³⁻. It adopts an anti-La₂O₃ type structure with hexagonal P3̄m1 space-group symmetry.²⁰ The structure is schematically illustrated in Fig. 1. In terms of ion arrangements, Mg₃Bi₂ incorporates two alternating layers: layer A consists of tetrahedrally coordinated Mg²⁺ cations and octahedral interstitial sites forming ‘covalent’ Mg₂Bi₂²⁻ sheets whereas layer B consists of octahedrally coordinated Mg²⁺ cations and tetrahedral interstitial sites. These interstitial sites are expected to play an important role in ionic diffusion²¹ as discussed further below.

Fig. 1 (a) Crystal structure of Mg₃Bi₂. Positions of tetrahedrally coordinated Mg (orange sphere, Tet) and octahedrally coordinated Mg (green sphere, Oct) are shown. Bismuth atoms sit on the line vertices (not shown for clarity). (b) Schematic illustration of the Mg₃Bi₂ structure showing the possible diffusion pathways involving the tetrahedral and the octahedral interstitial sites (‘Tet Int’ and ‘Oct Int’).

From DFT-based lattice relaxation, we can see that the PBE functional overestimates the experimentally determined a cell parameter²² at both the scalar relativistic and full relativistic levels of theory, whereas it marginally underestimates the c parameter at the scalar relativistic level, and overestimates it at the full relativistic level (Table 1). The HSE06 functional again overestimates the a lattice parameter, although to a smaller extent than PBE, whereas it underestimates the c lattice parameter. Finally, inclusion of SOC results in slight expansion of the cell which may arise from reduced electrostatic interactions (see below).

Table 1 DFT-predicted cell parameters of Mg₃Bi₂ using the PBE and HSE06 exchange–correlation functionals. SR and FR refer to scalar relativistic and full relativistic (i.e. explicit spin–orbit coupling) calculations, respectively. Experimental lattice constants are from Lazarev et al.²²

	PBE		HSE06		Expt
	SR	FR	SR	FR
a/Å	4.71	4.72	4.65	4.66	4.62
c/Å	7.40	7.44	7.35	7.36	7.41

---

Finally, we note that Mg₃Bi₂ undergoes a phase transition above 703 °C to a defective body-centered cubic structure similar to AgI, with excess Mg cations as the tetrahedral interstitials, resulting in Mg_1.5Bi.²³ As expected from the AgI structure, it shows superionic conduction of Mg²⁺ cations as determined through neutron diffraction.²⁴ In this study, however, we restrict the investigation to the room-temperature hexagonal phase which is more relevant to Mg-ion batteries.

3.2 Electronic structure

Despite being first reported in 1933 by Zintl,²⁰ the precise electronic structure and the band gap (E_g) of Mg₃Bi₂ has not been determined with either experiment or theoretical calculations. In a series of work, Ferrier and co-workers have speculated that Mg₃Bi₂ is semimetallic based on a conductivity-composition plot;^25,26 Lazarev et al. have assumed Mg₃Bi₂ is semiconducting with E_g = 0.1 eV;²² Watson et al. obtained E_g = 0.5 eV from Mg X-ray emission spectra.²⁷ Work on the amorphous Mg₃Bi₂ alloy have shown from conductivity measurements that it is semiconducting with a band gap of 0.15 eV.²⁸ These varying results may arise from the difficulty in preparing these materials stoichiometrically, due to the high vapor pressure of Mg. Previous literature describing electronic structure calculations is also sparse: Sedighi et al. concluded Mg₃Bi₂ is a semiconductor with E_g = 0.25 eV based on the Engel–Vosko Generalized Gradient Approximation (EV-GGA),²⁹ whereas Imai et al., Xu et al., and Zhang et al. concluded it is a semimetal based on pure GGA calculations.^30–32 However, all these works except the last one did not consider explicit spin–orbit coupling (SOC), which has been shown to play an important role in the electronic structure of compounds involving heavy atoms such as bismuth;^33,34 in addition, the well known problem of band gap underestimation in semilocal DFT necessitates the use of a higher level of theory to predict the accurate electronic structure of this material. Hence, we revisit the electronic structure of this compound using state-of-the-art electronic structure methods such as hybrid functionals and many-body perturbation theory in the G₀W₀ approximation, and also consider the effects of spin–orbit coupling.

First we look at the Bader charges of Mg and Bi in the structure, shown in Table 2. The HSE06 functional results in an increased ionicity in the system compared to the semilocal PBE values. In all cases, inclusion of SOC resulted in decreased ionicity which is consistent with the density of states data shown below.

Table 2 Scalar relativistic (SR) and full relativistic (FR) Bader charge analysis of Mg₃Bi₂ using the Perdew–Burke–Ernzerhof (PBE) and Heyd–Scuseria–Ernzerhof (HSE06) exchange–correlation functionals. Only the valence charge is calculated

	Bi	Mgoct	Mgtet
PBE SR	−2.10	+1.42	+1.39
PBE FR	−2.06	+1.39	+1.37
HSE06 SR	−2.21	+1.51	+1.45
HSE06 FR	−2.18	+1.48	+1.44

---

Fig. 2 shows the density of states (DOS) plots obtained for Mg₃Bi₂ using DFT within the semilocal PBE and hybrid HSE06 approximations. Both methods show the same trend comparing the results with and without SOC. As expected from the charge distribution of Mg₃Bi₂, the top of valence band is strongly dominated by the bismuth 6p contribution down to around −5 eV from the Fermi level and has a negligible Mg contribution. As expected from the literature, the relativistic contraction of the bismuth 6s states around −12 eV from the Fermi level results in a relatively large separation of around 5 eV between the top of 6s and the bottom of 6p states, which makes it chemically inactive (the so-called inert-pair effect).³⁵ The downshift in 6s energy by adding the full relativistic effect is around −0.15 eV; effects of similar magnitudes were observed in PbO₂, an important active material in lead–acid batteries.⁹

Fig. 2 Density of states plot of Mg₃Bi₂ using the PBE and HSE06 exchange–correlation functionals. Local atomic DOS projections inside the sphere defined by the Wigner–Seitz radii (1.63 and 1.52 Å for Bi and Mg, respectively) are also shown. All energies were referenced to the highest occupied state.

Fig. 3 shows the band structure of Mg₃Bi₂ using DFT within the semilocal PBE and hybrid HSE06 approximations, with and without the spin–orbit interaction. At the PBE level, Mg₃Bi₂ exhibits semimetallic behavior, consistent with earlier reports of a type-II nodal line semimetal at this level of theory.³² The inclusion of spin–orbit coupling opens a small gap in the nodal line, but the system remains semimetallic. At the HSE06 level, Mg₃Bi₂ is a semiconductor with the minimum band gap of 0.36 eV located at the Γ point. The inclusion of spin–orbit coupling reduces the band gap to 0.17 eV, but the system remains semiconducting. The G₀W₀ calculations show that the Γ-point band gap increases to 0.58 eV without spin–orbit coupling. We expect that spin–orbit coupling would reduce the band gap by an amount similar to that observed in HSE06, and thus Mg₃Bi₂ is also semiconducting at the G₀W₀ level of theory.

Fig. 3 Band structure plots of Mg₃Bi₂ using the PBE (left) and HSE06 (right) exchange–correlation functionals without (grey) and with (red) spin–orbit coupling (SOC). The Fermi level is located at zero energy.

We note that inclusion of SOC results in a reduced band gap and increased dispersion of the Bi 6p states. Previous reports on binary Sb and Te compounds have attributed this phenomenon to spin–orbit coupling in the anion p-orbitals: energy splitting of the degenerate p-states result in the j = 3/2 state being split upward in energy, and the j = 1/2 state being split downward in energy.³⁶ The net effect of this splitting is a reduction in the band gap.

3.3 Defect energetics

To investigate the possibility of different types of defects present in the Mg₃Bi₂ sample, formation energies ΔE_f of various stoichiometric and non-stoichiometric point defects were calculated ab initio. The results outlined in Table 3 are separated into three categories: antisite defects, Frenkel defects, and vacancy defects, which we discuss in sequence. Nonstoichiometric antisite and vacancy defects were created by removing or inserting relevant atoms from the stoichiometric cell while retaining the charge neutrality of the cell due to Mg₃Bi₂ being a semimetal at the PBE level of theory (i.e. charged defects are not possible under such a situation). For instance, a Mg vacancy defect would involve removal of a neutral Mg atom resulting in a composition of Mg₂₃Bi₁₆ (for a 40-atom supercell) and two electron holes on the Fermi level; similarly, a Mg substitutional defect would involve a Bi atom replaced for a Mg atom, resulting in Mg₂₅Bi₁₅ (again for a 40-atom supercell) and five excess electrons. Charge neutrality is automatically maintained in the antisite and Frenkel defects since they are stoichiometric.

Table 3 Scalar relativistic (SR) and full relativistic (FR) ab initio formation energies ΔE_f of various stoichiometric and non-stoichiometric defects in Mg₃Bi₂ using the PBE functional. Defect notations follow the convention of Kroger and Vink with neutral sign omitted for clarity.¹⁶ All calculations assumed non-charged defects (see text) with cell dimensions fixed to simulate a dilute limit. Vacancy defect energies are referenced to the respective metals. For Frenkel defects, two scenarios where the Mg sits on a nearby (nn) or far interstitial sites were considered. Asterisks (*) indicate that the resulting structure was unstable and reverted back to the starting structure. All values are quoted in electron-volts. Only some of the calculations were performed under the 120-atom supercell condition after an initial screening with the 40-atom supercell

	ΔEf
	40-atom supercell		120-atom supercell
	SR	FR	SR	FR
Antisite defects
BiMg(oct)	1.84
BiMg(tet)	2.01
MgBi	2.92
BiMg(oct) + MgBi	2.51	2.27
BiMg(tet) + MgBi	1.17	1.15
Frenkel defects
VMg(oct) + Mgi(oct) (nn)	0.87	0.79	0.89
VMg(oct) + Mgi(oct) (far)	0.89	0.81	1.04
VMg(oct) + Mgi(tet) (nn)	*	*	*	*
VMg(oct) + Mgi(tet) (far)	1.49	1.41	1.55
VMg(tet) + Mgi(oct) (nn)	*	*	*	*
VMg(tet) + Mgi(oct) (far)	1.07	0.98	1.10
VMg(tet) + Mgi(tet) (nn)	*	*	*	*
VMg(tet) + Mgi(tet) (far)	1.17	1.15	1.61
Vacancy defects
VMg(oct)	1.09	0.40	1.03	0.33
VMg(tet)	1.17	0.50	1.12	0.42
VBi	1.88	2.02

---

First, the formation of antisite defects are energetically unfavorable, with the Mg on a Bi site having the highest formation energy of 2.92 eV. This is most likely due to the large difference between their ionic radii. Despite the fact that no ionic radius for the Bi³⁻ ion is reported, the trend can still be explained in terms of the corresponding atomic radii of Mg and Bi (Mg 1.50 Å; Bi 1.60 Å), the difference being expected to get even larger as Mg is oxidized and Bi is reduced. This is also supported from the Mg–Bi bond lengths in Mg₃Bi₂ crystal: Mg_oct–Bi 3.21 Å versus Mg_tet–Bi 2.92 Å. Using the Shannon radii of 0.72 and 0.57 Å for Mg_oct and Mg_tet, respectively,³⁷ Bi³⁻ radii of 2.49 and 2.35 Å are obtained in each case. In addition, the unfavourable electrostatic interaction between the ions of same charge (Mg occupying the Bi lattice site is coordinated by Mg ions, and vice versa) reinforces the high energy cost in forming these types of defects.

The ΔE_f of Frenkel-type defects (creation of a vacancy plus an interstitial) are also shown to be relatively high for both nearby and separated vacancy–interstitial pairs; around 0.8–1.2 eV is required for their formation. However, performing structural relaxation on some of the starting guesses with the defect Mg ion sitting on an interstitial site that is adjacent (in the first coordination shell) to a Mg vacancy (specifically, the V_Mg(oct) + Mg_i(tet) (nn), V_Mg(tet) + Mg_i(tet) (nn), and V_Mgi(tet) + Mg_i(tet) (nn) cases) resulted in the structure reverting back to that of the pristine cell, with no defects. This indicates that these defects are energetically unstable and they would revert back to the original structure, if formed at all. These sites, as we will see later, may play an important role as energy maximum transition state sites in the Mg-ion diffusion pathway.

Perhaps the most surprising results are the energies of the vacancy defects: while Mg vacancy defects have around the same ΔE_f as Frenkel defects without SOC, their magnitudes are reduced significantly when SOC is included in the Hamiltonian. From the fully relativistic calculation on the 120-atom supercell with SOC included, the ΔE_f of octahedral vacancies V_Mg(oct) is found to be as low as 0.33 eV, with a slightly higher value of 0.42 eV for tetrahedral vacancies V_Mg(tet). This effect is likely to be connected to the enhanced shielding of Bi 6p-levels due to the SOC effect as explained above: as the formation of a neutral Mg vacancy should involve loss of electrons from the Bi, this deshielding should result in a significant reduction of the formation energy. Hence, considering their low formation energies, we conclude that vacancy defects are the dominant type of defects present in Mg₃Bi₂.

3.4 Mg migration kinetics

Having established the defect chemistry in Mg₃Bi₂, we now turn our attention to the Mg diffusion in this structure and study the migration barriers.

As Mg vacancies were shown to have low formation energy, especially with SOC included, we attempt to simulate the effect of vacancy diffusion following the creation of one octahedral Mg vacancy V_Mg(oct). This has lower ΔE_f than the tetrahedral vacancy V_Mg(tet) as illustrated in the previous section. Then, a nearby tetrahedral Mg atom is removed from its original position and placed on the guessed ‘transition state’. Finally, a hybrid eigenvector-following approach is used to find the transition state, followed by the search for a steepest descent path connecting the two corresponding minima.

The result presented in Fig. 4 clearly shows a small diffusion barrier of hopping between the Mg_oct and Mg_tet (path 1 on Fig. 1b), with 0.34 eV for the fully relativistic calculation and 0.43 eV for the scalar relativistic calculation. In contrast, the Mg_oct–Mg_oct diffusion barrier (path 2 on Fig. 1b) is around twice that of the Mg_oct–Mg_tet diffusion barrier, indicating that the Mg diffusion must occur via octahedral–tetrahedral exchange. This is in line with the conclusion from ²⁵Mg NMR studies of Liu et al.⁸ An alternative exchange mechanism involving path 3 on Fig. 1b was also investigated, but the transition structure searching resulted in the same transition state as for the path 1, clearly indicating the absence of a diffusion pathway along this line.

Fig. 4 Diffusion profile of Mg ion for selected pathways as illustrated in Fig. 1b. (a) Scalar relativistic, (b) full relativistic, with SOC included. Energies are referenced to the bulk energy of each 120-atom supercell with one V_oct defect.

With the ab initio calculations of the defect creation and activation energies, we now compare our results with available experimental data. As reported in the literature, bismuth can be cycled reversibly to form Mg₃Bi₂ in a Mg-ion battery.⁴ Previous work on magnesium ion conduction in bismuth anodes either used the Galvanostatic Intermittent Titration Technique (GITT) on an electrochemical cell, or ex situ variable-temperature (VT) ²⁵Mg NMR spectroscopy to probe Mg transport.^7,8 In this work, however, we restrict the discussion to the latter NMR result, due to the following limitation of GITT experiments: in the GITT measurement, the voltage response (resulting from the relaxation) of the cell after the application of a short current pulse is modelled with the diffusion equation to extract the diffusion coefficient D under the assumption that the (de)insertion reaction occurs via solid solution (i.e. the relaxation of the potential after the current pulse is a measure of ionic transport through a single phase). In situ X-ray diffraction results, however, have clearly shown that magnesiation occurs via a two-phase reaction between Bi and Mg₃Bi₂.⁸ Thus the diffusion coefficients extracted from GITT measurements must be treated with caution, since the relaxation phenomena under these conditions may include multiple contributions such as (i) redistribution of Mg ions due to the formation of metastable (possibly non-stoichiometric solid solution) phases formed under operating conditions (and on application of an overpotential), and (ii) redistribution of phase boundaries to minimise the interfacial energies, etc.

We now consider transport in two cases, the first in a completely stoichiometric material where we need to consider the energy associated with defect formation. In this case, substitutional and vacancy defects are not relevant because they result in stoichiometry changes. The second case explores diffusion in off-stoichiometry materials Mg_3−XBi₂ and Mg₃Bi_2+y, representing Mg vacancy and (excess) Bi substitutional defects, respectively. As no extrinsic Mg vacancies are expected in the stoichiometric Mg₃Bi₂ compound, the diffusion in this case must occur via a vacancy diffusion mechanism involving thermally generated Mg vacancies; hence the observed D would be the diffusion coefficient of Mg vacancies. In the latter case, nonstoichiometry dictates that extrinsic Mg vacancies must exist in the compound.

Table 4 compares the effective activation barriers E^eff_a obtained experimentally versus the E^eff_a estimated from the ab initio calculations using the scalar relativistic (SR) and full relativistic (FR) treatments. Looking at the SR case first, a large variation in E^eff_a is observed where a significantly lower barrier is predicted for the Mg jumps with an existing nearby extrinsic vacancy defect (i.e. case 2, Mg_3−XBi₂ SR): 0.43 eV. In case 1 (stoichiometric Mg₃Bi₂), where generation of a thermal Mg vacancy is required for diffusion, the effective diffusion barrier E^eff_a should now include this formation of vacancy. Naturally, we would expect an increase in E^eff_a due to this inclusion of vacancy formation energy: depending on the vacant site (tetrahedral or octahedral Mg vacancy), the E^eff_a is estimated to be 1.46 eV (octahedral) or 1.55 eV (tetrahedral).

Table 4 Effective migration barriers E^eff_a estimated through DFT and VT NMR techniques. VT NMR data are taken from Liu et al.⁸ For the NMR measurements, Mg₃Bi₂ prepared through electrochemical insertion and mechanical milling were considered. SR and FR refer to scalar relativistic and full relativistic calculations, as described in the text

Case	Diffusion process	E ^effa/eV
		SR	FR
DFT
Case 1	Formation and diffusion of Voct	1.46	0.67
Case 1	Formation and diffusion of Vtet	1.55	0.76
Case 2	Pre-formed vacancy diffusion	0.43	0.34
Experiment
	VT NMR Echem	0.71
	VT NMR ball-mill	0.19

---

We also observe that inclusion of spin–orbit coupling enhances the diffusion noticeably in all three cases, where a reduction of effective activation barrier is seen (e.g. for case 1, V_oct is 1.46 and 0.67 eV for SR and FR, respectively). This is mainly due to the vacancy formation energy, where the FR case results in a dramatically lowered formation energy.

Comparing these results to the experiment, the effective activation barrier E^eff_a in the fully relativistic case ignoring the vacancy creation (case 2) agrees with the E^eff_a of the ball-milled sample, determined through VT NMR. On the other hand, E^eff_a obtained for the electrochemically prepared samples (VT NMR Echem) agree well with the E^eff_a assuming a Mg vacancy creation plus Mg diffusion (case 1). Hence, we conclude that the primary diffusion mechanism in electrochemically prepared samples (measured ex situ, i.e. not during battery operating conditions) should involve vacancy creation, whereas the mechanism in mechanically prepared samples only involves the vacancy diffusion. This could be explained by method of sample preparation: the electrochemical Mg insertion process in this case produced samples closer to the thermodynamic equilibrium creating fewer vacancies, whereas mechanical milling is largely a high energy process resulting in more vacancies (for instance, vacancy formation in ZnO through milling was previously observed with HRTEM studies³⁸). This is also a known phenomenon in the synthesis of intermetallic phases: mechanical milling can provide excess energy to the material, which can be stored in the sample as atomistic disorders of which vacancies are one example.³⁹ Furthermore, as discussed above, Mg deficiency is likely due to preferential Mg sticking to the ball-mill components (jar, balls). The already present vacancies in mechanically prepared samples act as potential diffusion sites for adjacent Mg ions, enhancing their diffusion.

Extending this result to a battery under operating conditions, it is important to stress that the Mg_3−XBi₂ phases formed in situ may not be stoichiometric. The kinetics of Bi (de)magnesiation will depend on a number of factors which include the interfacial energies between the Bi and Mg₃Bi₂ phases, and transport of Mg in Bi, Mg₃Bi₂, and at the various interfaces. In particular, the mechanism of demagnesiation will depend strongly on the ease of vacancy formation in Mg₃Bi₂. By analogy with our work on lithium silicides⁴⁰ and of others (e.g. CuTi₂S₄ (ref. 41)), this energy of vacancy formation may even be responsible for setting (or strongly influencing) the overpotential observed on charge. Rapid ionic transport within a phase also has important implications for the mechanism of the electrochemical oxidation/reduction, because allows the process of Mg₃Bi₂ reduction and dissolution of Mg²⁺ into the electrolyte to be decoupled.⁴⁵ Reduction of Mg₃Bi₂ to Bi requires transport of Mg²⁺ away from the Mg₃Bi₂/Bi interface; Mg²⁺ transport can now happen within the Mg₃Bi₂ phase, allowing Mg²⁺ dissolution to occur anywhere on the surface of the Mg₃Bi₂ particle rather than at a Mg₃Bi₂/Bi/electrolyte phase boundary.

Finally, we note that Zintl-type A₃B₂ materials to which Mg₃Bi₂ belongs have been identified as potentially promising thermoelectric materials.^42,43 Mg mobility in these materials could have important implications on the thermoelectric power generation, and work is in progress on doping other atoms into this structure to enhance, or suppress, this mobility.

4 Conclusion

In conclusion, advanced electronic structure calculations show that spin–orbit coupling plays an important role in structure and dynamics of Mg₃Bi₂, a promising Mg-ion battery anode material. Inclusion of relativistic spin–orbit coupling also significantly lowers the formation energies of Mg vacancy defects, which is crucial to the apparent low Mg-ion diffusion barrier. Using an efficient single-ended hybrid eigenvector-following approach, we have calculated the Mg migration barriers involving relativistic spin–orbit coupling which are as low as 0.34 eV for the octahedral to tetrahedral diffusion. The calculated activation barriers are in good agreement with the previous experimental report using variable temperature ²⁵Mg NMR experiments for materials prepared electrochemically and via ball-milling. Stoichiometric materials show higher activation energies, since the activation energy involves both the cost of vacancy generation and transport. An understanding of Mg transport and the energetics of vacancy formation are important in understanding the mechanisms for demagnesiation of Mg₃Bi₂. This work sheds light on ways of improving the Mg diffusion in similar materials such as Sn, which have been shown to have good capacity but poor rate performance in Mg-ion batteries.⁴⁴

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Professor David Wales for helpful discussions. Via our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). Research was also carried out at the Center for Functional Nanomaterials, Brookhaven National Laboratory, which is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, under Contract No. DE-AC02-98CH10886. J. L. acknowledges Trinity College Cambridge for the graduate studentship. B. M. acknowledges support from the Winton Programme for the Physics of Sustainability, and from Robinson College, Cambridge, and the Cambridge Philosophical Society for a Henslow Research Fellowship.

References

1. R. Van Noorden, Nature, 2014, 507, 26.
2. P. Canepa, G. Sai Gautam, D. C. Hannah, R. Malik, M. Liu, K. G. Gallagher, K. A. Persson and G. Ceder, Chem. Rev., 2017, 117, 4287–4341.
3. R. Mohtadi and F. Mizuno, Beilstein J. Nanotechnol., 2014, 5, 1291–1311.
4. Y. Shao, M. Gu, X. Li, Z. Nie, P. Zuo, G. Li, T. Liu, J. Xiao, Y. Cheng, C. Wang, J.-G. Zhang and J. Liu, Nano Lett., 2014, 14, 255–260.
5. T. S. Arthur, N. Singh and M. Matsui, Electrochem. Commun., 2012, 16, 103–106.
6. F. Murgia, L. Stievano, L. Monconduit and R. Berthelot, J. Mater. Chem. A, 2015, 3, 16478–16485.
7. M. Ramanathan, A. Benmayza, J. Prakash, N. Singh and F. Mizuno, J. Electrochem. Soc., 2016, 163, A477–A487.
8. Z. Liu, J. Lee, G. Xiang, H. F. J. Glass, E. N. Keyzer, S. E. Dutton and C. P. Grey, Chem. Commun., 2017, 53, 743–746.
9. R. Ahuja, A. Blomqvist, P. Larsson, P. Pyykk and P. Zaleski-Ejgierd, Phys. Rev. Lett., 2011, 106, 018301.
10. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
11. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
12. P. E. Blchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
13. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
14. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2006, 124, 219906.
15. J. Hafner, J. Comput. Chem., 2008, 29, 2044–2078.
16. F. Krøger and H. Vink, in Solid State Physics, Elsevier, 1956, vol. 3, pp. 307–435.
17. D. Seymour, S. Chakraborty, D. S. Middlemiss, D. J. Wales and C. P. Grey, Chem. Mater., 2015, 27, 5550–5561.
18. L. J. Munro and D. J. Wales, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 3969–3980.
19. Y. Kumeda, D. J. Wales and L. J. Munro, Chem. Phys. Lett., 2001, 341, 185–194.
20. E. Zintl and E. Husemann, Z. Phys. Chem., 1933, 21, 138–155.
21. S. R. Elliott, The Physics and Chemistry of Solids, Wiley, 1998.
22. V. B. Lazarev, S. F. Marenkin, V. Y. Shevchenko and G. I. Goncharenko, Neorg. Mater., 1974, 10, 1615–1616.
23. A. C. Barnes, C. Guo and W. S. Howells, J. Phys.: Condens. Matter, 1994, 6, L467.
24. W. Howells, A. Barnes and M. Hamilton, Phys. B, 1999, 266, 97–99.
25. R. P. Ferrier and D. J. Herrell, Philos. Mag., 1969, 19, 853–868.
26. R. P. Ferrier and D. J. Herrell, J. Non-Cryst. Solids, 1970, 2, 278–283.
27. L. M. Watson, C. A. W. Marshall and C. P. Cardoso, J. Phys. F Met. Phys., 1984, 14, 113.
28. C. Sutton, Solid State Commun., 1975, 16, 327–330.
29. M. Sedighi, B. Arghavani Nia, H. Zarringhalam and R. Moradian, Eur. Phys. J. Appl. Phys., 2013, 61, 10103.
30. Y. Imai and A. Watanabe, J. Mater. Sci., 2006, 41, 2435–2441.
31. R. Xu, R. A. de Groot and W. van der Lugt, J. Phys.: Condens. Matter, 1993, 5, 7551–7562.
32. X. Zhang, L. Jin, X. Dai and G. Liu, J. Phys. Chem. Lett., 2017, 8, 4814–4819.
33. D. F. Shao, X. Luo, W. J. Lu, L. Hu, X. D. Zhu, W. H. Song, X. B. Zhu and Y. P. Sun, Sci. Rep., 2016, 6, 21484.
34. H.-J. Noh, H. Koh, S.-J. Oh, J.-H. Park, H.-D. Kim, J. D. Rameau, T. Valla, T. E. Kidd, P. D. Johnson, Y. Hu and Q. Li, Europhys. Lett., 2008, 81, 57006.
35. K. B. Yatsimirskii, Theor. Exp. Chem., 1995, 31, 153–168.
36. J. M. Crowley, J. Tahir-Kheli and W. A. Goddard, J. Phys. Chem. Lett., 2016, 7, 1198–1203.
37. R. D. Shannon, Acta Crystallogr. A, 1976, 32, 751–767.
38. D. Chen, Z. Wang, T. Ren, H. Ding, W. Yao, R. Zong and Y. Zhu, J. Phys. Chem. C, 2014, 118, 15300–15307.
39. C. Suryanarayana, in Pergamon Materials Series, Elsevier, 1999, vol. 2, pp. 49–85.
40. K. Ogata, E. Salager, C. Kerr, A. Fraser, C. Ducati, A. Morris, S. Hofmann and C. Grey, Nat. Commun., 2014, 5, 3217–3228.
41. H.-C. Yu, C. Ling, J. Bhattacharya, J. C. Thomas, K. Thornton and A. Van der Ven, Energy Environ. Sci., 2014, 7, 1760–1768.
42. V. Ponnambalam and D. T. Morelli, J. Electron. Mater., 2013, 42, 1307–1312.
43. J. Sun and D. J. Singh, J. Mater. Chem. A, 2017, 5, 8499–8509.
44. N. Singh, T. S. Arthur, C. Ling, M. Matsui and F. Mizuno, Chem. Commun., 2013, 49, 149–151.
45. N. Yamakawa, M. Jiang and C. P. Grey, Chem. Mater., 2009, 21, 3162–3176.

---

Footnote

† Current Address: Max Planck Institute for Chemical Energy Conversion, Stiftstr. 3436, Mülheim an der Ruhr 45470, Germany.

---