An Ab Initio Investigation on the Electronic Structure, Defect Energetics, and Magnesium Kinetics in Mg$_3$Bi$_2$

We present a comprehensive ab initio investigation on Mg$_3$Bi$_2$, a promising Mg-ion battery anode material with high rate capacity. Through combined DFT (PBE, HSE06) and $G_0W_0$ electronic structure calculations, we find that Mg$_3$Bi$_2$ 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.


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. 1 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 2+ , Al 3+ , and Ca 2+ have gained attention due to their low cost and high volumetric densities, features 1 arXiv:1807.10856v1 [cond-mat.mtrl-sci] 27 Jul 2018 important for large-scale applications such as electric vehicles and grid-based energy storage systems. 2 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 2+ 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. 3 While possible high-voltage MIB cathode candidates have been identified, 2 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. 3 In this respect, bismuth metal is a promising anode material with a low discharge voltage of 0.2 V versus Mg metal. [4][5][6] Bismuth alloys with magnesium to form an intermetallic Mg 3 Bi 2 phase, resulting in a theoretical capacity of 385 mAh/g. However, the most interesting feature of 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 25 Mg Nuclear Magnetic Resonance (NMR) spectroscopy to study Mg 3 Bi 2 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. 8 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 indepth ab initio investigation on the Mg 3 Bi 2 system, the end product of discharge process. The results are presented in four sections: in the first section, we discuss the structural characteristics of Mg 3 Bi 2 , 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 0 W 0 ) to calculate the electronic structure. In particular, the Greens function based G 0 W 0 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. 9 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 3 Bi 2 , 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.

Computational Details
All Density Functional Theory (DFT) calculations were performed with the VASP code 10,11 employing the projector-augmented wave (PAW) method. 12 Spin-polarized Perdew-Burke-Ernzerhof 2 (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. In addition, 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 −4 eV. Monkhorst-Pack k-point sampling of < 0.05Å -1 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. Cellular relaxations SCF cycles were converged to a 10 −6 eV limit. Γ-centered k-point sampling of < 0.03Å -1 was used.
Single-shot G 0 W 0 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 0 W 0 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 −8 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; 15 and a 'full relativistic' correction that explicitly includes the spin-orbit coupling term alongside the scalar relativistic correction.
The initial structure of Mg 3 Bi 2 was fully relaxed until the energy differences between the subsequent steps are converged to 10 −5 eV per cell and the forces are < 0.05 eV/Å -1 . 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). 16 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.

Transition State Searching
Transition state and steepest-descent pathway were found using the hybrid eigenvector-following approach as implemented in the OPTIM code. 17 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Å -1 . Then up to five LBFGS minimization steps are performed in the tangent space until the root-mean-square gradient is smaller than 10 −3 eVÅ -1 . 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 −3 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, HEF method can sample the path at an arbitrary step size. This allows us to capture the finer details of the energy variance along the pathway.

Crystal structure
Mg 3 Bi 2 is the last in a series of magnesium pnictide compounds Mg 3 X 2 (X=P, As, Sb, Bi) and is a Zintl-type compound with closed shell ions with formal charges of Mg 2+ and Bi 3 -. It adopts an anti-La 2 O 3 type structure with hexagonal P3m1 space-group symmetry. 20 The structure is schematically illustrated in Figure 1. In terms of ion arrangements, Mg 3 Bi 2 incorporates two alternating layers: layer A consists of tetrahedrally coordinated Mg 2+ cations and octahedral interstitial sites forming 'covalent' Mg 2 Bi 2 -2 sheets whereas layer B consists of octahedrally coordinated Mg 2+ cations and tetrahedral interstitial sites. These interstitial sites are expected to play an important role in ionic diffusion 21 as discussed further below.  From DFT-based lattice relaxation, we can see that the PBE functional overestimates the experimentally determined a cell parameter 22 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). Finally, we note that Mg 3 Bi 2 undergoes a phase transition above 703 • C to a defective bodycentered cubic structure similar to AgI, with excess Mg cations as the tetrahedral interstitials, resulting in Mg 1.5 Bi. 23 As expected from the AgI structure, it shows superionic conduction of Mg 2+ cations as determined through neutron diffraction. 24 In this study, however, we restrict the investigation to the room-temperature hexagonal phase which is more relevant to Mg-ion batteries.

