Cross-Diffusion Reaction-Diffusion Systems

• Cross-diffusion reaction-diffusion systems are models where the flux of one component is influenced by others through non-diagonal, state-dependent diffusion processes.
• They modify classical Turing instability criteria by introducing off-diagonal transport, relaxing diffusion contrast conditions and shifting instability thresholds via geometry and heterogeneity.
• These systems are central to understanding complex pattern selection, dynamic front propagation, and nonlinear phenomena in ecology, chemotaxis, and biomedical applications.

A cross-diffusion reaction-diffusion system is a reaction-diffusion model in which the flux of one component is influenced by other components, so the diffusion operator is non-diagonal, state-dependent, or both. In linear form this appears through a diffusion matrix such as $$D=\begin{pmatrix}1&d_v\ d_u&d\end{pmatrix}$$; in nonlinear formulations it is often written through chemical potentials,

$$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$

The subject therefore includes classical cross-diffusion in ecology and chemotaxis, triangular systems in which only one equation contains cross-diffusion, effective Markovian reductions of anomalous transport, and fast-reaction limits that generate cross-diffusive fluxes. Its central consequence is that transport is no longer a purely diagonal smoothing mechanism: it can become the destabilizing agent behind Turing, wave, and localized instabilities, and it can reorganize nonlinear pattern selection, front propagation, and long-time dynamics (Villar-Sepúlveda et al., 2024, Aymard, 2020, Baron et al., 2018).

1. Mathematical forms and diffusion structures

A convenient general formulation is

$$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$

with no-flux conditions $\nabla\mu_i\cdot n=0$. In this representation, linear diffusion is encoded by $d_i$, self-diffusion by $d_{ii}\phi_i^{\alpha_i}$, and cross-diffusion by the off-diagonal terms $d_{ij}\phi_j^{\beta_{ij}}$. This formulation is especially useful when diffusion is written in chemical-potential form and an energy law is available (Aymard, 2020).

Two-species models often expose the same structure more directly. In the activator-depleted Schnakenberg setting on an annulus, the nondimensional system is

$$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$

with diffusion matrix

$$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$

Here the off-diagonal entries $d_u,d_v$ are linear cross-diffusion coefficients, and the determinant condition is the basic well-posedness constraint (Yigit et al., 2024).

Nonlinear state-dependent cross-diffusion is frequently expressed through a Jacobian of fluxes. In the BVAM-type system with self- and cross-diffusion,

$$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$0

the diffusion matrix obtained by linearizing the chemical potentials is

$$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$1

The off-diagonal entries are now state-dependent and symmetric, so diffusion coefficients depend on amplitude as well as species identity (Aymard, 2024).

A major subclass is the triangular system. In the form

$$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$2

$$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$3

the first equation contains self-diffusion and cross-diffusion, while the second contains only self-diffusion. This “triangular” structure is important analytically because coupling is asymmetric: $$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$4 is affected by $$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$5, but $$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$6 is not cross-diffused by $$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$7 (Guerand et al., 2022).

2. Linear instability, Turing mechanisms, and the role of geometry

For a homogeneous steady state $$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$8, the linearized spectral problem is

$$\frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,\qquad \mu_i=\left(d_i+d_{ii}\phi_i^{\alpha_i}+\sum_{j\neq i}d_{ij}\phi_j^{\beta_{ij}}\right)\phi_i.$$9

where $$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$0 and $$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$1. Within this framework, a Turing instability means that at least one real eigenvalue satisfies $$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$2 for some $$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$3, while a wave instability means that a complex conjugate pair satisfies $$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$4 with nonzero imaginary part. The essential point is that cross-diffusion changes the matrix pencil $$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$5, so diffusion-driven instability is no longer governed by the classical diagonal-diffusion criteria (Villar-Sepúlveda et al., 2024).

Several model classes show that cross-diffusion is not merely a perturbation of ordinary diffusion. In the Lotka-Volterra predator-prey model

$$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$6

linear instability requires

$$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$7

In the terminology of that analysis, linear diffusion of the predator is stabilizing, while the cross-diffusion term is the only possible destabilizer. The bifurcation at onset is stationary, not oscillatory, because $$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$8 with zero imaginary part (Tulumello et al., 2013).

