Effect of diffusion constant on the morphology of dendrite growth in lithium metal batteries

An important problem in lithium batteries is the growth of dendrites or whiskers. The formation of dendrites can reduce battery efficiency and capacity and can also result in a short circuit, causing the battery to catch fire. There has therefore been considerable effort devoted to developing strategies to suppress or prevent dendrite growth.^1–4 The problem is interesting from a fundamental standpoint because the surface does not relax after deposition and, therefore, standard equilibrium theories are not appropriate.⁵ Although a number of strategies have been tried, including modifying the physical properties of electrolytes^6,7 and the metal–electrolyte interface,^8–10 the physics behind the growth mechanism is not well understood. In this work, we investigate the effect of the diffusion constant of lithium and the electric field on dendrite growth using stochastic dynamics simulations of a simple model.

Two strategies to limit dendrites focus on the nucleation at the surface and the dynamics of the metal ions. Suppressing nucleation is difficult because there is no nucleation barrier.¹¹ Furthermore, this has to be achieved for every charging cycle. The other strategy is to control the movement of metal ions during charging by controlling either the electric field or the diffusion of the lithium ions. It has been suggested that the growth of dendrites can be suppressed by increasing the transference number, i.e., the fraction of current carried by the cations.¹² 

In this work, we use computer simulations of a simple model to study the effect of the lithium ion diffusion constant, D, and applied electric potential, V_a, on the growth of dendrites. Only one diffusing lithium ion is modeled explicitly, and all other species, e.g., solvent and other ions, are replaced with a continuum. Following Ref. 13, the lithium ion is modeled as a hard sphere that can stick to the electrode surface or previously deposited ions. The local electric field in the simulations is obtained by solving the Laplace equation, and the system is evolved using stochastic dynamics. The ion is deposited if it is within a specified distance, and the process is repeated with the next ion after a deposition event. Our work is a three-dimensional version of a two-dimensional model investigated by Aryanfar et al.;^13–17 we show that there are significant differences between the two models. We find that the diffusion constant plays an important role in dendrite growth with a low diffusion constant, resulting in more growth. An increase in D also alters the structure of the dendrites from a tall broccoli structure to a flatter cauliflower structure.

We perform stochastic dynamics simulations to study the deposition of particles on a surface. The particles are modeled as hard spheres of diameter r₊ = 1.2 Å (we use Å as the unit of length in this paper). The simulation cell is a rectangular parallelepiped of dimensions 166.7 × 166.7 × 200, with a hard wall at z = 0. At the start of the simulation, a particle is released at a randomly chosen position along x and y, and z = 200. The value z = 200 is large enough that it does not affect the dendrite growth, and all calculations presented have a maximum height that is less than 200. The position of the particle is propagated using the method of Ermak,^18,19

where r(t) is the position of the particle at time t, Δt is the time step, D is the diffusion coefficient, μ is the mobility, g is a displacement vector with random direction and magnitude ∣g∣ = 1, and E is the electric field. We use parameter values suggested by Aryanfar et al.,¹³ and these are collected in Table I.

The position of the particle is propagated until it comes within 0.1 of the surface or an already deposited particle, at which point it sticks and does not move further. After the particle sticks, another particle is released at z = 200. The procedure is then repeated, one particle at a time, until a total of N = 10⁵ particles have been deposited. The entire process is repeated ten times, and the results presented are averages over all trajectories.

The electric field is calculated at the beginning of the simulation and updated after every ten particle trajectories are complete, i.e., ten particles have been deposited. For numerical efficiency, the effect of the diffusing particle on the field is ignored. The electric field is calculated from electrostatic potential, which is obtained by solving the three-dimensional Laplace equation

where Φ is the electrostatic potential.

The Laplace equation is solved by the finite-difference method where the simulation cell is divided into a grid with 100 grid points in each dimension (a total of 10⁶ grid points) and the field at each grid point is calculated. If i, j, and k denote the indices in the x, y, and z directions, respectively, and Φ(i, j, k) is the potential at the grid point, then the boundary conditions are Φ(i, j, 1) = V_a and Φ(i, j, 100) = V_c. We use periodic boundary conditions in the x and y directions, i.e., Φ(1, j, k) = Φ(100, j, k) and Φ(i, 1, k) = Φ(i, 100, k). In finite-difference form, the Laplace equation is

