{'fraction_non_alphanumeric': 0.04068551775350283, 'fraction_numerical': 0.01834952025482188, 'mean_word_length': 5.016223841993886, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dendrite formation in rechargeable lithium-metal batteries: Phase-field modeling using open-source finite element library

We describe a phase-field model for the electrodeposition process that forms dendrites within metal-anode batteries. We derive the free energy functional model, arriving at a system of partial differential equations that describe the evolution of a phase field, the lithium-ion concentration, and an electric potential. We formulate, discretize, and solve the set of partial differential equations describing the coupled electrochemical interactions during a battery charge cycle using an open-source finite element library. The open-source library allows us to use parallel solvers and time-marching adaptivity. We describe two- and three-dimensional simulations; these simulations agree with experimentally-observed dendrite growth rates and morphologies reported in the literature.


Introduction
Conventional lithium-ion batteries (LIB) that use intercalated graphite electrodes have dominated the rechargeable battery market over the past three decades; nevertheless, this battery technology is reaching its theoretical limit of 250 Wh/kg [1,2]. Therefore, our society has a strong need for new chemistries and designs to achieve commercial implementation of ultra-high energy density batteries of around 500 Wh/kg [3].
Currently, metallic lithium (Li) is the most prominent negative electrode material due to its combined high theoretical energy density (3860 mAh/g) and low reduction potential (-3.04 V vs. standard hydrogen electrode). However, despite recent successful operation reports for lithium metal batteries (LMBs), all reported short life cycles due to uneven and unstable lithium deposition leading to dendrite formation [4,5]. Therefore, developing a stable rechargeable lithium metal anode is crucial to attaining higher energy density rechargeable technologies, such as Li-air, Li-S, and Li-flow batteries [6,7].
Lithium dendrites form during the battery charge cycle and can disconnect from the anode, generating "dead lithium" compounds that stop participating in the electrochemical reaction, resulting in loss of Coulombic efficiency [8]. Ultimately, the uncontrolled growth of these lithium dendrites can cause the battery to short-circuit when they become large enough to contact the cathode electrode after penetrating through the separator [9,10,11].
The formation of dendrites in Li anodes has been noted to vary considerably from one system to the other. Bai et al. [38] noted the formation of mossy/whiskers-like morphologies at current densities lower than the intrinsic diffusion-limited current density (reaction limited regime, where lithium grows from the base). On the other hand, at higher current densities (greater than the diffusion-limited current density, i.e. transport or diffusion limited regime, where lithium grows from the tip), fractal dendritic morphology or finger-like structures are obtained. Transition from base-controlled mossy structure to tip-controlled dendrite growth was also observed experimentally when under limiting current condition. Further, various 2D numerical simulations have found the formation of spike-like or tree-like structures where dendrite growth was found to occur on multiple branches and in multiple directions. These types of formation, however, have not been observed widely experimentally due to instabilities at the branch structures and mass transfer rate limitations at these regions [39].
The experimental investigation of the lithium dendrite formation in rechargeable metal batteries is challenging [40]. Thus, the combined insights from experiment and simulation enhance our understanding of the mechanisms of dendrite formation and growth in lithium anodes [39,41,42]. Within these computational models, the thermodynamic ones include several surface nucleation models based on density functional theory (DFT) that explore the mechanisms of lithium electrodeposition and dendrite suppression. Considering lithium dendrite growth as an inherent property of the metal, these models invoke surface energies and diffusion barriers as keys parameters in determining the specific metal structures growing on electrodes [43,44,45]. Dynamic models describe the electrochemical propagation of dendrites; for example, the space-charge (migration-limited) models violate the electroneutrality at the dendrite tips leading to high overpotentials and the subsequent formation of ramified electrodeposits [46]. Furthermore, stress and inelastic deformation models simulate the dendrite's surface-tension mitigated growth [47,48]. Additionally, film growth models apply classical film-growth theory within the lithium anode battery [49], and diffusion-limited aggregation models (DLA) develop fractal shapes using random-walk statistics [17]. The latter successfully modeled extremely low currents without considering surface energy.
Moreover, several phase-field (diffuse-interface) models simulate the morphologic evolution of lithium electrodeposits as a result of reaction-driven phase transformation within the metal anode batteries and describe the dendrite morphologies observed experimentally for low and high applied current densities [50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67]. The phase-field method describes and predicts the mesoscale microstructural evolution by implicitly tracking the boundaries and interfaces using an auxiliary function (i.e., the order parameter). These models rigorously incorporate interfacial energy, interface kinetics, and curvature-driven phase boundary movement. The evolution of the phase-field variables satisfies local equilibrium [68] and free energy minimization [69], leading to nonlinear partial differential equations (PDE's). Different phase-field models of the electrochemical dendrite deposition describe the phase-field evolution using the Cahn-Hilliard equation [70,71], the Allen-Cahn equation [72], or a modified non-linear Allen-Cahn reaction equation [55]. In [73], the authors present a comprehensive review of the existing phase-field simulations in rechargeable batteries.
Despite continuous efforts, there are several open issues related to the evolution of dendritic patterns in lithium metal electrodes. Most published models solve two-dimensional problems to rationalize the dendritic patterns and suppress their growth [60,61,62,63,65,68]. Nevertheless, the lithium dendrite morphologies are inherently three dimensional [64,74,75]. Thus, simulation tools to study and understand the three-dimensional effects triggering pattern formation are essential. Its application includes enhancing our understanding of the differences between 2D and 3D simulations on the free space effect on the Li-ion diffusion to the anode surface, the surface anisotropy, the electric field distribution, the nuclei distribution, the interaction between neighboring dendrites; and assessing control strategies to prevent the dendrite formation (e.g., 3D porous current collectors hosts with lithiophilic sites) [76,77].
Nevertheless, few attempts to simulate three-dimensional lithium dendrite growth exist in the literature. For instance, [78] develops a 3D growth model for the electrode-electrolyte interface in lithium batteries in the presence of an elastic prestress; the model assesses the linearized stability of planar interface growth. Later, Jang et al. performed stochastic modeling of the three-dimensional effect of applied voltage and diffusion constant on the growth of lithium dendrites [79]. Recently, Yangyang et al. [65] presented a phase-field model to study the three-dimensional effect of the exchange current density on the electrodeposition behavior of lithium metal; however, without focusing on dendritic morphologies. The scarcity of three-dimensional phase-field results is evidence of the complexity and large computational cost involved in solving the highly non-linear set of equations describing coupled electrochemical interactions during a battery charge cycle.
In this paper, we present a 3D phase-field framework to describe the dendritic electrodeposition of lithium. This work constitutes a step toward the physics-based, quantitative models to rationalize hazardous threedimensional dendritic patterns needed to achieve the commercial realization of Li-metal batteries. Herein, we use several computational efficiency improvements to deliver the 3D simulations at a reasonable cost, such as time step adaptivity strategy, mesh rationalization, parallel computation, and balanced phase-field interface thickness to mesh resolution ratio.
We organize the paper as follows: Section 2 presents the basic equations describing the lithium-battery dendrite growth process. Then, we derive the weak variational statement and its numerical implementation in Section 3. We describe the time-adaptive strategy and the parameters that deliver convergence in Section 4 and detail its implementation in Section 5. We discuss numerical simulations of lithium-battery dendrites growth in Section 6, where we validate a two-dimensional model in terms of dendrite propagation rates and spatial distribution analysis of the system's variables in comparison with phase-field simulation results reported in the literature. Furthermore, we simulate three-dimensional spike-like lithium structures that grow under high current density [64,74,75] (fast battery charge); these structure's growth is dangerous for battery operation. We perform single and multiple nuclei numerical experiments to study the 3D distribution of the electric field and the lithium-ion concentration to understand the mechanism behind tip-growing lithium morphologies better. Finally, we draw conclusions in Section 7.

Governing partial differential equations
We model a battery cell composed of a solid pure Lithium metal anode and a binary liquid electrolyte. The variables of interest are the lithium concentration, ζ, and the electric potential, φ. Under an electro potential difference ∆φ, the battery charges and the electrolyte dissociates in Li + cation and A − anion species; Faradic reactions occur, the current passes through the electrode-electrolyte interface, the Li + cation gains an electron and deposits on the anode surface (Li + +e − → Li). Thus, the physical processes involved in the electrochemical deposition of Lithium are charge and mass transport. The Butler-Volmer kinetics (1) is a standard phenomenological model that describes the charge and mass transport at the electrode-electrolyte interface [80,81], known as the current-overpotential equation: where i is the current density and i 0 is the exchange current density, identified as an intrinsic kinetic parameter and assumed constant in this case [65]. The applied overpotential η (ζ, φ) is the total free energy change per charge transferred [55]. The first and second term in brackets represent the oxidation and reduction reactions, respectively, whereas α is the charge transfer coefficient that characterizes the forward and reverse reactions [82]. R, T, n, and F represent the gas constant, temperature, valence, and Faraday's constant, respectively. The lithium electrodeposition rate depends on the applied current density via a Faradic reaction [83], where λ represents the electrodeposited film thickness over a time t, and C s m is the site density of Lithium metal (inverse of molar volume) (e.g., see [55,84] for further details).
where we adjust the interface mobility L σ value (constant for each simulation) for the selected voltage charging the battery. These mobility values balance the phase-field interface energy and the electrochemical reaction contribution [84]. The electrochemical reaction kinetic coefficient is [ where the quantities of Li surface energy γ and the phase-field interfacial thickness δ P F relate to the model parameters according to δ P F = 2κ0 3γ [91]. The interpolation function h (ξ) smooths the diffuse interface in the current implementation; we use a sigmoid interpolation function [92,84], which has better computational efficiency than the popular choice of polynomial interpolation function [93]: where ς is a parameter that determines the interface thickness of the interpolation function; we use ς = 20, to adjust for interpolation between ξ = 0 and ξ = 1.

2.2.
Diffusion-migration process. A diffusion-migration equation describes the motion of charged chemical species (lithium-ion) in the fluid electrolyte. The temporal evolution of ζ + , satisfies the modified Nernst-Planck diffusion equation that describes the flux of Li-ions under the influence of a concentration gradient ∇ζ + and an electric field ∇φ [46]. We add the term C s m C0 ∂ξ ∂t to the Nernst-Planck diffusion equation to account for the eliminated lithium ion from the electrolyte solution, due to electrodeposition on the solid phase (metal electrode). Thus, the diffusion equation becomes: where we interpolate the effective diffusivity by , with D s and D l as the electrode and electrolyte diffusivities, respectively.
2.3. Electrostatic Potential. We account for the electrostatic potential distribution φ using the charge continuity equation [51] (10) where i is the current density and ρ = Σ i n i F ζ i is the charge density. Experimental observations support the assumption that space-charge effects do not affect the stability of electrodeposits [94]. Therefore, we ignore the effects of the double-layer structure and assume electroneutrality, that is, where we use Ohm's law in the continuity equation, i = σ E, and σ is the conductivity and E = −∇φ is the electric field. Thus, we arrive at a Poisson equation, including a source term to represent the charges that enter or leave the system due to the electrochemical reaction [59] (11) ∇ · σ eff (ξ) ∇φ = nF C s m ∂ξ ∂t .
We interpolate the effective conductivity using where σ s and σ l are the electrode and electrolyte phase conductivity, respectively. We collect the system of equations that models the physical process in the following box.

Box 1: Dendrite growth model
Lithium-battery dendrite growth process based on phase-field theory. Find Ξ = ξ, ζ + , φ fulfilling where V is the problem domain with boundary ∂V = ∂V N ∪ ∂V D , the subscript N and D related to the Neumann and Dirichlet parts, with outward unit normal n, and I is the time interval.

Weak variational formulation and space-time discretization
3.1. Weak formulation. We now state the weak variational formulation [95,96] of the system of equations where• = ∂• ∂t denotes the time derivative of the variable • = ξ, ζ + , u, w V = V wu dV expresses the inner product, and the functions on the right-hand side are We use standard finite element spaces where each function and its gradient are square integrable.

3.2.
Time semi-discretization. We use a second-order backward-difference (BDF2) time marching scheme with an adaptive time step size. BDF2 is an implicit time marching method that requires the solution at two previous time instants; the initial step uses a first-order backward-difference method (BDF1). BDF2 has second-order accuracy and numerically damps the highest frequencies of the solution, unlike the conservative Crank-Nicolson method [95]. Liao et al. [97] showed that BDF2 is an effective time integrator for the phase-field crystal model, especially when coupled with an adaptive time-step strategy.
We discretize the time interval into t 0 < t 1 < ... < t n < ... < t f and define the time-step size as , respectively, as ξ n ,ξ n , ζ +n and˙ ζ +n , and express the time increments as ∆• = • n+1 − • n for • = ξ,ξ, ζ + ,˙ ζ + . We use a second-order approximation of the time derivative at t n+1 as follows [98], Then using this definition and letting V v = (v, 0, 0), we define the ξ-residual as . We now approximate a (v, Ξ n+1 ) as a Taylor series expansion from Ξ n to obtain where O ∆ 2 represents neglected higher-order terms in the expansion and represents the directional Gâteaux derivative of the functional f in the direction •. Combining (15) and (16), we obtain a linear equation system to solve, that is Similarly, we define the weighting functions V w = (0, w, 0) and V p = (0, 0, p); the residuals R • with • = ζ + , φ and linearize the resulting residuals to obtain the linearized system of equations to update the Newton iteration. We collect the linearized system of equations that models the physical process in the following box.

Box 2: Linearized equation system
Discrete linearized equations for lithium-battery dendrite growth process based on phase-field theory.
For completeness, we summarize the Gâteaux derivatives in (19) in the following box.
3.3. Space discretization. We discretize the domain V ⊂ R 3 with Dirichlet boundary condition ∂V D and Neumann boundary condition ∂V N using finite elements. We express the domain as the union of nonoverlapping elements, We define continous piecewise polynomial functions over the discrete domain. In particular, we use linear functions over simplexes (i.e., triangles in 2D, tetrahedra in 3D).
Since the three variables share the spatial distribution of the boundary conditions, we use the shape functions N A (x) for each degree of freedom • n A at t n satisfying the Dirichlet boundary conditions to discretize each degree of freedom • = X, Z, Y corresponding to the variables ξ, ζ + , φ, respectively; thus, we have where S corresponds to the total number of unknowns in each solution variable. We define the weighting spaces using test spaces using the same functions but restricting those with support on ∂V D to be zero, thus where W corresponds to the total number of weighting functions for each variable. Using these spatial discretizations, we obtain the fully discrete algebraic problem:

Box 4: Matrix blocks
The mass matrix blocks are

Time-adaptive strategy
Time step adaptivity is highly useful in this problem, where the time step requirements vary significantly at different simulation stages. For example, initially, the simulation requires small time steps to achieve convergence during the development of the phase-field interface and thereon, their size grows a few orders of magnitude (i.e., from ∆t = 10 −6 [s] up to ∆t = 10 −1 [s]) depending on the parameters (e.g., mesh size, electrodeposition rate).
Following [99,100], we express the local truncation error for BDF2 using Taylor expansions as follows where u = ξ, ζ + , since φ's time derivative does not appear explicitly in the formulation. We use the solutions u n+1 , u n , u n−1 and u n−2 from the BDF2 scheme to estimate the truncation error of (26) using the third-order backward difference formula (BDF3) (27) ...
Thus, we estimate BDF2's local truncation error by substituting (27) into (26) to obtain Finally, we compute the weighted local truncation error as an error indicator [101] (29) where ρ abs and ρ rel are user-defined parameters that define the absolute and relative tolerances, respectively.
In our examples we set these parameters to ρ abs = ρ rel = 10 2 . The time error is computed using the maximum time error of time-dependent variables and the time-step adaptivity simply follows [102,103] where k is the time-step refinement level and ρ tol is a safety factor parameter set to 0.9 in our simulations. We summarize the time-adaptive scheme in Algorithm 1, where we define two tolerances tol max and tol min that limit the range of reduction or increments of the time-step size.

Implementation details
This section briefly discusses the implementation aspects and details the step-by-step calculations. We perform all the numerical experiments using the open-source FEniCS environment [104] using the FIAT package [105] to integrate numerically and assemble the matrices and vectors. We use four-node quadrilateral elements in 2D and eight-node hexahedral elements in 3D,; in all cases, we use consistent Gaussian quadrature. We use a passing interface package MPI4py [106,107,108,109] for parallelization and solve nonlinear using SNES combined with BiCGStab for each linear system, including a Nonlinear Additive Schwarz methods (NASM) for parallel solution [110]. Table 1 summarizes the parameters we use. We compute the time increment by gathering the unknown vectors at the master core, estimating the error for the current time step, adapting the time step size, and broadcasting its value to the other cores. This straightforward implementation is practical given the small number of processors used. We initialize the simulations using a first-order backward difference formula (BDF1) until • n−2 is different from null. We perform the simulations using a laptop with a 2.4 GHz processor with 8-core Intel Core i9 and 16 GB 2667 MHz DDR4 RAM, obtaining satisfactory results using a regular computer. Algorithm 1 sketches the time marching scheme, where rank is the core number, and master is the master core used as the communicator.
6.1. System layout and properties. Generally, the computational domain for a battery simulation comprises the anode and cathode regions as well as the space between the electrodes filled with electrolyte [111]. However, in most phase-field simulations of metal electrodeposition, including those we perform, the cathode region reduces to a current collector boundary condition on the electrolyte side of the domain (see Figure 2). We model a battery cell with a traditional sandwich architecture, undergoing a recharging process under fixed applied electro potential. We represent this cell as an l x × l y rectangular domain in 2D and as an l x × l y × l z hexagonal domain in 3D. The initial perturbation (dendrite nuclei) significantly impacts the simulation outcome. These nuclei are usually part of the problem's initial conditions. Some models include initial seeds (artificial deformation on the electrode surface) combined with surface anisotropy to propagate the dendrite shape [53,56,41]. Alternatively, other models combine surface anisotropy with Langevin noise at the phase-field interface to generate a random nucleation pattern [60,62,63,68]. In our experience, the shape of the initial seed and the noise levels have a major influence on the resulting dendrite morphology.
We study two different initial structures. The first structure consists of a 5µm-thick metal anode (l 0 = 5 [µm]), made up of pure Lithium, separated from the liquid electrolyte by a smooth interface, as Figure 2-a depicts. The initial condition drives from the equilibrium solution for a one-dimensional transition zone between solid (ξ = 1) and liquid (ξ = 0), where our variables ξ, ζ + , φ vary in the "x" spatial direction In the second structure we assume that artificial nucleation regions, ellipsoidal protrusions (seeds) with semi-axes r x , r y , r z , and center l 0x , l 0y , l 0z , exist at the surface of the anode undergoing electrodeposition (see Figure 2-b). This is a widely-used strategy in phase-field simulations of electrodeposition [58,59,62,63]; since it reduces the computational cost as the lithium metal is only able to electrodeposit on the same nuclei, enhancing dendrite growth and allowing for detailed study of its morphology. For the artificial nucleation case, we modify the initial condition formula, replacing "x" by the hyperbolic tangent argument, to account for a smooth transition between the solid seed (lithium metal anode) and the surrounding liquid electrolyte region. We assume the cell's electrolyte to be 1M LiPF 6 EC/DMC 1:1 volume ratio solution [61]; we compute the electrolyte's site density (C l m ) using its density (1.3 g/cm 3 ) and molar mass (90 [g/mol]) and electrode's site density C s m using the density (0.534 g/cm 3 ) and molar mass (6.941 [g/mol]) of pure Lithium [61]. We model the anode boundary as a Dirichlet boundary condition ξ = 1 at x = 0 for the phase-field parameter (solid electrode phase). In contrast, we use a no-flux Neumann boundary condition at x = l x (cathode) to allow the electrodeposition process to take place (ξ changing from 0 to 1) when the reaction front approaches the short-circuit condition [84]. For the Li-ion concentration, we apply Dirichlet boundary conditions setting ζ + = 0 and ζ + = 1 at the anode (x = 0) and cathode (x = l x ) boundaries, respectively. These boundary conditions allow the Li-ion flux into the battery (electrolyte side), ensuring that the Li deposited at the electrode-electrolyte interface equals the amount of Li + supplied to the electrolyte. Thus, these boundary conditions avoid the quick Li-ion depletion and keep the electrodeposition process running for the full simulation [84]. Finally, we apply periodic boundary conditions on the remaining faces; these are (x, 0, z), (x, l y , z), (x, y, 0), and (x, y, l z ). Figure 2 summarizes the boundary conditions we apply.
Additionally, Table 2 presents the phase-field model parameters. The normalization constants for length, time, energy and concentration scales are set as h 0 = 1 [µm], t 0 = 1 [s], E 0 = 2.5 × 10 6 J/m 3 , and C 0 = 1 × 10 3 mol/m 3 , respectively [61]. Table 2 also shows that although the electrode and electrolyte materials can exhibit Li/Li + dependent conductivities and diffusivities, their values are set to constants across each phase for simplicity [61,60]. 6.2. Lithium dendrite propagation rate: Consistency with experimental data. Validation of phasefield models with experimental results is a well-known challenge of electrodeposition simulations [60]; partly due to the computational cost of simulating the detailed lithium electrodeposition at the whole-cell scale [62]. One limiting factor is the domain size, defined by the inter-electrode separation distance in experimental cells, which ranges from 1 to 10mm [83,116,117].  We perform 2D simulations to directly compare our simulations against experimental results obtained from thin-cell geometries. We use 2D simulation results to validate the predicted lithium dendrite propagation rates, assuming that the electrodeposition rates are not significantly affected by the problem's size.
We achieve large simulation sizes through several computational efficiency improvements. We use a sigmoid interpolating function [84] instead of the popular polynomial interpolating function [58,56,59]. We use time-step size adaptivity and mesh mapping in the region of interest (close to the electrode-electrolyte interface). We also use parallel solution strategies and carefully select the phase-field interface thickness (δ P F ), given the interface thickness effect on the electrodeposition rate, which ultimately determines the evolution dynamics (motion) of the lithium electrodeposits [84].
We compare the simulated lithium dendrite propagation rate against the experimental measurements of Nishikawa et al. [83]. Nishikawa et al. studied the dendritic electrodeposition of lithium as a function of total charge (electric current passed over time) using microscope observations. They reported dendrite growth rates of 0.01 − 0.06 µms −1 [34] using a 1M LiPF 6 electrolyte, establishing upper and lower bounds for comparison with our simulation results. The wide range of values reported in many studies [29,83,118] may be caused by the uneven current distribution on the electrode surface [83].  We simulate a lithium electrodeposition process at an applied current density of i = 10 mA/cm 2 [83]. We use a 2D rectangular domain l x = 5000 [µm] long (distance betweenelectrodes), and l y = 50 [µm] wide (sufficient width to allow dendrite formation). We use a mesh 8, 000 × 140 elements, with mesh distribution that guarantees an x-spatial resolution of approximately 0. 35 [µm] (h < 4δ P F ) in the region of interest of the domain with square elements. We apply a charging electro potential of φ b = −0.45 [V ] to the cell (a common value in the literature [59,61,60]), well above the threshold overpotential for dendrite growth [61]. We perturb the initial electrode surface by (±0. 5 [µm]) in x-direction according to 0.4 sin (0.44y + π) e −0.004(y−0.5ly) 2 , consisting of a periodic function combined with a Gaussian distribution to promote dendrite formation at the center of the domain ( ly 2 = 25 [µm]). Figure 3 shows the initial conditions (t = 0 [s]) and the evolution of the system's variables ξ, ζ + , φ . We zoom in the phase-field variable ξ at each instant to emphasize the growth of the lithium dendrite and show the instantaneous dendrite length measurement. Figure 4 shows the evolution of the simulated lithium dendrite length over time (blue). While in the range of experimental data (green), the results are sensitive to differences between the transport and kinetic properties and those of the experiments [83] and our simulator. Therefore, further well-controlled experimental studies, with detailed characterization of the system parameters, combined with modeling will be necessary to improve the correlation [34].  . Comparison between 2D phase-field simulations of lithium dendrite growth at a current density of i = 10 mA/cm 2 (blue), experimental data taken from Nishikawa et al. [83] (green), and the Akolkar analytical model predictions [34] (red). 1D planar electrode phase field simulation for reference (black).
Additionally, Figure 4 includes analytical dendrite propagation results for the tip current density developed by Akolkar [34] for comparison. Akolkar's model calculates the dendrite propagation rate by analyzing various overpotentials that develop at the dendrite tip and the flat electrode surface, assuming a current density of i = 10 mA/cm 2 applied and a constant radius of r = 1 [µm] at the dendrite's tip. This simplification constitutes departs from our phase-field model, where the dendrite's tip radius is variable throughout the simulation. This difference may partly explain the higher growth rate obtained in our simulations. Furthermore, we also include 1D planar electrode phase-field results using an identical set-up as to the 2D simulation [84].
6.3. 2D Simulation of lithium dendrite growth. The following 2D numerical experiment compares our model's features to other phase-field simulations reported in the literature. This comparison verifies our work as a preliminary step to continue with a series of larger and more complex 3D simulations of lithium dendrite growth in metal anode batteries latter in this paper.
In this case, we use a 2D square domain set as 360 × 360 µm 2 , within the range of typical distances between electrodes in metal anode batteries [6]. We use a 500 × 720 structured mesh with a mapping in the x-direction to obtain square elements of size 0.5 × 0.5 µm 2 in the region of interest (see Section 6.4 for details). We apply a charging electro potential of φ b = −0.7 [V ] to the cell, with an artificial nucleation region with 3 protrusions (ellipsoidal seeds), equally spaced, growing from the left boundary, with semi-axes Figure 5 shows the evolution of our system's variables ξ, ζ + , φ at different time steps. The corresponding temporal-spatial distribution of the phase-field variable ξ depicts the dynamic morphological evolution of the simulated lithium dendrite.   Figure 5 shows the evolution of bush-like lithium dendrite formation, having morphological resemblance to the dendritic patterns reported in previous phase-field studies of dendrite growth [62,66]. The lithium dendrites grow from each nucleation site and move towards the opposite electrode, following the electric field ( E = −∇φ). As the dendrite grows, side branches appear, growing across the electrolyte, in agreement with lithium dendrite experiments [38]. The computed width range of lithium branches is between 5 and 13 [µm]. These bush-like morphologies grow fast across the electrolyte region and penetrate through porous separators, becoming potentially dangerous as they can produce battery short circuits [6]. The spatial distribution of the lithium-ion concentration ζ + , and electric potential φ, corresponding to the second and third rows of Figure 5 respectively, show that they follow the same trend reported in other phase-field models of lithium dendrite growth [62,56]. The lithium-ion concentration remains equal to the bulk concentration, ζ + = 1, through the electrolyte phase and decreases near the dendrite front (solid phase) due to the electrodeposition process [59]. The electric potential has a large gradient that spreads over the electrolyte phase and remains constant and equal to φ b at the electrode phase. As the lithium dendrite moves across the electrolyte region, the gradient of the electric potential (electric field) increases, in agreement with simulation results reported by Hong et al. [56].    6.4. 3D simulation of lithium dendrite growth: Single nucleus experiment. This section describes 3D phase-field simulations of lithium dendrite formation to study three-dimensional highly branched "spikelike" dendritic patterns, commonly observed experimentally. These patterns form under high current density loads, which correspond to fast battery charge [64,74,75]. We select a geometrical unit that characterizes a real cell structure [62,111]. We choose a computational domain of 80 × 80 × 80 µm 3 . Consequently, given the domain size, we expect growth rates up-to two orders of magnitude faster-than-normal, due to the short separation between electrodes l x = 80 [µm]. We use a 3D structured mesh with eight-node hexahedral elements. We distribute the mesh to focus the node's mapping on the area of interest (see Figure 7a). In particular, in the x-direction x r = (2/π) arcsin(x u ); where x u is the node's x coordinate normalized by l x , before mapping (uniform distribution), and x r is the node's mapped coordinate. The arcsin function transitions smoothly, inducing a relatively small variation in the mesh size within 50% of the physical domain where the lithium electrodeposit process occurs [62,63,67]. The mapping produces a finer mesh to properly capture the phase-field interface thickness (4 elements in the interface [84]) and the steepest gradients of ζ + and φ. We use a 120 3 tensor-product mesh with a mesh size of 0.4 [µm] in the region of interest (bottom half of the domain). We partition the mesh into eight processors. Figure 7a identifies each core with a different color, showing that the tensor-product mesh can efficiently allocate resources in the region of interest (more resources/colors allocated to the bottom half of the domain). As in the previous 2D simulation, we apply a charging electro-potential value of φ b = −0.7 [V ] to the cell, with an artificial nucleation region formed by a single protrusion (ellipsoidal seed), with its center at (0; l y /2; l z /2), and semi-axes 4 [µm] × 2 [µm] × 2 [µm] as Figure 7b shows. For visualization purpose, here we rotate the 3D figures to depict a vertical x-axis, showing an upright lithium dendrite growth, in agreement with the convention used for lithium dendrite experiments [83,74,75]. Figure 8 shows the morphological evolution of the simulated lithium dendrite (isosurface plot of the phasefield variable ξ). The simulation forms a spike-like, symmetric, and highly branched pattern, consisting of the main trunk and sets of four equal side branches growing in each horizontal direction. A result consistent with the body-centered cubic (bcc) crystallographic arrangement of lithium metal [62]. Our model includes the four folded surface anisotropy (6). In-situ optical microscopic investigations [74,75] report spike-like lithium dendrite formation. These dendritic patterns grow when high (over-limiting) current densities are applied to the cell (> 20 mA/cm 2 ) [6]. In such cases, the rate of lithium deposition overcomes the rate of solid-electrolyte interface formation, allowing the lithium deposit to grow almost free from the influence of the interface [6]. The enrichment of lithium-ion concentration appears in the vicinity of dendrite tips, reaching peak values of up to ζ + = 2.1, and triggering tip-growing dendritic lithium. Electrolyte regions with higher lithium-ion concentration ( ζ + > 1) are represented by orange volumes in Figure 8. It is worth mentioning that this phenomenon was previously reported by Hong et al. [56] in 2D simulations of lithium dendrite formation, proposing a compositionally graded electrolyte, which could potentially suppress the dendrite initiation.
We calculate the electric field distribution by differentiation of the resolved electric potential E = −∇φ h . The magnified view in Figure 9 shows how the electric field localizes in the vicinity of the dendrite tip, leading to an enriched concentration of the lithium-ion it induces due to the strong migration from the surrounding regions (see (9) and [56]). Figure 10 shows the performance of the time-adaptive scheme, throughout the 0.6 [s] of simulation. Starting with a small time-step of ∆t 0 = 10 −7 [s] to achieve convergence, followed by a rapid increase in size, until reaching a stationary value of about ∆t n+1 = 10 −3 [s]. The evolution of the weighted truncation error e n+1 (blue) stays close to the minimum tolerance limit (10 −7 ) during 0.3 [s] of simulation; beyond this point, the error estimate increases due to the acceleration of lithium dendrite propagation rate as the dendrites approach the positive electrode. The time-step size remains unchanged since the error estimate remains within the error bounds.
The construction of phase-field models satisfies an a priori nonlinear stability relationship, expressed as a time-decreasing free-energy functional; nevertheless, standard discrete approximations do not inherit this stability property. Thus, Figure 11 shows the evolution of the Gibbs free energy of the system Ψ, see (3), using our adaptive time integration scheme. We plot the total energy curve (black), as well as three additional energy curves that correspond to each one of its terms, namely, the Helmholtz (chemical) free energy V f ch dV (blue), surface energy V f grad dV (green), and electrostatic energy V f elec dV (red), as the figure indicates.   Figure 11. Energy time series for 3D spike-like dendrite growth simulation. The inset plots the increasing surface energy in smaller scale for better appreciation. Figure 11 shows that the total systems' discrete free energy does not increase with time. Therefore, we obtain discrete energy stable results using our second-order backward-difference (BDF2) time-adaptive marching scheme. Alternative, provably unconditionally stable second-order time accurate methods may deliver larger time-step sizes for phase-field models [119,120,121,122,99], but these are beyond the scope of this work. Additionally, the system's chemical and surface energies increase as the lithium surface area grows as time progresses. In parallel, the electrostatic energy decreases in time. This interaction is consistent with the electrodeposition process, where the system stores the applied electrostatic energy as electrochemical energy as the battery charges.
6.5. 3D Simulation of lithium dendrite growth: Multi-nuclei experiment. The lithium dendrite nucleation depends on local inhomogeneities that may arise from different causes, such as defects and impurities in the metal anode, imperfect contact between the electrode and electrolyte caused by the development of a solid-electrolyte interface, and variations in the local concentration or temperature in the electrolyte [61]. Given the random nature of the nucleation phenomenon, we need to deal with some degree of randomness and uncertainty when defining the artificial nuclei in the simulation. We study the simulations' sensitivity to the artificial nuclei size, shape, and proximity following 2D studies that show this  dependance [59,57,62]. Thus, the following 3D numerical experiment tests the simulation's sensitivity to the nuclei distribution and proximity by comparing four-nuclei morphological results with those previously obtained with a single artificial protrusion. The simulation setup is similar to the previous 3D experiment with a different nucleation arrangement. Therefore, we form the artificial nucleation region with four protrusions (ellipsoidal seeds), with semi-axes 4  Figure 12 shows the morphological evolution of the simulated lithium dendrite (isosurface plot of the phase-field variable ξ). In agreement with previous numerical experiments, the simulation forms a spike-like, symmetric, and highly branched pattern. The dendrite morphology consists of four main trunks growing from each nucleus, with pairs of orthogonal branches developing to the sides.
The dendrite growth does not occur perpendicular to the stack but at an angle of about 20 • between main trunks. The stacks seem to repel each other, showing morphological similarity with dendritic patterns observed in lithium experiments performed by Ding [74], see Figure 13. As a consequence of the reduced distance between electrodes (limited by the computational cost), the time scales we simulate are significantly shorter (∼ two orders of magnitude faster) than the experimentatl ones [83,74,75]. Figure 14a details the simulated dendritic morphology at time t = 0.5 [s], together with the spatial distribution of the lithium-ion concentration ζ + in the electrolyte region. No side branches form facing the center of the nucleation arrangement (between the dendrites), where the deposition process depletes the lithium-ion concentration (shown in blue). Figure 14a also shows how dendrite tips with higher lithium-ion concentration consume lithium ions from the dendrite valley [56], creating a shadow that inhibits branching growth in the spatial proximity of more developed adjacent dendrites. As a consequence, the side branches    do not grow dendrites facing the center of the nucleation arrangement (unfavorably oriented dendrites) due to the 3D interactions with other dendrites (adjacent to the main trunks) [89]. Therefore, 2D phase-field models of dendrite electrodeposition cannot capture this 3D phenomenon. Figure 14b shows the spatial distribution of the electric potential φ overlaid with the lithium dendrite at the same instant t = 0.5 [s]. In agreement with the 2D simulations of Section 6.3, the electric potential has a steep gradient over the electrolyte phase and its steepness increases as the lithium dendrite grows. Figure 15 shows the evolution of the electric field E = −∇φ and the enriched lithium-ion concentration ( ζ + > 1) in the vicinity of the dendrite tips (red volumes represent enriched concentration). As before, a larger electric field localizes in the vicinity of the dendrite tips, leading to a higher lithium-ion concentration with peak values of ζ + = 1.8. The electric migration forces overcome the concentration diffusion gradient and drive lithium cations from surrounding regions with lower concentrations (i.e., Li-ion depletion of valley regions) to accumulate around dendrite tips, triggering tip-growth with highly branched dendritic lithium.
We analyze the dendritic patterns' morphology evolution using single and multiple nuclei. Following Yufit et al. [117], we characterize the morphology by tracking the dendrites' volume-specific area ratio (µm 2 /µm 3 ) and the branch number evolution in time. Figure 16a plots the growth of the volume versus the surface   area for the 3D lithium patterns we simulate. The volume-specific area average ratios of 0.83 and 0.91 µm 2 /µm 3 of the single and multiple nuclei simulations, respectively; where a higher surface area/volume ratio is indicative of a more branched shape. Given the lack of experimental data available in the literature for quantitative characterization of the spike-like lithium morphologies, we rely on experimental results available for zinc dendrites. Yufit et al. [117] report values between 0.86 and 1.04 µm 2 /µm 3 for experimental formation of dendrites in zinc batteries. These values are slightly higher ratios than the ones we obtain for lithium. Nevertheless, these values are reasonable since zinc has a hexagonal crystal structure and tends to grow more branching (fractal) dendrites than lithium, which has a cubic crystal structure (less inherent anisotropy)and, thus, produces needle (spiky) dendrites [60]. Figure 16b compares the number of branches developed over time in each case. The simulatilons produce ratios of 48 and 72 branches per second [1/ [s]] for the single and multiple nuclei simulations, respectively. We compute the number of branches by visual inspection of the simulated morphologies, where we consider new protuberances as incipient branches under development. Direct comparison with dendrite experiments was impractical in this case, due to the faster dynamics of the simulation (small distance between electrodes). Figure 17 shows the performance of the time-adaptive scheme, throughout the 0.8 [s] of the simulation. Starting with a small time-step of ∆t 0 = 10 −7 [s] to initially achieve convergence, followed by a rapid increase in size, until reaching a stationary value of about ∆t n+1 = 10 −3 [s]. The weighted truncation error e n+1 (blue) stays close to the minimum tolerance limit (10 −7 ) during 0.4 [s] of the simulation, after which the estimated error grows due to the acceleration of lithium dendrite propagation rate as it approaches the positive electrode. The time-step size remains unchanged since the error estimate remains within the tolerances Figure 18 shows the systems' energy evolution. The total energy curve (black) does not increase with time, showing energy stable results. As before, the system's chemical (blue) and surface energies (green) increase with time as the lithium surface area increases, while the electrostatic energy (red) decreases. This energy reduction happens as the process transforms electrostatic energy and stores it as electrochemical energy (battery charge). Figure 18 shows that the surface energy for the four-nuclei simulation is almost four times larger than the single nucleus surface energy of Figure 11. The proportionately larger surface area due to the four seeds explains this scaling.

Conclusions
Using a phase-field model, we describe the dendrite formation process during metal (lithium) electrodeposition.
We present the weak variational form of the coupled equations describing electrochemical interactions during the battery charge cycle and its finite element implementation using an open-source library. We use a time-adaptive strategy and detail its implementation and parametrization. The numerical experiments demonstrate numerically that the adaptive strategy produces energy stable results when we use a minimum tolerance for the time-adaptive scheme of tol min = 10 −7 . We reduce the computational cost of simulating the detailed lithium electrodeposition at the scale of the whole cell (distance between electrodes of about l x = 5000 [µm]) by using time-step size adaptivity, mesh mapping to concentrate the mesh in the region of interest, parallel computation, and a careful selection of the phase-field interface thickness (δ P F ) considering the applied electro potential value of φ b = −0.45 [V ].
We compare our 2D simulation results (lithium dendrite propagation rate) with experimental measurements by Nishikawa et al. [83] at an applied current density of i = 10 mA/cm 2 . The growth rates we predict (∼ 0.04 µms −1 ) are within the range of experimental results (0.01 − 0.06 µms −1 [83]). Further, we compare the two-dimensional behavior of our model with previous phase-field simulations from the literature. The 2D evolution of the phase-field order parameter ξ under φ b = −0.70 [V ] applied voltage depicts the dynamic growth of "bush-like" lithium dendrites. The spatial distribution of the lithium-ion concentration ζ + , and electric potential φ, follow the same trend reported in [56,62,66].
We also perform 3D simulations of lithium dendrite formation to explain 3D highly branched "spikelike" dendritic patterns formed under high current density status (fast battery charge). We study the 3D morphological evolution of dendritic lithium under φ b = −0.70 [V ] applied voltage. We mimic the cubic crystal structure and surface anisotropy of lithium by using a 3D four-fold anisotropy model (6) of William et al. [88] to simulate crystal growth.
We analyze the sensitivity of the model to initial conditions by comparing results using different setups (single nucleus versus a four-nuclei arrangement) for accounting for the random nature of the nucleation process. We obtain spike-like, symmetric, and highly-branched patterns in both cases. We obtain volume to specific-area average ratios of 0.83 and 0.91 µm 2 /µm 3 for the single and multiple nuclei simulations, respectively; where a higher surface area/volume ratio indicates a more branched dendritic shape.
We study the 3D distribution of the electric field and lithium-ion concentration to better understand the mechanism behind spike-like dendrite growth. In both simulations, the electric field increases in the vicinity of the dendrite tips, increasing the lithium-ion concentration that peaks at ζ + = 2.1 for a single nucleus and 1.8 for multiple nuclei. Thus, electric migration forces overcome the concentration diffusion gradient, causing lithium cations to move from less concentrated surrounding regions (i.e., lithium-ion depletion of valley regions) and accumulate around dendrite tips, triggering tip-growing and highly branched dendritic lithium. In multiple-nuclei simulation we observe that the deposition process depletes the lithium-ion concentration at the center of the nucleation arrangement. Thus, 3D interactions with other dendrites creates a shadow that inhibits branching growth in the spatial proximity of more developed adjacent dendrites [89]. Therefore, 2D phase-field models of dendrite electrodeposition cannot capture this 3D phenomenon.
Finally, unlike 2D simulations, the 3D simulations performed have a smaller domain size due to the high computational cost. Consequently, the short separation between electrodes in the 3D model (l x = 80 [µm]) induce faster deposition rates. Thus, future work may explore different meshing strategies to simulate larger 3D domains to achieve realistic simulation time scales.
In conclusion, we demonstrate the proposed phase-field modeling framework's effectiveness in performing 2D and 3D simulations of dendrite formation in lithium metal batteries with a reasonable computational cost. Beyond lithium electrodeposition, the framework has the potential to model other metal deposits in metal-anode batteries, such as zinc anode batteries.

Acknowledgments
This work was supported by the Aberdeen-Curtin Alliance Scholarship. This publication was also made possible in part by the Professorial Chair in Computational Geoscience at Curtin University. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 777778 (MATHROCKS). The Curtin Corrosion Centre and the Curtin Institute for Computation kindly provide ongoing support.