DFT modelling of explicit solid–solid interfaces in batteries: methods and challenges

Kevin Leung
Sandia National Laboratories, MS 1415, Albuquerque, NM 87185, USA. E-mail: kleung@sandia.gov

Received 30th November 2019 , Accepted 4th February 2020

First published on 11th February 2020

Density Functional Theory (DFT) calculations of electrode material properties in high energy density storage devices like lithium batteries have been standard practice for decades. In contrast, DFT modelling of explicit interfaces in batteries arguably lacks universally adopted methodology and needs further conceptual development. In this paper, we focus on solid–solid interfaces, which are ubiquitous not just in all-solid state batteries; liquid-electrolyte-based batteries often rely on thin, solid passivating films on electrode surfaces to function. We use metal anode calculations to illustrate that explicit interface models are critical for elucidating contact potentials, electric fields at interfaces, and kinetic stability with respect to parasitic reactions. The examples emphasize three key challenges: (1) the “dirty” nature of most battery electrode surfaces; (2) voltage calibration and control; and (3) the fact that interfacial structures are governed by kinetics, not thermodynamics. To meet these challenges, developing new computational techniques and importing insights from other electrochemical disciplines will be beneficial.

I Introduction

High energy density batteries are instrumental to vehicle electrification and to grid storage which helps alleviate intermittency in solar and wind energy generation. Deeper understanding of existing materials and the search of new electrode and electrolyte materials are crucial for developing higher capacity, faster-charging, safer, and longer-lasting batteries.1 Battery interfaces are also universally acknowledged to be critical for governing battery rate capability and stability.2 Charge/discharge processes, as well as many side reactions that lead to degradation, self-discharge, and thermal runaway, initiate at interfaces.

Electronic structure Density Functional Theory (DFT) modelling of the crystalline interior of electrodes, for the purpose of predicting phase stability, equilibrium voltages, and lithium diffusion kinetics, has been widely practised for decades.3–7 It provides important guidance to experiments. In contrast, DFT modelling of explicit battery interfaces arguably needs further conceptual development and systematization. Here we distinguish explicit interfaces where two phases are in contact in the same simulation cell, from single-phase DFT calculations used to infer interfacial properties indirectly.8 The present work focuses on modelling methods, offers somewhat pedagogical discussions on the rationale behind the DFT approaches used in our group, and examines future directions and improvements. Simple material interfaces, mainly involving lithium metal, are used as illustrations, but the principles involved should be broadly applicable to other electrodes. Even simple materials exhibit complex interfaces that require substantial approximations; the nature and consequences of some key implicit approximations are highlighted. We restrict ourselves to vanishing current densities.9 This paper is intended to be a topical, critical overview, not a comprehensive review of battery interface DFT modelling.10–12

A useful illustration of the multiplicity of interfacial processes is the first charge of graphite anodes used in commercial liquid electrolyte-based lithium ion batteries (Fig. 1a). The “open circuit” potential (OCP) of pristine graphite in organic carbonates, e.g., ethylene carbonate, EC, mixed with salts and cosolvents, is ∼3 V vs. Li+/Li(s) (henceforth this is the reference used unless noted). We stress that OCP is measured by immersing the electrode in the electrolyte over experimental time scale, implying that the interface is well-equilibrated. As charging begins and the voltage drops, graphite acts as an inert electrode like those in electrochemical capacitors. Excess e go to the surface and their negative charges are compensated by an enhanced local concentration of cations in the electrolyte near the interface; the familiar electric double layer (EDL, effectively a dipole sheet) is formed. As the voltage approach ∼0.7 V, graphite acts as an electron emitter – possibly also as an electrocatalyst at its edge site functional groups – and the organic solvent is electrochemically reduced in what will be called “parasitic reactions.” This is because high energy batteries typically operate at voltages outside the electrochemical stability limit of electrolytes.13–15 (Counterions can be reduced at higher voltages, apparently more slowly and to a lesser extent.16–18) The passivating “solid electrolyte interphase” (SEI)10,13,19 grows from these decomposed electrolyte fragments, but e can tunnel or diffuse through the nascent film, possibly as part of Li atoms, until later charging cycles. The EDL, which sustains the applied voltage, may now partly reside in the solid film, not just in the liquid (Fig. 1b). Li+ does not intercalate into graphite sheets until 0.2–0.1 V (not counting edge sites). This small 0.1 V voltage window has seen the majority of DFT studies,20 in the bulk LixC6 region. By continuity, in this regime, the voltage must be a manifestation of both the graphite electronic structure which changes with x, and the EDL contribution. The SEI continues to grow/evolve in this window. When a LiC6 stoichiometry is attained, the fully charged graphite reverts to an inert electrode, albeit coated with SEI, that may plate lithium metal at negative voltages. Lithium plating is detrimental to battery life and safety. Related processes occur on transition metal oxide cathode surfaces, on which “cathode electrolyte interphase” (CEI) is formed (Fig. 1b); however, the voltage-dependent parasitic reactions there most likely involve species in contact with the surface, just like in water splitting, rather than long-range electron transfer like on the anode.

image file: c9cp06485k-f1.tif
Fig. 1 Schematics illustrating (a) first charge of a graphite anode in lithium ion batteries; (b) proposed voltage profile on battery with “dirty” electrode surfaces. “iSEI” and “oSEI” refer to inorganic and organic SEI components. The iSEI has nanometer thickness. The separator occupies the electrolyte region and is omitted.

Thus, as our previous overview emphasized,21 batteries interfaces embody the rich physics of multiple electrochemical devices, and DFT-based battery modelling benefits from borrowing from diverse computational disciplines. Electrified liquid/solid interfaces on passive, pristine carbon or platinum surfaces are relevant to electrochemical capacitors and electrocatalysis; these areas have well-established DFT modelling protocols.22–32 Transition metal oxide/liquid interfaces in batteries are in some ways related to electrochemical and photoelectrochemical water splitting.33–36 Electrodes covered with passivating solid films are relevant to breakdown of surface oxide films which protect metal surfaces against corrosion, and constructions similar to Fig. 1b have been invoked in corrosion studies.37,38,40

Passivating films on organic-solvent-based battery anode surfaces have inorganic (iSEI) and organic (oSEI) components (Fig. 1b),62 leading to multiple solid–solid interfaces. We argue that it is as urgent to model these comparatively neglected inorganic solid interfaces as the solid–liquid interfaces which have been the mainstay of computational electrochemistry. DFT simulations have been conducted on explicit battery solid–solid interfaces relatively recently,15,41–49 partly because all-solid-state batteries with ceramic- or sulfide-based solid electrolytes constitute a timely research area. Many aspects of solid–solid interface studies in liquid-electrolyte batteries are informed by and are transferrable to these solid electrolytes.15,44 The main difference is the solid film/electrolyte thickness - nanometers vs. microns. The much thicker solid electrolytes come with a wider “space charge” region50,69 which can potentially be modelled using continuum methods coupled to DFT currently applied to the Mott–Schottky layer at semiconductor interfaces.28 One theme common to both is that the charge carrying cation M+ (Li+ and Na+) is highly mobile by solid state conductivity standards (otherwise the material would not be used in batteries). Hence the total number of M atoms should vary with voltage at the interface as well as inside electrodes.

We identify three computational challenges somewhat unique to battery interfaces. One is the ubiquity of the aforementioned solid–solid interfaces. Solid surface films in which M+ can diffuse limit the utility of modelling techniques traditionally used in the liquid state like DFT-based molecular dynamics (also called ab initio molecular dynamics or AIMD).10,51–54 This is because of time scale mismatch: ionic motion and relaxation in solid occur many orders of magnitude slower than liquid state diffusion.

Another critical issue is voltage determination and control. First we stress that an assumption often invoked when modelling explicit battery interfaces, namely that voltages in interfacial simulation cells are determined by the lithium content/energetics, is incorrect. Instead, we adopt the definition used in other electrochemical disciplines22–30,55 and reconcile it with traditional battery modelling approaches. Using the correct definitions is necessary to establish equilibrium, and to construct simulation cells at overpotentials conditions.

Finally, battery interfaces are determined by kinetics, not thermodynamics. Inorganic electrode and electrolyte materials are often sintered for hours at ≥900 °C. A phase diagram, thermodynamics modelling approach suffices under those conditions.8 In contrast, interfaces that involve organic solvents are assembled at room temperatures where chemical reactions are slow. Even all-solid-state battery interfaces are usually created at modest (<200 °C) temperatures, where thermodynamic equilibrium does not always apply,56 although thermodynamic phase diagrams for bulk8 and surfaces57 provide the foundation for kinetics investigations. Thermodynamic (single-phase) calculations are elegant and systematic.8,12 Interfacial kinetics calculations tend to be idiosyncratic. They depend on the realism of the model interface, lattice matching, and other issues common to all interface modelling efforts. Here we focus on issues unique to electrochemical interfaces.

It must be stressed that a reaction with a large exothermicity (ΔG) does not necessarily exhibit a low barrier (ΔG*). One famous example is the Grotthuss mechanism for proton exchange in water:

H2O + OH → OH + H2O.
The left and right sides are the same, and ΔG is exactly zero. Thus the thermodynamic driving force is the lowest possible; a more unfavorable ΔG would stop the reaction from occurring. Yet the reaction is almost barrierless (ΔG* = 0).58 Hence kinetics and thermodynamics are not always correlated.

