An atomistically informed mesoscale model for growth and coarsening during discharge in lithium-oxygen batteries

Lithium-oxygen batteries have recently received much interest due to their potential for high energy density compared to conventional lithium-ion batteries.¹ Li-O₂ batteries are fundamentally different from lithium-ion batteries because they use chemical transformations that result in growth of a new bulk phase on the cathode as opposed to intercalation of lithium into the metal oxide cathode in the case of lithium-ion batteries. There are many challenges remaining for the development of Li-O₂ batteries, such as poor cycle life and large charge overpotentials, which require fundamental understanding of the discharge and charge processes. One key aspect of the discharge processes that is still poorly understood is the resulting morphology (e.g., toroid and film growth) that plays an important role in determining the charging behavior of the Li-O₂ battery. The discharge product in Li-O₂ batteries is generally composed of lithium peroxide, Li₂O₂,^1,2 although in some cases it may be a mixture of lithium peroxide and lithium superoxide, LiO₂.^3–6 

There have been various different mechanisms proposed for the growth of the discharge product in Li-O₂ batteries. Some researchers have postulated a surface-growth mechanism with the basic reaction step involving two-electron oxygen reduction and addition of two Li cations on the cathode surface.^7,8 Others have postulated that after oxygen reduction on the surface, subsequent reactions occur in the solution phase involving Li₂O₂ or LiO₂ followed by nucleation and growth of one of these species on the cathode surface.^9–12 Nazar et al. reported a Li₂O₂ formation mechanism by the disproportionation of LiO₂ in solution,¹³ which can account for the lithium peroxide nucleation and growth from solution. Bruce et al. also reported evidence for such a mechanism.¹⁴ A through-solution mechanism has also been proposed involving heterogeneous nucleation from a supersaturated solution of LiO₂.^12,6 The LiO₂ mechanism is supported by experimental evidence that LiO₂ is soluble in the electrolyte prior to deposition on the surface based on quartz microbalance studies.¹⁵ There is experimental and computational evidence that in some cases, there can be slow disproportionation, which results in a discharge product of a mixture of lithium peroxide and lithium superoxide.^5,16 The factors that control the morphology/microstructure of the discharge products are difficult to determine experimentally. Most of the computational studies have focused on the electronic structure of the discharge products, with only a few on mesoscale or continuum level modeling.^17,18

In this paper, we present an atomistically informed mesoscale model for deposition of the discharge product from the electrolyte in a lithium-air battery and coarsening between particles. Due to the fact that there are various different reaction mechanisms that may occur in a Li − O₂ battery depending on cell configuration, electrocatalyst, electrolyte, carbon substrate, etc., we have chosen one specific system to investigate that has been well studied experimentally as well as computationally at the electronic structure level.^3–6 This system is based on activated carbon as a cathode material and tetraglyme (TEGDME) as an electrolyte. Systematic experimental studies have revealed that the disproportionation of LiO₂ occurs after its deposition as a bulk material on the surface so that the discharge product is a mixture of LiO₂ and Li₂O₂. The model we develop addresses only the LiO₂ deposition and not its disproportionation, which will be reported in a future paper. Modeling the LiO₂ nucleation and growth process can provide important fundamental understanding of the growth mechanism occurring for activated carbon in this system as well for other Li-O₂ cathode systems involving deposition from solution.

The paper is organized as follows. In Sec. II, we give an overview of the reaction and deposition mechanism, its formulation as a set of reaction-diffusion and phase-field equations, a multiscale coupling scheme, and a sequence of steps for determining the particle-size evolution. Section III presents the pertinent input material properties, and some of the Density Functional Theory (DFT) calculations used to obtain them. Section IV discusses the numerical implementation of the phase-field model and its validation against the Gibbs-Thompson equation. In Sec. V, we apply the model to a particular set of experiments. Some implications of the model for other scenarios are discussed in Sec. VI followed in Sec. VII by a summary of our main conclusions.

The model is developed as follows: First, the proposed reaction and deposition mechanism is described. Second, we translate the conceptual model into a set of phase-field and reaction-diffusion equations. These equations are then applied to two coupled regimes distinguished by length scale. Finally, a scheme for the sequence of particle evolution including nucleation, growth, and coarsening is described.

A schematic of a mechanism for formation of LiO₂ in solution and deposition on a solid cathode is shown in Figure 1, with reactions numbered and written out in Table I. This mechanism involves reaction and diffusion in the electrolyte. The Li⁺ ions enter the system from the anode, which is not explicitly included in our model. Molecular O₂ enters the electrolyte through the porous cathode and as a reaction product from disproportionation, but is reduced at active Oxygen Reduction Reaction (ORR) sites on the cathode surface to form O₂⁻ anions. The Li⁺ and O₂⁻ ions diffuse in the electrolyte and react to form LiO₂ in the solution. This reaction is thermodynamically favored and should be a fast reaction.¹² The process continues until LiO₂ reaches saturation and deposits as a LiO_2(solid) phase on the surface.

