Electrochemical modeling of silicon in lithium-ion batteries using a multi-species, multi-reaction framework with atomistic insights

Abstract

Silicon is a promising anode material for lithium-ion batteries due to its high capacity and potential for fast charging. However, its electrochemical behavior is dominated by pronounced voltage hysteresis, particle-size-dependent voltage plateaus, and relaxation processes induced by hysteresis. Conventional Doyle–Fuller–Newman models cannot capture these phenomena. Here, we present a multi-species, multi-reaction framework that explicitly considers the multiphase lithium–silicon system by assigning an independent equilibrium potential to each phase, derived from modified Nernst equations and parameterized with experimental and atomistic data. The model captures both asymmetric lithiation and delithiation pathways as well as phase-fraction evolution in silicon half-cells. Quantitative comparison yields root-mean-square errors of 5.4–36.9 mV during constant-current and pulse protocols, corresponding to a relative root mean square error of 0.6–4.1% of the overall voltage window. Simulations further reveal that phase fractions continue to evolve during relaxation through thermodynamic redistribution of lithium between phases, governed by phase-specific equilibrium potentials and kinetics. This cross-dimensional approach enables a mechanistic representation of voltage hysteresis, providing a pathway toward improved state estimation, cell design, and battery management.

---

Broader context

Decarbonizing transportation and energy systems is essential to mitigating the accelerating impacts of climate change, including extreme weather events and rising global temperatures. Lithium-ion batteries are central to this transformation, but conventional graphite anodes limit their energy density. Silicon, with a theoretical capacity ten times higher and abundant in nature, is a key factor for next-generation batteries that power electric vehicles and renewable energy storage. However, its practical implementation is hampered by pronounced voltage hysteresis, complex phase transformations, and long relaxation times that remain poorly understood and difficult to model. Here, we introduce a multi-species, multi-reaction (MSMR) modeling framework that directly couples atomistic energetics with macroscopic electrochemistry, capturing the multiphase behavior of the Li–Si system. This approach reproduces with high quantitative accuracy the asymmetric lithiation and delithiation pathways, voltage hysteresis, and relaxation dynamics. By linking material-scale thermodynamics to electrode-scale performance, the model provides a physically coherent basis for optimizing Si-based anodes and contributes to the development of high-energy-density batteries.

1. Introduction

The growing demand for sustainable mobility and energy storage drives the development of lithium-ion batteries with higher energy densities and faster charging capabilities. Silicon (Si) is a highly promising anode material due to its theoretical capacity of 3579 mAh g⁻¹, abundance, and environmental benefits, providing nearly ten times the capacity of graphite.^1–5 Its higher average potential versus lithium further supports fast charging by allowing greater overpotentials.⁶ However, practical application is hindered by large volume changes during lithiation and delithiation, which induce mechanical stress, particle fracture, and electrode swelling.^5–9

Equally critical are electrochemical phenomena that remain poorly understood, including voltage hysteresis, extended relaxation times, and complex phase transformations. Voltage hysteresis in Si reaches about 200 mV, considerably larger than in graphite or lithium iron phosphate.^10–12 Its origins have been linked to multiple factors. Dreyer et al. emphasized the thermodynamic contributions of phase coexistence and distinct equilibrium potentials in two-phase systems such as lithium iron phosphate.¹¹ Bazant et al. advanced this understanding by developing a generalized kinetic theory based on nonequilibrium thermodynamics, introducing the concept of mosaic instability to describe multiparticle phase transformations.^13–15 Köbbing et al. proposed that mechanical effects, particularly the viscoelastic response of the solid electrolyte interphase, can explain slow voltage relaxation and hysteresis in Si nanoparticles.^16,17 Together, these perspectives underline the complex interplay of thermodynamics, kinetics, and mechanics in shaping the electrochemical response of Si.

Another line of explanation focuses on the multiphase nature of the Li–Si alloy system. Baker and Verbrugge introduced multi-species multi-reaction models (MSMR), which assign independent equilibrium potentials to each phase using modified Nernst equations.^18–21 Jiang et al. extended this framework to include Si-specific crystallization mechanisms.³ However, a spatially resolved electrochemical model capturing both phase coexistence and relaxation dynamics remains absent. Moreover, uncertainties regarding the dominant phases, crystalline-to-amorphous transformations, and the associated chemical potentials continue to hinder reliable parametrization. However, experimental studies on the Li–Si system provide essential benchmarks. Wen and Sharma identified multiple equilibrium voltage plateaus in Li–Si alloys at high temperature.^22,23 Weydanz et al. demonstrated that even at very low currents, the potential of electrochemical (de)lithiation deviates significantly from equilibrium.²⁴ Atomistic simulations, including density functional theory, molecular dynamics, Minima Hopping Method (MHM), and CALPHAD phase diagram calculations, confirm the intermetallic nature of the system and predict both crystalline and amorphous phases.^25–31 These findings call into question the direct transferability of modeling frameworks originally designed for solid-solution and systems undergoing spinodal decomposition that rely on phase-field descriptions via Cahn-Hilliard and similar equations.^32,33

From a practical perspective, voltage hysteresis complicates the determination of the true open-circuit potential, alters the energy balance, and increases heat generation.³⁴ To address these issues, empirical approaches such as Plett or Preisach-type hysteresis models have been proposed.^12,35,36 While these methods can reproduce hysteresis loops phenomenologically, they lack direct electrochemical significance and fail to capture particle-size effects, cutoff-voltage dependence, and state-of-charge-dependent relaxation.^37–39

A modeling framework that explicitly accounts for the phase diversity of Si offers the opportunity to overcome these limitations. By integrating experimental data with atomistic insights, such an approach can reproduce voltage hysteresis and associated relaxation phenomena while remaining grounded in physical parameters. Although it does not answer all open questions, this study deepens the understanding of Si (de)lithiation and provides a robust foundation for accurate state-of-charge and state-of-health estimation in advanced battery systems.

2. Fundamentals

2.1. Reactions in electrochemical models for porous electrodes

The conventional electrochemical description of lithium-ion batteries, as formulated in the Doyle–Fuller–Newman (DFN) model, treats the active material as being lithiated or delithiated by lithium ions originating from the liquid electrolyte, as expressed in eqn (1) (ref. 40):

AM + Li⁺ + e⁻ ⇋ AMLi (1)

Eqn (1) represents the intercalation or alloying of lithium into the active material (AM), with the kinetics of this process governed by the exchange current density j.

In eqn (2), ϕ_s denotes the potential in the solid phase and ϕ_l the potential in the electrolyte, with their difference defining the overpotential at the surface of the solid phase (η = ϕ_s − ϕ_l). The first exponential term in the Butler–Volmer expression corresponds to the anodic contribution, while the second represents the cathodic contribution. The parameters α and c are the charge transfer coefficients, describing how the activation energy barrier is distributed between anodic and cathodic reactions on the molecular level, thereby governing the sensitivity of the current density to changes in overpotential.

The reaction rate constant k₀ in eqn (3) determines the equilibrium exchange current density, which quantifies the rate of the reaction described in eqn (1) under equilibrium conditions. Its value depends on the concentrations of reactants and products, as expressed in eqn (3):

Here c_AM and c_Li⁺ denote the concentration of reactants and c_AMLi the concentration of products in the direction of lithiation. Since the direct determination of the rate constant k₀ is challenging, a representative parameter j₀ is used in applied electrochemical models to calculate the exchange current density. To capture changes in the equilibrium potential (E) during (de)lithiation, the state of lithiation (SoL) is typically introduced as shown in eqn (4) in conventional models.

In this relation, the open-circuit potential (OCP) lookup table defines the equilibrium potential of the active material vs. Li/Li⁺. Here c_max denotes the maximum lithium concentration that the material can accommodate, while c_x corresponds to the current lithium-ion concentration in the active material. The OCP function itself is typically obtained from experimental datasets recorded under slow (e.g., C/40) (de)lithiation, ensuring conditions close to equilibrium. A key distinction of the electrochemical model presented in this study is that the equilibrium potential is not imposed as an empirical SoC-voltage lookup table but is consistently calculated from the underlying phase behavior.

