aabills · Oct 22, 2026 · Oct 28, 2026 · Oct 28, 2026 · Oct 28, 2026 · Oct 28, 2026
@@ -7,6 +7,7 @@
 
 ## Features
 
+- Added `"coupled ion-electron transfer"` as an `"intercalation kinetics"` option, following Fraggedakis et al. (2021). It uses the existing reorganization energy and a new temperature-dependent `"... electrode coupled ion-electron transfer rate constant [A.m-2]"` per electrode, and its current saturates at large overpotential. ([#5993](https://github.com/pybamm-team/PyBaMM/pull/5993))
 - Replaced Codacy with open checks that also run locally. Ruff's `flake8-bandit` (`S`) rules replace Bandit, and CodeQL scans the Python and C++ code. `nox -s type-check` runs Pyrefly over `pybamm`, which now type-checks cleanly, and a new type error fails CI. clang-tidy reports solver findings to code scanning as an advisory check. ([#5974](https://github.com/pybamm-team/PyBaMM/pull/5974))
 
 ## Bug fixes

@@ -0,0 +1,7 @@
+Coupled ion-electron transfer
+=============================
+
+.. autoclass:: pybamm.kinetics.CoupledIonElectronTransfer
+    :members:
+
+.. footbibliography::
@@ -5,6 +5,7 @@ Kinetics
 
   base_kinetics
   butler_volmer
+  coupled_ion_electron_transfer
   diffusion_limited
   linear
   marcus

@@ -251,6 +251,17 @@ @article{Fan2022
   publisher = {Elsevier}
 }
 
+@article{Fraggedakis2021,
+  title     = {Theory of coupled ion-electron transfer kinetics},
+  author    = {Fraggedakis, Dimitrios and McEldrew, Michael and Smith, Raymond B and Krishnan, Yamini and Zhang, Yirui and Bai, Peng and Chueh, William C and Shao-Horn, Yang and Bazant, Martin Z},
+  journal   = {Electrochimica Acta},
+  volume    = {367},
+  pages     = {137432},
+  year      = {2021},
+  publisher = {Elsevier},
+  doi       = {10.1016/j.electacta.2020.137432}
+}
+
 @article{Gustafsson2020,
   doi       = {10.21105/joss.02369},
   year      = {2020},
@@ -904,6 +915,17 @@ @article{Yang2017
   doi       = {10.1016/j.jpowsour.2017.05.110}
 }
 
+@article{Zeng2014,
+  title     = {Simple formula for {Marcus--Hush--Chidsey} kinetics},
+  author    = {Zeng, Yi and Smith, Raymond B and Bai, Peng and Bazant, Martin Z},
+  journal   = {Journal of Electroanalytical Chemistry},
+  volume    = {735},
+  pages     = {77--83},
+  year      = {2014},
+  publisher = {Elsevier},
+  doi       = {10.1016/j.jelechem.2014.09.038}
+}
+
 @article{Zhao2018,
   title     = {Modeling the effects of thermal gradients induced by tab and surface cooling on lithium ion cell performance},
   author    = {Zhao, Yan and Patel, Yatish and Zhang, Teng and Offer, Gregory J},

@@ -160,6 +160,8 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 Model for intercalation kinetics. Can be "symmetric Butler-Volmer"
                 (default), "asymmetric Butler-Volmer", "linear", "Marcus",
                 "Marcus-Hush-Chidsey" (which uses the asymptotic form from Zeng 2014),
+                "coupled ion-electron transfer" (which uses the form from Fraggedakis
+                2021),
                 or "MSMR" (which uses the form from Baker 2018). A 2-tuple can be
                 provided for different behaviour in negative and positive electrodes.
             * "interface utilisation": str
@@ -353,6 +355,7 @@ def __init__(self, extra_options):
                 "linear",
                 "Marcus",
                 "Marcus-Hush-Chidsey",
+                "coupled ion-electron transfer",
                 "MSMR",
             ],
             "interface utilisation": ["full", "constant", "current-driven"],
@@ -1388,6 +1391,8 @@ def get_intercalation_kinetics(self, domain):
             return pybamm.kinetics.Marcus
         elif options["intercalation kinetics"] == "Marcus-Hush-Chidsey":
             return pybamm.kinetics.MarcusHushChidsey
+        elif options["intercalation kinetics"] == "coupled ion-electron transfer":
+            return pybamm.kinetics.CoupledIonElectronTransfer
         elif options["intercalation kinetics"] == "MSMR":
             return pybamm.kinetics.MSMRButlerVolmer
 

@@ -49,6 +49,65 @@ def __init__(self, param, domain, reaction, options=None, phase="primary"):
         domain_options = getattr(self.options, domain)
         self.size_distribution = domain_options["particle size"] == "distribution"
 
+    def _get_intercalation_concentrations(self, variables):
+        """
+        Get the electrolyte concentration, particle surface concentration and
+        temperature seen by the main intercalation reaction.
+
+        Parameters
+        ----------
+        variables : dict
+            The variables in the full model.
+
+        Returns
+        -------
+        c_e, c_s_surf, T : :class:`pybamm.Symbol`
+            Electrolyte concentration, particle surface concentration and
+            temperature, on matching domains.
+        """
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+        c_e = variables[f"{Domain} electrolyte concentration [mol.m-3]"]
+        T = variables[f"{Domain} electrode temperature [K]"]
+
+        # For "particle-size distribution" submodels, take distribution version
+        # of c_s_surf that depends on particle size.
+        if self.size_distribution:
+            c_s_surf = variables[
+                f"{Domain} {phase_name}particle surface "
+                "concentration distribution [mol.m-3]"
+            ]
+            # If all variables were broadcast (in "x"), take only the orphans,
+            # then re-broadcast c_e
+            if (
+                isinstance(c_s_surf, pybamm.Broadcast)
+                and isinstance(c_e, pybamm.Broadcast)
+                and isinstance(T, pybamm.Broadcast)
+            ):
+                c_s_surf = c_s_surf.orphans[0]
+                c_e = c_e.orphans[0]
+                T = T.orphans[0]
+
+                # as c_e must now be a scalar, re-broadcast to "current collector"
+                c_e = pybamm.PrimaryBroadcast(c_e, ["current collector"])
+            # broadcast c_e, T onto "particle size"
+            c_e = pybamm.PrimaryBroadcast(c_e, [f"{domain} {phase_name}particle size"])
+            T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"])
+        else:
+            c_s_surf = variables[
+                f"{Domain} {phase_name}particle surface concentration [mol.m-3]"
+            ]
+            # If all variables were broadcast, take only the orphans
+            if (
+                isinstance(c_s_surf, pybamm.Broadcast)
+                and isinstance(c_e, pybamm.Broadcast)
+                and isinstance(T, pybamm.Broadcast)
+            ):
+                c_s_surf = c_s_surf.orphans[0]
+                c_e = c_e.orphans[0]
+                T = T.orphans[0]
+        return c_e, c_s_surf, T
+
     def _get_exchange_current_density(self, variables):
         """
         A private function to obtain the exchange current density
@@ -65,54 +124,13 @@ def _get_exchange_current_density(self, variables):
         """
         phase_param = self.phase_param
         domain, Domain = self.domain_Domain
-        phase_name = self.phase_name
         domain_options = getattr(self.options, domain)
 
         c_e = variables[f"{Domain} electrolyte concentration [mol.m-3]"]
         T = variables[f"{Domain} electrode temperature [K]"]
 
         if self.reaction == "lithium-ion main":
-            # For "particle-size distribution" submodels, take distribution version
-            # of c_s_surf that depends on particle size.
-            domain_options = getattr(self.options, domain)
-            if domain_options["particle size"] == "distribution":
-                c_s_surf = variables[
-                    f"{Domain} {phase_name}particle surface "
-                    "concentration distribution [mol.m-3]"
-                ]
-                # If all variables were broadcast (in "x"), take only the orphans,
-                # then re-broadcast c_e
-                if (
-                    isinstance(c_s_surf, pybamm.Broadcast)
-                    and isinstance(c_e, pybamm.Broadcast)
-                    and isinstance(T, pybamm.Broadcast)
-                ):
-                    c_s_surf = c_s_surf.orphans[0]
-                    c_e = c_e.orphans[0]
-                    T = T.orphans[0]
-
-                    # as c_e must now be a scalar, re-broadcast to
-                    # "current collector"
-                    c_e = pybamm.PrimaryBroadcast(c_e, ["current collector"])
-                # broadcast c_e, T onto "particle size"
-                c_e = pybamm.PrimaryBroadcast(
-                    c_e, [f"{domain} {phase_name}particle size"]
-                )
-                T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"])
-
-            else:
-                c_s_surf = variables[
-                    f"{Domain} {phase_name}particle surface concentration [mol.m-3]"
-                ]
-                # If all variables were broadcast, take only the orphans
-                if (
-                    isinstance(c_s_surf, pybamm.Broadcast)
-                    and isinstance(c_e, pybamm.Broadcast)
-                    and isinstance(T, pybamm.Broadcast)
-                ):
-                    c_s_surf = c_s_surf.orphans[0]
-                    c_e = c_e.orphans[0]
-                    T = T.orphans[0]
+            c_e, c_s_surf, T = self._get_intercalation_concentrations(variables)
             # Get main reaction exchange-current density (may have empirical hysteresis)
             j0_option = getattr(domain_options, self.phase)["exchange-current density"]
             if j0_option == "single":

@@ -3,6 +3,7 @@
 from .butler_volmer import SymmetricButlerVolmer, AsymmetricButlerVolmer
 from .linear import Linear
 from .marcus import Marcus, MarcusHushChidsey
+from .coupled_ion_electron_transfer import CoupledIonElectronTransfer
 from .tafel import ForwardTafel  # , BackwardTafel
 from .no_reaction import NoReaction
 from .msmr_butler_volmer import MSMRButlerVolmer
@@ -15,6 +16,7 @@
 from .inverse_kinetics.inverse_butler_volmer import InverseButlerVolmer
 from .inverse_kinetics.inverse_linear import InverseLinear
 
-__all__ = ['base_kinetics', 'butler_volmer', 'diffusion_limited',
+__all__ = ['base_kinetics', 'butler_volmer', 'coupled_ion_electron_transfer',
+           'diffusion_limited',
            'inverse_kinetics', 'linear', 'marcus', 'msmr_butler_volmer',
            'no_reaction', 'tafel', 'total_main_kinetics']
@@ -0,0 +1,123 @@
+#
+# Coupled ion-electron transfer kinetics class
+#
+from __future__ import annotations
+
+import numpy as np
+
+import pybamm
+
+from .base_kinetics import BaseKinetics
+
+
+class CoupledIonElectronTransfer(BaseKinetics):
+    """
+    Intercalation rate law in which the ion and electron transfer in one concerted
+    step :footcite:t:`Fraggedakis2021`, using the asymptotic Marcus-Hush-Chidsey
+    integral of :footcite:t:`Zeng2014`. In units of :math:`RT/F`,
+
+    .. math::
+        j = k_0 \\sqrt{\\pi\\lambda}\\,\\mathrm{erfc}\\left(\\frac{\\lambda - \\sqrt{1 + \\sqrt{\\lambda} + \\eta_f^2}}{2\\sqrt{\\lambda}}\\right) \\left[\\frac{\\theta}{1 + e^{-\\eta_f}} - \\frac{\\tilde{c}_e (1 - \\theta)}{1 + e^{\\eta_f}}\\right],
+
+    where :math:`\\eta_f = \\eta_r + \\ln(\\tilde{c}_e (1 - \\theta) / \\theta)`,
+    :math:`\\theta = c_{s,surf}/c_{max}` and :math:`\\tilde{c}_e = c_e/c_{e,init}`.
+    Not available for a lithium metal electrode, which has no host lattice.
+
+    Parameters
+    ----------
+    param : parameter class
+        model parameters
+    domain : str
+        The domain to implement the model, either: 'Negative' or 'Positive'.
+    reaction : str
+        The name of the reaction being implemented
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
+    phase : str, optional
+        Phase of the particle (default is "primary")
+    """
+
+    def __init__(self, param, domain, reaction, options, phase="primary"):
+        if reaction == "lithium metal plating":
+            raise pybamm.OptionError(
+                "Coupled ion-electron transfer kinetics need a host lattice, so "
+                "they cannot be used for the lithium metal electrode of a half-cell. "
+                "Set 'intercalation kinetics' per electrode instead, e.g. "
+                "('Marcus-Hush-Chidsey', 'coupled ion-electron transfer') for a "
+                "positive working electrode"
+            )
+        if reaction != "lithium-ion main":
+            raise pybamm.OptionError(
+                "Coupled ion-electron transfer kinetics are only implemented for "
+                f"lithium-ion intercalation, not '{reaction}'"
+            )
+        super().__init__(param, domain, reaction, options, phase)
+        domain_options = getattr(self.options, self.domain)
+        if getattr(domain_options, self.phase)["exchange-current density"] != "single":
+            raise pybamm.OptionError(
+                "Coupled ion-electron transfer kinetics do not use the "
+                "exchange-current density parameter, so 'exchange-current density' "
+                "must be 'single'"
+            )
+        pybamm.citations.register("Fraggedakis2021")
+        pybamm.citations.register("Zeng2014")
+
+    def _get_exchange_current_density(self, variables):
+        c_e, c_s_surf, T = self._get_intercalation_concentrations(variables)
+        phase_param = self.phase_param
+        tol_c_s = pybamm.settings.tolerances["j0__c_s"]
+        tol_c_e = pybamm.settings.tolerances["j0__c_e"]
+        stoichiometry = pybamm.maximum(
+            pybamm.minimum(c_s_surf / phase_param.c_max, 1 - tol_c_s), tol_c_s
+        )
+        relative_c_e = pybamm.maximum(c_e, tol_c_e) / phase_param.c_e_init
+
+        # Kept for the rate law, which only receives j0 and the overpotential
+        self._stoichiometry = stoichiometry
+        self._relative_c_e = relative_c_e
+
+        # Linear-response slope, so j ~ j0 F eta_r / RT as for Butler-Volmer
+        F_RT = self.param.F / (self.param.R * T)
+        eta_f_equilibrium = self._formal_shift(stoichiometry, relative_c_e)
+        return (
+            relative_c_e
+            * (1 - stoichiometry)
+            * self._rate_factor(eta_f_equilibrium, F_RT, T)
+            * self._fermi(-eta_f_equilibrium)
+        )
+
+    @staticmethod
+    def _formal_shift(stoichiometry, relative_c_e):
+        return pybamm.log(relative_c_e * (1 - stoichiometry) / stoichiometry)
+
+    @staticmethod
+    def _fermi(x):
+        """Return 1 / (1 + exp(-x)) without overflow."""
+        return (1 + pybamm.tanh(x / 2)) / 2
+
+    def _rate_factor(self, eta_f, F_RT, T):
+        """Return k0 times the Fermi-weighted Marcus integral prefactor."""
+        lambda_T = F_RT * self.phase_param.mhc_lambda
+        k0 = self.phase_param.k0_coupled_ion_electron_transfer(T)
+        a = 1 + pybamm.sqrt(lambda_T)
+        arg = (lambda_T - pybamm.sqrt(a + eta_f**2)) / (2 * pybamm.sqrt(lambda_T))
+        return k0 * pybamm.sqrt(np.pi * lambda_T) * pybamm.erfc(arg)
+
+    def _match_domains(self, symbol, eta_r):
+        domain = self.domain
+        if symbol.secondary_domain == [
+            "current collector"
+        ] and eta_r.secondary_domain == [f"{domain} electrode"]:
+            symbol = pybamm.SecondaryBroadcast(symbol, [f"{domain} electrode"])
+        return symbol
+
+    def _get_kinetics(self, j0, ne, eta_r, T, u):
+        stoichiometry = self._match_domains(self._stoichiometry, eta_r)
+        relative_c_e = self._match_domains(self._relative_c_e, eta_r)
+
+        F_RT = self.param.F / (self.param.R * T)
+        eta_f = F_RT * eta_r + self._formal_shift(stoichiometry, relative_c_e)
+        oxidation = stoichiometry * self._fermi(eta_f)
+        reduction = relative_c_e * (1 - stoichiometry) * self._fermi(-eta_f)
+        return u * self._rate_factor(eta_f, F_RT, T) * (oxidation - reduction)
@@ -661,6 +661,16 @@ def j0(self, c_e, c_s_surf, T, lithiation=None):
             post_processor=regulariser,
         )
 
+    def k0_coupled_ion_electron_transfer(self, T):
+        """Coupled ion-electron transfer rate constant [A.m-2]"""
+        Domain = self.domain.capitalize()
+        inputs = {"Temperature [K]": T}
+        return pybamm.FunctionParameter(
+            f"{self.phase_prefactor}{Domain} electrode coupled ion-electron "
+            "transfer rate constant [A.m-2]",
+            inputs,
+        )
+
     def j0_stripping(self, c_e, c_Li, T):
         """Dimensional exchange-current density for stripping [A.m-2]"""
         tol = pybamm.settings.tolerances["j0__c_e"]

@@ -159,6 +159,22 @@ def test_kinetics_mhc(self):
         )
         self.run_basic_processing_test(options, parameter_values=parameter_values)
 
+    def test_kinetics_coupled_ion_electron_transfer(self):
+        options = {"intercalation kinetics": "coupled ion-electron transfer"}
+        parameter_values = pybamm.ParameterValues("Marquis2019")
+        parameter_values.update(
+            {
+                "Negative electrode reorganization energy [eV]": 0.2,
+                "Positive electrode reorganization energy [eV]": 0.2,
+                "Negative electrode coupled ion-electron transfer rate constant "
+                "[A.m-2]": 10,
+                "Positive electrode coupled ion-electron transfer rate constant "
+                "[A.m-2]": 10,
+            },
+            check_already_exists=False,
+        )
+        self.run_basic_processing_test(options, parameter_values=parameter_values)
+
     def test_kinetics_marcus(self):
         options = {"intercalation kinetics": "Marcus"}
         parameter_values = pybamm.ParameterValues("Marquis2019")

@@ -30,7 +30,7 @@
 'exchange-current density': 'single' (possible: ['single', 'current sigmoid'])
 'heat of mixing': 'false' (possible: ['false', 'true'])
 'hydrolysis': 'false' (possible: ['false', 'true'])
-'intercalation kinetics': 'symmetric Butler-Volmer' (possible: ['symmetric Butler-Volmer', 'asymmetric Butler-Volmer', 'linear', 'Marcus', 'Marcus-Hush-Chidsey', 'MSMR'])
+'intercalation kinetics': 'symmetric Butler-Volmer' (possible: ['symmetric Butler-Volmer', 'asymmetric Butler-Volmer', 'linear', 'Marcus', 'Marcus-Hush-Chidsey', 'coupled ion-electron transfer', 'MSMR'])
 'interface utilisation': 'full' (possible: ['full', 'constant', 'current-driven'])
 'lithium plating': 'none' (possible: ['none', 'reversible', 'partially reversible', 'irreversible'])
 'lithium plating porosity change': 'false' (possible: ['false', 'true'])