Modeling of Magnesium Intercalation into Chevrel Phase Mo 6 S 8 : Report on Improved Cell Design

A good understanding of the limiting processes in rechargeable magnesium batteries is key to develop novel high-capacity/ high-voltage cathode materials. Thereby, the performance of magnesium-ion batteries can strongly depend on the morphology of the intercalation cathode. Moreover, high mass loadings are essential for commercialization. In this work the influence of different mass loadings are studied in addition to the impact of the particle size distribution of the active material. Therefore, a detailed continuum model is developed, which is able to describe the complex intercalation of magnesium into a Chevrel phase (CP) cathode. The model considers the thermodynamics, kinetics and interplay of the two energetically different intercalation sites of Mo 6 S 8 , which results from its unique crystal structure, as well as the impact of the desolvation on the electrochemical reactions and possible ion agglomeration. Ideal combinations of mass loading and electrolyte concentration as well as the desired CP particle size are determined for the state-of-the-art magnesium tetrakis(hexafluoroisopropyloxy)borate Mg[B(hfip) 4 ] 2 electrolyte.


Introduction
One essential challenge of the current century is to meet the globally increasing energy demand on a sustainable basis.Therefore, renewable energy sources as well as appropriate energy storage technologies gain in importance.The latter includes the development of rechargeable high energy density batteries based on abundant elements.[6] Magnesium is not only very abundant [1][2][3] but also less prone to dendrite formation than many other metals. [7,8]Together with the bivalency of the magnesium cations the resulting possibility to safely use a metal anode enables batteries with high specific capacities.
][11][12] Consequently, the choice of suitable electrolytes and cathode materials is not straightforward.Reversible magnesium insertion with reasonable kinetics was possible for the first time with a Chevrel Phase (CP) cathode (Mo 6 S 8 ). [13]isadvantages of CP are its comparatively low specific capacity and potential against Mg metal.14À 27] Key to explain the electrochemical behavior of the CP cathode is its unique crystal structure, which is built from Mo 6 S 8 face-centered cubes.Thereby, the eight sulfur atoms occupy the corners of the cube whereas the six molybdenum atoms are located on the faces, forming an octahedron.These Mo 6 S 8 cubes are connected in a three-dimensional way, whereby two different types of cavities can be occupied by magnesium cations.Each of the cavities provides 6 equivalent lattice sites for magnesium intercalation, which often are referred to as inner and outer sites and which together form 3D transport channels.However, due to electronic and steric limitations only two of the 12 potential intercalation sites per Mo 6 S 8 unit are simultaneously accessible for magnesium intercalation, which comprise one of the inner and one of the outer sites.[30][31] The kinetics and dynamics of Mg insertion and transport in CP has been well studied by Aurbach et al. [15,19,20,[32][33][34][35][36][37] Moreover, it is known, that chloride anions play an important role for efficient magnesium intercalation into CP.The chlorides as well as the CP itself can catalytically facilitate the intercalation process. [15,23]][40][41] Thereby, the presence of chloride anions is also beneficial for the desolvation process. [1,9,10,42,43]][46][47][48][49][50][51][52] Recently, successful and reversible magnesium intercalation into CP was shown with this promising chloridefree magnesium borate electrolyte. [17]owever, the processes in the CP as well as at its surface and the resulting limitations are not yet fully understood.][55] For instance, smaller particles and thinner electrodes are favorable for the solid-state diffusion and for the ionic transport in the electrolyte, respectively.
The aim of this work is to provide a more detailed understanding of the influence of the CP morphology on the battery performance, whereby the effect of the particle size distribution as well as the effect of different mass loadings are analyzed.Therefore, a comprehensive continuum model for a CP cathode is developed, which explicitly considers, that magnesium intercalates into two sites with different thermodynamic and kinetic properties.The model is then applied to cells using the state-of-the-art Mg[B(hfip) 4 ] 2 /DME electrolyte.Moreover, the impact of the desolvation on the electron-transfer reaction [10] as well as possible ion agglomeration [56] are included.The simulation results are validated and supported by experimental data.

