Linear stability analysis of transient electrodeposition in charged porous media: suppression of dendritic growth by surface conduction

TL;DR

This paper analyzes how charged porous media influence electrodeposition stability, showing that negative surface charges can suppress dendritic growth by enabling surface conduction, with results validated against experimental data.

Contribution

It introduces a linear stability framework for transient electrodeposition in charged porous media, highlighting the stabilizing effect of negative surface charges and validating predictions with experiments.

Findings

Negative charges enable sustained overlimiting current via surface conduction.

Negative charges significantly reduce surface instabilities and dendritic growth.

Theoretical predictions align well with experimental data on copper electrodeposition.

Abstract

We study the linear stability of transient electrodeposition in a charged random porous medium, whose pore surface charges can be of any sign, flanked by a pair of planar metal electrodes. Discretization of the linear stability problem results in a generalized eigenvalue problem for the dispersion relation that is solved numerically, which agrees well with the analytical approximation obtained from a boundary layer analysis valid at high wavenumbers. Under galvanostatic conditions in which an overlimiting current is applied, in the classical case of zero surface charges, the electric field at the cathode diverges at Sand's time due to electrolyte depletion. The same phenomenon happens for positive charges but earlier than Sand's time. However, negative charges allow the system to sustain an overlimiting current via surface conduction past Sand's time, keeping the electric field bounded. Therefore, at Sand's time, negative charges greatly reduce surface instabilities and suppress dendritic growth, while zero and positive charges magnify them. We compare theoretical predictions for overall surface stabilization with published experimental data for copper electrodeposition in cellulose nitrate membranes and demonstrate good agreement between theory and experiment. We also apply the stability analysis to how crystal grain size varies with duty cycle during pulse electroplating.

Figures

Tables

Table 1. Table 1: Scales used for nondimensionalization.

Variables and parameters	Scale
$x$, $y$, $z$, $L_y$, $L_z$, $r$, $r_{\text{m}}$, $h$, $\lambda$, $\xi$	$L_x$
$t$	$\frac{L_x^2}{D_{\text{amb0}}}$ (diffusion time)
$c_{\pm }$	${\nu }_{\pm }{c}_0$
$\varphi$, ${\varphi }_{\text{e}}$, $\Delta {\varphi }^{\text{eq}}$, $E^{\Theta }$, $\Delta \varphi$, $\eta$	$\frac{k_{\text{B}}T}{e}$ (thermal voltage)
$D_{\pm 0}$	$D_{\text{amb0}}$
${\rho }_{\text{s}}$	$z_{+}{\nu }_{+}e{c}_0=-z_{-}{\nu }_{-}e{c}_0$
$F_{\pm }$	$\frac{{ϵ}_{\text{p}}{D}_{\text{amb0}}{\nu }_{\pm }{c}_0}{L_x}$
$u$, $u_x$, $v_{\text{I}}$, $v_{\text{In}}$	$\frac{{ϵ}_{\text{p}}{D}_{\text{amb0}}}{L_x}$
$J$	$J_{\text{lim}}$
$I$	$I_{\text{lim}}$
$j_0$, $J_{\text{F}}$	$\frac{J_{\text{lim}}}{{ϵ}_{\text{p}}}$
$k_0$	$\frac{J_{\text{lim}}}{e{ϵ}_p}$
$\gamma$	$\frac{L_x{k}_{\text{B}}T}{\Omega }$
$\mathscr{H}$, $k_y$, $k_z$, $k$	$\frac{1}{L_x}$
$\omega$	$\frac{D_{\text{amb0}}}{L_x^2}$ (reciprocal diffusion time)

Table 2. Table 2: Dimensionless parameters that are used for results in Sections  IV.2 ,  IV.3 ,  IV.4 and  IV.6 for a typical electrolyte in a typical nanoporous medium.

Dimensionless parameter	Description	Value
${\nu }_{+}$	Number of cations formed from complete dissolution of 1 molecule of neutral salt	$1$
$\nu$	Number of anions formed from complete dissolution of 1 molecule of neutral salt	$1$
$z_{+}$	Cation charge number	$1$
$z$	Anion charge number	$-1$
$D_{+}$	Cation diffusivity	$1$
$D$	Anion diffusivity	$1$
$n$	Number of electrons transferred in charge transfer reaction	$1$
$\alpha$	Charge transfer coefficient	$0.5$
$\text{Ca}=\gamma$	Capillary number, ratio of capillary constant to inter-electrode distance (Equation 22)	$8.74\times {10}^{-5}$
${\beta }_{\text{m}}$	Ratio of atomic volume of solid metal to reciprocal electrolyte concentration	$1.30\times {10}^{-4}$
${\beta }_{\text{D}}$	Weighted ratio of cation and anion diffusivities	$0.25$
${\beta }_{\text{v}}$	Ratio of ${\beta }_{\text{m}}$ to ${\beta }_{\text{D}}$	$5.20\times {10}^{-4}$
${\xi }_{+}$	Ratio of electrolyte concentration to standard cation concentration	$0.01$
$E^{\Theta }$	Standard electrode potential	$0$
$L_y$	Aspect ratio in $y$ direction	$100$
$L_z$	Aspect ratio in $z$ direction	$100$
${\rho }_{\text{s}}$	Ratio of background charge density to electrolyte concentration	$-1,-0.75,-0.5,-0.25$, $-0.05,0,0.05$
Da	Damköhler number, ratio of reaction rate to diffusion rate (Equation 14)	$0.1,1,10$

Equations

$${ϵ}_{\text{p}}\frac{\partial c_{\pm }}{\partial t}+\nabla \cdot F_{\pm }$$

$$\nabla \cdot J$$

$${\rho }_{\text{s}}=\frac{{\sigma }_{\text{s}}}{h_{\text{p}}}$$

$$J_{\textnormal{lim}}=\frac{2\mathopen{}\mathclose{{}\left(z_{+}-z_{-}}\right)e\epsilon_{\textnormal{p}}D_{+0}\nu_{-}c_{0}}{L_{x}}$$

$${\tilde{t}}_{\text{s}}=\frac{\pi }{16{\tilde{J}}_{\text{a}}^2},{\tilde{J}}_{\text{a}}>1.$$

$${\mu }_{\text{O}}$$

$${\mu }_{\text{e}}$$

$${\mu }_{\text{R}}$$

$$\Delta {\varphi }^{\text{eq}}=\frac{k_{\text{B}}T}{ne}\ln \frac{a_{\text{O}}{a}_{\text{e}}^n}{a_{\text{R}}}+E^{\Theta }-\frac{2\Omega \gamma \mathcal{H}}{ne},E^{\Theta }=\frac{{\mu }_{\text{O}}^{\Theta }+n{\mu }_{\text{e}}^{\Theta }-{\mu }_{\text{R}}^{\Theta }}{ne},$$

$$J_{\textnormal{F}}=nek_{0}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{\mu_{\ddagger}^{\textnormal{r,ex}}-\mu_{\textnormal{O}}-n\mu_{\textnormal{e}}}{k_{\textnormal{B}}T}}\right)-\exp\mathopen{}\mathclose{{}\left(-\frac{\mu_{\ddagger}^{\textnormal{r,ex}}-\mu_{\textnormal{R}}}{k_{\textnormal{B}}T}}\right)}\right]$$

$$\mu_{\ddagger}^{\textnormal{r,ex}}=k_{\textnormal{B}}T\ln\gamma_{\ddagger}^{\textnormal{r}}+\mathopen{}\mathclose{{}\left(1-\alpha}\right)\mathopen{}\mathclose{{}\left(z_{\textnormal{O}}e\phi-ne\phi_{\textnormal{e}}+\mu_{\textnormal{O}}^{\Theta}+n\mu_{\textnormal{e}}^{\Theta}}\right)+\alpha\mathopen{}\mathclose{{}\left(z_{\textnormal{R}}e\phi+\mu_{\textnormal{R}}^{\Theta}+2\Omega\gamma\mathscr{H}}\right).$$

$$J_{\textnormal{F}}=j_{0}\mathopen{}\mathclose{{}\left\{\exp\mathopen{}\mathclose{{}\left(-\frac{\alpha ne\eta}{k_{\textnormal{B}}T}}\right)-\exp\mathopen{}\mathclose{{}\left[\frac{\mathopen{}\mathclose{{}\left(1-\alpha}\right)ne\eta}{k_{\textnormal{B}}T}}\right]}\right\},\quad j_{0}=\frac{k_{0}ne\mathopen{}\mathclose{{}\left(a_{\textnormal{O}}a_{\textnormal{e}}^{n}}\right)^{1-\alpha}a_{\textnormal{R}}^{\alpha}}{\gamma_{\ddagger}^{\textnormal{r}}},$$

$$\Delta {\varphi }^{\text{eq}}=\frac{k_{\text{B}}T}{ne}\ln {\hat{c}}_{+}+E^{\Theta }-\frac{2\Omega \gamma \mathcal{H}}{ne},j_0=k_{0}ne{\hat{c}}_{+}^{1-\alpha }.$$

$$\text{Da}=\frac{e{ϵ}_{\text{p}}{k}_0}{J_{\text{lim}}}.$$

$$h\mathopen{}\mathclose{{}\left(\xi,t}\right)=h^{\mathopen{}\mathclose{{}\left(0}\right)}\mathopen{}\mathclose{{}\left(t}\right)+\epsilon\Re\mathopen{}\mathclose{{}\left[h^{\mathopen{}\mathclose{{}\left(1}\right)}\exp\mathopen{}\mathclose{{}\left(ik\cdot\xi+\omega t}\right)}\right]+\mathcal{O}\mathopen{}\mathclose{{}\left(\epsilon^{2}}\right)$$

$$\displaystyle c_{-}\mathopen{}\mathclose{{}\left(x,\xi,t}\right)$$

$$\displaystyle\phi\mathopen{}\mathclose{{}\left(x,\xi,t}\right)$$

$$\displaystyle\theta\mathopen{}\mathclose{{}\left(x=h}\right)$$

$$\displaystyle\frac{\partial\theta}{\partial x}\mathopen{}\mathclose{{}\left(x=h}\right)$$

$$\displaystyle\frac{\partial\theta}{\partial y}\mathopen{}\mathclose{{}\left(x=h}\right)$$

$$\displaystyle\nabla\theta\mathopen{}\mathclose{{}\left(x=h}\right)$$

$$\text{Ca}=\tilde{\gamma }=\frac{\Omega \gamma }{L_x{k}_{\text{B}}T},$$

$$\displaystyle c_{0,t}-D\mathopen{}\mathclose{{}\left[c_{0,xx}+z\mathopen{}\mathclose{{}\left(c_{0}\phi_{0,x}}\right)_{x}}\right]$$

$$\displaystyle\beta_{\textnormal{D}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(D-D_{+}}\right)c_{0,xx}+z_{+}D_{+}\rho_{\textnormal{s}}\phi_{0,xx}-\mathopen{}\mathclose{{}\left(z_{+}D_{+}-zD}\right)\mathopen{}\mathclose{{}\left(c_{0}\phi_{0,x}}\right)_{x}}\right]$$

$${\varphi }_{\text{e}}^{\text{a}}$$

$$\displaystyle-D\mathopen{}\mathclose{{}\left(-c_{0,x}-zc_{0}\phi_{0,x}}\right)$$

$$\hat{n}\cdot J_0$$

$$-{\dot{h}}_0$$

$$\displaystyle-D\mathopen{}\mathclose{{}\left(c_{0,x}+zc_{0}\phi_{0,x}}\right)$$

$$\hat{n}\cdot J_0$$

$${\dot{h}}_0$$

$$\displaystyle c_{0}\mathopen{}\mathclose{{}\left(t=0}\right)$$

$$J_{\text{a}}$$

$$\displaystyle=\mathopen{}\mathclose{{}\left.\beta_{\textnormal{D}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(D-D_{+}}\right)c_{0,x}+z_{+}D_{+}\rho_{\textnormal{s}}\phi_{0,x}-\mathopen{}\mathclose{{}\left(z_{+}D_{+}-zD}\right)c_{0}\phi_{0,x}}\right]}\right\rvert_{x=h_{0}^{\textnormal{a}}}.$$

$$\displaystyle D\mathopen{}\mathclose{{}\left\{c_{1}^{\prime\prime}-k^{2}c_{1}+z\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(c_{0}\phi_{1}^{\prime}+\phi_{0,x}c_{1}}\right)_{x}-k^{2}c_{0}\phi_{1}}\right]}\right\}=\omega c_{1},$$

$$\displaystyle\mathopen{}\mathclose{{}\left(D-D_{+}}\right)\mathopen{}\mathclose{{}\left(c_{1}^{\prime\prime}-k^{2}c_{1}}\right)+z_{+}D_{+}\rho_{\textnormal{s}}\mathopen{}\mathclose{{}\left(\phi_{1}^{\prime\prime}-k^{2}\phi_{1}}\right)-\mathopen{}\mathclose{{}\left(z_{+}D_{+}-zD}\right)\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(c_{0}\phi_{1}^{\prime}+\phi_{0,x}c_{1}}\right)_{x}-k^{2}c_{0}\phi_{1}}\right]=0,$$

$$\displaystyle c_{0,t}h_{1}^{\textnormal{a}}-D\mathopen{}\mathclose{{}\left[-c_{1}^{\prime}-z\mathopen{}\mathclose{{}\left(c_{0}\phi_{1}^{\prime}+\phi_{0,x}c_{1}}\right)}\right]=0,$$

$$\displaystyle\beta_{\textnormal{v}}j_{0,0}\mathopen{}\mathclose{{}\left(\hat{D}_{1}h_{1}^{\textnormal{a}}+\hat{D}_{2}c_{1}+\hat{D}_{3}\phi_{1}}\right)=\omega h_{1}^{\textnormal{a}},$$

$$\displaystyle\beta_{\textnormal{m}}\mathopen{}\mathclose{{}\left[-\mathopen{}\mathclose{{}\left(D-D_{+}}\right)c_{1}^{\prime}-z_{+}D_{+}\rho_{\textnormal{s}}\phi_{1}^{\prime}+\mathopen{}\mathclose{{}\left(z_{+}D_{+}-zD}\right)\mathopen{}\mathclose{{}\left(c_{0}\phi_{1}^{\prime}+\phi_{0,x}c_{1}}\right)}\right]=\omega h_{1}^{\textnormal{a}},$$

Full text

I Introduction

