Linear stability analysis of transient electrodeposition in charged porous media

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. In contrast, 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 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 pulse electroplating to analyze how the crystal grain size varies with duty cycle.


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 [1][2][3]. 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 [4][5][6][7] 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 [8]. 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 [2,[9][10][11][12][13][14][15], solidification [1][2][3]9] and morphogenesis [3,16]. Some particular examples of morphological stability analysis include the Saffman-Taylor instability (viscous fingering) [17][18][19][20], viscous fingering coupled with electrokinetic effects [21], the Mullins-Sekerka instability of a spherical particle during diffusion-controlled or thermally controlled growth [22] and of a planar interface during solidification of a dilute binary alloy [23,24], and control of phase separation using electro-autocatalysis or electro-autoinhibition in driven open electrochemical systems [25,26].

A. 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.
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 [103,104]. In practice, overlimiting current (OLC) beyond the electrodiffusion limit has been observed experimentally in ion-exchange membranes [105][106][107][108][109][110][111][112][113][114][115][116] and microchannels and nanochannels [117][118][119][120][121][122][123][124]. Depending on the length scale of the pores or channel, some possible physical mechanisms for OLC [125] are surface conduction [119][120][121][126][127][128], electroosmotic flow [129,130] and electroosmotic instability [8,131]. Some chemical mechanisms for OLC include water splitting [114,115] and current-induced membrane discharge [132].
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 [125]. 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 [127,128,133] and microchannels and nanochannels [119-121, 125, 126, 134]. 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.

B. Theories of pattern formation
Morphological stability analysis of electrodeposition is typically performed in the style of the pioneering Mullins-Sekerka stability analysis [22,23]. 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 [22]. 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 [23]. 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 [135].
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 [136][137][138]. In 1984, Aogaki and Makino extended their previous work to account for surface diffusion of adsorbed metal atoms under galvanostatic [139,140] and potentiostatic conditions [141,142]. 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 [143] and applied it to study how hydrogen adsorption affects these extracted parameters under galvanostatic conditions [144]. Later work by Barkey, Muller and Tobias in 1989 [145,146], and Chen and Jorne in 1991 [147] 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 [148]. 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 [149]. 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 [150]. The role of electrolyte additives in stabilizing electrodeposition was examined in the linear stability analysis performed by Bocarsly in 2002 and2003 [151-153], and McFadden et al. in 2003 [154]. 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 [155]. 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 [156,157]. 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 [158]. Tikekar, Archer and Koch then extended this work in 2016 by accounting for elastic deformations that provide further stabilization [159]. Building on Monroe and Newman's 2004 and 2005 work on interfacial deformation effects [156,157], 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 [160], and further generalized this work in the same year to account for anisotropy [161]. Natsiavas, Weinberg, Rosato and Ortiz in 2016 also investigated the stabilizing effect of prestress and showed good agreement between theory and experiment [162].
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 [163].
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 [164]; some experimental results were also obtained by Fleury, Chazalviel, Rosso and Sapoval in support of this model [165], and some of the numerical results of the original analysis were subsequently improved by Rosso, Chazalviel and Chassaing [166]. Via an asymptotic anal-ysis 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 [167]. Building on past theoretical and experimental work done on silver electrodeposition by Barton and Bockris [168], and zinc electrodeposition by Despic, Diggle and Bockris [169,170], Monroe and Newman investigated the propagation velocity and length of a dendrite tip that grows via Butler-Volmer kinetics [171]. 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 [172]. 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 [173] and how they are affected by temperature [174].