Kinetics is particularly relevant to the viability of conversion cathode materials59 and to stability against parasitic reactions.10 Currently DFT modelling of conversion cathodes is focused on elucidating metastable intermediate states,59 not on computing voltage-dependent reaction barriers;60 in electrocatalysis, both intermediates and barriers are known contribute to overpotentials.60 Finally, kinetics approaches may not identify the final products. To circumvent this, whenever possible, we have focused on elucidating the stability of proposed surface film components rather than start parasitic reaction calculations from pristine electrolyte and pristine electrode surfaces.15

II Methods and models

DFT calculations only deal with half cells because DFT is a ground state theory which supports a single Fermi level (EF). Hence we insert a vacuum region into Fig. 1b, generating two half cells (Fig. 2a). This construction is rigorous for liquid electrolyte systems61 and has been applied in DFT modelling work.31,32
image file: c9cp06485k-f2.tif
Fig. 2 (a) Inserting a vacuum region into the liquid electrolyte of Fig. 1b. This step is rigorous. (b) Reduced anode model and cathode models, introducing approximations. A metallic current collector is added at the backside of the cathode, (c) schematic of solid electrolytes with vacuum inserted. (d) Schematic of metal oxide-coated metal surface, relevant to corrosion. Arrows with electrons indicate formal paths for calculating image file: c9cp06485k-t1.tif; in physical systems it may be ions that cross the boundaries. M = lithium in this figure.

A. Models of thin-film-covered anode interfaces

Focusing on the anode, the multilayer SEI model62 suggests that inorganic products coat the active anode material, separating it from the thicker, amorphous, and porous organic/polymeric layer outside.73 The organic SEI is complex, and likely exhibits a low dielectric constant. Its micro- and even atomic-structures are subjects of current studies and debate.63–65 SEI chemical compositions also evolve with time.66,67

To make progress, we replace the outer organic SEI (oSEI), with vacuum (Fig. 2b). This model focuses on the anode active material/innermost inorganic SEI (iSEI) interface, expected to play a main role in blockage of e tunneling to the electrolyte. Our liquid-free, charge-neutral simulation cells resemble those we use to model all-solid-state battery interfaces (Fig. 2c).56 In principle, the omitted electrostatic contribution from a liquid electrolyte can be computed using classical force-field-based molecular dynamics, but this is not applicable when an organic SEI covers the inorganic film. (To give an order-of-magnitude estimate, the liquid EC/vacuum interface exhibits a ∼0.3 V surface potential,55 dwarfed in magnitude by solid–solid interface effects discussed below.) Another approach might be to replace the vacuum in Fig. 2b with implicit solvation.22,27,68 While noble metal–liquid interactions have been calibrated in implicit solvation models, oxide/water or oxide/organic liquid interactions have not. Currently no dielectric solvation model is universally available or uniformally implemented into every DFT code. To enhance reproducibility, we work with vacuum interfaces. The EDL resides entirely within the thin solid film regions in our models. Our approach is clearly a limiting-case approximation, to be improved upon in the future.

Our half-cell model for solid–electrolyte/electrode interfaces (Fig. 2c) has an unavoidable ambiguity: the absolute work function Φ (and hence the estimated voltage, see below) depends on the nature of the cleaved facet, unlike liquids which relax to universal liquid–vapor interfaces. Kelvin Probe Force Microscopy (KPFM) measurements of voltages in solids share this ambiguity.69 The voltage difference between cathode and anode does not depend on the cleaved surface provided the two vacuum–solid electrolyte interfaces are identical. Vacuum regions in fact exist in solid electrolytes, inside pores and cracks.70

B. Voltage definitions for metallic electrodes

Unique among electrochemical devices, battery electrodes can emit both electrons and the charge-carrying cation M+. The primary battery functions are charge and discharge, which involve M+ intercalation and exit from electrodes into the electrolyte, accompanied by e flow in the external circuit. However, as discussed above, e can also be emitted into (reduction) or captured from (oxidation) the electrolyte, accompanied by local M+ redistribution but not M+ transfer between electrodes.

To account for these separate functions, we define two voltages, image file: c9cp06485k-t2.tif and image file: c9cp06485k-t3.tif.55 The “ionic” voltage image file: c9cp06485k-t4.tif is the canonical definition in the battery DFT literature.3 It measures the energy gain by inserting an M atom, referenced to M(s) → M(atom):

image file: c9cp06485k-t5.tif(1)
where EnM is the total energy of simulation cell with an electrode with nM M atoms, μM is the M chemical potential in its bulk metal phase, and |e| is the electronic charge. The entropy change is small in solids and is generally neglected in eqn (1). image file: c9cp06485k-t6.tif is the widely used definition of voltages in computational electrochemistry; it is Φ = (EvacuumEF)/|e| modulo a constant, where Evacuum is the vacuum level in the same simulation cell that there Fermi level EF is computed. e is the prime mover in electrochemistry; potentiostats and voltmeters control and measure e energy (EF), not M content.

image file: c9cp06485k-t7.tif is notoriously difficult to define and control in DFT calculations. For liquid electrolytes associated with corrosion, fuel cells, and batteries, AIMD-based thermodynamic integration approaches have been applied to calibrate the voltage in condensed phase simulation cells without vacuum regions.18,22,25,26,55 (At least some of these works18,55 were influenced by the pioneering research of Sprik and coworkers.71,72) Such methods calculate half-cell electrochemical potentials via free energy changes; the predictions are solvent-dependent.

The majority of image file: c9cp06485k-t8.tif calculations rely on having a vacuum region (or a quasi-vacuum containing an implicit solvent) in the simulation cell. Vacuum is in effect the reference electrode. In periodic boundary DFT calculations, both finite temperature AIMD23 and T = 0 K configuration optimization calculations27,31 have used vacuum to compute absolute voltages. The Berry's phase method, under development for liquid electrolytes,30 can potentially circumvent the need for vacuum.

The relation between voltage referenced to Standard Hydrogen Electrode (SHE), and the work function (Φ), was famously quantified by Trassati for liquid electrolytes.61 It is Φ minus 4.44 eV divided by |e|. This relation has been incorporated into or used as calibration for implicit solvent formulations of voltage calculations.22,27,29 In cluster-based quantum chemistry and DFT calculations, where only a solute and sometimes a few solvent molecules are surrounded with implicit solvent, calibration is not possible, and the Trasatti relation is directly applied to relate electron affinities (EA) and ionization potentials (IP) to redox potentials.16,74 When the SHE reference is replaced with Li+/Li(s), the shift is changed from 4.44 eV to between 1.37 and 1.40 eV. The O(0.01 eV) differences are within DFT uncertainties.

A metallic electrode is just a large molecule where EA = IP, with the exception that different facets may have different Φ's which are reconciled by charge transfer and/or geometric considerations.75 The complex organic SEI (Fig. 2a) precludes calibration calculations here. Hence in this and our previous work, we simply assign

image file: c9cp06485k-t9.tif(2)
to electrodes in our overall charge-neutral simulation cells, just like for small molecules, and acknowledge the approximation in omitting the material outside the inorganic solid. Trasatti likely never intended the use of eqn (2) in vacuum, but one can imagine substituting a low dielectric electrolyte, like argon at low density. We also use this relation for pitting corrosion calculations, where water is optimized at T = 0 K. In the future, we will switch to sampling H2O at T = 300 K, which is the more rigorous approach. Unlike ref. 31, which also opens a vacuum gap in the frozen bulk liquid electrolyte region, we include at most a sub-monolayer of water molecules. Like ref. 18, 25 and 55 and unlike ref. 22 and 27, we control image file: c9cp06485k-t10.tif using atoms and e, not with effective medium approaches extrinsic to standard DFT.

C. i–equilibrium

We define “i–equilibrium” to be established when image file: c9cp06485k-t11.tif. Rigorously speaking, i–equilibrium is required for comparison with constant-voltage condition experiments. The exact i–equilibrium condition is affected by the accuracy of the DFT functional used, which is neglected herein. Since the two electrochemical potentials follow the relation image file: c9cp06485k-t12.tif implies a specific [small mu, Greek, macron]M+ in the electrolyte55 which can be used for future calibration calculations. When i–equilibrium is violated, the system is at an overpotential with respect to at least one process. Overpotentials often occur in batteries because M+ and e motions are generally separated by orders of magnitude in time scales. Experimentally, M+ = Li+ relaxation times in cathode materials can be minutes.76 For this reason, we also call image file: c9cp06485k-t13.tif and image file: c9cp06485k-t14.tif the “instantaneous” and “equilibrium” voltages, respectively. In DFT calculations, Li+ relaxation and motion are far more limited than in experiments because the time scale (if using AIMD) is much shorter, and there is no Li reservoir. Therefore DFT modelling of interfaces is even more susceptible to overpotentials. As long as the “electrode” (a large clump of atoms) is metallic, image file: c9cp06485k-t15.tif is instantaneously well defined, even if there are forces on atoms. image file: c9cp06485k-t16.tif is in contrast only well-defined when atomic forces are zero at zero temperature, or when the Boltzmann distribution is satisfied at finite T.

