SurfaceReactionBC: statistical factor for mixed isotopologue channels cannot be expressed

[RemDelaporteMathurin]

(generated by Claude because I didn't have time, will go back to it later)

Summary

Flagging a convention issue rather than a bug: SurfaceReactionBC has no way to express the statistical degeneracy of a mixed isotopologue channel, and the natural way to write H + D <-> HD is silently inconsistent with detailed balance. Same root cause as the InterfaceReaction discussion in #1222: the rate is the plain product of reactant concentrations, and stoichiometry is expressed only by repeating a species in reactant.

The physics

For recombination, the forward step carries a degeneracy that the reverse step does not: there are two ways to pick an (H, D) pair to recombine, but an HD molecule dissociates only one way. Writing the homonuclear channel as k_r c_H^2 absorbs the one-half that avoids double counting identical pairs, so consistency requires

k_r(HD) = 2 k_r(HH)

with equal k_d. This is what reproduces the classical P_HD^2 / (P_H2 P_D2) = 4, i.e. H2 : HD : D2 = 1 : 2 : 1 at equal H and D.

Nothing in the API expresses this. value_fenics = kd * P - kr * prod(reactants) and the flux is applied once per entry in reactant, so every structural knob (repeating a species, declaring the BC twice) multiplies the whole rate and scales forward and reverse together. The degeneracy is forward-only, so it can only go into k_r0.

Reproducer

A slab equilibrating with a gas that is itself at equilibrium: P_H2 = P_D2 = 1, P_HD = 2, so P_HD^2/(P_H2 P_D2) = 4. All rate constants equal except k_r0 of the mixed channel. The correct answer is c_H = c_D = sqrt(k_d/k_r) = 1, with every channel individually at rest.

import numpy as np
import festim as F


def run(k_r_HD):
    model = F.HydrogenTransportProblem()
    model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 101))
    mat = F.Material(name="mat", D_0=1, E_D=0)
    vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat)
    left = F.SurfaceSubdomain1D(id=1, x=0)
    model.subdomains = [vol, left]

    H, D = F.Species("H"), F.Species("D")
    model.species = [H, D]
    model.temperature = 500

    model.boundary_conditions = [
        F.SurfaceReactionBC(reactant=[H, H], gas_pressure=1.0,
                            k_r0=1.0, E_kr=0, k_d0=1.0, E_kd=0, subdomain=left),
        F.SurfaceReactionBC(reactant=[D, D], gas_pressure=1.0,
                            k_r0=1.0, E_kr=0, k_d0=1.0, E_kd=0, subdomain=left),
        F.SurfaceReactionBC(reactant=[H, D], gas_pressure=2.0,
                            k_r0=k_r_HD, E_kr=0, k_d0=1.0, E_kd=0, subdomain=left),
    ]

    model.settings = F.Settings(atol=1e-12, rtol=1e-12, final_time=100, stepsize=0.5)
    model.initialise()
    model.run()
    return H.post_processing_solution.x.array.max()


for label, k_r_HD in [("k_r(HD) = k_r(HH)", 1.0), ("k_r(HD) = 2 k_r(HH)", 2.0)]:
    print(f"{label:24s}  c_H = {run(k_r_HD):.6f}   (expected 1.0)")

k_r(HD) = k_r(HH)         c_H = 1.154701   (expected 1.0)
k_r(HD) = 2 k_r(HH)       c_H = 1.000000   (expected 1.0)

1.154701 = sqrt(4/3), a 15% overload. Worse than the number itself: at that composition the mixed channel has a nonzero net rate balanced by the two homonuclear ones, so the model sits in a permanent circulating exchange current, HD -> 1/2 H2 + 1/2 D2, at a gas composition that is already at chemical equilibrium. A user modelling isotope exchange would read that as a physical result.

Declaring the mixed BC twice ([H, D] and [D, H]) does not help, since it doubles the kd and kr terms together: it gives c_H = 1.224745, wrong in a different way.

Notes

• test/system_tests/test_2_isotopes_no_pressure uses k_r0 = 0.02 for HH and 0.01 for HD, which is the inverse of the relation above. It is harmless there because k_d0 = 0 removes the reverse step, so detailed balance never enters, but it does suggest the convention is currently unguided.
• The InterfaceReaction analogue was checked numerically: with equal forward constants a three-channel H/T set equilibrates to K_exch = 1, and to 4.000000 once k_HT+ = 2 k_HH+. Declaring the mixed channel twice leaves it at 1.

Possible responses

1. Document it. One line in the SurfaceReactionBC docstring stating that k_r0 must carry the statistical factor for mixed channels, with the factor 2 spelled out. Cheapest, and no existing input file changes meaning.

2. Divide the forward rate by the symmetry factor of the reactant multiset, prod(nu_i!), so that A + A gives 0.5 * kappa * c_A^2 and A + B gives kappa * c_A * c_B. This is the standard correction for counting each identical pair twice, and it makes equal rate constants across HH/HD/DD automatically correct:

   c_H2 = (kappa/2kd) c_H^2,  c_D2 = (kappa/2kd) c_D^2,  c_HD = (kappa/kd) c_H c_D
   => K = c_HD^2/(c_H2 c_D2) = 4
   

   The catch is that it halves the flux in the single-isotope case, from 2 k c^2 to k c^2, so every existing recombination coefficient shifts by a factor of 2. That collides with the pre-existing ambiguity in the literature between Gamma = K_r c^2 and Gamma = 2 K_r c^2, so it would need a release note and an explicit statement of which convention FESTIM adopts.

I would argue for (1) unless you want to settle the recombination-coefficient convention explicitly, in which case (2) is the cleaner end state. Happy to open a PR for (1).

(Edited: an earlier version of this issue also suggested passing stoichiometry as a dict, e.g. reactant={H: 1, D: 1}. That does not help. It carries the same information as the list form, since for an elementary reaction the rate-law exponents and the stoichiometric coefficients are the same numbers. The degeneracy is a separate constant on the forward term only.)