Turing-region preservation in matrix-oriented splitting methods for reaction-diffusion systems

Technical Summary: Turing-Region Preservation in Matrix-Oriented Splitting Methods

Problem Statement
Reaction–diffusion systems are fundamental to modeling self-organization in biology, ecology, and chemistry, particularly through the mechanism of diffusion-driven (Turing) instability. While theoretical analyses often focus on the continuous level or semi-discrete spatial approximations, practical numerical simulations rely on fully discrete dynamical systems. The paper argues that consistency of a numerical scheme with the continuous model as discretization parameters tend to zero is insufficient to guarantee the correct reproduction of patterns. In long-time integrations required to distinguish transient structures from robust asymptotic patterns, the time integrator can induce discrete dynamics that differ substantially from the continuous flow. Specifically, a scheme may anticipate instability on modes that are continuously stable, suppress modes that are continuously unstable, or generate spurious destabilization mechanisms. The central problem addressed is whether fully discrete maps induced by time integrators preserve the "Turing region"—the set of parameters and modes for which the system is unstable—of the underlying continuous reaction–diffusion system.

Methodology
The authors develop a matrix-oriented formulation for first-order splitting integrators applied to two-species reaction–diffusion systems on two-dimensional rectangular domains. By exploiting tensor-product spatial discretizations, the system is rewritten as a differential matrix equation, avoiding the explicit construction of large Kronecker-product matrices.

The study focuses on two families of first-order splitting schemes:

1. IE-S: Implicit Euler treatment of the diffusive flow combined with various local reaction solvers.
2. EX-S: Exact treatment of the diffusive flow (via matrix exponentials) combined with the same reaction solvers.

The reaction substeps utilize a set of first-order local integrators ($S \in \mathcal{S}$), including Symplectic Euler (SE), Adjoint Symplectic Euler (ASE), Poisson Euler (PE), their adjoints, and explicit variants (EV). The Gierer–Meinhardt activator–inhibitor system serves as the canonical benchmark.

The core analytical tool is a modal analysis based on the amplification matrices of the fully discrete map. The authors derive Jury conditions (stability criteria for discrete maps) for each scheme. These conditions are decomposed into three quantities ($J_1, J_2, J_3$):

• $J_1$: Corresponds asymptotically (and exactly for some schemes) to the continuous Turing polynomial $h(\mu)$. Its violation signifies the diffusion-driven instability.
• $J_2$: A purely discrete stability constraint with no continuous counterpart. Its violation indicates spurious numerical instability.
• $J_3$: Related to homogeneous mode stability.

The analysis separates the continuous Turing mechanism (carried by $J_1$) from purely discrete effects (carried by $J_2$).

Key Contributions and Results

1. Derivation of Discrete Jury Conditions: The paper derives explicit Jury conditions for the IE-S and EX-S families. It establishes that for first-order splitting schemes, the first Jury quantity $J_1$ generally reproduces the sign of the continuous Turing polynomial $h(\mu)$ asymptotically as the time step $h_t \to 0$. However, the second Jury quantity $J_2$ introduces purely discrete stability constraints.

2. Identification of Two Opposite Pathologies: Through the Gierer–Meinhardt model, the authors exhibit two distinct failure modes:

   • Spurious Pattern Formation (IMEX): In a continuous Turing-stable regime, the IMEX scheme (IE-EE) can violate the second Jury condition ($J_2 < 0$) while satisfying the first ($J_1 > 0$). This generates a stable, non-homogeneous numerical pattern that does not exist in the continuous model.
   • Spurious Stability (Adjoint-Symplectic): In a continuous Turing-unstable regime, the adjoint-symplectic families (IE-ASE/APE and EX-ASE/APE) may fail to violate the first Jury condition on the unstable modes ($J_1 > 0$), thereby suppressing the pattern that the continuous model predicts.

3. Exact Preservation by IMEX: A significant theoretical finding is that for the IMEX scheme (IE-EE), the first Jury quantity $J_1$ reproduces the sign of the continuous Turing polynomial $h(\mu)$ exactly, not just asymptotically. This makes IMEX a reliable detector of continuous Turing modes via $J_1$. However, the paper demonstrates that this exact preservation does not preclude spurious patterns generated by the violation of $J_2$.

4. Sufficient Conditions for Preservation: For the IMEX scheme, the authors derive an explicit time-step condition that guarantees the preservation of the continuous Turing region. This condition requires that in a Turing-stable regime, the time step must be small enough to ensure the dominance of $J_2$ over $J_1$ (specifically, $J_2 \ge J_1$), preventing spurious instability. In a genuine Turing regime, the homogeneous stability bound is the primary constraint.

5. Comparative Analysis of Families:

   • IE-S vs. EX-S: The qualitative behavior (green/orange/blue regions in parameter space) is governed primarily by the reaction substep rather than the diffusion treatment (implicit vs. exact).
   • SE/PE Superiority: The Symplectic Euler (SE) and Poisson Euler (PE) families (and their explicit variants) are identified as the most balanced first-order solvers. They avoid the $J_2$-driven spurious patterns of IMEX in stable regimes and do not suffer from the $J_1$-driven suppression of patterns seen in adjoint-symplectic families in unstable regimes.
   • Adjoint Families: The ASE/APE families are shown to be prone to suppressing continuous instabilities, a defect encoded in the reaction amplification matrix that persists even with exact diffusion treatment.

Significance and Claims
The paper frames Turing-region preservation as a qualitative property of the continuous flow, aligning with the principles of geometric numerical integration. The authors claim that consistency alone does not guarantee that the discrete and continuous Turing regions coincide at the time steps used in practice.

The primary contribution is methodological: before a splitting scheme is used to interpret numerical patterns, its discrete Turing region must be compared with the continuous one. The paper provides the algebraic criteria (Jury conditions) to perform this check. It does not claim to have found a single "optimal" scheme that guarantees Turing-region preservation for all time steps; rather, it highlights that each integrator induces its own discrete Turing region. The analysis shows that controlling the second Jury condition is necessary to prevent spurious patterns in stable regimes, while the first Jury condition's behavior determines whether genuine instabilities are detected or suppressed.

The work concludes that while the SE/PE family performs best among the tested first-order methods, the design of integrators that preserve the continuous Turing region as a geometric structure remains an open challenge for future research, potentially involving positivity-preserving or structure-preserving nonstandard integrators.