C. 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 [149,150,158]. For simplicity, we ignore bulk convection, electroosmotic flow, surface adsorption, surface diffusion of adsorbed species [139][140][141][142] and additional mechanical effects such as pressure, viscous stress and deformational stress [156,157,[159][160][161][162]. 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 [175][176][177][178] 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 [179], and also use the stability analysis as a tool for investigating the dependence of crystal grain size on duty cycle during pulse electroplating.

A. 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 [180]. The system under consideration is a binary asymmetric liquid electrolyte in a finite 3D charged random nanoporous medium where x ∈ [0, L x ], y ∈ [0, L y ] and z ∈ [0, L z ], 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 [179]. 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 [149,150,163], we pick a moving reference frame with a velocity u(t) = u x (t)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 p , τ , a p and σ s respectively. where ρ s is the volume-averaged background charge density.
We assume that there are no homogeneous reactions and all material properties such as p , τ , a p and σ 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 ±0 , where the 0 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 ±0 is related to the molecular (free solution) diffusivity D m ±0 by [181]. 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(t) = u x (t)e x is only significant in the equation describing electrode surface growth or dissolution [149,150,163], which we will discuss in Section II C. Under these assumptions, for ion transport, the Nernst-Planck equations describing species conservation, charge conservation equation and macroscopic electroneutrality constraint are given by where c ± , F ± , z ± , are the cation and anion concentrations, fluxes and charge numbers respectively, φ is the electrolyte electric potential, J is the electrolyte current density, h p = p ap is the effective pore size and ρ 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 ν ± , electroneutrality requires that z + ν + + z − ν − = 0. We will use the macroscopic electroneutrality constraint given by Equation 3 to eliminate c + as a dependent variable, therefore leaving c − and φ as the remaining dependent variables. For a classical system with an uncharged nanoporous medium, i.e., ρ s = 0, the maximum current density that the system can possibly attain under electrodiffusion is called the diffusion-limited or limiting current density J lim , which is given by [180] where c 0 is the neutral salt bulk concentration. The limiting current I lim is then given by I lim = J lim A where A = L y L z is the surface area of the electrode. Under galvanostatic conditions, when a current density J a larger than J lim is applied, the cation and anion concentrations at the cathode reach 0 and the electrolyte electric potential and electric field there diverge in finite time, which is called Sand's time [182] denoted as t s ; see [183] for a discussion of some subtlety associated with this transition time when J a is exactly equal to J lim . Defining the dimensionless Sand's timet s = D amb0 ts L 2 x and dimensionless applied current [103,180] is the ambipolar diffusivity of the neutral salt in the dilute limit and L 2 x D amb0 is the diffusion time scale,t s is given by [183] t s = π 16J 2 a ,J a > 1.
For galvanostatic conditions, t s is a critically important time scale because the formation of dendrites often occurs near or at t s , therefore it will be central to the linear stability analysis results discussed in Section IV.
Unlike the classical case of ρ s = 0, for ρ s < 0, the system can sustain an overlimiting currentJ 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 s and eventually reach a steady state, in stark contrast to the finite time divergence of the classical case at t s . On the other hand, for ρ 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 ρ s > 0 relative to ρ s = 0, at the cathode, the bulk electrolyte concentration vanishes and the electric field diverges earlier than t s ; in other words, ρ s > 0 effectively reduces t 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 max , to be smaller than J lim , and J max decreases as ρ s increases. In effect, a more positive ρ s decreases J lim , which thus leads to a smaller t s for a given J a according to Equation 5. Details of how to numerically compute J max are found in [180]; note that J max here corresponds to I BV max in [180].

B. 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 H is given by H = − 1 2 ∇ s ·n where ∇ s is the surface gradient operator andn is the unit normal that points outwards from the electrolyte [184]. 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 electron e with charge −1, R z R is the reduced species R with charge z R , and z O − n = z R because of charge conservation. If the reduced species is solid metal, i.e., z 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 µ i of the oxidized species O, electron e and reduced species R for i ∈ {O, e, R} are given by where the surface energy term 2ΩγH [149, 150, 156-158, 163, 171, 184] is included in µ R when the reduced species is solid metal (z R = 0) and the Θ superscript indicates standard state. The activity of species i is given by a i = γ iĉi where γ i is the activity coefficient of is the concentration of species i normalized by its standard concentration c Θ i . µ Θ i is the standard electrochemical potential of species i, φ e is the electrode electric potential, Ω = Mm ρm is the atomic volume of the solid metal where M m and ρ m are the atomic mass and mass density of the metal respectively, and γ is the isotropic surface energy of the metal/electrolyte interface. The quantity Ωγ k B T is the capillary constant that has units of length [22][23][24]. The interfacial electric potential difference ∆φ is defined as ∆φ = φ e − φ.
At equilibrium when µ O + nµ e = µ R , we obtain the Nernst equation where the "eq" superscript denotes equilibrium and E Θ is the standard electrode potential.
When the system is driven out of equilibrium so that µ O + nµ e = µ R , it generates a Faradaic current density J F that is given by [25,26,181] where k 0 is the overall reaction rate constant and µ r,ex ‡ is the excess electrochemical potential of the transition state for the Faradaic reaction. Using the Butler-Volmer hypothesis, consists of a chemical contribution k B T ln γ r ‡ , where γ 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 α, which is the charge transfer coefficient. Therefore, µ r,ex ‡ is given by Defining the overpotential η as 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 α c and α a respectively, as α c = α and α a = 1 − α such that α c + α a = 1. We note that our particular choice of µ r,ex ‡ in Equation 11 corresponds to choosing the "mechanical transfer coefficient" α m defined in [156] to be equal to α a , causing j 0 to not depend explicitly on 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 A, 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ĉ ± . Therefore, ∆φ eq and j 0 simplify to To compare the reaction and diffusion rates, we define the Damköhler number Da as the ratio of the Faradaic current density scale e p k 0 and limiting current density J lim given by When Da is large, i.e., Da 1, the system is diffusion-limited but when Da is small, i.e., Da 1, the system is reaction-limited.

C. Boundary conditions, constraints and initial conditions
We use "a" and "c" superscripts to denote the anode and cathode respectively, r = z + e −n · u(r = r a,c m ). For galvanostatic conditions in which we apply a current I a on the system, we requirén · J(r = r c m ) dS c =´−n · J(r = r a m ) dS a = I a to satisfy charge conservation whereas for potentiostatic conditions in which we apply an electric potential V on the cathode, we set where c 0 is the initial neutral salt bulk concentration [180], and x a m (t = 0) = 0 and x c m (t = 0) = L x , i.e., the anode and cathode are initially planar.

A. 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(y, z, t) where h is the electrode surface height. Given h, we can derive explicit expressions for surface variables such as the curvature H and normal interfacial velocity v In in terms of h and its spatial and temporal derivatives [184,185], which are provided in Section I of Supplementary Material.
For brevity, we let k = [k y , k z ] T and ξ = [y, z] T where k is the wavevector, and k y and k z are the wavenumbers in the y and z directions respectively. Therefore, k · ξ = k y y + k z z, where · 2 is the L 2 -norm and k 2 is the overall wavenumber, and the wavelength λ is given by λ = 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 where 1 is a dimensionless small parameter, the (0) and (1) superscripts denote the base and perturbed states respectively, (·) gives the real part of a complex number, h (1) is the complex-valued perturbation amplitude of the electrode surface height, and ω 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 φ are similarly given by where c − and φ (1) are the complex-valued perturbation amplitudes of anion concentration and electrolyte electric potential respectively.
To evaluate c − and φ and their gradients ∇c − and ∇φ at the interface at x = h, we require their Taylor series expansions around the base state interface at x = h (0) . Lettinĝ = exp(ik · ξ + ωt) and θ ∈ {c − , φ}, these expansions are given by [149,150,163] After substituting these perturbation expressions into the full model in Section II, we obtain the base and perturbed states by matching the O(1) and O( ) terms respectively.
The dispersion relation ω(k) is subsequently computed by solving these O(1) and O( ) equations. The growth rate ω is generally complex-valued and for a particular k value, there is an infinite discrete spectrum of ω values. However, for linear stability analysis, we are primarily interested in the maximum of the real parts of the ω values, which is denoted as max{ (ω)}, that corresponds to the most unstable eigenmode. If max{ (ω)} < 0, the perturbation decays exponentially in time and the base state is linearly stable. On the other hand, if max{ (ω)} > 0, the perturbation grows exponentially in time and the base state is linearly unstable. Lastly, if max{ (ω)} = 0, the perturbation does not decay nor grow and the base state is marginally stable.

B. Nondimensionalization
To make the equations more compact and derive key dimensionless parameters, in Table I, we define the scales that are used for nondimensionalizing the full model in Section II and the perturbation expressions in Section III A.L y andL z are the aspect ratios in the y and z directions respectively. For convenience, we also define parameters emerge from this nondimensionalization process, namely the Damköhler number Da =k 0 that is described earlier in Equation 14 and the capillary number Ca that is given by which is the ratio of the capillary constant Ωγ k B T [22][23][24] to the inter-electrode distance L x , andγ 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 (0) and (1) superscripts, which denote the base and perturbed states respectively, as 0 and 1 subscripts respectively. Similarly, we drop the 0 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, 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 C and III D below.
The equations for the base state are obtained by substituting the perturbation expressions in Section III A into the full model in Section II and matching terms at O(1). Equivalently, the base state is simply the full model specialized to 1D in the x direction with the curvaturerelated terms dropped, which only appear at O( ). Therefore, at O(1), the governing PDEs (partial differential equations) are given by 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 a 0 are given by where Since the unit normal at the cathode points in the opposite direction from that at the anode, the signs of the expressions involvingn at the cathode are opposite to that at the anode. Therefore, the boundary conditions at the cathode at x = h c 0 are given by We pick u x (x = h a 0 ) and u x (x = h c 0 ) such that the positions of the anode and cathode in the base state remain stationary, i.e.,ḣ a 0 =ḣ c 0 = 0. Therefore, u x = β vn · J 0 (x = h a 0 ) = −β vn · J 0 (x = h c 0 ) where the second equality automatically holds true because of charge conservation in the 1D O(1) 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 Sinceḣ a 0 =ḣ c 0 = 0, h a 0 (t) = 0 and h c 0 (t) = 1 at all t. For galvanostatic conditions in which we apply a current density J a on the system, we impose impose φ c e = 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 ρ s [180]. Specifically, c 0 , φ 0,x and their spatial derivatives can be analytically expressed in terms of the Lambert W function [186]. On the other hand, φ 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 a and ρ s and is found by numerically solving the algebraic Butler-Volmer equations given by Equations 27 and 30 with MATLAB's fsolve or fzero function.

D. Perturbed state
To derive the equations for the perturbed state at O( ), we substitute the perturbation expressions in Section III A into the full model in Section II and match terms at O( ). 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 O( ), the governing ODEs (ordinary differential equations) are given by 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 The boundary conditions at the anode at x = h a 0 are given by where theD 1 ,D 2 andD 3 parameters arê Because the unit normal at the cathode is in the opposite direction from that at the anode, the signs of the expressions involvingn at the cathode are opposite to that at the anode.
Hence, the boundary conditions at the cathode at x = h c 0 are given by where theĜ 1 ,Ĝ 2 andĜ 3 parameters arê The capillary number Ca = γ appears in theD 1 andĜ 1 parameters in the form of γ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 φ 0,x fields sets the scale for the critical wavenumber k c , which is the wavenumber at which the perturbation growth rate ω is 0 and the electrode surface is marginally stable.

E. 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 ∆x = 1 N −1 using second-order accurate finite differences [187]. 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 where Y, Z ∈ R (2N +2)× (2N +2) , v ∈ C 2N +2 , ω ∈ C and the second subscript in c 1,i and φ 1,i for i = 1, 2, . . . , 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 , φ 1 , h a 1 and h c 1 evaluated at the grid points, and the eigenvalue is the complex-valued growth rate ω. 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 [149].
The infinite eigenvalues are physically irrelevant to the linear stability analysis [188,189], 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 [190]. 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 [188,191,192] but they are more intrusive and require more extensive matrix manipulations as compared to the mapping technique [190] 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 O(N 3 ) while that for finding k ≤ N of them, where k = 1 in our case, scales as O(kN 2 ), the computational cost is dramatically reduced by a factor of O(N ) if we use an eigenvalue solver that can find subsets of eigenvalues and eigenvectors such as MATLAB's eigs solver.

F. Numerical implementation
The equations for the time-dependent base state in Section III C 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 E are then solved for using the eigs function in MAT-LAB R2018a. When the 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 [193,194], which is more robust than the eigs function. The colormaps used for some of the plots in Section IV are obtained from BrewerMap [195], which is a MATLAB program available in the MATLAB File Exchange that implements the ColorBrewer colormaps [196].

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 ρ s , Da and J a under galvanostatic conditions. ρ s = 0 corresponds to the classical case of an uncharged nanoporous medium while ρ s = 0 allows us to depart from this classical case and study its effects on the linear stability of the electrode surface. Experimentally, ρ s can be tuned via layer-by-layer deposition of polyelectrolytes [179,197,198] or tethered immobilized anions [76]. 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 a > 1 is applied on a classical system with ρ s = 0, as discussed in Section II A, the Sand's time t 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 s .
For the results discussed in Sections IV B, IV C, IV D and IV F below, we assume the following dimensional quantities for a typical electrolyte in a typical nanoporous medium: T = 298 K, M m = 6.941 g/mol (arbitrarily pick lithium metal) [199], ρ m = 0.534 g/cm 3 (arbitrarily pick lithium metal) [199], L x = 60 µm, L y = L z = 100L x = 6 mm, c 0 = 10 mM (initial neutral salt bulk concentration), c Θ + = 1 M = 10 3 mol/m 3 (standard concentration) and γ = 1 J/m 2 (typical surface energy of metal/electrolyte interface) [149]. Corresponding to these dimensional quantities, all dimensionless parameters that are kept constant for the results in Sections IV B, IV C, IV D and IV F are given in Table II.

A. 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 ω(k), we expect to see some local maxima or possibly just a single global maximum, which we denote as {k max , ω max }, where the electrode surface is maximally unstable. We also expect to see a critical wavenumber k c corresponding to ω = 0 where the electrode surface is marginally stable. When k is larger than k c , ω is always negative because the surface energy stabilizing effect always dominates when the wavenumber is sufficiently large. We note that k c is always greater than k max . Corresponding to k max and k c are the maximum wavelength λ max = 2π kmax and critical wavelength λ c = 2π kc respectively. In a porous medium, the characteristic pore size h c = 2d p , where d p is the pore diameter, sets a threshold or cutoff for overall electrode surface stabilization: we should observe stabilization if h c is smaller than λ c [158]. If h c is larger than λ c , then the most unstable eigenmode dominates the electrode surface growth with a growth rate of ω max and the characteristic length scale of this instability is λ max . Therefore, {k max , ω max } and k c are the most physically informative points of the dispersion relation. We now derive an approximation for the dispersion relation ω(k) that is valid at high values of k and will be useful for computing {k max , ω max } and k c quickly and accurately because k max and k 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 B.
When k is sufficiently large, at the cathode at x = h c 0 = 1, we expect k 2 c 1 to balance c 1 , and k 2 φ 1 to balance φ 1 in Equations 35 and 36 respectively. Therefore, k −2 is a small parameter multiplying the highest order spatial derivative terms c 1 and φ 1 , and the spatial profiles for c 1 and φ 1 form a boundary layer with characteristic thickness k −1 . Hence, as an ansatz for the boundary layer analysis, we assume where A and B are arbitrary constants that are determined from the boundary conditions at x = h c 0 = 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. Imposing the boundary conditions at x = h c 0 = 1, we obtain where we define α 1 = D − D + , α 2 = z + D + − zD and α 5 = α 2 c 0 − z + D + ρ s for brevity.
Approximate values of {k max , ω max } can be obtained by solving ω (k) = 0 and requiring ω (k) < 0 where the primes indicate total derivatives with respect to k. In addition, by solving ω(k) = 0, we can obtain approximate values of k 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 {k max , ω max } and k c , we first find a simpler and yet still accurate analytical approximation for k c , which can then used as an initial guess for numerically solving for {k max , ω max } using Equation 49 with MATLAB's fminbnd optimizer. Such an approximation can be obtained by assuming k c is large enough We observe that k c scales as Ca − 1 2 = γ − 1 2 , which is expected because the surface energy stabilizing effect appears in the form of γk 2 inĜ 1 in Equation 42, and this scaling agrees with that obtained in previous work done on linear stability analysis of electrodeposition [149,150,158,163,200].  Because we are mostly interested in the k max , ω max and k c points on the ω(k) curve, we plot them against N in Figure 3. We observe that the numerically computed k max , ω max and k 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 max and k 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.

C. Parameter sweeps
The base state anion concentration field c 0 , electrolyte electric potential field φ 0 and electric field E 0 = −φ 0,x possess salient features that are useful for understanding the linear qualitatively, the "total amount of instability" increases with t. For ρ s = −0.05 < 0, when compared to ρ s = 0 and ρ s = 0.05, the ω curve is the smallest at a given t because of a smaller base state electric field E 0 . The ω curve also remains bounded at all t and eventually reaches a steady state that is almost attained near t = 2t s because E 0 at the cathode behaves in the same fashion. In sharp contrast, for the classical case of ρ s = 0 near t s , the ω curve grows dramatically because of the rapidly increasing E 0 at the cathode, which eventually diverges at t s and in turn causes the ω curve to diverge at t s too. Compared to this classical case, for ρ s = 0.05 > 0, because E 0 at the cathode is larger at a given t and diverges earlier than t s , the ω curve accordingly grows even more rapidly at earlier times and diverges earlier than t 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 [150]. 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 ω but this effect is very insignificant because the application of an overlimiting current implies that the system is always diffusion- In the interest of space, plots of numerically computed (ω) against k for J a = 1 (limiting current) and J 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 a = 1, the trends observed for J a = 1 are qualitatively similar to our previous discussion for J a = 1.5, except that the ω values are smaller because a smaller applied current density results in a smaller electric field at the cathode. For J 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 ω curve remains bounded at all t and reaches a steady state eventually. Moreover, ω generally increases with Da, and this increase is especially pronounced when Da increases from 1 to 10. This is because as discussed in Section II B, the system becomes diffusion-limited when Da 1, thus increasing the ion concentration gradients and electric field at the cathode and resulting in a larger growth rate.
As discussed in Section IV A, at each t point, each ω curve exhibits a global maximum {k max , ω max } and a critical wavenumber k c , which is where the curve crosses the horizontal axis ω = 0. The {k max , ω max } and k c points provide a succinct way to summarize the most physically significant features of the ω(k) curve for all the parameter ranges we have explored thus far. Therefore, for ρ s ∈ {−0.05, 0, 0.05}, Da ∈ {0.1, 1, 10} and J a ∈ {0.5, 1, 1.5}, we plot numerically computed k max and ω max against t ts in Figure 6 and numerically computed k c against t ts in Figure 7. For J a ≥ 1, we observe that the k max and ω max curves diverge near t s for ρ s ≥ 0 but level off to constant values past t s for ρ s < 0, therefore these curves appear as if they are "fanning out". In contrast, for J a < 1, the k max and ω max curves level off past t s for all values of ρ s as the system eventually reaches a steady state when an underlimiting current is applied. The k c curves have the same qualitative shape as the k max curves except that they are larger, as expected. The effects of Da and J a on the k max , ω max and k 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 ρ s < 0 to determine how much additional stabilization a negative ρ s confers to the surface as it gets increasingly more negative. Subsequently, we plot numerically computed k max , ω max and k c against t ts for ρ s ∈ {−1, −0.75, −0.5, −0.25, −0.05}, Da = 1 and J a = 1.5 in Figure 8. While a more negative ρ s generally decreases k max , ω max and k 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 ω max remains positive even for ρ s = −1, albeit at a small value. In practice, it is probable that a sufficiently small and positive ω max value can be deemed to be small enough for considering an electrode surface "practically stable".

D. 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 max , ω max and k c against t ts for ρ s ∈ {−0.05, 0, 0.05}, Da = 1 and J a = 1.5 in Figure 9. In the interest of space, these plots for other values of Da and J a are provided in Figures 3 to 8 of Section VI of Supplementary Material. For all parameter ranges considered, the agreement between numerical and approximate values of k max , ω max and k c is excellent, giving us confidence that the approximations are useful for rapidly and accurately computing k max , ω max and k c . This confirms that k max and k 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 E and IV F that follow.   Bottom row: Plots of k c . In the legends, "num." refers to numerical solutions while "approx." refers to approximate solutions.

E. 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 [179] 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 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 For multistep electron transfer reactions, it is more convenient to work with ∆φ instead of η. Therefore, we rewrite J F in terms of ∆φ as where k c and k a are the cathodic and anodic rate constants respectively.
The reaction mechanism for copper electrodeposition and electrodissolution is given by [103,[201][202][203] Cu 2+ (aq) + e − Cu + (ads), 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 F in terms of η, we obtain where α 1 is the charge transfer coefficient of the first step.
Previously in Section II B, for a 1-step n-electron transfer metal electrodeposition reaction, the dimensionless forms of Equations 12 and 13 are given by By comparing Equations 57 and 58 with Equations 59 and 60, we set n = 2 and α = α 1 2 and replace γ with 2γ 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 currentvoltage relations, we have previously performed parameter estimation [180] for copper electrodeposition in a copper(II) sulfate (CuSO 4 ) electrolyte in cellulose nitrate (CN) membranes [179], which are a random nanoporous medium with well connected pores. The parameters that are estimated are ρ s , τ , Da, α 1 and p and their fitted values are provided in Table III in [180]. Other parameters specific to the copper electrodeposition reaction, CuSO 4 electrolyte and CN membranes used are also provided in Tables I and II in [180].
For the surface energy of the copper/electrolyte interface, we use dimensional γ = 1.85 J/m 2 given in Table I in [163].
For our analysis here, the specific experimental datasets that we focus on are labeled CN 2 (−) and CN 2 (+) in [180], which correspond to negatively and positively charged CN  To summarize the model predictions, we plot approximate dimensional values of λ c and λ max against the dimensional applied current I a in Figure 10. In the λ c plot in Figure 10(a), we also indicate the characteristic pore size h c of 0.5±0.1 µm, which is given by twice the pore diameter d p of 250±50 nm [179], in order to determine if the model predicts overall electrode surface stabilization. As discussed in Section IV A, we expect overall electrode surface stabilization if h c < λ c , which corresponds to the blue shaded region in the λ c plot. On the contrary, we expect overall electrode surface destabilization if h c > λ c , which corresponds to the red shaded region in the λ c plot, and the characteristic instability wavelength is 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.

F. 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 [204,205] or pulse charging of lithium metal batteries (LMBs) and lithium-ion batteries (LIBs) [206][207][208][209][210][211][212][213][214] 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 [204,205]. Using an applied direct current to perform silver electrodeposition under galvanostatic conditions, Aogaki experimentally observed that the crystal grain size decreases with time [137,138], which agrees well with theoretical predictions from linear stability analysis previously done by Aogaki and Makino [136]. 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 E, we generally expect the characteristic pore size h c to be larger than the critical wavelength λ c at high applied currents, therefore the electrode surface is unstable with a characteristic instability wavelength λ max . Because a pulse current is applied, λ max varies in time and hence, it would be useful to define an average λ max that averages out the effect of time. In this spirit, we define the average maximum wavenumber k max and the corresponding average maximum wavelengthλ max as where t f is the final time of the pulse and each maximum wavenumber k max is weighted by its corresponding maximum growth rate ω max . We expectλ 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 a with an "on" (charging) time of ∆t on , a "off" (relaxation) time of ∆t off , and a period T given by T = ∆t on + ∆t off . The duty cycle γ dc is given by γ dc = ∆ton T and the average applied current densityJ a over one period is given byJ a = J a,p γ dc where J a,p is the peak applied current density. Hence, for a particularJ a , a smaller γ dc implies a larger J a,p .
For the classical case of ρ s = 0, we pickJ a = 1 and ∆t on = 0.0125t s and vary γ dc from 0.2 to 1 (direct current) where the Sand's time t s is calculated based onJ a .J a , ∆t on and γ dc should be carefully chosen such that J 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 γ dc < 0.2 for our choice ofJ a = 1 and ∆t on = 0.0125t s cannot be numerically simulated. For ρ s = −0.05, we pickJ a = 1.5 and ∆t on = t s and vary γ dc from 0.1 to 1 (direct current) to drive the system at an overlimiting average applied current density.
We also fix Da = 1 for both cases and use Equations 49 and 50 to compute approximate values of k max and ω max . For these choices of parameters, as an illustrative example, we plot J a , approximate k max and approximate ω max against t for γ dc = 0.5 in Figure 11. We note that the large overshoot in k max at the beginning of each "on" cycle for ρ s = 0 is caused by the sharp rate of increase of the concentration gradients and electric field as J a rapidly increases from 0 in the "off" cycle to J a,p in the "on" cycle. Corresponding to these pulse waves, we plotλ max against γ dc in Figure 12. For both ρ s = 0 and ρ s = −0.05,λ max increases with γ dc , which agrees with the empirical observation that the crystal grain size generally decreases with applied current [204,205]. 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 [149,150,158]. 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 sur- face 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 [179]. 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 [134] 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 [215] and shock electrodialysis for water treatment [111][112][113].
We have made many assumptions and simplifications in the model presented, and relax- ing some of them offers opportunities for extending it in useful ways. First, we have ignored surface adsorption, surface diffusion of adsorbed species [139][140][141][142] and additional mechanical effects such as pressure, viscous stress and deformational stress [156,157,[159][160][161][162], 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 [216] that are used in solid state batteries. 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 [175][176][177][178], 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 [93][94][95], 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 [114,115] and current-induced membrane discharge [132] 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 [217,218] or asymmetric Marcus-Hush kinetics [219] 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.