An incorrect assumption often made45 (including in our early work54) is that simulation cells containing an interface is automatically at “open-circuit” voltage equal to image file: c9cp06485k-t17.tif. This ignores the possibility that i-equilibrium can be violated. In fact, image file: c9cp06485k-t18.tif is seldom discussed or reported when modelling batteries. The reason is likely historical. The majority of DFT calculations on electrode materials deal with the crystal interior (i.e., they are “single phase”). In the absence of interfaces, the average electrostatic potential in the simulation cell is undefined to a constant77 except in special cases. Hence the absolute EF, which depends the average potential, is also underfined – as is image file: c9cp06485k-t19.tif via eqn (2). The only logical recourse then is to assume i–equilibrium. When an electrode/electrolyte interface exists in the simulation cell, however, a potential drop occurs that can in principle be computed (Fig. 2). Other reasons interfaces may be out of equilibrium include the presence of high current densities, and electrode polarization whereby a liquid electrolyte responds too slowly to voltage changes. These are not addressed herein.

D. Grand canonical ensemble for e and Li+

Excess e surface densities on pristine electrode surfaces vary during charge and discharge of EDL capacitors, fuel cell electrodes, and other devices. As a result, mobile ions redistribute at interfaces and EDL's are modified. Grand canonical treatments of e, which permits fractional e compensated by an effective medium, have been developed for this purpose.22,29

In batteries, charge carriers M+ are also mobile; their concentrations evolve as voltage varies not only inside battery anodes and cathodes, but also at interfaces where there is often less crystalline order. To model changing M+ content calls for Grand Canonical Monte Carlo (GCMC) or related methods, in conjunction with electronic structure calculations. AIMD simulations, one of the mainstays in battery interfacial calculations, lack this capability. MC is time-independent and can in principle circumvent the orders-of-magnitude solid–liquid time scale mismatch. In the past, DFT/MC calculations were costly because the global DFT energy was recomputed after each single-atom motion in a MC trial.78 GCMC/DFT-like calculations are being implemented;79,80 they are promising for future studies of battery interfaces. Until such methods become widely available, the interfacial M+ content has to be varied manually or via high-throughput calculations, as it has been done in single phase battery electrode calculations.3,4

Many battery cathode materials are polaron-conducting transition metal oxides with finite band gaps. Rigorously speaking, image file: c9cp06485k-t20.tif is undefined such materials (unless surface/interface regions become metallic81); EF is pinned by defect/dopants levels not typically included in the simulation cell, and it cannot be said that all regions of the electrode are at the same voltage. Constant voltage methods which allow fractional e in DFT calculations22,29 cannot be used because fractional occupancy of transition metal d- or f-orbitals is unphysical in insulators. Here one can look to experiments for guidance. In practical batteries, the cathode is a composite with active oxides mixed with conductive components. By adding a metallic slab like lattice-matched Au(001) to the bottom of spinel LixMn2O4, one can recover metal-like EF behavior (Fig. 2b).43 Regarding the structure of cathode surface films, much less is known about thinner, hard-to-characterize CEI on battery cathode surfaces in liquid-electrolyte batteries. Further discussions of cathode interfaces will be left to future overviews.

E. Modelling kinetics and AIMD simulations

At least three general approaches have been applied to model interfacial kinetics in batteries: (1) unconstrained AIMD simulations at or near room temperature; (2) high temperature AIMD; (3) barrier predictions. (1) and (2) are both unconstrained and only differ in intention; the former is meant to directly mimic near-experimental and -device conditions while the latter is explicitly understood to artificially accelerate timescales. Because of the computational expense, current AIMD simulations seldom exceed 1 ns in trajectory lengths.52 So unconstrained AIMD trajectories52–54 are at least 1012 times shorter than experimental time scales. M+-Motion and bond-breaking reactions revealed in such AIMD trajectories are the initial steps. Even if AIMD predicts that the electrolyte is atomized,52 nucleation of these products to SEI solid phases takes more time. Logically, if a change is observed in an AIMD trajectory, it provides important insights about the mechanism. If no change occurs, it is possible AIMD trajectories are just not long enough. In general, AIMD trajectories should exceed relevant molecular relaxation times (e.g., those associated with rotational and intramolecular motion) to qualify as reflecting equilibrium or ‘open-circuit’ conditions.

Alternatively, parasitic reactions at interfaces have been examined using AIMD at substantially elevated temperatures.51 This accelerates reaction rates, but extrapolating predictions to target conditions (usually room temperature) is non-trivial. This approach can be extremely useful for discovering mechanisms that are subsequently quantified via barrier-finding calculations.

Reaction barrier calculations like umbrella sampling and metadynamics bypass trajectory length limits at the cost of having to choose pathways.18,82–84 In purely solid state calculations, zero-temperature nudged-elastic-band (NEB) calculations have become standard.85 From barrier calculations, we estimate mean reaction rates using the standard transition state theory expression

1/tave = ko[thin space (1/6-em)]exp(−ΔE*/kBT), (3)
where ΔE* is the activation energy and kBT is the thermal energy at room temperature. At finite temperatures, free energy barriers (ΔG*) should be used in eqn (3). We adopt ko = 1012 s−1, within a factor of ten of prefactors used in the literature. Any step in a proposed multi-step reaction mechanisms is considered viable if it is exothermic (ΔE < 0) and if tave is less than one hour – which translates into ΔE* < ∼1 eV.15

Rates may depend on voltages, while voltages are modified by atomic motion in finite sized simulation cells. This is because the work function Φ is altered by the effective dipole moment via

ΔΦ = 4πσd. (4)
σ is the surface density of a uniform point dipole sheet, d is the average dipole magnitude and includes screening effects, and atomic units are used. This is an “interface” effect but the ΔΦ magnitude does not decay with distance from the interface. M+ displacement or the transfer of an e following an electrochemical reaction changes d, Φ, and hence image file: c9cp06485k-t21.tif via eqn (4) in finite simulation cells. Under constant voltage conditions, image file: c9cp06485k-t22.tif should be zero over the course of a reaction. To achieve this condition, constant image file: c9cp06485k-t23.tif ensembles86 or extrapolation to large cell sizes87 have been applied in electrocatalysis calculations. Related methods have recently been applied to battery interfaces,88 although applications of such techniques to solid-solid interfaces remain rare.

In some cases, apparent degradation or “disordering” reactions can extend to micron depth beyond the interface89 – far beyond DFT length scales. Fitting reactive force fields can extend the time- and length-scale associated with the battery interface reaction zone.90 However, spin-dependent reactive potentials have not been devised for transition metal oxides. Standard modelling protocols like the “master equation” exist to deal with kinetics-controlled, multistep reactions in catalysis and gas phase reactions.91 So far the complexity of battery reactions has precluded their use.

To avoid challenging prospect of calculating all possible reaction mechanisms, whenever possible, we have focused on elucidating the stability of proposed surface film components.15 We call this the “surface film instability paradigm.” Recent application of artificial intelligence to determine locally optimal interfacial atomic structures48 can yield interfaces stable against bond-breaking. This promising approach may efficiently satisfy our stability criterion if quasi-kinetic constraints can be added.

F. Non-DFT and non-planewave methods

We briefly mention non-DFT electronic structure methods which are not our focus.119 Quantum mechanics/molecular mechanics (QM/MM) methods have accelerated solution phase reactions by confining DFT to a small “QM” spatial region.92,93 Since electrochemical calculations involve voltages which require metallic electrodes, the “QM” required may still be substantial. Tight-binding methods94 have been long proposed for aqueous phase electrochemistry but require more development; tight-binding treatment of organic species and materials associated with battery o-SEI has yet to be developed. We apply a plane-wave plus projector-augmented wave (PAW) basis set, but the discussions above are applicable to all DFT basis sets, including mixed-planewave/atomic basis.78

G. DFT details

All DFT calculations are conducted under T = 0 K ultrahigh vacuum (UHV) condition, using periodically replicated simulation cells and the Vienna Atomic Simulation Package (VASP) version 5.3.95–98 A 400 eV planewave energy cutoff and a 10−4 eV convergence criterion are enforced. Most calculations apply the PBE functional.99 In one case, HSE06 is used as a spot check.100–102 The standard dipole correction is applied.103

We consider Li metal(001)/LiF(001) interfaces as exemplars of Fig. 2b. The base model has a 17.28 × 17.28 × 36 Å3 simulation cell with stoichiometric Li294F144. We conduct a 10 ps AIMD trajectory at T = 400 K to equilibrate the interface, and then optimize the final MD configuration. The bottom two Li metal layers are held fixed throughout. When an O2 molecule is added to the surface, the c lattice constant is expanded by 4 Å. 2 × 2 Brillouin sampling is used. For the purpose of checking PBE results with the HSE06 functional in this model, Γ-point sampling is employed; decreasing the k-point grid changes image file: c9cp06485k-t24.tif by less than 0.12 V.

To make connection with corrosion research, Al(111)/α-Al2O3(0001) interfaces are also considered. No definitive metal–oxide interfacial structure is known in the experimental literature, and proposed DFT-configurations with subtle differences yield significant changes in Φ.104,105 This is unlike the weakly interacting and more robust Li-metal/LiF interface. Like ref. 40, we apply AIMD simulations to equilibrate this interface. Our AIMD is conducted on a primitive surface cell with dimensions 4.81 × 8.33 × 50 Å3 and Al70O54(H2O)6 stoichiometry. This is the size of a single surface unit cell, chosen to be commensurate with those of ref. 104 and 105, so that a direct comparison of the respective properties with those work can be made. We aim to provide more relaxation than the purely T = 0 K configuration relaxation of these two work. Our cell size is smaller than that in ref. 40. We do not claim our procedure yields the definitive structure of this complex interface, which likely depends on material preparation conditions.