Linear stability analysis is routinely applied to nonlinear systems to study how the onset of instability is related to system parameters and to provide physical insights on the conditions and early dynamics of pattern formation Langer (1980); Kessler et al. (1988); Cross and Hohenberg (1993). Some examples in hydrodynamics include the Orr-Sommerfeld equation that predicts the dependence on Reynolds number of the transition from laminar flow to turbulent flow Orr (1907a, b); Sommerfeld (1908); Orszag (1971) and the electroconvective instability that causes the transition of a quasiequilibrium electric double layer to an nonequilibrium one that contains an additional extended space charge region Zaltzman and Rubinstein (2007). Here, we focus on morphological stability analysis in which linear stability analysis is used to analyze morphological instabilities of interfaces formed between different phases observed in various diverse phenomena such as electrodeposition Kessler et al. (1988); Ben-Jacob and Garik (1990); Lopez-Tomas et al. (1995); Léger et al. (1999a); Sagués et al. (2000); Rosso et al. (2003); Schwarzacher (2004); Rosso (2007), solidification Langer (1980); Kessler et al. (1988); Ben-Jacob and Garik (1990); Cross and Hohenberg (1993) and morphogenesis Turing (1952); Cross and Hohenberg (1993). Some particular examples of morphological stability analysis include the Saffman-Taylor instability (viscous fingering) Saffman and Taylor (1958); Saffman (1986); Bensimon et al. (1986); Homsy (1987), viscous fingering coupled with electrokinetic effects Mirzadeh and Bazant (2017), the Mullins-Sekerka instability of a spherical particle during diffusion-controlled or thermally controlled growth Mullins and Sekerka (1963) and of a planar interface during solidification of a dilute binary alloy Mullins and Sekerka (1964); Sekerka (1965), and control of phase separation using electro-autocatalysis or electro-autoinhibition in driven open electrochemical systems Bazant (2013, 2017).

I.1 Stability of metal electrodeposition

We focus on electrodeposition as a specific example of an electrochemical system for which morphological stability has been widely researched both theoretically and experimentally. The fundamental aspect of electrodeposition concerns the inherent instability of the governing physics while the practical aspect is about applications such as electroplating of metals and charging of metal batteries. To elucidate the physics behind electrodeposition, in liquid electrolytes, the morphologies of electrodeposits formed and their transitions for metals such as copper, zinc and silver are particularly well studied Ben-Jacob and Garik (1990); Lopez-Tomas et al. (1995); Léger et al. (1999a); Sagués et al. (2000); Rosso et al. (2003); Rosso (2007); Schneider et al. (2017). Depending on conditions such as applied current, applied voltage and electrolyte concentration, a variety of morphological patterns such as diffusion-limited aggregation (DLA) patterns Witten and Sander (1981, 1983); Meakin (1983); Vicsek (1984); Brady and Ball (1984); Matsushita et al. (1984); Ben-Jacob et al. (1985); Sawada et al. (1986); Grier et al. (1986); Trigueros et al. (1991); Erlebacher et al. (1993), dense branching morphologies (DBM) Ben-Jacob et al. (1985); Sawada et al. (1986); Grier et al. (1986); Ben-Jacob et al. (1986, 1987); Grier et al. (1987); Ben-Jacob et al. (1988); Garik et al. (1989); Trigueros et al. (1991); Fleury et al. (1991a); Erlebacher et al. (1993); Elezgaray et al. (2000); Léger et al. (2000a) and dendritic structures Ben-Jacob et al. (1983, 1984a, 1984b, 1985); Sawada et al. (1986); Grier et al. (1986); Ben-Jacob et al. (1987, 1988); Trigueros et al. (1991) have been observed. Ion concentration fields Argoul et al. (1996); Léger et al. (1997, 1998, 1999b, 2000b); Rosso et al. (2002), electroconvection Fleury et al. (1992, 1993, 1994); Rosso et al. (1994); Huth et al. (1995), gravity-induced convection (buoyancy) Rosso et al. (1994); Huth et al. (1995); Chazalviel et al. (1996) and the presence of impurities Fleury et al. (1991b, a); Kuhn and Argoul (1993, 1994a, 1994b) have also been examined to determine their significant effects on morphology. While the range of possible morphologies of electrodeposits is diverse, for electroplating of metals, it is desirable to electrodeposit layers that are as uniform and homogeneous as possible.

Electrodeposition is also a critical process in the development of rechargeable/secondary lithium metal batteries (LMBs) that use lithium metal for the negative electrode. For negative electrodes that use lithium chemistry, because lithium metal has the lowest standard reduction potential ($-3.04\,\textnormal{V}$ vs. the standard hydrogen electrode), highest theoretical specific ($3860\,\textnormal{mAh}/\textnormal{g}$) and volumetric ($2061\,\textnormal{mAh}/\textnormal{cm}^{3}$) capacities, and lowest mass density ($0.53\,\textnormal{g}/\textnormal{cm}^{3}$) out of all possible candidates, it is the most promising choice for achieving high energy densities in LMBs Brandt (1994); Tarascon and Armand (2001); Aurbach et al. (2002); Xu (2004); Armand and Tarascon (2008); Cairns and Albertus (2010); Scrosati and Garche (2010); Etacheri et al. (2011); Xu et al. (2014); Li et al. (2014); Tu et al. (2015); Tikekar et al. (2016a); Choi and Aurbach (2016); Blomgren (2017); Lin et al. (2017a); Cheng et al. (2017); Lin et al. (2017b); Albertus et al. (2018); Lu et al. (2018). However, the charging of LMBs is equivalent to lithium electrodeposition at the negative electrode, which is an inherently unstable process that can cause the formation of dendrites that penetrate the separator and result in internal short circuits and thermal runaway during charging Brandt (1994); Tarascon and Armand (2001); Aurbach et al. (2002); Xu (2004); Armand and Tarascon (2008); Cairns and Albertus (2010); Scrosati and Garche (2010); Etacheri et al. (2011); Xu et al. (2014); Li et al. (2014); Tu et al. (2015); Tikekar et al. (2016a); Choi and Aurbach (2016); Blomgren (2017); Lin et al. (2017a); Cheng et al. (2017); Lin et al. (2017b); Albertus et al. (2018); Lu et al. (2018). This process has been especially well investigated in lithium polymer batteries that use a polymer electrolyte Brissot et al. (1998, 1999a, 1999b, 2001); Rosso et al. (2001); Teyssot et al. (2005); Rosso et al. (2006a); Stone et al. (2012). Detailed studies of various growth modes of lithium in liquid electrolytes during charging have been recently performed Bai et al. (2016); Kushima et al. (2017); Bai et al. (2018), which will aid in the development of better models for lithium electrodeposition. Modern lithium-ion batteries (LIBs) Brandt (1994); Tarascon and Armand (2001); Aurbach et al. (2002); Whittingham (2004); Xu (2004); Armand and Tarascon (2008); Bruce et al. (2008); Cairns and Albertus (2010); Scrosati and Garche (2010); Goodenough and Kim (2010); Etacheri et al. (2011); Goodenough and Park (2013); Li et al. (2014); Blomgren (2017); Lu et al. (2018) partially mitigate this problem of dendrite formation and propagation by employing lithium intercalating materials such as graphite for the negative electrode where only lithium ions and not reduced lithium atoms are involved in the intercalation reactions, which is also known as the “rocking chair technology” Tarascon and Armand (2001). Nonetheless, lithium plating still occurs in LIBs when they are charged at high rates or low temperatures Vetter et al. (2005); Goodenough and Kim (2010); Goodenough and Park (2013); Li et al. (2014); Palacín and Guibert (2016). Although the root causes of the widely publicized LIB failures in two Boeing 787 Dreamliners in January 2013 were not conclusively identified Williard et al. (2013), there is no doubt that safety is of paramount importance in both LMBs and LIBs, which requires a thorough understanding of dendrite formation.

For both electroplating of metals and charging of high energy density LMBs, it would be advantageous to perform them at as large a current or voltage as possible without causing dendrite formation. It is therefore important to understand the possible mechanisms for the electrochemical system to sustain a high current or voltage and how these mechanisms interact with the metal electrodeposition and LMB charging processes. In a neutral channel or porous medium containing an electrolyte, when ion transport is governed by diffusion and electromigration, which is collectively termed electrodiffusion, the maximum current that can be attained by the electrochemical system is called the diffusion-limited current Newman and Thomas-Alyea (2004); Bard and Faulkner (2000). In practice, overlimiting current (OLC) beyond the electrodiffusion limit has been observed experimentally in ion-exchange membranes Rubinstein and Shtilman (1979); Rösler et al. (1992); Krol et al. (1999a, b); Rubinshtein et al. (2002); Rubinstein et al. (2008); Deng et al. (2013); Schlumpberger et al. (2015); Schlumpberger (2016); Nikonenko et al. (2010, 2014); Strathmann (2010) and microchannels and nanochannels Kim et al. (2007); Yossifon and Chang (2008); Zangle et al. (2009, 2010a, 2010b); Nam et al. (2015); Schiffbauer et al. (2015); Sohn et al. (2018). Depending on the length scale of the pores or channel, some possible physical mechanisms for OLC Dydek et al. (2011) are surface conduction Zangle et al. (2009, 2010a, 2010b); Mani et al. (2009); Mani and Bazant (2011); Dydek and Bazant (2013), electroosmotic flow Yaroshchuk et al. (2011); Rubinstein and Zaltzman (2013) and electroosmotic instability Rubinstein and Zaltzman (2000); Zaltzman and Rubinstein (2007). Some chemical mechanisms for OLC include water splitting Nikonenko et al. (2010, 2014) and current-induced membrane discharge Andersen et al. (2012). In this paper, we focus on porous media consisting of pores with a nanometer length scale, therefore the dominant OLC mechanism is expected to be surface conduction Dydek et al. (2011). When a sufficiently large current or voltage is applied across a porous medium whose pore surfaces are charged, the bulk electrolyte eventually gets depleted at an ion-selective interface such as an electrode. In order to sustain the current beyond the electrodiffusion limit, the counterions in the electric double layers (EDLs) next to the charged pore surfaces migrate under the large electric field generated in the depletion region. This phenomenon is termed surface conduction and results in the formation and propagation of deionization shocks away from the ion-selective interface in porous media Mani and Bazant (2011); Dydek and Bazant (2013); Yaroshchuk (2012) and microchannels and nanochannels Zangle et al. (2009, 2010a, 2010b); Dydek et al. (2011); Mani et al. (2009); Nielsen and Bruus (2014). The deionization shock separates the “front” electrolyte-rich region, in which bulk electrodiffusion dominates, from the “back” electrolyte-poor region, in which electromigration in the EDLs dominates.

I.2 Theories of pattern formation

Morphological stability analysis of electrodeposition is typically performed in the style of the pioneering Mullins-Sekerka stability analysis Mullins and Sekerka (1963, 1964). The destabilizing effect arises from the amplification of surface protrusions by diffusive fluxes while the main stabilizing effect arises from the surface energy penalty incurred in creating additional surface area. The balance between these two effects, which is influenced by system parameters, sets a characteristic length scale or wavenumber for the surface instability. In 1963, by applying an infinitesimally small spherical harmonic perturbation to the surface of a spherical particle undergoing growth by solute diffusion or heat diffusion, Mullins and Sekerka derived a dispersion relation that related growth rates of the eigenmodes to particle radius and degree of supersaturation Mullins and Sekerka (1963). Similarly, in 1964, Mullins and Sekerka imposed a infinitesimally small sinusoidal perturbation on a planar liquid-solid interface during the solidification of a dilute binary alloy and obtained a dispersion relation relating the surface perturbation growth rate to system parameters such as temperature and concentration gradients Mullins and Sekerka (1964). In the spirit of the Mullins-Sekerka stability analysis, about 16 years later in 1980, Aogaki, Kitazawa, Kose and Fueki applied linear stability analysis to study electrodeposition with a steady-state base state in the presence of a supporting electrolyte, i.e., electromigration of the minor species can be ignored, and without explicitly accounting for electrochemical reaction kinetics Aogaki et al. (1980). Following up on this work, from 1981 to 1982, Aogaki and Makino changed the steady-state base state to a time-dependent base state under galvanostatic conditions while keeping other assumptions intact Aogaki and Makino (1981); Aogaki (1982a, b). In 1984, Aogaki and Makino extended their previous work to account for surface diffusion of adsorbed metal atoms under galvanostatic Aogaki and Makino (1984a, b) and potentiostatic conditions Aogaki and Makino (1984c, d). In the same year, Makino, Aogaki and Niki also used such a linear stability analysis to extract surface parameters of metals under galvanostatic and potentiostatic conditions Makino et al. (1984a) and applied it to study how hydrogen adsorption affects these extracted parameters under galvanostatic conditions Makino et al. (1984b). Later work by Barkey, Muller and Tobias in 1989 Barkey et al. (1989a, b), and Chen and Jorne in 1991 Chen and Jorne (1991) additionally assumed the presence of a diffusion boundary layer next to the electrode.

Subsequent developments in linear stability analysis of electrodeposition relaxed some assumptions made in the past literature and added more physics and electrochemistry. Butler-Volmer reaction kinetics was first explicitly considered by Pritzker and Fahidy in 1992 for a steady-state base state with a diffusion boundary layer next to the electrode Pritzker and Fahidy (1992). Also considering Butler-Volmer reaction kinetics with a steady-state base state, in 1995, Sundström and Bark used the Nernst-Planck equations for ion transport without assuming the existence of a diffusion boundary layer, numerically solved for the dispersion relation and performed extensive parameter sweeps over key parameters of interest such as surface energy and exchange current density Sundström and Bark (1995). Extending these two papers in 1998, Elezgaray, Léger and Argoul used a time-dependent base state under galvanostatic conditions, numerically solved for both the time-dependent base state and perturbed state to obtain the dispersion relation and demonstrated good agreement between their theory and experiments on copper electrodeposition in a thin gap cell Elezgaray et al. (1998). The role of electrolyte additives in stabilizing electrodeposition was examined in the linear stability analysis performed by Haataja, Srolovitz and Bocarsly in 2002 and 2003 Haataja and Srolovitz (2002); Haataja et al. (2003a, b), and McFadden et al. in 2003 McFadden et al. (2003). By demonstrating that the effects of the anode can be ignored under certain conditions when deriving the dispersion relation, BuAli, Johns and Narayanan in 2006 simplified Sundström and Bark’s analysis to obtain an analytical expression for the dispersion relation BuAli et al. (2006). In 2004 and 2005, Monroe and Newman included additional mechanical effects such as pressure, viscous stress and deformational stress to the linear stability analysis of electrodeposition, which provided more stabilization beyond that provided by surface energy Monroe and Newman (2004, 2005). For a steady-state base state, in 2014, Tikekar, Archer and Koch studied how tethered immobilized anions provide additional stabilization to electrodeposition by reducing the electric field at the cathode and, after making some approximations, derived analytical expressions for the dispersion relation for small and large current densities Tikekar et al. (2014). Tikekar, Archer and Koch then extended this work in 2016 by accounting for elastic deformations that provide further stabilization Tikekar et al. (2016b). Subsequently in 2018, Tikekar, Li, Archer and Koch looked at using polymer electrolyte viscosity to suppress morphological instabilities driven by electroconvection Tikekar et al. (2018). Building on Monroe and Newman’s 2004 and 2005 work on interfacial deformation effects Monroe and Newman (2004, 2005), Ahmad and Viswanathan identified a new mechanism for stabilization driven by the difference of the metal density in the metal electrode and solid electrolyte in 2017 Ahmad and Viswanathan (2017a), and further generalized this work in the same year to account for anisotropy Ahmad and Viswanathan (2017b). Natsiavas, Weinberg, Rosato and Ortiz in 2016 also investigated the stabilizing effect of prestress and showed good agreement between theory and experiment Natsiavas et al. (2016). Relaxing the usual assumption of electroneutrality, in 2015, Nielsen and Bruus performed linear stability analysis for a steady-state base state that accounts for the extended space charge region that is formed when the electric double layer becomes nonequilibrium at an overlimiting current Nielsen and Bruus (2015).

