Theory and Atomistic Simulation of Electrodeposition

Abstract

Electrodeposition is a fundamental process in electrochemistry, and has applications in numerous industries, such as corrosion protection, decorative finishing, energy storage, catalysis, and electronics. While there is a long history of using electrodeposition, its application for controlled nanostructure growth is limited. The establishment of an atomic-scale understanding of the electrodeposition process and dynamics is crucial to enable the controlled fabrication of metal nanoparticles and other nanostructures. Significant advancements in molecular simulation capabilities and the electronic structure theory of electrified solid-liquid interfaces bring theory closer to realistic applications, but a gap remains between realistic applications, theoretical understanding of dynamics, and atomistic simulation. In this review we briefly summarize the current state-of-the-art computational techniques available for the simulation of electrodeposition and electrochemical growth on surfaces, and identify the remaining open challenges.

Warwick-Chemistry]Department of Chemistry, University of Warwick, Coventry, CV4 7AL, United Kingdom [Nottingham]School of Chemistry, University of Nottingham, Nottingham, NG7 2RD, United Kingdom Warwick-Chemistry]Department of Chemistry, University of Warwick, Coventry, CV4 7AL, United Kingdom [Warwick-Physics]Department of Physics, University of Warwick, Coventry, CV4 7AL, United Kingdom

1 Introduction

Electrodeposition is the formation of solid structures on the surface of an electrode when an electrochemical potential is applied, and is a viable nanofabrication process alongside more established methods such as nanoimprint lithography ^{1, 2}, pH-driven precipitation ^{3, 4}, and directed assembly ^{5, 6, 7, 8, 9, 10, 11, 12}. Metal electrodeposition is a fundamental electrochemical process, and is an important method with many applications such as carbon dioxide reduction catalysis ¹³, water splitting ¹⁴, fuel cell applications ¹⁵, and materials for energy storage and conversion. Metal electrodeposition is also inherently used in various industrial applications, such as electroplating ^{16, 17, 18}, electrowinning ^{19, 20, 21}, and electrocatalysis ^{22, 23, 24}. The key challenge in metal electrodeposition is to control the structure, size, and stability of surface-adsorbed nanostructures on an atomistic scale, which in turn define the reactivity and electrochemical properties of the resulting materials. The establishment of an atomic-scale understanding of the electrodeposition process and dynamics is thus crucial to enable the controlled fabrication of metal nanostructures.

Experimental techniques to track, characterize and harness metal electrodeposition have vastly improved over recent years ^{25, 26, 27, 28, 29, 30}. The complementary use of microscopy approaches, surface spectroscopy methods, and electrochemical analysis provide unprecedented resolution at the nanoscale and, to a more limited extent, resolution in the time domain ²⁶. Simultaneously significant advancements have been made in molecular simulation capabilities and the electronic structure theory of electrified solid-liquid interfaces ^{31, 32}. Yet a large gap between realistic applications, theoretical understanding of dynamics, and atomistic simulation remains. Challenges that require further understanding include electrode-electrolyte interactions, electronucleation and growth mechanisms, and reaction rates and kinetics. Both theory and experiment face challenges when it comes to bridging this gap and reaching an atomic-level understanding of electrodeposition. Theoretical and computational studies must be able to simulate realistic models capable of replicating experimental conditions, accounting for factors such as the electrochemical potential and surface heterogeneity. On the other hand, model experimental studies should ideally be conducted under well-defined and idealized conditions (e.g. atomically-flat electrode interfaces and well-purified electrolytes) to allow for atomistic simulations and theoretical analyses to be applied ³³. The synergy between experiment and simulation has the potential to deeply enrich the field, as modeling methods can be refined once information about atomic structure is attained from experiment, while simulations can be used to make predictions that experiments can validate ³³. Given the rapid advancement in simulation methodologies and the increasing importance of metal electrodeposition, it is now timely to comprehensively review the vast range of modeling methods available for atomistic simulations. We hope that this review will guide future efforts in this field and help synergize theory and experiment.

Existing reviews tend to focus on specific aspects of computational electrochemistry, such as electron transfer processes ^{34, 35}, modelling methods ^{36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46}, solvation and solid-liquid interfaces ^{47, 48, 49, 50}, and the electrochemical double layer ^{51, 52}. A useful collection of computational electrochemistry reviews is presented in Koper et al. ⁵³, while Gamburg and Zangari ⁵⁴ is a recommended resource that covers the theory and practice of metal electrodeposition. However, existing reviews do not focus on metal electrodeposition in particular, nor do they discuss all the various aspects that need to be considered in atomistic simulations of complex electrochemical processes and environments.

In this review, we summarize the key concepts and aspects of electrodeposition, as well as the various state-of-the-art computational techniques that are available for the atomistic simulation of electrodeposition. Finally, we discuss open questions and challenges in the field.

2 Principles of Electrodeposition

The process of metal electrodeposition can be described in four steps (Figure 1): diffusion of metal cations through the solvent after the application of an electric current; adsorption of metal cations at the cathode surface via electron transfer reactions; migration of metal adatoms along the cathode surface; and the nucleation of larger metal nanostructures on the cathode surface.

Figure 1: A simplified diagram showing the process of metal electrodeposition onto a cathode surface. First, solvated metal cations electrodiffuse through the solvent (shown as water) towards the cathode (step 1). Once close enough to the cathode, electron transfer will occur that reduces the metal cation, resulting in electrosorption (step 2). The metal atom can then migrate along the cathode surface (step 3), either coming to rest at isolated sites on the cathode surface or coalescing with other metal atoms, resulting in the electronucleation of nanostructures (step 4). The nucleation process is in competition with dissolution of nanostructures, which can result in further surface migration. Solvated cations and anions are shown as ‘+’ and ‘$-$’, respectively, and hydrogen and oxygen atoms are colored white and red, respectively.

2.1 Electrodiffusion of Solvated Ions

There are several mathematical relationships that can be used to calculate useful macroscopic quantities in the context of electrodeposition. The Nernst-Planck equation, shown in Equation (1), is a generalization of Fick’s laws ⁵⁵ and describes the flux, $J_{i}$, of a charged particle $i$ (such as a precursor cation), due to diffusive and electrostatic forces.

$$J_{i}=-D_{i}\nabla c_{i}-z_{i}Fc_{i}\nabla\phi$$ (1)

where $D_{i}$ is the diffusion coefficient of $i$, $c_{i}$ is the concentration of $i$, $z_{i}$ is the charge of $i$, $F$ is the Faraday constant, and $\phi$ is the electric potential.

The distribution of charged particles within the system will also have an effect on the electric potential, and this relationship can be formulated via the Poisson equation, as shown in Equation (2).

$$\nabla^{2}\phi=-\dfrac{\rho}{\varepsilon}$$ (2)

where $\rho$ is the charge density and $\varepsilon$ is the permittivity of the medium.

The Gouy-Chapman equation, shown in Equation (3), can be used to describe the behavior of ions as they approach a charged surface, such as a cathode.

$$\sinh{\bigg{(}\dfrac{ze\phi}{4k_{\mathrm{B}}T}\bigg{)}}=\dfrac{c_{+}}{c_{-}}$$ (3)

where $z$ is the number of electrons transferred, and $F$ is the Faraday constant, $e$ is the elementary charge, $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ is the temperature, and $c_{+},c_{-}$ are the concentrations of cations and anions, respectively.

While Equations (1)–(3) provide many valuable insights into the diffusion of solvated ions during electrodeposition, they possess some limitations. The Nernst-Planck equation assumes a dilute solution and that the local electrostatic field is not affected by the ionic concentration, which is likely to not hold at high ionic concentrations. Furthermore, ion-ion interactions are neglected within the formulation, which can significantly influence the diffusive behavior of solvated ions. Additionally, Equation (1) was originally formulated for only binary solutions and would require modifications in order to be applicable for more complicated and realistic systems comprising many ionic species. Despite these limitations, Equation (1) has been widely used to study the time evolution and transport of charged species ^{56, 57, 58, 59, 60, 61} in experiments.

The Poisson equation assumes solvated ions are stationary and that their motion is not affected by changes in the electrostatic potential. In addition, the distribution of the charge density, $\rho$, in Equation (2) is assumed to be continuous. Both of these assumptions will not always be accurate in realistic systems, especially with charged metal cations interacting with solvent molecules. Despite this shortcoming, Equation (2) has been used to elucidate various processes during metal electrodeposition, including electrodiffusion, ionic migration and electroconvection ^{62, 63, 64, 65}.