Unlike previous DFT work, we add 3 H2O atop each surface Al3+ in an attempt to make surface Al ions 6-coordinated. After a 0.8 ps AIMD trajectory at T = 400 K, the system is cooled to T = 100 K in 0.15 ps. We find that one-third of the added water evaporates/detaches from the surface; they are removed. Another third dissociates into H+ and OH. The atomic configuration is then optimized at T = 0 K, resulting in a configuration where each surface Al is 4-coordinated, and the intact H2O per surface Al lies parallel to the surface, coordinated only to other surface O or OH groups. A 3 × 2 surface supercell is constructed from this primitive cell and is applied to examine oxygen vacancies with 2 × 2 × 1 Brillouin sampling.

The LiOH(001)/Li(111) interface is examined for kinetics purposes. Here the base model has a 13.53 × 14.33 × 40 Å3 simulation cell and a Li224O06H96 stoichiometry, and 2 × 2 × 1 Brillouin sampling is applied. In all cases Li and Al metals, softer than fluorides and oxides, are strained to match the fluoride or oxide supercell lattice constants. This changes the pure metal work function by ∼0.1 V. All simulation cells are charge neutral.

III Results

A. Lithium metal in vacuum

Consider a BCC lithium solid cleaved in the (001) direction in vacuum. The cohesive energy is that of Li-metal, and image file: c9cp06485k-t25.tif vs. Li+/Li as we remove one entire layer of Li atoms, according to eqn (1). (The small strain energy is neglected.) The PBE-predicted work function is Φ = 3.02 eV, close to the experimental value of 2.93 eV.106 From eqn (2), image file: c9cp06485k-t26.tif vs. Li+/Li(s). Hence the system is out of i–equilibrium, at overpotential against Li dissolution. If there were a Li reservoir and an electrolyte, and image file: c9cp06485k-t27.tif were enforced by a potentiostat, the Li slab would completely dissolve.

Fig. 3a reinforces this conclusion. A 8 Å-thick LiF slab with its (001) facet exposed is located 5 Å apart from the Li(001) slab. Neither the Li nor the LiF slab has a net charge or a net dipole moment when isolated. Fig. 4a depicts z-axis-resolved orbital positions showing no electronic overlap between the slabs separated by vacuum. The electrostatic potential profile (Fig. 4c) suggests that the 5 Å-separated slabs remain largely dipole free. image file: c9cp06485k-t28.tif, only weakly perturbed from that of isolated Li-metal slab. Once again, image file: c9cp06485k-t29.tif.

image file: c9cp06485k-f3.tif
Fig. 3 (a) Li(001) and LiF(001) slabs separated by ∼5 Å; (b) Li(001) and LiF(001) slabs in contact; (c) and (d) optimized configurations of panel (a) when O2 is added, and of panel (b) when 9 interlayer Li atoms and a O2 are added. Blue, pink, and red spheres depict Li, F, and O atoms, respectively. The +z direction is to the right.

image file: c9cp06485k-f4.tif
Fig. 4 (a and b) Local densities of state of valence electrons associated with Fig. 3a and b. Red circles depict Kohn–Sham orbitals on atoms which contribute to the wavefunction by more than 2%; the green line are Fermi levels and vacuum is at ΔE = 0.0 eV. (c and d) Electrostatic potentials associated with Fig. 3a and b, respectively; the left and right side vacuum levels correspond to coated and uncoated Li metal, separated by a discontinuity due to the artificial dipole layer.103

This finding cannot be over-emphasized: it is incorrect to assume that Li metal is always at i–equilibrium, i.e., at 0.0 V vs. Li+/Li(s), in an interfacial simulation cell. Experimentally, lithium is stable at or below 0.0 V; when stripping Li under overpotential conditions, Li also exists above 0.0 V. No conclusion about image file: c9cp06485k-t30.tif can be drawn from the existence of Li metal in a simulation cell. The same is true of other electrodes, including intercalation and conversion electrodes. The overpotential associated with plating and stripping of Mg anodes can approach 1 V.19 Au electrodes are featured in numerous DFT interfacial simulation cells.27 Reported voltages image file: c9cp06485k-t31.tif in that literature are correctly pegged to (−EF) modulated by the electric double layer,27 not assumed to be at the Au3+ + 3e → Au(s) half cell voltage. Distinguishing image file: c9cp06485k-t32.tif and image file: c9cp06485k-t33.tif establishes electrochemical equilibrium and permits the modelling of overpotential effects.

B. Surface films strongly affect image file: c9cp06485k-t34.tif

Fig. 3b depicts the configuration obtained after we put the LiF and Li surface slabs in contact. The system develops a large dipole moment, and Φ drops from the vacuum value of 3.02 eV to 1.54 eV. image file: c9cp06485k-t35.tif becomes 0.17 V (eqn (2)), much closer to image file: c9cp06485k-t36.tif than bare Li(001). The surface film brings it close to i–equilibrium – without adding surface charges or liquid electrolyte. Similar, significant “contact potential” drops upon adding an inorganic film, which reduces Φ (and hence image file: c9cp06485k-t37.tif via eqn (2)), have been noted in the corrosion,38 electronics,39 solid electrolyte,43,107 and Li–air battery108 modelling literature. They are partly caused by e density rearranging within each material at the interface, creating a dipole layer (eqn (4)). The reduction in Φ is reminiscent of the effect of a water monolayer on inert metal surfaces,109 except that the LiF film exhibits an induced, not permanent, dipole moment. Significant contact potentials can also be inferred from experiments,50,69 although they are conducted at less than atomic resolution.

We can further reduce image file: c9cp06485k-t38.tif to −0.01 V in Fig. 3b (i–equilibrium condition) by manually adding an interlayer of Li atoms directly underneath some of the F anions at the interface.44 Adding Li atoms in this way is found to be thermoneutral to within 0.05 eV after accounting for Li chemical potential (eqn (1)). It suggests that the energy/voltage landscape at the interface has many similar local minima. This interlayer approach, which reflects the grand canonical ensemble environment, is motivated by the “interfacial charge storage” idea.107 Fig. 4b and d depict the resulting orbital alignment and electrostatic potential. A finite band gap remains in the LiF region, and no LiF surface state is occupied.

Doubling the LiF slab thickness to 18 Å without the Li-interlayer atoms yields a ΔΦ slightly larger by 0.05 eV, leading to image file: c9cp06485k-t39.tif. The somewhat different image file: c9cp06485k-t40.tif partly reflects system size effect, but is likely within the uncertainty margin. The Li2O(111)/Li(001) interface yields a contact potential drop of roughly the same size; however, there is small amount of e leakage to the outer Li2O surface,44 which in reality would have led to further passivation by the organic materials outside. The PBE functional used herein has delocalization errors111 which may permit unphysical e leakage from Li metal into LiF. To check that this does not overestimate the contact potential, we also apply the more accurate HSE06 functional to the Li/8 Å-LiF interface. The resulting ΔΦ is in fact 10% larger than the PBE prediction, showing that our PBE calculations most likely do not exaggerate the contact potential drop.

The flatness of the valence band edge (Fig. 4b) emphasizes that we are at “potential-of-zero-charge”-like conditions. Our focus is the near-anode region, not the voltage profile over long length scales (Fig. 1b). Nevertheless, if in real-life batteries the lithium anode behaves like Fig. 4b, the voltage profile in Fig. 1b would remain flat in the liquid electrolyte region until it reaches the cathode EDL. This ignores the fact that LiF is not the only SEI component. Charges at defect sites and other factors may also introduce electric fields near the anode surface and creates an EDL, part of which may now reside at the solid–solid interface. The exact partitioning of solid state vs. liquid state EDL depends on the free energy cost of developing those EDL's. This little-explored EDL partitioning offers a potential design principle for improving battery interfaces.

C. How i-equilibrium and image file: c9cp06485k-t41.tif affect DFT predictions

image file: c9cp06485k-t42.tif is seldom reported or controlled in DFT modelling of battery interfaces. Some reported DFT interface simulations are likely out of i–equilibrium;45,54 they are at as-synthesized rather than battery cycling conditions with well-defined voltages. What are the consequences?

image file: c9cp06485k-t43.tif is crucial for determining the electrochemical equilibrium conditions where M+ insertion from liquid electrolyte into metallic electrode is not at an overpotential.55,88,94 The properties under these conditions, e.g., M+ insertion kinetics and barrier, is expected to strongly depend on image file: c9cp06485k-t44.tif. image file: c9cp06485k-t45.tif is also crucial for passive electrodes acting as e emitter. The correct image file: c9cp06485k-t46.tif value yields redox event onset at the right voltage. In this light, the edge-termination dependence of LiC6 on electrolyte redox behavior in our early work54 is likely the result of the graphite edge sites giving different image file: c9cp06485k-t47.tif, not computed at that time.