The Schnakenberg system with linear cross-diffusion makes the same point in a different way. The modified Turing conditions include

$$\forall i\in\{1,\dots,M\},\quad \frac{\partial \phi_i}{\partial t}-\nabla^2\mu_i=R_i,$$9

$\nabla\mu_i\cdot n=0$0

$\nabla\mu_i\cdot n=0$1

In that setting, cross-diffusion can produce Turing instability even when the activator diffuses as fast as or faster than the inhibitor, and pattern formation is possible even at $\nabla\mu_i\cdot n=0$2, where the classical system would not Turing-unstable. The asymmetry of the two cross-diffusion channels is explicit: increasing $\nabla\mu_i\cdot n=0$3 lowers the threshold $\nabla\mu_i\cdot n=0$4, whereas increasing $\nabla\mu_i\cdot n=0$5 raises it (Gambino et al., 2015).

Spatial heterogeneity and domain geometry modify the same instability logic rather than replacing it. In a one-dimensional heterogeneous cross-diffusion system with weak transport, a WKB ansatz shows that local Turing conditions are governed by

$\nabla\mu_i\cdot n=0$6

with pattern-forming regions determined by

$\nabla\mu_i\cdot n=0$7

and

$\nabla\mu_i\cdot n=0$8

The resulting instability is localized to subregions of the domain, rather than global (Gaffney et al., 2022). On a two-dimensional annulus, the Laplacian spectrum changes again, and the thickness $\nabla\mu_i\cdot n=0$9 enters directly into lower and upper bounds separating Turing-only behavior from Hopf/transcritical regimes (Yigit et al., 2024).

A common misconception is that the classical Turing picture survives unchanged once off-diagonal diffusion is added. The cited analyses instead show that cross-diffusion can be the sole destabilizer, can relax the usual diffusion-contrast requirement, and can shift instability thresholds through geometry, heterogeneity, and mode selection (Tulumello et al., 2013, Gambino et al., 2015, Gaffney et al., 2022).

3. Emergent cross-diffusion from anomalous transport and fast reactions

Cross-diffusion need not be introduced phenomenologically. One route is anomalous transport. In a reaction-subdiffusion system derived from a CTRW with Mittag-Leffler waiting times, the original equation is history dependent because reactions and transport are coupled through a memory term. The large-scale pattern-forming behavior can nevertheless be represented by an effective Markovian system

$d_i$0

with constitutive law

$d_i$1

The off-diagonal terms arise because the per-capita removal rate of one species depends on the concentration of the other. In this effective description, the original subdiffusive system and the Markovian cross-diffusion system have the same Turing instability and the same stationary patterns; when particles are short-lived, the transient dynamics are captured as well (Baron et al., 2018).

A second route is the fast-reaction limit. In a slow-fast competition model with two fast states $d_i$2 and a second species $d_i$3, the stiff switching term $d_i$4 enforces a local equilibrium

$d_i$5

As $d_i$6, the limit system is

$d_i$7

with

$d_i$8

This is a triangular cross-diffusion system: the $d_i$9-flux depends on $d_{ii}\phi_i^{\alpha_i}$0 through the equilibrium splitting, whereas the $d_{ii}\phi_i^{\alpha_i}$1-equation remains a linear heat equation with reaction (Brocchieri et al., 10 Mar 2025).

An analogous mechanism appears in the fast reversible reaction $d_{ii}\phi_i^{\alpha_i}$2. There, the three-species system converges to a two-variable cross-diffusion system in the conserved combinations $d_{ii}\phi_i^{\alpha_i}$3 and $d_{ii}\phi_i^{\alpha_i}$4, supplemented by the equilibrium constraint

$d_{ii}\phi_i^{\alpha_i}$5

The limiting system inherits an entropy structure, but the entropy is not the sum of the entropies of the reduced variables; it is inherited through the nonlinear equilibrium constraint (Daus et al., 2017).

These derivations imply that cross-diffusion is not restricted to explicit taxis or avoidance laws. It can emerge as an effective macroscopic transport law from memory effects or fast local equilibration (Baron et al., 2018, Brocchieri et al., 10 Mar 2025).

