Understanding Performance Limitations to Enable High Performance Magnesium-Ion Batteries

A mathematical model was developed to investigate the performance limiting factors of Mg-ion battery with a Chevrel phase (Mg x Mo 6 S 8 ) cathode and a Mg metal anode. The model was validated using experimental data from the literature [Cheng et al., Chem. Mater. , 26 , 4904 (2014)]. Two electrochemical reactions of the Chevrel phase with signiﬁcantly different kinetics and solid diffusion were included in the porous electrode model, which captured the physics sufﬁciently well to generate charge curves of ﬁve rates (0.1C–2C) for two different particle sizes. Limitation analysis indicated that the solid diffusion and kinetics in the higher-voltage plateau limit the capacity and increase the overpotential in the Cheng et al.’s thin (20- μ m) electrodes. The model reveals that the performance of the cells with reasonable thickness would also be subject to electrolyte-phase limitations. The simulation also suggested that the polarization losses on discharge will be lower than that on charge, because of the differences in the kinetics and solid diffusion between the two reactions of the Chevrel phase.

During the past 30 years, Li-ion batteries have come to dominate the market for use in electric vehicles (EVs) and portable electronic devices. However, to extend the driving range of EVs and to make the portable electronic devices more compact and lighter, the energy density of the batteries must be increased. Significant effort has been devoted to finding new battery chemistries with the potential to exceed the energy density of Li-ion batteries. [1][2][3][4] Among the novel approaches being studied, Mg-ion batteries are promising candidates, [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] partly because of the advantages of Mg metal anode. Although the reduction potential of Mg metal is similar to that of Li (-2.37 V vs. H + /H 2 for Mg compared with -3.0 V vs. H + /H 2 for Li), the volumetric capacity is significantly higher (3,833 mAh cm -3 for Mg compared with 2,046 mAh cm -3 for Li). In addition, Mg-ion batteries may be safer than Li-ion batteries, not only because Mg does not form hazardous dendrites when charged and discharged repeatedly 25 but also because a strong and stable oxide develops on the surface, protecting the battery if the metal is exposed to ambient conditions. Finally, the divalent charge on the Mg ion enables exploration of higher-capacity cathode materials that promise increases in the overall cell-energy density.
Although multivalent systems are promising from an energydensity perspective, they are limited by the slow transport of the divalent cation. This issue manifests itself in the electrolyte phase, where the Mg ion is thought to form ion pairs, 26,27 and in typical cathode materials, where the diffusion of ions can be fairly slow. 28,29 Furthermore, although not clearly explored in the literature, the sluggish movement of Mg ions can also manifest itself at the interface during charge transfer. To date, it is unclear which of these limitations (bulk cathode vs. bulk electrolyte vs. interfaces) limits the rate capability of Mg-ion systems most.
In 2000, Aurbach et al. introduced the first generation of rechargeable Mg-ion batteries 5 using Mg metal as the anode, Chevrel phase (Mo 6 S 8 ) as the cathode, and dichloro complex in tetrahydrofuran (DCC/THF) as the electrolyte. Two plateaus in the discharge and charge curves indicated that the Chevrel phase underwent dualphase transformation. Aurbach et al. demonstrated that the system has good cyclability, e.g., >2,500 cycles. Although the Mg-metal anode had a high specific capacity (2,205 mAh g −1 ), the low capacity of the Chevrel phase (122 mAh g −1 ) limited the overall cell capacity.
Since the first efforts of Aurbach et al., much research has been devoted to developing new cathode materials, 3,9,20,22,24,29 such as vana-dium oxide (V 2 O 5 ) 17 and manganese dioxide (MnO 2 ). 14 However, so far, no systematic study has been undertaken to understand the limitations in the systems imposed by the divalent nature of the working ion. The purpose of this paper is to understand the exact cause of the limitation in the Mg-ion system. We used the Chevrel-phase cathode as an example system to explore different limitations including the polarization loss from the bulk cathode, bulk electrolyte, and the two interfaces. We performed this analysis by developing a continuum model for the chemistry using the macrohomogeneous formulation pioneered by Newman and coworkers [30][31][32][33][34][35] and comparing it to experimental data obtained from a literature. 21 Reported literature values were used for many of the parameters for the various phenomena, and the model was compared with the data to ascertain the limitations in the system. The Appendix provides the detailed equations and properties used in this study. To the best of our knowledge, this paper is the first to examine the limitations in this promising next-generation battery system.
Our results indicate that although the Mg metal reaction can be reasonably fast in an all-phenyl-complex (APC) electrolyte, the system is limited by the polarization losses in the liquid electrolyte and both the bulk diffusion and interfacial charge-transfer resistance in the Chevrel phase during the reaction at the upper plateau. Finally, we show that the Chevrel system exhibits unique behavior, wherein the discharge rate capability can be higher than the charging rate capability due to the faster transport and charge-transfer process in the lower Mg insertion plateau.

