THE JOURNAL OF CHEMICAL PHYSICS 143, 224113 (2015)

An atomistically informed mesoscale model for growth and coarsening during discharge in lithium-oxygen batteries Michael J. Welland,1 Kah Chun Lau,1 Paul C. Redfern,1 Linyun Liang,2 Denyun Zhai,3 Dieter Wolf,1 and Larry A. Curtiss1,a)

1

Materials Science Division, Argonne National Laboratory, Argonne, Illinois 60439, USA Mathematics and Computer Science, Argonne National Laboratory, Argonne, Illinois 60439, USA 3 Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne, Illinois 60439, USA 2

(Received 21 July 2015; accepted 29 October 2015; published online 10 December 2015) An atomistically informed mesoscale model is developed for the deposition of a discharge product in a Li-O2 battery. This mescocale model includes particle growth and coarsening as well as a simplified nucleation model. The model involves LiO2 formation through reaction of O2− and Li+ in the electrolyte, which deposits on the cathode surface when the LiO2 concentration reaches supersaturation in the electrolyte. A reaction-diffusion (rate-equation) model is used to describe the processes occurring in the electrolyte and a phase-field model is used to capture microstructural evolution. This model predicts that coarsening, in which large particles grow and small ones disappear, has a substantial effect on the size distribution of the LiO2 particles during the discharge process. The size evolution during discharge is the result of the interplay between this coarsening process and particle growth. The growth through continued deposition of LiO2 has the effect of causing large particles to grow ever faster while delaying the dissolution of small particles. 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, although kinetic factors need to be included. The approach described in this paper synergistically combines models on different length scales with experimental observations and should have applications in studying other related discharge processes, such as Li2O2 deposition, in Li-O2 batteries and nucleation and growth in Li-S batteries. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4936410]

I. INTRODUCTION

Lithium-oxygen batteries have recently received much interest due to their potential for high energy density compared to conventional lithium-ion batteries.1 Li-O2 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-O2 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-O2 battery. The discharge product in Li-O2 batteries is generally composed of lithium peroxide, Li2O2,1,2 although in some cases it may be a mixture of lithium peroxide and lithium superoxide, LiO2.3–6 There have been various different mechanisms proposed for the growth of the discharge product in Li-O2 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 a)[email protected]

0021-9606/2015/143(22)/224113/13/$30.00

cathode surface.7,8 Others have postulated that after oxygen reduction on the surface, subsequent reactions occur in the solution phase involving Li2O2 or LiO2 followed by nucleation and growth of one of these species on the cathode surface.9–12 Nazar et al. reported a Li2O2 formation mechanism by the disproportionation of LiO2 in solution,13 which can account for the lithium peroxide nucleation and growth from solution. Bruce et al. also reported evidence for such a mechanism.14 A through-solution mechanism has also been proposed involving heterogeneous nucleation from a supersaturated solution of LiO2.12,6 The LiO2 mechanism is supported by experimental evidence that LiO2 is soluble in the electrolyte prior to deposition on the surface based on quartz microbalance studies.15 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−O2 battery

143, 224113-1

© 2015 AIP Publishing LLC

224113-2

Welland et al.

J. Chem. Phys. 143, 224113 (2015)

FIG. 1. Schematic of particle formation mechanism from LiO2 production in the electrolyte. The vertical axis is logarithmic. Solvated O2− is produced at active ORR sites on the cathode and diffuses into the electrolyte where it reacts with Li+ to form solvated LiO2. The LiO2 deposits as solid LiO2 particles that are formed on surface nucleation sites. Disproportionation occurs in the solid LiO2 (red) to produce Li2O2 (green). This figure is not meant to represent the exact distribution of LiO2 and Li2O2, which is more complex than shown.

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 LiO2 occurs after its deposition as a bulk material on the surface so that the discharge product is a mixture of LiO2 and Li2O2. The model we develop addresses only the LiO2 deposition and not its disproportionation, which will be reported in a future paper. Modeling the LiO2 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-O2 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 phasefield 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.

II. MODEL DEVELOPMENT

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. Proposed reaction mechanism for formation of discharge product

