First-principles studies on the atomistic properties of metallic magnesium as anode material in magnesium-based batteries

Rechargeable magnesium-ion batteries (MIBs) are a promising alternative to commercial lithium-ion batteries (LIBs). They are safer to handle, environmentally more friendly, and provide a five-time higher volumetric capacity (3832 mAh·cm-3) than commercialized LIBs. However, the formation of a passivation layer on metallic Mg electrodes is still a major challenge towards their commercialization. Using density functional theory, the atomistic properties of metallic magnesium, such as bulk, surface, and adsorption properties, were examined. Well-selected self-diffusion processes on perfect and imperfect Mg surfaces were investigated to better understand the initial surface growth phenomena. Subsequently, rate constants and activation temperatures of crucial diffusion processes on Mg(0001) and Mg(101!1) were determined, providing preliminary insights into the surface kinetics of metallic Mg electrodes.


INTRODUCTION
Rechargeable batteries are of outstanding importance in achieving global climate goals as they play a key role in electronic applications, in the mobility sector for the electrification of cars, or stationary grid energy storage systems for surplus wind and solar power. The most prominent representatives are Lithium-ion batteries (LIBs) [1][2][3][4][5] , which are facing controversy due to the use of critical raw materials such as lithium, cobalt, nickel, or graphite. [6] Metallic lithium anodes have a high theoretical capacity (3861 mAh·g -1 ), the lowest negative electrochemical potential (-3.04 V vs. standard hydrogen electrode (SHE)), and a low density (0.59 g·cm -3 ), however, are prone to form dendrites. [7] Dendrites are needle-shaped structures leading to short circuits, thermal runaways, and even battery failures, although promising approaches such as the electrostatic shielding of the lithium surface through additives have been proposed to tackle this problem. [8] Due to these safety concerns, graphite with intercalated lithium ions is still used as anode material in commercial LIBs, reducing its theoretical capacity by roughly one order of magnitude to 372 mAh·g -1 in present day commercial LIBs. [9] Therefore, rechargeable multivalent magnesium-ion batteries (MIBs) have emerged as a promising alternative, mainly due to the bivalency of the magnesium atom, which allows it to carry twice the charge resulting in a higher volumetric capacity (3833 mAh·cm -3 vs. 777 mAh·cm -3 of graphite anodes in commercial LIBs). [10] Additionally, magnesium has a low electrochemical potential (-2.37 V vs. SHE), is highly abundant, of low cost, and is environmentally friendly. These benefits are reflected in the exponentially increasing number of publications that have recently been dedicated to the topic of MIBs, as shown in Figure 1 from the publication by Li et al. [11] On the other hand, a great deal of effort is required to find suitable electrolytes that allow reversible magnesium deposition while at the same time providing a broad electrochemical window. [9,12] Conventional carbonate electrolytes cannot be used for rechargeable MIBs because the solid electrolyte interface (SEI) formed by decomposition products of the electrolyte is not permeable to magnesium-, in contrast to lithium-cations. [10,13] Proposed solutions are electrolytes based on Grignard reagents, Mg(TFSI)2 salts or ionic liquids, as well as the use of polymeric interlayers as artificial SEI. [13,14] Nevertheless, the structure and growth mechanism of SEI formation remains poorly understood due to its high complexity, highlighting the need for adequate theoretical modeling on different time and length scales. [15] The second controversial issue is whether dendrite formation is possible on metallic magnesium anodes. On the one hand, the tendency of occurrence is significantly lower than for the lithium equivalent, [16] yet dendrites have been demonstrated under harsh deposition conditions, as shown by the Banerjee group. [17] As a result, there have been requests to conduct experiments to determine critical current densities [18] and to define precise voltage windows [19] for dendrite-free regions. In a recently published study, the critical overpotential of dendrite formation on lithium, sodium, and magnesium was calculated using a grandcanonical density functional theory (DFT) approach, indicating a combination of high surface tension, low surface capacitance, and low potential of zero charge (PZC) as indicators for mitigated dendrite growth. [20] Another explanation for enhanced dendrite growth on lithium compared to magnesium surfaces includes the accumulation of negative excess charges at the peaks of protrusions of the lithium electrodes inducing a strong electric field that attracts further lithium cations. [21] In addition, lower step-down and terrace self-diffusion barriers of magnesium compared to lithium are considered responsible for the reduced dendrite growth. [22,23] In this work, periodic DFT calculations have been performed to determine the atomistic properties of magnesium and to explain the initial stages of surface growth. This technique has previously been applied for lithium [24] and sodium [25] in our group. For this purpose, the bulk, surface, and adsorption properties of magnesium were first calculated. Based on the obtained surface energies, the corresponding equilibrium shape of a magnesium crystal was determined within the framework of the Wulff construction. [26] Subsequently, a wide variety of two-and three-dimensional diffusion processes on perfect and imperfect Mg(0001) and Mg(101 ! 1) surfaces were investigated and analyzed concerning surface growth phenomena.
Finally, the obtained DFT results provide a data set for parametrizing a force field for metallic Mg or performing corresponding kinetic Monte-Carlo studies. [27,28]