Finally, the Gouy-Chapman equation possesses two major limitations. Firstly, all ions are assumed to be point charges within the system and have the same size. This assumption can result in inaccuracies with solvated metal cations, where particle size can play a part in the final cationic distribution near the cathode surface. Secondly, Equation (3) assumes that the charge at the cathode surface is uniformly distributed in a single plane, and thus neglects more complex situations such as surface heterogeneity and already-adsorbed species.

The limitations of these approaches can be mitigated and/or complemented by atomistic simulation methods, which can be used to provide detailed insights into dynamics and interactions that are not captured by such macroscopic models.

2.2 Electrosorption of Metal Cations

Faraday’s laws of electrolysis ⁶⁶, summarized in Equation (4), relate the total electrodeposited metal mass, $m$, with the net charge, $Q$, within an electrochemical cell:

$$m=\dfrac{QM}{zF}$$ (4)

where $M$ is the molar mass. Equation (4) is widely used to calculate various quantities, such as the amount of electrodeposited metal or the time required to achieve a given amount of electrodeposited metal. Faraday’s laws of electrolysis can be manipulated to describe a range of processes, including electroplating, metal corrosion and electrowinning ⁶⁷, and have been widely applied to experiments ^{68, 17}. However, Faraday’s law does not take into account the interaction between the metal adsorbate and the underlying surface, and only relates $m$ with $Q$. It therefore assumes that any electrodeposition that occurs is irreversible and effective.

From a more atomistic point of view, the electrosorption of species onto electrode surfaces is typically described using the same theories that are used to treat adsorption in the gas phase, such as $d$-band theory ⁶⁹. During the adsorption process, bonding and antibonding hybrid orbitals form between the $d$-band of the surface and the valence orbital of the adsorbate. The strength of the binding interaction is determined by the filling of antibonding states ^{69, 70, 71}. However, unlike gas-phase chemistry, where this filling is determined by the number of electrons, at a metal surface, the filling is governed by the energy of the antibonding states with respect to the Fermi level ⁷⁰ (the energy level at which the probability of an electronic state being occupied by an electron is 50%). As the energies of antibonding states will lie above those of $d$-states, the energy at which the center of the $d$-state distribution lies above the Fermi level can be used to gauge the strength of the binding interaction ⁷⁰. The higher the energy of the $d$-states with respect to the Fermi level, the higher in energy the antibonding states are and therefore the stronger the binding interaction ^{70, 72}. However, $d$-band theory only considers adsorbate-adsorbent interactions and does not take into account further interactions with the electrolyte or solvent molecules ⁷²; it is thus too simplistic to fully describe electrosorption dynamics ^{73, 74}.

To accurately model metal electrosorption, one must first consider the fundamental differences between adsorption at solid-gas and solid-liquid interfaces. The key difference, in the context of metal electrodeposition, is that once the electrode is negatively polarized, the hydrogen atoms in water molecules will be attracted to the electrode surface (due to their positive dipoles), and these molecules need to be displaced from the cathode surface before metal cations can adsorb. Adsorption isotherms are commonly used to describe adsorption and desorption processes, and determining which adsorption isotherm to use is dependent on the specific characteristics of the system and what assumptions are made by the isotherm for the intended purpose. The simplest example is the Langmuir adsorption isotherm ⁷⁵, given by Equation (5):

$$\dfrac{\theta}{1-\theta}=\exp\bigg{(}-\dfrac{\Delta G^{\minuso\prime}_{\mathrm{ads}}}{RT}\bigg{)}a_{\mathrm{ads}}$$ (5)

where $\theta$ is the fractional surface coverage (by metal atoms), $\Delta G^{\minuso\prime}_{\mathrm{ads}}$ is the standard electrochemical free energy of adsorption, $R$ is the ideal gas constant, $T$ is the temperature, and $a_{\mathrm{ads}}$ is the activity of the adsorbate in bulk solution. Equation (5) can be rewritten to take into account the applied electrode potential, $E$:

$$\dfrac{\theta}{1-\theta}=\exp\bigg{(}-\dfrac{\Delta G^{\minuso}_{\mathrm{ads}}+zF(E-E^{\minuso})}{RT}\bigg{)}a_{\mathrm{ads}}$$ (6)

where $\Delta G^{\minuso}_{\mathrm{ads}}$ is the standard free energy of adsorption and $E^{\minuso}$ is the standard potential of the adsorption reaction ⁷⁶. While the Langmuir adsorption isotherm provides a first approximation of electrosorption behavior, there are several inherent assumptions that limit its applicability. Firstly, it does not take into account the displacement of solvent molecules, which is a key component of electrosorption. It also assumes the adsorbent surface to be homogeneous, with all adsorbates being assumed to adsorb as a monolayer. Furthermore, all adsorption sites are assumed to be energetically equivalent and interactions between adjacent adsorbed species are neglected. Several extensions to the Langmuir adsorption isotherm exist to mitigate these assumptions ^{77, 78, 76}. For example, a Flory-Huggins-type isotherm can be used to also describe the displacement of solvent molecules during electrosorption ^{79, 80, 76}, while the Frumkin adsorption isotherm ⁷⁷ can be used to also account for interactions between adjacent adsorbed species, and the Freundlich adsorption isotherm ⁷⁸ is more appropriate for describing multilayer adsorption on heterogeneous surfaces.

Adsorption isotherms are widely used in experiments ^{81, 82, 83, 84, 85} to determine parameters such as the maximum adsorption capacity or the relationship between the quantity of particles in a metal deposit and their concentration in an electrolyte. While powerful tools, adsorption isotherms require prior knowledge of the characteristics of the system in some manner and do not take into account electrosorption on an atomistic scale. Atomistic simulations can, however, be used to address these limitations by including factors that are excluded within isotherms.

2.3 Electron Transfer during Electrodeposition

Electron transfer reactions are the process by which metal cations adsorb onto a given cathode surface. Such reactions are termed to be either inner-sphere or outer-sphere. Inner-sphere electron transfer is when the reaction occurs via a strong electronic interaction, such as a covalent bond, and the reactants become connected by a chemical bridge. In contrast, outer-sphere electron transfer (OS-ET) occurs between two ‘unconnected’ reactants, and the electron has to move from one reactant to the other over some space, typically at least a solvent layer from the cathode surface ⁸⁶ (in the context of electrodeposition). With transition metals, the mechanism of electron transfer for transition metals can be generally assumed to be outer-sphere as, in solution, transition metal cations can form complex coordination compounds with ligands, such as water molecules. The presence of such ligands can create steric hindrance as well as tightly bind the inner coordination sphere around the metal cation, which makes inner-sphere electron transfer less favorable. However, exceptions can arise depending on factors such as the choice of solvent.

Marcus theory is commonly used to describe the rate of electron transfer reactions ⁸⁷. For multivalent ions, Marcus theory suggests that the simultaneous transfer of multiple electrons is improbable, which implies that the cathodic deposition of multivalent metal cations should occur in a sequence of one-electron steps ^{87, 88, 89}, with the adsorption reaction being the final one-electron step i.e. for a metal element, M: $\textnormal{M}^{z+}+\textnormal{e}^{-}\rightarrow\textnormal{M}^{(z-1)+}$ until the neutral atom is reached (i.e. $z-1=0$), rather than $\textnormal{M}^{z+}+z\textnormal{e}^{-}\rightarrow\textnormal{M}$. Marcus theory has since been combined with a Newns-Anderson Hamiltonian ^{90, 91} (vide infra) to describe electrochemical processes by Schmickler and coworkers ^{88, 89}, and has been used to deduce that small monovalent metal cations, such as Ag⁺, are able to get close to the electrode surface without losing solvation energy as they fit into the water structure well, unlike larger multivalent cations ^{88, 92}. As the valency of metal cations increases, the potentials of mean force sharply increase on approach to the electrode ^{88, 89}. When divalent cations approach the electrode surface, they shed their secondary solvation shells, causing their free energy to rise. This means a close approach of divalent (and multivalent) cations to the electrode surface is energetically unfavorable ⁸⁹.