Effect of an inhomogeneous particle size distribution
It is well-known, that the size of the Mo 6 S 8 particles plays an important role for magnesium intercalation into CP. [16,55]specially for the double charged magnesium cation the solidstate diffusion in the host is the most critical step for decent charge storage kinetics.Smaller particles thereby shorten the diffusion path, which is beneficial for the battery performance since transport limitations due to slow solid-state diffusion are less pronounced.Practically, the active material always consists of multiple particle sizes.Therefore, one should choose multiple representative particle sizes to get a better description of the material properties and consequently better predictions by the simulations.However, due to the additional complexity of considering the detailed particle size distribution, which also makes the calculations much more time-intensive, many simulations are based on a single average particle size.This simplification is often reasonable when the particle sizes are normally and sufficiently narrow distributed.
The Mo 6 S 8 material used in our study can not be wellrepresented by a Gaussian distribution (Figure S1).Therefore, we specifically investigated the differences coming from describing the Mo 6 S 8 particles by an average size (5.90 μm) and considering multiple particle sizes.In the latter case the results from a two peak fit (48.3 vol% 9.75 μm and 51.7 vol% 1.26 μm) were used on the one hand and the whole particle size distribution (Table S2) represented by 29 discrete particle sizes was considered on the other hand.The comparison of the voltage profiles can be found in Figure 1a and b.The corresponding differential capacity analysis is shown in Figure S8 of the Supporting Information.
Important to notice is that the shape of the voltage curve resulting from simulations including two and 29 different (Table S2) particle sizes is very similar (Figure 1a,b).The only difference is a slightly lower capacity in the case of the fully considered particle size distribution.But overall, the simulations using two representative particle sizes are in good agreement with the simulations including the whole particle size distribution.Therefore, we conclude that two representative particle sizes are in the present case sufficient to predict the battery performance without unreasonably raising the calculation time compared to the simulations based on the detailed particle size distribution.
More eye-catching is the pronounced difference between the simulations considering multiple particle sizes and the one based on the mean particle size (Figure 1a,b).Thereby, not only a pronounced difference in the reachable capacity can be seen, but also the general shape of the voltage profiles differ significantly: The first voltage plateau is less pronounced in the case of one representative particle size and moreover, the voltages of the two plateaus are slightly lower.Figure 1c-f shows the corresponding amount of intercalated magnesium ions in the case of a mean particle size and two particle sizes and enables to get more insight into the observed differences.It can be seen that almost all inner and outer sites of the averaged sized particle (especially in the middle of the particle) are already fully filled during the first discharge (Figure 1c,d).This is also the case for the smaller particles of the two particle size simulation (Figure 1e,f).In contrast the bigger particles show significant transport limitations so that the intercalation into the middle of the particle is not complete during the first discharge.Interestingly, the magnesiation of the middle of the large particles continues during charge since the incomplete filling causes a concentration gradient which drives the solidstate diffusion into this undesired direction (Figure 1e,f).That changes as soon as enough magnesium is extracted, so that the magnesium concentration close to the surface is lower than in the middle of the particle.This effect, which finally intensifies trapping of magnesium, is only observed for the large particles (9.75 μm) and can't be seen for the averaged sized ones (Figure 1c,d).Therefore, the average particle sizes can provide a larger capacity during the first discharge (Figure 1a).Due to trapping of magnesium ions in the inner sites the capacity of the following charge is significantly reduced.This is also true for the second cycle (Figure 1b), whereby the charge and discharge capacity are very close since an almost fixed amount of trapped magnesium is reached after the first cycle.These simulation results are consistent with previous reports about the intercalation characteristics of CP. [33,57] The trapping in the inner sites also causes the first voltage plateau at around 1.1 V  S2).The corresponding occupancy for the inner (c,e) and the outer (d,f) intercalation site are shown for the simulations with the mean CP particle size (c,d) and the two CP particle sizes (e,f).Dashed vertical lines indicate discharge and charge phases, respectively.
to be less pronounced in the second and following cycles compared to the first one (Figure 1a,b).Figure 1c and e reveal that this undesired trapping is very pronounced for the average-sized particles as well as for the bigger particles of the simulations with two different particle sizes: After a magnesium ion is intercalated into the middle of the particle and occupies an inner site deintercalation becomes improbable at relevant applied currents (< 5 %).However, in the case of the smaller (1.26 μm) particles more than 60 % of the ions can be reversibly inserted into the inner sites of the middle of the particle (Figure 1e).Consequently compared to the mean-sized particles the trapping of magnesium is lower in the case of two particle sizes.Therefore, the capacity during the second cycle is significantly higher in the simulations with two representative particle sizes (Figure 1b).Moreover, it can be seen from Figure 1c and e that the inner sites at the surface of the average-sized particles are much faster occupied compared to the two different sized particles.Therefore, the first plateau is less pronounced for the simulations with a mean particle size (Figure 1a,b).
An interesting finding from Figure 1a,b is that in the case of multiple particle sizes an additional third plateau at ca. 0.75 V appears at the end of the discharge process, which is not observed for the simulations based on the mean particle size.Thereby, the position of the 3rd plateau in the voltage profile indicates, that it is caused by intercalation into the outer sites.This assumption is confirmed by simulations, which only allow intercalation into one of the two different sites and do not consider an exchange between the intercalation sites (Figure S9a,b).To get a better understanding about the different intercalation behavior of the inner and outer sites, the corresponding amount of intercalated magnesium ions is also analyzed and can be found in Figure S9c,d.
The beneficial effect of a reduced diffusion path can clearly be seen for both intercalation sites.In the smaller particles significantly more magnesium ions can reversibly intercalate into the middle of the particle.However, to explain the origin of the 3rd plateau the occupancy of the intercalation sites at the surface of the particle has to be examined, since this is the parameter that determines the OCV (Figure S9e,f) and therefore the cell voltage.In the case of the inner sites there is almost no difference in the amount of intercalated magnesium ions at the surface of the smaller and the bigger particle, so that a similar OCV and consequently a single voltage plateau is observed (Figure 1e and Figure S9a,c,e).In contrast, the ion insertion and extraction into and from the outer sites of the bigger particles is only initially equally fast as in the smaller particles but then gets significantly delayed (Figure 1f and Figure S9d).This means that the outer sites of the smaller particle, not only in the middle but also at the interface, are filled faster than in the case of the bigger particle.Consequently, the additional plateau represents the situation where the outer sites of the small particles are almost completely filled but intercalation can still take place into the bigger particles, so that the OCV of the two different sized particles differs significantly (Figure 9b,d,e).Since the solid-state diffusion is much slower on the inner sites (Table 1) an analogue situation does not occur in the case of these sites and the corresponding first voltage plateau.For the inner sites the transport limitations prevent a complete intercalation into the middle of the particle even for the 1.26 μm ones (Figure S9c).
Of course the possible exchange in between the sites complicates the intercalation behavior.However, the similar occupancy of the inner sites at the surface of the two different sized particles and a delayed filling of the outer sites at the surface of the bigger particles is also observed for simulations considering both intercalations sites as well as the exchange (Figure 1e,f), which finally leads to the occurrence of a third plateau in the overall voltage profile (Figure 1a,b).