Rate Behavior of Chevrel Phase
In this paper we used experimental data on the Chevrel system reported by Cheng et al. 21 Cheng et al. provide information on the rate capability for materials using two different particle sizes in APC/THF electrolyte. Although the "conventional" material was reported to have a wide particle size distribution ranging from 200 nm to 4 μm, the authors also prepared nanosized particles of uniform size (approximately 250 nm) and performed cycling experiments. Scanning electron microscopy (SEM) was used to verify the sizes and uniformity of the particles. 21 By using the APC in THF electrolyte instead of the DCC electrolyte, the electrochemical window could be widen up to 3.3 V. 20,36 Figure 1 (reproduced from Reference 21) presents the charge and discharge curves of two cells containing the Chevrel phase with different particle sizes. Figure 1a shows the data on the conventional material, and Figure 1b presents the data for the material with the smaller and uniform particles. Both figures show the typical decrease in capacity with increasing rate. However, the data sets illustrate the effect of particle size on the performance of the cathode. The smallerparticle-size materials exhibit overall higher capacity and better capacity utilization at high rates. In addition, the concentration polarization is smaller in the nanosized material. Figure 1 also shows that there is asymmetry between the charge and discharge with higher polarization losses on charge compared to discharge. The data in Figure 1 was used for model development and validation.
Note that Cheng et al. charged and discharged the cells from low to high C-rates using a constant current mode at fixed cutoff voltages. 21,37 Because the typical constant-current-constant-voltage (CC-CV) mode was not implemented, it is not clear that the cells were all fully charged before the discharge experiment was implemented. Therefore, it is possible that the decrease in discharge capacity with rate is simply a consequence of not charging the battery to its full capacity. To test this hypothesis, in Figure 2, we plot the ratio between the discharge and charge capacities obtained from Figures 1a and 1b. Figure 2 demonstrates that the ratio between the discharge and charge capacities is always close to unity even at high C-rates. The discharge capacity is almost the same as the charge capacity obtained on the previous charge. This suggests that the limitations observed in Figure 1 are all from the inability of the battery to be charged, rather than its inability to be discharged. In other words, if these cells were charged using CC-CV protocols, the capacity at high C-rates could be even higher. Because of this interesting behavior, in this article, we used the charging experiments to ascertain the limitations in the system. Figure 2 also suggests that these cells exhibit a small, but finite, side reactions; evident from the ratio being less than one. In addition, the ratio increases with increasing current density, a classic indication of side reactions wherein decreasing the time for reaction decreases the inefficiency. Levi et al. performed ex-situ X-ray diffraction studies of Mg 2 Mo 6 S 8 and observed that the following two reactions predominate at the high-and low-voltage plateaus observed in Figure 1: 28,29 Mg 2+ + 2e − + Mo 6 S 8 → MgMo 6 S 8 slow at high voltage; [1] Mg 2+ +2e − +MgMo 6 S 8 → Mg 2 Mo 6 S 8 fast at low voltage. [2] By executing slow-scan rate cyclic voltammetry, Levi et al. also demonstrated that the first reaction (Eq. 1) is much slower than the second (Eq. 2). The reaction kinetics differ because of the characteristic arrangement of the cation sites in the crystal structure of Mg 2 Mo 6 S 8 . The first stage of Mg insertion (Eq. 1) describes occupying the inner lattice sites. In the second reaction (Eq. 2), one Mg 2+ ion per formula unit is located at the inner sites and another is located at the outer site. The minimal distance between the inner sites, 0.9 Å, is only half of that between the inner and outer sites, 2.0 Å. Therefore, the two reactions exhibit different activation-energy barriers. 28,29 The characteristic structure of Mg 2 Mo 6 S 8 also appears to affect the ion mobility in the solid. In 2005, Levi et al. used the classical potentiostatic intermittent titration technique (PITT) to measure the solid diffusivity of Mg x Mo 6 S 8 as a function of Mg content, x. 38 Levi et al. reported that measurement of the diffusivity represented in Eq. 1 was not possible because the value was too low. The diffusivity of Eq. 2 was measured to be 2.0 × 10 -17 m 2 s -1 , which is one to two orders lower than the solid diffusivity of Li + in LiCoO 2 . 39 This result indicates that not only the surface kinetics but also the ion diffusion differ greatly in the two stages of Mg insertion.