Current density, $j$, is typically used as a quantitative measure of the rate at which electrons are transferred between solvated metal cations and the electrode surface during electrodeposition. Current density is typically used instead of just current as larger electrodes can carry more current. The normalization of current by the surface area of the electrode means processes such as electron transfer can be more easily compared on different electrodes.

The current density associated with a first-order reaction is given by Equation (7).

$$j=zFkc$$ (7)

where $k$ is the rate constant. Within general (electro)chemical rate theory and transition state theory, $k$ obeys the Arrhenius equation ^{93, 94} and is given as Equation (8).

$$k=\kappa\exp{\bigg{(}-\dfrac{\Delta G^{\ddagger}}{k_{\mathrm{B}}T}\bigg{)}}$$ (8)

where $\kappa$ is a potential-dependent pre-exponential factor and $\Delta G^{\ddagger}$ is the free energy barrier (sometimes referred to as the activation energy), which depends on the electrode potential. $\kappa$ accounts for the attempt frequency in transition state theory and may include other effects such as non-adiabatic corrections, solvent dynamics, and nuclear quantum effects ⁹⁵. It should be noted that the frequency factor, $\nu$, is often considered in lieu of $\kappa$; $\nu$ represents the frequency with which a single reactant passes from the reactant side to the product side, and is related to $\kappa$ via Equation (9), where $\delta x$ is the effective area of the region where the reaction occurs ⁹⁴.

$$\kappa\approx\nu\delta x$$ (9)

The reader is directed to He et al. ⁹⁴ for a more detailed review on the importance of the pre-exponential factor in electrochemistry.

Depending on the degree of electronic coupling between the metal nanostructure and the electrode, OS-ET is classified to be either adiabatic or non-adiabatic. Identifying the adiabaticity of OS-ET for different adsorbate-electrode pairs is of fundamental importance in order to optimize the efficiency and mechanism of electron transfer ⁹⁶. In the adiabatic regime, $\kappa$ is independent of the electron tunneling probability between the metal nanostructure and the electrode, and the rate of OS-ET is independent of the electrode material, assuming there exists a sufficiently strong electronic interaction between the adsorbate and the electrode ^{94, 97}. Using the early version of Marcus theory ⁸⁷ and collision theory, the value for the pre-exponential factor was typically set to be $10^{4}<\kappa<10^{5}$ $\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$. Other (and more elaborate) formulae for $\kappa$ have been proposed ^{98, 99, 100, 101, 102} but they typically give rise to the same values in practice ⁹⁴.

In the non-adiabatic regime, $\kappa$ is proportional to the density of electronic states near the Fermi level, as shown in Equation (10) where $V$ is the interaction potential and $\rho$ is the electron density.

$$\kappa\propto\left|V\right|^{2}\rho$$ (10)

Factors that can affect $V$ include electronic-vibrational interactions ¹⁰³, the surrounding environment, and the applied potential ^{104, 105}. For metal electrodes, it is usually enough to use the value at the Fermi level as electrons are mainly exchanged at this level ⁹⁴. OS-ET reactions on pure metal electrodes are often adiabatic ⁹⁶, while some doubt regarding their adiabaticity on other electrodes remains ^{96, 106, 107}.

2.4 Electronucleation

Once metal atoms have adsorbed onto an electrode surface, they can start to coalesce to form larger nanostructures in a process known as electronucleation. Individual metal adatoms may migrate along the electrode surface and coalesce to form metastable nanoclusters. Furthermore, larger, crystalline nanostructures can form if smaller structures rearrange and amalgamate together, while closely-spaced nanoclusters can also dissemble and feed atoms into existing nanostructures. In contrast, isolated metal adatoms that are not part of a nanostructure and do not move along the electrode surface might indicate the presence of point defect sites on the substrate surface ¹⁰⁸. The on-surface dynamics of electronucleation can thus be extremely rich and complex, and atomistic simulations must be able to account for the various thermodynamic and kinetic effects that play a defining part in the size distribution, growth rate, and the rate-determining steps of surface-adsorbed metal nanostructures.

When it comes to modeling the process of electronucleation, both classical and atomistic theories exist to describe the formation of stable nuclei ^{26, 109, 110, 111}. Classical nucleation theory relies on macroscopic physical quantities that are applicable to sufficiently large clusters such that their size, $n$ (number of atoms), can be considered a continuous variable. In this case, the Gibbs free energy of nucleation, $\Delta\overline{G}(n)$, is differentiable ^{110, 111}, which allows the rate of nucleation to be predicted ¹⁰⁹ using Equation (11):

$$\varrho=NZ\overline{J}\exp\bigg{(}-\dfrac{\Delta\overline{G}}{k_{\textnormal{B}}T}\bigg{)},$$ (11)

where $\varrho$ is the nucleation rate, $N$ is the number of nucleation sites, $Z$ is the Zeldovich factor (the probability that a nucleus at the top of the barrier for nucleation will progress to form a new phase rather than dissolve ¹⁰⁹), $\overline{J}$ is the rate at which atoms attach to the nucleus, $k_{\textnormal{B}}$ is the Boltzmann constant and $T$ is the temperature ¹⁰⁹. $Z$ is derived using the assumption that nuclei near the top of the barrier effectively diffuse along the radial axis. By statistical fluctuations, nuclei can either diffusively grow into a larger nucleus that will eventually form a new phase, or can lose atoms and diminish in size.

When it comes to experimental electrochemical measurements, techniques such as chronoamperometry (which is a method where the current decay over time can provide information regarding the rate of cationic reduction and deposition onto an electrode surface) provide mostly macroscopic information, from which nanoscopic behaviour such as nucleation rates can be inferred ^{26, 112, 113}. Such inferences, however, have been found to be inappropriate to describe the initial stages of nucleation where individual atoms and few-atom clusters are present ^{26, 25, 114, 115, 116, 117, 118} due to the assumptions within classical nucleation theory. Classical nucleation theory assumes the holding of the capillarity approximation, which treats the nucleus interior as a bulk incompressible liquid and assumes that the surface tension of a small liquid droplet is equal to the surface tension of a flat surface ^{119, 110, 111}. This has been shown to break down for small systems ^{120, 121, 122}. Furthermore, clusters are assumed to either grow or shrink via single-atom absorption or emission respectively, which places kinetic restrictions on the nucleation pathways ^{110, 111, 119}. This does not hold in reality as entire clusters can merge or fragment, and these kinetic pathways cannot be ignored. While improvements to classical nucleation theory do exist, such as dynamical nucleation theory ¹²³, mean-field kinetic nucleation theory ¹²⁴, coupled flux theory ^{125, 126, 127, 128} and diffuse interface nucleation theory ^{129, 130}, these have mostly been applied to describing the condensation of supersaturated vapors into the liquid phase and crystal nucleation studies rather than investigating metal electronucleation ^{119, 131}. Despite its shortcomings, classical nucleation theory is still a powerful theory and has been shown to be capable of qualitatively capturing nucleation thermodynamics and kinetics for many systems ¹¹⁹.

The growth of metal nanostructures can be controlled either by kinetics, diffusion, or a combination thereof. Diffusion control typically occurs when there is a low concentration of metal cations within aqueous solutions ¹³². The growth of a new-phase hemispherical metal nanostructure of radius $r$, under mixed kinetic-diffusion control, can be expressed using Equation (12). It should be noted that Equation (12) assumes that nucleation proceeds at moderate supersaturation levels and the aforementioned assumptions within classical nucleation theory are valid.

$$j=\dfrac{\exp(\alpha f\eta)-\exp(-\beta f\eta)}{\dfrac{1}{j_{0}}+\dfrac{r\exp(\alpha f\eta)}{zec_{0}D}}$$ (12)

where $j,j_{0}$ are the current density and the exchange current density at the electrolyte/cluster interface respectively; $\alpha,\beta$ are transfer coefficients such that $\alpha+\beta=1$; $\eta$ is the overpotential; $c_{0}$ is the bulk concentration of depositing ions; and $f=\nicefrac{{ze}}{{k_{\mathrm{B}}T}}$, where $e$ is the elementary electric charge ¹³². For purely kinetic-controlled growth, Equation (12) can be rearranged to form the Butler-Volmer equation ^{133, 132}, which can be seen in Equation (13); it can also be rearranged to form an equation for pure diffusion-controlled growth in terms of the diffusion coefficient, as can be seen in Equation (14).