Interpretation of the simulation results in comparison to experimental data
In addition to the comparison of simulations based on different approaches to include the particle size distribution (Figure 1), also the comparison to the corresponding experimental results (Figure 2) are revealing.Thereby, the mass loading of the CP cathode and the C-rate of the simulations were chosen accordingly to the experiment.In addition to the dischargecharge curves (Figure 1a and 2) a differential capacity analysis is shown in Figure S10.In contrast to the simulations the measurements show a significantly lower capacity for the first discharge than for the following ones (Figure 2).This implies that during the first discharge an unfavorable conditioning effect dominates over the beneficial situation that the intercalation host is fully empty and there is no trapping yet.In general, it is well-known that magnesium electrolytes including Mg[B(hfip) 4 ] 2 often require conditioning to activate the elec- trode surfaces. [17,48,58]Moreover, recently it was reported, that the B(hfip) 4 -anion can be oxidized at high voltages causing the formation of a cathode-electrolyte interphase, which facilitates the desolvation of Mg(DME) 3  2 + and therefore the intercalation into CP. [59]These kind of side reactions can cause low capacities as well as an asymmetry between charge and discharge during the first cycle(s).Since the model does not consider such complex activation processes, the comparison of the simulation results with the experimental data focuses on later cycles.However, the charge and discharge capacity already become very similar from the second cycle onwards (Figure 2), which is also observed in the simulations and indicates a fixed amount of trapping.Moreover, it can be seen that the capacity reached in the experiment is slightly lower than predicted by the simulations (Figure 1a,b and 2).Additionally, the first plateau is rarely observed in the experiment.These two findings imply that the solid-state diffusion in between the inner sites is significantly slower than assumed in the simulations.Therefore, a parameter study regarding the corresponding diffusion coefficient D 1 was performed (Figure S6a).The analysis shows that capacity and shape of the first plateau of experiment and simulation matches better when the diffusion between the inner sites is about one order of magnitude slower than initially assumed (Figure 2 and S6a).Since the parameterization of the diffusion coefficient was based on PITT measurements, which can only measure the effective transport parameter, more precisely the combination of ion hopping between same sites and an exchange between the different sites, it seems reasonable that the diffusion coefficient D 1 of the model, which only accounts for hopping between inner sites, should be lower than the experimental data.In contrast, D 2 of the outer sites should be very similar to the experimental value because the measurement represents the case where the inner sites are already occupied and impact of the exchange should be very small so that the effective transport should be dominated by the hopping between outer sites.This is also supported by a corresponding parameter study (Figure S6b) and the parameter optimization during validation of the model (Section S2.1.of the Supporting Information).However, the most interesting feature of the measured voltage curves is that there is an additional small plateau at about 0.75 V. Thanks to the simulations with different approaches to consider the particle size distribution (Figure 1) the extra plateau can be ascribed to the non-Gaussian distribution of Mo 6 S 8 particle sizes (Figure S1).All in all, the different simulations (Figure 1) and the comparison to experimental data (Figure 2) clearly show, that considering two different Mo 6 S 8 particle sizes is a simple and efficient hence effective way to include the morphological impact of multiple particle sizes on the performance of Mo 6 S 8 cathodes into the simulations.Therefore, all following simulations are based on two different particle sizes of the active material Mo 6 S 8 representing the corresponding share of active material in the electrode formulation.Note, that as seen in Figure 1a,b this approach slightly overestimates the reachable capacity.

Effect of particle size on the electrode performance
Since smaller Mo 6 S 8 particles show significantly less transport limitations by the slow solid-state diffusion of magnesium (Figure 1e,f and S12) they enable an improved battery performance, which is shown in Figure 3.For instance with big 9.75 μm particles only 49 % of the theoretical capacity can be reached from the second cycle onwards (Figure S11b and 3) since inaccessible sites and trapping in the middle of the particle are a big issue (Figure S1).In contrast, the usage of smaller 1.26 μm particles is extremely beneficial for the performance limiting solid-state diffusion (Figure S12) and enables to reach almost twice the capacity of the bigger 9.75 μm particles (Figure S11  and 3).The herein studied mixture of different particle sizes contains a significant amount of small particles (Figure S12), which are crucial for a sufficient battery performance.Nevertheless, Figure 1 and S11 clearly show that the sluggish solidstate diffusion still limits the reachable capacity.The analysis of the influence of different particle sizes on the battery performance shown in Figure 3 reveals that limitations by solid-state diffusion are minor for Mo 6 S 8 particles smaller than about 0.1-0.5 μm.Especially at higher, more practical C-rates a small  diffusion length is crucial to reach sufficient capacities (Figure 3).However, the larger surface area of smaller particles may not only be beneficial for the intercalation reaction but could also enhance undesired side reactions.Moreover, small particles might lead to unfavorable electrode morphologies causing transport limitations for magnesium ions in the electrolyte.Finally, processing of the material might be more challenging on an industrial scale.Nevertheless, small particle sizes can be considered to be key to a good performance of CP and other Mg intercalation cathodes (Figure 3).

Influence of the mass loading on magnesium intercalation into Chevrel Phase
High mass loadings are still an essential bottleneck for commercialization of magnesium-ion batteries.Therefore, we also investigate the impact of different mass loadings on the battery performance.

