An Ab Initio Investigation on the Electronic Structure, Defect Energetics, and Magnesium Kinetics in Mg$_3$Bi$_2$
Jeongjae Lee Bartomeu Monserrat Ieuan Seymour Zigeng Liu Sian Dutton Clare Grey

[NO \title GIVEN]
    [NO \author GIVEN]
    December 30, 2023
======================




We present a comprehensive ab initio investigation on Mg3Bi2, 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 Mg3Bi2 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.<cit.> 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 important for large-scale applications such as electric vehicles and grid-based energy storage systems.<cit.> 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.<cit.>

While possible high-voltage MIB cathode candidates have been identified,<cit.> 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.<cit.> In this respect, bismuth metal is a promising anode material with a low discharge voltage of 0.2 V versus Mg metal.<cit.> Bismuth alloys with magnesium to form an intermetallic Mg3Bi2 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.<cit.> Recently, our group has used ^25Mg Nuclear Magnetic Resonance (NMR) spectroscopy to study Mg3Bi2 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.<cit.>

 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 Mg3Bi2 system, the end product of discharge process. The results are presented in four sections: in the first section, we discuss the structural characteristics of Mg3Bi2, 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_0W_0) to calculate the electronic structure. In particular, the Greens function based G_0W_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.<cit.> 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 Mg3Bi2, 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.




§ METHODOLOGY



 §.§ Computational Details

All Density Functional Theory (DFT) calculations were performed with the VASP code<cit.> employing the projector-augmented wave (PAW) method<cit.>. Spin-polarized Perdew-Burke-Ernzerhof (PBE) and Heyd-Scuseria-Ernzerhof (HSE06) exchange-correlation functionals were adopted.<cit.> 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_0W_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_0W_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;<cit.> and a `full relativistic' correction that explicitly includes the spin-orbit coupling term alongside the scalar relativistic correction.

The initial structure of Mg3Bi2 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. BiMg(oct) denotes a neutral Bi sitting on an octahedral Mg site).<cit.> 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.<cit.> Details on the algorithm are available elsewhere.<cit.> 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.



§ RESULTS AND DISCUSSION




 §.§ Crystal structure

Mg3Bi2 is the last in a series of magnesium pnictide compounds Mg3X2 (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-La2O3 type structure with hexagonal P3̅m1 space-group symmetry.<cit.> The structure is schematically illustrated in Figure <ref>. In terms of ion arrangements, Mg3Bi2 incorporates two alternating layers: layer A consists of tetrahedrally coordinated  Mg^2+ cations and octahedral interstitial sites forming `covalent' Mg2Bi2^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<cit.> as discussed further below.
 
 
 
From DFT-based lattice relaxation, we can see that the PBE functional overestimates the experimentally determined a cell parameter<cit.> 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 <ref>). 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 Mg3Bi2 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.<cit.> As expected from the AgI structure, it shows superionic conduction of Mg^2+ cations as determined through neutron diffraction.<cit.> 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<cit.>, the precise electronic structure and the band gap (E_g) of Mg3Bi2 has not been determined with either experiment or theoretical calculations. In a series of work, Ferrier and co-workers have speculated that Mg3Bi2 is semimetallic based on a conductivity-composition plot;<cit.> Lazarev et al. have assumed Mg3Bi2 is semiconducting with E_g=0.1 eV;<cit.> Watson et al. obtained E_g=0.5 eV from Mg X-ray emission spectra.<cit.> Work on the amorphous Mg3Bi2 alloy have shown from conductivity measurements that it is semiconducting with a band gap of 0.15 eV.<cit.> 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 Mg3Bi2 is a semiconductor with E_g=0.25 eV based on Engel-Vosko Generalized Gradient Approximation (EV-GGA),<cit.> whereas Imai et al., Xu et al., and Zhang et al. concluded it is a semimetal based on pure GGA calculations.<cit.> 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.<cit.>; 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_0W_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 <ref>. 
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. 