The LiO_2(solid) phase is thermodynamically unstable with respect to dispropotionation to Li₂O_2(solid) and release of O₂. In a Li-O₂ cell, this O₂ will be released into the electrolyte, which can then be reduced again according to reaction 1.⁵ 

The problem is naturally divided into formation and deposition of solid LiO₂, and the transformation from solid LiO₂ to Li₂O₂ and O₂. The current work focuses on the first part, i.e., the production of LiO₂ particles as described by reactions 1-3. The solid state disproportionation reaction is addressed in a separate paper.¹⁹ This work therefore studies a possible process of product formation and particle coevolution during discharge. The disproportionation process is neglected as a first approximation to this system model in part due to the suggested slow reaction process of halflife ∼24 h,^3,5,6,16 comparable to the duration of the experiments. As will be shown below, this simple model explains the consistently observed trend that increasing the discharge current produces smaller, more abundant particles.²⁰ 

In the model developed in this work, the solid state particles are predicted to be spherical, although the experimentally observed shape of the discharge products are often toroidal. Other models exist which propose mechanisms of toroidal formation;¹⁷ however, no explicit driving force for a non-spherical shape is present in this work. The spherical particles in this work are considered representative for a torus of the same outer radius.

The model described above is a reaction-diffusion type for the production of LiO₂ in the electrolyte described by reactions 1 and 2. Reaction 3 describes the formation of a condensed phase through deposition, which can be well represented with a phase-field model. The evolution of the system can then be predicted based on thermodynamic self-consistency through consideration of the free-energy functional.

The numerical model implementation follows that of the work of Welland et al.²¹ in which species are assumed to be present everywhere in the domain, but to vanishing small quantities, on the order of 10⁻²³ M, where they would otherwise be considered “insoluble.” Chemical potentials of species are constant within the interfacial region, and concentrations vary smoothly between their equilibrium values in either phases.²² All species are solved for explicitly, and the local hydrostatic pressure, P, is used as a means to enforce constant total concentration.

The local phase, either liquid electrolyte or LiO_2(solid), is modeled by an Allen-Cahn type phase-field model.²³ The local phase is represented by an independent thermodynamic state variable ϕ, which varies smoothly in the range [0,1]. Pure solid LiO₂ and electrolyte phases correspond to 0 and 1, respectively, with a diffuse interface existing between the two. The local mole fraction of the electrolyte phase is calculated by the function²⁴ 

This model is applied to the case of dimethoxyethane (DME), which is representative of the TEGDME ether used in the experiments on activated carbon with which these results are compared. This formulation is general enough to be also applicable to other electrolytes with a suitable change in material properties.

Following the Allen-Cahn phase-field methodology, the interface energy between the condensed LiO₂ phase and the electrolyte, γ, is introduced into the Gibbs energy functional through the addition of the term $γ3dϕ21−ϕ2+3d∇ϕ2$⁠.²¹ The parameter d represents the thickness of the diffuse interface and is chosen to balance a realistic approximation of the physical interface width and the computational expense required to adequately resolve the interface.

Under conditions of constant temperature and ambient pressure, the appropriate Gibbs free-energy functional describing the system is²¹ 

where c_i and μ_i are the concentration and chemical potential of species i, respectively.

The species in solution, with concentrations c_i (where i = Li⁺, O₂⁻, or LiO₂; see Fig. 1 and Table I), are subject to diffusion and reaction. For each species, the following partial differential equations must be satisfied:

where υ_i and R₂ are the stoichiometric constants and reaction rate for reaction 2 in Table I. J_i is the diffusive flux of species i driven by spatial gradients in the chemical potential,

where D_i is the diffusion coefficient, R_g is the ideal gas constant, and T = 298.15 K is the temperature which is assumed to be room temperature in the course of the current work. The total concentration of all species is c = ∑_ic_i.

For this model, we consider the electrolyte to be an ideal solution, which is reasonable for small concentrations of dissolved Li⁺, O₂⁻, and LiO₂. In the LiO₂ particle, Li⁺, O₂⁻, and DME are insoluble and are represented by vanishingly small concentrations. The chemical potential is then given by²¹ 

where $μi,(amorph)0$ is the reference chemical potential in the solid phase, x_i = c_i/c is the mole fraction of species i, V_i is the specific volume, $qi=expμi,(amorph)0/μi,(elec)0$ is the partition coefficient, and $μi,(elec)0$ is the reference chemical potential of the electrolyte.

As discussed by Welland et al.,²¹ the local hydrostatic gauge pressure is introduced as a means to force the total concentration of all species to remain equal to the equilibrium site concentration, c_sites. In linearized form about P = 0, it is calculated as

where κ_T is the bulk modulus. As discussed in Sec. IV, interface energy increases the hydrostatic pressure of small particles, which implies c > c_sites.

The rate of reaction 2 is $R2=k2μLiO2−μLi+−μO2−$⁠, where k₂ is the associated rate constant and is equal for forward and backward reactions. Mass transport is therefore described by