Mathematical Model, Assumptions, and Material Properties
We applied the well-established macrohomogeneous battery simulation model to analyze the limitations of the Mg/Chevrel phase battery. Figure 3 shows the schematic of the Mg/Chevrel phase  21 it has Mg metal as the anode, Chevrel phase (Mo 6 S 8 ) as the cathode, and APC/THF as the electrolyte. While the Chevrel phase is a porous electrode, the metal anode is a planar interface. All the experiments, and therefore the model, are for room temperature operation. Thermal analysis was not addressed in the paper because of lack of extensive experimental data at different temperatures. This paper focused on the limitation analysis at room temperature only.
Newman and co-workers developed a model for the Li-ion battery cell based on the porous-electrode theory, the concentrated-solution theory, the Ohm's law, and the intercalation kinetics 30-35 that incor-porates three core physics concepts within the battery system: lithium transport in the electrolyte and solid, interfacial reaction kinetics, and reaction thermodynamics. The Appendix provides the detailed equations used in this study. Table I lists the material properties used in our model. When we were unable to find a material property in the literature, we used the most acceptable value of a similar chemistry. Specifically, for the electrolyte property, Mizrahi et al. reported the ionic conductivity of APC in THF at various temperatures. 36 Therefore, we could use the value for the electrolyte property at room temperature directly in our modeling. No values for the diffusion coefficient or transference number of electrolytes for APC in THF were found in the literature; we therefore utilized the data for a magnesium organohaloaluminate electrolyte (C 2 H 5 MgCl−((C 2 H 5 ) 2 AlCl) 2 /THF) reported by Benmayza et al. 40 We note that the data shown in Figure 1 were obtained on electrodes that were 20-μm thick, 21 much thinner than actual cells'. Under these conditions the electrolyte limitations are expected to be small in the cathode, therefore, we believe that the simulations that predict the experimental data are not affected by these assumptions. We examine the importance of the electrolyte at the end of the paper.
Special care was taken to describe the two reactions in the cathode. First, we assume that the bulk reaction in the cathode proceeded via intercalation. While the flat potential would suggest a phase change reaction (a consequence of the Gibbs phase rule), no information is available in the literature on the juxtaposition of the different phases. Because this paper uses information from Levi et al. where a similar assumption was made, the diffusion  coefficient and the exchange-current density used here should be viewed as average properties that capture both the transport phenomena and the phase movement.
As described above, both the kinetics and diffusivity in the solid differed greatly in the two reactions. To capture these experimental observations, we developed the following equations for the solid diffusivity and exchange-current density of the cathode: where s D and s i are smoothing factors that control the extent to which the step shape is smoothed, and y is the depth of discharge (DOD). The smoothing factors were fitting parameters as well; therefore, they were fitted using the voltage curves in Figure 1. The smoothing factors of the solid diffusivity, s D , and the exchange-current density, s i , do not have the theoretical basis to be the same; however, it turned out that a same value of 65 for both of s D and s i fit the voltage curves adequately and was used in the paper. The solid diffusivity and the exchange-current density values are presented in Figures 4a and 4b.
For the solid diffusivity of Mg x Mo 6 S 8 (Eq. 3), we were able to adopt the value reported by Levi et al. 38 for the fast reaction, D s,fast . Levi et al. measured D s,fast using the PITT method; however, as mentioned above, they could not measure the lower diffusivity, D s,slow , because the value was too small, caused by Mg ions trapping. 28,29,38 For the exchangecurrent density, i 0 , no reported value of Mg x Mo 6 S 8 in APC/THF (Eq. 4) was found in the literature; therefore, we used i 0,fast , i 0,slow , and D s,slow as fitting parameters in the model.