Without performing a linear stability analysis, some models focus on describing the initiation and subsequent propagation of dendrites. The classic work in this class of models is by Chazalviel in 1990 in which they used the Poisson’s equation for electrostatics, i.e., electroneutrality is not assumed, and showed that the initiation of ramified electrodeposits is caused by the creation of a space charge layer upon anion depletion at the cathode, the induction time for initiation is the time needed for building up this space charge layer, and the velocity of the ramified growth is equal to the electromigration velocity of the anions Chazalviel (1990); some experimental results were also obtained by Fleury, Chazalviel, Rosso and Sapoval in support of this model Fleury et al. (1990), and some of the numerical results of the original analysis were subsequently improved by Rosso, Chazalviel and Chassaing Rosso et al. (2006b). Via an asymptotic analysis of the Poisson-Nernst-Planck equations for ion transport, Bazant also showed that the velocity of the ramified growth is approximately equal to the anion electromigration velocity and estimated the induction time for the onset of ramified growth Bazant (1995). Building on past theoretical and experimental work done on silver electrodeposition by Barton and Bockris Barton and Bockris (1962), and zinc electrodeposition by Despic, Diggle and Bockris Despic et al. (1968); Diggle et al. (1969), Monroe and Newman investigated the propagation velocity and length of a dendrite tip that grows via Butler-Volmer kinetics Monroe and Newman (2003). By examining the thermodynamics and kinetics of heterogeneous nucleation and growth, which is assumed to proceed via the linearized Butler-Volmer equation valid for small overpotentials, Ely and García identified five different regimes of nucleus behavior Ely and García (2013). Assuming a concentration-dependent electrolyte diffusivity and the existence of a hemispherical dendrite “precursor” that grows with Tafel kinetics, Akolkar studied the subsequent propagation velocity and length of the dendrite Akolkar (2013) and how they are affected by temperature Akolkar (2014).

I.3 Contributions of this work

In this paper, we perform linear stability analysis of electrodeposition inside a charged random porous medium, whose pore surface charges can generally be of any sign, that is filled with a liquid electrolyte and flanked on its sides by a pair of planar metal electrodes. The linear stability analysis is carried out with respect to a time-dependent base state and focuses on overlimiting current carried by surface conduction. By doing so, we combine and generalize previous work done in Sundström and Bark (1995); Elezgaray et al. (1998); Tikekar et al. (2014). For simplicity, we ignore bulk convection, electroosmotic flow, surface adsorption, surface diffusion of adsorbed species Aogaki and Makino (1984a, b, c, d) and additional mechanical effects such as pressure, viscous stress and deformational stress Monroe and Newman (2004, 2005); Tikekar et al. (2016b); Ahmad and Viswanathan (2017a, b); Natsiavas et al. (2016). We expect surface diffusion of adsorbed species, which alleviates electrodiffusion limitations, and interfacial deformation effects to stabilize electrodeposition, hence our work here can be considered a worst-case analysis. The only electrochemical reaction considered here is metal electrodeposition, therefore in the context of LMBs and LIBs, electrochemical and chemical reactions between lithium and the electrolyte that cause the formation of the solid electrolyte interphase (SEI) layer Verma et al. (2010); Cheng et al. (2016); Peled and Menkin (2017); Wang et al. (2018) are not included. We first derive governing equations for the full model that consists of coupling ion transport with electrochemical reaction kinetics, followed by applying linear stability analysis on the full model via the imposition of sinusoidal spatial perturbations around the time-dependent base state. For the dispersion relation, we perform a boundary layer analysis on the perturbed state to derive an accurate approximation for it and a convergence analysis of its full numerical solution. To better understand the physics of the dispersion relation, we carry out parameter sweeps over the pore surface charge density, Damköhler number and applied current density under galvanostatic conditions. We also compare the numerical and approximate solutions of the maximum wavenumber, maximum growth rate and critical wavenumber in order to verify the accuracy of these approximations. Subsequently, we apply the linear stability analysis to compare theoretical predictions and experimental data for copper electrodeposition in cellulose nitrate membranes Han et al. (2016a), and also use the stability analysis as a tool for investigating the dependence of crystal grain size on duty cycle during pulse electroplating.

II Full model

II.1 Transport

The starting point for modeling ion transport is the leaky membrane model that is able to predict overlimiting current carried by surface conduction, which we have previously coupled with Butler-Volmer reaction kinetics for analyzing steady state current-voltage relations and linear sweep voltammetry Khoo and Bazant (2018). The system under consideration is a binary asymmetric liquid electrolyte in a finite 3D charged random nanoporous medium where x\in\mathopen{}\mathclose{{}\left[0,L_{x}}\right], y\in\mathopen{}\mathclose{{}\left[0,L_{y}}\right] and z\in\mathopen{}\mathclose{{}\left[0,L_{z}}\right], whose 2D projection is illustrated in Figure 1. We assume that the nanoporous medium is random with well connected pores such as cellulose nitrate membranes so that we can investigate macroscopic electrode-scale morphological instabilities Han et al. (2016a). The cations are electroactive and the anions are inert. Initially, at $t=0$, the anode surface is located at $x=0$ while the cathode surface is located at $x=L_{x}$. As is typical for linear stability analysis of electrodeposition Sundström and Bark (1995); Elezgaray et al. (1998); Nielsen and Bruus (2015), we pick a moving reference frame with a velocity u\mathopen{}\mathclose{{}\left(t}\right)=u_{x}\mathopen{}\mathclose{{}\left(t}\right)e_{x} that is equal to the velocity of the electrode/electrolyte interface so that the average positions of the dissolving anode and growing cathode remain stationary. For the porous medium, we denote its porosity, tortuosity, internal pore surface area per unit volume and pore surface charge per unit area as $\epsilon_{\textnormal{p}}$, $\tau$, $a_{\textnormal{p}}$ and $\sigma_{\textnormal{s}}$ respectively. We assume that there are no homogeneous reactions and all material properties such as $\epsilon_{\textnormal{p}}$, $\tau$, $a_{\textnormal{p}}$ and $\sigma_{\textnormal{s}}$ are constant and uniform. We also assume that dilute solution theory holds true, hence all activity coefficients are $1$ and the cation and anion macroscopic diffusivities $D_{\pm 0}$, where the [math] subscript indicates dilute limit, are constant, uniform and independent of concentrations. Accounting for corrections due to the tortuosity of the porous medium, the macroscopic diffusivity $D_{\pm 0}$ is related to the molecular (free solution) diffusivity $D_{\pm 0}^{\textnormal{m}}$ by $D_{\pm 0}=\frac{D_{\pm 0}^{\textnormal{m}}}{\tau}$ Ferguson and Bazant (2012). The assumption of dilute solution theory further implies that the convective flux in the moving reference frame is negligible and the effect of the moving reference frame velocity u\mathopen{}\mathclose{{}\left(t}\right)=u_{x}\mathopen{}\mathclose{{}\left(t}\right)e_{x} is only significant in the equation describing electrode surface growth or dissolution Sundström and Bark (1995); Elezgaray et al. (1998); Nielsen and Bruus (2015), which we will discuss in Section II.3. Under these assumptions, for ion transport, the Nernst-Planck equations describing species conservation, charge conservation equation and macroscopic electroneutrality constraint are given by

[TABLE]

where $c_{\pm}$, $F_{\pm}$, $z_{\pm}$, are the cation and anion concentrations, fluxes and charge numbers respectively, $\phi$ is the electrolyte electric potential, $J$ is the electrolyte current density, $h_{\textnormal{p}}=\frac{\epsilon_{\textnormal{p}}}{a_{\textnormal{p}}}$ is the effective pore size and $\rho_{\textnormal{s}}$ is the volume-averaged background charge density. Denoting the numbers of cations and anions that are formed from complete dissociation of $1$ molecule of neutral salt as $\nu_{\pm}$, electroneutrality requires that $z_{+}\nu_{+}+z_{-}\nu_{-}=0$. We will use the macroscopic electroneutrality constraint given by Equation 3 to eliminate $c_{+}$ as a dependent variable, therefore leaving $c_{-}$ and $\phi$ as the remaining dependent variables.

For a classical system with an uncharged nanoporous medium, i.e., $\rho_{\textnormal{s}}=0$, the maximum current density that the system can possibly attain under electrodiffusion is called the diffusion-limited or limiting current density $J_{\textnormal{lim}}$, which is given by Khoo and Bazant (2018)

[TABLE]

where $c_{0}$ is the neutral salt bulk concentration. The limiting current $I_{\textnormal{lim}}$ is then given by $I_{\textnormal{lim}}=J_{\textnormal{lim}}A$ where $A=L_{y}L_{z}$ is the surface area of the electrode. Under galvanostatic conditions, when a current density $J_{\textnormal{a}}$ larger than $J_{\textnormal{lim}}$ is applied, the cation and anion concentrations at the cathode reach [math] and the electrolyte electric potential and electric field there diverge in finite time, which is called Sand’s time Sand (1899) denoted as $t_{\textnormal{s}}$; see van Soestbergen et al. (2010) for a discussion of some subtlety associated with this transition time when $J_{\textnormal{a}}$ is exactly equal to $J_{\textnormal{lim}}$. If we have not assumed macroscopic electroneutrality and instead modeled electric double layers that can go out of equilibrium at high currents and voltages, then the electric field would be large but finite at and past $t_{\textnormal{s}}$ Zaltzman and Rubinstein (2007); Nielsen and Bruus (2015); Bazant et al. (2005); Chu and Bazant (2005). Defining the dimensionless Sand’s time $\tilde{t}_{\textnormal{s}}=\frac{D_{\textnormal{amb0}}t_{\textnormal{s}}}{L_{x}^{2}}$ and dimensionless applied current density $\tilde{J}_{\textnormal{a}}=\frac{J_{\textnormal{a}}}{J_{\textnormal{lim}}}$ where D_{\textnormal{amb0}}=\frac{\mathopen{}\mathclose{{}\left(z_{+}-z_{-}}\right)D_{+0}D_{-0}}{z_{+}D_{+0}-z_{-}D_{-0}} Newman and Thomas-Alyea (2004); Khoo and Bazant (2018) is the ambipolar diffusivity of the neutral salt in the dilute limit and $\frac{L_{x}^{2}}{D_{\textnormal{amb0}}}$ is the diffusion time scale, $\tilde{t}_{\textnormal{s}}$ is given by van Soestbergen et al. (2010)

[TABLE]

For galvanostatic conditions, $t_{\textnormal{s}}$ is a critically important time scale because the formation of dendrites often occurs near or at $t_{\textnormal{s}}$, therefore it will be central to the linear stability analysis results discussed in Section IV.

Unlike the classical case of $\rho_{\textnormal{s}}=0$, for $\rho_{\textnormal{s}}<0$, the system can sustain an overlimiting current $\tilde{J}_{\textnormal{a}}>1$ via surface conduction that is the electromigration of the counterions in the electric double layers (EDLs), which are formed next to the charged pore surfaces, under a large electric field generated in the depletion region next to the cathode. This additional surface conductivity thus enables the system to go beyond $t_{\textnormal{s}}$ and eventually reach a steady state, in stark contrast to the finite time divergence of the classical case at $t_{\textnormal{s}}$. On the other hand, for $\rho_{\textnormal{s}}>0$, the counterions in the EDLs are the inert anions instead of the electroactive cations, which contribute a surface current that flows in the opposite direction from that of the bulk current. Because of this “negative” surface conductivity conferred by $\rho_{\textnormal{s}}>0$ relative to $\rho_{\textnormal{s}}=0$, at the cathode, the bulk electrolyte concentration vanishes and the electric field diverges earlier than $t_{\textnormal{s}}$; in other words, $\rho_{\textnormal{s}}>0$ effectively reduces $t_{\textnormal{s}}$. A more quantitative way of explaining this is that the “negative” surface conductivity causes the maximum current density that can be achieved, which is denoted as $J_{\textnormal{max}}$, to be smaller than $J_{\textnormal{lim}}$, and $J_{\textnormal{max}}$ decreases as $\rho_{\textnormal{s}}$ increases. In effect, a more positive $\rho_{\textnormal{s}}$ decreases $J_{\textnormal{lim}}$, which thus leads to a smaller $t_{\textnormal{s}}$ for a given $J_{\textnormal{a}}$ according to Equation 5. Details of how to numerically compute $J_{\textnormal{max}}$ are found in Khoo and Bazant (2018); note that $J_{\textnormal{max}}$ here corresponds to $I_{\textnormal{max}}^{\textnormal{BV}}$ in Khoo and Bazant (2018).

II.2 Electrochemical reaction kinetics

In order to analyze how spatial perturbations of an electrode surface affect its linear stability, we need to account for the effects of surface curvature and surface energy in the electrochemical reaction kinetics model. The mean curvature of the electrode/electrolyte interface $\mathscr{H}$ is given by $\mathscr{H}=-\frac{1}{2}\nabla_{\textnormal{s}}\cdot\hat{n}$ where $\nabla_{\textnormal{s}}$ is the surface gradient operator and $\hat{n}$ is the unit normal that points outwards from the electrolyte Deen (2011). In this paper, we consider a charge transfer reaction that involves only the cations and electrons while the anions are inert. More concretely, we suppose the following charge transfer reaction consuming $n$ electrons $\textnormal{O}^{z_{\textnormal{O}}}+n\textnormal{e}^{-}\rightleftharpoons\textnormal{R}^{z_{\textnormal{R}}}$ where $\textnormal{O}^{z_{\textnormal{O}}}$ is the oxidized species O with charge $z_{\textnormal{O}}$, $\textnormal{e}^{-}$ is the electron e with charge $-1$, $\textnormal{R}^{z_{\textnormal{R}}}$ is the reduced species R with charge $z_{\textnormal{R}}$, and $z_{\textnormal{O}}-n=z_{\textnormal{R}}$ because of charge conservation. If the reduced species is solid metal, i.e., $z_{\textnormal{R}}=0$, as is the case in metal electrodeposition, the creation of additional electrode/electrolyte interfacial area results in a surface energy penalty that appears in the electrochemical potential of the reduced species. Therefore, the electrochemical potentials $\mu_{i}$ of the oxidized species O, electron e and reduced species R for i\in\mathopen{}\mathclose{{}\left\{\textnormal{O},\textnormal{e},\textnormal{R}}\right\} are given by

[TABLE]