Fig. 3c and d qualitatively illustrate this effect by placing a O2 molecule on to the outer LiF surfaces of Fig. 3a and b. This example draws on DFT modelling of O2 reduction on Al2O3-coated Al surfaces,40 relevant to the growth of a passivating film there, and sidesteps more complex organic molecule redox demonstrations.18,44 When LiF is coated on Li with an “interlayer Li” (Fig. 3b), image file: c9cp06485k-t48.tif, and an e readily transfers through the LiF film to the adsorbed O2. Bader analysis110 confirms that a −1.00|e| net charge exists on the oxygen species, i.e., it is a superoxide (O2). When the Li and LiF slabs are well-separated (Fig. 3a), image file: c9cp06485k-t49.tif is much less reductive (1.65 V), and the O2 weakly physisorbs instead. Magnetic moment analysis suggests that O2 is a triplet molecule. Bader analysis yields a −0.31|e| net charge. The non-integer value is unphysical and is likely due to DFT/PBE localization error111 for this space-separated configuration, which is challenging for DFT methods.

In the literature, it has been reported LiF-coated Li metal rapidly decomposes EC molecules in AIMD simulations.53,112 Fig. 3 suggests the reason is that LiF lowers image file: c9cp06485k-t50.tif to values which initiate electrochemical reduction. On bare Li metal surfaces, image file: c9cp06485k-t51.tif is significantly higher. Despite this, many organic species also rapidly decomposes on bare Li in AIMD simulations.51,52 But these reactions are likely chemical, featuring reactants in contact, rather than electrochemical (i.e., long range e transfer) in nature. We will return to this theme in Section III.F.

image file: c9cp06485k-t52.tif and electric fields associated with it should affect charged carriers transport properties.41,120 Here we use a negative example and examine how Li+ vacancy motion towards the metallic electrode affect computed voltages. Again we consider a Li metal slab as in Fig. 3b, with a LiF layer with twice the thickness (∼18 Å thick, not shown). Without Li-vacancy, image file: c9cp06485k-t53.tif. With a Li-vacancy 10 Å from the lithium metal surface, d increases and image file: c9cp06485k-t54.tif rises to 0.98 V. When the vacancy is 4 Å away from the surface, image file: c9cp06485k-t55.tif, still higher than the no-vacancy configuration. Hence the system is not at constant potential during Li-vacancy transport. This is an artifact of the finite surface area of the simulation cell (eqn (4)). The need to move atoms while keeping image file: c9cp06485k-t56.tif constant is known in electrocatalysis modelling.22,86,87,113 Battery calculations generally have larger lateral surface areas than electrocatalysis and smaller σ (eqn (4)), but d can be much larger unless M+ moves in a high-dielectric liquid electrolyte18,55,88 instead of a relatively low dielectric solid film. In general, if there are large changes in image file: c9cp06485k-t57.tif along a Li diffusion pathway, substantial corrections may be needed to recover the constant potential energy landscape.

D. Relation to atmospheric corrosion

Next we consider contact potential and electric field effects on the occupation of orbitals associated with defects. In battery modelling, related studies have been treated mostly at a flat-band level without explicitly considering electrified interfaces.44,118 Because more prior relevant work exists in the Al2O3 literature, we use the Al(111)/α-Al2O3(0001) interface40 relevant to atmospheric pitting corrosion of aluminum, as example. Of particular interest is the charge state of oxygen vacancies. The “point-defect model” (PDM),37 a widely used approach for predicting passivating oxide breakdown that leads to pitting, proposes doubly positively charged O-vacancies. When Al metal is present in a DFT simulation cell, the charges of defect states are determined by EF and image file: c9cp06485k-t58.tif.

Fig. 5a depicts the orbital alignment of this interface with an oxygen vacancy at position shown in Fig. 5c. By comparing with a slab without vacancies, we identify orbitals circled in green as those associated with the O-vacancy. They lie below the Fermi level EF, indicating that the vacancy is uncharged. Consistent with this, the valence band edge is relatively flat, showing little electric field. image file: c9cp06485k-t59.tif vs. SHE.

image file: c9cp06485k-f5.tif
Fig. 5 (a) Local densities of state of valence electrons associated with the Al(111)/α-Al2O3(0001) interface. The occupied states in the gap which correspond to a neutral oxygen vacancy are circled in green. (b) Same as (a), but with 2 H+ vacancies at the vacuum interface to create an electric field. (c) Atomic configuration of the slab. The green O-atom at z = 32.2 Å is removed to create a vacancy. Yellow spheres are Al atoms. Green line are Fermi levels and vacuum is at 0.0 eV.

To recover positively charged oxygen vacancies, we remove H atoms from two of the twelve surface AlOH groups. As expected, DFT predicts that the H-vacancies carry negative charges, which should be compensated by positive charges on Al metal surfaces in the charge-neutral simulation cell. This raises image file: c9cp06485k-t60.tif to a more oxidizing 0.59 V, and raises the local valence band edge near the surface (Fig. 5b). Now creating the same O-vacancy as before yields unoccupied orbital levels above EF (Fig. 5b), and a doubly positively charged defect which dovetails with PDM37 emerges. A more significant electric field can be inferred from the valence band edge in the range 24 < z < 30 Å.

The O-vacancy is periodically replicated in Fig. 5c. In reality, a spatial distribution of O-vacancies should exist over a range of z values. In pitting corrosion, adsorbed anions are likely to play the electrostatic role of negatively charged H vacancy in Fig. 5c. Adsorbed halide ions are indeed strongly linked with pitting corrosion onset. The image file: c9cp06485k-t61.tif's constructed for these slabs are too high compared with Al pitting onset voltages. Nevertheless, this slab model appears a promising platform for investigating voltage effects in Al2O3 passivating films. Hybrid DFT functionals and other oxide phases should also be considered using this approach. Another fruitful connection between battery and corrosion is the possible galvanic lithium corrosion.114

E. The Li(s)/LiOH interface, a kinetics case study

To illustrate the kinetics-controlled nature of battery interfaces, we consider lithium hydroxide on lithium metal surface. LiOH has occasionally been cited as a SEI component, both experimentally and theoretically.53,62 In ref. 53, AIMD simulations show that organic solvent decomposes on LiOH(010)-coated Li metal(001) surfaces; the LiOH film itself does not react. However, OH groups are expected to give off H2 gas at low voltages.62 Here we explicitly re-examine LiOH stability.

First we consider thermodynamics. At T = 0 K, both

2LiOH(s) + 2Li(s) → 2Li2O(s) + H2(g), (5)
LiOH(s) + 2Li(s) → Li2O(s) + LiH(s) (6)
are exothermic, by −1.74 and −1.70 eV, respectively. The reaction
H2(g) + 2Li(s) → 2LiH(s) (7)
is also exothermic by −1.65 eV. H2 gas release is favored by a ∼0.4 eV entropy gain at room temperature over that at T = 0 K, but this is insufficient to reverse the exothermicity of eqn (7). LiH62 has indeed been observed using cryo-TEM,115 although its presence may depend on electrolyte water content.116 Thermodynamically, LiOH solid is predicted to be unstable against lithium metal. Kinetically, both the LiOH film and isolated OH units, which come from reaction with water, are stable within <100 ps AIMD trajectory lengths.53 Here we turn to activation barrier (ΔE*) calculations to extrapolate to longer time scales.

Fig. 6a depicts an isolated LiOH unit optimized on the Li(001) surface. In Fig. 6b, the proton is detached to the Li surface, with a favorable ΔE = −1.29 eV released. NEB calculations reveal that ΔE* = 0.72 eV. From eqn (3), this translates into an average reaction time of 0.82 s. This is fast on battery time scales, but too slow to be observed in unconstrained AIMD simulations (<10−9 s in duration). We conclude that isolated LiOH is kinetically unstable on Li(001).

image file: c9cp06485k-f6.tif
Fig. 6 (a and b) Intact and decomposed LiOH unit on Li(001). (c) LiOH(010) slab on Li(001) surface focusing on interface region, showing all OH groups there coordinated to surface Li atoms. (d) Adding an extra Li (circled in yellow) at the interface. The arrow indicates the H atom to be transferred. (e) Same as (d), but with a hydrogen transferred to the surface. (f) Top of LiOH(010) slab with removal of a H atom from within yellow circle. The +z direction is to the right.

Next we consider the Li(001)/LiOH(010) interface. Its proposed existence assumes that adsorbed LiOH monomers can nucleate in sub-second time scales before they decompose. The lattice matching is such that all OH groups at the interface are coordinated to surface Li metal atoms (Fig. 6c), unlike LiF and Li2O.44 Breaking one of the three inequivalent OH bonds to form a H on the lithium metal surface releases ΔE = −0.79 eV, favorable for reaction. The barrier is ΔE* = 1.27 eV, which indicates a reaction exceeding experimental time scales. The other two inequivalent OH groups at the interface exhibit ΔE = −0.60 eV but ΔE* are similar at 1.20 eV. If one were to deposit Li vapor on defect-free LiOH(010), reactions are expected only at elevated temperatures.

However, the above calculations ignore the grand canonical nature of the interface during cycling. A metastable configuration can be created by adding a Li atom to the surface (Fig. 6d). ΔE′ is 0.22 eV after accounting for lithium chemical potential. From this configuration, breaking the O–H bond (Fig. 6e) releases ΔE = −1.07 eV and ΔE* = 0.76 eV. Accounting for ΔE′ but not zero point energy (ZPE) that further reduces ΔE* in proton transfer reactions,58 the energy at the transition state is 0.98 eV above that of Fig. 6c. From eqn (4), this is one reaction per 8 hours at T = 300 K. Without mapping out all possible reaction paths via the master equation,91 DFT/PBE calculations predict at least one pathway whereby that a pristine LiOH(010) slab should react with lithium metal on an experimentally observable time scale, likely forming Li2O on Li metal surfaces. We conclude that the LiOH SEI component seen in experiments are likely not in contact with Li metal. It is more likely lodged on other, more stable SEI components such as Li2O.