The deposition (or solidification) reaction 3 can be represented by the phase-field equation and included in the above system through the chemical potentials. The phase evolution equation is

where k₃ is the reaction kinetic term for deposition in Table I. Deposition then proceeds according to the equilibrium in Table I, as captured by the term ∑_i c_i∂μ_i/∂ϕ. The second term, which comes from the surface energy in Eq. (2), drives the phase to minimize its surface to volume ratio, producing a spherical geometry. In using a single kinetic term with two driving forces, we imply that the kinetics of both processes is the same. While this may not be true in practice, it is a useful first approximation.

Although subject to a single set of governing equations, the system evolves on two distinct length scales. As schematically illustrated in Fig. 1, microstructure evolution occurs on the scale of hundreds of nanometers on the cathode surface; by contrast, reaction-diffusion occurs on the scale of millimeters in the electrolyte between the cathode and anode. It would be computationally prohibitive to simulate the microstructure over the extent of the entire electrolyte. Therefore, the scheme in Fig. 1 is divided into two computational regimes, one for the microstructure on the cathode surface and one for the electrolyte above that surface as illustrated in Fig. 2. The electrolyte regime considers Eqs. (7a)–(7d) for O₂⁻, Li⁺, LiO₂, and DME reacting and diffusing in one dimension normal to the cathode surface. The microstructure regime only considers the phase-field evolution and interdiffusion of LiO₂ and DME via Eqs. (7c), (7d), and (8). The microstructure regime is simulated in either 1 or 3 dimensions as required for the scenario of interest.

The two regimes share a common boundary: the cathode surface in the electrolyte and the top of the microstructure regime. At this boundary, the two regimes are coupled strongly by the concentration of LiO₂. In the electrolyte, solvated LiO₂ is produced through the reaction of O₂⁻ produced at the cathode surface, the diffusion of the anion away from the surface and the reaction with Li⁺ from the electrolyte; the solvated LiO₂ subsequently diffuses away from the surface. Concurrently, microstructure evolution is governed by the local concentration of LiO₂ at the cathode, and any growth or shrinkage produces a sink of LiO₂ in the electrolyte. The concentration of LiO₂ at the cathode surface, and its rate of removal from solution is the means by which the two regimes are coupled.

Since in reality a large number of particles is being produced, for which simulation would be impractical, the microstructure simulations effectively consider only one nucleation site, although in some cases with multiple particles. The change in particle size in one microstructure simulation is therefore scaled by the nucleation site density, n_s, and a surface scaling parameter, s, which accounts for the difference in area between the total surface area of the powdered cathode and the simulated flat surface. The surface-scaling factor is fit experimentally as described in Sec. V A 1.

Computationally, this coupling between the two regimes is implemented with a Dirichlet boundary condition for c_LiO2 at the top surface in the microstructure simulations and an outward flux (sink) of LiO₂, $JLiO2s$⁠, at the cathode surface for the electrolyte reaction-diffusion model. The iterative scheme is as follows: given a guess for $JLiO2s$⁠, reaction-diffusion in the electrolyte is simulated and a new value of c_LiO2 at the cathode is determined. This value is used in the microstructure model to evolve the particle size and shape. The new guess of $JLiO2s$ is calculated from the growth rate of those particles by

where V_M is the volume of the microstructure regime. Iterations continue until consistent values of $JLiO2s$ and c_LiO2 are found. Instead of using Eq. (9) directly to calculate the next guess of $JLiO2s$⁠, a Newton method is employed which increases efficiency and robustness of the iterative process. Iterations are halted when the sink varies less than 0.1% between iterations, which was found to typically occur after 3-4 iterations.

The following three-stage scheme is proposed to model the nucleation, growth, and coarsening process for an initially pristine cathode, as shown in Fig. 3. As an approximate treatment of nucleation, we specify the critical nucleus size and use this to determine the nucleation supersaturation threshold.

• Stage 1: LiO₂ accumulates in the electrolyte until the nucleation supersaturation threshold is reached and nucleate on the cathode, indicated by point A.

• Stage 2: During stage 2, the nuclei past the critical size start to grow and new nuclei are still able to form. As the particles grow, they deplete the electrolyte of LiO₂, leading to a maximum in the concentration at point B. Eventually, the concentration drops below the nucleation threshold and formation of the new nuclei will cease at point C.

• Stage 3: As the concentration falls, small particles become unstable and dissolve while larger particles continue growing, referred to as coarsening. The concentration of LiO₂ in the electrolyte is being partially replenished by continued discharge and the dissolution of small particles. Ultimately, particles will grow large enough to not be affected by surface energy, and thus can be considered bulk phases.

An additional point D in stage 3 may be considered as a reference point for our later analysis of experimental data. Given the time, maximum particle radius, and minimum stable particle radius at point D, the attachment kinetics and surface-scaling factor may be fit, as will be demonstrated below (see Sec. V).

Few material properties for this system are experimentally available due to LiO₂ being thermodynamically unstable to disproportionation. Here, DFT calculations of material properties will be employed in these cases, which also provides insight into important mechanisms. Alternately, approximations based on similar measured properties will be used.