where the surface energy term $2\Omega\gamma\mathscr{H}$ Sundström and Bark (1995); Elezgaray et al. (1998); Monroe and Newman (2003, 2004, 2005); Deen (2011); Tikekar et al. (2014); Nielsen and Bruus (2015) is included in $\mu_{\textnormal{R}}$ when the reduced species is solid metal ($z_{\textnormal{R}}=0$) and the $\Theta$ superscript indicates standard state. The activity of species $i$ is given by $a_{i}=\gamma_{i}\hat{c}_{i}$ where $\gamma_{i}$ is the activity coefficient of species $i$ and $\hat{c}_{i}=\frac{c_{i}}{c_{i}^{\Theta}}$ is the concentration of species $i$ normalized by its standard concentration $c_{i}^{\Theta}$. $\mu_{i}^{\Theta}$ is the standard electrochemical potential of species $i$, $\phi_{\textnormal{e}}$ is the electrode electric potential, $\Omega=\frac{M_{\textnormal{m}}}{\rho_{\textnormal{m}}}$ is the atomic volume of the solid metal where $M_{\textnormal{m}}$ and $\rho_{\textnormal{m}}$ are the atomic mass and mass density of the metal respectively, and $\gamma$ is the isotropic surface energy of the metal/electrolyte interface. The quantity $\frac{\Omega\gamma}{k_{\textnormal{B}}T}$ is the capillary constant that has units of length Mullins and Sekerka (1963, 1964); Sekerka (1965). The interfacial electric potential difference $\Delta\phi$ is defined as $\Delta\phi=\phi_{\textnormal{e}}-\phi$. At equilibrium when $\mu_{\textnormal{O}}+n\mu_{\textnormal{e}}=\mu_{\textnormal{R}}$, we obtain the Nernst equation

[TABLE]

where the “eq” superscript denotes equilibrium and $E^{\Theta}$ is the standard electrode potential. When the system is driven out of equilibrium so that $\mu_{\textnormal{O}}+n\mu_{\textnormal{e}}\neq\mu_{\textnormal{R}}$, the system generates a Faradaic current density $J_{\textnormal{F}}$ that is given by Bazant (2013, 2017); Ferguson and Bazant (2012)

[TABLE]

where $k_{0}$ is the overall reaction rate constant and $\mu_{\ddagger}^{\textnormal{r,ex}}$ is the excess electrochemical potential of the transition state for the Faradaic reaction. Using the Butler-Volmer hypothesis, $\mu_{\ddagger}^{\textnormal{r,ex}}$ consists of a chemical contribution $k_{\textnormal{B}}T\ln\gamma_{\ddagger}^{\textnormal{r}}$, where $\gamma_{\ddagger}^{\textnormal{r}}$ is the activity coefficient of the transition state for the Faradaic reaction, and a convex combination of the electrostatic energies, surface energies (only for the reduced species) and standard electrochemical potentials weighted by $\alpha$, which is the charge transfer coefficient. Therefore, $\mu_{\ddagger}^{\textnormal{r,ex}}$ is given by

[TABLE]

Defining the overpotential $\eta$ as $\eta=\frac{\mu_{\textnormal{R}}-\mu_{\textnormal{O}}-n\mu_{\textnormal{e}}}{ne}=\Delta\phi-\Delta\phi^{\textnormal{eq}}$, we rewrite $J_{\textnormal{F}}$ as

[TABLE]

where $j_{0}$ is the exchange current density. In this form, we can identify the cathodic and anodic charge transfer coefficients, which are denoted as $\alpha_{\textnormal{c}}$ and $\alpha_{\textnormal{a}}$ respectively, as $\alpha_{\textnormal{c}}=\alpha$ and $\alpha_{\textnormal{a}}=1-\alpha$ such that $\alpha_{\textnormal{c}}+\alpha_{\textnormal{a}}=1$. We note that our particular choice of $\mu_{\ddagger}^{\textnormal{r,ex}}$ in Equation 11 corresponds to choosing the “mechanical transfer coefficient” $\alpha_{\textnormal{m}}$ defined in Monroe and Newman (2004) to be equal to $\alpha_{\textnormal{a}}$, causing $j_{0}$ to not depend explicitly on $\mathscr{H}$.

In this paper, we assume that the only charge transfer reaction occurring is metal electrodeposition that happens via the electrochemical reduction of cations in the electrolyte to solid metal on the electrode. The activity of solid metal is $1$ while we assume that the activity of electrons is $1$. In addition, like in Section II.1, we assume that dilute solution theory is applicable, therefore the activity coefficients of the cation, anion and transition state for the Faradaic reaction are $1$ and we replace activities of the cation and anion with their normalized concentrations $\hat{c}_{\pm}$. Therefore, $\Delta\phi^{\textnormal{eq}}$ and $j_{0}$ simplify to

[TABLE]

To compare the reaction and diffusion rates, we define the Damköhler number Da as the ratio of the Faradaic current density scale $e\epsilon_{\textnormal{p}}k_{0}$ and limiting current density $J_{\textnormal{lim}}$ given by

[TABLE]

When Da is large, i.e., $\textnormal{Da}\gg 1$, the system is diffusion-limited but when Da is small, i.e., $\textnormal{Da}\ll 1$, the system is reaction-limited.

II.3 Boundary conditions, constraints and initial conditions

We use “a” and “c” superscripts to denote the anode and cathode respectively, r=\mathopen{}\mathclose{{}\left[x,y,z}\right]^{\mathrm{T}} to denote the position vector, and r_{\textnormal{m}}^{\textnormal{a,c}}=\mathopen{}\mathclose{{}\left[x_{\textnormal{m}}^{\textnormal{a,c}},y_{\textnormal{m}}^{\textnormal{a,c}},z_{\textnormal{m}}^{\textnormal{a,c}}}\right]^{\mathrm{T}} to denote the positions of the anode and cathode. We ground the anode at all times, i.e., $\phi_{\textnormal{e}}^{\textnormal{a}}=0$. Because the cations are electroactive, we impose no-flux boundary conditions for the cation flux on all boundaries except the anode and cathode where Faradaic reactions involving the cations occur. On the other hand, because the anions are inert, we impose no-flux boundary conditions for the anion flux on all boundaries. At the anode and cathode, we require the conservation of charges across the electrode/electrolyte interfaces. In summary, the boundary conditions are given by \hat{n}\cdot F_{-}\mathopen{}\mathclose{{}\left(r=r_{\textnormal{m}}^{\textnormal{a,c}}}\right)=0, \hat{n}\cdot J\mathopen{}\mathclose{{}\left(r=r_{\textnormal{m}}^{\textnormal{a,c}}}\right)=\epsilon_{\textnormal{p}}J_{\textnormal{F}}^{\textnormal{a,c}}, \hat{n}\cdot F_{-}\mathopen{}\mathclose{{}\left(r=r_{\textnormal{other}}}\right)=0 and \hat{n}\cdot J\mathopen{}\mathclose{{}\left(r=r_{\textnormal{other}}}\right)=0 where $r_{\textnormal{other}}$ refers to all boundaries except the anode and cathode.

The velocity of the electrode/electrolyte interface $v_{\textnormal{I}}^{\textnormal{a,c}}$ is defined as $v_{\textnormal{I}}^{\textnormal{a,c}}=\epsilon_{\textnormal{p}}\frac{\mathrm{d}r_{\textnormal{m}}^{\textnormal{a,c}}}{\mathrm{d}t}$ and its normal component $v_{\textnormal{In}}^{\textnormal{a,c}}$ is given by $v_{\textnormal{In}}^{\textnormal{a,c}}=\hat{n}\cdot v_{\textnormal{I}}^{\textnormal{a,c}}=\epsilon_{\textnormal{p}}\hat{n}\cdot\frac{\mathrm{d}r_{\textnormal{m}}^{\textnormal{a,c}}}{\mathrm{d}t}$. Because the growth (dissolution) of the electrode surface is caused by metal electrodeposition (electrodissolution), taking into account the moving reference frame velocity u\mathopen{}\mathclose{{}\left(t}\right)=u_{x}\mathopen{}\mathclose{{}\left(t}\right)e_{x}, the normal interfacial velocity is related to the normal current density by v_{\textnormal{In}}^{\textnormal{a,c}}=-\frac{\Omega\hat{n}\cdot J\mathopen{}\mathclose{{}\left(r=r_{\textnormal{m}}^{\textnormal{a,c}}}\right)}{z_{+}e}-\hat{n}\cdot u\mathopen{}\mathclose{{}\left(r=r_{\textnormal{m}}^{\textnormal{a,c}}}\right) and therefore, \epsilon_{\textnormal{p}}\hat{n}\cdot\frac{\mathrm{d}r_{\textnormal{m}}^{\textnormal{a,c}}}{\mathrm{d}t}=-\frac{\Omega\hat{n}\cdot J\mathopen{}\mathclose{{}\left(r=r_{\textnormal{m}}^{\textnormal{a,c}}}\right)}{z_{+}e}-\hat{n}\cdot u\mathopen{}\mathclose{{}\left(r=r_{\textnormal{m}}^{\textnormal{a,c}}}\right).

For galvanostatic conditions in which we apply a current $I_{\textnormal{a}}$ on the system, we require \int\hat{n}\cdot J\mathopen{}\mathclose{{}\left(r=r_{\textnormal{m}}^{\textnormal{c}}}\right)\,\mathrm{d}S^{\textnormal{c}}=\int-\hat{n}\cdot J\mathopen{}\mathclose{{}\left(r=r_{\textnormal{m}}^{\textnormal{a}}}\right)\,\mathrm{d}S^{\textnormal{a}}=I_{\textnormal{a}} to satisfy charge conservation whereas for potentiostatic conditions in which we apply an electric potential $V$ on the cathode, we set $\phi_{\textnormal{e}}^{\textnormal{c}}=V$. For initial conditions, we set c_{-}\mathopen{}\mathclose{{}\left(t=0}\right)=\nu_{-}c_{0}-\frac{\rho_{\textnormal{s}}+\mathopen{}\mathclose{{}\left\lvert\rho_{\textnormal{s}}}\right\rvert}{2z_{-}e}\equiv\beta_{1} where $c_{0}$ is the initial neutral salt bulk concentration Khoo and Bazant (2018), and x_{\textnormal{m}}^{\textnormal{a}}\mathopen{}\mathclose{{}\left(t=0}\right)=0 and x_{\textnormal{m}}^{\textnormal{c}}\mathopen{}\mathclose{{}\left(t=0}\right)=L_{x}, i.e., the anode and cathode are initially planar.

III Linear stability analysis

III.1 Perturbations and linearization

Linear stability analysis generally involves imposing a spatial perturbation around a base state, keeping constant and linear terms of the perturbed state, and determining the dispersion relation that relates the growth rate of the perturbation to its wavenumber or wavelength. For electrodeposition specifically, the objective is to impose a spatial perturbation on a planar electrode surface and determine the effects of key system parameters on the linear stability of the surface in response to this perturbation. In this paper, we will choose a time-dependent base state, therefore the dispersion relation is also time-dependent. In 3D, the electrode/electrolyte interface can be written explicitly as x=h\mathopen{}\mathclose{{}\left(y,z,t}\right) where $h$ is the electrode surface height. Given $h$, we can derive explicit expressions for surface variables such as the curvature $\mathscr{H}$ and normal interfacial velocity $v_{\textnormal{In}}$ in terms of $h$ and its spatial and temporal derivatives Deen (2011); Crank (1987), which are provided in Section I of Supplementary Material. For brevity, we let k=\mathopen{}\mathclose{{}\left[k_{y},k_{z}}\right]^{\mathrm{T}} and \xi=\mathopen{}\mathclose{{}\left[y,z}\right]^{\mathrm{T}} where $k$ is the wavevector, and $k_{y}$ and $k_{z}$ are the wavenumbers in the $y$ and $z$ directions respectively. Therefore, $k\cdot\xi=k_{y}y+k_{z}z$, k^{2}=\mathopen{}\mathclose{{}\left\lVert k}\right\rVert_{2}^{2}=k_{y}^{2}+k_{z}^{2} where \mathopen{}\mathclose{{}\left\lVert\cdot}\right\rVert_{2} is the $L^{2}$-norm and \mathopen{}\mathclose{{}\left\lVert k}\right\rVert_{2} is the overall wavenumber, and the wavelength $\lambda$ is given by \lambda=\frac{2\pi}{\mathopen{}\mathclose{{}\left\lVert k}\right\rVert_{2}}. For brevity again, we write the overall wavenumber as $k$ and it is obvious from context whether $k$ refers to the wavevector or overall wavenumber. The perturbation that will be imposed is sinusoidal in the $y$ and $z$ directions given by

[TABLE]

where $\epsilon\ll 1$ is a dimensionless small parameter, the \mathopen{}\mathclose{{}\left(0}\right) and \mathopen{}\mathclose{{}\left(1}\right) superscripts denote the base and perturbed states respectively, \Re\mathopen{}\mathclose{{}\left(\cdot}\right) gives the real part of a complex number, h^{\mathopen{}\mathclose{{}\left(1}\right)} is the complex-valued perturbation amplitude of the electrode surface height, and $\omega$ is the complex-valued growth rate of the perturbation. In response to such an electrode surface perturbation, we assume that the perturbations to $c_{-}$ and $\phi$ are similarly given by

[TABLE]

where c_{-}^{\mathopen{}\mathclose{{}\left(1}\right)} and \phi^{\mathopen{}\mathclose{{}\left(1}\right)} are the complex-valued perturbation amplitudes of anion concentration and electrolyte electric potential respectively.

To evaluate $c_{-}$ and $\phi$ and their gradients $\nabla c_{-}$ and $\nabla\phi$ at the interface at $x=h$, we require their Taylor series expansions around the base state interface at x=h^{\mathopen{}\mathclose{{}\left(0}\right)}. Letting \hat{\epsilon}=\epsilon\exp\mathopen{}\mathclose{{}\left(ik\cdot\xi+\omega t}\right) and \theta\in\mathopen{}\mathclose{{}\left\{c_{-},\phi}\right\}, these expansions are given by Sundström and Bark (1995); Elezgaray et al. (1998); Nielsen and Bruus (2015)

[TABLE]

After substituting these perturbation expressions into the full model in Section II, we obtain the base and perturbed states by matching the \mathcal{O}\mathopen{}\mathclose{{}\left(1}\right) and \mathcal{O}\mathopen{}\mathclose{{}\left(\epsilon}\right) terms respectively. The dispersion relation \omega\mathopen{}\mathclose{{}\left(k}\right) is subsequently computed by solving these \mathcal{O}\mathopen{}\mathclose{{}\left(1}\right) and \mathcal{O}\mathopen{}\mathclose{{}\left(\epsilon}\right) equations. The growth rate $\omega$ is generally complex-valued and for a particular $k$ value, there is an infinite discrete spectrum of $\omega$ values. However, for linear stability analysis, we are primarily interested in the maximum of the real parts of the $\omega$ values, which is denoted as \max\mathopen{}\mathclose{{}\left\{\Re\mathopen{}\mathclose{{}\left(\omega}\right)}\right\}, that corresponds to the most unstable eigenmode. If \max\mathopen{}\mathclose{{}\left\{\Re\mathopen{}\mathclose{{}\left(\omega}\right)}\right\}<0, the perturbation decays exponentially in time and the base state is linearly stable. On the other hand, if \max\mathopen{}\mathclose{{}\left\{\Re\mathopen{}\mathclose{{}\left(\omega}\right)}\right\}>0, the perturbation grows exponentially in time and the base state is linearly unstable. Lastly, if \max\mathopen{}\mathclose{{}\left\{\Re\mathopen{}\mathclose{{}\left(\omega}\right)}\right\}=0, the perturbation does not decay nor grow and the base state is marginally stable. The imposition of the boundary conditions at $r=r_{\textnormal{other}}$ described in Section II.3 results in the quantization of the set of valid wavenumbers, which is explained in detail in Section IIIC of Supplementary Material.