where Δx, Δy, and Δz are grid spacings in the three dimensions. We use an iterative method to solve the Laplace equation. Starting with an initial guess for the grid points (usually a previous case or zero if none is available), we update the values of the potential, starting from k = 2. For the deposited particles, the potential at the nearest grid point is set to zero. The procedure is assumed to be converged if the sum of absolute differences in potential in successive iterations is less than 10⁻⁹. Once the potential has been calculated, the field, E, is obtained from a finite-difference gradient,

The average height $H$ is calculated by dividing the xy plane into 50 × 50 subsections, using the largest z component of particles in the bin as the height of the bin, and then averaging over all bins. The density profile ρ(z) is calculated by binning in the z-direction with a bin size of Δz.

The Minkowski–Bouligand fractal dimension is calculated using the box counting method.²⁰ The simulation cell is discretized into equally sized cubes with edge ɛ, and N(ɛ) is the number of cubes, which contains at least one deposited particle. The fractal dimension F is given by

We choose 1.667 Å ≤ ɛ ≤ 166.7 Å because N(ɛ) displays power law behavior in this (scaling) regime. Figure 1 shows that linear behavior is observed over this range, and we extract F from a fit to the five highest values of ɛ.

The electric field tip effect does not play a significant role in the growth of dendrites. When there are dendrites in the system, the solution of the Laplace equation shows that the field vectors point toward the tips of the dendrite. A previous study on two-dimensional systems suggested that this focusing, referred to as the tip effect, is the major factor in dendrite growth.¹³ We investigate the tip effect by performing simulations where the electric field is not updated during the simulation, i.e., the field is calculated once in the beginning and not updated thereafter. When the electric field is not updated, it is not affected by the growing dendrite except at the dendrite itself. There is, therefore, no focusing of the field vectors near the dendrite and, therefore, no tip effect. Therefore, not updating the field is a way of artificially turning this tip effect off. In our simulations, we find that turning off the tip effect has only a small effect on dendrite growth, unlike the case in two dimensions.

Figure 2(a) depicts the average height of the deposited layer as a function of the number of deposited particles and shows that that the tip effect does not have a qualitative effect on dendrite growth although the growth is slightly faster when the tip effect is turned on. Part (b) shows that the density profile is similar with the tip effect on or off, although, since the dendrites are slightly higher with the tip effect on, the density profile has a correspondingly longer range. This is true for all the values of diffusion constant studied, and two are depicted in the figure. Figure 2 shows that the diffusion constant has a significant effect on the dendrite growth, with a much greater height of the dendrite when the diffusion constant is low. At the lower value of the diffusion constant, there is a plateau in the density profile caused by dendrites of constant density. Although the tip effect is not significant, all the simulations reported in this work include the tip effect since this is a more realistic representation of the system.

The growth of dendrites is weaker with stronger electric potentials or higher diffusion coefficients. Figures 3(a) and 3(b) show the variation of ⟨H⟩ with the diffusion constant and applied potential, respectively. The height decreases roughly linearly with D/D₀ except for the lowest value of V_a. The height decreases in a roughly logarithmic fashion with V_a, indicating a sharp decrease for low values of V_a. The weaker growth with a higher diffusion coefficient can be understood in terms of the ion moving more before it sticks and therefore finding spots closer to the surface. A stronger electric potential drives the ion closer to the surface, resulting in a more dense dendrite layer, which therefore does not grow as high.

The structure of dendrites is not affected by the electric potential. Figure 4 shows the density profiles for two values of the electric potential and various values of the diffusion coefficient. For a given value of the diffusion constant [compare curves in parts (a) and (b) of the same color], the shape of the density profile is not sensitive to V_a. When there is a plateau, the only effect is to increase the length of the plateau, which corresponds to a greater ⟨H⟩. For a given V_a, however, the shape of the density profile is sensitive to the diffusion coefficient.