4. Nonlinear pattern selection, fronts, oscillations, and chaos

Near threshold, cross-diffusion systems are commonly reduced by multiple-scales analysis to amplitude equations. For simple critical eigenspaces, the generic reduction is the cubic Stuart-Landau equation

$d_{ii}\phi_i^{\alpha_i}$6

with $d_{ii}\phi_i^{\alpha_i}$7 corresponding to supercritical bifurcation and $d_{ii}\phi_i^{\alpha_i}$8 to subcritical bifurcation. In subcritical regimes, one must pass to a quintic normal form such as

$d_{ii}\phi_i^{\alpha_i}$9

In the Schnakenberg cross-diffusion model, this distinction organizes hysteresis, finite-amplitude patterned branches, and the transition from rolls to hexagons in resonant two-dimensional settings (Gambino et al., 2015).

Mode multiplicity and resonance determine the pattern family. In a competitive Lotka-Volterra system with nonlinear cross-diffusion on a rectangle, a simple critical eigenvalue leads to rolls or squares; a double non-resonant eigenvalue yields coupled Landau equations and mixed-mode states such as supersquares; and a double resonant eigenvalue yields quadratic-cubic amplitude systems supporting rolls and hexagons, together with bistability and hysteresis (Gambino et al., 2012). The same taxonomy reappears in the Schnakenberg model, where numerical simulations confirm rolls, squares, rhombi, rectangles, hexagons, and mixed-mode patterns (Gambino et al., 2015).

A recurrent subtlety is that stationary Turing onset does not preclude later temporal complexity. In the predator-prey model, linear theory predicts only a stationary Turing bifurcation, but numerical simulations in a large portion of the subcritical zone show oscillating patterns, a secondary Hopf bifurcation of the patterned branch, torus bifurcation, period doubling, and eventually chaos (Tulumello et al., 2013). A closely related sequence appears in the BVAM-type model augmented with self- and cross-diffusion: stable Turing patterns lose stability through a Hopf bifurcation, generating oscillating Turing patterns and periodic orbits; further variation of the control parameter produces period doubling and strange attractors, with phase portraits studied in both $d_{ij}\phi_j^{\beta_{ij}}$0-space and $d_{ij}\phi_j^{\beta_{ij}}$1-space (Aymard, 2024).

Cross-diffusion also alters front propagation. A minimal two-species model without self-diffusion,

$d_{ij}\phi_j^{\beta_{ij}}$2

admits exact monotone traveling fronts when the reaction term is quartic and the profile satisfies a quadratic reduction-of-order ansatz. The wave speed is uniquely fixed by

$d_{ij}\phi_j^{\beta_{ij}}$3

Depending on the root structure of the quartic, the front may resemble the Fisher-KPP stability pattern; however, the constructed monotone fronts are reported to be dynamically unstable and evolve into oscillatory or pulse-like structures (Aldurayhim et al., 2020).

5. Entropy structures, solvability theories, and segregation

A central analytical theme is the search for entropy or Lyapunov structures. For systems

$d_{ij}\phi_j^{\beta_{ij}}$4

a formal entropy identity has the form

$d_{ij}\phi_j^{\beta_{ij}}$5

and positivity of $d_{ij}\phi_j^{\beta_{ij}}$6 is the key coercivity mechanism. This viewpoint is strong enough to treat mixed convex/concave cross-diffusion, including the case $d_{ij}\phi_j^{\beta_{ij}}$7, $d_{ij}\phi_j^{\beta_{ij}}$8, where one cross-diffusion is convex and the other concave (Desvillettes et al., 2014).

For triangular population-dynamics systems, entropy and duality methods yield global weak solutions. In the class

$d_{ij}\phi_j^{\beta_{ij}}$9

$$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$0

the global existence theory covers low-regularity nonnegative initial data and includes the triangular SKT model as the case $$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$1 (Trescases, 2015). A more regular triangular theory proves that the unique local smooth solution given by Amann theorem is actually global and remains in

$$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$2

under explicit exponent restrictions (Guerand et al., 2022).