$$j=j_{0}\big{(}\exp({\alpha f\eta})-\exp({-\beta f\eta})\big{)}$$ (13)

$$j=\dfrac{zec_{0}D}{r}\big{(}1-\exp(-f\eta)\big{)}$$ (14)

In contrast, atomistic nucleation theories can be applied to clusters so small that $n$ is no longer continuous, as is the case with first-order phase transitions at high supersaturation ($\Delta\overline{\mu}$) levels resulting in $\Delta\overline{G}(n)$ being non-differentiable ^{110, 111}. Figure 2 shows the intrinsic differences between classical and atomistic theories of nucleation, as well as the (dis)continuity of cluster size. Atomistic theories allow for high $\Delta\overline{\mu}$ levels to be modeled and have been validated against experimental studies ¹¹¹. But despite the existence of such theories, much remains unclear regarding the initial stages of electronucleation and the role of the atomic-scale structure of the electrode surface. In this regard, explicit atomistic simulations can play an important role in elucidating the initial processes and mechanisms (vide infra).

Figure 2: Schematics showing the differences between classical (blue) and atomistic (red) nucleation theories. Shown in (a) is the relationship between the free energy of nucleation, $\Delta\overline{G}$, and the cluster size of a cluster ($n$ atoms) for a constant $\Delta\overline{\mu}$. In classical/atomistic theories, $\Delta\overline{\mu}$ is small/large enough to result in $\Delta\overline{G}(n)$ being differentiable/non-differentiable. In both theories, $n=n_{c}$ (the supercritical cluster size) is the size at which $\Delta\overline{G}(n)$ has a global maximum. Shown in (b) is the relationship between supersaturation, $\Delta\overline{\mu}$, and $n_{c}$ according to both theories. In classical/atomistic nucleation theory, $\Delta\overline{\mu}$ is low/high enough such that $n_{c}$ is continuous/discontinuous. Also shown in (b) are visualizations of how both theories treat clusters: atomistic theories are applicable for small clusters, whereas classical theories can be used for large clusters which are treated as spheres.

3 Simulation of Electrochemical Reaction Conditions

In this section, we summarize the various methodologies employed to simulate electrochemical reaction conditions. This includes an accurate account of the electrode potential, the electrolyte, and the electrode surface.

3.1 The Electrode Potential

The electrode potential can greatly affect the reaction thermodynamics and kinetics during electrodeposition, and it thus becomes essential to accurately model the potential during atomistic simulations. Electrochemical experiments are typically performed under a constant potential and are referenced against well-defined reference electrodes. However, modeling a constant chemical potential within atomistic simulations is quite challenging, and various schemes have been developed to control the applied potential ^{134, 135, 48}. The current schemes to model the applied electrode potential can be roughly split into three categories ¹³⁶: classical forcefields ¹³⁷, finite-field methods ^{138, 139, 140, 141} which can be used alongside forcefields or density functional theory (DFT), and grand-canonical ensembles (GCEs) with an electronic structure method ^{142, 143} such as DFT ^{144, 145}. The first two methods require the full cell, while forcefields, in particular, do not treat electrons and can only treat electrostatics. Grand-canonical treatments, however, only require a half-cell to be modelled and can vary the number of electrons. The difference between full- and half-cell approaches is illustrated in Figure 3.

Classical forcefield and finite-field approaches typically describe the applied electrode potential in terms of the inner potential difference between two electrodes, and do not take into account the electronic structure ^{146, 147, 148, 138, 139, 149, 150, 151, 137, 140, 141, 136}. In both of these approaches, two electrode interfaces need to be simulated in order to build an electrode potential difference within the cell to enforce charge neutrality. However, as the properties of only the working electrode are of interest, the second electrode acts as a passive counter electrode. While the treatment of the additional electrode can be justified for forcefield potentials, it becomes computationally intractable for electronic structure methods, which computationally scale more severely than classical forcefields ¹⁵² with the number of atoms.

Forcefields come in a variety of flavors, including non-reactive and reactive, and often possess a trade-off between accuracy, via highly-detailed parameterizations, and transferability, ie. the number of materials to which they apply. Forcefields are widely used to simulate the energetics and dynamics of large systems by describing their nuclear interactions with a set of empirical equations rather than allowing for the electronic degrees of freedom that quantum mechanical methods explicitly specify. Recently, machine learning (ML) has emerged as a method to parameterize classical forcefields by supplying training data from quantum mechanical methods. Most ML-generated forcefields only account for local interactions though, which will not suffice for electrodeposition simulations, where long-range interactions have an important part to play, unless they are cleverly integrated alongside long-range interactions ^{153, 154, 155, 156, 157, 158, 159, 160, 161}.

In contrast, GCEs only require one electrode to be modeled and thus enable half-cells to be simulated. This is achieved by fixing the Fermi level (the energy level at which there is a 50% probability of it being occupied by an electron) of the electrode (which is equal to the chemical potential of electrons, $\tilde{\mu}_{e}$), while the chemical potential of the electrolyte, $\tilde{\mu}_{s}$, is dependent on the electrolyte solution and its concentration ³². In order to have a well-defined treatment of an electrochemical solid-liquid interface, the electronic structure method (typically DFT) needs to be part of a GCE with a fluctuating number of electrons and ions (at a given temperature), as depicted in Figure 3, rather than the more common canonical ensemble, where the number of particles is constant but the chemical potentials are allowed to fluctuate ³². Practically, this is most elegantly achieved by fixing the Fermi level and allowing the number of electrons to fluctuate during a calculation ³². While stable algorithms for this method have been developed ^{31, 162}, the fluctuating number of electrons can cause difficulties with convergence in calculations ³². An alternative methodology is to perform calculations at several points with a constant number of electrons and then interpolating to the desired Fermi level ^{163, 164, 165, 166, 167, 168, 169, 170}.

Figure 3: Comparison of (a) full cell and (b) half-cell simulation approaches. Shown in (b) is the interface between the solvent and the electrode in a grand-canonical ensemble, which only requires a half-cell to be simulated. The chemical potentials of the solvent/electrolyte and electrons are fixed to $\tilde{\mu}_{s}$ and $\tilde{\mu}_{e}$, respectively.

Changes in the Fermi level directly correspond to changes in $\tilde{\mu}_{e}$, which is obtained by changing the charge state of the electrode. This presents a problem for simulations as electrochemical systems are usually partially periodic, but such systems need to be charge neutral. Various methodologies have been proposed to address this, including the introduction of a homogeneous background charge, correction schemes ^{164, 171, 172}, joint DFT ^{173, 168, 31}, and modified Poisson-Boltzmann implicit solvation models ^{174, 175, 176} (vide infra).

An important quantity that is often used to analyze how changes in the electrode potential impact the electrochemical system is the electrostatic potential. The redistribution of electrons and ions in the system needed to maintain charge neutrality after changes in the electrode potential is reflected by changes to the electrostatic potential profile, particularly in the double layer region near the electrode. Figure 4 shows a typical electrostatic potential profile; at the cathode surface, the electrostatic potential typically exhibits oscillatory behavior due to the periodic arrangement of atoms and the resulting alternating regions of positive and negative charge density. In the bulk solution, the electrostatic potential tends to zero as the overall charge distribution is around neutral due to the lack of nearby charged surfaces and the cancellation of potentials from solvated ions. In the double layer, due to the presence of counterions close to the cathode, there is a sharp drop in the electrostatic potential often referred to as the Stern potential, which typically exhibits a large gradient near the cathode surface and decays towards zero as the potential extends towards the bulk solvent.

Figure 4: Schematic of an example electrostatic potential ($\phi$) profile (black) within the cathode, the double layer, and the bulk solvent.

Fermi level-fixed GC-DFT has been widely used to model electrochemical thermodynamics and kinetics. However, this is not suitable for outer-sphere reactions, semiconductors or two-electrode systems as in these systems, the Fermi level typically lies within the band gap rather than within the conduction or valence bands. To address the shortcomings of Fermi level-fixed GC-DFT, constant inner potential DFT has recently been proposed, which utilizes the local electrode inner potential as the thermodynamic parameter for the electrode potential, rather than the global Fermi level ¹³⁶. Both GC-DFT variants have been shown to provide identical results for metallic electrodes ¹³⁶, but differences can arise for semiconducting metal oxide–water interfaces ¹⁷⁷.