2.2. Mass transport and current densities

Electrochemical models not only describe interfacial reactions but also account for the transport of species and the spatial distribution of current within the electrode. Accordingly, the electronic current density in the solid conductive components of the electrode, including the active material and the conductive additives, is described by Ohm's law (eqn (5)):

i_s = −σ_s∇ϕ_s (5)

In applied models (σ_s) is treated as an effective conductivity of the porous electrode, reflecting the combined contributions of active material, conductive additive, and binder.

In the electrolyte phase, where lithium ions act as the primary charge carriers, the ionic current density is described by eqn (6):

In eqn (6), the first term accounts for migration, while the second term accounts for the diffusion contribution to the current density.

The current density interrelates with the mass transport of lithium-ions, according to Fuller et al. (eqn (7)).⁴⁰

At the interface of active material and electrolyte, the mass flow of lithium ions is described according to eqn (8)

The flux of lithium-ions entering or leaving the geometry is also dependent on the surface area, which is why eqn (8) contains the multiplication by the geometry-specific surface area to obtain the molar flow rate. In terms of ideal spherical particles or phase volumes, ε_s,eff. stands for volume fraction, which, in combination with the prefactor three and the particle radius r_s corresponds to the effective specific surface area. Furthermore, in eqn (8), n stands for the number of electrons participating, while ν designates the stoichiometric coefficients in the reaction.

Inside the solid particles, the transport of lithium ions is considered in the conventional DFN model to be purely concentration-dependent and potential-independent. Therefore, the intraparticle mass transport is given by Fick's diffusion eqn (9):

N_s = −D_s∇c_s. (9)

2.3. Thermodynamic equilibrium of the electrochemical reaction

From an electrochemical standpoint, a macroscopic reaction between electrochemical species, that is, lithium-ions and SiLi^x (x denotes the corresponding SiLi phase, as indicated by roman numerals in Fig. 1A), can be traced back to the dynamic behavior of a redox couple (eqn (10)):

Ox + e⁻ ⇋ red (10)