Positive and negative effects of the loading on galvanostatic cycling
The voltage as well as the electrolyte concentration profiles of the second galvanostatic cycle at 0.1 mA cm À 2 for cathodes with different mass loading are shown in Figure 4.As seen in Figure 4a, a higher capacity is accessible when the CP cathode has a higher mass loading.Thereby, the first plateau is significantly longer and, additionally, the voltage of the plateaus is higher for higher mass loadings.The main reason for this is that the solid-state transport limitations, which usually limit the reachable capacity, are mitigated by less demanding operation conditions (Figure S14).This can be seen in the inverse correlation between the C-rate and the mass loading for a certain current density [Eq.( 16)], which shows that there is more active material involved in providing the desired current density, so that the local current density at the surface of each Mo 6 S 8 particle is smaller, which diminishes limitations by solid-state diffusion.Therefore, a higher mass loading enables a larger amount of magnesium ions to reversibly intercalate (Figure S14), which in turn enables to reach a higher capacity.This effect can be seen most prominently for the inner sites of the smaller particles (Figure S14c).In the case of a low mass loading (< 5 mg cm À 2 ) more than 50 % of the magnesium ions are trapped in the middle of the 1.26 μm particles, whereas for higher mass loadings (> 15 mg cm À 2 ) about 90 % of the ions can be extracted from the inner sites in the middle of the small particles (Figure S14c).This enhanced reversibility of magnesium intercalation into the inner sites is also the reason why the first plateau is significantly more pronounced for higher mass loadings (Figure 4a).Moreover, for the lowest analyzed mass loading of 1 mg cm À 2 the transport limitation by sluggish solid-state diffusion even shows for the outer sites of the small particle.As revealed in Figure S14d they can't be filled completely anymore, which can also clearly be seen in a missing third plateau at 0.75 V (Figure 4a).However, higher mass loadings can also become problematic for the battery performance since they also cause more pronounced concentration gradients in the electrolyte (Figure 4b).Therefore, a higher mass loading is only beneficial as long as the current density is small enough so that no detrimental transport limitations occur in the electrolyte.As seen in Figure 4b this is the case for 0.1 mA cm À 2 , where ion depletion and accumulation is still rather small.Even for a very high mass loading (25 mg cm À 2 ) the maximum electrolyte concentration is significantly below the critical concentration for ion agglomeration of about 0.35 M [56] .Consequently, the positive effect of an enlarged mass loading dominates.
In conclusion, an enormous impact of the mass loading on the solid-state diffusion and therefore on the battery performance can be observed, when a constant current is applied.In contrast, Figure S15a,c shows that the different mass loadings can provide a similar capacity at a given C-rate, which is expected since the local current density at the surface of the individual Mo 6 S 8 particles and therefore the solid-state diffusion limit is identical.However, the concentration gradients in the electrolyte become significantly more pronounced for the higher mass loadings even at the rather small C-rate of C/10 (Figure S15b).This is the reason, why the overpotential for the electrochemical reaction becomes slightly higher for high mass loadings, which results in a lower voltage of the two plateaus as seen in Figure S15a.

Interpretation of experimental data
The typical changes of the voltage profiles caused by different mass loadings as implied by Figure 4a can also be observed in the experimental data shown in Figure 5.As predicted by the simulations (Figure 4a and S13), the first plateau is significantly more pronounced for a higher mass loading and the additional plateau appears at about 0.75 V (Figure 5).These two observations indicate that the limiting effect of the solid-state diffusion is less pronounced for the CP cathode with the higher mass loading.However, this positive effect on magnesium intercalation does not lead to an improved capacity (Figure 5).Since the transport in the electrolyte should be sufficiently fast at 0.1 mA cm À 2 (Figure 4b), there seems to be another critical effect, which limits the reachable capacity.This could be side reactions at the cathode-electrolyte interface, but also electrical contact issues might arise in thicker electrodes.Since magnesium cations can't intercalate into uncontacted CP particles, this ‚dead' fraction of the active material means a loss of specific capacity and would result in effectively higher local currents.Nevertheless, the measured capacity for the low mass loading (2.75 mg cm À 2 ) is very close to the one predicted by the simulations (Figure 5a and S13).Consequently, the additional limitation only occurs for high mass loadings.Contact issues are more probable for thicker cathodes, since the distance to the current collector becomes longer.Moreover, an uneven distribution of conductive carbon black can cause that more CP particles loose contact and can't be accessed for magnesium intercalation.But from Figure 5 it is also observed that the conditioning takes significantly longer for thicker CP cathodes.Therefore, an extended study of the interfacial processes between Mo 6 S 8 and the Mg[B(hfip) 4 ] 2 /DME electrolyte is needed.Moreover, the role of conductive additive and the cathode manufacturing process need to be investigated.However, this is beyond the scope of this work.

Ideal combination of mass loading and electrolyte concentration
Altogether, it becomes clear that the cathode mass loading significantly impacts the magnesium transport in the CP particles as well as in the electrolyte.On the one hand, a high mass loading creates less demanding local conditions, which mitigates solid-state diffusion limitations and enables a battery operation at higher current densities.But at the same time an increase of the mass loading also leads to more pronounced concentration gradients in the electrolyte.Since magnesium electrolytes suffer from a poor salt solubility, higher mass loadings can easily cause ion accumulation and/or depletion over/under critical concentrations, especially at higher C-rates (Figure S16).To avoid a failure due to precipitation or shortage of magnesium ions, an upper limit for the CP mass loading depending on the electrolyte concentration and the C-rate or current density, under which the battery will be operated, is determined and the results can be found in Figure 6c, S16 and S17.
Consequently, the choice of mass loading is always a compromise between delivering high capacities and current densities by mitigating solid state diffusion limitations and avoiding transport limitations in the electrolyte.For a certain Crate the mass loading can be maximized by choosing the appropriate electrolyte salt concentration.The ideal combination of mass loading and electrolyte concentration for different C-rates can be found in Figure 6a as well as in Figure S17.For instance, the CP mass loading of a battery operating at C/10 should not exceed 51 mg cm -1 and the corresponding ideal salt concentration, which enables to avoid both precipitation and magnesium depletion, is found to be 0.25 M (Figure 6a and  S17a).Furthermore, it is proposed that also ion agglomeration happening at even lower concentrations (c max = 0.35 M) has a negative effect on battery performance, [56] .If this regime has to be avoided the mass loading and electrolyte concentration should both be significantly lower.In this case, the ideal battery for operating at C/10 is filled with a 0.18 M electrolyte and the maximum CP mass loading is 42 mg cm -1 (Figure 6a and S17a).Important to notice is, that the ideal electrolyte concentration slightly increases for operation at higher C-rates while the maximum mass loading becomes significantly lower.For cycling at 1 C critical transport limitations in the electrolyte already occur for CP mass loadings of about 14 mg cm -1 (Figure 6a and S17c).Figure 6b summarizes the ideal combinations of mass loading and electrolyte concentration in dependence of the C-rate and indicates the reachable current density.It is observed that the change of the ideal mass loading and concentration is significantly more pronounced for C-rates between 0.1 and 0.5 C, whereas they are very similar in the range of 0.5-1 C.Moreover, it becomes clear that CP cathodes with a smaller mass loading are able to provide significantly higher current densities without experiencing critical transport limitations in the electrolyte (Figure 6b).
All in all, the basis for magnesium batteries with high cathode mass loadings is a good ionic transport in the electrolyte.Therefore, the choice of electrolyte is crucial.However, poor salt solubility and/or ionic conductivity are issues for most magnesium electrolytes that are compatible with both electrodes.