METHODOLOGY
Periodic DFT calculations were performed with the plane-wave based Vienna ab initio simulation package (VASP). [29,30] The core electrons were described using the projector augmented wave (PAW) [31] method of Blöchl as implemented in VASP. [32] Exchange correlation effects were described within the generalized gradient approximation (GGA) of Perdew, Burke, and Enzerhof (PBE). [33] Additionally, the Bayesian error estimation functional with van der Waals correlation (BEEF-vdW) was applied to determine the standard deviation of the investigated quantities. [34] A plane-wave cutoff energy of 400 eV was employed after detailed convergence studies (SI Figure 1), where Mg 2s electrons were described as valence electrons. Following the scheme of Monkhorst and Pack, a k-point mesh density of at least 0.14 Å -1 was applied for all total energy calculations. [35] Thermal smearing of one-electron states was allowed using the Gaussian smearing method (σ = 0.1 eV) to determine the partial occupancies of each orbital and acquire faster convergence with respect to the number of kpoints. The electronic self-consistent field (SCF) was considered converged when the total energy difference was less than 10 -5 eV, and the norms of all the forces were smaller than 10 -3 eV·Å -1 . Surface energies were calculated on symmetrical (1×1) cells with a minimum slab thickness of 30 Å and a minimum vacuum region of 20 Å (SI Table I). A single atom in the middle of the slab was fixed to avoid net translations. [24] Adsorption energies, dimer interaction, diffusion barriers, and frequencies were generally calculated on converged slabs with 6 surface layers, where the two uttermost layers were fixed to emulate bulk properties. If it became necessary, vicinal slabs were used. The exact surface characteristics are summarised in SI Table II. The minimum energy path (MEP) was determined with the climbing image nudged elastic band (CI-NEB) method as implemented in the transition state tools for VASP (VTST) with generally 7 images between the initial states separated by a spring constant of 5.0 eV·Å -2 . [36,37] As an initial guess for the MEP the image-dependent pair potential (IDPP) was applied. [38,39] To confirm the transition state and to determine rate constants, we calculated vibrational frequencies with the dynamical matrix method from the VTST package.
For this purpose, the electronic and ionic convergence criteria were increased to 10 -8 eV and 10 -8 eV·Å -1 , respectively, and the diffusing as well as the nearest neighbor atoms along the diffusion pathway were symmetrically displaced by 0.007 Å in each spatial direction.
The cohesive energy !"# is defined as the energy gain an isolated atom receives when it is embedded in a certain crystal structure. It was calculated by subtracting the energies of the isolated atoms $%"& from the energy of the bulk crystal structure '()* , divided by the total number of atoms in the bulk structure: Bulk moduli were calculated based on the jellium equation of state. [40] In order to be able to compare the theoretical cohesive energy with the experimental reference, the zero-point vibrational energy +,-. was estimated as follows: [41] +,-. = 9 8 / 0 (2) / corresponds to the Boltzmann constant and 0 to the Debye temperature, which is 400 K for magnesium. [41] The surface energy is defined as the surface excess free energy per unit area and was calculated by subtracting a multiple of the bulk energy • '()* from the energy of the slab 1)$' , divided by the area of the surface times two (because of the symmetric slab configuration): All surface energies were emloyed to carry out the Wulff construction. [26] The area fractions of the facets were calculated using the Python package WulffPack for a nanocrystal with 5000 atoms. [42] Adsorption energies $2 were obtained by subtracting the energy of the pure surface 1)$' plus the energy of an isolated magnesium atom 345$%"& from the energy of the relaxed system %"% containing the slab with the adsorbed atom: The interaction energy 67% was calculated by subtracting two times the adsorption energy of a monomer $2 &"7"&89 from the adsorption energy of a dimer $2 26&89 : The Ehrlich-Schwoebel barrier .: was determined for all step-down diffusion processes. [43,44], [45] It corresponds to the additional energy barrier that a diffusing atom has to overcome when descending a step and was calculated by subtracting the terrace self-diffusion barrier $ ;"9 (>899$!8) from the barrier of the step-down process $ ;"9 (:%8@52"A7) : .: = $ ;"9 (:%8@52"A7) − $ ;"9 (>899$!8) The reaction rate @C> at room temperature ( = 293.15 ) was calculated for each diffusion process by means of transition state theory as given by the Arrhenius equation: where the pre-exponential factor was determined based on the Einstein approximation [46] with the vibrational frequencies of the diffusing adatom and the nearest surface atoms along the migration pathway for the initial and the transition state. $ equals the activation energy and / is the Boltzmann constant.
The activation temperature $ above which a process runs at the rate Γ was calculated with the equation used by Bogicevic et al. [47] Above $ , a process is considered to be 'activated' following the formula: Since the rate Γ depends on the experimental growth rate, Γ = 1 s 5D was set, which is valid for a deposition rate of about 0.001 − 0.1 ML • s 5D . As indicated, $ corresponds to the activation energy, / to the Boltzmann constant, and to the pre-exponential factor. For some ) with very small activation barriers (<0.02 eV), the pre-exponential factor could not be determined due to imaginary frequencies in the TS. In these cases, a value of = 5.0 • 10 DE s 5D was assumed, which corresponds to the average value of the pre-exponential factor, usually in the range between 10 DE − 10 DF s 5D .

a. Bulk properties
The bulk properties of magnesium, including lattice constants, binding energies, and bulk moduli, were calculated for the following crystal structures with the PBE and the BEEF-vdW functional: hexagonal-closed-packed (hcp), double-hexagonal-closed-packed (dhcp), bodycentered-cubic (bcc), face-centered-cubic (fcc), simple cubic (sc), β-tungsten (a15) and diamond (dia). The calculated values for both functionals for all crystal structures and a comparison with experimentally obtained data are shown in Table I. Magnesium crystallizes in accordance with experimental studies in a hcp structure under ambient conditions. The calculated cohesive energy of 1.50 eV·atom -1 from the PBE functional is in almost perfect agreement with the experimental value (Table I) and with other theoretical studies (SI Table   III). The BEEF-vdW functional slightly underestimates the cohesive energy, but it is within its estimated error. Moriarty et al. first predicted a phase transformation from the hcp to the bcc structure at around 50 GPa from generalized pseudopotential theory [48][49][50] , which Olijnyk and Holzapfel later confirmed by X-ray diffraction (XRD). [51] This is contradicted by the study of Erandonea et al., who assumed the dhcp instead of the bcc structure at high pressures and temperatures. [52] Moriarty et al. also predicted a second phase transformation from the bcc to the fcc structure at even higher pressures [49] , which could not be observed experimentally. [53] Table I. Calculated and experimentally observed physical constants for different bulk phases of magnesium. The lattice constants and are given in Å, bulk moduli in GPa, and the cohesive energy in eV·atom -1 . For the BEEF-vdW functional, the calculated standard deviation of the cohesive energy is given in parenthesis. The experimental results were taken from Ref. [54]. The influence of the zero-point vibrational energy on the experimental values is given in parentheses. The values were adapted from Alchagirov et al. for , , and and calculated with equation (2) for . [41] Bulk

b. Surface properties
The surface energy is regarded as a measure for the stability of a surface and determines the equilibrium shape of a crystal. [55] In agreement with literature [56][57][58] , Mg(0001) is the most stable crystal surface, followed by Mg(101 ! 0)A and Mg(101 ! 1), which differ in our calculations by less than 1 meV·Å -2 (Table II). Mg(0001) lies in the basal plane of the hcp unit cell, which is characterized by a particularly high packing density, resulting in thermodynamic stabilization. [54] At this point, we refer to the work of Tang et al., who established a model to predict the relative stability of magnesium surfaces based on the number of broken basal and non-basal bonds. [59] They also showed that Friedel oscillations are especially present on lowindex Mg surfaces, influencing surface relaxation. On the other hand, on high index surfaces, relaxation is driven by charge depletion and charge smoothing effects. In an embedded atom method (EAM) potential study, Mg (101 ! 1) was considered to be the second most stable surface, which is not the case according to our DFT calculations. [60] Caution is required for  Interestingly, and as already shown in literature [56,57] , the Wulff construction ( Figure 1 Table III. Other publications also assign Mg(112 ! 0) area fractions in the Wulff construction, which was not observed in our case. [60,62] It is important to mention that even over a wide potential range, the shape of the magnesium crystal behaves almost constant, which suggests that in the working potential range of magnesium batteries, our calculated surface area fractions are still valid. [63] The calculated area fractions emphasize the need to study the diffusion properties of all three surfaces collectively, rather than looking at one surface termination individually.   . [56] We also want to reference two EAM potentials for magnesium by Liu et al. [64] and Sun et al. [65] However, in the potential of Liu et al., only the bridge position on Mg(0001) is stable, which is not the case according to our DFT calculations. [60] We also want to point to the section Diffusion barriers, which further clarifies why fcc sites are energetically favored over hcp sites, although this contradicts the nature of the hcp lattice, which should follow the ABABAB-stacking sequence.
All calculated adsorption energies are consistent with previously published theoretical results (on PBE level). [22,56,66,67] The error range for the BEEF-vdW functional varies between 0.10 and 0.31 eV, indicating a strong dependence on the chosen functional. However, all PBE results are within the BEEF-vdW standard deviation, and the energy ratios are equal for both functionals.

Adsorption site
This work Other works

Dimer and trimer self-diffusion on Mg(0001) and Mg(101 ! 1)
The formation of polyatomic clusters represents the first stage of island formation. Here the question arises if attractive or repulsive interactions between adatoms are at play, resulting in a uniform distribution of the adsorbates or the formation of dimer, trimer, or cluster structures.
Our calculations reveal only negative or negligible (|Eint| < 0.03 eV) interaction energies between two Mg adsorbates and an interaction sequence given by Eint, Table VI, Table VII and Table VIII. The corresponding dimer sites are given in Figure 3. The differences in interaction energies may be related directly to the number of nearest surface atoms of the adsorbates following the opposite order (3Mg(0001) On Mg (0001)

This work
Other works [22], [56] [22] * not stable; Transition into fcc-hcp configuration.  [56] 0 ⟷ 2 0.01 0.01 (0.02) *Since the adsorption site of atom 4 is not stable outside of the dimer configuration, the adsorption energy for this atom was assumed to be equivalent to an fcc position on Mg(0001) to calculate the interaction energy. As the interaction energy studies ( Figure 5, Table IX      Step, kink and corner self-diffusion on Mg(0001) and Mg(101 ! 1) After the formation of initial polyatomic clusters, we would like to discuss further surface growth, especially island formation, by studying additional processes of the terrace-step-kink (TSK) model (step-edge, step-vacancy, kink, inner-corner, and outer-corner). [70] The TSK model describes the diffusion processes taking place during surface growth and discusses several growth mechanisms. [71] As described previously, islands on Mg(0001) are more stable on hcp sites starting from a heptamer. Thus, only processes at hcp steps, kinks, and corners were investigated, which consist of more than 7 atoms. A detailed overview of the processes studied may be found in Table XI, Figure 8, and Figure 10 for Mg(0001) and Table XII, and

Mg(0001) Path Pathway
Step-edge Initially, we want to discuss the process of a free atom approaching the step-edge on Mg(0001). As we have already seen with similar dimer and trimer processes, the forward barrier is negligibly small, but the detachment barrier is many times higher [S0/3 ↔ S1 On Mg(101 ! 1), it is noticeable that an approaching atom hast to overcome higher forward barriers of at least 0.14 up to 0.44 eV, depending on the direction of the step. However, this behavior is expected since the barriers for terrace self-diffusion on Mg(101 ! 1) are already at least 0.30 eV. The ratio between forward and reverse barriers is between 1-4 times, significantly lower than for Mg(0001) (25-40 times).
Step-edge diffusion on Mg (0001) 1)). Fractal islands are formed when an adatom arrives at an edge, but relaxation to a more favorable position is hindered. [71] Looking at the diffusion along the steps on Mg(101 ! 1), it is noticeable that depending on the The reason might be due to different coordinations of the diffusing atom in the initial/final sites (hollow vs. fcc-like sites). Step-edge: Step-edge: Step-vacancy: Step-vacancy: (101 & 1). The values are given in eV.

Ehrlich-Schwoebel-barrier and upper step self-diffusion on Mg(0001) and Mg(101 ! 1)
So far, we have only considered 2-dimensional processes. However, to get a complete picture of surface growth, 3-dimensional processes such as upper-step diffusion followed by its descent are particularly important. The activation energies of the respective processes (stepdown, step-down (dimer), upper-step) as well as the calculated Ehrlich-Schwoebel barriers for Mg(0001) and Mg(101 ! 1) are summarized in Table XIII and Table XIV. The schematic illustrations of the processes are given in Figure 11 and Figure 12. [ " ] Ex. Ehrlich-Schwoebel barriers indicating a uniform and smooth growth due to a high rate of interlayer mass transport. [71] Likewise, the descent is clearly preferred to the ascent following the observed increase in coordination partners. All step-down barriers are in agreement with the literature. [23,62] Going one step further, we wanted to investigate the influence of a second adatom on the is fully complete. [71] A further interesting case is the exchange process at the [1 ! 1 ! 20]-step. In this process, the barriers for ascent and descent are very close together. However, it is important to note that the second atom follows the diffusing atom into the fcc position above the exchange atom during descent. This arrangement maintains the dimer conformation until beyond the TS (SI Figure 11) and lowers the activation energy of the process to half of the corresponding exchange process at the [112 ! 0]-step. (101 & 1). The values are given in eV.

[12 & 10]B-step
Step-down: [ " ] e. Activation temperature of diffusion processes on Mg(0001) and Mg(101 ! 1) In addition to the diffusion barriers, the pre-exponential factors of the Arrhenius equations were calculated via the Einstein approximation to determine the activation temperatures of the processes studied on Mg(0001) and Mg(101 ! 1). The obtained activation temperatures at which the respective diffusion processes start to contribute to surface growth are shown in Figure 13 and Figure 14.
On Mg(0001), processes such as terrace self-diffusion, the formation of dimers and trimers, and the step-down process do not have significant barriers and are already activated at temperatures below 25 K. It must be noted that dimer and trimer formation is highly favored compared to the respective separation processes due to a very high interaction energy. At this point, we would like to address the hypothesis that lower terrace self-diffusion and step-down barriers on Mg(0001) vs. Li(100) are responsible for reduced dendrite formation. [23] We were able to confirm lower barriers for Mg(0001). However, the corresponding processes for metallic lithium become activated at 16 and 43 K, respectively. [24] At room temperature, these processes should run almost equally and have no significant influence on dendrite formation.
We believe that to adequately explain dendrite growth, a holistic and more exhaustive model needs to be developed that considers the influence of the SEI interface, the stability of different surface terminations, the diffusion barriers on the respective surfaces, the applied potential, the deposition conditions, and the electrolyte used. It remains an open question to what extent this is already possible with the current theoretical tools. However, the data presented in this work could serve as a reference for future studies where these additional effects are successively accounted for.
At a temperature below 100 K on Mg(0001), a variety of further diffusion processes become active such as dimer and trimer propagation, diffusion along an edge, the formation of 60° and 120° corners, switching sides at 240° and 300° corners and kink incorporation. Especially the activated corner-crossing barriers contribute to the growth of compact rather than fractal islands. [71] Interestingly, both the dimer step-down and step-up processes become possible in this temperature range, contrary to uniform and smooth surface growth.
On Mg(101 ! 1), we found that the activation temperatures for diffusion processes are generally higher. Below 100 K, it is only possible for an adatom to step-down, form a dimer, or diffuse towards an edge. Terrace diffusion within and across a channel and diffusion along an edge on Mg(101 ! 1) becomes possible at elevated temperatures above 100 K.
In the same temperature range, processes in which the coordination number of the diffusing atom is reduced become available on Mg(0001), such as in the 120° corner and kink breaking processes. Above 200 K, the dimer and trimer separation processes, edge evaporation, and 60° corner separation are also enabled, as well as the step-up process. Despite that, the statement still applies that bond breakages are extremely unlikely for magnesium and that the increase of the coordination number always proceeds preferentially. For example, kink evaporation and the formation of a step-vacancy are not possible at room temperature for Mg(0001). Nevertheless, activated evaporation and separation processes are necessary for island ripening (redistribution of mass). [71] We assume many small islands form in the lowtemperature range, and new adsorbates join existing ones. At higher temperatures, island ripening should begin, and smaller islands dissolve at the expense of larger islands.
On Mg(101 ! 1), all processes studied can occur at ambient conditions. In the temperature range between 200 and 300 K, dimer and trimer separation, edge evaporation, and the exchange step-up process are activated. In addition, concerted dimer propagation proceeds at about 200 K.  At this point, we would like to refer to the Supporting Information, which summarizes all calculated data. SI Table V and SI Table VI show the pre-exponential factors, activation energies, room-temperature rate constants, and activation temperatures for the terrace selfdiffusion processes for the PBE and BEEF-vdW functionals, respectively. Corresponding schematic representations and energy profiles are provided in SI Figure 4. In addition, overviews of all investigated diffusion processes on Mg(0001) are available in SI Figure 6 and on Mg(101 ! 1) in SI Figure 8. The corresponding pre-exponential factors, activation energies, room temperature rate constants, activation temperatures, and energy profiles may be found in SI Table VIII

CONCLUSION
This work aimed to determine and summarize the atomistic properties of magnesium and discuss the initial stages of surface growth in a possible magnesium-ion battery. Therefore, bulk, surface, adsorption, and diffusion properties of magnesium were calculated at the level of density functional theory. We are aware that in a real battery environment, the operation conditions, the potential, the electrolyte, charge-and discharge products, and especially the SEI environment might have a tremendous impact on the atomistic properties of magnesium.
Nevertheless, we still believe that our studies already provide informative insights into the actual processes taking place. The reason for the stacking fault might be Friedel oscillations causing an increased charge density pocket at the fcc site for smaller magnesium clusters. [67] Terrace self-diffusion of single atoms runs almost barrier-free on Mg(0001) but with higher activation energies on Mg(101 ! 1). When adsorbates collide on the surfaces, initial dimer, trimer, and cluster structures are formed, which is highly favored due to high attractive interactions between the adsorbates. The small cluster structures can be seen as seeds for further island growth, which shape on Mg(0001) could be either triangular with only <1 ! 1 ! 20>steps as proposed by Jäckle et al. [23] or hexagonal with alternating <112 ! 0>-and <1 ! 1 ! 20>-steps at elevated temperatures [72] . However, the island pattern should be compact rather than fractal due to low step-edge and corner crossing barriers. [71] On Mg(101 ! 1), the diffusion barriers are in general higher compared to Mg(0001), but the ratios for bond formation and breakage are indicating a uniform and smooth surface growth due to a high rate of interlayer mass transport at low deposition rates. [71] However, a second adatom significantly increases the descent barriers for hopping and exchange processes, fostering the aggregation of new islands on top of existing islands if the adatom meets another adatom faster than it reaches the step-edge.