Electronic structure
Despite being first reported in 1933 by Zintl, 20 the precise electronic structure and the band gap (E g ) of Mg 3 Bi 2 has not been determined with either experiment or theoretical calculations. In a series of work, Ferrier and co-workers have speculated that Mg 3 Bi 2 is semimetallic based on a conductivity-composition plot; 25, 26 Lazarev et al. have assumed Mg 3 Bi 2 is semiconducting with E g = 0.1 eV; 22 Watson et al. obtained E g = 0.5 eV from Mg X-ray emission spectra. 27 Work on the amorphous Mg 3 Bi 2 alloy have shown from conductivity measurements that it is semiconducting with a band gap of 0.15 eV. 28 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 3 Bi 2 is a semiconductor with E g = 0.25 eV based on Engel-Vosko Generalized Gradient Approximation (EV-GGA), 29 whereas Imai et al., Xu et al., and Zhang et al. concluded it is a semimetal based on pure GGA calculations. [30][31][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 0 W 0 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 3 Bi 2 using the Perdew-Burke-Ernzerhof (PBE) and Heyd-Scuseria-Ernzerhof (HSE06) exchange-correlation functionals. Only the valence charge is calculated. Figure 2 shows the density of states (DOS) plots obtained for Mg 3 Bi 2 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 3 Bi 2 , 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). 35 The downshift in 6s energy by adding the full relativistic effect is around -0.15 eV; effects of similar magnitudes were observed in PbO 2 , an important active material in lead-acid batteries. 9 Figure 3 shows the band structure of Mg 3 Bi 2 using DFT within the semilocal PBE and hybrid HSE06 approximations, with and without the spin-orbit interaction. At the PBE level, Mg 3 Bi 2 exhibits semimetallic behavior, consistent with earlier reports of a type-II nodal line semimetal at this level of theory. 32 The inclusion of spin-orbit coupling opens a small gap in the nodal line, but the system remains semimetallic. At the HSE06 level, Mg 3 Bi 2 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 0 W 0 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 3 Bi 2 is also semiconducting at the G 0 W 0 level of theory.
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. 36 The net effect of this splitting is a reduction in the band gap.

Defect energetics
To investigate the possibility of different types of defects present in the Mg 3 Bi 2 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  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 3ion 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 3 Bi 2 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, 37 Bi 3radii 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 Mg(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 Mg-ion diffusion.
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, ∆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 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 3 Bi 2 .

Mg migration kinetics
Having established the defect chemistry in Mg 3 Bi 2 , 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.  Table 3: Scalar relativistic (SR) and full relativistic (FR) ab initio formation energies ∆E f of various stoichiometric and non-stoichiometric defects in Mg 3 Bi 2 using the PBE functional. Defect notations follow the convention of Kroger and Vink with neutral sign omitted for clarity. 16 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.
The result presented in Figure 4 clearly shows a small diffusion barrier of hopping between the Mg oct and Mg tet (Path 1 on Figure 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 Figure 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 25 Mg NMR studies of Liu et al. 8 An alternative exchange mechanism involving Path 3 on Figure  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.
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 3 Bi 2 in a Mg-ion battery. 4 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) 25 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 3 Bi 2 . 8 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-x Bi 2 and Mg 3 Bi 2+y , representing Mg vacancy and (excess) Bi substitutional defects, respectively. As no extrinsic Mg vacancies are expected in the stoichiometric Mg 3 Bi 2 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-x Bi 2 SR): 0.43 eV. In Case 1 (stoichiometric Mg 3 Bi 2 ), 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). 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, SR with V oct 1.46 eV versus FR V oct 0.67 eV. This is mainly due to the vacancy formation energy, where the    8 For the NMR measurements, Mg 3 Bi 2 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.
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 38 ). 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. 39 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-x Bi 2 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 3 Bi 2 phases, and transport of Mg in Bi, Mg 3 Bi 2 , and at the various interfaces. In particular, the mechanism of demagnesiation will depend strongly on the ease of vacancy formation in Mg 3 Bi 2 . By analogy with our work on lithium silicides 40 and of others (e.g. CuTi 2 S 4 41 ), this energy of vacancy formation may even be responsible for setting (or strongly influencing) the overpotential observed on charge.
Finally, we note that Zintl-type A 3 B 2 materials to which Mg 3 Bi 2 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.

Conclusion
In conclusion, advanced electronic structure calculations show that spin-orbit coupling plays an important role in structure and dynamics of Mg 3 Bi 2 , 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 25 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 3 Bi 2 . 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. 44