III.2 Nondimensionalization

To make the equations more compact and derive key dimensionless parameters, in Table 1, we define the scales that are used for nondimensionalizing the full model in Section II and the perturbation expressions in Section III.1. $\tilde{L}_{y}$ and $\tilde{L}_{z}$ are the aspect ratios in the $y$ and $z$ directions respectively. For convenience, we also define \beta_{\textnormal{D}}=-\frac{z_{-}D_{-0}}{2\mathopen{}\mathclose{{}\left(z_{+}D_{+0}-z_{-}D_{-0}}\right)} (weighted ratio of cation and anion diffusivities), $\beta_{\textnormal{m}}=\nu_{+}c_{0}\Omega$ (ratio of atomic volume of solid metal to reciprocal electrolyte concentration), $\beta_{\textnormal{v}}=\frac{\beta_{\textnormal{m}}}{\beta_{\textnormal{D}}}$ and $\xi_{+}=\frac{\nu_{+}c_{0}}{c_{+}^{\Theta}}$ (ratio of electrolyte concentration to standard cation concentration), and note that \hat{c}_{+}=\xi_{+}\mathopen{}\mathclose{{}\left(\tilde{c}_{-}-\tilde{\rho}_{\textnormal{s}}}\right). Two important dimensionless parameters emerge from this nondimensionalization process, namely the Damköhler number $\textnormal{Da}=\tilde{k}_{0}$ that is described earlier in Equation 14 and the capillary number Ca that is given by

[TABLE]

which is the ratio of the capillary constant $\frac{\Omega\gamma}{k_{\textnormal{B}}T}$ Mullins and Sekerka (1963, 1964); Sekerka (1965) to the inter-electrode distance $L_{x}$, and $\tilde{\gamma}$ is the dimensionless isotropic surface energy of the metal/electrolyte interface.

To avoid cluttering the notation, we drop tildes for all dimensionless variables and parameters, and all variables and parameters are dimensionless in the following sections unless otherwise stated. We also rewrite the \mathopen{}\mathclose{{}\left(0}\right) and \mathopen{}\mathclose{{}\left(1}\right) superscripts, which denote the base and perturbed states respectively, as [math] and $1$ subscripts respectively. Similarly, we drop the [math] subscript for diffusivities and the $-$ subscript for anion-related variables and parameters. As shorthand, we use subscripts to denote partial derivatives with respect to $x$ (with the exception that $u_{x}$ denotes the $x$ component of the moving reference frame velocity $u$), $y$, $z$ and $t$, primes to denote total derivatives with respect to $x$, and an overhead dot to denote the total derivative with respect to $t$. All equations for the dimensionless full model are provided in Section II of Supplementary Material. Details for deriving the dimensionless equations for the base and perturbed states are provided in Section III of Supplementary Material, and we summarize them in Sections III.3 and III.4 below.

III.3 Base state

The equations for the base state are obtained by substituting the perturbation expressions in Section III.1 into the full model in Section II and matching terms at \mathcal{O}\mathopen{}\mathclose{{}\left(1}\right). Equivalently, the base state is simply the full model specialized to 1D in the $x$ direction with the curvature-related terms dropped, which only appear at \mathcal{O}\mathopen{}\mathclose{{}\left(\epsilon}\right). Therefore, at \mathcal{O}\mathopen{}\mathclose{{}\left(1}\right), the governing PDEs (partial differential equations) are given by

[TABLE]

where the first PDE is the Nernst-Planck equation describing species conservation of anions and the second PDE is the charge conservation equation. The boundary conditions at the anode at $x=h_{0}^{\textnormal{a}}$ are given by

[TABLE]

where \eta_{0}=\phi_{\textnormal{e}}-\phi_{0}-\frac{1}{n}\ln\mathopen{}\mathclose{{}\left[\xi_{+}\mathopen{}\mathclose{{}\left(c_{0}-\rho_{\textnormal{s}}}\right)}\right]-E^{\Theta}, j_{0,0}=\textnormal{Da}n\mathopen{}\mathclose{{}\left[\xi_{+}\mathopen{}\mathclose{{}\left(c_{0}-\rho_{\textnormal{s}}}\right)}\right]^{1-\alpha} and \hat{n}\cdot J_{0}=\beta_{\textnormal{D}}\mathopen{}\mathclose{{}\left[-\mathopen{}\mathclose{{}\left(D-D_{+}}\right)c_{0,x}-z_{+}D_{+}\rho_{\textnormal{s}}\phi_{0,x}+\mathopen{}\mathclose{{}\left(z_{+}D_{+}-zD}\right)c_{0}\phi_{0,x}}\right]. Since the unit normal at the cathode points in the opposite direction from that at the anode, the signs of the expressions involving $\hat{n}$ at the cathode are opposite to that at the anode. Therefore, the boundary conditions at the cathode at $x=h_{0}^{\textnormal{c}}$ are given by

[TABLE]

where \eta_{0}=\phi_{\textnormal{e}}-\phi_{0}-\frac{1}{n}\ln\mathopen{}\mathclose{{}\left[\xi_{+}\mathopen{}\mathclose{{}\left(c_{0}-\rho_{\textnormal{s}}}\right)}\right]-E^{\Theta}, j_{0,0}=\textnormal{Da}n\mathopen{}\mathclose{{}\left[\xi_{+}\mathopen{}\mathclose{{}\left(c_{0}-\rho_{\textnormal{s}}}\right)}\right]^{1-\alpha} and \hat{n}\cdot J_{0}=\beta_{\textnormal{D}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(D-D_{+}}\right)c_{0,x}+z_{+}D_{+}\rho_{\textnormal{s}}\phi_{0,x}-\mathopen{}\mathclose{{}\left(z_{+}D_{+}-zD}\right)c_{0}\phi_{0,x}}\right].

We pick u_{x}\mathopen{}\mathclose{{}\left(x=h_{0}^{\textnormal{a}}}\right) and u_{x}\mathopen{}\mathclose{{}\left(x=h_{0}^{\textnormal{c}}}\right) such that the positions of the anode and cathode in the base state remain stationary, i.e., $\dot{h}_{0}^{\textnormal{a}}=\dot{h}_{0}^{\textnormal{c}}=0$. Therefore, u_{x}=\beta_{\textnormal{v}}\hat{n}\cdot J_{0}\mathopen{}\mathclose{{}\left(x=h_{0}^{\textnormal{a}}}\right)=-\beta_{\textnormal{v}}\hat{n}\cdot J_{0}\mathopen{}\mathclose{{}\left(x=h_{0}^{\textnormal{c}}}\right) where the second equality automatically holds true because of charge conservation in the 1D \mathcal{O}\mathopen{}\mathclose{{}\left(1}\right) base state. Physically, $u_{x}$ is equal to the velocity of the growing planar cathode/electrolyte interface or the dissolving planar anode/electrolyte interface in the base state. The initial conditions are given by

[TABLE]

Since $\dot{h}_{0}^{\textnormal{a}}=\dot{h}_{0}^{\textnormal{c}}=0$, h_{0}^{\textnormal{a}}\mathopen{}\mathclose{{}\left(t}\right)=0 and h_{0}^{\textnormal{c}}\mathopen{}\mathclose{{}\left(t}\right)=1 at all $t$. For galvanostatic conditions in which we apply a current density $J_{\textnormal{a}}$ on the system, we impose

[TABLE]

For potentiostatic conditions in which we apply an electric potential $V$ on the cathode, we impose $\phi_{\textnormal{e}}^{\textnormal{c}}=V$.

The equations for the time-dependent base state cannot generally be solved analytically, therefore we would have to solve them numerically. However, at steady state, the base state admits semi-analytical solutions for any $\rho_{\textnormal{s}}$ Khoo and Bazant (2018). Specifically, $c_{0}$, $\phi_{0,x}$ and their spatial derivatives can be analytically expressed in terms of the Lambert W function Corless et al. (1996). On the other hand, $\phi_{0}$ is known semi-analytically because it can be analytically expressed in terms of the Lambert W function up to an additive constant, which is a function of $J_{\textnormal{a}}$ and $\rho_{\textnormal{s}}$ and is found by numerically solving the algebraic Butler-Volmer equations given by Equations 27 and 30 with MATLAB’s $\mathtt{fsolve}$ or $\mathtt{fzero}$ function.

III.4 Perturbed state

To derive the equations for the perturbed state at \mathcal{O}\mathopen{}\mathclose{{}\left(\epsilon}\right), we substitute the perturbation expressions in Section III.1 into the full model in Section II and match terms at \mathcal{O}\mathopen{}\mathclose{{}\left(\epsilon}\right). One important outcome is that the curvature-related terms appear as functions of $k^{2}$ because they are associated with second-order spatial partial derivatives in the $y$ and $z$ directions. At \mathcal{O}\mathopen{}\mathclose{{}\left(\epsilon}\right), the governing ODEs (ordinary differential equations) are given by

[TABLE]

where the first ODE describes the perturbation in species conservation of anions and the second ODE describes the perturbation in charge conservation. For brevity, we define \alpha_{3}=-\alpha\exp\mathopen{}\mathclose{{}\left(-\alpha n\eta_{0}}\right)-\mathopen{}\mathclose{{}\left(1-\alpha}\right)\exp\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(1-\alpha}\right)n\eta_{0}}\right]. The boundary conditions at the anode at $x=h_{0}^{\textnormal{a}}$ are given by

[TABLE]

where the $\hat{D}_{1}$, $\hat{D}_{2}$ and $\hat{D}_{3}$ parameters are

[TABLE]

Because the unit normal at the cathode is in the opposite direction from that at the anode, the signs of the expressions involving $\hat{n}$ at the cathode are opposite to that at the anode. Hence, the boundary conditions at the cathode at $x=h_{0}^{\textnormal{c}}$ are given by

[TABLE]

where the $\hat{G}_{1}$, $\hat{G}_{2}$ and $\hat{G}_{3}$ parameters are

[TABLE]

The capillary number $\textnormal{Ca}=\gamma$ appears in the $\hat{D}_{1}$ and $\hat{G}_{1}$ parameters in the form of $\gamma k^{2}$, which is the source of the surface stabilizing effect that arises from the surface energy penalty incurred in creating additional surface area. The competition between this surface stabilizing effect and the surface destabilizing effect arising from the $c_{0}$, $c_{0,x}$ and $\phi_{0,x}$ fields sets the scale for the critical wavenumber $k_{\textnormal{c}}$, which is the wavenumber at which the perturbation growth rate $\omega$ is [math] and the electrode surface is marginally stable.

III.5 Discretization of perturbed state

Without making further approximations, the equations for the perturbed state do not admit analytical solutions, thus we have to resort to numerical methods to solve them. To do so, the equations for the perturbed state are spatially discretized over a uniform grid with $N$ grid points and a grid spacing $\Delta x=\frac{1}{N-1}$ using second-order accurate finite differences LeVeque (2007). Details of this discretization process are provided in Section IV of Supplementary Material. In summary, the discretized equations can be written as a generalized eigenvalue problem given by

[TABLE]

where Y,Z\in\mathbb{R}^{\mathopen{}\mathclose{{}\left(2N+2}\right)\times\mathopen{}\mathclose{{}\left(2N+2}\right)}, $v\in\mathbb{C}^{2N+2}$, $\omega\in\mathbb{C}$ and the second subscript in $c_{1,i}$ and $\phi_{1,i}$ for $i=1,2,\ldots,N$ denotes the grid point index. In the context of a generalized eigenvalue problem, the eigenvector $v$ consists of the complex-valued amplitudes $c_{1}$, $\phi_{1}$, $h_{1}^{\textnormal{a}}$ and $h_{1}^{\textnormal{c}}$ evaluated at the grid points, and the eigenvalue is the complex-valued growth rate $\omega$. Although $Y$ is non-singular, the time-independent terms in the equations for the perturbed state introduce rows of zeros in $Z$, therefore $Z$ is singular and the generalized eigenvalue problem cannot be reduced to a standard eigenvalue problem. Specifically, $Y$ is non-singular with rank $2N+2$ while $Z$ is singular with rank $N$, and the total number of eigenvalues is $2N+2$.

Because $Z$ is singular with rank $N$, there are $N$ finite eigenvalues and $N+2$ infinite eigenvalues. This mathematical property is not always consistently noted in past literature on linear stability analysis of electrodeposition, although Sundström and Bark did mention that $N$ different eigenvalues are obtained with $N$ grid points that give rise to $2N+2$ equations without explicitly stating that the other $N+2$ eigenvalues are infinite Sundström and Bark (1995). The infinite eigenvalues are physically irrelevant to the linear stability analysis Valério et al. (2007); Kawano et al. (2013), therefore we would want to focus on solving for the finite eigenvalues. This can be achieved by mapping the infinite eigenvalues to other arbitrarily chosen points in the complex plane via simple matrix transformations Goussis and Pearlstein (1989). Details of how these transformations are performed are given in Section IV of Supplementary Material. There are methods for directly removing the infinite eigenvalues such as the “reduced” method Gary and Helgason (1970); Peters and Wilkinson (1970); Valério et al. (2007) but they are more intrusive and require more extensive matrix manipulations as compared to the mapping technique Goussis and Pearlstein (1989) that we use.

The modified generalized eigenvalue problem that results from these transformations can then be solved using any eigenvalue solver. For linear stability analysis, we only need to find the eigenvalue with the largest real part instead of all the finite eigenvalues. Since the time complexity of finding all the eigenvalues typically scales as \mathcal{O}\mathopen{}\mathclose{{}\left(N^{3}}\right) while that for finding $k\leq N$ of them, where $k=1$ in our case, scales as \mathcal{O}\mathopen{}\mathclose{{}\left(kN^{2}}\right), the computational cost is dramatically reduced by a factor of \mathcal{O}\mathopen{}\mathclose{{}\left(N}\right) if we use an eigenvalue solver that can find subsets of eigenvalues and eigenvectors such as MATLAB’s $\mathtt{eigs}$ solver.

III.6 Numerical implementation

The equations for the time-dependent base state in Section III.3 are numerically solved using the finite element method in COMSOL Multiphysics 5.3a. The eigenvalue with the largest real part and its corresponding eigenvector from the generalized eigenvalue problem for the perturbed state in Section III.5 are then solved for using the $\mathtt{eigs}$ function in MATLAB R2018a. When the $\mathtt{eigs}$ function occasionally fails to converge for small values of the wavenumber $k$, we use Rostami and Xue’s eigenvalue solver based on the matrix exponential Rostami and Xue (2018); Xue , which is more robust than the $\mathtt{eigs}$ function. The colormaps used for some of the plots in Section IV are obtained from BrewerMap Cobeldick , which is a MATLAB program available in the MATLAB File Exchange that implements the ColorBrewer colormaps Harrower and Brewer (2003).

IV Results

