Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

An ab initio investigation on the electronic structure, defect energetics, and magnesium kinetics in Mg3Bi2

Jeongjae Lee a, Bartomeu Monserrat b, Ieuan D. Seymour a, Zigeng Liu a, Siân E. Dutton b and Clare P. Grey *a
aDepartment of Chemistry, Lensfield Road, Cambridge CB2 1EW, UK. E-mail: cpg27@cam.ac.uk
bCavendish Laboratory, J. J. Thomson Avenue, Cambridge CB3 0HE, UK

Received 22nd December 2017 , Accepted 25th July 2018

First published on 28th July 2018


Abstract

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 G0W0 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.


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.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 Mg2+, Al3+, and Ca2+ 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.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 Mg2+ 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–6 Bismuth alloys with magnesium to form an intermetallic Mg3Bi2 phase, resulting in a theoretical capacity of 385 mA h g−1. 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 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.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 in-depth ab initio investigation on the Mg3Bi2 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 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 G0W0) to calculate the electronic structure. In particular, the Greens function based G0W0 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 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.

2 Methodology

2.1 Computational details

All Density Functional Theory (DFT) calculations were performed with the VASP code10,11 employing the projector-augmented wave (PAW) method.12 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−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. SCF cycles were converged to a 10−6 eV limit. A Γ-centered k-point sampling of <0.03 Å−1 was used.

Single-shot G0W0 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 G0W0 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 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).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.

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.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, 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

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 Mg2+ and Bi3−. It adopts an anti-La2O3 type structure with hexagonal P[3 with combining macron]m1 space-group symmetry.20 The structure is schematically illustrated in Fig. 1. In terms of ion arrangements, Mg3Bi2 incorporates two alternating layers: layer A consists of tetrahedrally coordinated Mg2+ cations and octahedral interstitial sites forming ‘covalent’ Mg2Bi22− sheets whereas layer B consists of octahedrally coordinated Mg2+ cations and tetrahedral interstitial sites. These interstitial sites are expected to play an important role in ionic diffusion21 as discussed further below.
image file: c7ta11181a-f1.tif
Fig. 1 (a) Crystal structure of Mg3Bi2. 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 Mg3Bi2 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 parameter22 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 Mg3Bi2 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.22
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 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 Mg1.5Bi.23 As expected from the AgI structure, it shows superionic conduction of Mg2+ 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.

3.2 Electronic structure

Despite being first reported in 1933 by Zintl,20 the precise electronic structure and the band gap (Eg) 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;25,26 Lazarev et al. have assumed Mg3Bi2 is semiconducting with Eg = 0.1 eV;22 Watson et al. obtained Eg = 0.5 eV from Mg X-ray emission spectra.27 Work on the amorphous Mg3Bi2 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 Mg3Bi2 is a semiconductor with Eg = 0.25 eV based on the 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–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 G0W0 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 Mg3Bi2 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 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).35 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.9


image file: c7ta11181a-f2.tif
Fig. 2 Density of states plot of Mg3Bi2 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 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.32 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 G0W0 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 G0W0 level of theory.


image file: c7ta11181a-f3.tif
Fig. 3 Band structure plots of Mg3Bi2 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.36 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 Mg3Bi2 sample, formation energies ΔEf 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 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. 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 ΔEf of various stoichiometric and non-stoichiometric defects in Mg3Bi2 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
Δ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
[thin space (1/6-em)]
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
[thin space (1/6-em)]
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 Bi3− 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,37 Bi3− 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 ΔEf 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 VMg(oct) + Mgi(tet) (nn), VMg(tet) + Mgi(tet) (nn), and VMgi(tet) + Mgi(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 ΔEf 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 ΔEf of octahedral vacancies VMg(oct) is found to be as low as 0.33 eV, with a slightly higher value of 0.42 eV for tetrahedral vacancies VMg(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 Mg3Bi2.

3.4 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 VMg(oct). This has lower ΔEf than the tetrahedral vacancy VMg(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 Mgoct and Mgtet (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 Mgoct–Mgoct diffusion barrier (path 2 on Fig. 1b) is around twice that of the Mgoct–Mgtet 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.8 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.


image file: c7ta11181a-f4.tif
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 Voct 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 Mg3Bi2 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) 25Mg 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 Mg3Bi2.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 Mg3−XBi2 and Mg3Bi2+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 4 compares the effective activation barriers Eeffa obtained experimentally versus the Eeffa 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 Eeffa is observed where a significantly lower barrier is predicted for the Mg jumps with an existing nearby extrinsic vacancy defect (i.e. case 2, Mg3−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 Eeffa should now include this formation of vacancy. Naturally, we would expect an increase in Eeffa due to this inclusion of vacancy formation energy: depending on the vacant site (tetrahedral or octahedral Mg vacancy), the Eeffa is estimated to be 1.46 eV (octahedral) or 1.55 eV (tetrahedral).

Table 4 Effective migration barriers Eeffa estimated through DFT and VT NMR techniques. VT NMR data are taken from Liu et al.8 For the NMR measurements, Mg3Bi2 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
[thin space (1/6-em)]
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, Voct 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 Eeffa in the fully relativistic case ignoring the vacancy creation (case 2) agrees with the Eeffa of the ball-milled sample, determined through VT NMR. On the other hand, Eeffa obtained for the electrochemically prepared samples (VT NMR Echem) agree well with the Eeffa 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 studies38). 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 Mg3−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 silicides40 and of others (e.g. CuTi2S4 (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 Mg3Bi2 reduction and dissolution of Mg2+ into the electrolyte to be decoupled.45 Reduction of Mg3Bi2 to Bi requires transport of Mg2+ away from the Mg3Bi2/Bi interface; Mg2+ transport can now happen within the Mg3Bi2 phase, allowing Mg2+ dissolution to occur anywhere on the surface of the Mg3Bi2 particle rather than at a Mg3Bi2/Bi/electrolyte phase boundary.

Finally, we note that Zintl-type A3B2 materials to which Mg3Bi2 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 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.44

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

Footnote

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

This journal is © The Royal Society of Chemistry 2018