Influence of ion agglomeration on magnesium intercalation into Chevrel Phase
The previous simulations clearly showed, that quite high concentration gradients can occur in Mg-CP cells and that an increase of the CP mass loading and/or the C-rate can easily cause the ion concentration to locally exceed the critical concentration for ion agglomeration (Figure 4, 6, S15 and S17).
To get an impression of the impact of the clustering, the model was expanded by a description for ion agglomeration in Mg[B(hfip) 4 ] 2 /DME. [56]The corresponding results are shown in Figure 7 and S18.
Interestingly, even at small current densities, where only neglectable concentrations of ion agglomerates are formed in the electrolyte (c Cluster < 10 À 4 moL L -1 ), minor differences can already be observed between simulations with and without considering clustering (Figure 7a and S18a).The plateaus are shifted to slightly lower voltages (Figure 7a) and the concentration gradient in the cell, especially at the cathode side, gets slightly more pronounced (Figure S18a).Nevertheless, these difference are very small at C/10 and can still be considered as neglectable.In contrast the effect of clustering is getting very prominent at higher current densities (Figure 7b and S18b).At a current density of 1 mA cm À 2 (= 0.67 C) for instance the cluster concentration already reaches up to almost 0.2 M locally.The resulting shortage of free ions clearly impacts the battery performance.On the one hand the bulky ion agglomerates hamper the ionic transport in the electrolyte, which can be seen as self-amplifying effect since this causes even more pronounced concentration gradients in the cell (Figure S18b).Consequently, the formation of ion pairs and prenucleation clusters causes peaks in local concentration surpassing the solubility limit.Therefore, this effect possibly promotes precipitation, which in the end might be fatal for the battery performance.On the other hand clustering hinders the electrochemical reaction at the electrode-electrolyte interface, which leads to an increase in overpotentials (Figure 7b).Moreover, the presence of ion agglomerates causes an increased slope of the two voltage plateaus and the first plateau is significantly less pronounced (Figure 7b).All in all, the adverse effects of ion aggregation can clearly be seen in the simulations (Figure 7).Nevertheless, even with a significant amount of clusters in the electrolyte the battery performance can still be satisfying.The shape of the voltage profile can act as indicator for the occurrence of ion agglomeration, which further implies that concentrations close to the solubility limit are present during cycling.Such an observation could possibly serve as an indicator for precipitation and be exploited for battery management.Altogether it is best for the battery performance, when ion agglomeration can be avoided.Therefore, the mass loading and electrolyte concentration should be chosen according the the desired Crate or current density of battery operation (Figure 6b).

Conclusion
In summary, our proposed intercalation model considers the thermodynamics and kinetics of multiple intercalation sites.The model provides important insights into transport limitations in Chevrel phase particles as well as in the Mg[B(hfip) 4 ] 2 / DME electrolyte.A detailed understanding of these limitations is key to improve the performance of magnesium-ion batteries.We found that the efficiency of the solid-state transport, which is usually limiting the reachable capacity, is determined by the particle size of the active material on the one hand, whereby Chevrel phase particles smaller than 0.5 μm are desireable to avoid transport limitations over a wide range of practical current densities.One the other hand a high mass loading is beneficial for solid-state transport of magnesium ions in the Chevrel phase cathode.However, high mass loadings can also cause critical concentration gradients in the electrolyte and both salt precipitation and total magnesium depletion can become an issue.Therefore, mass loading and electrolyte concentration have to be balanced and chosen depending on the operation condition.Hence, our detailed study on the ideal combination of those three parameters can act as a useful guide for battery design.
Moreover, it becomes apparent, that the cathode morphology can not only determine the amount of trapped magnesium in the crystal structure and therefore the practical capacity but also the general shape of the discharge curves.Latter, is also affected by the occurrence of ion agglomerates, which can act as additional indicator for critical transport limitations in the electrolyte.
All in all, the morphology of the intercalation cathode significantly impacts the performance of magnesium ion batteries.This work is a first step towards a better understanding of the limiting factors and simulation-based optimization of magnesium-ion batteries.To better capture all the effects of the cathode morphology, further research on improving the battery design should include an extension of the herein presented pseudo-2D model to a 3D-resolved description of the Chevrel phase cathode.Moreover, future work should address the role of the solvent/electrolyte and the corresponding kinetic description of magnesium intercalation into Chevrel phase into more detail including simulations of cyclic voltammograms.

Electrolyte synthesis
The magnesium tetrakis(hexafluoroisopropyloxy)borate (MgBOR) salt was synthesized by a reaction between magnesium borohydride Mg(BH 4 )2 and hexafluoroisopropanol in dimethoxyethane (DME) following previous work. [60]Liquid electrolyte solution was prepared by dissolving proper amount of MgBOR salt in DME (0.3 M), the concentration is based on the molecular weight of Mg[B(hfip) 4 ] 2 •3 DME.

Characterization of Chevrel Phase
To determine the sizes of the Mo 6 S 8 particles, a dispersion of the Mo 6 S 8 in isopropyl alcohol (IPA) was studied by dynamic light scattering (DLS) with a Microtrac -Bluewave particle analyzer.The porosity of calendared CP cathodes with two different mass loadings (12 mg cm À 2 and 20 mg cm À 2 ) was determined by mercury porosimetry with a Micromeritics MicroActive AutoPore V 9600, whereby a pure current collector (carbon coated nickel) without any CP coating was analyzed as control sample.

