SurfaceReactionBC: statistical factor for mixed isotopologue channels cannot be expressed

_(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.

```python
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.)*

Start with the SurfaceReactionBC implementation and run the supplied isotope reproducer, then inspect test/system_tests/test_2_isotopes_no_pressure. Confirm whether the project wants a documentation-only convention or a changed rate calculation; done requires that choice to be explicit and covered by regression tests.

Written by the indexing model from the issue text.