A schematic of a mechanism for formation of LiO2 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 O2 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 O2− anions. The Li+ and O2− ions diffuse in the electrolyte and react to form LiO2 in the solution. This reaction is thermodynamically favored and should be a fast reaction.12 The process continues until LiO2 reaches saturation and deposits as a LiO2(solid) phase on the surface. The LiO2(solid) phase is thermodynamically unstable with respect to dispropotionation to Li2O2(solid) and release of O2. In a Li-O2 cell, this O2 will be released into the electrolyte, which can then be reduced again according to reaction 1.5 The problem is naturally divided into formation and deposition of solid LiO2, and the transformation from solid LiO2 to Li2O2 and O2. The current work focuses on the first part, i.e., the production of LiO2 particles as described by reactions 1-3. The solid state disproportionation reaction is addressed in a separate paper.19 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.20 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

TABLE I. Reactions leading to the formation of LiO2 particles with subsequent disproportionation to Li2O2(solid) as illustrated in Fig. 1. Reaction number 1 2 3 4

Equation

Process

O2 + e− ↔ O−2 O−2 + Li+ ↔ LiO2 LiO2 ↔ LiO2(solid) 2LiO2(solid) ↔ O2 + Li2O2(solid)

ORR Association Deposition Disproportionation

224113-3

Welland et al.

J. Chem. Phys. 143, 224113 (2015)

of toroidal formation;17 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.

1. Rate-equation model for reaction and diffusion

The species in solution, with concentrations ci (where i = Li+, O2−, or LiO2; see Fig. 1 and Table I), are subject to diffusion and reaction. For each species, the following partial differential equations must be satisfied:

B. Mesoscale model for reaction, diffusion, and deposition

The model described above is a reaction-diffusion type for the production of LiO2 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.21 in which species are assumed to be present everywhere in the domain, but to vanishing small quantities, on the order of 10−23 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.22 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 LiO2(solid), is modeled by an Allen-Cahn type phase-field model.23 The local phase is represented by an independent thermodynamic state variable φ, which varies smoothly in the range [0,1]. Pure solid LiO2 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 function24 p (φ) = φ3[6φ2 − 15φ + 10].

(1)

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 LiO2 phase and the electrolyte, γ, is introduced into the   Gibbs energy functional through the addition of the term γ d3 φ2[1 − φ]2 + 3d|∇φ|2 .21 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 is21 G=

  V



 3 2 ci µi + γ φ [1 − φ]2 + 3d|∇φ|2 , i d

(2)

where ci and µi are the concentration and chemical potential of species i, respectively.

∂ci = −∇·Ji + νi R2, ∂t

(3)

where υi and R2 are the stoichiometric constants and reaction rate for reaction 2 in Table I. Ji is the diffusive flux of species i driven by spatial gradients in the chemical potential, Ji = −

ci Di ∇µi , RgT

(4)