Figure <ref> shows the density of states (DOS) plots obtained for Mg3Bi2 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 Mg3Bi2, 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).<cit.>. The downshift in 6s energy by adding the full relativistic effect is around -0.15 eV; effects of similar magnitudes were observed in PbO2, an important active material in lead-acid batteries.<cit.>

Figure <ref> shows the band structure of Mg3Bi2 using DFT within the semilocal PBE and hybrid HSE06 approximations, with and without the spin-orbit interaction. At the PBE level, Mg3Bi2 exhibits semimetallic behavior, consistent with earlier reports of a type-II nodal line semimetal at this level of theory.<cit.> The inclusion of spin-orbit coupling opens a small gap in the nodal line, but the system remains semimetallic. At the HSE06 level, Mg3Bi2 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_0W_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 Mg3Bi2 is also semiconducting at the G_0W_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.<cit.> 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 Mg3Bi2 sample, formation energies Δ E_f of various stoichiometric and non-stoichiometric point defects were calculated ab initio. The results outlined in Table <ref> 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 Mg3Bi2 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 Mg23Bi16 (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 Mg25Bi15 (again for a 40-atom supercell) and five excess electrons on the Fermi level. Charge neutrality is automatically maintained in the antisite and Frenkel defects since they are stoichiometric.

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^3- 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 Mg3Bi2 crystal: Mgoct–Bi 3.21 Å versus Mgtet–Bi 2.92 Å. Using the Shannon radii of 0.72 and 0.57 Å for Mgoct and Mgtet, respectively,<cit.> Bi^3- 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_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 Mg3Bi2.


		


 §.§ Mg migration kinetics

Having established the defect chemistry in Mg3Bi2, 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 Figure <ref> clearly shows a small diffusion barrier of hopping between the Mg_oct and Mg_tet (Path 1 on Figure <ref>b), 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 <ref>b) 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 ^25Mg NMR studies of Liu et al.<cit.> An alternative exchange mechanism involving Path 3 on Figure <ref>b 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 Mg3Bi2 in a Mg-ion battery.<cit.> 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) ^25Mg NMR spectroscopy to probe Mg transport.<cit.> 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 Mg3Bi2.<cit.> 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-xBi2 and Mg3Bi_2+y, representing Mg vacancy and (excess) Bi substitutional defects, respectively. As no extrinsic Mg vacancies are expected in the stoichiometric Mg3Bi2 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 <ref> compares the effective activation barriers E_a^eff obtained experimentally versus the E_a^eff 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_a^eff 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-xBi2 SR): 0.43 eV. In Case 1 (stoichiometric Mg3Bi2), where generation of a thermal Mg vacancy is required for diffusion, the effective diffusion barrier E_a^eff should now include this formation of vacancy. Naturally, we would expect an increase in E_a^eff due to this inclusion of vacancy formation energy: depending on the vacant site (tetrahedral or octahedral Mg vacancy), the E_a^eff 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 FR case results in a dramatically lowered formation energy.

Comparing these results to the experiment, the effective activation barrier E_a^eff in the fully relativistic case ignoring the vacancy creation (Case 2) agrees with the E_a^eff of the ball-milled sample, determined through VT NMR. On the other hand, E_a^eff obtained for the electrochemically prepared samples (VT NMR Echem) agree well with the E_a^eff 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<cit.>). 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.<cit.> 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-xBi2 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 Mg3Bi2 phases, and transport of Mg in Bi, Mg3Bi2, and at the various interfaces. In particular, the mechanism of demagnesiation will depend strongly on the ease of vacancy formation in Mg3Bi2. By analogy with our work on lithium silicides<cit.> and of others (e.g. CuTi2S4<cit.>), 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_3B_2 materials to which Mg3Bi2 belongs have been identified as potentially promising thermoelectric materials.<cit.> 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 Mg3Bi2, 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 ^25Mg 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 Mg3Bi2. 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.<cit.>




§ ACKNOWLEDGEMENT

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. 



§ SUPPORTING INFORMATION

Research data in support of this manuscript will be made available on Apollo, the University of Cambridge repository.

rsc