Thus AIMD simulations of liquid or solid electrolyte decomposition on Li metal surfaces are extremely useful for giving qualitative insights about possible reaction mechanisms, especially at the initial stages of interfacial degradation reactions. But they currently cannot rule out that reactions can occur at >1 ns time scales.

F. Voltage and electric field effects on kinetics

This LiOH(010) example permits an in-depth exploration of the effect of image file: c9cp06485k-t62.tif on ΔE and ΔE* associated with parasitic reactions. Fig. 6c, d, and e are consistent with image file: c9cp06485k-t63.tif, respectively. Unlike the LiF interface (Fig. 3a), all O atoms at the surface are coordinated to Li metal atoms, and attempts at adding an interlayer like in ref. 44 fail to lower image file: c9cp06485k-t64.tif. To achieve i–equilibrium image file: c9cp06485k-t65.tif requires cations outside the LiOH slab to reduce d via eqn (4). Here we take an opposite approach. Motivated by our Al2O3 analysis (Fig. 5), we create a H vacancy at the LiOH outer surface (Fig. 6f). As expected, DFT calculations yield a negatively charged H-vacancy in the neutral cell, which makes d more positive in eqn (4) and increases image file: c9cp06485k-t66.tif from 0.42 V to 2.17 V while image file: c9cp06485k-t67.tif remains 0.0 V. The system is now at a significant overpotential with respect to Li dissolution from the anode. The valence band edge shown in Fig. 7 suggests a significant electric field of ∼0.22 V Å−1.
image file: c9cp06485k-f7.tif
Fig. 7 Local densities of state of valence electrons associated with Fig. 6f with surface H vacancy. For color key, see Fig. 4.

Remarkably, ΔE and ΔE* are almost unchanged. They are now −1.09 eV and +0.77 eV, within 0.02 eV of the values at image file: c9cp06485k-t68.tif. This is despite the fact that Bader charge analysis suggests the H transferred to the surface now has a −2.1|e| charge (i.e., it is a hydride). In electrocatalysis modelling, |e| times changes in image file: c9cp06485k-t69.tif (designated as “U”) is often added post-processing to ΔE on metal surfaces in reaction steps that involve e transfer, and ΔE* have larger field-dependence.60 In Fig. 6, an added image file: c9cp06485k-t70.tif should emerge in ΔE only if a charged species transverse the entire solid state EDL. The lack of change in ΔE* is however unusual, and may reflect strong metal surface screening.

Other parasitic reactions on Li metal we have considered, like those at the Li/LiPON56 (oxy-phosphorus nitride inorganic polymers) and Li/Li2CO315 interfaces, exhibit similarly small voltage dependences until the field strengths imposed by localized defects are sufficient to initiate different degradation mechanisms that involve long-range e transfer.15 Experimentally, it is known that storage (without applied voltage) and cycling conditions (applying voltages) yield different behavior in metal anodes; the dendrite growth in particular depends on charging rate and overpotentials applied.117 If our anecdotal DFT evidence is true in general, it suggests that the difference between storage and cycling conditions does not stem from primary reaction rates, but originates from long range e transfer to defects which may occur in unison with strictly chemical reactions at the interface, and/or Li+ transport effects. At other electrode interfaces, electric field- and voltage-effects remain to be distinguished and elucidated.

IV Conclusions and outlook

In this work, we apply simple examples with explicit model electrode/surface film interfaces to illustrate some key challenges specific to electronic structure DFT modelling of interfaces in high energy density batteries. These include voltage calibration/control, the thin-film covered (“dirty”) nature of electrodes in monovalent cation (M+) batteries where M = Li or Na, and the fact that interfaces are kinetics-controlled. We find that contact potential at metal-solid electrolyte interfaces, which can only be computed with explicit interfaces in DFT models, cannot be neglected in understanding the voltage profile across interfaces. The ionic (“open circuit”) voltage image file: c9cp06485k-t71.tif, which implies equilibrium conditions, and electronic (“instantaneous”) voltage image file: c9cp06485k-t72.tif, must be distinguished. These two quantities determine whether interfaces are at equilibrium, and enable modelling of interfaces at overpotential conditions. image file: c9cp06485k-t73.tif is shown to govern long range transfer and electric fields that modify transport properties. It appears to play little role in chemical reactions where reactants are in direct contact with lithium metals anodes. However, the roles of voltage and electric fields on reactions at other battery interfaces need further clarification.

The surface-film covered nature of electrodes is not an inconvenience; it is a major part of the interfacial physics of batteries, and should be regarded as a new basic science paradigm and opportunity. In contrast, modelling pristine metal electrode surfaces of M+ batteries, is mainly relevant to the initial stages of battery assembly. Unconstrained ab initio molecular dynamics (AIMD) simulations cannot conclusively prove the stability of interfacial components because trajectory lengths are 12 orders of magnitude shorter than battery operational time scales; barrier-finding or related computational techniques are needed to understand kinetics-controlled interfaces. The M+ content across the entire interface follows the grand canonical ensemble. Cathode interfaces are challenging because of the polaron-conducting nature of cathode materials. Borrowing from other branches of computational electrochemistry like corrosion, supercapacitors, and electrocatalysis presents good opportunities for advancing the current state of battery interface modelling.

We have avoided exploring a wide variety of battery electrode/electrolyte materials and battery concepts such as lithium–sulfur (Li–S), lithium–oxygen (Li–O2), and conversion cathode materials on purpose. The latter areas have been subjects of modelling overviews, while a critical review of modelling methods has been lacking. The unifying modelling concepts discussed in this work focus on metallic conductor electrodes and reasonably well-characterized surface films or solid electrolytes. These should be broadly applicable to all battery electrodes like Na anodes, lithiated graphite, and lithiated silicon. One exception is multivalent (e.g., Mg) batteries, where surface films have been avoided to prevent slow multivalent ion transport. Li–S battery apply Li-metal cathodes and often carbon-based cathodes that trap sulfur; these electrodes are metallic and contain surface CEI films too. Carbon cathode/Li2O2 interfaces in Li–air batteries also fit this category. All-solid state with metallic anodes would be prototype systems to apply these computational ideas, but the thicker solid electrolyte compared to surface films in liquid electrolyte batteries may require special treatment (e.g., those devised for DFT modelling of Mott–Schottky contacts) to address the larger space charge region. Computing overpotential effects on kinetics in kinetics-limited conversion cathode materials may be a particularly fruitful area to apply our approach; the calculations of reaction barriers at non-overpotential conditions require careful turning of the electronic voltage at the interface. Interfaces between transition metal oxide cathodes and liquid electrolytes may present the most significant challenges. While the non-metallic nature in many such oxides can be circumvented with by adding a computational “current collector” at the bottom, the very thin surface film (“CEI”) covering oxide surfaces, which has mostly unknown atomic structures, hinders assignment of voltages and predictions of degradation and transport reaction kinetics. In the final analysis, at this and all other battery interfaces, the voltage is an imposed, DFT-constructed quantify; the question is whether model interfacial structures that provide the imposed voltages yield electric fields that reflect experimental conditions.

Conflicts of interest

There are no conflicts to declare.