The diffusion coefficients and solubilities of Li⁺, O₂⁻, LiO₂, and DME are summarized in Table II. For Li⁺ and O₂, material properties are obtained from Lu et al.²⁵ for the case of 0.1 M Li⁺ in DME. The diffusion coefficient of O₂⁻ is assumed equal to that of O₂. As a first approximation, the diffusion coefficients D_LiO2 and D_DME are assumed equal to that of Li⁺. Moreover, the diffusion coefficients are assumed to be the same in both the electrolyte and the solid phase; however, the latter is inconsequential since Li⁺, O₂⁻, and DME are only present in vanishingly small concentrations.

In our simulations, the densities of both phases are held at the value of LiO₂. This simplification avoids the need to consider flow of the electrolyte upon growth or shrinkage of the LiO₂ particle and the corresponding advective currents for species transport, greatly reducing the computational complexity.

The rate of production of LiO₂ in the electrolyte from the applied discharge current is governed by the rate constant k₂ and by the equilibrium concentrations of Li⁺, O₂⁻, and LiO₂ in the electrolyte. It is thought that this reaction proceeds very quickly and quickly consumes all the O₂⁻.¹² Therefore, we assign k₂ = 10⁴ mol s/m⁻³ and $CO2−DME=10−10$⁠, which is sufficient during the simulations to immediately react O₂⁻ with Li⁺.

Table III summarizes other material properties, including those derived from DFT calculations. The nucleation-site density is an estimate derived from counting the number of particles observed after 300 mA h/g of discharge in a Li-O₂ cell based on an activated-carbon cathode;³ this value may differ for other systems.

To compute the density and bulk modulus of bulk amorphous LiO₂, the total energy vs the total equilibrium simulation-cell volume, E(V), shown in Fig. 4 is used to fit the equation-of-state (EOS) of the amorphous LiO₂ structure. The static E(V) dataset includes the self-consistent total energies calculated using the VASP DFT code²⁶ with periodic boundary conditions for different volumes with isotropic volume changes while keeping the atomic positions fixed. The calculations were spin-polarized and carried out using the gradient-corrected exchange-correlation functional of Perdew, Burke, and Ernzerhof (PBE) with plane-wave basis sets up to a kinetic energy cutoff of 400 eV. The PAW method was used to represent the interaction between the core and valence electrons, with the Kohn-Sham valence states (i.e., 2s for Li, 2s 2p for O) expanded into plane-wave basis sets. The convergence criterion of the total energy was set to be within 1 × 10⁻⁵ eV for the K-point integration, and all the atoms and geometries were optimized until the residual forces became less than 1 × 10⁻² eV/Å within a single K-point grid (i.e., Γ-point). The initial configuration of the DFT-predicted LiO₂ bulk amorphous structure is obtained by 54 randomly packed LiO₂ molecules inside a cubic simulation cell before undergoing the full geometry optimization (including both the atomic coordinates and cell lattices). From the relaxed geometry, the atomic density of the amorphous LiO₂ structure is found to be ∼2.19 g/cm³, slightly smaller than the predicted atomic density of crystalline Li₂O₂ crystal of ∼2.37 g/cm³ reported in our previous study.²⁷ By using the quasi-harmonic Debye model,²⁸ the computed E(V) data were fit to the third-order Birch-Murnagham EOS. This yields the calculated bulk modulus κ_T for use in Eq. (6) to be ∼85.82 GPa.

The solubility of LiO₂ in the electrolyte in equilibrium with a flat interface (i.e., radius = ∞), $cLiO2∞$⁠, is one of the two factors which mediates mass exchange between particles, which must be included in the mesoscale model of growth and ripening. Unfortunately, there is no reported solubility for LiO₂ in any electrolyte, although there is experimental evidence for its solubility from quartz microbalance studies.¹⁵ Thus, the only available solubility information for LiO₂ comes from computational studies. In a separate study,³⁰ we have calculated solvation free energies of molecular LiO₂ in a DME solvent using various explicit and implicit solvent models, as well as ab initio molecular dynamics. Best estimates for the LiO₂ solvation energy from these calculations along with the calculated lattice energy of LiO₂ were used to determine the solubility of bulk LiO₂. The best estimate for the solubility of LiO₂ is about 1 mM, although this has very large uncertainties because small differences in the calculated solvation energy lead to large differences in the solubility. For example, a difference of 2 kcal/mol in the solvation energy translates to more than 2 orders of magnitude difference in the solubility. Unfortunately, the solvation energies are difficult to calculate even to 2 kcal/mol in accuracy. In practice, we initially use this value for the solubility of LiO₂ obtained from our DFT calculations in our mesoscale model. Subsequently, a value for the solubility will be derived from fits to reproduce experimental data (see Sec. V); however, the solubility from these fits turns out to be similar to that predicted from the DFT calculations.