For arbitrary numbers of species, global renormalized solutions are available when the diffusion matrix is non-diagonal and generally neither symmetric nor positive semi-definite, provided an entropy-compatible structural condition such as weak cross-diffusion or detailed balance holds. In that setting, the reaction terms may satisfy reversible mass-action kinetics and need not obey any growth condition; the solution concept is renormalized precisely to absorb the lack of direct $$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$3 control on the reactions (Chen et al., 2017).

Not all solvability results rely on full entropy. The glioblastoma model with densities $$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$4 and $$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$5,

$$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$6

has only a partial entropy structure. The one-dimensional existence proof instead combines this partial structure with fully implicit time discretization, variational arguments, fixed points, and weighted $$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$7-type estimates that control degeneracy at $$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$8 (Burger et al., 2017).

A distinct line of work concerns segregation. In the one-dimensional system

$$u_t=\Delta_r u+d_v\Delta_r v+\gamma f(u,v),\qquad v_t=d_u\Delta_r u+d\Delta_r v+\gamma g(u,v),$$9

$$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$0

a variational splitting scheme combining reaction ODEs and optimal-transport JKO steps yields weak solutions and preserves segregation for initially segregated data when $$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$1, even in the presence of vacuum (Carrillo et al., 2017).

6. Energetic numerics, computational methods, and model applications

The energetic formulation is not only analytical; it also guides discretization. For symmetric cross-diffusion with $$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$2 and $$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$3, the free energy

$$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$4

satisfies

$$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$5

A discrete counterpart, combined with Heun reaction steps and Strang splitting,

$$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$6

was implemented in FreeFem++ with adaptive mesh refinement and used to study a modified Gray-Scott system. The reported patterns include ring-shaped travelling waves and complex Turing-like spot/maze patterns, and geometry was found to influence local and final pattern structure (Aymard, 2020).

For SKT-type systems with self- and cross-diffusion, a variable-step nonlinear operator-splitting method based on modified Douglass-Gunn factorization provides a second-order scheme in space and time. The paper reports that the overall work per step is $$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$7, with only six $$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$8 vectors needed in memory, and establishes invertibility and stability under a CFL-type restriction (1901.10069).

Spectral continuation methods have become important for time-periodic and chaotic regimes. In the BVAM-type cross-diffusion model on a periodic interval, a Fourier spectral method in space, RK4 in time, Newton-Krylov solvers for equilibria and periodic orbits, and Floquet multiplier computations were used to track Hopf points, period-doubling bifurcations, and strange attractors. The methodological emphasis there is that high accuracy is necessary to catch periodic orbits and perform their linear stability analysis via Floquet multipliers (Aymard, 2024).

Applications span several scientific domains. Population-dynamics models include competitive Lotka-Volterra systems, triangular SKT systems, and predator-prey systems in which predator motion is biased by prey gradients (Gambino et al., 2012, Trescases, 2015, Tulumello et al., 2013). Chemotaxis and epidemiology appear in Keller-Segel-type formulations and in a spatio-temporal Ross-Macdonald malaria model whose cross-diffusion was designed to induce prescribed Turing or wave instabilities (Villar-Sepúlveda et al., 2024). Biomedical applications include glioblastoma growth with size exclusion (Burger et al., 2017), localized patterns in heterogeneous cross-diffusion systems motivated by bacterial chemotaxis and ecological interactions (Gaffney et al., 2022), annular cross-diffusion patterns interpreted through hypoxic-tumour cross sections (Yigit et al., 2024), and a pandemic model in which the death density obeys

$$D=\begin{bmatrix}1&d_v\ d_u&d\end{bmatrix}, \qquad d-d_ud_v>0.$$9

so that the spatial redistribution of deaths is driven by the diffusion of the infected population rather than by $d_u,d_v$0 itself (Cherniha et al., 2020).

Taken together, these developments show that cross-diffusion reaction-diffusion systems are not a single model class but a transport principle. They encompass explicit off-diagonal flux laws, nonlinear chemical-potential formulations, effective descriptions of anomalous motion, and singular limits of multiscale reaction systems. Their characteristic mathematical features are non-diagonal transport, altered instability criteria, rich amplitude dynamics, and a strong dependence on entropy structure, geometry, and asymptotic regime.