We thank Mira Todorova, Minoru Otani, Axel Gross, Harald Oberhofer, and De-en Jiang for discussions during the 2019 Telluride Workshop on Electrochemical Interfaces. We also thank Quinn Campbell and Catalin Spataru for commenting on the manuscript. Above all, we gratefully acknowledge Michiel Sprik's pioneering research. The Li/LiOH work is funded by Nanostructures for Electrical Energy Storage (NEES), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DESC0001160. The Li/LiF work is funded by the Laboratory Directed Research and Development Program at Sandia National Laboratories. The Al/Al2O3 work is funded by the Advanced Strategic Computing (ASC) Program. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the document do not necessarily represent the views of the U.S. Department of Energy or the United States Government.


  1. J. Liu, et al., Nat. Energy, 2017, 4, 180–186 CrossRef.
  2. Basic Research Needs for Next Generation Electrical Energy Storage, http://science.energy.gov/community/resources/report.
  3. M. K. Aydinol, A. F. Kohan, G. Ceder, K. Cho and J. Joannopoulos, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 1354–1365 CrossRef CAS.
  4. G. Ceder, G. Hautier, A. Jain and S. P. Ong, MRS Bull., 2011, 36, 185–191 CrossRef CAS.
  5. J. A. Dawson, P. Canepa, T. Famprikis, C. Masquelier and M. S. Islam, J. Am. Chem. Soc., 2018, 140, 362–368 CrossRef CAS PubMed.
  6. P. Canepa, J. A. Dawson, G. S. Gautam, J. M. Statham, S. C. Parker and M. S. Islam, Chem. Mater., 2018, 30, 3019–3027 CrossRef CAS.
  7. T. Famprikis, P. Canepa, J. A. Dawson, M. S. Islam and C. Masquelier, Nat. Mater., 2019, 18, 1278–1291 CrossRef CAS PubMed.
  8. Y. Zhu, X. He and Y. Mo, ACS Appl. Mater. Interfaces, 2015, 7, 23685–23693 CrossRef CAS PubMed.
  9. A. Baskin and D. Prendergast, J. Electrochem. Soc., 2017, 164, E3834–E3447 CrossRef.
  10. See A. Wang, S. Kadam, H. Li, S. Shi and Y. Qi, npj Comput. Mater., 2018, 4, 15 CrossRef , for a comprehensive review.
  11. O. Borodin, X. Ren, J. Vatamanu, A. v. W. Cresce, J. Knap and K. Xu, Acc. Chem. Res., 2017, 50, 2886–2894 CrossRef CAS PubMed.
  12. K. T. Butler, G. S. Gautam and P. Canepa, npj Comput. Mater., 2019, 5, 19 CrossRef.
  13. K. Xu, Chem. Rev., 2014, 114, 11503–11618 CrossRef CAS PubMed.
  14. Here “electrochemical stability” is defined as the empirically-determined voltage range over which interfaces are stable. It involves a mixture of thermodynamics, kinetics, and electron transfer considerations15.
  15. K. Leung, F. Soto, K. Hankins, P. B. Balbuena and K. L. Harrison, J. Phys. Chem. C, 2016, 120, 6302–6313 CrossRef CAS.
  16. L. Xing, O. Borodin, G. D. Smith and W. Li, J. Phys. Chem. A, 2011, 115, 13896–13905 CrossRef CAS PubMed.
  17. O. Borodin, M. Olguin, C. E. Spear, K. W. Leiter and J. Knap, Nanotechnology, 2015, 26, 354003 CrossRef PubMed.
  18. K. Leung, Phys. Chem. Chem. Phys., 2015, 17, 1637–1643 RSC.
  19. Exceptions include metal anodes for divalent cations (e.g., Mg2+) which must avoid solid-state SEI that drastically slow down cation transport. see, e.g., H. D. Yoo, I. Shterenberg, Y. Gofer, G. Gershinsky, N. Pour and D. Aurbach, Energy Environ. Sci., 2013, 6, 2265–2279 RSC.
  20. K. Persson, Y. Hinuma, Y. S. Meng, A. Van der Ven and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 125416 CrossRef.
  21. K. Leung, J. Phys. Chem. C, 2013, 117, 1539–1547 CrossRef CAS.
  22. J. Haruyama, T. Ikeshoji and M. Otani, Phys. Rev. Mater., 2018, 2, 095801 CrossRef.
  23. N. Bonnet, T. Morishita, O. Sugino and M. Otani, Phys. Rev. Lett., 2012, 109, 266101 CrossRef PubMed.
  24. I. Hamada, O. Sugino, N. Bonnet and M. Otani, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 155427 CrossRef.
  25. S. Surendralal, M. Todorova, M. W. Finnis and J. Neugebauer, Phys. Rev. Lett., 2018, 120, 246801 CrossRef CAS PubMed.
  26. M. Todorova and J. Neugebauer, Phys. Rev. Appl., 2014, 1, 014001 CrossRef.
  27. K. Letchworth-Weaver and T. A. Arias, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 075140 CrossRef.
  28. Q. Cambell and I. Dabo, Phys. Rev. B, 2017, 96, 205308 CrossRef.
  29. N. G. Hörmann, O. Andreussi and N. Marzari, J. Chem. Phys., 2019, 150, 041730 CrossRef PubMed.
  30. T. Dufils, G. Jeanmairet, B. Rotenberg, M. Sprik and M. Salanne, Phys. Rev. Lett., 2019, 123, 195501 CrossRef CAS PubMed.
  31. C. D. Taylor, S. A. Wasileski, J.-S. Filhol and M. Neurock, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 165402 CrossRef.
  32. M. Otani, I. Hamada, O. Sugino, Y. Morikawa, Y. Okamoto and T. Ikeshoji, J. Phys. Soc. Jpn., 2008, 77, 024802 CrossRef.
  33. J. Cheng and M. Sprik, Phys. Chem. Chem. Phys., 2012, 14, 11245–11267 RSC.
  34. C. Zhang, J. Hutter and M. Sprik, J. Phys. Chem. Lett., 2019, 10, 3871–3876 CrossRef CAS PubMed.
  35. B. Wen, W. J. Yin, A. Selloni and L. M. Liu, J. Phys. Chem. Lett., 2018, 9, 5281–5287 CrossRef CAS PubMed.
  36. T. A. Pham, Y. Ping and G. Galli, Nat. Mater., 2017, 16, 401–408 CrossRef CAS PubMed.
  37. C. Y. Chao, L. F. Lin and D. D. Macdonald, J. Electrochem. Soc., 1981, 128, 1187–1194 CrossRef CAS.
  38. X.-X. Yu and L. Marks, Corrosion, 2019, 75, 152–166 CrossRef CAS.
  39. M. Kondo and T. Matsushita, ACS Omega, 2019, 4, 13426–13434 CrossRef CAS PubMed.
  40. D. Costa, T. Ribeiro, F. Mercuri, G. Pacchioni and P. Marcus, Adv. Mater. Interfaces, 2014, 1, 1300072 CrossRef.
  41. K. C. Santosh, R. C. Longo, K. Xiong and K. Cho, J. Electrochem. Soc., 2014, 161, F3104–F3110 CrossRef.
  42. N. D. Lepley and N. A. W. Holzwarth, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 214201 CrossRef.
  43. K. Leung and A. Leenheer, J. Phys. Chem. C, 2015, 119, 10234–10246 CrossRef CAS.
  44. K. Leung and K. L. Jungjohann, J. Phys. Chem. C, 2018, 121, 20188–20196 CrossRef.
  45. H. M. Tang, Z. Deng, Z. N. Lin, Z. B. Wang, I. H. Chu, C. Chen, Z. Y. Zhu, C. Zheng and S. P. Ong, Chem. Mater., 2018, 30, 163–173 CrossRef CAS.
  46. J. Haruyama, K. Sodeyama, L. Han, K. Takada and Y. Tateyama, Chem. Mater., 2014, 26, 4248–4255 CrossRef CAS.
  47. J. Haruyama, K. Sodeyama and Y. Tateyama, ACS Appl. Mater. Interfaces, 2017, 9, 286–292 CrossRef CAS PubMed.
  48. B. Gao, R. Jalem, Y. Ma and Y. Tateyama, Chem. Mater., 2020, 32, 85–96 CrossRef CAS.
  49. M. Sumita, Y. Tanaka, M. Ikeda and T. Ohno, J. Phys. Chem. C, 2015, 119, 14–22 CrossRef CAS.
  50. K. Yamamoto, Y. Iriyama, T. Asaka, T. Hirayama, H. Fujita, C. A. J. Fisher, K. Nonaka, Y. Sugita and Z. Ogumi, Angew. Chem., Int. Ed., 2010, 49, 4414–4417 CrossRef CAS PubMed.
  51. H. Yildirim, J. B. Haskins, C. W. Bauschlicher and J. W. Lawson, J. Phys. Chem. C, 2017, 121, 28214–28234 CrossRef CAS.
  52. B. V. Merinov, S. V. Zybin, S. Naserifar, S. Morozov, J. Oppenheim, W. A. Goddard, J. Lee, J. H. Lee, H. E. Han, Y. C. Choi and S. H. Kim, J. Phys. Chem. Lett., 2019, 10, 4577–4586 CrossRef CAS PubMed.
  53. E. P. Kamphaus, S. Angarita-Gomez, X. Qin, M. Shao, M. Engelhard, K. T. Mueller, V. Murugesan and P. B. Balbuena, ACS Adv. Mater. Interfaces, 2019, 11, 31467–31476 CrossRef CAS PubMed.
  54. K. Leung and J. L. Budzien, Phys. Chem. Chem. Phys., 2010, 12, 6583–6586 RSC.
  55. K. Leung and C. M. Tenney, J. Phys. Chem. C, 2013, 117, 24224–24235 CrossRef CAS.
  56. K. Leung, A. J. Pearse, A. A. Talin, E. J. Fuller, G. W. Rubloff and N. A. Modine, ChemSusChem, 2018, 11, 1956–1969 CrossRef CAS PubMed.
  57. S.-H. Yoo, M. Todorova and J. Beugebauer, Phys. Rev. Lett., 2018, 120, 066101 CrossRef CAS.
  58. M. E. Tuckerman, D. Marx and M. Parrinello, Nature, 2002, 417, 925–929 CrossRef CAS PubMed.
  59. M. Fu, Z. Yao, X. Ma, H. Dong, K. Sun, S. Hwang, E. Hu, H. Gan, Y. Yao, E. A. Stach, C. Wolverton and D. Su, Nano Energy, 2019, 63, 103882 CrossRef CAS.
  60. K. Chan and J. K. Norskov, J. Phys. Chem. Lett., 2015, 14, 2663–2668 CrossRef.
  61. S. Trasatti, J. Electroanal. Chem., 1986, 209, 417–428 CrossRef.
  62. A. Schechter, D. Aurbach and H. Cohen, Langmuir, 1999, 15, 3334–3342 CrossRef CAS.
  63. L. Wang, A. Menakath, F. D. Han, Y. Wang, P. Y. Zavalij, K. J. Gaskell, O. Borodin, D. Iuga, S. P. Brown, C. S. Wang, K. Xu and B. W. Eichhorn, Nat. Chem., 2019, 11, 789–796 CrossRef CAS PubMed.
  64. R. Jorn, R. Kumar, D. P. Abraham and G. A. Voth, J. Phys. Chem. C, 2013, 117, 3747–3761 CrossRef CAS.
  65. M. J. Boyer and G. S. Hwang, Phys. Chem. Chem. Phys., 2019, 40, 22449–22455 RSC.
  66. Z. Zhuo, P. Lu, C. Delacourt, R. Qiao, K. Xu, F. Pan, S. J. Harris and W. Yang, Chem. Commun., 2018, 54, 814–817 RSC.
  67. G. M. Veith, M. Doucet, J. K. Baldwin, R. L. Sacci, T. M. Fears, Y. Wang and J. F. Browning, J. Phys. Chem. C, 2015, 119, 20339–20349 CrossRef CAS.
  68. S. Ringe, H. Oberhofer and K. Reuter, J. Chem. Phys., 2017, 146, 134103 CrossRef.
  69. H. Musada, N. Ishida, Y. Ogata, D. Ito and D. Fujita, Nanoscale, 2017, 9, 893–898 RSC.
  70. E. J. Cheng, A. Sharafi and J. Sakamoto, Electrochim. Acta, 2017, 223, 85–91 CrossRef CAS.
  71. C. Adriaanse, J. Cheng, V. Chau, M. Sulpizi, J. VandeVondele and M. Sprik, J. Phys. Chem. Lett., 2012, 3, 3411–3415 CrossRef CAS PubMed.
  72. F. Costanzo, M. Sulpizi, R. G. Della Valle and M. Sprik, J. Chem. Theory Comput., 2008, 4, 1049 CrossRef CAS PubMed.
  73. In contrast, multivalent batteries with M2+ and M3+ charge carriers must use electrolyte that does not form surface films, because these ions moves slowly through solid films. These battery concepts requires substantial development.
  74. P. Peljo and H. H. Girault, Energy Environ. Sci., 2018, 11, 2306–2309 RSC.
  75. N. W. Ascroft and N. D. Mermin, Solid State Physics, CBS Publishing, Tokyo, 1976 Search PubMed.
  76. A. Barai, K. Uddin, M. Dubarry, L. Somerville, A. McGordon, P. Jennings and I. Bloom, Prog. Energy Combust. Sci., 2019, 72, 1–31 CrossRef , and references therein.
  77. S. W. de Leeuw, J. W. Parram and E. R. Smith, Proc. R. Soc. London, 1980, 373, 27–56 CAS.
  78. M. J. McGrath, J. I. Siepmann, I. F. W. Kuo, C. J. Mundy, J. VandeVondele, J. Hutter, F. Mohamed and M. Krack, ChemPhysChem, 2005, 6, 1894–1901 CrossRef CAS PubMed.
  79. R. B. Wexler, T. Qiu and A. M. Rappe, J. Phys. Chem. C, 2019, 123, 2321–2328 CrossRef CAS.
  80. S. Kasamatsu and O. Sugino, J. Phys.: Condens. Matter, 2019, 31, 085901 CrossRef PubMed.
  81. I. Scivetti and G. Teobaldi, J. Phys. Chem. C, 2015, 119, 21358–21368 CrossRef CAS.
  82. See, e.g., Z. Jiang, K. Klyukin and V. Alexandrov, Phys. Chem. Chem. Phys., 2017, 19, 14897–14901 RSC.
  83. Y. Okuno, K. Ushirogata, K. Sodeyama, G. Shukri and T. Tateyama, J. Phys. Chem. C, 2019, 123, 2267–2277 CrossRef CAS.
  84. R. Benedek, J. Phys. Chem. C, 2017, 121, 22049–22053 CrossRef CAS.
  85. G. Henkelman, B. P. Uberuaga and H. Jonsson, J. Chem. Phys., 2000, 113, 9901–9904 CrossRef CAS.
  86. G. Kastlunger, P. Lindgren and A. A. Peterson, J. Phys. Chem. C, 2018, 122, 12771–12781 CrossRef CAS.
  87. M. Van den Bossche, E. Skulason, C. Rose-Petruck and H. Jonsson, J. Phys. Chem. C, 2019, 123, 4116–4124 CrossRef CAS.
  88. T. Ohwaki, T. Ozaki, Y. Okuno, T. Ikeshoji, H. Imai and M. Otani, Phys. Chem. Chem. Phys., 2018, 20, 11586–11591 RSC.
  89. Z. Wang, D. Santhanagopalan, W. Zhang, F. Wang, H. L. Xin, K. He, J. Li, N. J. Dudney and Y. S. Meng, Nano Lett., 2016, 16, 3760–3767 CrossRef CAS PubMed.
  90. D. Bedrov, G. D. Smith and A. C. T. Van Duin, J. Phys. Chem. A, 2012, 116, 2978–2985 CrossRef CAS.
  91. X. H. Li, A. W. Jasper, J. Zador, J. A. Miller and S. J. Klippenstein, Proc. Combust. Inst., 2017, 36, 219–227 CrossRef CAS.
  92. S. Dohm, E. Spohr and M. Korth, J. Comput. Chem., 2017, 38, 51–58 CrossRef CAS.
  93. T. Fujie, N. Takenaka, Y. Suzuki and N. Nagaoka, J. Chem. Phys., 2018, 149, 044113 CrossRef PubMed.
  94. Y. Li and Y. Qi, Energy Environ. Sci., 2019, 12, 1286–1295 RSC.
  95. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169 CrossRef CAS PubMed.
  96. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50 CrossRef CAS.
  97. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775 CrossRef CAS.
  98. J. Paier, M. Marsman and G. Kresse, J. Chem. Phys., 2007, 127, 024103 CrossRef PubMed.
  99. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868 CrossRef CAS PubMed.
  100. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215 CrossRef CAS.
  101. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2006, 124, 219906 CrossRef.
  102. O. A. Vydrov, J. Heyd, A. V. Krukau and G. E. Scuseria, J. Chem. Phys., 2006, 125, 074106 CrossRef PubMed.
  103. J. Neugebauer and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 16067–16080 CrossRef CAS PubMed.
  104. E.-G. Kim and J.-L. Bredas, Org. Electron., 2013, 14, 569–574 CrossRef CAS.
  105. X.-G. Wang, J. R. Smith and A. Evans, Phys. Rev. Lett., 2002, 89, 286102 CrossRef PubMed.
  106. Recall the DFT/PBE result quoted is for strained lithium, in CRC Handbook for Chemistry and Physics, CRC Press, 2012, ch. 12 Search PubMed.
  107. L. Fu, C.-C. Chen, D. Samuelis and J. Maier, Phys. Rev. Lett., 2014, 112, 208301 CrossRef.
  108. A. C. Luntz, V. Viswanathan, J. Voss, J. B. Varley, J. K. Norskov, R. Scheffler and A. Speidel, J. Phys. Chem. Lett., 2013, 4, 3494–3499 CrossRef CAS.
  109. S. Schnur and A. Gross, Catal. Today, 2011, 165, 129–137 CrossRef CAS.
  110. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354–360 CrossRef , note that charge decomposition schemes like Bader are approximate.
  111. A. J. Cohen, P. Mori-Sanchez and W. T. Yang, Science, 2008, 321, 792–794 CrossRef CAS.
  112. J. Alvarado, M. A. Schroeder, M. H. Zhang, O. Borodin, E. Gobrogge, M. Olguin, M. S. Ding, M. Gobet, S. Greenbaum, Y. S. Meng and K. Xu, Mater. Today, 2018, 21, 341–353 CrossRef CAS.
  113. M. E. Bjorketun, Z. H. Zeng, R. Ahmed, V. Tripkovic, K. S. Thygesen and J. Rossmeisl, Chem. Phys. Lett., 2013, 555, 145–148 CrossRef CAS.
  114. D. Lin, Y. Liu, Y. Li, Y. Li, A. Pei, J. Xie, W. Huahng and Y. Cui, Nat. Chem., 2019, 11, 382–389 CrossRef CAS PubMed.
  115. M. J. Zachman, Z. Tu, S. Choudhury, L. A. Archer and L. F. Kourkoutis, Nature, 2018, 560, 345–351 CrossRef CAS.
  116. C. Fang, J. Li, M. Zhang, F. Yang, J. Z. Lee, M.-H. Lee, J. Alvarado, M. A. Schroeder, Y. Yang, B. Lu, Ni. Williams, M. Ceja, L. Yang, M. Cai, J. Gu, K. Xu, X. Wang and Y. S. Meng, Nature, 2019, 572, 511–515 CrossRef CAS PubMed.
  117. M. Jackle, K. Helmbrecht, M. Smits, D. Stottmeister and A. Gross, Energy Environ. Sci., 2018, 11, 3400–3407 RSC.
  118. H. Tian, Z. Liu, Y. Z. Ji, L. Q. Chen and Y. Qi, Chem. Mater., 2019, 31, 7351–7359 CrossRef CAS.
  119. A. A. Franco, A. Ricci, D. Brandell, C. Frayret, M. Gaberscek, P. Jankowski and P. Johansson, Chem. Rev., 2019, 119, 4569–4627 CrossRef CAS PubMed.
  120. R. C. Longo, L. E. Camacho-Forero and P. B. Balbuena, J. Mater. Chem. A, 2019, 7, 8527–8539 RSC.

This journal is © the Owner Societies 2020