To estimate the thermodynamic stability of the interface between amorphous LiO₂ and the electrolyte, we used data from previously published ab initio molecular dynamics (AIMD) simulations.⁶ Based on the DME electrolyte (i.e., ∼0.5M LiCF₃SO₃/DME) that has an interface with the LiO₂ amorphous slab having a thickness of ∼1.2 nm (i.e., 525 atoms in a simulation cell), we optimized the amorphous LiO₂/electrolyte interface and found it to be stable with an interfacial energy of approximately −64.5 meV/Ų.⁶ From the AIMD simulation, the solvated LiO₂ amorphous surface is found to be stabilized by DME and LiCF₃SO₃ molecules at the interfaces.

The calculations showed significant charge transfer and strong adsorption of DME onto the LiO₂ surface. Surface adsorption can produce negative surface energy without destabilizing it;³¹ however, the current phase-field model does not include surface adsorption, and a negative surface energy would destabilize the surface as seen in Eq. (8). The problem is that DFT is more physically comprehensive in that multi-physics phenomena such as surface adsorption are included naturally, whereas the current phase-field model is derived in such a way that these phenomena are not included.

As a pragmatic solution, we proceed with the phase-field simulation based on an approximate positive surface energy corresponding to that of Li₂O₂, [0001] of 1.136 J/m².⁷ The effects of this uncertainty will be discussed in Sec. V.

The interface attachment kinetic factor k₃ plays a critical role in the current work due to its impact on the coarsening process. Ripening between particles of different sizes can broadly be categorized as being limited by either speed of diffusion of matter between particles or the kinetics of attachment of solvated mass onto the particle. In the case of diffusion-controlled Ostwald ripening in the absence of growth, Lifshitz-Slyozov-Wagner^32,33 theory predicts that the rate of change of the average radius $R$ is given by

Given the above material properties, a particle with a radius of 100 nm would grow at a rate of almost 15.5 μm/h. Since a particle with a radius of ∼125 nm was seen to exist at 4.8 h in an activated-carbon Li-O₂ cell,³ this is clearly not the case even considering large uncertainty in $cLiO2∞$⁠.³ Therefore, we consider that k₃ is relatively slow, and coarsening is attachment-kinetics controlled.

The low value of k₃ may be related to both the solvation shell around LiO₂ in solution and the binding of DME molecules on the LiO₂ surface. Deposition of LiO₂ involves removal of the solvation shell and penetration of the adsorbed layer. This parameter is fit to reproduce experimental observations as described in Sec. V.

The model is implemented using the FEniCS package,^34–38 an open source Finite Element Method package, and configured to use the PETSc^39,40 library for the linear solver backend. Quadratic Lagrange elements were used at a resolution of d/3.5. An adaptive time stepping algorithm was used with an absolute tolerance of 0.05 on the maximum change in solution vector. The systems of nonlinear equations were solved using Newton’s method with a GMRES Krylov solver. The Additive Schwarz preconditioner with overlap of 2 was used with the Incomplete Lower Upper factorization solver with 2 levels of fill solver on each core.

The microstructure phase-field model described in Sec. III A can be validated against the classical Gibbs-Thompson theory before being applied to more complex scenarios. The effect of positive surface energy is to compress the particle inwards mechanically, modifying its thermodynamic state from that of a flat surface, i.e., an infinite radius of curvature.⁴¹ From the Gibbs-Thompson theory, it is known that the internal pressure of the particle of radius R differs from the hydrostatic electrolyte pressure far from the interface resulting in a change to the chemical potential,

where μ^R is the chemical potential between the inside of the particle, σ is the surface energy, and μ^∞ is far from the interface. Assuming an ideal solution, this implies

A range of radii was simulated with a 1D, spherically symmetric implementation of Eqs. (7c), (7d), and (8) with d = 10 nm. Fig. 5 compares the simulation results with the classical Gibbs-Thompson effect in Eq. (11) and shows excellent agreement. Fig. 6 shows simulated concentration profiles for a 199 nm sphere, indicating the pressure change inside the particle, and the constituent variations across the interface.

We now consider application of the phase-field model described in Secs. III and IV to results reported for a Li − O₂ battery based on activated carbon³ and some additional experiments that have been carried out for this paper to help distinguish growth from coarsening.

In the Li-O₂ battery considered in the experiments powdered activated carbon was loaded onto a 1 cm diameter disk to a density of approximately 1.6 mg/cm². The density of activated carbon is 2 mg/cm³, which implies a layer approximately 8 μm thick depending on particle packing. The activated carbon has surface area of approximately 3840 m²/g which is mostly attributed to nanopores, but includes microscopic surface features and the surface area of a range of particle sizes.³ The disk has 6 drops of DME placed on top, which corresponds to ∼3.8 mm thick layer of DME between the cathode and anode using the properties in Table III. A constant discharge current of 0.1 mA h/g is applied during discharge.

Examination of particle spacing in SEM images³ shows the nucleation site density to be approximately 12 sites/μm². We assume a critical nucleus radius for our simplified nucleation model of 10 nm, which seems physically reasonable.

