---
title: Cross-Diffusion Reaction-Diffusion Systems
url: https://www.emergentmind.com/topics/cross-diffusion-reaction-diffusion-system
type: topic
---

# Cross-Diffusion Reaction-Diffusion Systems

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 [2409.06860] [2001.10069] [1811.04054].

## 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 [2001.10069].

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 [2412.20097].

Nonlinear state-dependent cross-diffusion is frequently expressed through a Jacobian of fluxes. In the BVAM-type system with self- and cross-diffusion,
\[
\mu_1=d_1u_1+d_{11}u_1^3+d_{12}u_2^2u_1,\qquad
\mu_2=d_2u_2+d_{22}u_2^3+d_{12}u_1^2u_2,
\]
the diffusion matrix obtained by linearizing the chemical potentials is
\[
D(u_1,u_2)=
\begin{pmatrix}
d_1+3d_{11}u_1^2+d_{12}u_2^2 & 2d_{12}u_1u_2\\
2d_{12}u_1u_2 & d_2+3d_{22}u_2^2+d_{12}u_1^2
\end{pmatrix}.
\]
The off-diagonal entries are now state-dependent and symmetric, so diffusion coefficients depend on amplitude as well as species identity [2412.17076].

A major subclass is the triangular system. In the form
\[
\partial_t u-\Delta\Big((d_u+d_{\alpha}\langle u\rangle^\alpha+d_{\beta}\langle v\rangle^\beta)u\Big)
=u(r_u-r_a\langle u\rangle^a-r_b\langle v\rangle^b),
\]
\[
\partial_t v-\Delta\Big((d_v+d_{\gamma}\langle v\rangle^\gamma)v\Big)
=v(r_v-r_c\langle v\rangle^c-r_d\langle u\rangle^d),
\]
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: \(u\) is affected by \(v\), but \(v\) is not cross-diffused by \(u\) [2202.10256].

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

For a homogeneous steady state \(P\), the linearized spectral problem is
\[
\det(J-\mu D-\lambda I)=0,\qquad \mu=k^2,
\]
where \(J=\partial_{\mathbf u}f(P)\) and \(D=D(P)\). Within this framework, a Turing instability means that at least one real eigenvalue satisfies \(\lambda>0\) for some \(\mu>0\), while a wave instability means that a complex conjugate pair satisfies \(\Re(\lambda)>0\) with nonzero imaginary part. The essential point is that cross-diffusion changes the matrix pencil \(J-\mu D\), so diffusion-driven instability is no longer governed by the classical diagonal-diffusion criteria [2409.06860].

Several model classes show that cross-diffusion is not merely a perturbation of ordinary diffusion. In the Lotka-Volterra predator-prey model
\[
\partial_t u=\Gamma u(r-\gamma u-v)+\nabla^2u,\qquad
\partial_t v=\Gamma v(-1+u)+d_{21}\nabla\cdot(v\nabla u)+d_2\nabla^2v,
\]
linear instability requires
\[
q<0\quad\Longleftrightarrow\quad d_{21}>\frac{\gamma d_2}{\det(J)}.
\]
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 \(\lambda(k_c)=0\) with zero imaginary part [1311.1748].

The Schnakenberg system with linear cross-diffusion makes the same point in a different way. The modified Turing conditions include
\[
d-d_ud_v>0,
\]
\[
dJ_{11}+J_{22}-d_uJ_{12}-d_vJ_{21}>0,
\]
\[
\bigl(dJ_{11}+J_{22}-d_uJ_{12}-d_vJ_{21}\bigr)^2
-4(d-d_ud_v)(J_{11}J_{22}-J_{12}J_{21})>0.
\]
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 \(d=1\), where the classical system would not Turing-unstable. The asymmetry of the two cross-diffusion channels is explicit: increasing \(d_v\) lowers the threshold \(d_c\), whereas increasing \(d_u\) raises it [1501.04890].

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
\[
B_0(x)=D^{-1}(x)J(x),
\]
with pattern-forming regions determined by
\[
\operatorname{tr}(J(x))<0,\qquad \det(J(x))>0,
\]
and
\[
\operatorname{tr}(B_0(x))>0,\qquad 
\operatorname{tr}(B_0(x))^2-4\det(B_0(x))>0.
\]
The resulting instability is localized to subregions of the domain, rather than global [2210.10155]. On a two-dimensional annulus, the Laplacian spectrum changes again, and the thickness \(\rho=b-a\) enters directly into lower and upper bounds separating Turing-only behavior from Hopf/transcritical regimes [2412.20097].

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 [1311.1748] [1501.04890] [2210.10155].