The particles penetrate closer to the electrode surface when the diffusion constant is high. To estimate the penetration, we define a layer density, ρ_ν(N), as follows: The simulation cell is divided into ten layers in the z-direction [see Fig. 5(a)], and the number density of particles in each layer is tabulated as a function of the number particles deposited. ρ_ν(N) is the number density of particles in layer ν when N particles have been deposited. This is shown in Fig. 5(b). The density in the first layer is always higher when the diffusion coefficient is higher. As a consequence, the ρ_ν(N) has a lower value in higher layers since the total number of particles N is fixed in this comparison. Note also that the layer density plateaus in the lower layers, because these layers are saturated, and additional particles cannot penetrate. The dendrites are clearly not close packed or space filling. This leads to different structures under varying conditions as discussed below.

The number of neighbors for each particle is sensitive to the diffusion coefficient. To quantify the neighbors, we define a coordination number, N_c, as the number of particles between r₊ and 1.5 r₊ of a given particle, averaged over all particles. Both the position of the peak and the width of the distribution depend on the diffusion constant at fixed V_a [Fig. 6(a)], but the distribution is insensitive to the applied potential [panel (b)]. This emphasizes that the local environment of the deposited particles is sensitive to the diffusion coefficient but not the applied potential.

The dendrites have a fractal dimension, F, which varies between 2.72 and 2.85 depending on D/D₀ and V_a. This may be compared to the fractal dimension of a broccoli (F ∼ 2.7) and cauliflower (F ∼ 2.8).²¹ We are not aware of any significance of this fractal dimension although dimensions in this range are common for vegetation.²² The dendrites grow from the substrate upward, such as a tree with a trunk and branches. A dendrite with a high fractal dimension (closer to three) is more space filling and hemispherical, such as a cauliflower, while that with a lower fractal dimension is more conical, such as a broccoli. Figure 7 depicts the fractal dimension as a function of the diffusion coefficient and applied potential and shows that the most conical dendrites occur for low values of D and V_a. For high values of V_a, the fractal dimension is insensitive to the diffusion coefficient, but the variation is more marked for low values of V_a. Similarly, there is a stronger dependence on V_a for low values of D. The broccoli and cauliflower shapes of the dendrites are shown as configurations in Fig. 8.

We present a simple model to study dendrite formation at a surface. The model has two parameters, the applied potential and diffusion coefficient, which we vary to study their impact on the structure and growth of dendrites. We find that the diffusion coefficient is a far more important parameter than the potential and strong dendrite growth is seen when the diffusion coefficient is low. This is contrary to conclusions drawn in a two-dimensional version of the model where the applied potential was found to be more important.¹³ A suggested strategy to limit dendrite growth batteries is therefore to increase the Li⁺ diffusion coefficient.

We speculate that when the diffusion constant is high, the lithium ion can avoid sticking to a growing dendrite and penetrate closer to the surface. This is also true when the applied voltage is high. When both are low, there is a smaller electrical driving force pulling the ion to the surface, and the low diffusion constant makes it more likely for the ion to find a tip of a growing dendrite. This could also explain the difference in behavior between two-dimensional and three-dimensional systems. In two dimensions, the free space becomes non-percolating at a low area fraction of the dendrite, and this prevents the diffusion of new ions to the surface even when the diffusion constant is high. In three dimensions, the free space is never close to the percolation threshold for diffusion. Since the ions cannot penetrate the dendrite in two dimensions, the tip effect plays a significant role in driving the ions to growing dendrites. The physics is therefore different in two vs three dimensions.

We find a non-linear dependence of both average height $H$ and fractal dimension F on the diffusion coefficient and applied potential. If one of these variables is sufficiently high, dendrite growth is suppressed. However, they contribute in a synergistic fashion to dendrite growth when both are low, and the structure of the dendrites is different depending on these parameters. Under low diffusion and applied potential conditions, the dendrites are thinner and grow in a cone-like fashion with a fractal dimension similar to broccoli. This prevents penetration of new particles into the growing structure. Under high diffusion and applied potential conditions, the dendrites are more hemispherical with a fractal dimension similar to cauliflower because particles have penetrated into the growing structure.

The model presented here has only two physical parameters relevant to a real system. It would be interesting to test these conclusions with more realistic atomistic models of battery electrolytes and electrodes. Another direction would be to investigate the effect of surface roughness on dendrite growth. A combination of models at multiple scales would help further our understanding of battery electrolytes and perhaps lead to design principles important to applications.

This study was supported by the U.S. Department of Energy, Basic Energy Sciences (Contract No. DE-SC0017877).

The data that support the findings of this study are available from the corresponding author upon reasonable request.