In the reported results for the activated carbon Li-O₂ battery,³ discharge was stopped at several capacities: 150, 300, 500, 750, and 1000 mA h/g corresponding to 2.4, 4.8, 8, 12, and 16 h, respectively. We first consider experiments stopped at or before 500 mA h/g in the context of the scheme described in Sec. II D. Longer discharge times are considered later.

Initially, the cathode surface is free of LiO₂ and the electrolyte is saturated in Li⁺ but without dissolved $O2−$ or LiO₂. As discharge proceeds, the concentration of LiO₂ in the electrolyte at the cathode follows the simulation results presented in Fig. 7.

During stage 1, the discharge current introduces $O2−$ into the system at the cathode surface, which immediately reacts with Li⁺ and forms LiO₂. The concentration of LiO₂ at the cathode surface rises until the nucleation supersaturation level, $cLiO210 nm$⁠, is reached. We consider that the cusp in the discharge voltage profile corresponds to nucleation, which occurs at ∼0.5 h.

Since no particles have formed yet, only the reaction-diffusion model in the electrolyte needs to be simulated. The cathode surface concentrations are shown in Fig. 7, and the concentration profiles at the end of stage 1 are shown in Fig. 8. At 0.46 h (point A), the simulated concentration of LiO₂ is 13.4 mM, which corresponds to $cLiO2∞=2.6mM$ from Eq. (12). This calculation compares well with the value calculated by DFT of 1 mM.

Calculating the bulk solubility in this manner depends on the surface energy via Eq. (12). Due to the uncertainty in σ discussed in Sec. III, it is useful to quantify the magnitude of this effect. We calculate that σ = 1.136% ± 10% implies $cLiO2∞=2.62%±16%mM$⁠.

The assumption of a 10 nm critical particle size affects the calculated bulk solubility quite strongly due to the exponent in Eq. (12). However, a more complex nucleation model is beyond the scope of this work. Nonetheless, good agreement with the DFT calculation supports this value as a pragmatic model and allows us to proceed.

During stage 2, LiO₂ particles nucleate on the cathode surface and become sinks for solvated LiO₂. Since the growth rate is proportional to the surface area, the largest LiO₂ particles, i.e., the first to have nucleated, will contribute most to the LiO₂ sink. Therefore, this stage is simulated with the microstructure being represented by a spherically symmetric model of one of the largest particles beginning with radius 10 nm at 0.46 h. The interface width d is set to 1 nm as a reasonably realistic estimate of the interface width.

Points B and C occur at ∼0.8 h and ∼1.1 h, respectively, as shown in Fig. 7. The concentration profile in the electrolyte at 0.79 is shown in Fig. 8. The radius of the simulated particle is shown in Fig. 9.

In stage 3, as the concentration of LiO₂ in the electrolyte drops, small particles start to dissolve and add their mass into the electrolyte. This mass contribution is relatively small due to small particle size and simulation of these particles dissolving would be computationally prohibitive. Stage 3 is therefore subdivided by the reference point D corresponding to a SEM image for the activated carbon results³ at 4.8 h (300 mA h/g). This image was chosen since there is an observable decrease in the density of particles on this SEM image from the previous image. Since 50 nm is the minimum particle radius we can reliably distinguish from surface features on these images, we consider this point to be the time at which 50 nm particles dissolve. Following point D, particles that are of an appreciable size dissolve and contribute significantly to coarsening. This coarsening is captured in stage 3b with four particles.

A single spherical particle is simulated continuing directly from stage 2, while disregarding the coarsening resulting from small particle dissolution. At point D, we consider 50 nm particles to begin to dissolve, which implies $cLiO2=cLiO250 nm$ at the cathode surface. The largest particles are ∼150 nm as seen in the inset of Fig. 7. Given the size of the largest nanoparticles at 4.8 h and the expected concentration $cLiO250 nm$⁠, we may fit the attachment kinetic term, k₃, and surface scaling parameter, s, to 2 × 10⁻¹¹ m³/(Js) and 32, respectively.

The fit value of the surface-scaling factor is reasonable considering the thickness of the activated carbon layer and the approximate carbon particle size (∼1–10 μm) in these experiments.³ The current work assumes this ratio to remain the same for all samples, but other experiments may employ different particle size distribution or structures.⁴² 

After point D in Fig. 7, the supersaturation in the electrolyte is no longer sufficient to sustain particles of radius less than 50 nm. Such particles therefore dissolve and contribute their mass back into the electrolyte. To simulate this, the microstructure model is solved in 3D with differently sized particles. A periodic array of particles on the surface is considered with different radii as shown in Fig. 10. This approximation reduces the required computational problem to four particles, for which only 1/8th of the particle is required to be simulated.

For continuity with stage 3a, and to study a suitable particle size range, particles of initial radii 150, 100, 50, and 20 nm are simulated. The results at the beginning and end of stage 3b, 4.8 and 8 h, respectively, are shown in Fig. 11 where the volume for ϕ > 0.5 has been made translucent. An animation showing all the time steps is Fig. S1.⁴³ The concentration at the surface is shown in Fig. 7 and shows a slight perturbation as the 20 nm particle dissolves. The concentration profile in the electrolyte at the end of stage 3b is shown in Fig. 8.