Because of the large number of dimensionless parameters present, the parameter space is too immense to be explored thoroughly in this paper. Instead, the key dimensionless parameters that we focus on and vary are $\rho_{\textnormal{s}}$, Da and $J_{\textnormal{a}}$ under galvanostatic conditions. $\rho_{\textnormal{s}}=0$ corresponds to the classical case of an uncharged nanoporous medium while $\rho_{\textnormal{s}}\neq 0$ allows us to depart from this classical case and study its effects on the linear stability of the electrode surface. Experimentally, $\rho_{\textnormal{s}}$ can be tuned via layer-by-layer deposition of polyelectrolytes Han et al. (2014, 2016a); Hammond (2004) or tethered immobilized anions Tu et al. (2015). Da is very sensitive to the specific reactions considered and varies significantly in practice. We focus on galvanostatic conditions instead of potentiostatic conditions because when an overlimiting current $J_{\textnormal{a}}>1$ is applied on a classical system with $\rho_{\textnormal{s}}=0$, as discussed in Section II.1, the Sand’s time $t_{\textnormal{s}}$ provides a time scale at which the electric field at the cathode diverges that causes the perturbation growth rate to diverge too. This allows us to focus the linear stability analysis on times immediately before, at and immediately after $t_{\textnormal{s}}$.

For the results discussed in Sections IV.2, IV.3, IV.4 and IV.6 below, we assume the following dimensional quantities for a typical electrolyte in a typical nanoporous medium: $T=298\,\textnormal{K}$, $M_{\textnormal{m}}=6.941\,\textnormal{g}/\textnormal{mol}$ (arbitrarily pick lithium metal) Rumble (2017), $\rho_{\textnormal{m}}=0.534\,\textnormal{g}/\textnormal{cm}^{3}$ (arbitrarily pick lithium metal) Rumble (2017), $L_{x}=60\,\mu\textnormal{m}$, $L_{y}=L_{z}=100L_{x}=6\,\textnormal{mm}$ ($L_{x}\ll L_{y}=L_{z}$ to model a thin gap cell that reduces effects of gravity-induced convection (buoyancy) Elezgaray et al. (1998)), $c_{0}=10\,\textnormal{mM}$ (note that $c_{0}$ here is the dimensional initial neutral salt bulk concentration, not the dimensionless base state anion concentration), $c_{+}^{\Theta}=1\,\textnormal{M}=10^{3}\,\textnormal{mol}/\textnormal{m}^{3}$ (standard concentration) and $\gamma=1\,\textnormal{J}/\textnormal{m}^{2}$ (typical surface energy of metal/electrolyte interface) Sundström and Bark (1995). Corresponding to these dimensional quantities, all dimensionless parameters that are used for the results in Sections IV.2, IV.3, IV.4 and IV.6 are given in Table 2.

IV.1 Approximations

At the heart of the linear stability analysis is the competition between the destabilizing effect that arises from the amplification of surface protrusions by diffusive fluxes in a positive feedback loop and the stabilizing effect that arises from the surface energy penalty incurred in the creation of additional surface area. Therefore, in the dispersion relation \omega\mathopen{}\mathclose{{}\left(k}\right), we expect to see some local maxima or possibly just a single global maximum, which we denote as \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\}, where the electrode surface is maximally unstable. We also expect to see a critical wavenumber $k_{\textnormal{c}}$ corresponding to $\omega=0$ where the electrode surface is marginally stable. When $k$ is larger than $k_{\textnormal{c}}$, $\omega$ is always negative because the surface energy stabilizing effect always dominates when the wavenumber is sufficiently large. We note that $k_{\textnormal{c}}$ is always greater than $k_{\textnormal{max}}$. Corresponding to $k_{\textnormal{max}}$ and $k_{\textnormal{c}}$ are the maximum wavelength $\lambda_{\textnormal{max}}=\frac{2\pi}{k_{\textnormal{max}}}$ and critical wavelength $\lambda_{\textnormal{c}}=\frac{2\pi}{k_{\textnormal{c}}}$ respectively. In a porous medium, the characteristic pore size $h_{\textnormal{c}}=2d_{\textnormal{p}}$, where $d_{\textnormal{p}}$ is the pore diameter, sets a threshold or cutoff for overall electrode surface stabilization: we should observe stabilization if $h_{\textnormal{c}}$ is smaller than $\lambda_{\textnormal{c}}$ Tikekar et al. (2014). If $h_{\textnormal{c}}$ is larger than $\lambda_{\textnormal{c}}$, then the most unstable eigenmode dominates the electrode surface growth with a growth rate of $\omega_{\textnormal{max}}$ and the characteristic length scale of this instability is $\lambda_{\textnormal{max}}$. Therefore, \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} and $k_{\textnormal{c}}$ are the most physically informative points of the dispersion relation. We now derive an approximation for the dispersion relation \omega\mathopen{}\mathclose{{}\left(k}\right) that is valid at high values of $k$ relative to $k_{\textnormal{c}}$ and will be useful for computing \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} and $k_{\textnormal{c}}$ quickly and accurately because $k_{\textnormal{max}}$ and $k_{\textnormal{c}}$ tend to be large. The approximation is also useful for verifying the full numerical solution at high $k$, which will be discussed in Section IV.2.

When $k$ is sufficiently large compared to $k_{\textnormal{c}}$, at the cathode at $x=h_{0}^{\textnormal{c}}=1$, we expect $k^{2}c_{1}$ to balance $c_{1}^{\prime\prime}$, and $k^{2}\phi_{1}$ to balance $\phi_{1}^{\prime\prime}$ in Equations 35 and 36 respectively. Therefore, $k^{-2}$ is a small parameter multiplying the highest order spatial derivative terms $c_{1}^{\prime\prime}$ and $\phi_{1}^{\prime\prime}$, and the spatial profiles for $c_{1}$ and $\phi_{1}$ form a boundary layer with characteristic thickness $k^{-1}$. Hence, as an ansatz for the boundary layer analysis, we assume

[TABLE]

where $A$ and $B$ are arbitrary constants that are determined from the boundary conditions at $x=h_{0}^{\textnormal{c}}=1$. By assuming such an ansatz, the cathode is effectively decoupled from the anode and the perturbation growth rate is entirely dependent on the boundary conditions at the cathode. The validity of this ansatz is corroborated by our observations that $h_{1}^{\textnormal{c}}\gg h_{1}^{\textnormal{a}}$ generally in our numerical simulations, especially at large values of $k$, which was also observed by Sundström and Bark Sundström and Bark (1995). Imposing the boundary conditions at $x=h_{0}^{\textnormal{c}}=1$, we obtain

[TABLE]

where we define $\alpha_{1}=D-D_{+}$, $\alpha_{2}=z_{+}D_{+}-zD$ and $\alpha_{5}=\alpha_{2}c_{0}-z_{+}D_{+}\rho_{\textnormal{s}}$ for brevity.

Approximate values of \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} can be obtained by solving \omega^{\prime}\mathopen{}\mathclose{{}\left(k}\right)=0 and requiring \omega^{\prime\prime}\mathopen{}\mathclose{{}\left(k}\right)<0 where the primes indicate total derivatives with respect to $k$. In addition, by solving \omega\mathopen{}\mathclose{{}\left(k}\right)=0, we can obtain approximate values of $k_{\textnormal{c}}$. However, this process is tedious because the first term inside the braces in Equation 49 is a rational function that consists of polynomials in $k$ of relatively high degrees. Specifically, after multiplying the numerator and denominator of this term by $k$, it becomes a rational function with a numerator that is a polynomial in $k$ of degree $4$ and a denominator that is a polynomial in $k$ of degree $2$. Therefore, for the purpose of quickly approximating \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} and $k_{\textnormal{c}}$, we first find a simpler and yet still accurate analytical approximation for $k_{\textnormal{c}}$, which can then used as an initial guess for numerically solving for \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} using Equation 49 with MATLAB’s $\mathtt{fminbnd}$ optimizer. Such an approximation can be obtained by assuming $k_{\textnormal{c}}$ is large enough that $\hat{G}_{2}c_{1}\ll\hat{G}_{1}h_{1}$ and $\hat{G}_{3}\phi_{1}\ll\hat{G}_{1}h_{1}$ and then setting $\omega=0$ in Equation 42, thus resulting in

[TABLE]

We observe that $k_{\textnormal{c}}$ scales as $\textnormal{Ca}^{-\frac{1}{2}}=\gamma^{-\frac{1}{2}}$, which is expected because the surface energy stabilizing effect appears in the form of $\gamma k^{2}$ in $\hat{G}_{1}$ in Equation 42, and this scaling agrees with that obtained in previous work done on linear stability analysis of electrodeposition Sundström and Bark (1995); Elezgaray et al. (1998); Krishnan et al. (2002); Tikekar et al. (2014); Nielsen and Bruus (2015).

IV.2 Convergence analysis

Before analyzing the physical significance of the linear stability analysis results, we would want to first establish the accuracy and convergence of the full numerical solution of the dispersion relation \omega\mathopen{}\mathclose{{}\left(k}\right). To this end, we perform a numerical convergence analysis in which we examine the convergence of the numerical solution as the number of grid points $N$ increases. At the same time, we also compute the approximate \omega\mathopen{}\mathclose{{}\left(k}\right) given by Equation 49 because we expect the numerical and approximate solutions to agree well at high values of $k$; this therefore provides another way of checking the accuracies of both the numerical and approximate solutions.

To demonstrate how the numerical dispersion relation \omega\mathopen{}\mathclose{{}\left(k}\right) changes with $N$, we fix $\textnormal{Da}=1$ and $J_{\textnormal{a}}=1.5$ (overlimiting current) and plot numerically computed \Re\mathopen{}\mathclose{{}\left(\omega}\right) against $\log_{10}k$ for \rho_{\textnormal{s}}\in\mathopen{}\mathclose{{}\left\{-0.05,0,0.05}\right\} and N\in\mathopen{}\mathclose{{}\left\{251,501,1001,2001,4001}\right\} at specific $\frac{t}{t_{\textnormal{s}}}$ values in Figure 2. As expected, the numerical solutions converge quickly as $N$ increases from $251$ to $4001$. For $\rho_{\textnormal{s}}=0$ and $\rho_{\textnormal{s}}=0.05$ at small values of $k$, when the value of $N$ is small at $251$ or $501$, we observe that there are anomalously large values of \Re\mathopen{}\mathclose{{}\left(\omega}\right) that vanish at larger values of $N$. This is because when $N$ is too small, the grid is not sufficiently fine to accurately resolve the base state variables, in particular the rapidly increasing electric field at the cathode near $t_{\textnormal{s}}$, thus leading to an overestimation of the destabilizing effect caused by electrodiffusion and an underestimation of the stabilizing effect caused by surface energy. The numerical and approximate solutions also expectedly agree well with each other at large values of $k$ and this agreement improves as $N$ increases, thus confirming that the approximations are accurate at high $k$.

Because we are mostly interested in the $k_{\textnormal{max}}$, $\omega_{\textnormal{max}}$ and $k_{\textnormal{c}}$ points on the \omega\mathopen{}\mathclose{{}\left(k}\right) curve, we plot them against $N$ in Figure 3. We observe that the numerically computed $k_{\textnormal{max}}$, $\omega_{\textnormal{max}}$ and $k_{\textnormal{c}}$ curves rapidly level off and converge to constant values as $N$ increases. The numerical and approximate solutions also agree very well as $N$ increases, which is expected because $k_{\textnormal{max}}$ and $k_{\textnormal{c}}$ are large and the approximations are accurate at high $k$. As a compromise between numerical accuracy and computational time, we pick $N=1001$ for all numerical and approximate solutions computed in the following sections.

IV.3 Parameter sweeps

The base state anion concentration field $c_{0}$, electrolyte electric potential field $\phi_{0}$ and electric field $E_{0}=-\phi_{0,x}$ possess salient features that are useful for understanding the linear stability analysis results. We focus on galvanostatic conditions under an overlimiting current $J_{\textnormal{a}}>1$ because as explained in Section II.1, doing so provides us with Sand’s time $t_{\textnormal{s}}$ as a time scale at which the bulk electrolyte is depleted at the cathode. Depending on the sign of $\rho_{\textnormal{s}}$, the $c_{0}$, $\phi_{0}$ and $E_{0}$ fields behave differently at $t=t_{\textnormal{s}}$ and beyond. Fixing $\textnormal{Da}=1$ and $J_{\textnormal{a}}=1.5$ (overlimiting current), we plot $c_{0}$, $\phi_{0}$ and $E_{0}$ against $x$ for various $\frac{t}{t_{\textnormal{s}}}$ values for \rho_{\textnormal{s}}\in\mathopen{}\mathclose{{}\left\{-0.05,0,0.05}\right\} in Figure 4. For $\rho_{\textnormal{s}}=-0.05$, because the system can go beyond $t_{\textnormal{s}}$ and eventually reach a steady state, we show plots up to $t=2t_{\textnormal{s}}$. For $\rho_{\textnormal{s}}=0$, since $\phi_{0}$ and $E_{0}$ at the cathode diverge at $t_{\textnormal{s}}$, which cause the numerical solver to stop converging, we can only show plots up to $t=0.95t_{\textnormal{s}}$. For $\rho_{\textnormal{s}}=0.05$, because $\rho_{\textnormal{s}}>0$ effectively reduces $t_{\textnormal{s}}$ as discussed in Section II.1, we show plots up to $t=0.85t_{\textnormal{s}}$.

For $\rho_{\textnormal{s}}=-0.05<0$, the distinguishing features of running the system at an overlimiting current carried by surface conduction are the anion depletion region at the cathode and the bounded and constant electric field $E_{0}$ in this depletion region after $t=t_{\textnormal{s}}$. Because the anion concentration gradient almost vanishes in the depletion region, the current in this region is predominantly not carried by electrodiffusion but by electromigration of the counterions in the electric double layers (EDLs) under the aforementioned bounded and constant electric field $E_{0}$, i.e., surface conduction. Moreover, because of this additional surface conductivity, when compared to $\rho_{\textnormal{s}}=0$ and $\rho_{\textnormal{s}}=0.05$, $E_{0}$ is always smaller at all $x$ for a given $t$. On the other hand, for the classical case of $\rho_{\textnormal{s}}=0$, $E_{0}$ at the cathode quickly increases near $t_{\textnormal{s}}$ and eventually diverges at $t_{\textnormal{s}}$. Relative to this classical case, for $\rho_{\textnormal{s}}=0.05>0$, $E_{0}$ is always greater at all $x$ for a given $t$ and eventually diverges at the cathode earlier than $t_{\textnormal{s}}$ because of the “negative” surface conductivity conferred by the positive background charge as discussed in Section II.1.