3.2 The Electrolyte

Before solvated cations are electrodeposited onto an electrode surface, they must diffuse through the solvent. Atomistic simulations must therefore ensure they appropriately describe the electrolyte, which includes the solvent itself and charged ions that, upon application of an electrochemical potential, form an electrochemical double-layer that modifies the electrostatic potential above the electrode surface ³².

Figure 5: Schematic representations of explicit, hybrid, and implicit solvation environments for a metal nanocluster adsorbed onto an electrode surface within a simulation box. The electrode shown is an atomistic model of a copper (111) surface. Analytical representations of the solvent are shown in blue, while hydrogen, oxygen, adsorbed metal, and copper atoms are colored white, red, gold, and brown, respectively.

Figure 5 shows the various ways the solvent environment can be described for an electrode surface. Within atomistic simulations, the solvent environment can be described either by using explicit molecules, an implicit (analytical) model, or a hybrid of the two. In the former case, the explicit modeling of solvent molecules is a more ‘realistic’ and physically meaningful description of the system. However, explicit solvation models are typically more computationally expensive than implicit models and in order to model a physically meaningful system, many explicit solvent molecules need to be included within the model, and these can contribute to over 90% of the atoms within a modeled system ¹⁷⁸. Using quantum mechanical methods to model solute-solute, solvent-solvent, and solute-solvent interactions can therefore quickly become computationally intractable with explicit solvent models comprising many hundreds of atoms. In order to reduce computational costs, empirical molecular mechanical forcefields have become a popular choice to treat interactions within the system. However, care must be taken during the parameterization of such forcefields in order to not sacrifice the accuracy that comes with ab initio approaches for the sake of computational tractability. In this regard, machine-learned interatomic potentials (MLIPs) offer a lot of promise ^{179, 180}, assuming the existence of appropriate training data (vide infra).

In most implicit solvent models, the solvent is treated as an electrostatic continuum with predefined dielectric and interfacial properties ^{181, 182, 183, 178}, and the solute is placed inside a cavity within the continuum. In some cases, a dependency on distance from the solute can be included, and it is also possible to introduce a dependence on the rate of a particular process, whereby the response of the solvent varies for fast and slow processes. When using implicit solvent models, a number of choices need to be made, including the shape and size of the cavity that the solute will be placed within. Examples of simple cavity shapes are spherical and ellipsoidal, while more complex ones that can be generated algorithmically such as van der Waals surfaces (based on the van der Waals radii ¹⁸⁴ of atoms), Lee-Richards molecular surfaces ¹⁸⁵, Connolly surfaces ^{186, 187}, and surfaces based on charge density isosurfaces ¹⁸⁸ are also possible cavities.

One of the advantages of implicit solvent models is that the number of interacting particles and the number of degrees of freedom within a system are significantly reduced. Implicit solvation models are thus typically computationally cheaper than explicit models and are therefore a good practical choice for computationally demanding studies ^{178, 189, 190}. The reduced computational cost also means quantum mechanical methods become more tractable and can be used to treat the solute more accurately than typical molecular mechanical methods. However, the lack of an explicit atomistic description of the solvent can result in numerous interactions, such as hydrogen bonds (both with the solvent and within the solute), being neglected, overstabilized salt bridges, incorrect ion distribution ¹⁹¹, and unphysical sampling ^{178, 192, 193}. Furthermore, the implicit description of the solvent also introduces an artificial boundary between the solute and solvent, and this interface needs to be treated carefully.

Typically, the formation of the electrochemical double-layer is modeled using modified Poisson-Boltzmann implicit solvation models, which present the simplest level of fixed-potential GC-DFT for solid-liquid interfaces. Such models come in many varieties, including: linear, non-linear, and ion-size effects. With modified Poisson-Boltzmann solvers, charge neutrality can be maintained via the ionic distribution in the double-layer, though this constraint is not automatically fulfilled in non-linear models ¹⁹⁴ due to an over-simplification of equations and the presence of the cavity exclusion function ³². This means that when combining modified Poisson-Boltzmann solvers with a quantum mechanical method for charged periodic systems, charge neutrality must be enforced either using Lagrange multipliers ¹⁹⁴, a solvated jellium with a constant-charge background tempered by an implicit solvent ¹⁹⁵, or by treating these systems using ‘metallic’ boundary conditions ¹⁶³. With a solvated jellium model, a constant Fermi level is iteratively obtained and only selectively altered with respect to the eigen-levels of the adsorbed solvent ¹⁹⁵.

Finally, hybrid solvation models that seek to reap the benefits of both explicit and implicit solvation models can also be used to effectively model the electrolyte. In terms of the modeling, hybrid models can be facilitated using embedded cluster models (vida infra) and methodologies such as hybrid quantum mechanics/molecular mechanics (QM/MM). Here, the system can be partitioned into three regions, where the central region is treated using QM and comprises the solute and some explicit solvent molecules. This central part is embedded within a second layer and contains more explicit solvent molecules but is treated using MM. Finally, both aforementioned regions are embedded within a third layer, where the solvent is described using implicit models and represents the bulk solvent. This hybrid approach allows the local region of interest to be modeled with the accuracy of QM without subjecting the entire system to the typically higher computational costs that come with QM. However, care must be taken to ensure the various interfaces between the regions are modeled appropriately, and some studies have also observed a dependency on the number of added explicit solvent molecules ¹⁹⁶.

3.3 The Electrode Surface

3.3.1 The Electronic Structure Method

In order to simulate metal electrodeposition, atomistic models need to include an appropriate description of the electrode surface onto which metals will be deposited. The extended surface model and the choice of quantum mechanical method thus need to correctly account for the rich diversity of interactions that are present at adsorbate-electrode interfaces. Such interactions include hybrid organic-inorganic, long-range van der Waals, and long-range electrostatic interactions of charged species.

Kohn-Sham DFT ^{144, 145} is one of the most commonly used electronic structure methods to describe extended surfaces ¹⁹⁷ and materials ¹⁹⁸. Within DFT, increasingly accurate density-functional approximations are being developed ¹⁹⁹ that can represent the energetics and electronic structure of complex materials. To describe electrified solid-liquid interfaces efficiently, a pragmatic selection of well-tested density-functional approximations that balance computational efficiency and predictive accuracy is required ^{198, 197}. For example, a posteriori long-range dispersion corrections, such as the Grimme ^{200, 201, 202, 203, 204} and Tkatchenko-Scheffler ^{205, 206, 207, 208, 209, 210, 211, 212} families, to generalized gradient approximations have been shown to provide a reliably accurate representation of adsorption structures and energetics. However, considering the existing limitations in the quantitative experimental characterization of electrochemical systems in general, and the kinetics and dynamics of electrodeposition in particular, the accuracy of existing density-functional approximations and dispersion correction schemes is currently not the limiting factor in atomistic electrodeposition simulations. Nevertheless, the choice of quantum mechanical method used remains key for the accurate modeling of processes at electrode surfaces.

However, the intrinsic computational scalability of DFT and even beyond-DFT, such as wavefunction methods and many-body perturbation theory, provides a challenge for systems comprising more than a few hundred atoms, which can often be the case for complex electrode surfaces. DFT typically has a formal $\mathcal{O}(N^{3})$ scaling, while higher-level ab initio wavefunction methods can formally scale beyond $\mathcal{O}(N^{6})$ (where $N$ is the number of basis functions). Semi-empirical tight-binding methods, such as density functional tight-binding ²¹³, have grown in popularity as they comprise a good compromise between computational cost and accuracy and have been successfully employed to study the electrodeposition of metals ²¹⁴. However, tight-binding parameterizations are typically developed for a particular subset of elements for a specific purpose, which means they possess low transferability. Furthermore, only a few reliable parameterizations therefore currently exist for metal-organic interfaces ^{215, 216, 217, 214}, which can be problematic for describing metal electrodeposition on carbonaceous electrodes, such as graphite/graphene and boron-doped diamond, or even for multimetallic systems. For these reasons, DFT is still typically more popular than tight-binding methods. However, newer semi-empirical methods, such as xTB, have sought to correct for this drawback by employing a global and element-specific parameterization ^{218, 219}, rather than the element pair-specific parameters employed within density functional tight-binding.