where Di is the diffusion coefficient, Rg 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 = i ci . For this model, we consider the electrolyte to be an ideal solution, which is reasonable for small concentrations of dissolved Li+, O2−, and LiO2. In the LiO2 particle, Li+, O2−, and DME are insoluble and are represented by vanishingly small concentrations. The chemical potential is then given by21 ) p + PVi , + RT ln (x i ) − RT ln 1 − p + qi (

µi =

µ0i

(5)

where µ0i,(amorph) is the reference chemical potential in the solid phase, x i = ci /c is the mole i, Vi is the ( fraction of species ) specific volume, qi = exp µ0i,(amorph)/µ0i,(elec) is the partition coefficient, and µ0i,(elec) is the reference chemical potential of the electrolyte. As discussed by Welland et al.,21 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, csites. In linearized form about P = 0, it is calculated as P = κT [c − csites] ,

(6)

where κT is the bulk modulus. As discussed in Sec. IV, interface energy increases the hydrostatic pressure of small particles, which implies c > csites.   The rate of reaction 2 is R2 = k2 µLiO2 − µLi+ − µO− , 2 where k2 is the associated rate constant and is equal for forward and backward reactions. Mass transport is therefore described by ∂cLi+ = ∂t ∂cO2 = ∂t ∂cLiO2 = ∂t ∂cDME ∂t

− ∇ · JLi+ − R2,

(7a)

− ∇ · JO2 − R2,

(7b)

− ∇ · JLiO2 + R2,

(7c)

= − ∇ · JDME.

(7d)

224113-4

Welland et al.

J. Chem. Phys. 143, 224113 (2015)

2. Phase-field model for deposition

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 δG ∂φ = −k3 ∂t δφ     ∂ µi 6  = −k3 ci +γ φ 1 − 3φ + 2φ2 − 3d∇2φ , i ∂φ d (8) where k 3 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 ci ∂ µ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.

C. Multiscale coupling

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 O2−, Li+, LiO2, 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 LiO2 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 LiO2. In the electrolyte, solvated LiO2 is produced through the reaction of O2− produced at the cathode surface, the diffusion of the anion away from the surface and the reaction with Li+ from the electrolyte; the solvated LiO2 subsequently diffuses away from the surface. Concurrently, microstructure evolution is governed by the local concentration of LiO2 at the cathode, and any growth or shrinkage produces a sink of LiO2 in the electrolyte. The concentration of LiO2 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

FIG. 2. Schematic of coupling between microstructure and electrolyte regimes with a common boundary at the left boundary of the electrolyte and the top of the microstructure. The electrolyte regime shows the fluxes of O−2 and LiO2 at the cathode and an illustrative LiO2 concentration profile. The microstructure regime shows a LiO2 particle growing into the DME from the solid to shaded red circles. The concentration of LiO2 at the boundary in the electrolyte is passed to the microstructure as a fixed value. The corresponding particle growth rate (rate removal of LiO2 from solution) is passed to the electrolyte as an outward flux of LiO2.

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 cLiO2 at the top surface in the microstructure simulations and an s outward flux (sink) of LiO2, JLiO , at the cathode surface for 2 the electrolyte reaction-diffusion model. The iterative scheme s is as follows: given a guess for JLiO , reaction-diffusion in the 2 electrolyte is simulated and a new value of cLiO2 at the cathode is determined. This value is used in the microstructure model s to evolve the particle size and shape. The new guess of JLiO 2 is calculated from the growth rate of those particles by  ∂ cLiO2dVM s JLiO2 = n s s , (9) ∂t where VM is the volume of the microstructure regime. s and cLiO2 are Iterations continue until consistent values of JLiO 2 found. Instead of using Eq. (9) directly to calculate the next s guess of JLiO , a Newton method is employed which increases 2 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. D. Scheme for nucleation, growth, and coarsening

The following three-stage scheme is proposed to model the nucleation, growth, and coarsening process for an initially

224113-5

Welland et al.

J. Chem. Phys. 143, 224113 (2015) TABLE II. Material properties for species i in the LiO2(solid) and DME electrolyte phases. Diffusion coefficients are from Ref. 20. The solubilities of Li+, O2−, and DME in the LiO2(solid) phase are vanishingly small as discussed in the text, and the LiO2 species is the molar density predicted by DFT calculations. The solubilities in the DME electrolyte phase are discussed in DME is fit to the particular experiment. the text, and C LiO 2

Di (cm2/s) CDME (M) i LiO Ci 2(M)

FIG. 3. Concentration of LiO2 in the electrolyte at the cathode vs time during the proposed particle evolution scheme. The various stages and points are described in the text.

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: LiO2 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 LiO2, 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 LiO2 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).

III. MATERIAL PROPERTIES

Few material properties for this system are experimentally available due to LiO2 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+, O2−, LiO2, and DME are summarized in Table II. For Li+ and O2,

Li+

O 2−

LiO2

DME

1.0 · 10−5 0.1 10−23

4.0 · 10−5 10−10 10−23

1 · 10−5 10−3 56.2

1 · 10−5 56.2 10−23

material properties are obtained from Lu et al.25 for the case of 0.1 M Li+ in DME. The diffusion coefficient of O2− is assumed equal to that of O2. As a first approximation, the diffusion coefficients DLiO2 and DDME 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+, O2−, and DME are only present in vanishingly small concentrations. In our simulations, the densities of both phases are held at the value of LiO2. This simplification avoids the need to consider flow of the electrolyte upon growth or shrinkage of the LiO2 particle and the corresponding advective currents for species transport, greatly reducing the computational complexity. The rate of production of LiO2 in the electrolyte from the applied discharge current is governed by the rate constant k 2 and by the equilibrium concentrations of Li+, O2−, and LiO2 in the electrolyte. It is thought that this reaction proceeds very quickly and quickly consumes all the O2−.12 Therefore, we assign k 2 = 104 mol s/m−3 and CODME = 10−10, which is − 2 sufficient during the simulations to immediately react O2− 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-O2 cell based on an activated-carbon cathode;3 this value may differ for other systems. A. Density and bulk modulus

