# Reactant-bimolecular instability without D-unstable cores

This standalone example implements *Reactant-bimolecular mass-action instability without D-unstable cores*. Original PDF SHA-256: `12479870603ad4aafcc2a7cef8d0c32f22a78f07e210ba83c4f904c249063dbf`.

A child selection assigns distinct consuming reactions to a subset of species. Its square stoichiometric matrix is D-nonunstable when every positive diagonal column scaling has eigenvalues with nonpositive real parts. That structural property holds for **all 25 children** of this example, yet the full mass-action Jacobian has an unstable complex pair. A different operating point of the same network is strictly stable. The example separates the reusable kinetic model, structural certificates, operating conditions and numerical trajectories.

## Run and reuse

```sh
python -m pip install -r requirements.txt
python -m unittest -v
python example.py --output outputs
```

Use Python 3.11 or later. Editable inputs at the top set circulation ratio T, inverse D concentration L, throughput s, dilution, local perturbation size/duration and sweep bounds. Defaults reproduce the manuscript's dimensionless witness. The default run writes exact source, spectral and child certificates; a numerical operating grid; moving-equilibrium perturbation data; local trajectory data; two figure pairs; and source/output hashes. It does not run Lean.

```python
from example import OperatingFamily, Quartic, DilutionCompletion

family = OperatingFamily()
stable = family.model(T=2, L=1)
unstable = family.model(T=100, L=100)
print(Quartic(stable.jacobian()).classification())
print(Quartic(unstable.jacobian()).classification())
reactor = DilutionCompletion(unstable, delta='1/2')
print(reactor.feed_concentrations)
```

`Reaction` stores integer complexes. `Source` owns the stoichiometric matrices, monomials and equilibrium-rate reconstruction. `MassAction` supplies exact symbolic fields/Jacobians and a numerical stiff-integration interface. `OperatingFamily` supplies the paper's complete positive stationary-flux cone and a closed-form equilibrium for independently edited rates. `ChildLattice`, `Quartic` and `DilutionCompletion` are separate analysis/composition components. There is no shared theory framework to configure.

## The kinetic model

The reactions, in rate-vector order, are:

1. `A+D -> 2C`
2. `B+C -> empty`
3. `2C -> D`
4. `2D -> A+4B+4C`
5. `empty -> A`
6. `empty -> D`

Each reaction consumes at most two molecules, and no species appears on both sides of a reaction. The fourth reaction produces nine molecules. Homodimer event rates are `k*C**2` and `k*D**2`; the factor two for consumption comes from stoichiometry. These are deterministic mass-action conventions, not unordered-pair stochastic propensities.

All positive stationary fluxes are `v=s*(T,4,T,1,T-1,2)` with `T>1`, `s>0`. For any positive composition x, rate reconstruction uses `k_j=v_j/product(x_i**Y_ij)`. For `x*=(1,1,1,1/L)`, this becomes `k=(s*T*L,4*s,s*T,s*L**2,s*(T-1),2*s)`. Thus changing L also changes the rates needed to keep the selected composition stationary; it is not a concentration perturbation at fixed rates. Throughput s rescales time.

Default `T=L=100`, `s=1` gives rates `(10000,4,100,10000,99,2)`, equilibrium `(1,1,1,1/100)` and flux `(100,4,100,1,99,2)`. The exact characteristic polynomial is

`p(z)=z**4+10908*z**3+185600*z**2+1280000*z+64000000`.

Its coefficients are positive but the third Hurwitz determinant is `-5025252352000000`. The exact shifted-root construction evaluates that determinant for `p(z+t)` at t=2 and t=4, giving opposite signs. At an intermediate zero, `omega=sqrt(C(t)/A(t))` yields a root `t+i*omega`. This certifies a positive real part between 2 and 4 without relying on numerical eigenvalues. Numerical evaluation locates the pair at approximately `2.99974676 +/- 15.68952573 i`.

`Quartic.classification` uses the paper's positive-coefficient criterion. If its hypotheses do not hold, it returns `unknown`; a failed sufficient check is not treated as stability. For positive coefficients, positive H gives strict stability, negative H gives a complex-pair instability, and H=0 gives an imaginary-pair boundary. The shifted-root certificate in the output is explicitly the fixed manuscript witness, even if editable operating inputs are changed.

## All children, including boundary cases

`ChildLattice` enumerates every injective consuming assignment, including the empty assignment. Absent species are padded with a decoupled diagonal -1, which adds negative eigenvalues under positive column scaling. For each of the 25 patterns it independently derives the scaled quartic and its Hurwitz polynomial in positive variables `(a,b,c,d)`.