Machine-learned interatomic potentials (MLIPs) are another methodology that can offer high computational efficiency and perform calculations at an ab initio-level of accuracy, if the appropriate training data is supplied. Currently, most MLIPs do not sufficiently capture long-range electrostatics and thus will likely not suffice for electrochemical simulations, which require long-range interactions to also be described, though some studies have successfully captured long-range effects within short-ranged MLIPs^{153, 154, 155, 156, 157, 158, 159, 160, 161}. It should be noted though, that while there exists no universally best MLIP for every problem (as any pair of optimization algorithms should be equivalent when their performances are averaged over all possible problems ²²⁰), for practical purposes, care must be taken to ensure the atomic environment is appropriately represented by the machine learning architecture.

Quantum mechanical methods can be used to elucidate many key interactions and processes in electrodeposition. For example, DFT (or an alternative method) can be used to conduct single-point energy calculations or geometry optimizations, the latter resulting in a structure that is an energetic minimum on the potential energy surface. It should be noted that the relaxed structure could be a local minimum rather than a global minimum; various algorithmic processes exist that can be used in conjunction with quantum mechanical methods to identify global and local minima ^{221, 222, 223}. Quantum mechanical methods such as DFT and even semi-empirical methods ²²⁴ typically also permit the calculation of the electronic density of states, which can be used to elucidate various interactions. For example, if a copper cation is a few angstroms away from the cathode surface, its projected electronic density of states will be centered at the Fermi level and will contain one electron in the 4$s$ orbital; if it has been deposited, it will have a much broader density of states ⁸⁹. Such interactions can also be extracted from some machine learning frameworks, which have been developed to predict the electronic properties; such approaches retain full access to the electronic structure at forcefield-level efficiency while capturing QM in an analytically differentiable representation ²²⁵.

3.3.2 Structural Model of the Electrode Surface

For homogeneous electrode surfaces which are atomically flat, periodic boundary conditions can be used as an effective method to model extended surfaces and interfaces ¹⁹⁷. Here, the surface can be represented as a repeated slab, which is trivially defined using a unit cell that is infinitely repeated in three directions. For a surface, it is important to ensure the unit cell possesses a large enough vacuum in the $z-$direction to avoid periodic images from interacting with each other ¹⁹⁷. The unit cell should also be large enough to ensure calculations do not suffer from finite size effects ¹⁹⁷; this can be assessed by performing convergence tests on different unit cell sizes. The reader is directed to Hofmann et al. ¹⁹⁷ for a more detailed review on the repeated slab approach and practical considerations.

However, the surfaces of many electrode materials, particularly semiconducting electrodes, often possess a high degree of heterogeneity and can include structural defects such as point defects, dopants, dislocations, and lower coverages. In fact, defect sites within electrode surfaces have been shown to anchor and stabilize metal nanostructures ^{108, 226}. This stronger interaction between the nanostructure and the defect site can increase the reactivity exhibited by the supported nanostructure ^{227, 228, 229, 230, 231} in catalytic reactions. Atomistic simulations thus need to be able to account for the heterogeneity of such electrode surfaces. Modeling isolated defects at extended surfaces can be challenging with periodic boundary conditions, as isolated defects within the unit cell acquire a periodicity which can be unphysical. Furthermore, some doped electrode materials typically possess a relatively low dopant concentration e.g. the dopant concentration within boron-doped diamond electrodes is typically around 0.1% ²³². To model such an electrode surface using periodic boundary conditions would necessitate a very large unit cell, which can be computationally intractable with higher-level theories.

The challenges associated with the periodic representation of defects can be overcome by creating truncated cluster models. However, this removes the long-range properties of any bulk material and such calculations can be plagued by spurious finite size effects. Embedded cluster calculations are a viable alternative to periodic slab calculations as they acknowledge the intrinsic locality of surface defect chemistry and permit isolated defects to be modeled whilst breaking translational periodicity. Embedded cluster models have been treated using a variety of approaches, such as QM/QM ²³³, MM/MM and QM/MM, though care must be taken to ensure the embedded region is truncated appropriately, and the interface between the embedded and embedding regions is treated correctly. Such embedded cluster models are also typically computationally cheaper, which makes the application of higher-level theories more straightforward.

Figure 6: Schematic of how electrode surfaces can be represented within atomistic simulations. For well-defined pristine areas, the electrode can be modeled as a repeated slab via a unit cell to define periodic boundary conditions. For local areas of interest, such as a defect (shown as a yellow X), the surface can be modeled using an embedded cluster approach and partitioned into different regions, shown in blue and green. The electrode surface shown is an atomistic model of an oxygen-terminated polycrystalline boron-doped diamond electrode ^{26, 234, 108}.

4 Simulation of Electrodeposition Processes

Here, we review current atomistic simulation methods and how they capture the elementary steps of electrodeposition described in Figure 1.

4.1 Electrodiffusion of Solvated Ions

The most common atomistic simulation method to model diffusion is molecular dynamics (MD). MD simulations make use of interaction potentials, and the dynamics of particles (ions, atoms or molecules) are described by numerically integrating Newton’s equations of motion, with the forces being computed as the derivatives of the interaction potential. MD simulations that run for long enough should be able to describe the dynamical and structural properties at finite temperatures, as well as the thermodynamic equilibrium. Furthermore, each atomistic trajectory resulting from an MD run allows for the complex mechanisms that drive (electro)chemical processes to be identified. MD simulations, however, cannot easily reach timescales beyond nanoseconds and, without any form of enhanced sampling, cannot describe rare events such as electrodeposition, which can occur over a timescale of seconds. Furthermore, by treating atoms classically, MD alone can only calculate the positions, accelerations and forces of atoms. As protons and electrons cannot be explicitly modeled, key electrochemical processes such as proton conduction cannot be simulated by MD alone. Ab initio MD has become popular as it allows for electronic interactions to be included within simulations. Unlike classical MD, ab initio MD uses a first-principles method, typically DFT, to calculate forces by solving the time-dependent Schrödinger equation. While ab initio MD can be more accurate than classical MD, the additional complexity introduced by taking into account the QM interactions results in a much higher computational cost.

An alternative to MD are Monte Carlo simulations, which randomly sample potential configurations of a system rather than solve Newton’s equations of motions. Monte Carlo simulations can therefore be used to gain insights into the probabilistic nature of metal cation motion, as well as other quantities such as the distribution of electrodeposited metal atoms. However, they cannot directly extract dynamical observables connected to time-correlated functions, such as diffusion rates. Kinetic Monte Carlo ²³⁵ simulations analyze the time evolution of a system based on the rates of different elementary processes. In fact, kinetic Monte Carlo simulations are particularly useful for studying rare events and processes that occur over large timescales, which are typically not computationally feasible with conventional MD or DFT, making them very applicable for the description of diffusion phenomena during electrodeposition ^{236, 237, 238, 239, 240}.

4.2 Electron Transfer Reactions

Commonly, model Hamiltonians are used to describe interactions between solvated species and the electrode surface and to represent the electronic structure of the system. This is achieved via the inclusion of terms that represent the electronic states of the system and by defining the energy levels of electrons within various species. Other terms that are included within the model Hamiltonian account for Coulombic repulsion (for electron-electron interactions during electron transfer) and hybridization terms (for the overlap between orbitals on the adsorbed metal species and the electrode surface) which can influence the probability of electron transfer occurring.

The Newns-Anderson model Hamiltonian ^{90, 91} is a very popular framework that has been used to describe electronic interactions during atomistic simulations of electrodeposition. This framework reduces the complexity of the system while capturing the essential components of electron transfer. However, there are several limitations of the Newns-Anderson model Hamiltonian, such as the absence of electron correlation and the assumption that the electronic coupling between the electrode and the electrolyte is constant. Solvent effects are also not accounted for, rendering the Newns-Anderson Hamiltonian inappropriate for describing electron transfer processes in highly polarized solvents, and the adiabatic approximation is assumed to be valid, where atomic positions are assumed to not instantaneously change with electronic states during electron transfer, which is inappropriate for describing non-adiabatic effects such as coupled electron-proton transfer or dynamically induced electron transfer.