## 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
\[
\frac{\partial \rho_i}{\partial t}
=\frac{\partial^2}{\partial x^2}\big[\hat D_i(\rho_1,\rho_2)\rho_i\big]+f_i
=\frac{\partial}{\partial x}\left[D_{i1}\frac{\partial \rho_1}{\partial x}
+D_{i2}\frac{\partial \rho_2}{\partial x}\right]+f_i,
\]
with constitutive law
\[
\hat D_i(x,t)=\frac{\sigma_i^2}{\eta_i^{\alpha_i}
\left\{\frac{R_i^-[x,t]}{\rho_i(x,t)}\right\}^{1-\alpha_i}}.
\]
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 [1811.04054].

A second route is the fast-reaction limit. In a slow-fast competition model with two fast states \(u_a^\varepsilon,u_b^\varepsilon\) and a second species \(v^\varepsilon\), the stiff switching term \(\varepsilon^{-1}Q\) enforces a local equilibrium
\[
Q(u_a^\ast,u_b^\ast,v)=0,\qquad u_a^\ast+u_b^\ast=u.
\]
As \(\varepsilon\to0\), the limit system is
\[
\partial_t u-\Delta(A(u,v))=F_u(u,v),\qquad
\partial_t v-\Delta(d_v v)=F_v(u,v),
\]
with
\[
A(u,v)=d_a u_a^\ast(u,v)+d_b u_b^\ast(u,v).
\]
This is a triangular cross-diffusion system: the \(u\)-flux depends on \(v\) through the equilibrium splitting, whereas the \(v\)-equation remains a linear heat equation with reaction [2503.07156].

An analogous mechanism appears in the fast reversible reaction \(A\rightleftharpoons B+C\). There, the three-species system converges to a two-variable cross-diffusion system in the conserved combinations \(u_1+u_2\) and \(u_1+u_3\), supplemented by the equilibrium constraint
\[
q_1(u_1)=q_2(u_2)q_3(u_3).
\]
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 [1710.03590].

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 [1811.04054] [2503.07156].

## 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
\[
\frac{dA}{dT}=\sigma A-LA^3,
\]
with \(L>0\) corresponding to supercritical bifurcation and \(L<0\) to subcritical bifurcation. In subcritical regimes, one must pass to a quintic normal form such as
\[
\frac{dA}{dT}=\bar{\sigma}A-\bar{L}A^3+\bar{R}A^5.
\]
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 [1501.04890].

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 [1211.4412]. The same taxonomy reappears in the Schnakenberg model, where numerical simulations confirm rolls, squares, rhombi, rectangles, hexagons, and mixed-mode patterns [1501.04890].

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 [1311.1748]. 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 \((u_1(x=0,t),u_2(x=0,t))\)-space and \((E,\dot E)\)-space [2412.17076].

Cross-diffusion also alters front propagation. A minimal two-species model without self-diffusion,
\[
u_t=f(u)-v+D_vv_{xx},\qquad v_t=-D_u u_{xx},
\]
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
\[
c=6\kappa^3D_uD_v.
\]
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 [2004.07962].

## 5. Entropy structures, solvability theories, and segregation

A central analytical theme is the search for entropy or Lyapunov structures. For systems
\[
\partial_t u_i-\Delta\big[a_i(U)u_i\big]=r_i(U)u_i,
\]
a formal entropy identity has the form
\[
\frac{d}{dt}\int_\Omega \Phi(U)
+\sum_{j=1}^d\int_\Omega 
\left\langle \partial_jU,\;D^2\Phi(U)D(A)(U)\partial_jU\right\rangle
=\text{reaction terms},
\]
and positivity of \(D^2\Phi(U)D(A)(U)\) is the key coercivity mechanism. This viewpoint is strong enough to treat mixed convex/concave cross-diffusion, including the case \(\gamma_2>1\), \(0<\gamma_1<1/\gamma_2\), where one cross-diffusion is convex and the other concave [1410.7377].