Preparation of Mo 6 S 8 cathode and mg anode
Chevrel Phase (CP) was provided by American NEI Corporation, USA.CP cathodes were prepared by a standard wet coating process.A Mo 6 S 8 slurry was prepared using 80-90 wt % CP, 5-10 wt % carbon black (Super C45, Imerys) and 5-10 wt % poly(vinyl difluoride) (PVDF) binder (Solef ® 5130, Solvay).The carbon black was dispersed into a PVDF solution (8 wt % dissolved in N-methylpyrrolidone solvent) using a highspeed disperser (Dispermat, VMA-Getzmann GmbH).Then the CP powder was added to the mixture and further dispersed until an homogeneous slurry was obtained, which was casted onto carbon-coated aluminium and nickel current collector foil (Custom cell and Gelon) with mass loadings of 2, 8, 12 and 20 mg cm À 2 .The electrodes were dried at 80 °C for 12 hours under dynamic vacuum and then calendered to 30 % porosity.
High purity mg metal foil (99.9 %, 100 μm) was provided by Gelon LIB group.Prior to cell assembly, the mg foil was scratched carefully on both sides in the glovebox to remove the native surface layer.

Electrochemical measurements
2025 and 2032 type coin cells consisted of one layer of pure magnesium metal foil anode, single-side coated cathode and a vacuum dried glass fiber membrane (Whatman, GF/C, 260 μm) in between wetted with 80 μL of electrolyte solution.The assembly of the coin cells was done in a high-purity argon (99.999 %) filled glovebox under ideal conditions (O 2 < 1 ppm, H 2 O < 1 ppm).Galvanostatic charge-discharge experiments were performed in a voltage range of 0.4-1.6 and 0.5-2.5 V vs. Mg/Mg 2þ and at a current rate of 0.1 C with an Arbin battery cycling unit under 25 °C.

Continuum model
The presented magnesium-ion battery model is a typical 1 + 1D (pseudo-2D) model, which describes the transport in the liquid electrolyte [Eq.(S2)] and in the solid particles of the CP [Eq.(7)] as well as the electrochemical reactions at the electrode-electrolyte interface [Eq. ( 2)].The required potentials of the electrodes and the electrolyte are obtained by solving the corresponding charge balance equations [Eqs.(S1) and (S3)] [61,62] Moreover, the unique crystal structure of the Chevrel phase cathode, which provides two possible sites with different thermodynamics and kinetics for magnesium intercalation, is explicitly considered.A schematic illustration of the pseudo-2D model for a magnesium-ion battery with a CP cathode is shown in Figure 8.However, the herein described model is of more general nature and can easily be used to describe other intercalation materials, where multiple types of cavities can be occupied. [63,64]Also the transfer to blended cathode materials consisting of two or more active intercalation compounds is straightforward.

Differentiation of individual contributions to the open circuit potential (OCV)
The open circuit potential (OCV) represents the thermodynamics of an intercalation material.Thereby, the Gibbs free energy DG 0 i of the general intercalation reaction Mg 2þ sol þ #i þ 2e À Ð Mg i is connected to the chemical potentials μ 0 of all reactants as well as to the corresponding OCV U 0 i [Eq.( 1)].
whereby z þ describes the number of elementary charges per intercalating ion, which is equal to the number of transferred electrons and F is the Faraday constant.
The employed continuum model differentiates between two different vacancies in the material for the intercalation magnesium.Therefore, the individual contributions of the sites U 0 i to the overall OCV need to be extracted from the measured data.The here chosen approach is based on a substitutional lattice model and can be applied in a very general context to energy storage materials. [63,65]Following previous work on lithium-ion battery cathode materials, [63] the process is transferred to cathode materials with multiple-electron intercalation reactions.
Instead of considering only a single reaction for ion-intercalation, we distinguish the individual intercalation reactions into the different vacancies of the host material.Note, that this choice needs to be made based on information about the host material in order to link a physical meaning to the assigned sites.Since the Nernst equation describes the ideal open-circuit potential of an electrochemical reaction, Equation ( 2) is chosen as an expression for the open circuit voltage of a site i.
Therein, U 0 i;ref denotes the standard potential of the vacancy, x i ¼ c i =c max;i the fraction of sites occupied by the intercalated atoms and X i the total available fraction with 0 < x i < X i .R and T are the ideal gas constant and the temperature, respectively, and ω i incorporates deviations from ideality in the equation and can be linked to an activity coefficient. [63]Assuming equilibrium conditions for all sites U 0 :¼ U 0 1 ¼ U 0 2 ¼ ::: and inverting Equation ( 2), one obtains an expression for the total amount of occupied sites x. x This expression can then be used to split an overall OCV into the different OCV-parts related to the individual sites by employing a fit-algorithm to given OCV-data like the method of least squares.Note that this method has to be applied to the expression x U 0 ð Þ, as an explicit expression for the inverse of Equation ( 3) is not straightforward.
The above described procedure can also be applied by assuming any other function, that describes the chemical potential of a chemical reaction in a reasonable way and incorporates parameters ω i resembling deviations of an ideal shape.It should also be noted, that the main goal of this ansatz is to capture the characteristic features of an OCV and associate them with physical processes like specific vacancies in the host material's lattice structure.Therefore, discrepancies of the resulting differentiated OCV from measured data are expected as the procedure uses only a minimal amount of basis functions for the fitting-method.