Extensions to the Newns-Anderson model Hamiltonian have been developed to explicitly include key microscopic parameters in a single formulation ^{97, 241, 242, 243, 244}. For example, inspired by Marcus theory ⁸⁷, the Schmickler-Newns-Anderson Hamiltonian ^{242, 241} was developed to describe electrochemical electron transfer by incorporating the effect of the solvent. Recently, the Schmickler-Newns-Anderson Hamiltonian has been further modified to account for the electrostatic interaction between the electrode and the redox couple ^{245, 246, 95} and is expressed as the sum of four terms:

$$\mathcal{H}=\mathcal{H}_{\mathrm{el}}+\mathcal{H}_{\mathrm{sol}}+\mathcal{H}_{\mathrm{int}}+\mathcal{H}_{\phi}$$ (15)

where $\mathcal{H}_{\mathrm{el}}$ is the electronic contribution and for reactant and electrode surface orbitals ($a$ and $k$, respectively), can be expressed as:

$$\mathcal{H}_{\mathrm{el}}=\varepsilon_{a}(d)\hat{n}_{a}+\sum_{k}\varepsilon_{k}\hat{n}_{k}+\sum_{k}\big{(}V_{k}(d)c_{k}^{+}c_{a}+V_{k}^{*}(d)c_{a}^{+}c_{k}\big{)}$$ (16)

where $\varepsilon_{a}$ is the electronic energy of the redox couple, $d$ is the distance between the metal cation and the electrode, $\hat{n}$ is the operator for the occupation number of the redox orbital, and $V_{k}$ is the interaction parameter that characterizes the strength of electronic interactions ^{245, 95}. The last term in Equation (16) accounts for electron transfer between the metal cation and the electrode surface, with $c^{+}$ and $c$ denoting creation and annihilation operators, respectively ²⁴⁵. As the metal cation approaches the electrode surface, $\varepsilon_{a}$ will shift towards the Fermi level of the electrode ^{245, 244, 247}. This phenomenon is referred to as Fermi level pinning and occurs due to the stabilization of the cation by the electric field of the cathode and the interaction with the electronic states of the electrode.

The solvent contribution for classical nuclei, $\mathcal{H}_{\mathrm{sol}}$ in Equation (15), is given by:

$$\mathcal{H}_{\mathrm{sol}}=\dfrac{1}{2}\hbar\omega(p^{2}+q^{2})$$ (17)

where $\hbar$ is the Planck constant, $\omega$ is the solvent frequency, $p$ is the solvent momentum, and $q$ is the solvent reorganization coordinate. The interaction energy, $\mathcal{H}_{\mathrm{int}}$ in Equation (15), between the solvent and the reactant linearly depends on the solvent coordinate and the coupling strength, $g$:

$$\mathcal{H}_{\mathrm{int}}=(z-n)\hbar\omega gq$$ (18)

Finally, the electrostatic interaction between the redox couple and the electrode surface, $\mathcal{H}_{\phi}$ in Equation (15), can be expressed as:

$$\mathcal{H}_{\phi}=(z-n)\phi(d)$$ (19)

where $z\phi(d)$ represents a repulsive Coulombic interaction between the redox couple center and the electrode and thus possesses an opposite sign to the electrostatic potential obtained from electronic structure calculations, where the electrostatic potential is calculated from the perspective of electrons ⁹⁵. The right-hand side of Equation (19) reflects the change in Coulombic interaction due to reduction and as $\phi(d)$ increases, the stabilization due to reduction increases.

The extended Hamiltonian in Equation (15) contains information about both pre-exponential factor and the free energy barrier in Equation (8), and can thus be used to calculate both adiabatic and non-adiabatic OS-ET rates ²⁴¹. As atomistic simulations can be used to obtain all the quantities that enter Equation (15), it becomes much easier to gain a physical and chemical understanding of OS-ET kinetics.

4.3 Electronucleation

Simulating the amalgamation of surface-adsorbed metal atoms to form metal nanostructures is computationally expensive due to the large number of degrees of freedom and the sheer number of nucleation pathways possible. Some studies investigate the atom-by-atom growth of nanostructures via single-atom absorption, but as mentioned above, this does not hold in reality as entire clusters can merge or fragment. However, analyzing how adsorption energies, cohesive energies, and interatomic distances change as a function of nanostructure size can still provide many valuable insights ^{228, 248}.

However, for a given atom count, metal nanostructures can exhibit numerous metastable geometries. Experimentally, the determination of the geometric ground state of surface-adsorbed nanoclusters is very difficult, and simulations are much more promising for providing insights. However, the reliable identification and optimization of surface-adsorbed metal nanostructures is particularly challenging due to the structural complexity and the large number of degrees of freedom, such as the number of possible metastable geometries, adsorption site, and surface coverage ¹⁹⁷. Algorithmic approaches such as basin- ^{221, 222} and minima-hopping ²²³ can be used to identify global and local minima. Recently, nested sampling has been extended to calculate coverage-temperature adsorbate phase diagrams by incorporating all relevant configurational contributions to the free energy ²⁴⁹. However, an additional complexity is introduced when attempting to realize the configuration under certain electrochemical conditions.

Isaev et al. ¹³² sought to study the time dependence of the radius of a single hemispherical nanocluster electrodeposited onto an electrode. In particular, they proposed models of formation and diffusion-controlled growth of a new-phase nanocluster for potentiostatic electrodeposition, cyclic voltammetry, and galvanostatic electrodeposition. The nanocluster growth rate during galvanostatic deposition was found to be much lower than under potentiostatic condition due to the drop in the overpotential that occurs after formation of the nanocluster. However, at larger currents, multiple nucleation will occur and increasing the concentration of the depositing cations would cause the number of clusters to decrease and each cluster size to increase.

5 Capabilities of Software Packages

In order to simulate electrodeposition, users should seek to use established software packages that can account for key electrochemical variables such as solvation effects. The choice of software should be influenced by the task at hand. A good software package should be well-documented and easy to use. Furthermore, softwares should possess efficient computational scalability across multiple cores as atomistic simulations of complex processes such as electrodeposition can be computationally expensive. To this end, several software packages have impressive parallelism of calculations over central processing units (CPUs) and have started to support acceleration of calculations with graphical processing units (GPUs) ^{250, 251, 252, 253, 254}, which can be useful for improving the computational tractability of large systems and complex simulations. Many of the software packages mentioned below are also currently undergoing development to enable integration with machine learning models. This section seeks to provide a non-exhaustive list of software packages for atomistically simulating the various aspects of metal electrodeposition. For a more beginners’ introduction to computational modeling techniques, see Ismail et al. ²⁵⁵.

While there are numerous software packages that contain DFT implementations; David Sherrill et al. ²⁵⁶ is a useful resource that lists the several electronic structure packages that exist. Solid-state codes for condensed-phase systems are generally more applicable than molecular codes. Examples of such codes include (in alphabetical order): the wavelet code BigDFT ²⁵⁷, the plane-wave pseudopotential code CASTEP ²⁵⁸, CP2K ²⁵⁹, the all-electron numeric atomic orbital code FHI-aims ²⁶⁰, the Python-based projector augmented-wave code GPAW ^{261, 262, 254}, NWChem ^{263, 264}, the linear-scaling code ONETEP ²⁶⁵, the Python-based PySCF ^{266, 267}, Quantum ESPRESSO ^{268, 269}, TURBOMOLE ²⁷⁰, and VASP ²⁷¹. There are many factors to consider when selecting which DFT software package to use. One important choice is the type of basis set, used to represent the Kohn-Sham orbitals. Plane wave approaches rely on delocalized basis functions but require a special treatment for the core electrons, which can be treated using pseudopotentials or projector-augmented waves ¹⁹⁷. In contrast, atom-centered basis functions are centered on the nuclei and more natively describe core electrons ¹⁹⁷. However, it should be noted that simulations employing plane wave and atom-centered basis functions should give the same results in the converged limit ²⁷². If tight-binding methods are desired in preference to DFT, then software packages such as DFTB+ ²⁷³ and xTB ^{218, 219} can be used to perform atomistic simulations, though care must be taken when selecting parameter sets.

