Field–Noyes Reaction–Diffusion Model

• Field–Noyes reaction–diffusion model is a semilinear parabolic system that describes chemical concentrations with polynomial reaction terms and homogeneous Neumann conditions.
• The model uses analytic semigroup theory integrated with Galerkin finite elements and exponential Runge–Kutta time discretization to manage nonsmooth initial data.
• It bridges classical reaction–diffusion analysis with energetic variational frameworks, providing sharp error bounds and insights into temporal order reduction near t=0.

The Field–Noyes reaction–diffusion model is the initial-boundary value problem studied as a simplified macroscopic description of the Belousov–Zhabotinsky reaction, in which three unknowns $u_1,u_2,u_3$ represent chemical concentrations on a convex polygonal domain $\Omega\subset\mathbb{R}^d$, $d=1,2,3$, with homogeneous Neumann boundary conditions and positive parameters $a_1,a_2,a_3,\lambda,\delta,\rho,c$ (Zhang et al., 21 Jul 2025). In the formulation analyzed in "Exponential Runge-Kutta Galerkin finite element method for a reaction-diffusion system with nonsmooth initial data" (Zhang et al., 21 Jul 2025), the model is treated as a semilinear parabolic system with polynomial reaction terms and diffusion, while "Field Theory of Reaction-Diffusion: Mass Action with an Energetic Variational Approach" (Wang et al., 2020) provides a broader reaction–diffusion framework that is not a Field–Noyes paper but is explicitly described as relevant to the Belousov–Zhabotinsky / Oregonator / Field–Noyes family.

1. Definition of the model

The Field–Noyes reaction–diffusion system considered in (Zhang et al., 21 Jul 2025) is

$$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$

The three components are chemical concentrations, and the paper does not introduce a further nondimensionalization beyond this form (Zhang et al., 21 Jul 2025). The reaction terms are polynomial in the components and include the quadratic couplings $u_1u_3$ and $u_1^2$. This polynomial structure is central to the later low-regularity analysis because the nonlinear bounds are derived directly for the specific Field–Noyes nonlinearity.

A semigroup formulation is obtained by defining

$A=\mathrm{diag}\{A_1,A_2,A_3\},$

where $A_1,A_2,A_3$ are the realizations in $L^2$ with homogeneous Neumann boundary conditions of

$\Omega\subset\mathbb{R}^d$0

and by writing the nonlinear term as

$\Omega\subset\mathbb{R}^d$1

With $\Omega\subset\mathbb{R}^d$2, the PDE becomes

$\Omega\subset\mathbb{R}^d$3

and the mild solution is expressed through

$\Omega\subset\mathbb{R}^d$4

This abstraction is not merely notational. It is the basis for the analytic semigroup estimates and fractional-domain arguments used throughout the sharp error theory.

2. Fractional regularity and nonsmooth initial data

A defining feature of the analysis in (Zhang et al., 21 Jul 2025) is the assumption that the initial data satisfy only

$\Omega\subset\mathbb{R}^d$5

More precisely,

$\Omega\subset\mathbb{R}^d$6

and

$\Omega\subset\mathbb{R}^d$7

Accordingly, the regime $\Omega\subset\mathbb{R}^d$8 is treated as one of nonsmooth initial data, because classical second-order time-integrator analyses typically assume much stronger regularity, such as $\Omega\subset\mathbb{R}^d$9, $d=1,2,3$0, or sufficient smoothness to guarantee bounded $d=1,2,3$1, $d=1,2,3$2, and bounded derivatives of the nonlinearity along the solution (Zhang et al., 21 Jul 2025).

The regularity bounds stated in the paper are

$d=1,2,3$3

and

$d=1,2,3$4

These estimates show explicitly that if $d=1,2,3$5, derivatives can blow up as $d=1,2,3$6. In particular, for low $d=1,2,3$7, $d=1,2,3$8, $d=1,2,3$9, and derivatives of $a_1,a_2,a_3,\lambda,\delta,\rho,c$0 are singular at the initial time. This is the mechanism behind temporal order reduction: the paper states that the temporal order is not always $a_1,a_2,a_3,\lambda,\delta,\rho,c$1, but rather $a_1,a_2,a_3,\lambda,\delta,\rho,c$2 (Zhang et al., 21 Jul 2025).

The functional-analytic setting is built in

$a_1,a_2,a_3,\lambda,\delta,\rho,c$3

with $a_1,a_2,a_3,\lambda,\delta,\rho,c$4 self-adjoint, negative definite, and generating the analytic semigroup $a_1,a_2,a_3,\lambda,\delta,\rho,c$5. The norm equivalence

$a_1,a_2,a_3,\lambda,\delta,\rho,c$6

and the semigroup estimates

$a_1,a_2,a_3,\lambda,\delta,\rho,c$7

are then used to turn low initial regularity into explicit time-singular factors. This suggests that the model’s numerical analysis is inseparable from the parabolic smoothing properties of the diffusion operator.

3. Galerkin finite elements and exponential Runge–Kutta discretization