To compute the density and bulk modulus of bulk amorphous LiO2, 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 LiO2 structure. The static E(V) dataset includes the self-consistent total energies calculated using the VASP DFT code26 with periodic boundary conditions for different volumes with TABLE III. General material properties. Density LiO2(solid)a Molar weight LiO2 Bulk modulus LiO2(solid)a Density DME a From

DFT calculations described in text.

2.19 g/cm3 38.94 g/mol 85.82 GPa 0.8683 g/cm3

224113-6

Welland et al.

FIG. 4. The total energy, E(V), vs. volume; the inset shows the simulation cell used for the DFT calculations on bulk amorphous LiO2 (containing 54 LiO2 units) in DFT calculation. The minimum of the E(V) corresponds to the free energy formation 2.37 eV/Li, slightly less stable than the LiO2 crystalline structure (i.e., 2.61 eV/Li) as reported in our previous study.29

J. Chem. Phys. 143, 224113 (2015)

solvent using various explicit and implicit solvent models, as well as ab initio molecular dynamics. Best estimates for the LiO2 solvation energy from these calculations along with the calculated lattice energy of LiO2 were used to determine the solubility of bulk LiO2. The best estimate for the solubility of LiO2 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 LiO2 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. C. Surface energy between LiO2 and DME

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−5 eV for the K-point integration, and all the atoms and geometries were optimized until the residual forces became less than 1 × 10−2 eV/Å within a single K-point grid (i.e., Γ-point). The initial configuration of the DFT-predicted LiO2 bulk amorphous structure is obtained by 54 randomly packed LiO2 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 LiO2 structure is found to be ∼2.19 g/cm3, slightly smaller than the predicted atomic density of crystalline Li2O2 crystal of ∼2.37 g/cm3 reported in our previous study.27 By using the quasi-harmonic Debye model,28 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.

To estimate the thermodynamic stability of the interface between amorphous LiO2 and the electrolyte, we used data from previously published ab initio molecular dynamics (AIMD) simulations.6 Based on the DME electrolyte (i.e., ∼0.5M LiCF3SO3/DME) that has an interface with the LiO2 amorphous slab having a thickness of ∼1.2 nm (i.e., 525 atoms in a simulation cell), we optimized the amorphous LiO2/electrolyte interface and found it to be stable with an interfacial energy of approximately −64.5 meV/Å2.6 From the AIMD simulation, the solvated LiO2 amorphous surface is found to be stabilized by DME and LiCF3SO3 molecules at the interfaces. The calculations showed significant charge transfer and strong adsorption of DME onto the LiO2 surface. Surface adsorption can produce negative surface energy without destabilizing it;31 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 multiphysics 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 Li2O2, [0001] of 1.136 J/m2.7 The effects of this uncertainty will be discussed in Sec. V.

B. Solubility of LiO2

The solubility of LiO2 in the electrolyte in equilibrium ∞ with a flat interface (i.e., radius = ∞), cLiO , is one of the 2 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 LiO2 in any electrolyte, although there is experimental evidence for its solubility from quartz microbalance studies.15 Thus, the only available solubility information for LiO2 comes from computational studies. In a separate study,30 we have calculated solvation free energies of molecular LiO2 in a DME

D. Attachment kinetics

The interface attachment kinetic factor k3 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-Wagner32,33 theory predicts that the rate of change of the average radius ⟨R⟩ is given by

224113-7

Welland et al.

J. Chem. Phys. 143, 224113 (2015)

∞ 2 ∂⟨R⟩ 8γcLiO2VLiO2 DLiO2 1 . = ∂t 9Rg T 3⟨R⟩2

(10)

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-O2 cell,3 this is clearly not the case even ∞ 3 considering large uncertainty in cLiO . Therefore, we consider 2 that k 3 is relatively slow, and coarsening is attachment-kinetics controlled. The low value of k 3 may be related to both the solvation shell around LiO2 in solution and the binding of DME molecules on the LiO2 surface. Deposition of LiO2 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.

IV. IMPLEMENTATION AND VALIDATION

The model is implemented using the FEniCS package,34–38 an open source Finite Element Method package, and configured to use the PETSc39,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.41 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, µiR