Results and Discussion
Fitting results and extrapolation to larger particle sizes.-As noted above, Cheng et al. used CC-only protocol, without a CV step to fully-charge the battery. As shown in Figure 2, their data suggests that the drop in capacity at higher rates originates from limitations on charge and that the subsequent discharge results in complete insertion of Mg in the lattice. In other words, while the Mg content at the beginning of charge varies from rate to rate, the Mg content at the discharge end (x ≈ 2) is fairly consistent at all rates. Therefore, simulations were conducted on charge in order to ensure that the initial condition was consistent at different C-rates.
Figures 5a and 5b shows model-experiment comparisons for the two particle-sizes across five C-rates, from 0.1 C to 2 C. See Appendix for the model details. The figures show that the model was capable of predicting the rate behavior fairly accurately. Figure 5a shows our methodology and success in fitting the model with the three parameters (the solid diffusivity for the process in Eq. 1, the exchange-current densities of cathode for the processes in Eqs. 1 and 2) and two smoothing factors of the solid diffusivity and the exchange-current density. A systematic approach was taken to extract the five unknown parameters to ensure accuracy of the prediction. Previous modeling experience provides a pathway to extract the unknown parameters from the data, by comparing it to the mathematical model. During charging, lower exchange-current density moves the voltage higher without impacting the overall capacity, unless the voltage hits the cutoff voltage before complete insertion. The solid diffusivity, on the other hand, has a direct effect on the utilization while also increasing the voltage. Therefore, we first adjust the diffusivity to fit the charge capacity and then adjust the exchange-current density until the charge voltage is similar to that in the experiment. We perform this fit at one C-rate and predict the behavior at other rates without changing the parameters. Finally, we adjust the smoothing factors so that the voltage increase from one plateau to another has a similar qualitative shape to the experiment. A small smoothing factor results in a sharp change between the two plateaus while a larger factor results in a gradual change. We consider the smoothing function to be minimally important for the results.
This approach shows a remarkable ability to predict the experimental data. Despite the large deviations in capacity and voltage drop  ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 207.241.231.82 Downloaded on 2018-07-20 to IP between 0.1 C and 2 C, the model was able to predict the behavior while capturing characteristics such as the decreasing capacity of the higher voltage plateau at higher C-rates. The simulation results match the total capacities and operating voltages at five different C-rates as well as it captures how the capacities and voltages of high and low plateau change at different rates.
In addition, the model was able to successfully extrapolate to predict behavior for cells with larger particles, as shown in Figure  5b. It should be emphasized that this prediction was achieved with no other changes to the material properties in the model. Figure 5b shows that the model still matches the most of the total capacities and operating voltages well. This was possible because we implemented different transport and kinetic characteristics of two reactions, shown in Eq. 1 and Eq. 2. Note that the SEM image obtained by Cheng et al. 21 reveals that the traditional cathode material exhibits highly non-uniform particle-sizes. Without introducing the complexity of multiple particle sizes, we use single-average size of d = 600 nm and find that the model provides a very reasonable prediction. While the error at high C-rates for the larger particle appears large, we believe that these differences are due to the use of single-particle size. At higher C-rates, smaller particles can be discharged fully, providing more charging capacity than what would be predicted by a singleparticle size. While inclusion of multiple sizes would provide better predictions, the early qualitative understanding of the limitations in Mg-ion batteries is easily achieved without this added complexity. This feature will be the subject of future publications.