The spatial approximation in (Zhang et al., 21 Jul 2025) uses a linear Galerkin finite element method. For a regular triangulation $a_1,a_2,a_3,\lambda,\delta,\rho,c$8 of $a_1,a_2,a_3,\lambda,\delta,\rho,c$9 with maximal mesh size $$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$0,

$$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$1

The $$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$2-projection $$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$3 is defined by

$$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$4

The discrete operator $$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$5 is defined through

$$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$6

and the Ritz operator $$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$7 satisfies

$$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$8

The paper notes that the bilinear form is written as

$$\left\{ \begin{aligned} &\frac{\partial u_1}{\partial t}=a_1\Delta u_1+\lambda^{-1}(\rho u_3-u_1u_3+u_1-u_1^2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_2}{\partial t}=a_2\Delta u_2+u_1-u_2, &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_3}{\partial t}=a_3\Delta u_3+\delta^{-1}(-\rho u_3-u_1u_3+cu_2), &&\text{in } \Omega\times(0,T],\ &\frac{\partial u_1}{\partial n(x)}=\frac{\partial u_2}{\partial n(x)}=\frac{\partial u_3}{\partial n(x)}=0, &&\text{on } \partial\Omega\times(0,T],\ &u_1(x,0)=u_{1,0}(x),\quad u_2(x,0)=u_{2,0}(x),\quad u_3(x,0)=u_{3,0}(x), &&\text{in } \Omega . \end{aligned} \right.$$9

but also states that this appears to be a typographical inconsistency; the intended form is the diagonal sum corresponding componentwise to the three equations (Zhang et al., 21 Jul 2025). That observation matters because the rest of the theory treats the operator as block-diagonal.

The semidiscrete problem is

$u_1u_3$0

with

$u_1u_3$1

The corresponding discrete semigroup is $u_1u_3$2, and a crucial discrete norm equivalence is

$u_1u_3$3

together with inverse and projection estimates such as

$u_1u_3$4

Time integration is performed by a second-order exponential Runge–Kutta method. Given $u_1u_3$5, the inner stage is

$u_1u_3$6

The update is

$u_1u_3$7

where

$u_1u_3$8

The paper emphasizes that the method consists of linear Galerkin finite elements in space and second-order exponential Runge–Kutta in time (Zhang et al., 21 Jul 2025). The stiff linear part is handled exactly through $u_1u_3$9, $u_1^2$0, and $u_1^2$1, while the nonlinear part is treated through stage values and linear interpolation over each time interval. A plausible implication is that this architecture is particularly well matched to diffusion-dominated semilinear systems with polynomial kinetics.

4. Sharp error analysis for low-regularity data

The central result of (Zhang et al., 21 Jul 2025) is a fully discrete error theory that remains valid when the initial data belong only to $u_1^2$2, $u_1^2$3. The principal estimates are

$u_1^2$4

and

$u_1^2$5

These bounds combine spatial finite element error and temporal exponential Runge–Kutta error. Their distinctive feature is the presence of the singular prefactors in $u_1^2$6, which quantify the near-initial-time loss caused by nonsmooth data.

For the semidiscrete approximation, the solution satisfies

$u_1^2$7

and the spatial error estimate is

$u_1^2$8

Thus the $u_1^2$9 spatial error is $A=\mathrm{diag}\{A_1,A_2,A_3\},$0, and the $A=\mathrm{diag}\{A_1,A_2,A_3\},$1 spatial error is $A=\mathrm{diag}\{A_1,A_2,A_3\},$2.

For temporal error, the paper first proves a general $A=\mathrm{diag}\{A_1,A_2,A_3\},$3-estimate,

$A=\mathrm{diag}\{A_1,A_2,A_3\},$4

with

$A=\mathrm{diag}\{A_1,A_2,A_3\},$5

Optimizing over $A=\mathrm{diag}\{A_1,A_2,A_3\},$6 yields the sharp estimates

$A=\mathrm{diag}\{A_1,A_2,A_3\},$7

and

$A=\mathrm{diag}\{A_1,A_2,A_3\},$8

The key message is explicit in the paper: if $A=\mathrm{diag}\{A_1,A_2,A_3\},$9, time convergence is reduced to approximately $A_1,A_2,A_3$0; if $A_1,A_2,A_3$1, the method recovers its full second order (Zhang et al., 21 Jul 2025). The analysis therefore differs fundamentally from smooth-data proofs that assume bounded high derivatives and conclude uniform $A_1,A_2,A_3$2 local defects. Here the local defect is singular near $A_1,A_2,A_3$3, and the proof tracks that singularity term by term through the error recursion.

5. Nonlinear estimates, semigroup methods, and proof architecture

The proof strategy in (Zhang et al., 21 Jul 2025) combines analytic semigroup smoothing, fractional Sobolev product estimates, finite element semigroup approximation, and induction-based stability arguments. A central device is the admissible index set

$A_1,A_2,A_3$4

for which the product estimate

$A_1,A_2,A_3$5

holds. This estimate is transferred to the discrete setting and leads to bounds of the form

$A_1,A_2,A_3$6

as well as corresponding estimates for $A_1,A_2,A_3$7 and $A_1,A_2,A_3$8.

