{'fraction_non_alphanumeric': 0.04187192118226601, 'fraction_numerical': 0.019334975369458128, 'mean_word_length': 5.142980654922728, 'pattern_counts': {'":': 0, '<': 2, '<?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': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Three-dimensional experimental-scale phase-field modelling of dendrite formation in rechargeable lithium-metal batteries

We perform phase-field simulations of the electrodeposition process that forms dendrites within metal-anode batteries including anisotropic representation. We describe the evolution of a phase field, the lithium-ion concentration, and an electric potential, during a battery charge cycle, solving equations using time-marching algorithms with automatic time-step adjustment and implemented on an open-source finite element library. A modified lithium crystal surface anisotropy representation for phase-field electrodeposition model is proposed and evaluated through different numerical tests, exhibiting low sensitivity to the numerical parameters. Change of dendritic morphological behaviour is captured by a variation of the simulated inter-electrode distance. A set of simulations are presented to validate the proposed formulation, showing their agreement with experimentally-observed lithium dendrite growth rates, and morphologies reported in the literature.


Introduction
Global energy demand continues to rise due to industrial activity and the world's population expansion, with an average growth rate of about 1% to 2% per year since 2010 (pre-Covid19 pandemic levels) [1]. The increasing consumption of non-renewable energy reserves, such as coal, gas, and oil [2], and awareness of climate change [3,4], have triggered a steep growth in renewable energy sources (6% average annual growth worldwide over the past decade) [5], along with an urgent need for the development of improved energy storage systems [6]. Globally, around one-quarter of our electricity comes from renewables, which include hydropower, wind, solar, biomass, ocean energy, biofuel, and geothermal [1].
New chemistry and designs, such as metal anode batteries, are under active research to achieve an energy density of 500 Wh/kg and manufacturing costs lower than $100/kWh [7]. Despite enormous efforts, today's highest energy density remains below 400 Wh/kg, with an average growth rate of about 5% per year since 1970 [8]. As the energy density limitation (299 Wh/kg) of conventional lithium-ion batteries based on intercalated graphite anode cannot meet the current market demand, researchers are refocusing on lithium metal batteries (LMBs) [9]. LMBs can achieve ultra-high energy densities by avoiding the use of a graphite lattice to host Li + (intercalation process), as illustrated by the comparative schematic of Figure 1. The graphite material (host) drastically reduces the energy density of conventional Li-ion batteries by adding extra weight to the battery pack that does not participate in the electrochemical reaction [8,10].
(a) Li-ion Battery.   Although lithium metal has been an attractive anode alternative in rechargeable batteries since the early 1970s, its commercialization has been hindered due to several shortcomings. The greatest challenge to achieving the commercial realization of lithium-metal batteries is related to their stability and safety [10]. These issues are closely linked to lithium anode: dendrite formation due to the uneven deposition of lithium, dead lithium formed after dendrites breakage, formation of unstable solid electrolyte interphase (SEI), and volume expansion of the metal anode. Additionally, these mechanisms interact, causing synergistic detrimental effects [11].
Dendrite formation in LMBs is the consequence of lithium's uneven deposition, associated with thermodynamic and kinetic factors, such as the inhomogeneous distribution of Li-ion concentration and electric potential on the electrode surface. Furthermore, the morphology of the electrodeposited lithium is influenced by different factors such as the magnitude and frequency of the applied current density, electrolyte concentration, temperature, pressure, ion transport, and mechanical properties in the electrolyte [12,13]. Understanding dendrite formation in LMBs combines theory, experiment, and computation [14]. Within computational research, based on classical chemical reaction kinetics, the phase-field method is a suitable modeling technique for studying mesoscale (µm) electro-kinetic phenomena, such as dendritic electrodeposition of lithium, with a reasonable computational cost. Phase-field (diffuse-interface) models can simulate the morphology evolution of Li electrodeposit due to reaction-driven phase transformation within metal anode batteries and rationalize morphology patterns of dendrites observed experimentally, under both low and high applied current densities [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44].
While progress has occurred in phase-field modeling of lithium dendrites in recent years, there are still several issues related to the evolution of dendritic patterns in lithium metal electrodes that remain unresolved [37].
The fundamental failure mechanism of lithium anode remains unclear and controversial [11]. A significant effort has gone into using 2D models to rationalize 3D dendritic patterns observed experimentally [42]. Furthermore, various strategies exist to suppress Li dendrites' growth and weaken side reactions. Some of these strategies address the battery operating conditions, including pulse charging lithium dendrite suppression [40] and control of internal temperature [27]. Other alternatives focus on the electrode (anode), such as modeling of 3D conductive structured lithium metal anode [32,43], and low porosity and stable SEI structure [29]. Besides, other approaches center on the electrolyte management and separator design, proposing a compositionally graded electrolyte [26], dendrite suppression using flow field (forced advection) [31], the study of separator pore size inhibition effect on lithium dendrite [44].
Given the inherent 3D nature of lithium dendrite morphologies [30,45,46], it is critical to develop phase-field models to understand the impact of 3D effects on triggering the formation of these patterns. Nevertheless, few papers attempt to simulate the full 3D lithium dendrite growth process using phase-field models. For instance, Mu at al. [36] performed parallel three-dimensional phase-field simulations of dendritic lithium evolution under different electrochemical states, including charging, suspending, and discharging states. Later, Liu et al. [39] presented a phase-field model to study the three-dimensional effect of exchange current density on electrodeposition behavior of lithium metal; however, without focusing on dendritic morphologies. Recently, Arguello et al. [42] presented 3D phase-field simulations using an open-source finite element library, to describe hazardous three-dimensional dendritic patterns in LMBs. The authors used time step adaptivity, mesh rationalization, parallel computation, and balanced phase-field interface thickness to mesh resolution ratio. The high computational cost of simulating the detailed lithium electrodeposition is a well-known challenge that has limited the domain size of phase-field simulations [25,28]. Thus, higher-than-normal dendrite growth rates were reported in the literature for 3D phase-field modeling of dendrite growth, due to the short separation between electrodes used [36,42]. Recently, experimental observations by Chae et al. [47] have revealed a change in the lithium deposition behavior and morphology from "hazardous" needle-and moss-like dendritic structures to "safer" morphologies (smooth and round shaped surface) as interelectrode spacing increases. Therefore, simulating the dendrite formation at the experimental scale has significant practical relevance.
In this paper, we present a 3D phase-field model that describes the anisotropic dendritic electrodeposition of lithium, building on [41,42]. This work constitutes a step toward the physics-based, quantitative models to rationalize hazardous three-dimensional dendritic patterns needed to achieve the commercial realization of Li-metal batteries [42]. We organize the paper as follows: Section 2 presents the basic equations describing the lithium-battery dendrite growth process and details its implementation where we introduce a modified representation of the surface anisotropy of lithium metal. Section 3 describes the system layout and properties, together with the implementation of symmetric boundary conditions for a detailed study of symmetric dendritic patterns [46,30]. We discuss numerical simulations of spike-like lithium-battery dendrites growth in Section 4, where we analyze the sensitivity of the simulation results for a series of spatial resolutions and phase-field interface thickness. Section 5 evaluates the behavior of the surface anisotropy representation model for metal anode battery simulations through different numerical tests. We compare thedendritic patterns with the results obtained in preceding simulation work [42]. Towards the goal of experimental-scale simulations, we test the modified surface anisotropy model under a larger interelectrode distance in Section 6. The lithium electrodeposition behavior changes as we increase the interelectrode distance. We study and describe the lithium dendrite propagation rates and morphologies obtained under different charging voltages. Finally, we draw conclusions in Section 7.

Formulation & implementation
2.1. Surface anisotropy representation for phase-field electrodeposition models. In this section, we present a modified representation of the 3D surface anisotropy of lithium crystal. We start by considering the surface energy expression, following [41,42]: f grad = 1 2 κ (ξ) (∇ξ) 2 , where its variational derivative (surface anisotropy of lithium crystal) is: δf grad δξ = κ (ξ) ∇ 2 ξ, consistent with most recent phase-field models of dendritic electrodeposition [23,24,25,26,29,37]. However, a more accurate representation of δf grad δξ may include an additional term, as originally derived by Kobayashi [48] for 2D crystal growth. In 3D, we use the variational derivative version derived by George & Warren [49] to simulate the surface anisotropy of crystal growth (1) δf grad δξ = δ δξ where a 2 = κ (ξ) is the three-dimensional surface anisotropy or gradient coefficient. The first term after the last equality remains the same as in the previous surface anisotropy expression; however, we add a second term (derivation due to κ variable with ξ). We calculate the partial derivative ∂a ∂ ∂ξ ∂x i in (1). We express the 3D surface anisotropy coefficient (four-fold anisotropy) [49,42], as: where x 1 = x, x 2 = y, and x 3 = z; κ 0 relates to the Lithium surface tension γ; and δ aniso is the strength of anisotropy [50,51]. Thus, we apply the quotient derivative rule to (2) and arrive at the partial derivative expression we use in (1); subsequently: where n i = ∂ξ ∂x i ||∇ξ|| for i = 1, 2, 3. Although its extensive use in phase-field models of crystal growth (solidification) [49,52], only Wang et al. [53] apply these models to a 2D phase-field simulation of dendrite growth in the recharging process of zinc-air batteries. This limited use of this known model is because it induces minor morphological changes in 2D electrodeposition process; compare the similarity of the 2D dendritic morphologies reported by Wang et al. [53] including the additional surface anisotropy term, and Zhang et al. [23] not using it. However, as we show later, its effect is crucial when modeling 3D dendritic growth.

2.2.
Governing equations. Based on phase-field theory, the governing equations for dendrite growth in lithium-metal batteries were discussed in [41,42]. Thus, herein, we only summarize this problem's basic equations, modifying just the phase-field equation for inclusion of the surface anisotropy model (see [41,42] for further details and references). Additionally, we include a symbols list with descriptions and associated units for reference.
The lithium-metal batteries dendrite problem using phase field formulation is simulated by the following set of equations:Find Ξ = ξ, ζ + , φ fulfilling where ξ is the phase-field order parameter, ζ + is the lithium-ion concentration, and φ is the electric potential; V is the problem domain with boundary ∂V = ∂V N ∪ ∂V D , the subscript N and D related to the Neumann and Dirichlet parts (see 2), with outward unit normal n, and I is the time interval.
2.3. Implementation details. Following [42], we 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 (FEniCS environment) [54]. We use eight-node (tri-linear) hexahedral elements. We use a message passing interface package MPI4py [55,56,57,58] for parallelization and solve nonlinear using SNES from PETSc [59]. 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 (see [42] for further details). Figure 2. System layout & boundary conditions for artificial nucleation simulations as defined in [42]. Reproduced with Journal's permission

System layout & properties
Generally, the computational domain for a battery simulation comprises the anode and cathode regions and the space between the electrodes filled with an electrolyte [60]. 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 × l z hexagonal domain. Figure 2 summarizes the boundary conditions we apply [42]. The initial structure consists of artificial nucleation regions as part of the problem's initial condition. Ellipsoidal protrusions (seeds) with semi-axes r x , r y , r z , and center 0, l 0y , l 0z , exist at the surface of the anode undergoing electrodeposition (see Figure 2); seeding is a widely-used strategy in phasefield simulations of electrodeposition [23,24,28,29], which 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 [42]. The initial condition for a transition zone between solid electrode (ξ = 1) and liquid electrolyte (ξ = 0), vary in the x spatial direction according to: [65]. For the artificial nucleation case, we modify the initial condition formula, replacing "x" by h 0 x rx 2 + y−l0 y ry 2 + z−l0 z rz 2 − 1 within the hyperbolic tangent argument, to account for a smooth transition between the solid seed (lithium metal anode) and the surrounding liquid electrolyte region [42]. We assume the cell's electrolyte to be 1M LiPF 6 EC/DMC 1:1 volume ratio solution [26]. The electrode phase is a pure solid (neglecting any solid phase nanoporosity). Furthermore, Table 2 presents the phase-field model parameters [42]. The normalization constants for length, time, energy and concentration scales are and C 0 = 1 × 10 3 mol/m 3 , respectively [26]. Table 3 provides a list of the simulation settings and numerical parameters used in each numerical test presented in this work. For reference, the numerical tests have been enumerated in the order in which they appear in this paper.
3.1. Symmetric boundary conditions. The symmetric nature of spike-like lithium morphology [46,30,42] allows us to reduce the computational cost by using symmetry boundary conditions to model only one-quarter of the domain. Thus, we split the domain in four, and apply Neumann boundary conditions (∇ξ · n = 0; ∇ ζ + · n = 0; ∇φ · n = 0) to those boundaries facing the center of the domain as depicted in Figure 3. Therefore, we reduce the size of our computational domain to 25% (l x , l y /2, l z /2). We verify our strategy by comparing the previous 3D simulation result using the whole domain (see [42]) and those obtained using symmetric boundary conditions (see Figure 4d). Thus, the symmetric boundary conditions reduce the computational cost, which allows us to use finer meshes in the sensitivity analysis.
4. Phase-field interface thickness to mesh resolution ratio: A sensitivity analysis The phase-field interface thickness significantly affects the simulated reaction rates [41]. Wider interfaces (larger δ P F ) increase the reactive area in the simulation, which induces faster electrodeposition rates. Numerical evidence shows that 1D interface-thickness-independent growth (convergent results) are possible well before reaching the physical nanometer interface width [66,41].
In this section, we perform a sensitivity analysis to study possible mesh-induced effects on the simulated 3D dendrite morphology, propagation rates (dendrite's height vs time), electrodeposition rates (dendrite's volume vs time), and energy levels. We compare 3D simulation results for different spatial resolutions and phase-field interface thicknesses. We use the anisotropy model from [42] to enable the comparison with the results reported therein.
We select a geometric unit that characterizes a real cell structure [28,60]. 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 those that occur in physical scale cells under the same applied voltage, due to the short separation between electrodes l x = 80 [µm] [42]. Figure 4 presents a collection of 3D spike-like lithium dendrite morphologies (isosurface plot of the phase-field variable ξ = 0.5), obtained by varying the phase-field interface thickness (δ P F = 1, 1.5 and 2 [µm]), and mesh sizes (h = 0.5, 0.375 and 0.25 [µm]) (Tests 1 to 6). Thus, we test different combinations of phase-field interface thickness to mesh resolution ratios (R = δ P F /h = 3 to 8). We compare dendrite morphologies at the moment they reach a height of H = 45 [µm]. Figure 4 shows spike-like patterns that exhibit morphological similarity. Each consists of the main trunk and four side branches growing in each horizontal direction (as a result of the body-centered cubic (bcc) crystallographic arrangement of lithium metal [28]). Smaller phase-field interface thicknesses produce more slender dendritic morphologies (cf. Figures 4f (δ P F = 1 [µm]) and 4c (δ P F = 2 [µm])). We use Paraview's mean curvature measurement [67] to analyse the dendrite's tip radius (isosurface plot of the phase-field variable ξ = 0.5). Thus, the measured dendrite's tip radius in Figure 4f In addition, Figure 4 shows that increasing the resolution ratio delivers thicker dendritic morphologies (R = δ P F /h), keeping the phase-field interface thickness constant (more elements at the interface). We compare the morphologies in the first row of Figure 4 (δ P F = 2 [µm]), against those on the second row of Figure 4 (δ P F = 1.5 [µm]); in both cases finer mesh resolutions (higher R) lead to less slender and less branched dendritic morphologies. Figure  5 shows the evolution of the Gibbs free energy [42]. We plot the total energy curve for different simulation set-ups (phase-field interface thickness δ P F and mesh resolution ratio R), as the figure indicates. Figure 5 shows that in all cases, the systems' discrete free energy does not increase with time (adaptivity delivers discrete energy stable results). We obtain a maximum energy difference of about 9% (t = 0.6 [s]) between the simulations with the maximum (δ P F = 2 [µm] & R = 8) and minimum (δ P F = 1 [µm] & R = 4) total energy levels (see Figure 5). In addition, those dendrites sharing a similar level of total energy (represented in green, orange, and purple) exhibit closer morphological resemblance, as the figure shows.  Figure 5. Comparison of energy time series for 3D spike-like dendrite growth simulations, for different phase-field interface thickness to mesh resolution ratios (R = δ P F /h), under φ b = −0.7 [V ] charging potential. Dendrite morphologies at height H = 45 [µm] for reference (colours by phase-field interface thickness and mesh size). Tests 1 to 6. Figure 6 shows the effect of the phase-field interface thickness and mesh resolution ratio on the dendrite's propagation rate (H vs t). Smaller phase-field interface thickness produces significantly faster propagation rates. For example, simulation using smaller interface thickness (δ P F = 1 [µm]) exhibits up-to 100% faster growth rates than those obtained under larger interface thickness (δ P F = 2 [µm]). Figure 6 shows that slower dendrite propagation rates occur as we increase the mesh resolution ratio (R = δ P F /h), keeping the phasefield interface thickness constant. The inset in Figure 6 plots the maximum Li-ion concentration surrounding the dendrite's tips for δ P F = 2 [µm] where the enriched Li-ion concentration decreases as we increase the mesh resolution R = δ P F /h (more accurate solution), leading to slower propagation rates (H vs t).
Time [s] Figure 6. Comparison of 3D spike-like dendrite propagation rate, for different phasefield interface thickness to mesh resolution ratios ( charging potential. Dendrite morphologies at common height H = 45 [µm] for reference (colours by phase-field interface thickness and mesh size). The inset shows maximum Liion concentration as a function of time for different R = δ P F /h ratios, using the same phase-field interface thickness (δ P F = 2 [µm]). Tests 1 to 6. Figure 7 shows the dendrite's volume evolution as a proxy of the overall electrodeposition rate (volume of lithium metal deposited over time). The effect of the phase-field interface thickness and mesh resolution ratio on the overall electrodeposition rate is less significant (percentage-wise) than it is for the dendrite's propagation rate (dendrite's height over time). For example, Figure 7 shows a maximum electrodeposition rate difference of less than 20% (volume vs time) between the fastest and slowest simulation results. Thus, faster dendrite's propagation rates occur for smaller phase-field interface thicknesses due to the lithium metal being deposited / spread over a smaller surface area (more slender dendritic morphologies), rather than differences in the overall electrodeposition rate (minor effect).
This analysis shows (see Figure 7) that for phase field interface thickness 2 [µm] or smaller, the simulated electrodepostion rate (volume of lithium metal deposited over time) is relatively insensitive to the numerical parameters (δ P F and R). On the other hand, the simulated dendrite propagation rate shows stronger numerical dependencies (see Figure 6), affecting the level of realism of our results. Thus, propagation predictions presented here should only be taken as a comparison indicator between numerical tests, as we work towards smaller phase-field interface thickness to increase the accuracy of our simulations.
. Comparison electrodeposited volume for 3D spike-like dendrite growth simulation, for different phase-field interface thickness to mesh resolution ratios ( charging potential. Tests 1 to 6.

3D simulations using modified surface anisotropy
We evaluate the performance of the surface anisotropy representation model for metal anode battery simulations (see Section 2.1). We perform numerical tests to gain insight into the benefits of this modification compared with the results previously obtained in Section 4 and the preceding 3D simulation work [42]. These studies consist of 3D phase-field simulations of lithium dendrite formation during battery charge state to explore three-dimensional highly branched "spike-like" dendritic patterns, commonly observed experimentally. These patterns form under high current density loads, which correspond to fast battery charge [30,45,46]. 5.1. Comparison of simulated patterns: Surface anisotropy model. We study the performance of the modified surface anisotropy representation (4) using a 3D numerical experiment (Test 7) and comparing the resulting morphologies with those obtained for single nucleus simulations using the standard anisotropy representation (see Figure 4d). We use the simulation setup of Section 3. Figure 8 shows the morphological evolution of the simulated lithium dendrite and the enriched lithium-ion concentration ( ζ + > 1) in the vicinity of the dendrite tips, with peak values of ζ + = 2.3. In agreement with previous numerical experiments, the simulation forms a spike-like, symmetric, and branched pattern. The dendrite morphology consists of the main trunk and sets of four equal orthogonal branches developing to the sides. The side branches grow up to 18[µm] long (60% longer than in previous simulations), and 5 to 10[µm] width. Furthermore, the side branches growth is not perpendicular to the main truck but at an angle of about 25 • to 50 • , with a separation of about 4 to 8[µm] between branches. These results show improved morphological similarity with dendritic patterns observed in lithium experiments performed by Tatsuma et al. [46] (within 10% to 20% difference compared to dendrite's measured features) . Figure 9a plots the evolution of the surface energy for the 3D lithium patterns we simulate, revealing equivalent energy levels (less than 4% difference) when compared against the results previously obtained using the initial, non-modified, anisotropy representation. Thus, the numerical experiment demonstrates that the modified anisotropy representation did not significantly affect the surface energy. Additionally, Figure 9b characterizes the morphology by tracking the dendrites' volume-specific area (µm 2 /µm 3 ). We compute the volume-specific area average ratios of 0.83 and 0.78 µm 2 /µm 3 for the single nucleus and modified anisotropy simulations, respectively. The slightly lower surface area/volume ratio of the modified anisotropy representation (-6%) indicates the dendrite growth has fewer branches (fewer but larger branches).   5.2. Mesh orientation effect for different surface anisotropy representations. We further compare the behavior of the standard (Test 9) and modified anisotropy representation behavior (Test 10) by studying the mesh orientation's effect on each simulated pattern. So far, the simulations results use structured meshes aligned with the Cartesian axes. Unlike previous 3D simulations, we now proceed to redistribute the mesh (node's mapping) by performing a 25 • rotation around the x-axis, as depicted in Figure 10.
Thus, here we test the dendrite's sensitivity to the mesh orientation. Figure 11 compares the dendrite morphologies using the standard anisotropy representation, using Cartesian (Test 8), as well as 25 • rotated mesh distribution (Test 9) (see Figure 10)    its center and connecting the two farthest points of the contour (see Figure 12). We compare dendrite morphologies at the moment they reach a height of H = 45 [µm]. The analysis of the dendrite's side branches (horizontal slices) in Figure 12 reveals that the standard anisotropy representation is more sensitive to the orientation of the mesh. For example, orientation analysis in Figure 12a depicts dendrite's rotation angles of around θ = 23 • , evidently aligned with the 25 • of rotation imposed to the mesh. The side branches, due to the modified surface anisotropy representation (4) exhibit significantly smaller rotations of about θ = 4 • , using the same simulation conditions (see Figure 12b). Thus, the modified anisotropy model shows reduced sensitivity with respect to the mesh. Additionally, mesh refinement reduces the simulations' sensitivity to the mesh orientation.

3D Orientation of lithium crystal:
A surface anisotropy-based strategy. Given the random nature of the nucleation process, we need to deal with some degree of randomness and uncertainty when determining the preferred growth direction of the dendrite's crystal in the battery. The orientation of the crystal, determined by the orientation of the surface anisotropy, will direct the preferred direction of growth of the lithium dendrite. Thus, we adapt this well-known crystal growth model for solidification in [52] to electrodeposition dendrite growth. We define a material system of coordinates ( x, y, z), in which each axis corresponds to the 100 direction of a cubic lattice. The following coordinate transformation T is used between the coordinate systems of (x, y, z) and ( x, y, z): where θ x , θ y , and θ z are the rotation angles around the x, y, and z axes, respectively. We use (6) to compute the gradient of the phase-field variable (∇ξ) and use it in the surface anisotropy expression (1). Therefore, we can assign random values to each of the rotation angles (θ x , θ y , θ z ) to control the preferred growth direction of the lithium dendrite and side branches. We test the proposed strategy (Test 10) by applying a θ x = 35 • rotation (θ y = θ z = 0 • ) to the lithium surface anisotropy when using a Cartesian mesh. Figure 13 shows the simulated spike-like lithium dendrite morphology after t = 0.7 [s]. The top-view analysis 13a reveals that this rotation resulted in a dendrite rotation of about θ = 31 • under the applied anisotropy angle.

5.4.
Mesh size effect for different surface anisotropy representations. Following Section 4, we study the spatial sensitivity of the modified surface representation under mesh refinement (Test 12). Given spikelike lithium dendrite symmetry (see Figure 8), we use symmetry condition of Section 3.1 to reduce the computation cost and improve the mesh resolution. Thus, we model only one-quarter of the domain, using a 200 × 100 × 100 tensor-product mesh with a mesh spacing of 0. 25 [µm] in the region of interest (bottom half of the domain). The node mapping produces a finer mesh to properly capture the smaller phase-field interface thickness δ P F = 1 [µm] (4 elements in the interface [41]) and the steepest gradients of ζ + and φ. Figure 14 shows the simulated lithium dendrite (isosurface plot of the phase-field variable ξ = 0.5). The simulation forms a spike-like and highly branched pattern. We calculate the electric field distribution by differentiating the resolved electric potential E = −∇φ. Figure 14 shows how the electric field localizes in the vicinity of the dendrite tip [42]. In agreement with previous numerical experiments [42], the electric field distribution leads to an enriched concentration of the lithium-ion it induces due to the strong migration from the surrounding regions [21]. The dendrite growth does not occur perpendicular to the stack but at an angle of about 22 • between main trunks (cf. Figure 13)). The stacks repel each other, similar to previous multiple-nuclei results [42].  Figure 15 compares the effect of the mesh resolution and phase-field interface thickness on the simulated morphologies, with and without the presence of the modified surface anisotropy term. For the modified surface anisotropy representation, smaller phase-field interface thickness (δ P F ) and finer mesh resolution (h) lead to more branched and detailed dendritic patterns. However, in the standard case, finer mesh resolution leads to less branched micorstructures (see Figure 15c).
Despite the morphological differences mentioned above, the computed lithium electrodeposition average rate in this case (10,800 µm 3 /s ) is within analogous simulation results under coarser mesh resolution (9,150 µm 3 /s , see Figure 9b), as well as simulation result using the standard anisotropy representation, using coarse and fine mesh options (10,100 to 12,400 µm 3 /s , see Figure 7). These results show that using the modified surface anisotropy representation is robust in relation to the rate of electrodeposition (volume of lithium metal deposited over time), showing relatively low sensitivity to numerical parameters of our choice (δ P F and R). In practice, the amount of dendritic lithium directly reduces the Coulombic efficiency of the battery [68]. Therefore, we envisage a future application of our model in evaluating Coulombic efficiency reduction due dendite's formation in rechargeable lithium batteries.

Experimental-scale 3D simulations of lithium dendrite formation
This section evaluates the performance of the modified surface anisotropy model (see Section 2.1) in experimental-scale interelectrode distances. We map the nodal distribution concentrating the nodes in the region of interest, inspired by experimental and simulation results. The increased domain size affects the lithium electrodeposition behavior by increasing the interelectrode distance. We discuss the lithium dendrite propagation rates and morphologies for different charging voltages.
6.1. Meshing strategy for experimental-scale 3D simulations. The high computational cost of detailed 3D simulations of lithium electrodeposition at the whole-cell scale is a well-known challenge of electrodeposition simulations [25,28,42]. A limiting factor is the domain size, which imposes practical restrictions on the 3D simulations. Previously, we chose a domain size of (80 × 80 × 80 µm 3 ) that ensures the simulation volume at an affordable computational cost. But this short domain induces (electrode separation of l x = 80 [µm]) dendrite growth rates that are two orders of magnitude higher than those observed experimentally [42]. A detailed analysis of lithium dendrite experiments reveals that, despite the interelectrode separation distance in experimental cells, which ranges from 1 to 10[mm] [69,70,71], the lithium dendrites effectively occupy up to 20% of the interelectrode space. Thus, we focus on this area of interest, the region/volume of the experimental cell where lithium dendrites develop, near the anode surface. Furthermore, previous simulation results show that the spatial distribution of the variables in the bulk region (outside the area of interest) exhibit either constant values, such as ξ = 0 and ζ + = 1, or small electric potential gradients ∇φ. This weak variation indicates that only a few elements may adequately capture the bulk behavior. At the same time, we assign most computational resources to the area presenting the steepest gradients of ξ, ζ + and φ, representing a small portion of the whole domain.
This section applies the modified anisotropy representation in 3D simulations targeting experimental time and length scales. We describe a simple meshing strategy that exploits the aforementioned distribution by combining uniform node's mapping in the portion of the physical domain where the lithium electrodeposition process occurs (finer and regular mesh), with an exponential increment of the mesh size as we move away from the electrode into the electrolyte's bulk region. Thus, we use a 3D structured mesh with eight-node hexahedral elements. Within the bulk region, in particular in the x-direction x r = 2 j ×x u with j = 1, 2, ..., n; 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 exponential function transitions smoothly by doubling the element size when moving away from the area of interest into the bulk region. This focussed-mesh distribution in the area of interest and subsequent stretching allow us to achieve experimental interelectrode distances in only a few additional elements. Consequently, although the detailed portion of our domain (l xu ) remains the same (80 × 80 × 80 µm 3 ), we are now able to avoid simulations with higher-than-normal dendrite's growth rates, by achieving experimental interelectrode distances (l x up-to 5000 [µm]).
Thus, we select a geometrical unit that characterizes a real cell structure [28,60]. We choose a computational domain of 5000 × 80 × 80 µm 3 . Figure 2 summarizes the boundary conditions we apply. Lateral dimensions remain unchanged in this case (l y = l z = 80 µm 3 ), which along with periodic boundary conditions applied on the lateral faces, generates a 80 [µm] × 80 [µm] nucleation arrangement surrounding the simulated morphology (neighbouring dendrites). The implemented approach constitutes a more realistic alternative than modeling a single isolated dendrite [69,70]. Furthermore, neighboring dendrites act as a barrier (charge repulsion effect) that limits the side development of the simulated electrodeposit beyond the domain's boundaries.
We use a 180 × 100 × 100 tensor-product mesh, with a mesh size of 0.4 [µm] in the region of interest (l xu ). We partition the mesh into eight processors. Figure 16 identifies each core with a different color, showing that the tensor-product mesh can efficiently allocate resources in the region of interest (l xu ).  [42]. We modify the initial condition, by introducing a constant electric potential gradient in the liquid electrolyte region, from φ = φ b at the electrode-electrolyte interface, to φ = 0 at x = l x (cathode), which corresponds to the experimental observations by Nishida et al. [70]. They measured the initiation periods (time transient) for dendrite precursors to start to grow (become visible under an optical microscope) between 4 to 140 s [70]; shorter initiation times occur under larger applied current density. Therefore, sufficiently developed dendrite nuclei may take several seconds to appear, depending on the electrodeposition conditions. This time is sufficient for developing the electric potential gradient in the electrolyte. In addition, the initial conditions for ξ and ζ + remain the same, with a smooth transition between the solid seed (lithium metal anode) and the surrounding liquid electrolyte region [42]. Figure 17 shows the morphological evolution of the simulated lithium dendrite (isosurface plot of the phasefield variable ξ) under φ b = −0.7 [V ] charging potential (Test 13). This setup yields realistic simulation time scales due to the larger interelectrode distance we employ [69]. Stationary propagation rates (dendrite's tip speed) of around 0.2 [µm/s] are reached after 70 s of simulation (see Figure 20a). The simulated growth rates are larger than those reported by Nishikawa et al. [69] in experimental measurements of lithium dendrite growth in 1M LiPF 6 electrolyte (0.06 [µm/s]) due in part to the higher (almost double) applied current density in our model. Our results are within the range of lithium dendrite growth rates reported by Nishida et al. [70] (0.25−0.55 µms −1 ) using a different electrolyte type (LiTFSI). Therefore, further well-controlled experimental studies, with detailed characterization of the system parameters, combined with modeling will be necessary to improve the correlation [72]. Unlike previous 3D simulations, forming spike-like patterns [42], this case yields less branched, blunt tip, finger-like morphologies. The dendrite morphology consists of four main trunks growing from each nucleus, with pairs of orthogonal branches developing to the sides. The observed morphological difference is a consequence of the spatial distribution of the electrostatic potential in the electrolyte (φ). Although the applied electric potential remains the same (φ b = −0.7 [V ]) the larger interelectrode distance results in a significantly different electric field distribution ( E = −∇φ). The electric field surrounding the electrodeposit region can be 60 times smaller than in previous simulations lowering the current density (consistent with the change ratio in the interelectrode distance 5000 [µm] /80 [µm] = 62.5). This weaker current density results in a weaker action of the electric migration forces over the distribution of lithium ions in the electrolyte. Thus, lithium ions are less prone to accumulate around dendrite tips due to the counteracting influence of the concentration diffusion gradient, producing less branched and blunt morphologies.
The lower electric field effect is in agreement with experimental observations by Chae et al. [47], where a variation of the separation between the electrodes revealed a considerable difference in the electrochemical deposition of lithium (experiments under 1 mA/cm −2 applied current density). Chae et al. [47] observed that the lithium deposition behavior and morphology changed from "hazardous" needle-and moss-like dendritic structures to "safer" morphologies (smooth and round shaped surface) as interelectrode spacing increases. The variation of lithium deposition behavior was ascribed to a difference in the Li-ion concentration distribution. Thus, when under shorter interelectrode separation (< 500 [µm]), lithium electrodeposition occurs closer to the high Li-ion concentration regions (formed by the release of Li-ions from the counter electrode), producing a non-uniform directional deposition of lithium. Sharp dendritic structures can grow and penetrate porous separators, which are potentially dangerous as they can create a short battery circuit [73]. On the other hand, larger electrode separations (2000 and 4000 [µm]) lead to a more uniform deposition, without any angular edges or sharp tips, due to lower Li-ion concentration and electric potential gradients [47]. Although the current density applied in the present simulation is lower  than in previous numerical examples, it remains well above the limiting current density of the system, about i lim = 2 mA/cm −2 [74]. For over-limiting current densities applied to the cell, 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 SEI [73]. Figure 18 shows that the fully developed lithium dendrite morphology (t = 210 [s]) resides within the region of interest (well-resolved portion of the domain: ≤ 80 [µm]), which makes-up only 1.6% of the whole domain. Although the system size is similar to previous 3D simulations presented in this work (5,400,000 degrees of freedom), the temporal evolution is slower, increasing the computational time to four times.
Next, we present a 3D phase-field simulation of lithium dendrite formation under more negative applied voltage φ b = −1.4 [V ] (Test 14). We use the setup of the previous experiment, with the sole difference of the applied voltage φ b . We adjust the interfacial mobility parameter L σ to the newly applied electro potential to achieve the right balance between the phase-field interface energy term and the electrochemical reaction contribution (see Table 3) [41]. Figure 19 depicts the evolution of the lithium dendrite (ξ isosurface). As in the previous experimental-scale case, we obtain realistic simulated time scale, with stationary dendrite propagation rates of about 0.4 [µm/s] (see Figure 20b). The higher propagation rate in this case is due to the higher applied current density (from φ b = −0.7 to −1.4 [V ] charging potential), which agrees with experimental results, where higher current densities produce faster electrodeposition and dendrite propagation rates [61,69,75]. Furthermore, we see that computed dendrite propagation rates are within the range of lithium dendrite experiments reported by Nishida et al. [70].
The simulation produces a spike-like, symmetric, and highly branched pattern, with morphological resemblance to previous dendritic deposits obtained under shorter interelectrode distance [42]. The microstructure consists of four main trunks growing from each nucleus, with pairs of orthogonal branches        developing to the sides. Figure 21 depicts the evolution of the simulated spike-like dendritic morphology, together with the spatial distribution of the lithium-ion concentration ζ + in the electrolyte region. Figure 21 shows the evolution of the Li-ion concentration profile extending over 400 [µm] in the stack direction (x); where the deposition process depletes the lithium-ion concentration close to the electrode (shown in blue). This behavior contrasts with smaller-scale simulations presented earlier in this work and in previous work [42], where Li-ion concentration enriches the dendrite tips due to large electric migration forces (see Figure 9). This dendrite-tip enrichment can happen in a close-to-short-circuit condition (short interelectrode distance). Nevertheless, our simulations indicate that lower electro-potential gradients, such as those obtained under experimental-scale interelectrode distances, do not seem to generate high Li-ion concentration around the dendrite tips (competition between electric migration forces and Li-ion concentration diffusion gradient). This observation is in agreement with experimental measurements of Li-ion surface concentration by Nishida et al. [70], where the concentration of Li-ion near the electrode surface was reduced from 1 M (initially) to less than 0.1 M, after a few tens of seconds of electrodeposition, depending on the experiment's conditions.  We obtain apparent morphological differences from the previous dendritic lithium electrodeposition simulation under the experimental-scale domain (compare with Figure 17). Although in this case, the electric field ( E = −∇φ) surrounding the electrodeposit region remains low relative to previous simulations with shorter interelectrode separation ∼ 30 times smaller; the larger charging voltage (φ b = −1.4 [V ]) induces a spike-like and highly branched dendrite (over-limiting current density condition). This result agrees with previous two-dimensional phase-field studies investigating the effect of the applied voltage on the electrodeposit's morphological structure. Increasing the applied voltage produces faster dendrite formation with the tip splitting phenomenon [24], changing from a needle or finger-like structure to a tip splitting or spike-like pattern [29]. The reactive term of the phase-field equation, see (4), is exponentially affected by the electric potential through η a = φ − E Θ . The applied voltage increases the degree of polarization on the electrode, affecting the deposition and accumulation of lithium on the anode surface, which leads to changes in the morphology of lithium dendrites [29]. One verifies this by inspecting the Li-ion concentration gradient ∇ ζ + in the electrolyte region surrounding the dendrites morphologies.  (Figure 22b). Electrolyte regions with higher lithium-ion concentration gradients ( ∇ ζ + > 0.005) are represented with red volumes. Thus, higher lithium-ion concentration gradients appear in the vicinity of the dendrites' tips and side branches in Figure 22b leading to a spike-like, highly branched dendritic lithium (resembling the previously observed electric-migration versus Li-ion diffusion gradient competition happening here at a smaller scale). In contrast, Figure 22a, under lower applied voltage, only presents higher lithium-ion concentration gradients in the vicinity of upper tips of the dendrite triggering vertical and less branched growth. Therefore, the spike-like lithium morphologies forming under over-limiting current density (fast battery charge) [30] can occur either using a large electric field ( E = −∇φ) surrounding the electrodeposit region (close-to-short-circuit condition) or under a large applied voltage φ b (fast battery charge). This strong forcing produces strong electric migration forces, 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.
Following [71], we characterize the morphology by tracking the dendrites' volume-specific area (µm 2 /µm 3 ). Figure 23 compares the growth of the deposited volume versus the surface area for the 3D spike-like lithium pattern we simulate (short interelectrode separation vs experimental-scale results). Despite differences in the time and length scales between these simulations, we obtain similar volumespecific area average ratios; 0.91 and 0.97 µm 2 /µm 3 , for smaller-scale and experimental-scale simulations, respectively. The higher surface area/volume ratio indicates a more branched shape in the experimentalscale simulation. Both cases are within the volume-specific results reported for experimental formation of dendrites in zinc batteries (0.86 and 1.04 µm 2 /µm 3 ) [71] (the literature lacks experimental data for quantitative characterization of the spike-like lithium morphologies).
Remark 1. The similar area/volume average ratios between the dendritic microstructures formed using the experimental-scale simulation domain and the deposition patterns obtained under the short interelectrode distance setup (close-to-short-circuit condition) opens the possibility of using small-scale (lower-cost) 3D simulations. For example, the earlier ones in this work may be a useful testing tool to assess and adjust different 3D strategies before moving into more expensive, well-resolved larger-scale 3D simulations.   Figure 25 shows the behaviour of the time-adaptive scheme, throughout the 150 [s] of the simulation. Starting with a small time-step of ∆t 0 = 10 −8 [s] to initially achieve convergence, followed by an increase in size, until reaching a stationary value of about ∆t n+1 = 0.05 [s] (almost two orders of magnitude larger than previous simulations under smaller interelectrode distance [42]). The weighted truncation error e n+1 (blue) stays close to the minimum tolerance limit (10 −9 ) during the whole simulation. The estimated error does not grow exponentially as in previous cases [42] since the lithium dendrite remains far away from the positive electrode (propagation rate does not accelerate).
Standard discrete approximations do not inherit the a priori nonlinear stability relationship satisfied by phase-field models, expressed as a time-decreasing free-energy functional (see, e.g., [76,77,78] for discussions on energy stable time-marching methods). Therefore, we study the energetic evolution of our system. Figure 26 shows the evolution of the Gibbs free energy of the system Ψ, using our adaptive time integration scheme for the experimental-scale phase-field simulation. 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 26 shows that the total systems' discrete free energy does not increase with time. Thus, we obtain discrete energy stable results in experimental-scale simulations using our secondorder backward-difference (BDF2) time-adaptive marching scheme [42], although the method is not provably stable energetically. Moreover, while the system's surface and chemical energies grow as the area of the lithium deposit increase, the electrostatic energy decreases in time. This behavior, previously observed in smaller scale 3D simulations, is consistent with the electrodeposition process, where the system stores the applied electrostatic energy as electrochemical energy as the battery charges. The inset in Figure 26 shows that the surface energy of the fully developed pattern is almost four times larger than the surface energy computed in [42] for the smaller-scale simulation. The proportionately four-times larger surface area in the experimental-scale case (see Figure 23) explains this scaling.

Conclusions
We use phase-field modeling to investigate the electrodeposition process that forms dendrites within lithium-metal batteries (LMB). These simulations of lithium dendrite formation explain 3D highly branched "spike-like" dendritic morphologies formed under high current density (fast battery charge). These dendritic patterns grow fast across the electrolyte region and penetrate through porous separators, becoming hazardous for battery operation [73]. We analyze the dendrite formation in domains with various sizes using both, short (80 [µm]) and experimental-scale (5000 [µm]) interelectrode separation. The use of symmetry boundary conditions is adequate to exploit the symmetric nature of the spike-like lithium morphologies [46], reducing the computational cost down to 25% of the original requirement.
Through a resolution sensitivity analysis, we asses the mesh-induced effect on the simulated 3D dendrite morphology, propagation rates (dendrite's height vs time), electrodeposition rates (dendrite's volume vs time), and energy levels. The simulated electrodeposition rate (volume of lithium metal deposited over time) is the least sensitive to the numerical parameters of our choice (δ P F and R), while dendrite's propagation rate shows the strongest sensitivity. These results have practical significance since the amount of dendritic lithium produced during charge is directly linked to the reduction in Coulombic efficiency of the battery [68], and the growth rate is related to battery short-circuit predictions [79]. Therefore, future work may evaluate the Coulombic efficiency reduction due to dendrite formation in rechargeable lithium batteries.
We mimic lithium's cubic crystal structure and surface anisotropy by using a 3D four-fold anisotropy model of [49] to simulate crystal growth. We implement a modified 3D representation of the surface anisotropy. This modified model improves the simulation results being less sensitive to the mesh orientation. Furthermore, we introduce a surface anisotropy-based strategy that deals with randomness and uncertainty when determining the preferred growth direction of the dendrite crystal in the battery. Well-resolved simulations showed that the modified model preserves the robustness in the rate of lithium electrodeposition. We test the modified 3D surface anisotropy representation at experimental-scale interelectrode distances (higher computational cost). We simulate two charging voltages (φ b = −0.7 [V ] and φ b = −1.4 [V ]), revealing details about the mechanism behind spike-like dendrite growth at experimental scale. Furthermore, we verify measured morphological parameters, such as simulated dendrite propagation rates, volume-specific area, and side branching rates, within the reported ranges for experimental electrodeposition of spike-or tree-like metal dendrites.
Unlike simulations using shorter interelectrode separation, we observe no enrichment of Li-ion concentration surrounding the dendrite morphology at experimental scale ( ζ + < 1). However, electric migration forces continue to cause lithium cations to move from less concentrated surrounding regions and accumulate around dendrite tips (identified as higher lithium-ion concentration gradients ∇ ζ + > 0.005), triggering spike-growing and highly branched dendritic lithium in the case of φ b = −1.4 [V ] charging potential. In contrast, under 50% lower applied voltage (φ b = −0.7 [V ]), high lithium-ion concentration gradients are only present in the vicinity of the upper tips of the dendrite, triggering vertical and less branched growth, with smoother and rounder surface shapes [47].
Thus, our analysis at the experimental scale confirms what was previously observed under smaller-scale simulations: dendrite formation is connected to the competition between the lithium cation diffusion and electric migration forces, generating an uneven distribution of Li + on the electrode surface [42]. This fact gives insight into inhibition strategies focusing on enhancing the diffusion of lithium ions to achieve a more uniform concentration field on the anode surface, leading to lower dendrite formation propensity [80,81,82,83,84,85,86,87,88,89,90,91,92,93].
Given our understanding of the process, in future work we may add other physical aspects to the simulation; our 3D phase-field model coupled with additional fields will allow us to gain insight into other aspects of dendrite formation and assess some of the proposed strategies for dendrite suppression. Thus, strategies from 2D phase-field models available in the literature could be followed; for example, the current model does not consider heat transfer to simulate the thermal effect during the lithium dendrite growth process. Thermal-induced ion-diffusion may allow us to study dendrite suppression under high operating temperatures [27,40]. Also, the contribution from transport (forced advection) will allow us to study the effect of electrolyte hydrodynamics on the dendrite morphology in flow batteries [31,94], and electrochemical-mechanical phase-field models to study the role of stress in lithium dendrites [30,28]. We will also develop provably unconditionally stable second-order time accurate methods that may deliver larger time-step sizes for phase-field models [76,77,95,96,78], adaptive mesh refinement strategies [97], and improvement of the parallel computation efficiency [97,36].