−

µ∞ i

2Vi σ , = R

(11)

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 ( ) x iR 2Vi σ = exp . x∞ Rg T R i

(12)

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.

FIG. 5. The change in LiO2 solubility in the electrolyte of a particle with radius R, compared to a flat interface represented by R = ∞. Above—the change in chemical potential compared to the classical solution in Eq. (11). Below—the corresponding ratio of mole fractions of LiO2 considering an ideal solution.

FIG. 6. Simulated concentration and pressure profiles vs radius for a spherically symmetric particle of radius 200 nm. Large variations between the particle phase (left) and electrolyte (right) show large variations in species concentrations and the gauge pressure.

V. APPLICATION TO EXPERIMENTS

We now consider application of the phase-field model described in Secs. III and IV to results reported for a Li−O2 battery based on activated carbon3 and some additional experiments that have been carried out for this paper to help distinguish growth from coarsening. In the Li-O2 battery considered in the experiments powdered activated carbon was loaded onto a 1 cm diameter disk to a density of approximately 1.6 mg/cm2. The density of activated carbon is 2 mg/cm3, which implies a layer approximately 8 µm thick depending on particle packing. The activated carbon has surface area of approximately 3840 m2/g which is mostly attributed to nanopores, but

224113-8

Welland et al.

includes microscopic surface features and the surface area of a range of particle sizes.3 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 images3 shows the nucleation site density to be approximately 12 sites/µm2. We assume a critical nucleus radius for our simplified nucleation model of 10 nm, which seems physically reasonable. A. Particle evolution during discharge

In the reported results for the activated carbon Li-O2 battery,3 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. 1. Discharge up to first 500 mA h/g capacity

Initially, the cathode surface is free of LiO2 and the electrolyte is saturated in Li+ but without dissolved O−2 or LiO2. As discharge proceeds, the concentration of LiO2 in the electrolyte at the cathode follows the simulation results presented in Fig. 7. a. Stage 1. During stage 1, the discharge current introduces O−2 into the system at the cathode surface,

J. Chem. Phys. 143, 224113 (2015)