The concentration of LiO₂ along the x-axis in Fig. 11, between the 150 and 50 nm particles, is shown in Fig. 12 with the region between particles enlarged in the inset. In the inset, the gradient between particles can be seen implying direct communication between the particles, as opposed to a mean-field approximation.

There is a discrepancy in particle sizes at the end of stage 3b (8 h), where the experiment for the activated carbon³ shows particles of almost 250 nm, whereas the simulation only produces particles slightly less than 200 nm. This may be attributed to having only considered four particles in the simulation of stage 3b, a choice partly determined because of the large computational expense of simulating more. Only one of the four particles completely disappears during the simulated process, whereas SEM images show disappearance of almost 80% of the particles.³ The additional particles missing from the simulation would contribute their matter to the larger particles too, which would result in larger radii in Fig. 9.

As part of this study, we carried out additional experiments for longer discharge times with the activated carbon Li-O₂ cell using the same procedure and setup described in Ref. 3. Two samples were allowed to discharge up to 750 and 1000 mA h/g; these are shown in Fig. 13, along with another image of a 500 mA h/g sample. The maximum particle size is observed to increase up to approximately 500 nm after 16 h, while the number of particles decreases.

The initial condition for this simulation, i.e., the particle sizes and concentrations in the microstructure and electrolyte regimes, is taken from the end of stage 3b. The radii change over time is shown in Fig. 14, indicating that the larger particles continue to grow while the smaller ones initially grow and then start to dissolve. The concentration of LiO₂ at the cathode shown in Fig. 15 reveals a steady decrease over time. This behavior is consistent with the description of the behavior in stage 3 in that the surface LiO₂ concentration continues to drop due to particle growth, which eventually begins to dissolve the smaller particle. An animation of this coarsening process is given in Fig. S2.⁴³ 

Predictions of the particle size at 16 h only reach ∼250 nm in Fig. 15, while the SEM images are approximately 500 nm in Fig. 13. As discussed previously, this discrepancy may be a result of not considering enough small particles dissolving which would increase the growth of the larger particles.

Additional experiments for resting without discharge for the activated carbon Li-O₂ cell were also carried out using the same setup described in Ref. 3. Three samples in Fig. 16 were charged to 500 mA h/g and then allowed to rest for 5, 10, and 24 h before being disassembled, cleaned, and imaged. Interestingly, the particle radii do not appear to change significantly over the time of the experiment.

Simulations continued directly from the end of stage 3b as above but without the applied O₂⁻ flux. The concentration of LiO₂ in the electrolyte decreased due to deposition, without being replenished from the discharge current. The particle radii are plotted over time in Fig. 17 and show that the smaller particle immediately starts to dissolve. An animation of this process is given in Fig. S3.⁴³ Over time, the larger particle begins to consume the smaller; however, this process is slow and results in only a 19 nm increase in radius over a 24 h period. The decrease in the rate of large particle growth is due to the lower concentration of LiO₂ at the cathode as seen in Fig. 7.

The model provides some insight into the interpretation of the experiment. First, the rate of growth of the largest particles decreases substantially since the LiO₂ in the system is not being replenished by the discharge current. Small particles such as the initially 50 nm one quickly dissolve, but then the larger particles would ripen slowly. It should be noted that the current model does not consider disproportionation, reaction 4 in Table I, which would result in shrinkage of volume upon phase change, and counter the particle growth rate to some degree. The reason for the difference between discharge and rest is discussed in Sec. V C.

The interplay between growth and coarsening may be elucidated by comparing discharge beyond 500 mA h/g and resting (see Secs. V A 2 and V B, respectively). These experiments and simulations have the same conditions and differ only in the application of the discharge current. Considering the resting case first, the following analysis is proposed.

The concentration of LiO₂ at the cathode surface corresponds to a particular particle radius. Particles that are larger than this radius will grow, and those that are smaller than it will shrink. During resting, this concentration is determined by the evolution of the particle size distribution alone, which, in the current case, takes time due to the relatively slow attachment kinetics. In Fig. 15, after an initial transient, a group equilibrium is established. While the particles coarsen, the group equilibrium evolves accordingly.

When the discharge current is applied, the LiO₂ concentration is raised above the group equilibrium for the particular particle size distribution. Smaller particles are therefore stabilized and can even grow slightly as seen in Fig. 14. As the particles get larger, their surface area increases and they are able to deplete the electrolyte of LiO₂ significantly compared to the discharge rate. The effect of the current to raise the concentration therefore decreases with particle size and the concentration of LiO₂ drops towards the group equilibrium as seen for the discharge case in Fig. 15. When this happens, the smallest particles are no longer stabilized and will dissolve as seen for particle 2 in Fig. 14.