These estimates are the core nonlinear tools for low-regularity error analysis. They are needed because the nonlinearity interpolation defect and the inner-stage contribution must be controlled even when time derivatives of the exact solution are singular. The temporal local defect is analyzed by decomposing the error into three terms,

$A_1,A_2,A_3$9

where $L^2$0 propagates previous-step outer-stage error, $L^2$1 represents the inner-stage error contribution, and $L^2$2 is the interpolation defect of the exact nonlinearity (Zhang et al., 21 Jul 2025).

For the interpolation defect, the paper uses Taylor expansion,

$L^2$3

and derives weighted bounds such as

$L^2$4

and

$L^2$5

The novelty identified by the paper is not the choice of exponential Runge–Kutta itself, but the fact that the analysis allows $L^2$6, $L^2$7, quantifies temporal order reduction continuously in $L^2$8, handles the $L^2$9D case, and derives sharp $\Omega\subset\mathbb{R}^d$00 and $\Omega\subset\mathbb{R}^d$01 bounds for the Field–Noyes nonlinearity specifically (Zhang et al., 21 Jul 2025).

6. Relation to energetic variational reaction–diffusion theory, experiments, and limitations

The direct numerical-analysis paper (Zhang et al., 21 Jul 2025) does not formulate the Field–Noyes system through a thermodynamic variational principle. However, (Wang et al., 2020) provides a general reaction–diffusion framework with mass-action-type kinetics that the paper explicitly describes as relevant to the Belousov–Zhabotinsky / Oregonator / Field–Noyes family. In that framework, the starting point is an energy-dissipation law,

$\Omega\subset\mathbb{R}^d$02

with equilibrium determined by the free energy and kinetics determined by the dissipation functional. For reaction–diffusion systems, the template is

$\Omega\subset\mathbb{R}^d$03

The same paper states that this is not a Field–Noyes/Oregonator paper, but that it provides a thermodynamically grounded variational framework that can be specialized to models of the Field–Noyes class (Wang et al., 2020). This suggests a conceptual bridge: the polynomial reaction–diffusion equations used in the direct Field–Noyes analysis can be viewed within a broader distinction between equilibrium structure and nonequilibrium kinetics.

The numerical experiments in (Zhang et al., 21 Jul 2025) are carried out for the 2D problem on

$\Omega\subset\mathbb{R}^d$04

with

$\Omega\subset\mathbb{R}^d$05

and with the same initial profile for all three components,

$\Omega\subset\mathbb{R}^d$06

Four nonsmooth initial data are tested: $\Omega\subset\mathbb{R}^d$07

$\Omega\subset\mathbb{R}^d$08

$\Omega\subset\mathbb{R}^d$09

$\Omega\subset\mathbb{R}^d$10

A reference solution is computed using

$\Omega\subset\mathbb{R}^d$11

The observed spatial rates are approximately order $\Omega\subset\mathbb{R}^d$12 in $\Omega\subset\mathbb{R}^d$13 and approximately order $\Omega\subset\mathbb{R}^d$14 in $\Omega\subset\mathbb{R}^d$15 for all four initial data. The measured temporal rates at $\Omega\subset\mathbb{R}^d$16 are about $\Omega\subset\mathbb{R}^d$17 for case (1), about $\Omega\subset\mathbb{R}^d$18 for case (2), and about $\Omega\subset\mathbb{R}^d$19 for cases (3) and (4), in both $\Omega\subset\mathbb{R}^d$20 and $\Omega\subset\mathbb{R}^d$21. The paper also examines errors at the first time step $\Omega\subset\mathbb{R}^d$22, where the singular behavior is strongest, and reports that the numerical data largely confirm the theory, especially for initial data (2)–(4). For the sign-function data, the observed rates are slightly better than the generic estimate; the paper explains this by a more favorable multiplication structure, since products of sign functions retain special structure and give sharper nonlinear estimates than the generic fractional-Sobolev product bounds (Zhang et al., 21 Jul 2025).

Several limitations are stated explicitly. The error bounds are weighted by powers of $\Omega\subset\mathbb{R}^d$23, reflecting singular behavior near $\Omega\subset\mathbb{R}^d$24; a small $\Omega\subset\mathbb{R}^d$25 is unavoidable in several rates; the proof depends on delicate fractional multiplication estimates; no positivity-preserving theorem is proved; no invariant-region preservation result is established; and implementation details for evaluating matrix exponentials or $\Omega\subset\mathbb{R}^d$26-functions are not developed (Zhang et al., 21 Jul 2025). The paper also states that it does not prove discrete positivity for concentrations and does not focus on structure preservation. In parallel, (Wang et al., 2020) notes that its own numerical example does not present pattern formation, spiral waves, fronts, bistability, or oscillations, and is therefore more formal and constitutive than phenomenological for Field–Noyes-type wave phenomena. Taken together, these statements delimit the current scope: rigorous low-regularity convergence for a specific Field–Noyes reaction–diffusion system on one side, and a general energetic scaffold for mass-action reaction–diffusion systems on the other.