The code checks nonnegative polynomial coefficients and the required strict positivity of the first and third characteristic coefficients. This is an exact all-positive-scalings certificate via the paper's quartic exclusion lemma, not a sample of diagonal matrices. The saved JSON includes expanded polynomial expressions and monomial/coefficient data. A generic changed network can produce an inconclusive certificate; the code does not pretend to be a complete D-stability decision algorithm. Enumeration has an explicit budget and stops without an absence-of-core claim if exceeded.

Four children are singular. They are D-nonunstable, **not strictly D-stable**. Every original child is also checked to be a principal restriction of one of the four maximal selections. Empty children are retained in the count. The implementation rechecks finite algebraic certificates; it does not rebuild the formal theorem or replace its spectral lemmas.

## Operating map and the fast D mechanism

At `(T,L)=(2,1)` the same source is strictly stable, with characteristic coefficients `(24,124,288,256)` and H=626688. At T=100 and every L>=100, the exact polynomial `H(100,100+u)` has strictly negative coefficients, establishing an unstable ray. The figure's remaining shading and contour come from numerical evaluation, not an interval-certified bifurcation diagram.

Eliminating the fast D coordinate by a Schur complement gives a three-dimensional reduced matrix. Its cubic Hurwitz expression changes sign at `T*=(41+sqrt(2577))/4`, about 22.941. This is a **large-L threshold**, not the finite-L boundary. The reduction combines reaction sensitivities through the shared stationary flux; it is not a three-species child matrix. All 25 original structural certificates remain unchanged across the operating map.

The short-time illustration starts along the real part of a numerically computed mode, with a deterministic phase convention and perturbation size `1e-7`. Radau and BDF integrations are compared, and the nonlinear displacement is compared with `epsilon*Re(exp(lambda*t)*w)`. The growing spiral shows departure from an equilibrium. Neither it nor the unstable pair proves a Hopf bifurcation, a sustained oscillation, boundedness or persistence. If edited parameters select a stable mode, the routine can instead illustrate decay; its numerical eigenvalue and errors are always reported.

## Independent rate changes and a flow reactor

For the six original positive rate constants, the complete cone determines `s=k5/2` and `T=1+k4/s` (zero-based indices). Then

`D*=sqrt(s/k3)`, `C*=sqrt(s*T/k2)`, `B*=4*s/(k1*C*)`, `A*=s*T/(k0*D*)`.

`from_rates` uses this direct consequence of the paper's cone to construct the unique positive equilibrium and verifies stationarity exactly. The perturbation CSV changes each rate separately by +/-0.1%, moves the equilibrium accordingly and reports numerical spectral abscissae. These twelve probes are not a certificate for an entire six-dimensional box. The paper's rate-space robustness result is qualitative and local.

`DilutionCompletion` composes the first four reactions with four feeds and four first-order losses. With dilution delta, adjusted inflows are the old feed vector plus `delta*x*`, so the vector field is `F_delta(x)=F(x)+delta*(x*-x)` and the Jacobian is `J-delta*I`. Default delta=1 gives feed concentrations `(100,1,1,201/100)`. For the manuscript witness, every `0<delta<2` remains unstable by the exact root margin; with throughput s the corresponding sufficient bound is `delta<2*s`. Other edited operating points receive a fresh quartic classification.

The completed twelve-reaction source has 121 child selections. For each, the implementation identifies assigned loss reactions, checks their diagonal -1 columns, verifies the retained original child block, and transports its certificate through the resulting triangular matrix. Added inflows cannot be selected. This argument covers this completion; adding arbitrary reactions or changing complexes can introduce new children and requires fresh analysis.

## Scope and validation

These are abstract-species kinetics. The active reactions have no positive mass vector: `2C -> D` and `2D -> A+4B+4C` would require `mass_A+4*mass_B=0`. Feeding and dilution do not repair that obstruction. The model is not an atom-resolved chemical mechanism, and the two-sided bimolecular question remains open because the fourth product complex has nine molecules.

Concentrations and time are dimensionless. For chosen scales C0 and tau0, a reaction of reactant molecularity mu has physical rate `khat/(tau0*C0**(mu-1))`; this is a unit conversion, not experimental calibration. Stochastic pair propensities require a separate convention conversion.

Seven test groups verify source semantics, the complete flux cone, exact characteristic/shift identities, all 25 certificates and maximal coverage, the unstable ray and Schur reduction, all 121 completion transports with independent scaled determinant checks, moving-equilibrium reconstruction, the mass obstruction and local solver comparisons. No Lean proof is rerun. The paper's general minimality results—one reactant molecule cannot support an unstable positive equilibrium, and at least four species are needed for this counterexample—are explained by the manuscript, not inferred from this finite example.

Licensing: MIT is proposed, pending the owner's decision. No new license grant is made by this example.