Limitations of Mg-ion batteries.-In this section, we implemented
the developed model to analyze the limitations of these battery systems. By manipulating the system properties in the model, we could examine why a given battery does not reach its theoretical capacity or open-circuit voltage (OCV) when operated at a specific C-rate. 41,42 In general, higher C-rates decrease capacity and increase overpotentials. Our analysis provided qualitative information about which factors would have the greatest effects on battery performance. The information provides the direction to focus for better Mg-ion system developments. Figure 6 shows the limitations of the Mg/Mo 6 S 8 systems. The bullets in Figure 5 represent the experimental data obtained by Cheng et al. 21 for small particles at 0.5C. The blue line shows the simulation results. The black line is the OCV. The gap between the blue and black lines represents the additional energy required for charging and the capacity loss from overpotentials.
The lower voltage plateau in Eq. 1, where the reaction is more facile, shows almost no loss of capacity; there is, however, an approximately 13-mV overpotential in the electrolyte phase, as indicated by the gap between the blue and black line at ∼1.1 V. In the slow reaction area (higher voltage plateau in Eq. 2), the capacity loss was greater. The overpotentials generated by the solid diffusion and the kinetic overpotential of the cathode were much greater in the slow reaction.
The fast kinetics at the anode, however, results in negligible surface overpotential.
Discharge simulations.- Figures 1 and 2 demonstrate that the charge and discharge capacities of the Mg/Mo 6 S 8 system remain similar even at high C-rates. Although this system exhibits good discharge characteristics, the capacity might be increased even further by operating in a CC-CV mode. The modeling results support this conclusion. We performed one hypothetical simulation assuming the initial depth of discharge (DOD) to be zero. All the material properties remained the same as previous simulations; however, current sign was reversed. Figure 7 presents the modeling results, which indicate that the system would approach its theoretical capacity at a low C-rate as well as even at a high C-rate (2C), which agrees with the experimental observations presented in Figures 1a and 1b. The capacity will be close to the theoretical value regardless of the discharge rate. Figure 8 provides a possible explanation for the large discharging capacities, shown in Figure 7. Figure 8a demonstrates what would happen in the cathode particles while charging. The particles were filled with Mg ions when the charge starts. The cell was first fully discharged, y = 1; therefore, both the surface kinetics and solid diffusion are fast (as shown in Figure 4 and Eqs. 3 and 4). Because of the fast kinetics and solid diffusion, the normalized concentration uniformly decreased to y = 0.5 in the entire particle (0 ≤ r/r p ≤ 1). However, after this stage, the kinetics as well as the diffusion become slow; therefore, the surface will be depleted even though cyclable Mg ions remain inside the particles. This effect indicates that it would be extremely difficult for the Mg ions inside the particle, close to the core, to leave the particle. This condition would increase the surface overpotential and decrease the capacity. The OCV is determined by the surface concentration; therefore, the charging voltage will reach the cutoff voltage if the surface is depleted. The Mg ions, however, are still in the particles, which explain why the charge capacity is much less than the theoretical capacity. Figure 8b shows the discharge case. The situation is reversed to that of the charging case. When the discharge begins, the cathode particle is almost empty (y = 0), and the kinetics and solid diffusion are slow (as described in Figure 4 and Eqs. 3 and 4). With time, however, the normalized concentration exceeds 0.5, and then, both the kinetics and diffusion become fast; therefore, the entire particle will be filled with Mg ions. When the DOD reaches y = 0.5, the kinetics speed up, and the Mg ions are delivered rapidly to the inside. Therefore, the surface concentration does not increase rapidly, and the capacity will be near the theoretical value.  Effects of realistic cell thickness.-The data by Cheng et al. was obtained using cathode electrodes that were 20-μm thick, 21 much thinner than practical thickness used in actual batteries. For thin electrodes, the APC/THF electrolyte might appear to be an ideal electrolyte, based on Figure 6, which shows only minor limitations. A simulation performed on a hypothetical cell, however, revealed that such an electrolyte could also be a source of significant limitations. We performed simulations by increasing the thickness to 80 μm (typical for high energy Li-ion cells) to provide a fair comparison to Li systems currently in use. Figure 9 demonstrates that the electrolyte phase exhibited more pronounced limitations at the increased (but realistic) thickness of 80 μm when the other characteristics remained the same. The increased thickness of the cathode increased the overpotential in the liquid phase and decreased the capacity significantly. These effects result from the ionic resistance and concentration polarization in the porous electrode. More research in electrolytes with high transport properties are needed to ensure that these limitations are eliminated when operating cells under realistic designs.