Kinetic model for intercalation
[68][69] When the electrolyte contains multiple electroactive species or intercalation can happen into multiple sites or materials, the total current density at the electrode-electrolyte interface i se is given by the sum of the individual interface reactions [Eq.( 4)].
Thereby, the Butler-Volmer equation is used to calculate the interface current caused by intercalation into each of the different lattice sites of the CP cathode [Eq.( 5)].
|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } whereby k BV;i denotes the corresponding rate constant for intercalation into site i and a BV;i is the so-called symmetry factor, which enables to consider that the cathodic and anodic reaction can have different contributions to the overpotential h s;i .The latter is defined as deviation of the experimentally observed potential from the corresponding standard half-cell reduction potential given by the OCV [Eq.(6)].Details about the determination of the individual OCVs of the two different intercalation sites (U 0 i ) from the overall OCV can be found in the previous section (Differentiation of individual contributions to the open circuit potential (OCV)).
In the case of multivalent cations like magnesium it is known, that the desolvation plays an important role during the deposition or intercalation process, [1,9,[38][39][40][41] which would require advanced models to describe the electron-transfer kinetics.However, a careful parameterization of the Butler-Volmer equation, especially regarding the symmetry factor α BV,i , also enables to consider desolvation effects.Since findings from earlier work [10] enable a well-founded parameterization of the Butler-Volmer equation, it provides a sufficient estimate of the interface kinetics for modeling of magnesium intercalation into CP.The integration and development of more advanced models for the intercalation kinetics will be subject of future work.

Solid-state diffusion and exchange between the sites
The effective transport of magnesium in CP results from the combination of two different mechanisms (Figure 8 and Eq. ( 7)): On the one hand there is hopping between energetically similar sites (1 $ 1 and 2 $ 2), which can be described as solid-state diffusion Ñs;i ¼ D i rc i with the diffusion coefficient D i for the corresponding site i.On the other hand the magnesium can hop from one type of site to the other type of intercalation site (1 For describing this exchange @c ex i @t between the different sites the following equilibrium has to be considered whereby Mg 2þ i stands for an occupied and #i for an unoccupied site i.The corresponding rate equation for the exchange is based on the rate constants k for the forward (k 1!2 ) and backward (k 2!1 ) reaction as well as on the concentration of the reactants c [Eq. ( 8)].
Thereby, the concentration of unoccupied sites is given by c #i ¼ c max;i À c i , whereby the maximum concentration of intercalated magnesium c max;i is determined by the stoichiometry, the molar mass M and the density 1 of the intercalation material [Eq.(9)].Due to its crystal structure the maximum magnesium concentration is equal for both types of intercalation sites in CP.
Moreover, the two rate constants k 1!2 and k 2!1 are coupled to each other by the equilibrium constant K of the exchange reaction, which is determined by the different energetics of the inner (1) and outer (2) sites [Eq.(10)].
Combining Equations ( 8) -( 10) results in the following expression for the exchange rate [Eq.( 11)]: The standard Gibbs free energy of the exchange reaction DG 0 exchange , which is required to calculate the equilibrium constant K [Eq.(10)], is given by the sum of the stoichometric weighted chemical potentials μ 0 of all reactants [Eq.13)].
The combination of Equations ( 1), ( 10) and ( 12) enables to express the equilibrium constant K of the exchange in terms of the OCVs of the two different intercalation sites [Eq.( 13)].
Since the inner sites are known to be energetically more favorable, the equilibrium constant will be smaller than 1.However, the OCV of an intercalation material is dependent on its state of charge (SOC) and therefore on the concentration of intercalated magnesium.The determination of the OCVs of the two individual sites in the CP (U 0 1 and U 0 2 ) can be done by a fit to experimental data, which is described in detail in the previous section (Differentiation of individual contributions to the open circuit potential (OCV)).

Kinetics of Mg plating and stripping
Analogue to the intercalation reaction in Equation ( 5), the electrontransfer kinetics for magnesium deposition and dissolution at the metal anode are described by a thoroughly parametrized Butler-Volmer equation, [10] whereby the amplitude i 0 in Equation ( 5) is adapted for the case of a conversion electrode

Transport model
The governing equations for the ionic transport in the electrolyte and electrode and electrolyte potentials based on mass and charge conservation and are summarized in the Supporting Information [Eqs.(S1)-(S3) in the Supporting Information].This general transport theory for porous systems [68] is coupled to the above described model for the transport of magnesium in CP particles [Eq.(7)] and the Butler-Volmer equation for the interface reactions [Eq.( 5)].
The resulting equation system is simplified for an isothermal process (T ¼ 298:15 K) and solved for the electrolyte concentration c e , the concentration of magnesium in the two different intercalation sites of the CP cathode c s,1 , c s,2 as well as for the electrochemical potential of the electrolyte f e and the electric potentials of the two electrodes F s .

Ion Agglomeration
Since the CP cathode is an additional porous material it can be expected that the concentration gradients in the cell will be significantly higher compared to symmetric magnesium cells especially at higher mass loadings.[72][73][74] Therefore, we finally combine the herein presented model for CP cathodes with a model from earlier work [56] , which is able to consistently consider ion agglomeration in the electrolyte.Thereby, the effects of clustering on both the transport in the electrolyte as well as on the interface reaction are described.Note, that the influence of ion agglomeration on Mg intercalation into CP is studied as last step after understanding the impact of the cathode architecture on the battery performance and unless otherwise stated all previous simulations do not explicitly consider clustering.

Model parameters
In general the model parameters are chosen based on the literature on CP as well as our own experimental data.A summary of all parameters related to the morphology of the glass fiber separator, the transport properties of the electrolytes and the electron-transfer kinetics at the magnesium metal anode is given in Table S3 of the Supporting Information and all parameters related to the CP cathode are summarized in Table 1.
It was found, that the porosity of CP coatings with different mass loadings is very similar (51.6-51.7 %, Table S1), which indicates, that the calendaring process during cathode production does not result in a different electrode microstructure.Consequently, it can be seen as a general property of our cathode material that the mass loading does not significantly impact the morphology, porosity and effective transport parameters.Moreover, the sizes of the Mo 6 S 8 particles do not follow a Gaussian distribution (Figure S2).Therefore, the mean particle size might not be an ideal input parameter for the model, especially since it is well-know that the particle size plays an important role for the battery performance. [16,55]However, considering the detailed particle size distribution (Table S2) leads to much more time-intensive calculations.As alternative, the results from a two-peak fit (Table 1) are tested to effectively describe the impact of the bimodal active material particle size distribution.
The intercalation mechanism into CP from chloride-containing electrolytes was already studied from first principles and it was found that the stripping of the chloride anion from the coordination shell of the magnesium cation is facilitated by the Mo 6 S 8surface. [23]This finding suggests, that the surface of the CP cathode can catalytically assist the intercalation reaction.Moreover, it seems that the electronic screening effect also supports the desolvation process of the magnesium cation.However, it is very likely that the magnesium cation needs at least one free coordination site to approach the CP surface and benefit from its catalytic activity.Therefore, we assume that similar to magnesium deposition the initial desolvation of the magnesium cation will be the ratedetermining step for intercalation into CP. [10]Consequently, the intercalation rate constant k BV;i as well as the symmetry factor a BV;i can assumed to be similar for both intercalation sites of the Mo 6 S 8 .
Values for k BV and α BV including the for the electron transfer required partial desolvation of Mg 2 + cations have been determined in our previous work. [10]These values are used as initial best guess for the intercalation kinetics of magnesium into CP (Table 1).To additionally get more insights about the impact of the intercalation kinetics a parameter study regarding the rate constants k BV;i and symmetry factors a BV;i was done.The corresponding results can be found in Section S1.2.3.and Figure S4 of the Supporting Information.
Moreover, different values for the exchange rate constant k 2!1 were analyzed and it was found that the differences especially from the second cycle onwards are very small (Section S1.2.2. and Figure S2).Based on this parameter study we assume that 10 -7 m 3 s -1 mol -1 is a good estimate for the exchange rate constant and use this value for all following simulations (Table 1).Regarding the solid-state diffusion in CP extensive studies can be found in literature. [15,16,19,36]It was found that the diffusion during intercalation into the inner sites is about two orders of magnitude slower than for the second intercalation step.For our simulations we use diffusion coefficients determined by potentiostatic intermittent titrations (PITT, Table 1). [19]nally, the OCV fitting [Eq.( 3)] was done with data from magnesium insertion into CP at elevated temperature (T = 60 °C) [31] , since kinetic limitations due to the electron-transfer reaction and the solid-state diffusion can be neglected.Therefore, the measured chronopotential curve will solely be determined by the different thermodynamics of the intercalation sites.The results are shown in Figure S5.
Note, that simulations in the main part of this paper are done with parameters from the literature (Table 1) using the assumptions pointed out in the paragraph above.Moreover, we compared the simulations to our experimental data and additionally determined an optimized set of parameters for the magnesium intercalation into CP by a least-square fit to the measured data (Table S4).The corresponding voltage profiles during galvanostatic cycling (Figure S7) and a detailed discussion on the optimized model parameters can be found in Section S2.1.of the Supporting Information.
Commercialization of magnesium batteries with practical energy densities requires much higher mass loadings than usually tested in lab cells.A higher mass loading m act results in a thicker cathode L ca [Eq.( 14)], which means longer pathways for the ionic transport.Consequently, transport limitations play a more important role for the battery performance.Therefore, we want to analyze the impact of different mass loadings on the magnesium-ion battery performance.The mass loading m act thereby defines the thickness of the cathode L ca [Eq.( 14)].
L ca ¼ m act e act � 1 act (14)   whereby e act describes the volume fraction of the active material Mo 6 S 8 and 1 act is the active material density.The volume fraction of CP e act is determined by the electrode composition and its porosity ε e (Eq.( 15) and Table 1).
w j denotes the mass fraction of the electrode component j and 1 j is the corresponding material density (Table 1).Moreover, the mass loading also determines the C-rate in galvanostatic cycling experiments [Eq.( 16)].
C-rate ¼ i � A C theo � m act (16)   whereby i stands for the applied current density, A for the electrode area and is the theoretical specific capacity of the CP.
The additional parameterization required for simulations considering ion aggregation is based on earlier work, whereby a cluster size of z cluster ¼ 1, describing the agglomerate Mg[B(hfip) 4 ] 2 , and the corresponding equilibrium constant for the cluster formation K cluster ¼ 0:0276 are used as basis for the calculations. [56]

Figure 1 .
Figure 1.Comparison of the voltage profiles of the first (a) and second (b) galvanostatic cycle at C/10 resulting from simulations based on 12.24 mg cm À 2 CP and a mean particle size (5.90 μm), two different particle sizes (48.3 vol % 9.75 μm and 51.7 vol % 1.26 μm) as well as the detailed particle size distribution (TableS2).The corresponding occupancy for the inner (c,e) and the outer (d,f) intercalation site are shown for the simulations with the mean CP particle size (c,d) and the two CP particle sizes (e,f).Dashed vertical lines indicate discharge and charge phases, respectively.

Figure 2 .
Figure 2. Galvanostatic cycling of a Mg-ion battery with a 12.24 mg cm À 2 CP cathode at C/10.

Figure 3 .
Figure 3. Logarithmic representation of the practical capacity after the second galvanostatic cycle for different Mo 6 S 8 particle sizes at different Crates.

Figure 4 .
Figure 4. Comparison of the voltage (a) and electrolyte concentration (b, left side = cathode) profiles of the second galvanostatic cycle at 0.1 mA cm À 2 resulting from simulations with different CP mass loadings.

Figure 5 .
Figure 5. Galvanostatic cycling of a Mg-ion battery with a 2.75 mg cm À 2 (a) and a 8.25 mg cm À 2 (b) CP cathode at 0.1 C.

Figure 6 .
Figure 6.a) Local ion depletion and accumulation for different CP mass loadings and electrolyte concentrations: Critical local concentrations for ion aggregation (0.35 M), precipitation (0.5 M) and magnesium depletion (0.01 M) during galvanostatic cycling at different C-rates.Combinations of mass loading and electrolyte concentration, which don't cause critical local electrolyte concentrations are indicated in green (lighter just for 0.1 C, darker for both C-rates), while the ideal combinations are marked with green and yellow circles.b) Ideal combination of mass loading and electrolyte concentration for providing different C-rates and avoiding critical ion accumulation and depletion in the electrolyte.The corresponding current densities are indicated.

Figure 7 .
Figure 7. Simulations with and without considering ion agglomeration: Comparison of the voltage profiles of the second galvanostatic cycle with a 12.24 mg cm À 2 CP cathode at 0.15 mA cm À 2 (0.1 C, a) and 1.0 mA cm À 2 (0.67 C, b).The corresponding concentration profiles can be found in Figure S18.

Figure 8 .
Figure 8. Schematic illustration of the 1 + continuum model for a Mg / Mg 2 Mo 6 S 8 battery.