which immediately reacts with Li+ and forms LiO2. The concentration of LiO2 at the cathode surface rises until 10 nm the nucleation supersaturation level, cLiO , is reached. We 2 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 reactiondiffusion 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 LiO2 is ∞ 13.4 mM, which corresponds to cLiO = 2.6 mM from Eq. 2 (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 ∞ = 2.62% ± 16% mM. cLiO 2 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. b. Stage 2. During stage 2, LiO2 particles nucleate on the cathode surface and become sinks for solvated LiO2. Since the growth rate is proportional to the surface area, the largest LiO2 particles, i.e., the first to have nucleated, will contribute most to the LiO2 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

FIG. 7. Simulated concentrations of Li+ and LiO2 in the electrolyte at the cathode vs time following the scheme illustrated in Fig. 3. Horizontal lines indicate R , in equilibrium with a particle of radius R. From t = 0, the concentration rises past the bulk solubility limit c ∞ until it reaches point A where solubility, c LiO LiO2 2 the nucleation size limit is reached (assumed to be 10 nm particles for this case based on a cusp in the discharge voltage profile). At point B, the LiO2 deposition rate overcomes the LiO2 production rate and a maximum is predicted. At point C, the super-saturation has fallen below the threshold required to support the critical nucleation size. Particles continue to grow and deplete the LiO2 in the electrolyte past point D, which is where particles of ∼50 nm notably begin to disappear on SEM images for the activated carbon. Insets are reprinted with permission from Zhai et al., J. Am. Chem. Soc. 135, 15364 (2013). Copyright 2013 American Chemical Society.

224113-9

Welland et al.

J. Chem. Phys. 143, 224113 (2015)

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.

FIG. 8. Simulated concentration profiles of O−2 , Li+, and LiO2 in the electrolyte relative to the cathode at: 0.46 h, point A, the point at which the critical nucleation size is reached; 0.79 h, point B, the maximum in LiO2 concentration; 4.8 h, point D, the reference point; and 8.0 h, the end of stage 3, as indicated in Fig. 7.

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. c. Stage 3. In stage 3, as the concentration of LiO2 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 results3 at 4.8 h (300 mA h/g). This image was chosen since there

FIG. 9. Particle radii versus time during stages 2 and 3. Stage 3b considers the co-evolution of four particles in order to capture coarsening.

d. Stage 3a. 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 50 nm implies cLiO2 = cLiO at the cathode surface. The largest 2 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 50 nm concentration cLiO , we may fit the attachment kinetic term, 2 k3, and surface scaling parameter, s, to 2 × 10−11 m3/(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.3 The current work assumes this ratio to remain the same for all samples, but other experiments may employ different particle size distribution or structures.42 e. Stage 3b. 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.43 The

FIG. 10. Idealized schematic of particle size distribution on surface.

224113-10

Welland et al.

J. Chem. Phys. 143, 224113 (2015)

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.3 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. 2. Longer discharge

FIG. 11. Simulated LiO2 concentration at (a) the beginning (4.8 h) and (b) the end (8 h) of stage 3b. Particle radii are initially 20, 50, 100, and 150 nm. The electrolyte, determined by φ > 0.5, is made translucent. The concentration of LiO2 inside the particles exceeds the normal density due to compression by surface energy.

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 LiO2 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 carbon3 shows particles of almost 250 nm, whereas the simulation

FIG. 12. Simulated concentration of LiO2 along the x-axis of Fig. 11 between the initially 150 and 50 nm particles. The region between particles is enlarged in the inset showing the gradient between particles, and the lowering of the solubility with increasing particle size.

As part of this study, we carried out additional experiments for longer discharge times with the activated carbon Li-O2 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 LiO2 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 LiO2 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.43 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. B. Evolution during resting

Additional experiments for resting without discharge for the activated carbon Li-O2 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 O2− flux. The concentration of LiO2 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.43 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

224113-11

Welland et al.

(a)

J. Chem. Phys. 143, 224113 (2015)

(b)

(c)

FIG. 13. SEM images of samples at longer discharge capacities. As capacity increases, there is an increase in maximum particle radius and decrease in particle abundance. (a) 8 h (500 mA h/g), (b) 12 h (750mAh/g) and (c) 16 h (1000 mA h/g) .

to the lower concentration of LiO2 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 LiO2 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. C. Interplay between growth and coarsening

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 LiO2 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

FIG. 14. Simulated particle radii and their average while continuing to discharge from 500 mA h/g to 1000 mA h/g.

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 LiO2 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 LiO2 significantly compared to the discharge rate. The effect of the current to raise the concentration therefore decreases with particle size and the concentration of LiO2 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.

VI. IMPLICATIONS OF MODEL TO OTHER SCENARIOS

In some experiments, the discharge product deposits as a coherent film instead of isolated particles.36 Such a transition may be considered in the current work as a consequence of an increased discharge rate.18 A higher discharge current would result in faster production of LiO2. In the current work, the rate of LiO2 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

FIG. 15. Concentration of LiO2 at the cathode following stage 3 (500 mA h/g), for the discharging and resting experiments.

224113-12

Welland et al.

(a)

J. Chem. Phys. 143, 224113 (2015)

(b)

(c)

FIG. 16. SEM images of particles that are allowed to rest in the electrolyte following discharge to 500 mA h/g. The maximum particle size and abundance appear not to change qualitatively. (a) 5 h, (b) 10 h and (c) 24 h.

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 O−2 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 Li2O2.14 An estimate of Li2O2 solubility calculated in the same way as that of LiO2 is 6.02 × 10−19M,30 a factor of ∼10−16 lower than that calculated for LiO2. According to Eq. (10), if particle coarsening were controlled by Li2O2 diffusing in the electrolyte, the growth rate would be limited to ∼3 × 10−16 µ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 Li2O2 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 evidence44–46 that the solubility of Li2O2 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 LiO2, assuming Li2O2 is formed in solution by some disproportionation mechanism. Alternatively, if Li2O2 is formed by LiO2 disproportionation in the solid phase and it is as soluble like LiO2 as suggested by experiment, it could also participate in the coarsening.

FIG. 17. Simulated particle radii while resting with no further discharge, following discharge to 500 mA h/g.

VII. SUMMARY AND CONCLUSIONS

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-O2 battery. The highly idealized model considers LiO2 formation through reaction of O−2 and Li+ in the electrolyte, which deposits on the cathode surface when the LiO2 concentration reaches the supersaturation limit in the electrolyte. Our model incorporates a rate-equation approach for capturing the production of LiO2 in the electrolyte by a diffusion-reaction process and a phase-field model that describes the formation of a solid microstructure of LiO2 through deposition. Together the two approaches capture the diffusion-controlled evolution of the system in two distinct length-scales regimes: whereas reactiondiffusion 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 LiO2. 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 LiO2 being thermodynamically unstable to disproportionation, here most of the input properties come from DFT calculations, such as LiO2 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 LiO2 solubility. The following conclusion can be drawn from the results of this study for the case of nucleation and growth of LiO2 in a Li-O2 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 LiO2 particles during the discharge process.

224113-13

Welland et al.

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 LiO2 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 LiO2 disproportionation, can capture the evolving size and number of LiO2 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 Li2O2 deposition in Li-O2 batteries and nucleation and growth in Li-S batteries. ACKNOWLEDGMENTS

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. 1P.

G. Bruce, S. A. Freunberger, L. J. Hardwick, and J.-M. Tarascon, Nat. Mater. 11, 19 (2012). 2J. Lu, Y. Lei, K. C. Lau, X. Luo, P. Du, J. Wen, R. S. Assary, U. Das, D. J. Miller, J. W. Elam, H. M. Albishri, D. A. El-Hady, Y.-K. Sun, L. a. Curtiss, and K. Amine, Nat. Commun. 4, 2383 (2013). 3D. Zhai, H. H. Wang, J. Yang, K. C. Lau, K. Li, K. Amine, and L. A. Curtiss, J. Am. Chem. Soc. 135, 15364 (2013). 4J. Yang, D. Zhai, H.-H. Wang, K. C. Lau, J. A. Schlueter, P. Du, D. J. Myers, Y.-K. Sun, L. A. Curtiss, and K. Amine, Phys. Chem. Chem. Phys. 15, 3764 (2013). 5D. Zhai, H. H. Wang, K. C. Lau, J. Gao, P. C. Redfern, F. Kang, B. Li, E. Indacochea, U. Das, H. H. Sun, H. J. Sun, K. Amine, and L. A. Curtiss, J. Phys. Chem. Lett. 5, 2705 (2014). 6D. Zhai, K. C. Lau, H.-H. Wang, J. Wen, D. J. Miller, J. Lu, F. Kang, B. Li, W. Yang, J. Gao, E. Indacochea, L. A. Curtiss, and K. Amine, Nano Lett. 15, 1041 (2015). 7J. S. Hummelshoøj, a. C. Luntz, and J. K. Nørskov, J. Chem. Phys. 138, 034703 (2013). 8V. Viswanathan, K. S. Thygesen, J. S. Hummelshøj, J. K. Nørskov, G. Girishkumar, B. D. McCloskey, and A. C. Luntz, J. Chem. Phys. 135, 214704 (2011). 9S. H. Oh, R. Black, E. Pomerantseva, J.-H. Lee, and L. F. Nazar, Nat. Chem. 4, 1004 (2012). 10X. Ren, S. S. Zhang, D. T. Tran, and J. Read, J. Mater. Chem. 21, 10118 (2011).

J. Chem. Phys. 143, 224113 (2015) 11R.

R. Mitchell, B. M. Gallant, Y. Shao-Horn, and C. V. Thompson, J. Phys. Chem. Lett. 4, 1060 (2013). 12J. Lu, L. Cheng, K. C. Lau, E. Tyo, X. Luo, J. Wen, D. Miller, R. S. Assary, H.-H. Wang, P. Redfern, H. Wu, J.-B. Park, Y.-K. Sun, S. Vajda, K. Amine, and L. A. Curtiss, Nat. Commun. 5, 4895 (2014). 13R. Black, S. H. Oh, J.-H. Lee, T. Yim, B. Adams, and L. F. Nazar, J. Am. Chem. Soc. 134, 2902 (2012). 14L. Johnson, C. Li, Z. Liu, Y. Chen, S. a. Freunberger, P. C. Ashok, B. B. Praveen, K. Dholakia, J.-M. Tarascon, and P. G. Bruce, Nat. Chem. 6, 1091 (2014). 15S. Schaltin, G. Vanhoutte, M. Wu, F. Bardé, and J. Fransaer, Phys. Chem. Chem. Phys. 17, 12575 (2015). 16U. Das, K. C. Lau, P. C. Redfern, and L. A. Curtiss, J. Phys. Chem. Lett. 5, 813 (2014). 17N. B. Aetukuri, B. D. McCloskey, J. M. García, L. E. Krupp, V. Viswanathan, and A. C. Luntz, Nat. Chem. 7, 50 (2014). 18B. Horstmann, B. Gallant, R. Mitchell, W. G. Bessler, Y. Shao-Horn, and M. Z. Bazant, J. Phys. Chem. Lett. 4, 4217 (2013). 19L. Liang, “A phase field model for a LiO -Li O transformation in a lithium2 2 2 air battery” (unpublished). 20S. Lau and L. A. Archer, Nano Lett. 15, 5995 (2015). 21M. J. Welland, D. Wolf, and J. E. Guyer, Phys. Rev. E 89, 012409 (2014). 22S. G. Kim, W. T. Kim, and T. Suzuki, Phys. Rev. E 60, 7186 (1999). 23S. Allen and J. Cahn, Acta Metall. 20, 423 (1972). 24W. J. Boettinger, J. A. Warren, C. Beckermann, and A. Karma, Annu. Rev. Mater. Res. 32, 163 (2002). 25Y.-C. Lu, D. G. Kwabi, K. P. C. Yao, J. R. Harding, J. Zhou, L. Zuin, and Y. Shao-Horn, Energy Environ. Sci. 4, 2999 (2011). 26G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996). 27K. Lau, D. Qiu, X. Luo, J. Greeley, L. A. Curtiss, J. Lu, and K. Amine, Energies 8, 529 (2015). 28A. Otero-de-la-Roza, D. Abbasi-Pérez, and V. Luaña, Comput. Phys. Commun. 182, 2232 (2011). 29K. C. Lau, L. A. Curtiss, and J. Greeley, J. Phys. Chem. C 115, 23625 (2011). 30L. Cheng, “Computational studies of solubility of LiO and Li O in aprotic 2 2 2 solvents” (unpublished). 31A. Mathur, P. Sharma, and R. C. Cammarata, Nat. Mater. 4, 186 (2005). 32I. Lifshitz and V. Slyozov, J. Phys. Chem. Solids 19, 35 (1961). 33C. Wagner, Z. Elektrochem 65, 581 (1961). 34A. Logg, K.-A. Mardal, and G. N. Wells, Automated Solution of Differential Equations by the Finite Element Method (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012). 35A. Logg, G. N. Wells, and J. Hake, in Automated Solutions of Differential Equations by the Finite Element Method, edited by A. Logg, K.-A. Mardal, and G. N. Wells, Lecture Notes in Computational Science and Engineering (Springer, 2012), Vol. 84. 36M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells, ACM Trans. Math. Softw. 40, 1 (2014). 37R. C. Kirby, ACM Trans. Math. Softw. 30, 502 (2004). 38R. C. Kirby and A. Logg, ACM Trans. Math. Softw. 32, 417 (2006). 39S. Balay, J. Brown, K. Buschelman, D. W. Gropp, D. Kaushik, M. Knepley, L. McInnes, and F. B. Smith, http://www.mcs.anl.gov/petsc. 40S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, in Model Software Tools Science Computing, edited by E. Arge, A. M. Bruaset, and H. P. Langtangen (Birkhäuser Boston, Boston, MA, 1997), pp. 163–202. 41C. H. P. Lupis, Chemical Thermodynamics of Materials (North-Holland, 1983). 42B. D. Adams, C. Radtke, R. Black, M. L. Trudeau, K. Zaghib, and L. F. Nazar, Energy Environ. Sci. 6, 1772 (2013). 43See supplementary material at http://dx.doi.org/10.1063/1.4936410 for animations of simulation results. 44T. K. Zakharchenko, A. Y. Kozmenkova, D. M. Itkis, and E. A. Goodilin, Beilstein J. Nanotechnol. 4, 758 (2013). 45S. J. Visco, V. Y. Nimon, A. Petrov, K. Pridatko, N. Goncharenko, E. Nimon, L. De Jonghe, Y. M. Volfkovich, and D. A. Bograchev, J. Solid State Electrochem. 18, 1443 (2014). 46D. Zheng, H. S. Lee, X. Q. Yang, and D. Qu, Electrochem. Commun. 28, 17 (2013).