Conclusions
In this study, we developed a mathematical model that matches the charge curves of Mg/Mo 6 S 8 battery systems to identify the limitations of this system. Although models specific to multi-phase systems (such as core-shell 42 or phase-field model [43][44][45] ) may have been better modeling tools for this system, classical Dualfoil modeling [30][31][32][33][34][35] was sufficiently accurate to capture the core physics: the results fit the experimental data obtained by Cheng et al. 21 for five C-rates and two particle sizes. The different kinetic and transport properties in the two plateau's, as reported in the literature, suggest that the charge capability of this system should be more limited compared to the discharge capability. Analysis of Cheng et al.'s Mo 6 S 8 system 21 (20-μm thick) indicates that the slow phase represents the greatest source of energy loss for the battery. This study demonstrated the importance of examining the surface kinetics of the next-generation-cathode materials.
The study also showed that, in cells with practical electrode thickness, the APC/THF electrolyte might be a major limiting factor of the performance. The additional analysis revealed that an improved electrolyte will be essential to achieve good performance for the nextgeneration-multivalent-battery systems. equation, one initial condition and two boundary conditions are necessary: where c s,0 is the initial concentration of solid, d is the diameter of the particle, i n is the transfer current, and F is the Faraday constant. Eq. A3 is the symmetric boundary condition at the center of the particle, and Eq. A4 is the flux boundary condition at the particle surface.
Mass balance in electrolyte.-This part is applied to the separator and the cathode region with the governing equation where ε e is the porosity, c e is the concentration in electrolyte, D e,eff is the effective diffusion coefficient, t 0 + is the transference number, a V is the specific area, and x is the coordinate in the thickness direction (see Figure 3). The effective diffusion coefficient can be calculated based on the porosity and Bruggeman coefficient (D e,eff = ε Bruggeman D e ). Eq. A5 also requires one initial condition and two boundary conditions: where c e,0 is the initial concentration of the electrolyte; L pos and L sep are the thicknesses of the positive electrode and separator, respectively (see Figure 3); and i app is the applied current density. Eq. A7 is the no flux boundary condition at the positive electrode surface, and Eq. A8 is the flux boundary condition at the negative electrode surface.
Charge balance in solid phase.-Ohm's law describes the electronic phase: where σ eff is the effective conductivity of solid and φ s is the solid potential. Two boundary conditions are essential: Eq. A10 is the electronic current boundary condition at the interface between the positive electrode and current collector, and Eq. A11 represents the electrical insulation between the positive electrode and separator.
Charge balance in electrolyte.-A modified Ohm's law describes the ionic phase: where κ eff is the effective conductivity of the electrolyte, φ e is the electrolyte potential, R is the gas constant, f +-is the mean molar activity coefficient, and T is the temperature. The activity coefficient is assumed to be constant in this work. [A14] Eq. A13 represents the physical situation at the interface between the positive electrode and current collector: no ions can be transferred. Eq. A14 describes the anode surface: all the kinetic reactions will be ionic current.
Butler-Volmer kinetics.-For the reaction rate, the Butler-Volmer kinetic expression is implemented: where α a and α c are the anodic and cathodic transfer coefficients, respectively, and U is the open-circuit potential. In this work, the discharge and charge curves presented by Aurbach et al. 5 were averaged to obtain the open-circuit potential. The overpotential of the cycling voltage of Reference 5 was less than the overpotential of the lowest C-rate of Reference 21. Although Reference 5 does not provide information of the cycling C-rates, the lower overpotential suggests that the cell should be closer to its equilibrium state. Therefore, we adopted the OCV from Reference 5 in this paper. The exchange-current density of the positive electrode is defined in Eq. 4.