In this work, we consider the following type of reaction diffusion systems
where is the concentration of -th species, are diffusion coefficients, and are nonlinear reaction terms for the chemical reaction
Such a type of reaction-diffusion systems can be found in many mathematical models in chemical engineering, biology, soft matter physics and combustion theory, see [10, 29, 30, 33, 37, 46, 51, 52, 53, 58, 60] for examples.
Numerical simulation for the reaction-diffusion system (1) turns out to be very challenging, due to the stiffness brought by the reaction term. Moreover, a naive discretization to (1) may fail to preserve the positivity and the conservation property in the original system . To overcome these difficulties, many numerical methods have been developed to solve reaction kinetics and reaction-diffusion systems [5, 23, 32, 62], including some operator splitting approaches [8, 15, 25, 26, 62].
It has been discovered that for certain form of reaction-diffusion systems, in which the reaction part describes the reversible chemical reaction satisfying the law of mass action with detailed balance condition, the whole system admits an energy-dissipation law, which opens a door of developing structure-preserving numerical schemes. In more details, under certain conditions, which will be specified in the next section, the reaction-diffusion system (1) can be reformulated as a combination of two generalized gradient flows (with different patterns) for a single free energy [39, 59]. Since the reaction and diffusion parts of the original system dissipate the same free energy, it is natural to use an operator splitting approach to develop an energy stable scheme for the whole system. Based on this variational structure, a first order accurate operator splitting scheme has been constructed in a recent work , with the variational structure theoretically preserved for the numerical solution. In this approach, since the physical free energy is in the form of logarithmic functions of the concentration , a linear function of reaction trajectories , the positivity-preserving analysis of the numerical scheme at both stages has been established. Similar to the analysis in a recent article 
for the Flory-Huggins Cahn-Hilliard flow, an implicit treatment of the nonlinear singular logarithmic term is crucial to theoretically justify its positivity-preserving property. A more careful analysis reveals that, the convex and the singular natures of the implicit nonlinear parts prevent the numerical solutions approach the singular limiting values, so that the positivity-preserving property is available for the density variables of all the species. A detailed convergence analysis and error estimate have also been reported in a recent work. However, it is a not trivial task to develop a second order accurate operator splitting scheme based on this idea. In fact, most existing works of second order energy stable scheme for gradient flows are multi-step algorithms, based on either modified Crank-Nicolson or BDF2 temporal discretization, and a multi-step approximation to the concave terms is usually needed to ensure both the unique solvability and energy stability. On the other hand, a single step, second order approximation has to be accomplished at each stage in the operator splitting approach, while a theoretical justification of positivity-preserving and energy stability turns out to be very challenging.
In this article, we propose and analyze a second order accurate operator splitting scheme for the reaction-diffusion system with the detailed balance condition. Following the energetic variational formulation, the splitting scheme solves the reaction trajectory equation of at the reaction stage, and solves the diffusion equation for in the diffusion stage. To overcome the above-mentioned difficulties, we make use of a numerical profile created by the first order convex splitting algorithm, which is proved to be a second order accurate approximation to the physical quantity at time step , to construct a second order approximation to the mobility part. Then an application of modified Crank-Nicolson formula leads to a second order approximation to the mobility function at the intermediate time instant . Meanwhile, the physical energy does not contain any concave part in the reaction-diffusion system, so that a single step, modified Crank-Nicolson method leads to a second order accurate algorithm. In addition, an artificial second order Douglas-Dupont-type regularization term , in the form of , is added in the chemical potential, to ensure the positivity-preserving property. The energy stability is derived by a careful energy estimate, because of the choice in the modified Crank-Nicolson approximation. These techniques lead to a second order accurate, positivity preserving and energy stable algorithm in the reaction stage.
In the diffusion stage, an exact integrator, so called exponential time differencing (ETD) method is applied if the diffusion coefficients are constant. Such an ETD method solves the diffusion stage equation exactly (by keeping the finite difference spatial discretization), so that both the positivity-preserving and energy stability are ensured. If the diffusion coefficients are nonlinear, we have to apply a similar idea as in the reaction stage: a predictor-corrector approach in the mobility approximation and a modified Crank-Nicolson algorithm for the chemical potential. In either case, both the positivity-preserving and energy stability could be theoretically justified for the numerical solution in the diffusion stage. Finally, a combination of the numerical algorithms at both stages by the Strang splitting approach leads to a second-order accurate, structure preserving scheme for the original reaction-diffusion system.
The rest of this article is organized as follows. The energetic variational approach is reviewed in Section 2, for the reaction-diffusion systems with the detailed balance condition. Subsequently, the second-order operator splitting scheme is presented in Section 3. The positivity-preserving and energy stability analyses will be provided at each stage as well. Some numerical results will be presented in Section 4, to demonstrate the performance of the second order operator splitting scheme.
2 Review of the energetic variational approach for reaction-diffusion systems
In this section, we briefly review the energetic variational approach for reaction-diffusion systems with detailed balance, which will be the foundation of the second order operator splitting scheme developed in the next section. We refer interested readers to [41, 59] for more detailed descriptions.
The energetic variational approach (EnVarA) [22, 27, 40], which is inspired by the seminal works of Rayleigh  and Onsager [48, 49], provides a systematic way to derive the dynamics of the system from a prescribed energy-dissipation law. In more details, an energy-dissipation law, which comes from the first and second law of thermodynamics, can be written as
for an isothermal closed system, where is the total energy, including both the kinetic energy and the Helmholtz free energy , and is the energy dissipation rate which is equal to the entropy production in the process. The energy-dissipation law, along with the kinematics of employed variables, describe all the physics and the assumptions in the system. Starting with an energy-dissipation law, the EnVarA derives the dynamics of the systems through two variational principles, the Least Action Principle (LAP) and the Maximum Dissipation Principle (MDP). The LAP, which states the equation of motion for a Hamiltonian system can be derived from the variation of the action functional , with respect to the flow maps, gives a unique procedure to derive the conservative force for the system. In the MDP, variation of the dissipation potential , which equals to in the linear response regime, with respect to the rate (such as velocity), gives the dissipation force for the system. In turn, the force balance condition leads to the evolution equation to the system
In this formulation, the energy-dissipation law, along with the kinematics of state variables,describes all the physics and the assumptions for a given system. The energetic variational approach has been successfully applied to build up many mathematical models , including systems with chemical reactions [59, 60]; it has also provided a guideline of designing structure-preserving numerical schemes for systems with variational structures [41, 44, 45], etc.
2.1 Reaction kinetics
Consider a system with species and reversible chemical reactions given by
, the concentrations of all species. The variable vectorsatisfies the reaction kinetics
where is the reaction rate for the chemical reaction, and is the stoichiometric coefficients. From (4), it is noticed that
In turn, one can define linearly independent conserved quantities for the reaction network. In the classical chemical kinetics, is determined by the law of mass action (LMA), which states that the reaction rate is directly proportional to the product of the reactant concentrations, i.e.,
in which and are the forward and backward reaction constants for the -th reaction.
where the first part stands for the entropy, and is the internal energy associated with each species. In general, depends on and , and the choice of determines the equilibrium of the system. We assume that is a constant throughout this paper. Moreover, it has been shown that the reaction kinetics (4) along with the law of mass action (6) admits a Lyapunov function if there exists a strictly positive equilibrium point , satisfying
The condition is known as the detailed balance condition. Within , one can define the Lyapunov function as
It can be noticed that and are related through
To transform the reaction kinetics into a variational frame, it is important to introduce another state variable , known as the reaction trajectory [50, 59], or the extent of reaction [12, 36]. The -th component of corresponds to the number of -th reaction that has happened by time in the forward direction. For any initial condition , the value of can be represented in terms of as the following formula
This equation can be viewed as the kinematics of a reaction kinetics, which embodies the conservation properties (5). In particular, the positivity of requires a constraint on :
Subsequently, the reaction rate can be defined as , known as the reaction velocity . In the framework of the EnVarA, we can describe the reaction kinetics through the energy-dissipation law in terms of and :
where is the rate of energy dissipation in the chemical reaction process. Unlike mechanical systems, the rate of energy dissipation for reaction kinetics may not be quadratic in terms of , since the system is often far from equilibrium [4, 14]. For a general nonlinear energy dissipation
one can specify
which turns out to be the chemical affinity, and is the chemical potential of th species. The chemical affinity is the driving force of the chemical reaction [12, 13, 36], and the dissipation makes a connection between the reaction rate and the chemical affinity. A typical choice of is given by
One can derive the law of mass action by taking . Since near an equilibrium, we see that
In turn, the energy-dissipation law (12) becomes an -gradient flow in terms of .
The reaction kinetics can be viewed as a generalized gradient flow, with a nonlinear mobility in terms of the reaction trajectory. Hence, it is expected that the numerical techniques for -gradient flows can be applied to reaction kinetics.
2.2 Reaction-diffusion systems
One can extend the energetic variational formulation for reaction kinetics to reaction-diffusion system with detailed balance, which is the foundation of the operator splitting scheme developed in the next section. For a reaction-diffusion system with species and reactions, the concentration satisfies the kinematics
where is the average velocity of each species by its own diffusion, represents various reaction trajectories involved in the system, with being the stoichiometric matrix as defined in section 2.1. The quantities and can be obtained through an energy-dissipation law [6, 59]
which leads to a reaction-diffusion equation. Here is the free energy given by (7), and and are dissipations for the mechanical and reaction parts, respectively. One key point is that the reaction and diffusion parts of the system dissipate the same free energy. To derive the reaction diffusion equation (1), could be taken as
and could be taken as
The energetic variational approach could be applied to the reaction and diffusion parts, respectively, so that the “force balance equation” is obtained for the chemical and mechanical subsystems. Formally, a direct computation implies that
which in turn gives
In particular, a linear reaction-diffusion system can be obtained by choosing :
Other choices of can result in some porous medium type nonlinear diffusion equation 
where is the concentration-dependent diffusion coefficient.
In this formulation, the reaction part is reformulated in terms of reaction trajectories , and the reaction and diffusion parts impose different dissipation mechanisms for the same physical energy.
3 The second-order operator splitting scheme
In the section, we construct a second-order operator splitting scheme to a reaction-diffusion system based on the energetic variational formulation outlined in the last section, in which the numerical discretization for the reaction part is applied to the reaction trajectory in the reaction space, while the numerical method for the diffusion part is designed to the concentration in the species space. To illustrate the idea, we focus on a case with one reversible detailed balance reaction, given by
where and are constants. Moreover, we assume that the reaction-diffusion system satisfies the energy-dissipation law (20). Numerical schemes for systems involving multiple reversible reactions can be constructed in the same manner.
To simplify the numerical description, the reaction-diffusion equation (24) can be rewritten as
where is a reaction operator and a diffusion operator. Throughout this section, the computational domain is taken as with a periodic boundary condition, and with being the spatial mesh resolution throughout this section; a computational domain with other boundary condition or numerical mesh could be analyzed in a similar fashion. In addition, the discrete free energy is defined as follows, with the given spatial discretization:
where denotes the discrete inner product.
Following the second-order Strang splitting formula , the numerical solution can be obtained through three stages. Given with , we update via the following three stages.
Stage 1. First, we set and solve the reaction trajectory equation, subject to the initial condition , with a second-order, positivity-preserving, energy-stable scheme, with the temporal step-size . An intermediate numerical profile is updated as
Stage 2. Starting with the intermediate variable , we solve the diffusion equation by a second-order, positivity-preserving and energy-stable scheme with the temporal step-size to obtain .
Stage 3. We set and repeat the stage 1, i.e., solving the reaction trajectory equation, subject to the initial condition with the temporal step-size to obtain . The numerical solution at is updated as
More details of the numerical algorithms at each stage will be provided in the following subsections.
3.1 Second-order algorithm for reaction kinetics
We first develop a second order algorithm for the reaction stage, which only needs to be constructed in a point-wise sense. The discrete free energy can be reformulated in terms of at each mesh point, denoted by
For simplicity of presentation, we omit the grid index throughout this subsection. Following the earlier discussions, for a given initial condition , the reaction trajectory equation is given by
where is the nonlinear mobility that takes the form , with is the stoichiometric vector, and is the chemical potential associated with -species. Similar to an gradient flow, a second-order algorithm for the reaction trajectory equation (31) can be constructed through a Crank-Nicolson type discretization
where is a suitable approximation to the chemical affinity, , at , is an approximation to , which needs to be independent on . The primary difficulty is focused on the construction of and , to ensure the unique solvability, as well as the positivity of and .
First, we use a first-order scheme to obtain a rough “guess” to , denoted by , as a numerical solution to
in the admissible set. This first-order scheme was proposed in , while the unique solvability and the positivity preserving property have been proved. With at hand, we introduce . Although (33) corresponds to a first order truncation error, we see that is a second order approximation to , locally in time, due to the term in the denominator. In turn, becomes a second order approximation to . To approximate , we apply the idea of discrete variational derivative method [21, 24]. More specifically, the following function is introduced
as a second-order approximation to . In fact, it is also known as the discrete variation of .
With the combined arguments, the second-order algorithm is constructed as
The term is added for the theoretical analysis the positivity-preserving property. This term is artificial, and it will not effect the second order accuracy in the temporal discretization.
This algorithm can be reformulated as an optimization problem
where , and
is a function that measures the “distance” between and . An explicit form of is not available. On the other hand, we can prove that admits a unique minimizer in the admissible set. More precisely, the following theorem is valid.
To facilitate the proof of this result, the following smooth functions are introduced, for fixed :
By a direct calculation, it is straightforward to prove the following results, which will be used in the proof of Theorem 3.1.
For any fixed , we have: (1) for any ; (2) is convex in terms of ; (3) There exists between and such that ; (4) Since increases in terms of , we have for any .
Now we can proceed into the proof of Theorem 3.1.
Recall the minimization problem (36), and it is clear that is a strictly convex function over . We only need to prove that the minimizer of over could not occur on the boundary of , so that a minimizer corresponds to a numerical solution of (35) in .
The following closed domain is considered in the analysis:
A careful calculation indicates that, for any , the following bounds are satisfied
i.e., or . Since is a bounded, compact set, there exists a (may not unique) minimizer of over . Moreover, we have to prove that, such a minimizer could not occur on the boundary points in , if is sufficiently small, by using the singular property of logarithmic function approaches to .
Without loss of generality, the minimization point is assumed to be . A direct calculation gives
Next we show that is bounded, so that we can choose sufficiently small with
which leads to a contradiction since there will be such that
To derive a bound for , we notice that
where and is a constant. Since is an increasing function of for any , the following inequality is valid:
in which is sufficiently small such that . Similarly, we have
with sufficiently small such that . Hence,
Following the same argument, the following inequality could be derived:
where is a constant. So we can choose small enough such that , which leads to the contradiction inequality (43).
Using similar arguments, if , we can prove that
Then can be chosen to be sufficiently small such that , which leads to a contradiction. Meanwhile, if , we will have .
As a result, the global minimum of over could only possibly occur at an interior point, if is sufficiently small. In turn, there is a minimizer , in the interior region of , of , so that . In other words, has to be the numerical solution of (35), provided that is sufficiently small. Therefore, the existence of a “positive” numerical solution is proved. In addition, since is a strictly convex function over , the uniqueness of this numerical solution follows from a standard convexity analysis. The proof of Theorem 3.1 is finished.
The energy stability of the numerical scheme (35) is stated below.
For a given , the numerical solution to (35) satisfies the energy-dissipation estimate
Multiplying both side of (35) by and rearranging terms yields
In the derivation of the above inequality, the following fact has been used:
which comes from the monotonic property of the logarithmic function.
Without the additional term , the discrete energy dissipation law (52) is an exact time discretization to the continuous energy-dissipation law, which is the advantage of the discrete variational derivative method. It is crucial to add this term to establish the positivity-preserving property