We now examine the dispersion relation \omega\mathopen{}\mathclose{{}\left(k}\right) by plotting numerically computed \Re\mathopen{}\mathclose{{}\left(\omega}\right) against $k$ for various $\frac{t}{t_{\textnormal{s}}}$ values for \rho_{\textnormal{s}}\in\mathopen{}\mathclose{{}\left\{-0.05,0,0.05}\right\}, \textnormal{Da}\in\mathopen{}\mathclose{{}\left\{0.1,1,10}\right\} and $J_{\textnormal{a}}=1.5$ (overlimiting current) in Figure 5. In Figure 5, $\rho_{\textnormal{s}}$ increases from left to right and Da increases from bottom to top. Generally for all the parameters considered, the $\omega$ curve, in particular the $k_{\textnormal{max}}$, $\omega_{\textnormal{max}}$ and $k_{\textnormal{c}}$ points, increases and “moves in the northeast direction” as $t$ increases; qualitatively, the “total amount of instability” increases with $t$. For $\rho_{\textnormal{s}}=-0.05<0$, when compared to $\rho_{\textnormal{s}}=0$ and $\rho_{\textnormal{s}}=0.05$, the $\omega$ curve is the smallest at a given $t$ because of a smaller base state electric field $E_{0}$. The $\omega$ curve also remains bounded at all $t$ and eventually reaches a steady state that is almost attained near $t=2t_{\textnormal{s}}$ because $E_{0}$ at the cathode behaves in the same fashion. In sharp contrast, for the classical case of $\rho_{\textnormal{s}}=0$ near $t_{\textnormal{s}}$, the $\omega$ curve grows dramatically because of the rapidly increasing $E_{0}$ at the cathode, which eventually diverges at $t_{\textnormal{s}}$ and in turn causes the $\omega$ curve to diverge at $t_{\textnormal{s}}$ too. Compared to this classical case, for $\rho_{\textnormal{s}}=0.05>0$, because $E_{0}$ at the cathode is larger at a given $t$ and diverges earlier than $t_{\textnormal{s}}$, the $\omega$ curve accordingly grows even more rapidly at earlier times and diverges earlier than $t_{\textnormal{s}}$. Therefore, by bounding the electric field at the cathode, the presence of a negative background charge confers additional stabilization to the system beyond what is provided by surface energy effects, although it does not completely stabilize the system as there are still regions of positive growth rate in the dispersion relation. On the other hand, for the classical case of zero background charge, the system rapidly destabilizes near Sand’s time and ultimately diverges at Sand’s time because of the diverging electric field at the cathode, which is also demonstrated in Elezgaray et al. (1998). Relative to this classical case, the presence of a positive background charge destabilizes the system even further by generating an electric field at the cathode that is larger at a given time and diverges earlier than Sand’s time, resulting in higher growth rates at earlier times and in finite time divergence earlier than Sand’s time.

We observe that increasing Da generally increases $\omega$ but this effect is very insignificant because the application of an overlimiting current implies that the system is always diffusion-limited regardless of what Da is. Hence, in this regime of diffusion-limited electrodeposition under an overlimiting current, specific details of the electrochemical reaction kinetics model are not important in influencing the dispersion relation as long as the model includes the surface energy stabilizing effect, which typically occurs in the functional form of $\gamma k^{2}$.

In the interest of space, plots of numerically computed \Re\mathopen{}\mathclose{{}\left(\omega}\right) against $k$ for $J_{\textnormal{a}}=1$ (limiting current) and $J_{\textnormal{a}}=0.5$ (underlimiting current) are not shown here but are given in Figures 1 and 2 in Section V of Supplementary Material respectively. Since the system is still always diffusion-limited for $J_{\textnormal{a}}=1$, the trends observed for $J_{\textnormal{a}}=1$ are qualitatively similar to our previous discussion for $J_{\textnormal{a}}=1.5$, except that the $\omega$ values are smaller because a smaller applied current density results in a smaller electric field at the cathode. For $J_{\textnormal{a}}=0.5$, because the applied current density is underlimiting, Sand’s time is not defined and at the cathode, the bulk electrolyte concentration does not vanish and the electric field does not diverge at any $t$. Therefore, the $\omega$ curve remains bounded at all $t$ and reaches a steady state eventually. Moreover, $\omega$ generally increases with Da, and this increase is especially pronounced when Da increases from $1$ to $10$; this effect was also observed by Sundström and Bark Sundström and Bark (1995) who focused their analysis on underlimiting currents. This increase in $\omega$ is not directly caused by $E_{0}$ because $E_{0}$ does not change appreciably despite the increase in Da (refer to Figures 3 to 5 in Section VI of Supplementary Material). Rather, as discussed in Section II.2, the system becomes diffusion-limited when $\textnormal{Da}\gg 1$, causing the surface perturbations to destabilize faster.

As discussed in Section IV.1, at each $t$ point, each $\omega$ curve exhibits a global maximum \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} and a critical wavenumber $k_{\textnormal{c}}$, which is where the curve crosses the horizontal axis $\omega=0$. The \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} and $k_{\textnormal{c}}$ points provide a succinct way to summarize the most physically significant features of the \omega\mathopen{}\mathclose{{}\left(k}\right) curve for all the parameter ranges we have explored thus far. Therefore, for \rho_{\textnormal{s}}\in\mathopen{}\mathclose{{}\left\{-0.05,0,0.05}\right\}, \textnormal{Da}\in\mathopen{}\mathclose{{}\left\{0.1,1,10}\right\} and J_{\textnormal{a}}\in\mathopen{}\mathclose{{}\left\{0.5,1,1.5}\right\}, we plot numerically computed $k_{\textnormal{max}}$ and $\omega_{\textnormal{max}}$ against $\frac{t}{t_{\textnormal{s}}}$ in Figure 6 and numerically computed $k_{\textnormal{c}}$ against $\frac{t}{t_{\textnormal{s}}}$ in Figure 7. For $J_{\textnormal{a}}\geq 1$, we observe that the $k_{\textnormal{max}}$ and $\omega_{\textnormal{max}}$ curves diverge near $t_{\textnormal{s}}$ for $\rho_{\textnormal{s}}\geq 0$ but level off to constant values past $t_{\textnormal{s}}$ for $\rho_{\textnormal{s}}<0$, therefore these curves appear as if they are “fanning out”. In contrast, for $J_{\textnormal{a}}<1$, the $k_{\textnormal{max}}$ and $\omega_{\textnormal{max}}$ curves level off past $t_{\textnormal{s}}$ for all values of $\rho_{\textnormal{s}}$ as the system eventually reaches a steady state when an underlimiting current is applied. The $k_{\textnormal{c}}$ curves have the same qualitative shape as the $k_{\textnormal{max}}$ curves except that they are larger, as expected. The effects of Da and $J_{\textnormal{a}}$ on the $k_{\textnormal{max}}$, $\omega_{\textnormal{max}}$ and $k_{\textnormal{c}}$ values, which are previously discussed in the context of the dispersion relation, are also clearly reflected in Figures 6 and 7.

In an effort to make the electrode surface less unstable at overlimiting current, we focus on $\rho_{\textnormal{s}}<0$ to determine how much additional stabilization a negative $\rho_{\textnormal{s}}$ confers to the surface as it gets increasingly more negative. Subsequently, we plot numerically computed $k_{\textnormal{max}}$, $\omega_{\textnormal{max}}$ and $k_{\textnormal{c}}$ against $\frac{t}{t_{\textnormal{s}}}$ for \rho_{\textnormal{s}}\in\mathopen{}\mathclose{{}\left\{-1,-0.75,-0.5,-0.25,-0.05}\right\}, $\textnormal{Da}=1$ and $J_{\textnormal{a}}=1.5$ in Figure 8. While a more negative $\rho_{\textnormal{s}}$ generally decreases $k_{\textnormal{max}}$, $\omega_{\textnormal{max}}$ and $k_{\textnormal{c}}$, it is clear that there are diminishing returns to the amount of additional stabilization achieved. It also appears that complete stabilization is not possible as $\omega_{\textnormal{max}}$ remains positive even for $\rho_{\textnormal{s}}=-1$, albeit at a small value. In practice, it is probable that a sufficiently small and positive $\omega_{\textnormal{max}}$ value can be deemed to be small enough for considering an electrode surface “practically stable”, but experiments that measure and correlate $\omega_{\textnormal{max}}$ with observations of metal growth need to be performed in order to determine this threshold $\omega_{\textnormal{max}}$ value.

IV.4 Comparison between numerical and approximate solutions

To illustrate how well the approximations given by Equations 49 and 50 work for the parameter ranges considered, we plot numerical and approximate values of $k_{\textnormal{max}}$, $\omega_{\textnormal{max}}$ and $k_{\textnormal{c}}$ against $\frac{t}{t_{\textnormal{s}}}$ for \rho_{\textnormal{s}}\in\mathopen{}\mathclose{{}\left\{-0.05,0,0.05}\right\}, $\textnormal{Da}=1$ and $J_{\textnormal{a}}=1.5$ in Figure 9. In the interest of space, these plots for other values of Da and $J_{\textnormal{a}}$ are provided in Figures 6 to 11 of Section VII of Supplementary Material. For all parameter ranges considered, the agreement between numerical and approximate values of $k_{\textnormal{max}}$, $\omega_{\textnormal{max}}$ and $k_{\textnormal{c}}$ is excellent, giving us confidence that the approximations are useful for rapidly and accurately computing $k_{\textnormal{max}}$, $\omega_{\textnormal{max}}$ and $k_{\textnormal{c}}$. This confirms that $k_{\textnormal{max}}$ and $k_{\textnormal{c}}$ are large enough that Equations 49 and 50, which have assumed that $k$ is sufficiently large, are accurate for approximating them. We will therefore use Equations 49 and 50 extensively in Sections IV.5 and IV.6 that follow.

IV.5 Application to copper electrodeposition

We now apply linear stability analysis to the specific case of copper electrodeposition and electrodissolution and compare it with experimental data Han et al. (2016a) to determine how well theory agrees with experiment. Because copper electrodeposition involves the overall transfer of two electrons that are transferred one at a time in a serial manner, we need to first derive the overall expression for the Faradaic current density $J_{\textnormal{F}}$.

Assuming that the activity of electrons is $1$ and dilute solution theory is applicable, for a $n$-electron transfer reaction, the dimensionless forms of Equations 12 and 9 are given by

[TABLE]

For multistep electron transfer reactions, it is more convenient to work with $\Delta\phi$ instead of $\eta$. Therefore, we rewrite $J_{\textnormal{F}}$ in terms of $\Delta\phi$ as

[TABLE]

where $k_{\textnormal{c}}$ and $k_{\textnormal{a}}$ are the cathodic and anodic rate constants respectively.

The reaction mechanism for copper electrodeposition and electrodissolution is given by Newman and Thomas-Alyea (2004); Mattsson and Bockris (1959); Bockris and Enyo (1962); Brown and Thirsk (1965)

[TABLE]

where (aq), (ads) and (s) indicate aqueous, adsorbed and solid respectively. The first step is assumed to be the rate-determining step while the second step is assumed to be at equilibrium. Applying Equation 53 to each step, noting that the activity of solid metal is $1$ and rewriting $J_{\textnormal{F}}$ in terms of $\eta$, we obtain

[TABLE]

where $\alpha_{1}$ is the charge transfer coefficient of the first step.

Previously in Section II.2, for a $1$-step $n$-electron transfer metal electrodeposition reaction, the dimensionless forms of Equations 12 and 13 are given by

[TABLE]

By comparing Equations 57 and 58 with Equations 59 and 60, we set $n=2$ and $\alpha=\frac{\alpha_{1}}{2}$ and replace $\gamma$ with $2\gamma$ in the original set of equations in order to adapt the linear stability analysis for copper electrodeposition.

By carrying out nonlinear least squares fitting on experimental steady state current-voltage relations, we have previously performed parameter estimation Khoo and Bazant (2018) for copper electrodeposition in a copper(II) sulfate ($\textnormal{CuSO}_{4}$) electrolyte in cellulose nitrate (CN) membranes Han et al. (2016a), which are a random nanoporous medium with well connected pores. The parameters that are estimated are $\rho_{\textnormal{s}}$, $\tau$, Da, $\alpha_{1}$ and $\epsilon_{\textnormal{p}}$ and their fitted values are provided in Table III in Khoo and Bazant (2018). Other parameters specific to the copper electrodeposition reaction, $\textnormal{CuSO}_{4}$ electrolyte and CN membranes used are also provided in Tables I and II in Khoo and Bazant (2018). For the surface energy of the copper/electrolyte interface, we use dimensional $\gamma=1.85\,\textnormal{J}/\textnormal{m}^{2}$ given in Table I in Nielsen and Bruus (2015).

For our analysis here, the specific experimental datasets that we focus on are labeled $\textnormal{CN}_{2}(-)$ and $\textnormal{CN}_{2}(+)$ in Khoo and Bazant (2018), which correspond to negatively and positively charged CN membranes respectively with a dimensional electrolyte concentration $c_{0}$ of $100\,\textnormal{mM}$. We will drop the $2$ subscript for brevity. The morphologies of the electrodeposited copper films, which are visualized by EDS (energy dispersive X-ray spectroscopy) maps, for these $\textnormal{CN}(-)$ and $\textnormal{CN}(+)$ membranes at $2000\,\textnormal{s}$ for dimensional applied currents $I_{\textnormal{a}}$ of $15\,\textnormal{mA}$, $20\,\textnormal{mA}$ and $25\,\textnormal{mA}$ are given in Figure 10(a) that consists of magnifications of EDS maps taken from Figures 6(a) to 6(f) of Han et al. (2016a). At $15\,\textnormal{mA}$, the copper films for both $\textnormal{CN}(+)$ and $\textnormal{CN}(-)$ membranes appear to be uniform and stable. However, at $20\,\textnormal{mA}$ and $25\,\textnormal{mA}$, the film for $\textnormal{CN}(+)$ becomes very unstable and roughens more as the applied current increases. It is difficult to determine quantitatively the instability wavelength using the relatively low resolution EDS maps but it is probably much smaller than $5\,\mu\textnormal{m}$. In contrast, for $\textnormal{CN}(-)$, the film still remains uniform and stable at $20\,\textnormal{mA}$ but slightly destabilizes and roughens at $25\,\textnormal{mA}$ with an instability wavelength probably on the order of $5\,\mu\textnormal{m}$. In summary, the onset of overall electrode surface destabilization occurs at $20\,\textnormal{mA}$ for $\textnormal{CN}(+)$ with an instability wavelength of much smaller than $5\,\mu\textnormal{m}$ and at $25\,\textnormal{mA}$ for $\textnormal{CN}(-)$ with an instability wavelength of about $5\,\mu\textnormal{m}$.

Because the EDS maps are taken at $2000\,\textnormal{s}$, which is much longer than the diffusion times for $\textnormal{CN}(-)$ and $\textnormal{CN}(+)$ of $41.8\,\textnormal{s}$ and $40.9\,\textnormal{s}$ respectively, we assume that the system is at steady state. This assumption allows us to use the semi-analytical expressions for the base state variables in Khoo and Bazant (2018), which we have previously discussed in Section III.3, to compute approximate values of \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} and $k_{\textnormal{c}}$ using Equations 49 and 50. The $\textnormal{CN}(-)$ dataset has a dimensional limiting current of $18.2\,\textnormal{mA}$ while the $\textnormal{CN}(+)$ dataset has a dimensional maximum current, which we have discussed in Section II.1, of $16.9\,\textnormal{mA}$. Therefore, for $\textnormal{CN}(-)$, the three applied currents of $15\,\textnormal{mA}$, $20\,\textnormal{mA}$ and $25\,\textnormal{mA}$ correspond to underlimiting, slightly overlimiting and overlimiting currents respectively. On the other hand, for $\textnormal{CN}(+)$, the model eventually diverges and does not admit a steady state when the applied current $I_{\textnormal{a}}$ is above the maximum current $I_{\textnormal{max}}$, therefore we can only obtain finite values of \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} and $k_{\textnormal{c}}$ for the applied current of $15\,\textnormal{mA}$ while the model predicts infinite values of \mathopen{}\mathclose{{}\left\{k_{\textnormal{max}},\omega_{\textnormal{max}}}\right\} and $k_{\textnormal{c}}$ for the applied currents of $20\,\textnormal{mA}$ and $25\,\textnormal{mA}$ due to finite time divergence of the system. Other fitted key dimensionless parameters include $\rho_{\textnormal{s}}\approx-0.01$ and $\textnormal{Da}\approx 2.50$ for $\textnormal{CN}(-)$ and $\rho_{\textnormal{s}}\approx 0.236$ and $\textnormal{Da}\approx 0.473$ for $\textnormal{CN}(+)$.