Fig. 1 (A) Binary Li–Si phase diagram (adapted from Braga et al.) (ref. 25), highlighting the intermetallic nature of the system with distinct phase regions (I–V) and sharp stoichiometric transitions. (B) Equilibrium potentials of Li–Si alloys as a function of lithiation degree, compiled from experimental measurements (Sharma & Seefurth (ref. 22); Wen & Huggins (ref. 23)) and theoretical studies (Chevrier & Dahn (ref. 30); Valencia-Jaime et al. (ref. 29); Olou'ou Guifo et al. (ref. 26)). Thin lines represent OCPs from literature data, while dots represent OCPs obtained from atomistic simulations in this study. The thick lines correspond to the OCP functions implemented in the constant-current electrochemical simulations. Arrows indicate the stoichiometric composition corresponding to the atomistic structures shown in (C) and (D). (C) Atomistic model of crystalline Li–Si (Li₆₄Si₆₄), obtained using the combined Matlantis atomic-level AI simulator and ReaxFF-based molecular dynamics workflow. (D) Corresponding amorphous Li–Si model (Li₆₄Si₆₄), generated through the atomistic workflow.

The energy recoverable from redox potential is given by contrasting the potential of oxidized and reduced species expressed through ΔG, the Gibbs energy (eqn (11)).⁴¹

For macroscopic electrodes, the chemical potential is determined by the probability of states via the Boltzmann distribution, resulting in the activity-dependent Nernst equation.

The challenge of eqn (12) is to evaluate the activity coefficients. In eqn (13), the Boltzmann distribution describes the concentration ratio of the investigated redox couple species. The prevailing concentration of species can be fundamentally described by the Fermi level (eqn (14)) in combination with the Fermi–Dirac distribution (eqn (15)). For Si electrodes, the two characteristic energy levels may represent the energy difference between a bulk Si phase and the lithium ions in the electrolyte, but also the difference between two distinct Si phases with differing energies.

The associated energy states depend on the relative concentrations of oxidized and reduced species. Their distribution can be expressed through the Fermi–Dirac relation, where the total concentration of a redox couple is defined as c_tot = c_ox + c_red. By combining this relation with probability density functions, the occupancy of energy states and thus the equilibrium behavior of a given redox couple can be quantified.

The validity of this description (eqn (10)–(15)) is restricted to systems that remain ideally homogeneous during charge transfer. Under such conditions, the redox couple is consistently defined by a single pair of oxidized and reduced species throughout (de)lithiation. For Si electrodes undergoing (de)lithiation, the prerequisite of an invariable redox couple is not fulfilled, as the Li–Si phase diagram (Fig. 1A) shows the sequential formation of multiple intermetallic phases, each with its own characteristic equilibrium potential. To capture this behavior, the electrode reaction must be expressed as a cascade of redox transformations (eqn (16)–(20)):

Ox^I + e⁻ ⇋ red^I (16)

Red^I ≙ ox^II (17)

Ox^II + e⁻ ⇋ red^II (18)

Red^II ≙ ox^x (19)

Ox^x + e⁻ ⇋ ox^x+1 (20)

In this cascade, each product serves as the precursor for the next transformation step, reflecting the sequential (de)lithiation of Si. This principle schematically corresponds to the sharp vertical phase boundaries of the Li–Si phase diagram, which distinguish it from the continuous solid-solution behavior observed in materials such as lithium–iron–phosphate or spinel oxides. Equilibrium potentials of the electrochemical species linked to specific material phases can be obtained from atomistic calculations, for example, density functional theory (DFT) or molecular dynamics (MD), by evaluating their relative energies with respect to a lithium reference. In this way, atomistic simulations define the reference states required for electrochemical modeling and establish a direct link between microscopic energetics and macroscopic electrode behavior.

2.4. Atomistic parametrization with phase standard potentials

To establish a consistent data basis for the MSMR framework and to assess the validity of literature-derived equilibrium potentials, atomistic simulations were conducted for both crystalline and amorphous Li–Si phases. These simulations provide quantitative reference states that bridge atomistic energetics with macroscopic electrode behavior and form the foundation for the phase-specific equilibrium potentials used in the present model. The resulting potentials are summarized and contrasted with literature data in Fig. 1B. Thin lines represent OCPs compiled from previous experimental and theoretical studies. In contrast, dots correspond to the OCPs obtained from the atomistic calculations performed in this work. The thick lines indicate the OCP data implemented in the constant-current electrochemical model. An excerpt of the atomistic configurations used to calculate these potentials is displayed in Fig. 1C (crystalline Li₆₄Si₆₄) and Fig. 1D (amorphous Li₆₄Si₆₄). Both structures correspond to a 50 : 50 Li–Si stoichiometry and were generated using a combined Matlantis/ReaxFF workflow. The crystalline model exhibits the ordered intermetallic arrangement characteristic of Li–Si alloys, whereas the amorphous structure reveals a disordered atomic network typical of the glassy phase. Both the crystalline and amorphous Li–Si phases exhibit similar equilibrium potentials, which are in good agreement with previously reported experimental and theoretical data. The associated formation energies corresponding to the illustrated structures and their OCPs are provided in SI, Fig. S1.

3. Electrochemical model

The model presented in this study is based on the multiphase nature of the Li–Si system, as outlined in the previous section. Lithium insertion into Si proceeds through a sequence of distinct intermetallic phases with confined miscibility, and each phase exhibits specific electrochemical properties. To capture this behavior, every Li–Si phase is treated as an independent redox couple with its own equilibrium potential and kinetic parameters.

The overall electrode reaction can be expressed as a cascade of transformations (eqn (16)–(20)) and is summarized by:

SiLi^x + Li + e⁻ ⇋ SiLi^x+1 (21)

This formulation does not imply that each lithium-ion insertion triggers a phase change, but rather that transformations occur once a critical SoL is reached. The equilibrium potential of each phase is given by a modified Nernst relation including a non-ideality parameter w:

where c_max,SiLi^x denotes the maximum lithium concentration in the phase SiLi^x, and c_SiLi^x the corresponding current lithium content.

Each SiLi^x phase is represented as a spherical volume, whose surface area contributes to the local exchange current density (Fig. 2B). The governing transport equations of the DFN framework (eqn (1)–(9)) are spatially resolved in the solid and electrolyte domains (Fig. 2C). Phase-specific charge transfer is described by the Butler–Volmer relation (eqn (2)), with a concentration-independent but reaction-path-dependent exchange current density j₀. The non-ideality parameter w is likewise implemented as a factor that switches between lithiation and delithiation.

Fig. 2 (A) Schematic electrochemical model, capturing the multiphase nature of the Li–Si system with distinct phases coupled via the electrolyte. (B) Representation of individual SiLi^x phases as spherical volumes, where the particle surface determines the local exchange current density. (C) Integration of the multiphase description into the DFN model, with spatial resolution of charge and mass transport in the solid and electrolyte domains.

To capture transient phase evolution, the average phase-specific SoL of the electrode is used to calculate the fraction of each phase f_SiLi^x:³

Here, c_avg.,SiLi^x denotes the average lithium concentration of phase x and c_max,SiLi^x the corresponding maximum concentration. The latter is determined from the atomic fraction of lithium in each phase relative to the most lithiated reference phase as described in the SI.

The theoretical total capacity of the electrode is given by the maximum at% Li/Si ratio (SiLi^max) in combination with the active material loading (m_Si), molar mass (M_Si) and electrode surface (A) (eqn (26))

The theoretical maximum capacity is subsequently partitioned according to the stoichiometric limits of each individual phase.

Comparison of the calculated Q_max with experimental values provides a measure of the extent to which the active Si is utilized for lithium storage within the investigated voltage and stoichiometry window (Table 1). Further details on how the capacity partitioning and concentration values are assigned to the individual phases are provided in the SI.

Table 1 Electrochemical model parameter set for constant-current (de)lithiation of a 100 nm-Si-electrode

Parameter	Model	100 nm Si-Coin Cell	Comment/ref.
C/40 current	±0.595 A m^−2	±0.616 A m^−2	Measured
Charge transfer coefficient (α)	0.5		Ref. 3
Diffusion constant electrolyte (Dl)	2.2–4.2 × 10^−12 m^2 s^−1		Ref. 52
Diffusion constant Si (Ds)	1 × 10^−12 m^2 s^−1		Ref. 52
Electrode conductivity (σs)	20 S m^−1		Assumed
Electrode length	41.3 µm	41.3 µm	Measured
Electrode surface (A)		1.13 cm^2	Measured
Electrolyte conductivity (σl)	0–0.95 S m^−1		Ref. 52
Electrolyte Li^+ concentration (cl,ref)	1000 mol m^−3	1000 mol m^−3	Assumed
Electrolyte transport number (t+)	0.12–0.37		Ref. 52
Exchange current density SiLi^0 (j^00)	Lith. = 0.19; delith. = 0.425 A m^−2		Fitted
Exchange current density SiLi^I (j^I0)	Lith. = 0.875; delith. = 0.362 A m^−2		Fitted
Exchange current density SiLi^II (j^II0)	Lith. = 1; delith. = 7.5 A m^−2		Fitted
Li/Si ratio SiLi^0	0–1.32 (0–57 at% Li)		Fig. 1
Li/Si ratio SiLi^I	1.32–3.75 (57–78 at% Li)		Fig. 1
Li/Si ratio SiLi^II	3.75–3.7525 (78–79 at% Li)		Fig. 1
Mass loading Si (mSi)	1.29 mg cm^−2	1.29 mg cm^−2	Measured
Max. capacity (loading)	23 809 mAh m^−2		Measured
Max. capacity SiLi^0	16 300 mAh m^−2		Calculated
Max. capacity SiLi^I	30 010 mAh m^−2		Calculated
Max. capacity SiLi^II	30.87 mAh m^−2		Calculated
Max. Li concentration SiLi^0	97 791 mol m^−3		Calculated
Max. Li concentration SiLi^I	180 020 mol m^−3		Calculated
Max. Li concentration SiLi^II	185.21 mol m^−3		Calculated
Max. Li concentration of Si	278 000 mol m^−3		Calculated
Max. Li/Si ratio SiLi^max	3.7525 (79 at% Li)		Calculated
Molar mass Si (MSi)	28 g mol^−1		Ref. 60
Phase transition parameter (w^0)	Lith. = 1.24; delith. = 2.10		Fitted
Phase transition parameter (w^I)	Lith. = 1.25; delith. = 4.7		Fitted
Phase transition parameter (w^II)	Lith. = 1; delith. = 4.7		Fitted
Porosity (εl)	84%	85%	Measured
Separator length/thickness	420 µm	420 µm	Measured
Si capacity utilization (Siutil)	0.51		Calculated
Solid content electrode (εs)	16%	15%	Calculated
Sphere radius SiLi^x	50 nm	50 nm	Assumed
Standard potentials	375 mV		Fig. 1
Standard potentials	180 mV		Fig. 1
Standard potentials	53 mV		Fig. 1
Temperature (T)	25 °C	25 °C	Measured
Vol% SiLi^x (electrode) (εs,eff.)	10%	8.1%	Calculated
Vol% binder (electrode)	3.6%	4.0%	Calculated
Vol% con. add. (electrode)	2.4%	2.9%	Calculated

---

4. Physicochemical model implications

As shown in Fig. 2A, the current model captures the multiphase nature of the Si–Li system by modeling charge and mass transport, connecting different Si phases through the electrolyte. The intraparticle mass and charge transport indicated by the arrow inside the particle is currently neglected. Furthermore, a constant volume is assumed for each Si phase regardless of the degree of lithiation.

4.1. Polarization-dependent reference exchange current density

The polarization-dependent reference exchange current densities (j_{0,lith/delith}) imply an asymmetry in the resulting exchange current density (j), voltage relation, as described by the Butler-Volmer equation. Fig. 3A demonstrates that the asymmetry in exchange current and potential, resulting from a polarization-dependent j_{0,lith/delith} can also be emulated through dissimilar values of the charge transfer coefficients α_lith and α_delith. From a physicochemical interpretation, this asymmetry reflects different activation barriers for oxidation and reduction, suggesting that the reaction proceeds via distinct transition states or energy pathways. Such differences may arise from variations in the local structure of the SiLi^x system, solvation environment, or structural reorganization energies along the respective reaction paths. Consequently, the kinetics of lithiation and delithiation are no longer mirror-symmetric, which has implications for overpotential behavior, reaction reversibility, and overall electrode performance. Similar models, such as those proposed by Verbrugge et al., have previously implemented this asymmetric reaction pathway in their models by introducing polarization-dependent variations in the reference exchange current densities.⁴²

Fig. 3 (A) Illustration of exchange current density and potential relations induced by polarization-dependent reference exchange current densities or charge transfer coefficients. (B) Conode representation of phase equilibria, with phase fractions determined via the lever rule. (C) Predicted concentration profiles of oxidized and reduced species derived from the Fermi–Dirac distribution, considering the effect of the non-ideality parameter w. (D) Resulting influence of w on the voltage profile, demonstrating its role in capturing deviations from ideal Nernstian behavior during (de)lithiation.

4.2. Phase-specific equilibrium potentials

The local phase-specific equilibrium potential, determined by the Nernst equation, reflects the interplay of the oxidized and reduced species, acting as reactants and products during (de)lithiation. Following the conode representation in Fig. 3B, the phase fraction of the existing phases can be determined by the lever rule. At point ❶, only pure SiLi⁰ exists, whereas the concentration of SiLi^I is absent. The product-to-reactant ratio would correspond to . In contrast, at point ❸, only products exist, but the reaction has already consumed all reactants . Thereby, the fully delithiated state (❶) would correspond to the upper energy level of the Fermi–Dirac distribution of this redox couple, whereas the fully lithiated state (❸) would correspond to the lower energy level. In this analogy, the Fermi energy would correspond to (❷). In Fig. 3C, the predicted concentration profiles for the species of the redox couple are plotted over a potential sweep according to eqn (15) (one reaction). The equal population of oxidized and reduced species refers in this analogy to the reference or Fermi potential ❷. Because the species concentration equation (eqn (15)) and the Nernst-like equation (eqn (22)) correspond to each other in their transformations, the parameter w can be considered in eqn (15) as well. The resulting influence of the w parameter on the concentration profiles of the reduced and oxidized species can also be seen in Fig. 3C.

In contrast, Fig. 3D illustrates the effect of w on the voltage profile during (de)lithiation. From a physicochemical viewpoint, the first assumption would be that parameter w represents a fudge factor influencing the Fermi–Dirac distribution and its relation of species concentration to equilibrium potential. However, modifying a fundamental distribution function in such a way would be a questionable interpretation. One straightforward but incomplete explanation could be that w represents canonical non-ideality, accounting for the deviation between concentrations and activities in the system (eqn (12)). A more physically grounded interpretation of the parameter w can be obtained by relating it to the regular-solution interaction parameter Ω. In the regular-solution theorie, non-ideality enters the electrochemical potential as an additive enthalpic term , whereas the simplified expression used in this work incorporates non-ideality through a multiplicative scaling of the entropic logarithmic term.^43,44 Further information about the relation between regular solution theory and the MSMR framwork are provided in the SI.

5. Experimental and numerical framework

5.1. Silicon electrode preparation

The Si electrodes were prepared aiming for a target solid content of 60 mass% active material (nano-Si, d₅₀ ≈ 10 nm, d₅₀ ≈ 40 nm, and d₅₀ ≈ 100 nm, Merck), 20 mass% conductive additives (KS6L and C65 in a 2.25 : 1 volume ratio), and 20 mass% binder. The masses of KS6L and C65 were adjusted according to their densities (2.26 g cm⁻³ and 1.8 g cm⁻³, respectively) to ensure the specific volume fraction.

Initially, 445 mg of KS6L and 160 mg of C65 were mixed into 15.3 g of a 4 mass% LiPAA solution (prepared by neutralizing polyacrylic acid, M_v ≈ 459 kDa, Sigma-Aldrich, with 6 M LiOH to pH 7) using a magnetic stirrer (VELP Multistirrer) at 1000 rpm overnight under a closed lid to minimize water evaporation. Subsequently, 1.95 g of nano-Si was added and mixed for 8 h. Before coating, the mixture was homogenized in a SpeedMixer (DAC 150 SP, Hauschild) for 5 min at 2000 rpm. Coating was conducted via the doctor blade method (MTI Corporation, MSK-AFA-HC100) on a 25 µm thick copper foil (MTI) with a wet coating thickness of 80 µm and a feed rate of 15 mm s⁻¹. Drying of the electrodes was conducted in a vacuum oven at 80 °C and 60 mbar overnight.

5.2. Cell assembly and electrochemical characterization

Circular electrode discs with a diameter of 12 mm were punched from the uncalendared Si electrode using a precision punch tool (EL-Cut, EL-Cell GmbH). The mass of each electrode was recorded prior to cell assembly, corresponding to an average active material mass-loading of 0.73 mg cm⁻² for the d₅₀ = 10 nm Si material, 1.21 mg cm⁻² for the d₅₀ = 40 nm and 1.29 mg cm⁻² for the d₅₀ = 100 nm Si. Glass microfiber filters (GF/F, Whatman) were used as separators and punched into 18 mm discs (Seal Punch, Hoffmann Group). The separators were dried overnight under vacuum at 60 °C and transferred into an argon-filled glovebox (MBraun; O₂, H₂O <0.1 ppm) for cell assembly. Coin cells (CR2032 type) were assembled in a two-electrode configuration, using lithium metal (Sigma-Aldrich) as the counter and reference electrode, and the Si electrode as the working electrode. The electrolyte consisted of 1 M LiPF₆ in a 1 : 1 (v/v) mixture of ethylene carbonate (EC) and dimethyl carbonate (DMC) (Sigma Aldrich). For assembly, a 0.5 mm stainless-steel spacer was first placed into the positive case to ensure electrical contact with the working electrode, which was positioned directly above the spacer. The GF/F separator was then placed on top of the working electrode and soaked with 150 µL of electrolyte. A lithium metal chip with a 15.6 mm diameter (MSE PRO) was positioned centrally on the separator, followed by a second spacer and a wave spring to ensure mechanical stability and reliable electrical contact. The stack was enclosed with a negative case and hermetically sealed using a coin cell crimping device (HSMCC-H10, Wellcos Corporation). After a minimum wetting time of 12 h, the cells underwent formation and subsequent cycling in an Arbin cycler unit at 25 ± 1 °C. The C/40 constant-current (de)lithiation was performed following three preliminary charge–discharge adjustments, which served as formation cycles, to ensure a discharge duration of approximately 40 h within the voltage window of 1000–100 mV. For the pulse-relaxation protocol, the cells (de)lithiated five times for 8 h each at a C/40 current, with each step followed by a 36 h relaxation period.

5.3. Atomistic simulation

Atomistic simulations were performed using the ReaxFF reactive forcefield method as implemented in the Amsterdam Modeling Suite (AMS) simulation software and the Matlantis AI simulator, a web-based atomic-level simulation tool, to evaluate the electrochemical properties of crystalline and amorphous Li–Si phases.^45–47 Structural prototypes for the crystalline Li–Si compounds were obtained from the Materials Project database.⁴⁸ Geometry optimization of these bulk structures was carried out in Matlantis using the PFP v7.0.0 model in the CRYSTAL U0 mode estimator. This structural relaxation was performed using the FIRE optimizer until all atomic forces converged below 0.01 eV Å⁻¹, yielding minimum-energy configurations.⁴⁹ The OCP of each phase was then calculated directly from the total energies of these optimized bulk structures. Preliminary amorphous structures were generated using molecular-dynamics simulations based on the ReaxFF reactive force field as implemented in the AMS software. To obtain amorphous Li–Si configurations, crystalline Li–Si systems were heated from 300 K to 3000 K over 12.5 ps, equilibrated at 3000 K for 12.5 ps, and subsequently cooled back to 300 K over 12.5 ps at a constant rate. The resulting amorphous structures were then further relaxed in Matlantis, which was also used to compute their energies. All molecular-dynamics steps were conducted using the ReaxFF reactive force-field methodology implemented in the AMS, with the Li–Si ReaxFF forcefield parameter set developed by Ostadhossein et al.⁵⁰ Additional details regarding parameter selection, simulation time steps, and potential fitting are provided in the SI (Tables S1–S4).

5.4. Electrochemical simulation

The electrochemical simulations were carried out using COMSOL Multiphysics 6.3®, specifically employing the Battery Design Module to resolve the coupled electrochemical and transport processes.^51,52 Model parameterization was based on the experimental and atomistic data described in sections 2 and 6. To enable consistent comparison between measured and simulated results, the MATLAB LiveLink interface was applied for data handling, post-processing, and automated evaluation of the simulation output.

6. Results and discussion

6.1. Si-electrodes and microstructures

To ensure comparability between electrodes of different Si particle sizes, a single slurry formulation was selected that enabled reliable electrode fabrication without requiring any changes in composition. Fig. 4 visualized through scanning electron micrographs the Si active material as powders (Fig. 4A, C and E) and in the manufactured electrodes (Fig. 4B, D and F). For both the plain powders and the electrode, the increase in particle size (left to right) is apparent. Clearly differentiable from the 100 nm Si active material is in Fig. 4F the conductive additive C65 due to its significantly smaller particle size of approximately 20 nm.⁵³ Since the particle size difference between C65 and the Si active material is not very large in the 10 nm and 40 nm electrodes, traces of C65 cannot be clearly distinguished visually. However, in the 10 nm electrodes, the KS6L conductive graphite can be identified by its particle size difference from the Si nanoparticles, as shown in the SI, Fig. S2.

Fig. 4 Scanning electron micrographs of Si powders and corresponding electrodes. (A, C and E) Pristine Si powders with average particle sizes of 10 nm, 40 nm, and 100 nm. (B, D and F) Cross-sections of electrodes fabricated from the respective powders. The increase in particle size is evident from left to right. In the 100 nm electrode (F), the smaller carbon black particles (∼20 nm) can be distinguished from the Si active material.

6.2. Electrochemical properties of Si-half cells

In Fig. 5A, the first two cycles for the electrode using Si-nanoparticles as active material are shown. Evident is the distinct difference in the OCP of Si in the first cycle, apparent through the plateau-like voltage attributed to the initial amorphization of the Si active material.^54,55 The SoC is evaluated by integrating charge over time, whereby the SoC for charge and discharge was determined independently to compensate for the Coulombic efficiency, which was on average 74 ± 12% for the first cycle and 92 ± 12% for the second across all electrodes. Omitting the loss of lithium-ion host capacity is necessary for the SoC estimation in this study since charge loss mechanisms have not been considered in the electrochemical model. In Fig. 5B, the capacity of a C/40 (de)lithiation is normalized to the Si content in the electrode for the 10 nm, 40 nm, and 100 nm Si material. Using the specific capacity of graphite as a reference, the superior specific capacity of Si as an active material in lithium-ion batteries becomes evident. Among the investigated particle sizes, the 40 nm material displayed the highest delithiation capacity of 2179 mAh g⁻¹, followed by the 100 nm active material with 2021 mAh g⁻¹ and the 10 nm Si electrode with 1070 mAh g⁻¹. The delithiation capacities followed the same trend as the lithiation capacities.

Fig. 5 Electrochemical characterization of Si half-cells with 10 nm, 40 nm, and 100 nm particle sizes. (A) First two cycles aiming for C/40, showing the plateau-like voltage during the first lithiation associated with the amorphization of crystalline Si. (B) Normalized capacities for lithiation and delithiation, demonstrating the superior specific capacity of Si compared to graphite (ref. 59), with highest values for the 40 nm material. (C) Charge and discharge voltage profiles within 1000–100 mV, with 100% SoC defined at 100 mV cutoff with accurate C/40 current. (D) Differential capacity plots (dQ/dV), with main (de)lithiation features (❶–❻) assigned to distinct phase transformations and contributions from carbon additives.

Even higher lithiation capacities could potentially be achieved if lithiation were carried out below 100 mV vs. Li/Li⁺ reference. However, further lithiation was intentionally avoided, since the high SoL and the associated phase transformation coupled with volumetric expansion would impose severe mechanical stress on the material, which in turn would compromise the cycling stability. In Fig. 5C, the charge and discharge profiles of Si half-cells are shown for a voltage window between 1000 and 100 mV, where the capacity at the cutoff voltage of 100 mV is defined as 100% SoC. The charge and discharge times for all cells were 40.98 ± 1.1 h, indicating that the applied currents were experimentally adjusted to enforce comparable cycling durations and overpotentials. Applied C/40 currents normalized to the active material mass were 26.3 µA mg⁻¹ for the 10 nm electrode, 51.3 µA mg⁻¹ for the 40 nm electrode, and 47.6 µA mg⁻¹ for the 100 nm electrode. In Fig. 5D, the corresponding differential capacity plots of the C/40 cycle are shown. The mass loading of the associated electrodes was 0.80 mg for 10 nm, 1.45 mg for 40 nm, and 1.46 mg for 100 nm. Based on the 10 nm electrode (0.80 mg, 21.2 µA) as reference, the expected C/40 currents would scale to 38.6 µA for 40 nm and 38.6 µA for 100 nm.

In contrast, the applied currents were substantially higher: 74.1 µA for 40 nm (≈95% above expectation) and 69.7 µA for 100 nm (≈81% above expectation). Despite the identical electrode formulation, the disproportionate relation between electrode mass loading and applied current across all three particle sizes indicates a particle-size-dependent reference exchange current density. This behavior can be linked to the reaction rate constant and thus to a particle-size-dependent activation energy that should be further explored. In addition, structural factors such as particle agglomeration or electrode morphology may further influence the effective electrochemically active surface area, thereby contributing to the observed deviations.

Across all particle sizes, the differential capacity analysis (Fig. 5D) revealed that the main lithiation and delithiation peaks occurred within nearly identical voltage ranges, indicating that the fundamental reaction pathway is not significantly affected by particle size. However, the peak intensity increased systematically with particle size: the 10 nm electrode displayed the weakest peaks, the 40 nm electrode showed an intermediate intensity, and the 100 nm electrode exhibited the most pronounced peaks. The relative peak intensities also reflect the capacity contribution of the corresponding redox reactions; hence, the smaller peaks in the 10 nm electrode directly explain its lower overall capacity compared to the 40 nm and 100 nm electrodes. This trend may further be associated with a particle-size-dependent activation energy of the underlying redox reactions.

The individual peaks can be assigned to distinct electrochemical processes (see symbols ❶–❻ in Fig. 5D). The first lithiation peak (❶) corresponds to the transformation of amorphous Si (SiLi⁰) into the SiLi^I which is often denoted as Li₂Si phase or Li–Si phase in the literature.^55,56 The second lithiation feature (❷) indicates the completion of the SiLi^I formation process and the onset of the SiLi^II phase transformation. The third lithiation feature (❸) marks the further transformation of the SiLi^I phase into the SiLi^II phase, which corresponds to Li_3.5Si according to Graf et al.⁵⁵ Upon delithiation, the first peak observed (❹) is assigned to the deintercalation of the carbon additives.⁵⁵ The subsequent fifth peak (❺) represents the back-transformation of SiLi^II into SiLi^I, and finally, the sixth peak (❻) corresponds to the lithium delithiation of SiLi^I into amorphous SiLi⁰.

The first lithiation peak, corresponding to the SiLi^I phase transformation has been discussed in the literature with different proposed stoichiometries. While some studies suggest a clear composition, such as Li₂Si,²² a broader survey (Fig. 1B) places the phase and associated potentials of the SiLi^I between a stoichiometric range approximately between 50 at% and 75 at% Li. In the present model, the exact stoichiometry of SiLi^I is less critical than its relation to the maximum Li/Si ratio attainable, as this determines the fraction of the total capacity attributed to this phase. In general, the stoichiometries of the SiLi^x phases must be interpreted in relation to the highest achievable Li/Si ratio (at% Li in SiLi^max). Several values have been reported: Graf et al. suggest 3.5 for cutoff voltages below 170 mV, while Jiang and Chan et al. propose 3.75 without considering a possible super-alloyed Li–Si phase.^3,55,57 Braga et al. predict the existence of highly lithiated phases, such as Li₂₁Si₅, which would correspond to a maximum Li/Si ratio of 4.2.²⁵ Moreover, first-principles simulations often focus on a single type of morphology (e.g., crystalline or amorphous), which may limit their applicability to real electrode structures. Taken together, the consolidated literature indicates a feasible range of 3.00–4.20 for the maximum Li/Si ratio, in agreement with the modeling assumptions employed in this work.

The key outcome of the differential capacity analysis is that the cutoff potential, and thereby the degree of lithiation, directly governs the existence of different Li–Si phases. At lower anode potentials, phases with higher Li content become thermodynamically accessible. Lithiation must therefore be regarded as a cascade of reactions, in which the formation of one phase provides the reactant state for the subsequent transformation. This cascade principle is also reflected in the modeling approach of this study.

An important contribution to this line of interpretation was made by Berg et al., who presented the electrode potential as a function of the lithiation degree.⁵⁸ In their study, lithiation down to a cutoff potential of 170 mV versus Li/Li⁺ was considered.⁵⁸ Within this voltage window, the assumption was made that no phases beyond SiLi^I are formed, such that the observed capacity originates exclusively from the amorphous Li₂Si-like phase, which reaches its maximum capacity at 170 mV versus Li/Li⁺. For the electrodes investigated in this study, the corresponding voltage at which SiLi^I phase formation is completed, as represented by feature (❷) in Fig. 5D. However, in this study, lithiation is extended to a lower cutoff potential of 100 mV versus Li/Li⁺. As a result, the formation of the SiLi^II phase must be considered, since the additional voltage range below 170 mV enables the incorporation of more lithium into the Si structure. Beyond this, it is further assumed that continued lithiation initiates the formation of the SiLi^II phase, representing the most Li-rich state accessible within the investigated potential window.³ Consequently, the lithiation process observed here cannot be described solely by the SiLi⁰ → SiLi^I phase transformation, but must instead be treated as a sequence of consecutive phase transformations (SiLi⁰ → SiLi^I → SiLi^II), each contributing a distinct capacity fraction and serving as precursors for the subsequent reaction step.

6.3. Model vs. constant-current (de)lithiation of a 100 nm-Si-electrode

To demonstrate the applicability of the model presented in section 3 for simulating constant-current (de)lithiation processes, the parameterization was performed in a physics-based manner and tailored to accurately reproduce the electrochemical behavior of the fabricated 100 nm Si electrode. A constant-current C/40 (de)lithiation cycle with 98% Coulombic efficiency from the 100 nm-Si-electrode served as the ground truth. Table 1 summarizes the experimentally determined values together with the assumed and fitted parameters implemented in the simulation. In Fig. 1B, the standard potentials assigned to the respective material phases are indicated by thick lines, highlighting their proximity to standard potentials predicted by atomistic/first-principles simulation, literature and experiments in their respective atomic ratios. It is evident that the standard potentials of SiLi⁰ and SiLi^II lie within the standard deviation of the values obtained from the data consolidation of the phase-specific equilibrium potential. In contrast, the standard potential of SiLi^I deviates slightly from these predictions.

Fig. 6A compares the experimentally measured voltage–SoC data with the simulation results, showing a reasonable agreement between model and experiment. The absolute voltage error over SoC is shown in Fig. 6B for lithiation and delithiation. The major deviations occur at the low- and high-SoC regions (<2% and >98%) with an overall root mean square error (RMSE) of 29.0 mV for lithiation and 11.2 mV for delithiation, corresponding to 3.2% and 1.2% of the 900 mV voltage window, respectively. In contrast, when restricting the analysis to the 2–98% SoC range, the RMSE is reduced to only 5.4 mV for lithiation and 7.6 mV for delithiation, corresponding to 0.6% and 0.8%. The discrepancy at the SoC boundaries originates from the fact that a different phase transition parameter (w) would be required in these regions to meet the experimental data with more accuracy, reflecting lithiation-level-dependent reaction kinetics. Improved model fidelity could therefore be achieved by incorporating lithiation-dependent phase transformations, particle size distributions, and particle-size-dependent phase transition parameters. Furthermore, instead of rigid volumes assigned to each material phase, an enhanced model would simulate the dynamic growth and shrinkage of phases.

Fig. 6 Model validation against constant-current cycling of a 100 nm Si electrode at C/40. (A) Comparison of experimental and simulated voltage–SoC curves, showing good agreement. (B) Absolute voltage error between experiment and simulation, with RMSE values of 29.0 mV (lithiation) and 11.2 mV (delithiation); major deviations occur at SoC edges (<2% and >98%). (C) Evolution of phase fractions during cycling, highlighting the sequential transformation cascade and the asymmetric character of lithiation vs. delithiation. (D) State of lithiation (SoL) of each phase as a function of time, illustrating preferential lithiation of the phase with the highest equilibrium potential and the stepwise cascade of phase transitions.

Further insights into the phase fraction evolution during electrochemical alloying can be drawn from Fig. 6C, where the 40 h mark indicates current reversal and the onset of delithiation. During lithiation, the fraction of pristine SiLi⁰ progressively decreases and vanishes shortly after 20 h. At this point, the lithiation of the SiLi^I phase reaches its maximum, marking the complete conversion of SiLi⁰ into SiLi^I, while not-significant amounts of SiLi^II phase are present. Subsequent lithiation then initiates the transformation of SiLi^I into SiLi^II, as reflected by the decline of SiLi^I and the concurrent increase of SiLi^II phase fraction. During delithiation, the reaction cascade is reversed, with SiLi^II first converting into SiLi^I followed by the reappearance of delithiated and amorphous SiLi⁰. However, the lithiation and delithiation pathways are not fully symmetric. This asymmetry is evident from the distinct SiLi^I fraction peak around 55 h, which displays a lower intensity during delithiation compared to lithiation. In addition, the re-emergence of SiLi⁰ phase fraction is apparent as a non-mirror-symmetric graph to a vertical line drawn at 40 h, underscoring the asymmetric character of the underlying processes. From a modeling perspective, this asymmetric behavior is triggered by the use of a distinct phase transition (w^x) parameters and exchange current densities (j₀^x) for lithiation and delithiation (i.e.: polarization-dependent), which results in the experimentally observed voltage hysteresis.

Fig. 6D corroborates these observations by depicting the SoL of each SiLi^x phase. A steep ascent or descent in SoL indicates that a relatively large amount of Li is accommodated in, or released from, the respective phase volume for a time interval. Schematically, the ongoing lithiation and delithiation mechanism is presented in Fig. 8. According to the Nernst equation, each phase lowers its equilibrium potential upon lithiation. Thus, lithiation proceeds preferentially in the phase with the highest equilibrium potential versus Li/Li⁺, initially SiLi⁰ (Fig. 8A). As lithiation progresses, the potential of SiLi⁰ decreases until it is surpassed by that of the neighboring phase, leading to the preferential lithiation of SiLi^I (Fig. 8B). This is evident in Fig. 6D as a steep increase in the lithiation level of SiLi^I. The potential-driven thermodynamic driving force previously described, likewise, governs the lithiation of the SiLi^II phases and the corresponding phase transformations (Fig. 8C). During delithiation, the analogy holds with reversed polarization: the phase with the lowest equilibrium potential delithiates first, and its potential increases in Nernstian fashion until another phase becomes thermodynamically favored for delithiation (Fig. 8E–G).

6.4. Model vs. C/40 pulse-relaxation protocol with 8 h load and 36 h relaxation of a 100 nm-Si-electrode

Fig. 7A presents the experimentally determined voltage-SoC curves of the pulse-relaxation experiment for the 100 nm electrode, together with the corresponding simulation results. The model parameters in Fig. 6 and Table 1 have been adjusted for the simulation of the pulse-relaxation protocol with the experimental data from the corresponding electrode and cell, as shown in SI, Table S5. During relaxation, the voltage clearly tends towards the midpoint between the enveloping lithiation and delithiation OCP curves. The vertical appearance of the relaxation segments arises from the fact that the cell SoC remains constant while the voltage evolves during relaxation. Fig. 7B displays the corresponding voltage deviations between experiment and simulation over the SoC dimension. For the SoC values at which relaxation occurred (20%, 40%, 60%, 80%, and 100%), the last recorded data points of both experiment and simulation were compared. The resulting RMSE ranged from 36.9 mV for lithiation to 33.2 mV for delithiation, corresponding to a relative RMSE of 4.1% and 3.6% of the 900 mV voltage window. As already observed for constant-current cycling, the largest discrepancies between experiment and simulation occurred in the SoC edge zones (<2% and >98%).

Fig. 7 Comparison of experiment and simulation for the C/40 pulse-relaxation protocol of a 100 nm Si electrode (8 h load, 36 h relaxation). (A) Voltage–SoC curves showing convergence of relaxation segments toward the midpoint between lithiation and delithiation OCPs. (B) Voltage deviations between experiment and simulation at the onset of relaxation, with RMSE values of 36.9 mV for lithiation and 33.2 mV for delithiation. (C) Transient voltage response from experiment and simulation. (D) Absolute error over 400 h, with an overall RMSE of 27.0 mV. Phase fraction (E) and SoL (F) evolution during relaxation, demonstrating ongoing lithium redistribution between phases toward thermodynamically favored equilibrium states.

In Fig. 7C, the transient voltage profiles from experiment and simulation are compared directly, while Fig. 7D illustrates the corresponding error over a total time of 400 h. In Fig. 7D the RMSE of experiment and simulation corresponded to 27.0 mV over 400 h (relative RMSE of 3.0%). Finally, Fig. 7E and F provide insights into the processes occurring during the respective 36 h relaxation periods. The model predicts that the phase fractions continue to evolve or decay during relaxation, accompanied by the insertion and extraction of lithium ions. The interaction between the individual material phases becomes evident during the relaxation periods following lithiation (8–44 h, 52–88 h, 96–132 h, 140–176 h, and 184–220 h) as well as during those following delithiation (228–264 h, 272–308 h, 316–352 h, and 360–396 h). Throughout these intervals, mass transport of lithium-ions proceeds toward the thermodynamically favored phases, i.e., those with higher equilibrium potential, until the assimilation of the respective phase equilibrium potential is achieved.

This ongoing mass transport during relaxation is schematically represented in Fig. 8D for lithiation and Fig. 8H for delithiation. In particular, the SiLi^II phase releases lithium into the SiLi⁰ and SiLi^I phases, as their associated reference equilibrium potentials are higher. For the current parametrization, the model predicts that up to a SoC of 40% (corresponding to 52 h), the existence of the SiLi^II phase is thermodynamically so unfavored that it vanishes completely during relaxation by transferring its lithium ions to the other phases. During the relaxation period after 60% SoC (yellow segment, left border at 96 h), the phase fraction of SiLi^II continues to decrease but stabilizes once inter-phase transport processes are completed and phase potential differences are equalized.

Fig. 8 Schematic illustration of the phase cascade mechanism during cycling and relaxation of Si electrodes. (A–C) Lithiation proceeds stepwise, with lithium insertion occurring preferentially in the phase with the highest equilibrium potential until surpassed by that of the next phase. (D) During relaxation after lithiation, lithium redistributes from thermodynamically unfavored to favored phases. (E–G) Delithiation occurs in reverse order, beginning with the phase of lowest equilibrium potential and proceeding toward higher potentials. (H) Relaxation after delithiation drives further lithium redistribution until the phase equilibrium potentials are balanced.

At the onset of relaxation after reaching 60% SoC (96 h), the phase fraction amounts to approximately 7% SiLi⁰, 81% SiLi^I, and 12% SiLi^II. By the end of the relaxation period (132 h), these fractions shift to about 2% SiLi⁰, 90% SiLi^I, and 8% SiLi^II. Accordingly, the phase fraction of SiLi⁰ and SiLi^II decrease (black and blue arrow) during the associated relaxation period, whereas the phase fraction of SiLi^I increased (red arrow). During the increase of the SiLi^I fraction, its SoL decreases (12 → 8.3%), while the SoL of SiLi⁰ (93 → 98%) increases. This apparent paradox is resolved because a larger volume fraction is assigned to the SiLi^I phase at the end of relaxation, which lowers its effective SoL. Although the absolute amount of lithium remains nearly unchanged, it becomes distributed over an increased phase volume, thereby reducing the lithium concentration and resulting in a lower SoL, according to eqn (4).

7. Conclusions

In this work, we introduced an electrochemical modeling framework for Si anodes that explicitly accounts for the multiphase nature of the Li–Si system. By integrating atomistic and thermodynamic insights into a modified MSMR approach, the model successfully reproduces experimentally observed voltage hysteresis, SoC-dependent relaxation behavior, and phase evolution during lithiation and delithiation.

Comparison with Si half-cell data demonstrated that the model captures the dissimilar electrochemical potential behavior of the lithiation and delithiation reaction pathways, including non-symmetric phase transformations and fraction dynamics. The incorporation of phase-specific equilibrium potentials and reaction-path-dependent exchange current densities enables the model to reproduce voltage hysteresis and long relaxation times. These phenomena remain elusive for conventional DFN-type models.

Quantitative evaluation underlines the predictive capability of the framework. For the constant-current charge/discharge experiments of the 100 nm Si electrode, the RMSE across the SoC dimension amounted to 29.0 mV for lithiation and 11.2 mV for delithiation (relative RMSE 3.2% and 1.2%). For the dynamic pulse-relaxation protocol (C/40 with 8 h load and 36 h relaxation), the deviations were 36.9 mV for lithiation and 33.2 mV for delithiation, corresponding to a relative RMSE of 4.1% and 3.6% across the SoC dimension. In contrast, the error over the entire 400 h time domain was 27 mV (relative RMSE of 3.0%). These results highlight the model's ability to not only reproduce steady-state behavior but also capture transient relaxation phenomena with high accuracy.

The analysis of pulse-relaxation experiments further revealed that intra- and interphase lithium redistribution continues during rest periods, governed by thermodynamic driving forces toward phases with higher equilibrium potentials. This validates the importance of coupling atomistic-level parameters with electrode-scale simulations to bridge the gap between fundamental material behavior and macroscopic electrochemical performance. As a next step, future models could abandon the assumption of static phase-assigned volumes to enable continuous phase transformations and thus achieve an even more realistic representation of lithiation and delithiation dynamics.

Overall, the presented modeling framework provides a robust platform for studying the interplay between phase transformations, voltage hysteresis, and lithium transport in Si anodes. Beyond deepening the mechanistic understanding of the Li–Si system, this approach offers a pathway to plausibly and partially validate atomistic and first-principles predictions at the electrode level. Most importantly, the proposed methodology achieves a true cross-dimensional linkage between atomistic insights and macroscopic electrode-scale behavior, thereby supporting the design of next-generation batteries and advanced battery management strategies.

Author contributions

N. P.: methodology, investigation, data curation, validation, writing – original draft, writing – review and editing. O. Q.: methodology, validation, writing – review and editing. S. A.: validation, writing – review and editing. S. B.: investigation, writing – review and editing. J. M.: validation, writing – review and editing. S. S.: methodology, validation, writing – review and editing. V. P.: visualization, writing – review and editing.

Conflicts of interest

The authors have no conflict of interest.

Data availability

Data and models are available via https://doi.org/10.5281/zenodo.17600064.

Supplementary information (SI) is available. Content: exchange current density at equilibrium, relation between regular-solution parameter Ω and model parameter w, atomistic simulation details plus Li-Si potential/structure tables, electrochemical model assumptions and phase-capacity/concentration assignment, Si-electrode SEM micrographs with KS6L graphite, ReaxFF and Matlantis settings, pulse-relaxation model parameters, formation-energy plot for crystalline/amorphous Li-Si, supporting references. See DOI: https://doi.org/10.1039/d6eb00061d.

Acknowledgements

The authors thank Armin Asheri for the valuable scientific exchange and insightful discussions that have greatly contributed to this work. They also express their gratitude to Johannes Sperber for his support in addressing modeling challenges and to Justus Speichermann for his organizational assistance throughout the research. We acknowledge the support of the state government of Saarland within the framework of the EnFoSaar project, funded by the Saarland Transformation Program for Research and Knowledge Transfer (Transformationsprogramm Forschung und Wissenstransfer Saar).

References

1. L. Sun, Y. Liu, J. Wu, R. Shao, R. Jiang, Z. Tie and Z. Jin, Small, 2022, 18, e2102894.
2. G. Zhu, W. Luo, L. Wang, W. Jiang and J. Yang, J. Mater. Chem. A, 2019, 7, 24715–24737.
3. Y. Jiang, G. Offer, J. Jiang, M. Marinescu and H. Wang, J. Electrochem. Soc., 2020, 167, 130533.
4. M. T. McDowell, S. W. Lee, W. D. Nix and Y. Cui, Adv. Mater., 2013, 25, 4966–4985.
5. M. P. Bonkile, Y. Jiang, N. Kirkaldy, V. Sulzer, R. Timms, H. Wang, G. Offer and B. Wu, J. Power Sources, 2024, 606, 234256.
6. K. Nishikawa, J. Moon and K. Kanamura, J. Power Sources, 2016, 302, 46–52.
7. A. J. Louli, J. Li, S. Trussler, C. R. Fell and J. R. Dahn, J. Electrochem. Soc., 2017, 164, A2689–A2696.
8. Z.-L. Xu, X. Liu, Y. Luo, L. Zhou and J.-K. Kim, Prog. Mater. Sci., 2017, 90, 1–44.
9. O. Von Kessel, A. Avdyli, D. Vrankovic and K. P. Birke, J. Power Sources, 2024, 610, 234582.
10. A. Van Der Ven, K. A. See and L. Pilon, Battery Energy, 2022, 1, 20210017.
11. W. Dreyer, J. Jamnik, C. Guhlke, R. Huth, J. Moškon and M. Gaberšček, Nat. Mater., 2010, 9, 448–453.
12. D. Wycisk, M. Oldenburger, M. G. Stoye, T. Mrkonjic and A. Latz, J. Energy Storage, 2022, 52, 105016.
13. M. Z. Bazant, Acc. Chem. Res., 2013, 46, 1144–1160.
14. D. Burch, Ph.D., Massachusetts Institute of Technology, 2009.
15. R. B. Smith and M. Z. Bazant, J. Electrochem. Soc., 2017, 164, E3291–E3310.
16. L. Köbbing, A. Latz and B. Horstmann, Adv. Funct. Mater., 2023, 34, 2308818.
17. L. Köbbing, Y. Kuhn and B. Horstmann, ACS Appl. Mater. Interfaces, 2024, 16, 67609–67619.
18. M. W. Verbrugge, D. R. Baker, X. Xiao, Q. Zhang and Y.-T. Cheng, J. Phys. Chem. C, 2015, 119, 5341–5349.
19. D. R. Baker and M. W. Verbrugge, J. Electrochem. Soc., 2018, 165, A3952–A3964.
20. M. Verbrugge, D. Baker, B. Koch, X. Xiao and W. Gu, J. Electrochem. Soc., 2017, 164, E3243–E3253.
21. D. R. Baker, M. W. Verbrugge and X. Xiao, J. Appl. Phys., 2017, 122, 165102.
22. R. A. Sharma and R. N. Seefurth, J. Electrochem. Soc., 2019, 123, 1763–1768.
23. C. J. Wen and R. A. Huggins, J. Solid State Chem., 1981, 37, 271–278.
24. W. J. Weydanz, M. Wohlfahrt-Mehrens and R. A. Huggins, J. Power Sources, 1999, 81, 237–242.
25. M. H. Braga, A. Dębski and W. Gąsior, J. Alloys Compd., 2014, 616, 581–593.
26. S. B. Olou'ou Guifo, J. E. Mueller, D. Henriques and T. Markus, Phys. Chem. Chem. Phys., 2022, 24, 9432–9448.
27. S. B. Olou'ou Guifo, J. E. Mueller and T. Markus, J. Phys. Chem. C, 2024, 128, 4891–4904.
28. H. Valencia, P. Rapp, M. Graf, J. Mayer and H. A. Gasteiger, J. Electrochem. Soc., 2024, 171, 120507.
29. I. Valencia-Jaime, R. Sarmiento-Pérez, S. Botti, M. A. L. Marques, M. Amsler, S. Goedecker and A. H. Romero, J. Alloys Compd., 2016, 655, 147–154.
30. V. L. Chevrier and J. R. Dahn, J. Electrochem. Soc., 2009, 156, A454.
31. G. A. Tritsaris, K. Zhao, O. U. Okeke and E. Kaxiras, J. Phys. Chem. C, 2012, 116, 22212–22216.
32. H. D. Deng, N. Jin, P. M. Attia, K. Lim, S. D. Kang, N. Kapate, H. Zhao, Y. Li, M. Z. Bazant and W. C. Chueh, ACS Nano, 2024, 18, 2210–2218.
33. P. Bai, D. A. Cogswell and M. Z. Bazant, Nano Lett., 2011, 11, 4890–4896.
34. M. A. Roscher, O. Bohlen and J. Vetter, Int. J. Electrochem., 2011, 2011, 1–6.
35. C. P. Graells, M. S. Trimboli and G. L. Plett, 2020 IEEE Transportation Electrification Conference & Expo (ITEC), 2020.
36. J. Schmitt, I. Horstkötter and B. Bäker, J. Power Sources Adv., 2024, 26, 100139.
37. M. T. McDowell, I. Ryu, S. W. Lee, C. Wang, W. D. Nix and Y. Cui, Adv. Mater., 2012, 24, 6034–6041.
38. X. H. Liu, L. Zhong, S. Huang, S. X. Mao, T. Zhu and J. Y. Huang, ACS Nano, 2012, 6, 1522–1531.
39. K. Pan, F. Zou, M. Canova, Y. Zhu and J.-H. Kim, J. Power Sources, 2019, 413, 20–28.
40. M. Doyle, T. F. Fuller and J. Newman, J. Electrochem. Soc., 1993, 140, 1526–1533.
41. P. Peljo, C. Villevieille and H. H. Girault, Energy Environ. Sci., 2025, 18, 1658–1672.
42. M. Verbrugge, D. Baker and X. Xiao, J. Electrochem. Soc., 2015, 163, A262–A271.
43. P. Ombrini, A. Vasileiadis and M. Wagemaker, J. Electrochem. Soc., 2025, 172, 090524.
44. P. Ombrini, M. Z. Bazant, M. Wagemaker and A. Vasileiadis, npj Comput. Mater., 2023, 9, 148.
45. AMS GUI 2025.1, SCM, Amsterdam, The Netherlands, https://www.scm.com.
46. Matlantis® (PFP-based atomic simulation platform). https://www.matlantis.com. Preferred Networks, Inc., Tokyo, Japan.
47. S. Takamoto, C. Shinagawa, D. Motoki, K. Nakago, W. Li, I. Kurata, T. Watanabe, Y. Yayama, H. Iriguchi, Y. Asano, T. Onodera, T. Ishii, T. Kudo, H. Ono, R. Sawada, R. Ishitani, M. Ong, T. Yamaguchi, T. Kataoka, A. Hayashi, N. Charoenphakdee and T. Ibuka, Nat. Commun., 2022, 13, 2991.
48. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, APL Mater., 2013, 1, 011002.
49. E. Bitzek, P. Koskinen, F. Gahler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
50. A. Ostadhossein, E. D. Cubuk, G. A. Tritsaris, E. Kaxiras, S. Zhang and A. C. van Duin, Phys. Chem. Chem. Phys., 2015, 17, 3832–3840.
51. COMSOL Multiphysics® v. 6.3. https://www.comsol.com. COMSOL AB, Stockholm, Sweden.
52. Battery Design Module®. COMSOL Multiphysics® v. 6.3. https://www.comsol.com. COMSOL AB, Stockholm, Sweden.
53. K. Pfeifer, S. Arnold, Ö. Budak, X. Luo, V. Presser, H. Ehrenberg and S. Dsoke, J. Mater. Chem. A, 2020, 8, 6092–6104.
54. W.-J. Zhang, J. Power Sources, 2011, 196, 877–885.
55. M. Graf, C. Berg, R. Bernhard, S. Haufe, J. Pfeiffer and H. A. Gasteiger, J. Electrochem. Soc., 2022, 169, 020536.
56. E. Feyzi, M. R. Anil Kumar, X. Li, S. Deng, J. Nanda and K. Zaghib, Next Energy, 2024, 5, 100176.
57. M. K. Chan, C. Wolverton and J. P. Greeley, J. Am. Chem. Soc., 2012, 134, 14362–14374.
58. C. M. Berg, R. Morasch and H. A. Gasteiger, J. Electrochem. Soc., 2025, 172, 050516.
59. M. Naguib, M. Kurtoglu, V. Presser, J. Lu, J. Niu, M. Heon, L. Hultman, Y. Gogotsi and M. W. Barsoum, Adv. Mater., 2011, 23, 4248–4253.
60. C. Pereira-Nabais, J. Światowska, A. Chagnes, F. Ozanam, A. Gohier, P. Tran-Van, C.-S. Cojocaru, M. Cassir and P. Marcus, Appl. Surf. Sci., 2013, 266, 5–16.

---