For triangular population-dynamics systems, entropy and duality methods yield global weak solutions. In the class
\[
\partial_t u-\Delta_x\big[(d_u+d_u u^a+d_v v^c)u\big]
=u(r_u-r_au^a-r_bv^b),
\]
\[
\partial_t v-\Delta_x\big[(d_v+d_v v^c)v\big]
=v(r_v-r_cv^c-r_du^d),
\]
the global existence theory covers low-regularity nonnegative initial data and includes the triangular SKT model as the case \(a=b=c=d=1\) [1503.07468]. A more regular triangular theory proves that the unique local smooth solution given by Amann theorem is actually global and remains in
\[
C^{\frac{2+\nu}{2},2+\nu}([0,\infty)\times\overline\Omega)\cap C^\infty((0,\infty)\times\overline\Omega)
\]
under explicit exponent restrictions [2202.10256].

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 \(L^1\) control on the reactions [1711.01463].

Not all solvability results rely on full entropy. The glioblastoma model with densities \(p\) and \(m\),
\[
\partial_t p=D_\alpha(1-\rho)\Delta p+R_p(p,m),\qquad
\partial_t m=D_\nu\big((1-\rho)\Delta m+m\Delta\rho\big)+R_m(p,m),
\quad \rho=p+m,
\]
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 \(H^1\)-type estimates that control degeneracy at \(\rho=1\) [1710.03970].

A distinct line of work concerns segregation. In the one-dimensional system
\[
\partial_t \rho=\partial_x\!\left(\rho\,\partial_x(\rho+\eta)\right)+\rho F_1(\rho,\eta)+\eta G_1(\rho,\eta),
\]
\[
\partial_t \eta=\partial_x\!\left(\eta\,\partial_x(\rho+\eta)\right)+\eta F_2(\rho,\eta)+\rho G_2(\rho,\eta),
\]
a variational splitting scheme combining reaction ODEs and optimal-transport JKO steps yields weak solutions and preserves segregation for initially segregated data when \(G_1\equiv G_2\equiv0\), even in the presence of vacuum [1711.05434].

## 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_{ij}=d_{ji}\) and \(\beta_{ij}=\beta_{ji}=2\), the free energy
\[
E(\phi_1,\dots,\phi_M)=
\int_\Omega\left(
\sum_{i=1}^M d_{ii}\frac{\phi_i^{\alpha_i+2}}{\alpha_i+2}
+d_i\frac{\phi_i^2}{2}
+\sum_{j>i}d_{ij}\frac{\phi_i^2\phi_j^2}{2}
\right)
\]
satisfies
\[
\frac{d}{dt}E=
-\sum_{i=1}^M\|\nabla\mu_i\|_{L^2(\Omega)}^2+\sum_{i=1}^M(R_i,\mu_i).
\]
A discrete counterpart, combined with Heun reaction steps and Strang splitting,
\[
\phi^{n+1}=A_{\Delta t/2}B_{\Delta t}A_{\Delta t/2}(\phi^n),
\]
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 [2001.10069].

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 \(O(N^2)\), with only six \(N\times1\) 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 [2412.17076].

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 [1211.4412] [1503.07468] [1311.1748]. 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 [2409.06860]. Biomedical applications include glioblastoma growth with size exclusion [1710.03970], localized patterns in heterogeneous cross-diffusion systems motivated by bacterial chemotaxis and ecological interactions [2210.10155], annular cross-diffusion patterns interpreted through hypoxic-tumour cross sections [2412.20097], and a pandemic model in which the death density obeys
\[
v_t=d_2\Delta u+k(t)u,
\]
so that the spatial redistribution of deaths is driven by the diffusion of the infected population rather than by \(v\) itself [2012.13452].

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.

Source: https://www.emergentmind.com/topics/cross-diffusion-reaction-diffusion-system