To summarize the model predictions, we plot approximate dimensional values of $\lambda_{\textnormal{c}}$ and $\lambda_{\textnormal{max}}$ against the dimensional applied current $I_{\textnormal{a}}$ in Figure 10. In the $\lambda_{\textnormal{c}}$ plot in Figure 10(a), we also indicate the characteristic pore size $h_{\textnormal{c}}$ of $0.5\pm 0.1\,\mu\textnormal{m}$, which is given by twice the pore diameter $d_{\textnormal{p}}$ of $250\pm 50\,\textnormal{nm}$ Han et al. (2016a), in order to determine if the model predicts overall electrode surface stabilization. As discussed in Section IV.1, we expect overall electrode surface stabilization if $h_{\textnormal{c}}<\lambda_{\textnormal{c}}$, which corresponds to the blue shaded region in the $\lambda_{\textnormal{c}}$ plot. On the contrary, we expect overall electrode surface destabilization if $h_{\textnormal{c}}>\lambda_{\textnormal{c}}$, which corresponds to the red shaded region in the $\lambda_{\textnormal{c}}$ plot, and the characteristic instability wavelength is $\lambda_{\textnormal{max}}$. Comparing the $\lambda_{\textnormal{c}}$ plot with our previous discussion of the onset of overall electrode surface destabilization suggested by the copper film morphologies observed experimentally, we see that the model generally agrees well with experiment; the only disagreement is at $I_{\textnormal{a}}=20\,\textnormal{mA}$ where the model predicts destabilization for $\textnormal{CN}(-)$, which has $\lambda_{\textnormal{c}}=0.352\,\mu\textnormal{m}$, while the EDS map of the copper film at this applied current shows that the film appears to be stable. Nonetheless, this disagreement in theory and experiment is relatively minor since $\lambda_{\textnormal{c}}=0.352\,\mu\textnormal{m}$ for $\textnormal{CN}(-)$ at $I_{\textnormal{a}}=20\,\textnormal{mA}$ is only slightly smaller than the mean of $h_{\textnormal{c}}$ of $0.5\,\mu\textnormal{m}$ and is almost equal to the lower bound of $h_{\textnormal{c}}$ of $0.4\,\mu\textnormal{m}$. In addition, the model predicts $\lambda_{\textnormal{max}}=0$ (because $k_{\textnormal{max}}\rightarrow\infty$ when $I_{\textnormal{a}}>I_{\textnormal{max}}$) at $I_{\textnormal{a}}=20\,\textnormal{mA}$ and $I_{\textnormal{a}}=25\,\textnormal{mA}$ for $\textnormal{CN}(+)$ while it predicts $\lambda_{\textnormal{max}}=1.09\,\mu\textnormal{m}$ at $I_{\textnormal{a}}=25\,\textnormal{mA}$ for $\textnormal{CN}(-)$. These model predictions of $\lambda_{\textnormal{max}}$ qualitatively agree well with the experimentally observed instability wavelengths at these applied currents that we have previously discussed. Therefore, in conclusion, the theory agrees reasonably well with experimental data, especially given that many assumptions and simplifications are made in the model.

IV.6 Pulse electroplating and pulse charging

For many electrochemical applications such as electroplating and charging of metal batteries, which is equivalent to electrodeposition at the metal negative electrode, it is desirable to operate them as quickly as possible at a high current without causing the formation of dendrites that short-circuit the system. To delay or prevent the formation of dendrites, it is common to perform pulse electroplating of metals Devaraj et al. (1990); Chandrasekar and Pushpavanam (2008) or pulse charging of lithium metal batteries (LMBs) and lithium-ion batteries (LIBs) Li et al. (2001); Purushothaman et al. (2005); Purushothaman and Landau (2006); Zhang (2006); Hussein and Batarseh (2011); Shen et al. (2012); Savoye et al. (2012); Mayers et al. (2012); Aryanfar et al. (2014) so that there is sufficient time between pulses for the concentration gradients and electric field in the system to relax. For pulse electroplating of metals, it has been empirically observed that the crystal grain size generally decreases with applied current Devaraj et al. (1990); Chandrasekar and Pushpavanam (2008). Using an applied direct current to perform silver electrodeposition under galvanostatic conditions, Aogaki experimentally observed that the crystal grain size decreases with time Aogaki (1982a, b), which agrees well with theoretical predictions from linear stability analysis previously done by Aogaki and Makino Aogaki and Makino (1981). With all these considerations in mind, we apply our linear stability analysis with a time-dependent base state as a tool to investigate how pulse electroplating protocols with high average applied currents, which are inherently time-dependent, affect the linear stability of the electrode surface and the crystal grain size for both zero and negative pore surface charges.

Based on the results in Section IV.5, we generally expect the characteristic pore size $h_{\textnormal{c}}$ to be larger than the critical wavelength $\lambda_{\textnormal{c}}$ at high applied currents, therefore the electrode surface is unstable with a characteristic instability wavelength $\lambda_{\textnormal{max}}$. Because a pulse current is applied, $\lambda_{\textnormal{max}}$ varies in time and hence, it would be useful to define an average $\lambda_{\textnormal{max}}$ that averages out the effect of time. In this spirit, we define the average maximum wavenumber $\bar{k}_{\textnormal{max}}$ and the corresponding average maximum wavelength $\bar{\lambda}_{\textnormal{max}}$ as

[TABLE]

where $t_{\textnormal{f}}$ is the final time of the pulse and each maximum wavenumber $k_{\textnormal{max}}$ is weighted by its corresponding maximum growth rate $\omega_{\textnormal{max}}$. We expect $\bar{\lambda}_{\textnormal{max}}$ to be on the same order of magnitude as the the crystal grain size that is observed experimentally. As a simple example, we suppose that the pulse electroplating protocol is a periodic pulse wave $J_{\textnormal{a}}$ with an “on” (charging) time of $\Delta t_{\textnormal{on}}$, a “off” (relaxation) time of $\Delta t_{\textnormal{off}}$, and a period $T$ given by $T=\Delta t_{\textnormal{on}}+\Delta t_{\textnormal{off}}$. The duty cycle $\gamma_{\textnormal{dc}}$ is given by $\gamma_{\textnormal{dc}}=\frac{\Delta t_{\textnormal{on}}}{T}$ and the average applied current density $\bar{J}_{\textnormal{a}}$ over one period is given by $\bar{J}_{\textnormal{a}}=J_{\textnormal{a,p}}\gamma_{\textnormal{dc}}$ where $J_{\textnormal{a,p}}$ is the peak applied current density. Hence, for a particular $\bar{J}_{\textnormal{a}}$, a smaller $\gamma_{\textnormal{dc}}$ implies a larger $J_{\textnormal{a,p}}$.

For the classical case of $\rho_{\textnormal{s}}=0$, we pick $\bar{J}_{\textnormal{a}}=1$ and $\Delta t_{\textnormal{on}}=0.0125t_{\textnormal{s}}$ and vary $\gamma_{\textnormal{dc}}$ from $0.2$ to $1$ (direct current) where the Sand’s time $t_{\textnormal{s}}$ is calculated based on $\bar{J}_{\textnormal{a}}$. $\bar{J}_{\textnormal{a}}$, $\Delta t_{\textnormal{on}}$ and $\gamma_{\textnormal{dc}}$ should be carefully chosen such that $J_{\textnormal{a,p}}$ is not too high to deplete the bulk electrolyte at the cathode during the “on” cycle so that the system does not diverge at any point in time; this explains why $\gamma_{\textnormal{dc}}<0.2$ for our choice of $\bar{J}_{\textnormal{a}}=1$ and $\Delta t_{\textnormal{on}}=0.0125t_{\textnormal{s}}$ cannot be numerically simulated. For $\rho_{\textnormal{s}}=-0.05$, we pick $\bar{J}_{\textnormal{a}}=1.5$ and $\Delta t_{\textnormal{on}}=t_{\textnormal{s}}$ and vary $\gamma_{\textnormal{dc}}$ from $0.1$ to $1$ (direct current) to drive the system at an overlimiting average applied current density. We also fix $\textnormal{Da}=1$ for both cases and use Equations 49 and 50 to compute approximate values of $k_{\textnormal{max}}$ and $\omega_{\textnormal{max}}$. For these choices of parameters, as an illustrative example, we plot $J_{\textnormal{a}}$, approximate $k_{\textnormal{max}}$ and approximate $\omega_{\textnormal{max}}$ against $t$ for $\gamma_{\textnormal{dc}}=0.5$ in Figure 11. We note that the large overshoot in $k_{\textnormal{max}}$ at the beginning of each “on” cycle for $\rho_{\textnormal{s}}=0$ is caused by the sharp rate of increase of the concentration gradients and electric field as $J_{\textnormal{a}}$ rapidly increases from [math] in the “off” cycle to $J_{\textnormal{a,p}}$ in the “on” cycle. Corresponding to these pulse waves, we plot $\bar{\lambda}_{\textnormal{max}}$ against $\gamma_{\textnormal{dc}}$ in Figure 12. For both $\rho_{\textnormal{s}}=0$ and $\rho_{\textnormal{s}}=-0.05$, $\bar{\lambda}_{\textnormal{max}}$ increases with $\gamma_{\textnormal{dc}}$, which agrees with the empirical observation that the crystal grain size generally decreases with applied current Devaraj et al. (1990); Chandrasekar and Pushpavanam (2008). The ability to experimentally impart a negative pore surface charge to the nanoporous medium therefore enables pulse electroplating at overlimiting currents for electrodepositing a large amount of charge at a high rate and tuning the desired crystal grain size.

V Conclusion

We have derived the full model that couples the leaky membrane model for ion transport, which is capable of predicting overlimiting current due to surface conduction, with Butler-Volmer reaction kinetics, which describes the metal electrodeposition reaction, and performed linear stability analysis on it with respect to a time-dependent base state. The volume-averaged background charge density can generally be of any sign. As a result, we have generalized previous work on linear stability analysis of electrodeposition carried out in Sundström and Bark (1995); Elezgaray et al. (1998); Tikekar et al. (2014). We then performed a boundary layer analysis on the perturbed state in order to derive an accurate approximation for the dispersion relation and a convergence analysis to verify the accuracy and convergence of the full numerical solution of the dispersion relation. By performing parameter sweeps over the volume-averaged background charge density, Damköhler number and applied current density under galvanostatic conditions, we have concluded that a negative background charge significantly stabilizes the electrode surface instability, although it does not completely stabilize it, while a positive background charge further destabilizes this instability. We have also verified that the approximations for the maximum wavenumber, maximum growth rate and and critical wavenumber are very accurate, and applied them to demonstrate good agreement between theory and experimental data for copper electrodeposition in cellulose nitrate membranes Han et al. (2016a). Lastly, we have employed the linear stability analysis as a tool to analyze the dependence of the crystal grain size on duty cycle in pulse electroplating. These results demonstrate the predictive power and robustness of the theory despite its simplicity. Although detailed analysis of the Poisson-Nernst-Planck-Stokes equations for transport in a microchannel by Nielsen and Bruus Nielsen and Bruus (2014) reveals that the leaky membrane model for surface conduction is at best a rough approximation of the real system, the good agreement between theory and experiment that we have demonstrated suggests that the model is applicable in similar electrochemical systems using charged membranes such as shock electrodeposition for information storage applications Han et al. (2016b) and shock electrodialysis for water treatment Deng et al. (2013); Schlumpberger et al. (2015); Schlumpberger (2016).

We have made many assumptions and simplifications in the model presented, and relaxing some of them offers opportunities for extending it in useful ways. First, we have ignored surface adsorption, surface diffusion of adsorbed species Aogaki and Makino (1984a, b, c, d) and additional mechanical effects such as pressure, viscous stress and deformational stress Monroe and Newman (2004, 2005); Tikekar et al. (2016b); Ahmad and Viswanathan (2017a, b); Natsiavas et al. (2016), which confer additional stabilization to the electrode surface. Adding these physics and chemistry to the model are likely to result in finite values of the maximum wavenumber, maximum growth rate and critical wavenumber near and at Sand’s time under an overlimiting current for zero and positive background charges respectively, as opposed to diverging in our current model. The inclusion of these additional mechanical effects will also extend the applicability of the model to solid electrolytes Liu et al. (2017) that are used in solid state batteries. Deformation of the porous medium caused by metal growth inside the pores also results in a porosity that varies in both time and space whose effects would be interesting to study. Second, in order to apply the linear stability analysis to lithium metal batteries (LMBs), we would also need to model the solid electrolyte interphase (SEI) layer Verma et al. (2010); Cheng et al. (2016); Peled and Menkin (2017); Wang et al. (2018), which will certainly increase the complexity of the model but also make it more predictive. Incorporating these two aforementioned extensions into the model may help explain recent experimental studies of lithium growth that have demonstrated that competing SEI reactions and stress effects lead to root growth before Sand’s time or below limiting current Bai et al. (2016); Kushima et al. (2017); Bai et al. (2018), which is different from tip growth of dendrites under transport limitation that we have focused on in this paper. Third, other chemical mechanisms for overlimiting current such as water splitting Nikonenko et al. (2010, 2014) and current-induced membrane discharge Andersen et al. (2012) may be present. These effects are typically highly nonlinear and therefore, we expect them to significantly influence the transient base state and linear stability analysis. Fourth, we should consider the effects of coupling nucleation, which is fundamentally a nonlinear instability unlike spinodal decomposition that is a linear instability, to the current model. Specifically, nucleation may affect the transient base state during initial and early reaction-limited surface growth and create surface roughness on the scale of the characteristic nucleus size, which may in turn influence overall electrode surface stabilization or destabilization when the system reaches transport limitation near or at Sand’s time. Fifth, an interesting and useful generalization of the reaction model would be to use the symmetric Marcus-Hush-Chidsey kinetics Bai and Bazant (2014); Zeng et al. (2014) or asymmetric Marcus-Hush kinetics Zeng et al. (2015) instead of Butler-Volmer kinetics for modeling electron transfer reactions, which would afford us the reorganization energy as a key system parameter whose influence on the linear stability of the electrode surface can be investigated.

Acknowledgements.

E. Khoo acknowledges support from the National Science Scholarship (PhD) funded by Agency for Science, Technology and Research, Singapore (A*STAR). H. Zhao and M. Z. Bazant acknowledge support from the Toyota Research Institute through the D3BATT Center on Data-Driven-Design of Rechargeable Batteries.