Most DFT software packages include a wide range of density-functional approximations that can be used for calculations, with many of them having interfaces to the LIBXC ¹⁹⁹ library of functionals. In addition, most softwares also include dispersion corrections, with many of them interfacing with standalone programs such as DFT-D4 ²⁷⁴ or libraries such as libMBD ²¹² or libvdwxc ²⁷⁵. Many packages now include implementations of implicit solvation models that can be combined alongside QM methods such as DFT, which is key for simulating electrodeposition.

Quantum mechanical software packages that allow for a grand-canonical treatment of the system include CP2K ²⁵⁹, GPAW ^{262, 261, 254}, Quantum ESPRESSO ^{268, 269}, and VASP ²⁷¹, among others. DFT calculators supported within the popular Python-based ASE ²⁷⁶ can also be turned into GC-DFT calculators by changing only around 20 lines of code, as shown by Melander ¹⁴³. ASE itself provides interfaces to many software packages that can be used as calculators for many purposes such as energy calculations, geometry optimizations, and other complex algorithmic processes such as nudged elastic band for transition path calculations and global structure searches ²⁷⁶. Recently, the Python-based matscipy package was published, which is interoperable with ASE and includes a matscipy.electrochemistry module to permit electrochemical simulations that describe the motion and spatial distribution of charged species within an external electric field ²⁷⁷. The matscipy.electrochemistry module includes tools that sample discrete coordinate sets from continuum fields and apply steric corrections to avoid any overlap between species ²⁷⁷. It can therefore be used to generate molecular coordinates that resemble ionic distributions within an electrochemical double layer which can be used as the initial configuration for atomistic simulations.

Many software packages also have the ability to interface with others, which can be particularly useful for embedded cluster calculations, for example. Softwares such as the aforementioned Python-based ASE ²⁷⁶ and ChemShell ²⁷⁸ implement a modular approach to embedded cluster models, where a variety of QM and MM codes can be called upon to evaluate the energies and gradients of the various partitioned regions. ChemShell, scriptable using either the Tcl or Python programming languages ²⁷⁸, can call upon numerous QM and MM codes. The recent Atomic Simulation Interface (ASI), which has a native C-style API, also defines a set of functions that support both classical and ab initio MD, and hybrid QM/MM calculations ²⁷⁹. ASI provides an efficient way to import and export large arrays from electronic structure codes, such as the Hamiltonian, overlap, and density matrices, that are typically monolithic and can be used to develop universal and interoperable codes without sacrificing efficiency for portability ²⁷⁹.

Large-scale molecular mechanical modeling of materials are most efficiently achieved with established softwares such as DL_POLY ²⁸⁰, the Python-based HOOMD-blue ²⁸¹, and LAMMPS ²⁸². LAMMPS in particular is well-equipped to conduct electrochemical simulations, comprising implementations of packages such as ELECTRODE ²⁸³, which enable computationally efficient simulations of electrochemical reactions, electrode potentials, and ion transport within an electrochemical cell. MD simulations have become increasingly popular as they are computationally well-suited to capitalize upon the hardware speedups bestowed by Moore’s law and the resulting large-scale parallelism over both CPUs and GPUs in recent years ²⁸². Efficient MD codes can therefore be used to simulate systems from atomic to mesoscale length scales at picosecond to microsecond timescales ²⁸². Within MD, the computational cost for models with short-range interactions has formal linear scaling i.e. $\mathcal{O}(N)$, and long-range Coulombic interactions have formal $\mathcal{O}\big{(}N\log(N)\big{)}$ scaling ²⁸².

6 Applications and Challenges

New experimental techniques are emerging that can provide valuable data to support and validate atomistic simulations. For example, electron microscopy methods have been a popular choice to study metal electrodeposition due to the high resolution they offer ^{284, 285, 286, 287, 26}. Scanning transmission electron microscopy has been shown to be capable of dynamically visualizing the early stages of electronucleation for metals such as gold, with structural resolution on the atomic scale and time resolution defined by the sequential analysis of short electrodeposition runs (several milliseconds of deposition per run) ²⁶. Transmission electron microscopy has also been used to study the electrodeposition of metals such as lithium, silver, nickel and platinum at submicroscopic resolution ^{284, 285, 288, 286, 287, 289}. Several studies also report the use of scanning electron microscopy to investigate metal electrodeposition ^{290, 286, 291}. While liquid-cell transmission electron microscopy has made lots of progress in monitoring dynamic electrochemical systems ^{292, 293}, it has limited resolution due to factors such as electron-beam-induced gas bubble formation and electron scattering in the liquid ²⁶; in contrast, ex situ aberration-corrected scanning transmission electron microscopy is not only capable of resolving single atoms but can be used to quantify the number of atoms within a particle ^{294, 26}.

In contrast, scanning probe methods such as scanning tunnelling microscopy and atomic force microscopy generate images of surfaces using a physical probe that scans the sample ²⁹⁵. Both atomic force microscopy ²⁹⁶ and scanning tunnelling microscopy ^{297, 298, 299, 300, 301, 302, 303} are capable of analyzing the influence of current on the electrodeposited structure at submicroscopic resolutions. Scanning electrochemical cell microscopy has also been used to study the initial electronucleation stages and mobility of metals such as platinum ³⁰⁴, copper ³⁰⁵, iron ³⁰⁶ and silver ²⁵, and has gained much attention ^{307, 308, 309} due to its ability to routinely operate at submicroscopic scales ^{305, 310, 311}. Other microscopy techniques such as surface plasmon resonance microscopy ³¹² and dark-field scattering microscopy ^{313, 314} also exist and have been used to investigate metal electrodeposition. Surface plasmon resonance microscopy is a non-intrusive optical technique that is highly sensitive to nanoscopic objects, while dark-field scattering microscopy has a relatively simple experimental setup that can directly probe the plasmonic properties of individual structures and allow for correlation with electron microscopy ^{313, 314}. Their performance, however, is restricted by the very small field of view which increases the difficulty in acquiring quantitative data on electronucleation ^{312, 313}. Wide-field surface plasmon resonance microscopy, however, removes this constraint and allows for the growth of hundreds of nuclei to be tracked simultaneously at a reasonable time resolution (${\sim}$$1\text{\,}\mathrm{s}$) ³¹².

Advancements in the theory and atomistic simulation of metal electrodeposition have helped to elucidate many phenomena. However, several open challenges remain that hinder the full realization of accurate models. One primary challenge is the precise modeling of the electrochemical interface. The dynamic nature of this interface requires sophisticated simulation methods capable of capturing these nuances. Additionally, incorporating quantum mechanical effects into large-scale simulations without incurring prohibitive computational costs is an ongoing challenge.

To model electrodeposition, the electrosorption, electron transfer, on-surface migration and nucleation processes need to be captured. So far, no simulation approach exists that captures all these aspects atomistically. Bridging the gap between the femtosecond timescales of electron transfer and atomic motions, and the multisecond timescales of electrodeposition processes demands the development of advanced algorithmic approaches. Understanding electronucleation and growth mechanisms at the atomic level, particularly the role of surface defects and heterogeneities, pose additional challenges. The effects of solvent molecules and ions, including solvation dynamics and ion pairing, further complicate these simulations. Continuum and atomistic approaches need to be combined to create multiscale approaches that can be to used to tackle such challenges. It is evident that MLIPs and other ML surrogate models coupled with electronic structure theory ³¹⁵ will have to play a crucial role in developing a unified multi-scale modelling description of electrodeposition.

Both theory and experiment will be needed to answer the open questions in the field. In terms of electronucleation, this includes consideration of non-classical growth pathways such as surface migration, aggregation and coalescence of small nanoclusters ^{25, 118}; formation of metastable clusters into crystalline nanoparticles ²⁶; and nucleation and dissolution events that occur before stable nuclei form ³¹⁶. Other open questions ^{317, 111} include why the measured number of nuclei is higher than the calculated number of active sites and why single atoms are so stable. With the emergence of improved simulation methods and experiments with high spatial and temporal resolution, there is promise for these open questions to be answered in the near future.

This work, in part, builds on S.C.’s doctoral thesis ²⁴⁸. This work was funded by the EPSRC Centre for Doctoral Training in Diamond Science and Technology [EP/L015315/1], the Research Development Fund of the University of Warwick, Wellcome Leap as part of the Quantum for Bio Program, the UKRI Future Leaders Fellowship programme [MR/S016023/1 and MR/X023109/1], and a UKRI Frontier research grant [EP/X014088/1].