In some experiments, the discharge product deposits as a coherent film instead of isolated particles.³⁶ Such a transition may be considered in the current work as a consequence of an increased discharge rate.¹⁸ A higher discharge current would result in faster production of LiO₂. In the current work, the rate of LiO₂ production is balanced by the rate of deposition at point B in Fig. 3 before overwhelming it. With higher discharge rates, particles that are growing and nucleating may not have time to coarsen thus resulting in smaller particles for the same capacity.^18,42 At a high enough rate, they may coalesce before point B is reached in which case a film on the cathode would be a likely consequence. Such a film would eventually block the ORR sites in Fig. 1, reducing the incoming flux of $O2−$ and eventually shutting down the reaction.

It has also been found with many cathode materials that the nanoparticles forming on the cathode are composed of Li₂O₂.¹⁴ An estimate of Li₂O₂ solubility calculated in the same way as that of LiO₂ is 6.02 × 10⁻¹⁹M,³⁰ a factor of ∼10⁻¹⁶ lower than that calculated for LiO₂. According to Eq. (10), if particle coarsening were controlled by Li₂O₂ diffusing in the electrolyte, the growth rate would be limited to ∼3 × 10⁻¹⁶ μm/h, and no coarsening would be expected to be observed on the time scale of these experiments. However, this does not rule out the possibility of a different transport mechanism for Li₂O₂ between particles, such as migration along the cathode surface or reaction at the particle-electrolyte interface into Li bearing species with higher solubility. In addition, there is some experimental evidence^44–46 that the solubility of Li₂O₂ may be much higher than the predictions. This may be due to impurities in the electrolyte or surface effects. If so, this would affect the growth rate significantly and it could be similar to LiO₂, assuming Li₂O₂ is formed in solution by some disproportionation mechanism. Alternatively, if Li₂O₂ is formed by LiO₂ disproportionation in the solid phase and it is as soluble like LiO₂ as suggested by experiment, it could also participate in the coarsening.

We have developed a mesoscale model for the production of a condensed lithium superoxide microstructure as the discharge product on the cathode in a Li-O₂ battery. The highly idealized model considers LiO₂ formation through reaction of $O2−$ and Li⁺ in the electrolyte, which deposits on the cathode surface when the LiO₂ concentration reaches the supersaturation limit in the electrolyte. Our model incorporates a rate-equation approach for capturing the production of LiO₂ in the electrolyte by a diffusion-reaction process and a phase-field model that describes the formation of a solid microstructure of LiO₂ through deposition. Together the two approaches capture the diffusion-controlled evolution of the system in two distinct length-scales regimes: whereas reaction-diffusion occurs on the scale of millimeters in the electrolyte between the cathode and anode, microstructure formation and evolution occur on the scale of hundreds of nanometers on the cathode surface. The two regimes share a common boundary: the cathode surface in the electrolyte, and the top of the microstructure regime. At this boundary, the two regimes are coupled strongly by the concentration of LiO₂. A three-stage scheme for particle nucleation, growth, and coarsening is developed, based on a simplified nucleation model.

The mesoscale model requires as input key materials properties, which can be either extracted from experiments or determined by electronic structure calculations. In practice, since few material properties for this system are experimentally available due to LiO₂ being thermodynamically unstable to disproportionation, here most of the input properties come from DFT calculations, such as LiO₂ solubility in DME and elastic properties of the particles. These calculations also provide insight into important atomic-level mechanisms. Moreover, by using our model to fit some experimental data, we were able to extract as of yet unknown material properties, such as the attachment kinetic parameter, surface scaling parameter, and bulk LiO₂ solubility.

The following conclusion can be drawn from the results of this study for the case of nucleation and growth of LiO₂ in a Li-O₂ battery.

1. Our model predicts that coarsening, in which large particles grow and small ones disappear, has a substantial effect on the size distribution of the LiO₂ particles during the discharge process.

2. The size evolution during discharge is the result of the interplay between this coarsening process and particle growth. This growth through continued deposition of LiO₂ has the effect of causing large particles to grow faster and delaying the dissolution of small particles.

3. The predicted size evolution is consistent with experimental results for a previously reported cathode material based on activated carbon during discharge and when it is at rest.

4. The model, even without consideration of more complex microstructural features including LiO₂ disproportionation, can capture the evolving size and number of LiO₂ particles, but not their shape, e.g., toroid formation.

The mesoscale model developed here should have applications in studying other related discharge process such as Li₂O₂ deposition in Li-O₂ batteries and nucleation and growth in Li-S batteries.

Work supported by UChicago Argonne, LLC under Contract No. DE-AC02-06CH11357, the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), the Materials Sciences and Engineering Division, the Computational Materials and Chemical Sciences Network (CMCSN) project on Computational Microstructure Science (phase field modeling); and by the Joint Center for Energy Storage Research, a DOE-BES Energy Innovation Hub (electronic structure computations). We gratefully acknowledge the computing resources provided on Fusion and Blues, high-performance computing clusters operated by the Laboratory Computing Resource Center at Argonne National Laboratory.