We develop numerical methods for stochastic reaction-diffusion systems based on approaches used for fluctuating hydrodynamics (FHD). For hydrodynamic systems, the FHD formulation is formally described by stochastic partial differential equations (SPDEs). In the reaction-diffusion systems we consider, our model becomes similar to the reaction-diffusion master equation (RDME) description when our SPDEs are spatially discretized and reactions are modeled as a source term having Poisson fluctuations. However, unlike the RDME, which becomes prohibitively expensive for an increasing number of molecules, our FHD-based description naturally extends from the regime where fluctuations are strong, i.e., each mesoscopic cell has few (reactive) molecules, to regimes with moderate or weak fluctuations, and ultimately to the deterministic limit. By treating diffusion implicitly, we avoid the severe restriction on time step size that limits all methods based on explicit treatments of diffusion and construct numerical methods that are more efficient than RDME methods, without compromising accuracy. Guided by an analysis of the accuracy of the distribution of steady-state fluctuations for the linearized reaction-diffusion model, we construct several two-stage (predictor-corrector) schemes, where diffusion is treated using a stochastic Crank–Nicolson method, and reactions are handled by the stochastic simulation algorithm of Gillespie or a weakly second-order tau leaping method. We find that an implicit midpoint tau leaping scheme attains second-order weak accuracy in the linearized setting and gives an accurate and stable structure factor for a time step size of an order of magnitude larger than the hopping time scale of diffusing molecules. We study the numerical accuracy of our methods for the Schlögl reaction-diffusion model both in and out of thermodynamic equilibrium. We demonstrate and quantify the importance of thermodynamic fluctuations to the formation of a two-dimensional Turing-like pattern and examine the effect of fluctuations on three-dimensional chemical front propagation. By comparing stochastic simulations to deterministic reaction-diffusion simulations, we show that fluctuations accelerate pattern formation in spatially homogeneous systems and lead to a qualitatively different disordered pattern behind a traveling wave.

## I. INTRODUCTION

While deterministic reaction-diffusion models have been successfully applied to explain various spatiotemporal phenomena such as pattern formation, and to gain insight into nonequilibrium transitions, it is now widely appreciated that spatiotemporal fluctuations in the concentration of chemical species play an essential role. Such internal or thermodynamic fluctuations, which arise from both reaction and diffusion processes, have molecular origin; microscopically, those processes occur through the movement and collision of individual molecules under thermal fluctuations. Hence, the deterministic macroscopic description eventually fails at smaller scales where the fluctuations are significant, and a stochastic mesoscopic description is needed. Examples include fluctuation-induced instabilities,^{1} reversal of direction of front propagation,^{2} violation of the law of mass action,^{3} long-time tails in kinetics,^{4} emergence of new steady states^{5} and patterns,^{6} acceleration of pattern formation,^{7} enhanced induction time for ignition,^{8} and the onset of homogeneous oscillations.^{9} Due to a small number of proteins involved in cellular functions,^{10} processes in cell biology are good examples^{11–14} where the stochastic reaction-diffusion description provides an indispensable modeling tool.^{15,16}

A microscopic picture of reaction-diffusion, dating back to the work of Smoluchowski,^{17} assumes that molecules undergo independent Brownian motions and reactions can occur only when two molecules are close to each other. Based on this picture, the particle-based approach to simulate a reaction-diffusion system tracks the trajectories of diffusing molecules and uses the intermolecular distance to determine whether a reaction occurs. Exact sampling of the Smoluchowski model can be performed by first-passage kinetic Monte Carlo type algorithms;^{18–20} approximate reactive Brownian dynamics (BD) using a fixed time step size forms another class of algorithms.^{21,22} While molecular schemes, such as molecular dynamics (MD) and direct simulation Monte Carlo (DSMC), can be used for reaction-diffusion problems,^{23,24} they are computationally even more expensive. Hybrid methods combining particle and coarse-grained descriptions, either using operator splitting^{25} or domain decomposition,^{26,27} have also been proposed.

For the mesoscopic description of a reactive system, the master equation approach is commonly used. For a well-mixed (i.e., spatially homogeneous) system, the time evolution of the system (i.e., the number of molecules of each chemical species) is described by the chemical master equation (CME). Exact sampling of the CME can be performed by the stochastic simulation algorithm (SSA) of Gillespie,^{28} whereas the tau leaping method^{29} can be employed as an approximate algorithm with a given time step size. Several variants of these methods have been proposed.^{30,31} For a spatially inhomogeneous system, the time evolution of the system is commonly described by the reaction-diffusion master equation (RDME), which is also known as the multivariate master equation.^{32,33} In this approach, the system is divided into homogeneous subsystems or cells and the number of molecules of each chemical species in each cell is tracked. Changes in the molecule numbers occur either through hopping events of a molecule between adjacent cells or through chemical reactions within a cell. Hopping events correspond to diffusive transport and are treated as first-order reactions. Since the RDME is a spatial extension of the CME, exact sampling of the RDME can be performed by SSA-type algorithms,^{34,35} which are called inhomogeneous SSA (ISSA).

While conceptually simple and still widely used,^{6,7,36} the traditional approach of solving the RDME by ISSA has the computational issue that the method becomes prohibitively slow as the number of molecules or cells increases. Since the cell volume should be chosen sufficiently small to ensure homogeneity over each cell, large, finely resolved grids are required for two- or three-dimensional problems. As the spatial resolution increases, the time interval between successive events becomes very short due to rapid diffusive transfer, and hopping events greatly outnumber reaction events, which slows down ISSA.^{30} Several approaches have been proposed to improve the performance of stochastic sampling of the RDME, such as the next subvolume method^{35} and its parallel simulation version.^{37} Various implementations of the tau leaping method in a spatial context,^{38–40} and the time-dependent propensity for diffusion method,^{41} have also been proposed. A more aggressive approach to reduce the computational cost by avoiding the sampling of the individual diffusion events is to split diffusion and reaction in each time step and to treat diffusion in a more efficient manner. Various sampling methods for diffusion have been proposed, including the Gillespie multi-particle method,^{42} the multinomial simulation algorithm,^{43} the adaptive hybrid method on unstructured meshes,^{44,45} and the diffusive finite-state projection algorithm.^{46,47}

In this paper, we propose a numerical algorithm for stochastic reaction-diffusion systems based on approaches used for fluctuating hydrodynamics (FHD). To incorporate the effects of thermal fluctuations in a fluid, in FHD one assumes that the dynamics of the fluid can be described by the usual hydrodynamic equations (e.g., the Navier–Stokes equations), augmenting each dissipative flux with a stochastic flux.^{48} Those stochastic fluxes are modeled by spatiotemporal Gaussian white noise (GWN) and the resulting governing equations are written as stochastic partial differential equations (SPDEs). FHD was originally developed for equilibrium fluctuations by Landau and Lifshitz^{49} and its validity has been justified for nonequilibrium systems^{50} through the theory of coarse graining.^{51} For further discussion of FHD compared to MD, see Ref. 52. Various extensions and generalizations of FHD theory have been developed and successfully applied to fluctuation-induced phenomena; see Ref. 53 and references therein. Recent work by the authors has focused on FHD models of hydrodynamic transport^{54,55} in binary fluid mixtures,^{56,57} multiphase flows,^{58} multispecies fluid mixtures,^{59,60} multispecies reactive mixtures,^{61} and electrolytes.^{62}

Compared to our previous work,^{61} where the coupling effects of fluid hydrodynamic and chemical fluctuations have been investigated, here we focus on reaction and diffusion and neglect all other hydrodynamic processes (advection, viscous dissipation, thermal conduction, and cross term effects). Rather than using a Langevin description (i.e., based on Gaussian fluctuations) of chemistry, which is only valid in the limit of vanishing fluctuations,^{61} here we employ a more accurate description of reactions based on Poisson fluctuations. As pointed out in Ref. 51, even though a formal SPDE description is employed, the actual interpretation of FHD always requires the notion of a *coarse-graining* over a certain length scale. The FHD equations are discretized using a finite volume approach^{63,64} that represents the solution in terms of the average over cells, which provides an effective coarse-graining. Therefore, reactions can be treated in a similar manner to the RDME approach when the SPDEs are spatially discretized and integrated in time using SSA or a weakly second-order tau leaping method.^{65,66} Recent relevant work by others includes Ref. 67 in which the FHD approach has been applied to reaction-diffusion systems. However, only fluctuations arising from diffusion have been considered (i.e., no fluctuations from chemical reactions) and modeled as additive noise. The FHD approach has been also applied to concentration fluctuations in a ternary liquid mixture in equilibrium^{68} and the Model H equations for binary mixtures.^{69}

The key difference between the FHD and RDME descriptions lies in the more efficient treatment of fast diffusion. A number of approximate numerical methods for the RDME^{42–47} are based on operator splitting using first-order Lie or second-order Strang splitting.^{70} In Appendix A we review and discuss in more detail a split scheme that uses multinomial diffusion sampling^{71} for diffusion and SSA for reactions. These RDME-based schemes use a time step size $\Delta t$ comparable to the hopping time scale $\tau h=\Delta x2/(2dD)$ with *d* being the spatial dimension, $\Delta x$ being the grid spacing, and *D* being a typical diffusion coefficient. Even though $\tau h$ is much larger than the mean duration between successive events in ISSA, using $\Delta t$ comparable to $\tau h$ is still very restrictive for large *D* or small $\Delta x$. In our FHD formulation, we treat diffusion implicitly using backward Euler or Crank–Nicolson, so that the time step size can be significantly larger (e.g., an order of magnitude larger for a given accuracy tolerance) than the hopping time scale. Since the time steps used in RDME simulations are already (usually an order of magnitude or more) larger than those in BD simulations, our approach allows even larger time step size compared to particle-based methods.

While the development of numerical schemes for stochastic reaction-diffusion systems described by spatiotemporal GWN dates back to the 1990s,^{72,73} much of the prior work has not been guided by numerical analysis or extensive experience from deterministic computational fluid dynamics (CFD). With the help of well-established techniques for numerical solution of PDEs and SPDEs, we construct numerical schemes in a systematic manner to ensure that accuracy is maintained for a large time step size. To this end, we employ two-stage (i.e., predictor-corrector) Runge–Kutta temporal integrators.^{64,74} Rather than using operator splitting, we treat reaction and diffusion together in each stage in a manner that is second-order weakly accurate for general linearized FHD equations. The construction of these schemes is guided by a stochastic accuracy analysis of the (static) structure factor for linearized FHD.^{63,64} The structure factor is the steady-state spectrum of the concentration fluctuations, i.e., the covariance matrix in Fourier space, see Eq. (8). We apply the techniques in Refs. 63 and 64 to predict the discrete structure factors for our scheme and compare them to analytical predictions of the continuum structure factors for our model in the linearized setting.

The FHD approach inherently outperforms the RDME approach as the number of molecules per cell increases in exactly the same way that multinomial diffusion outperforms diffusion by hopping, or tau leaping outperforms SSA. In fact, the computational cost of FHD methods does not significantly change as the magnitude of the fluctuations changes. This is an obvious advantage of the FHD approach since the macroscopic limit cannot be efficiently simulated by the RDME approach. However, the validity of the FHD approach cannot be taken for granted when there are only a small number of molecules in each cell, since in FHD the number of molecules in each cell is a continuous real-valued variable, rather than a discrete non-negative integer variable as in the RDME. We investigate this issue carefully and propose techniques to improve the accuracy of the FHD description for the case of a small number of molecules per cell, making our numerical schemes robust even for large fluctuations. In particular, we develop a spatial discretization that significantly mitigates non-negativity of the species number densities and closely reproduces the Poisson thermodynamic equilibrium distribution for the number of molecules in a cell. For numerical examples considered in this paper, we show that the mean number of molecules in a cell can be as low as 10. However, if one is specifically interested in systems with only a small number of molecules per cell, one should use an integer-based description like RDME. Moreover, in very dilute cases, a particle-based description like BD is actually fastest since most cells will have essentially no molecules in them. However, for practical stochastic simulation of reaction-diffusion systems, where the populations of chemical species may have different orders of magnitude, this kind of robustness is required; even if there are a large number of molecules in a cell, some species may have a small number of molecules.

The rest of the paper is organized as follows. Section II presents the background for our approach, including the FHD description of reaction-diffusion systems and linearized analysis in a Gaussian approximation. Section III explains how the FHD reaction-diffusion equations can be spatially discretized using a finite-volume approach. Section IV presents temporal integrators for the spatially discretized equations that handle diffusion using existing FHD techniques and treat reactions using SSA or second-order tau leaping. Section V presents simulation results of several reaction-diffusion systems. In Section V A, for testing and validation of our numerical schemes, we use a one-species Schlögl model.^{75,76} In Section V B, to compare our methods to each other and to RDME-based methods, we study two-dimensional Turing-like pattern formation in the three-species Baras–Pearson–Mansour (BPM) model.^{77,78} In Section V C, to demonstrate the ability of our approach to scale to larger systems, we present numerical simulation results for three-dimensional front propagation in a two-species model.^{7} In Section VI, we offer some concluding remarks and suggest future research directions.

## II. BACKGROUND

In Section II A, we present the continuous-time continuous-space FHD description of reaction-diffusion systems. Here, we assume that fluctuations in chemistry are described by GWN (i.e., Langevin type). A more accurate description of chemistry based on Poisson fluctuations is incorporated in the continuous-time discrete-space description in Section III. In Sections II B and II C, we introduce the structure factor and the Schlögl reaction-diffusion model, respectively. As one of the criteria for the development and analysis of numerical schemes, later in the paper, we investigate how accurately a numerical scheme produces the structure factor for the Schlögl model.

In this section, we introduce several GWN vector and scalar random fields and denote them by $\U0001d4e9(\bm{x},t)=(Z1(\bm{x},t),\u2026,Zd(\bm{x},t))$ and $Z(\bm{x},t)$, respectively, with additional superscripts to distinguish the different fields. We assume that any two distinct processes are independent and that the noise intensity of each process is normalized: $\u27e8Zj(\bm{x},t)Zj\u2032(\bm{x}\u2032,t\u2032)\u27e9=\delta jj\u2032\delta (\bm{x}\u2212\bm{x}\u2032)\delta (t\u2212t\u2032)$ and $\u27e8Z(\bm{x},t)Z(\bm{x}\u2032,t\u2032)\u27e9$ = $\delta (\bm{x}\u2212\bm{x}\u2032)\delta (t\u2212t\u2032)$. Similarly, we denote GWN vector and scalar random processes by $W(t)$ and $W(t)$, respectively, and assume $\u27e8Wj(t)Wj\u2032(t\u2032)\u27e9=\delta jj\u2032\delta (t\u2212t\u2032)$ and $\u27e8W(t)W(t\u2032)\u27e9$ $=\delta (t\u2212t\u2032)$.

### A. FHD description

We consider a reaction-diffusion system having *N*_{s} species undergoing *N*_{r} reactions in *d*-dimensional space. By denoting the number density of species *s* by $ns(\bm{x},t)$, the equations of FHD for $\bm{n}(\bm{x},t)=(n1(\bm{x},t),\u2026,nNs(\bm{x},t))$ are written formally as the SPDEs^{61}

where *D*_{s} is the diffusion coefficient of species *s*, $ar(\bm{n})$ is the propensity function indicating the rate of reaction *r*, and $\nu sr$ is the stoichiometric coefficient of species *s* in reaction *r*. In the macroscopic limit of vanishing fluctuations, Eq. (1) approaches the deterministic reaction-diffusion PDE (law of large numbers),

We explain below how the diffusion and reaction parts are obtained by considering the diffusion-only (i.e., no-reaction) and reaction-only (i.e., well-mixed) cases.

#### 1. Diffusion

The diffusion-only SPDE

can be justified by considering a microscopic system where each molecule *i* undergoes independent Brownian motion,

In this *formal* derivation,^{79} one defines the instantaneous number density field $ns(\bm{x},t)=\u2211i\delta (\bm{x}\u2212\bm{x}s,i(t))$ and uses Ito’s rule to obtain Eq. (3). This equation can also be obtained from the diffusion portion of the general multispecies FHD equations,^{60,61} given by nonequilibrium statistical mechanics,^{80} by assuming a dilute solution and considering only solute species. In addition, the linearized version of Eq. (3) can be obtained from the multivariate master equation model (i.e., diffusion by hopping) near the macroscopic limit.^{32} However, although relations to those equations reaffirm Eq. (3) near the macroscopic limit, it is important to note that Eq. (3) is formally exact even in the case where fluctuations are large, since it is simply a rewriting of Eq. (4), in a representation in which the particle numbering (identity) is lost.^{81}

We note that the FHD equations (1) and (3) are not mathematically well-defined because the solution needs to be interpreted as a distribution (or a generalized function), and the square root of a distribution is not well-defined in general. The linearized FHD equation does not suffer from such an issue and is well-defined; the problems arise due to the multiplicative noise in Eq. (3). However, even though Eq. (3) is ill-defined, it is formally consistent with the law of large numbers (given by the deterministic reaction-diffusion equation (2)), the central limit theorem (given by the linearized FHD equation (12)), and large deviation theory for a collection of Brownian walkers. In this sense, Eq. (3) is a meaningful *representation* of the physical model that is useful in constructing *well-defined* mesoscopic descriptions via spatial discretization of the formal SPDEs. Compared to obtaining a mesoscopic model by directly coarse-graining the microscopic model, the spatial discretization of the SPDE is easier in general and can be done in a systematic manner.^{82}

#### 2. Reaction

To see how the reaction part of Eq. (1) is obtained, consider a well-mixed system with volume $\Delta V$. By assuming that the time evolution of $\mathit{n}(t)$ follows the CME, we express the change over the infinitesimal time interval *dt* as follows:^{30}

where $P(m)$ denotes a Poisson random variable having mean *m*. Note that Eq. (5) is equivalent to the CME if interpreted in the Ito sense. The specific form of the chemical rate function $ar(\bm{n})$ that we use in this work is described in Section III C. Henceforth, we will formally write Eq. (5) in the differential form,

The chemical Langevin equation (CLE)^{85} is obtained under the assumption that the mean number of reaction occurrences is large.^{30} That is, the assumption enables one to replace $P(m)$ by a Gaussian random variable having the same mean and variance, to give the CLE,

Since reaction is assumed to be local, the reaction part of Eq. (1) is obtained from the spatial extension of Eq. (7).

One of the important conclusions of our previous work^{61} was that the Langevin description (7) is not consistent with equilibrium statistical mechanics. Alternative formulations based on a Langevin diffusion description^{86,87} that are consistent at thermodynamic equilibrium fail to correctly model relaxation toward equilibrium.^{61} Instead, in order to correctly capture both small fluctuations and large deviations in equilibrium and non-equilibrium contexts, one must retain a description of chemical reactions as a Markov jump process. That is, one must describe reactions using a stochastic differential equation driven by Poisson rather than Gaussian noise.

### B. Structure factor

The structure factor is the steady-state spectrum of the concentration fluctuations,

i.e., the variance of the Fourier mode of the number density of species *s*,

Here we have assumed a periodic domain of volume *V* and defined $\delta n^s,\bm{k}=n^s,\bm{k}\u2212\u27e8n^s,\bm{k}\u27e9$, where the brackets $\u27e8\u27e9$ denote the equilibrium average. Here we derive an analytic expression of the structure factor from the linearized FHD equation. We assume that there is only one species, *N*_{s} = 1, and drop the subscript *s* for species, to write Eq. (1) as

where

and we have expressed fluctuations arising from all reactions by a single GWN field $\x90(R)$. At a spatially uniform stable steady state, $n(\bm{x},t)$ fluctuates around mean number density $n\xaf\u2261\u27e8n\u27e9$, where $a(n\xaf)=0$ and $a\u2032(n\xaf)<0$. The linearization of Eq. (10) around this equilibrium state is given by the central limit theorem,

where $r=\u2212a\u2032(n\xaf)>0$ is the effective reaction rate near equilibrium and $\Gamma \xaf=\Gamma (n\xaf)$.

The Fourier transform of Eq. (12) gives

Since Eq. (13) has the form of the Ornstein–Uhlenbeck equation,^{32} the structure factor is easily obtained as

where $\u2113=D/r$ denotes the penetration depth. From Eq. (14), we observe that there are two limiting cases. In the small wave number limit $k\u2113\u226a1$, $S(\bm{k})$ becomes $\Gamma \xaf/r$ and does not depend on diffusion. In fact, the result $S(\U0001d7ce)=\Gamma \xaf/r$ is also obtained from the CME assuming the whole system is well-mixed. On the other hand, in the large wave number limit $k\u2113\u226b1$, $S(\bm{k})$ becomes $n\xaf$, which is the result for the diffusion-only system. Hence, fluctuations are reaction-dominated at a length scale larger than $\u2113$ and are diffusion-dominated at a length scale smaller than $\u2113$.

We also observe that if the system is in detailed balance at its steady state, i.e., it is in thermodynamic equilibrium, then $\Gamma \xaf=n\xafr$ and $S(\bm{k})=n\xaf$, consistent with a product Poisson distribution with mean number density $n\xaf$. Therefore, in true thermodynamic equilibrium, the statistics of the fluctuations are independent of any kinetic parameters, as they must be according to equilibrium statistical mechanics.^{88} In particular, the presence of the reactions does not change the Poisson statistics of the state of thermodynamic equilibrium. In Section V A, we use this property to judge the quality of numerical schemes.

### C. Schlögl model

The Schlögl model^{75,76} is given by the chemical reactions for species X,

Hence, we have *N*_{s}=1, *N*_{r}=4, $\nu 1=\nu 3=1$, $\nu 2=\nu 4=\u22121$, *a*(*n*)=*k*_{1}*n*^{2}−*k*_{2}*n*^{3}+*k*_{3}−*k*_{4}*n*, and $\Gamma (n)=12(k1n2+k2n3+k3+k4n)$. Due to the cubic nonlinearity of *a*(*n*), the well-mixed system exhibits several kinds of distributions depending on the values of the rate constants. If detailed balance is satisfied, that is, $k1neq2=k2neq3$ and *k*_{3}= *k*_{4}*n*_{eq}, the system is in thermodynamic equilibrium and the distribution follows Poisson statistics with mean number density *n*_{eq}. Otherwise, depending on the number of real roots of *a*(*n*)= 0, the system exhibits a monostable distribution (for a single positive root) or a bistable distribution (for three positive roots).^{76}

The structure factor of the spatially extended Schlögl model can be calculated from Eq. (14). As expected from the fact that the equilibrium distribution of the system follows Poisson statistics, $S(\bm{k})=neq$ in the case of thermodynamic equilibrium. Note, however, that having a monostable distribution does *not* imply thermodynamic equilibrium. The structure factor of the out-of-equilibrium monostable case is not flat but exhibits a transition near $k\u2113\u223c1$. For the bistable case exhibiting metastability, the linearized theory is still applicable if one looks at fluctuations around one of the two peaks. However, in this work, we focus on the equilibrium and out-of-equilibrium cases where *a*(*n*) has a single positive root.

## III. SPATIAL DISCRETIZATION

In this section, we discuss spatial discretization of the FHD equation using a finite-volume approach^{63,64} that converts the SPDE into stochastic ordinary differential equations (SODEs) for the cell number density $ns,\bm{i}(t)$. We develop numerical schemes to solve these SODEs in Section IV. In Section III A, we first discretize the diffusion-only SPDE (3). In Section III B, we add reactions and present the continuous-time discrete-space description of the reaction-diffusion system. In Section III C, we discuss techniques to handle a small number of molecules per cell.

For simplicity, in this paper we only consider periodic systems. However, our methods can be straightforwardly generalized to standard types of physical boundary conditions (Dirichlet, Neumann, or Robin). In particular, since chemistry is local and does not require boundary conditions, one can rely on methods we have developed in prior work without chemistry; see, for example, the discussion in Ref. 89.

### A. Diffusion-only case

Due to the lack of regularity of $(s(D)(\bm{x},t)$ in Eq. (3), pointwise values of $ns(\bm{x},t)$ are not physically meaningful. Hence, we consider instead the spatial average of $ns(\bm{x},t)$ over a cell. We partition the system domain $L1\xd7L2\xd7\cdots \xd7Ld$ into cells of volume $\Delta V=\Delta x1\cdots \Delta xd$ and denote the cell number density of species *s* in cell $\bm{i}=(i1,\u2026,id)$ as

We denote the face of a cell using the index ** f**. If two contiguous cells have indices

**and $\bm{i}+\bm{e}j$ (with $\bm{e}j$ being the unit vector along the**

*i**j*-axis), the face

**shared by the cells is denoted by $\bm{i}+12\bm{e}j$.**

*f*To obtain a spatial discretization of Eq. (3) that ensures discrete fluctuation-dissipation balance,^{63,64} we use the standard second-order discrete Laplacian operator for the deterministic diffusion part $D\u22072sns$ and introduce a staggered grid for the stochastic diffusive flux term $\u2207\u22c5(2Dsns(s(D))$, see Fig. 1. For *d* = 1, a formally second-order spatial discretization of Eq. (3) is written as

The spatial average of $ns(\bm{x},t)$ over the interval of length $\Delta x$ around face $i\xb112$ is approximated by $n\u223cs,i\xb112(t)$, whereas that of $\x90s(\bm{x},t)$ is modeled by $1\Delta VWs,i\xb112(t)$. To close the equation, $n\u223cs,i\xb112(t)$ is approximated by an average of *n*_{s,i}(*t*) and $ns,i\xb11(t)$, that is, $n\u223cs,i\xb112=n\u223c(ns,i,ns,i\xb11)$. Natural candidates for the averaging function $n\u223c(n1,n2)$ would be the Pythagorean means: the arithmetic, geometric, and harmonic means. We choose a modified arithmetic average for $n\u223c(n1,n2)$ described in Section III C, for reasons detailed in Appendix B.

Generalization of the spatial discretization (17) to higher dimensions is straightforward. For each face, a GWN process $Ws,\bm{f}$ is defined and $n\u223cs,\bm{f}(t)$ is calculated from the cell number densities of the two cells sharing the face by using the averaging function $n\u223c(n1,n2)$, see Fig. 1. By introducing notations $ns(t)\u2261{ns,\bm{i}(t)}$, $Ws(t)\u2261{Ws,\bm{f}(t)}$, and $n\u223cs(t)\u2261n\u223c(ns(t))$, we express the resulting SODEs for ${ns,\bm{i}(t)}$ as

with the understanding that $\u2207d2$ denotes the standard (2*d* + 1)-point discrete Laplacian operator, and $\u2207d\u22c5$ denotes a discrete divergence operator.

### B. Reaction-diffusion system

By combining Eqs. (6) and (18), we obtain the spatial discretization of the reaction-diffusion FHD equations as a system of Ito SODEs,

In Eq. (1), fluctuations in the reaction rate are modeled as GWN, while in Eq. (19), we assume Poisson fluctuations. Since the latter fluctuations are consistent with discrete nature of reactions based on the CME, the description in Eq. (19) is physically more accurate. In fact, it has been shown that the CLE description can give physically incorrect results since it is *not* consistent with a Gibbs–Boltzmann or Einstein equilibrium distribution, even for the case of a single well-mixed cell.^{61} As shown above, the inclusion of Poisson fluctuations for reaction, however, requires the notion of a mesoscopic cell and thus can be realized only after the SPDE is spatially discretized.

The choice of appropriate cell size is a delicate issue for the RDME and FHD descriptions. An upper bound on the cell size is given by the penetration depth due to the underlying assumption that each cell is homogeneous and reactions occur within a cell. In fact, there is not only an upper bound of the cell size for a valid description but also a lower bound. This can be seen by considering the fact that bimolecular reactions would become increasingly infrequent as the cell size decreases.^{90,91} Several criteria for choosing the cell size have been proposed based on physical arguments^{90,92,93} and mathematical analysis.^{94} For a small value of the cell size, corrections in the rate constants of bimolecular reactions have been proposed.^{90,95,96} However, these corrections do not fix the underlying problem which comes from the fact that reactions are treated as a purely local process with no associated spatial length scale. In microscopic (particle) models of reaction-diffusion such as the Smoluchowski^{17} model or the Doi model,^{90} a microscopic reactive distance appears and controls the reaction rate for diffusion-limited reactions. By introducing a reactive distance into the model, and relaxing the restriction that a reaction should occur among the molecules in the same cell, a modified convergent RDME having well-defined limiting behavior for small cell size can be developed^{97} and could be combined with our FHD description of diffusion. The dependence of stochastic Turing patterns on the grid size has been also investigated.^{98}

### C. Maintaining non-negative densities

The spatially discretized FHD equation (18) or (19) is well defined but suffers from two issues that we now address. First, the number of molecules in a cell (i.e., $ns,\bm{i}\Delta V$) is not an integer. Second, the cell number density can become negative. When there are a small number of molecules per cell in the system, the behavior of the FHD description (19) depends sensitively on the averaging function $n\u223c$ and the propensity functions *a*_{r} that appear in the multiplicative noise terms. Hence, we carefully modify the form of $n\u223c$ and *a*_{r} for negative or very small densities in order to greatly reduce the chances of producing future negative densities.

In Section V A, we demonstrate that the arithmetic mean produces more accurate results for the equilibrium distribution than the other Pythagorean means. Based on the analysis given in Appendix B, we use the following modification to the arithmetic mean:

where

is a smoothed Heaviside function. The smoothed Heaviside function *H*_{0} is introduced to ensure the continuity of $n\u223c$ at *n*_{1} = 0 or *n*_{2} = 0. As explained in Appendix B, this averaging function guarantees non-negativity for the diffusion-only system (18) in the continuous-time description. In our simulations, we find this modification greatly reduces the occurrence of negative density while closely matching the true equilibrium distribution, noting that in our formulation the stochastic diffusive flux is *continuously* turned off at $n1\u22640$ or $n2\u22640$. We also note that the smoothing is based on the number of molecules in a cell and if both cells have at least one molecule (i.e., $ni\Delta V\u22651$), $n\u223c$ becomes exactly the arithmetic mean. As shown in Appendix B, the local modification near *n* = 0 does not cause any noticeable unphysical behavior for $ni\Delta V\u22651$. In Section V, we demonstrate that our numerical schemes based on Eq. (20) work very well even for a small number of molecules per cell.

For the propensity functions $ar(\bm{n})$, we use the following correction to the law of mass action, which is usually included in the RDME description: if the deterministic rate expression contains $ns2$ (or $ns3$, $\cdots $), then replace it by $ns(ns\u22121\Delta V)$ (or $ns(ns\u22121\Delta V)(ns\u22122\Delta V)$, $\cdots $). With this correction, at thermodynamic equilibrium, the mean reaction rate becomes equal to the one calculated from the deterministic rate expression with the mean number density. This can be seen from the fact that if $ns\Delta V$ follows Poisson statistics with mean $n\xafs\Delta V$, $\u27e8ns(ns\u22121\Delta V)\u27e9=n\xafs2$ and $\u27e8ns(ns\u22121\Delta V)(ns\u22122\Delta V)\u27e9=n\xafs3$.

When reactions are combined with an FHD treatment of diffusion, number densities are no longer restricted to non-negative integers and special treatment is required to make reaction rates non-negative and physically sensible for small numbers of molecules. In this work, we evaluate the rate $ar(\bm{n})$ by using continuous-range number densities ** n** (i.e., without trying to round $\bm{n}\Delta V$ to integers) and ensure that each term in the rate of each reaction is non-negative. For example, we take the rate expression of the Schlögl model (see Section II C) to be

where $n+=max(n,0)$. We note that more mathematically justified algorithms have been proposed to handle reactions in regard to negative densities using operator splitting and exact solutions of reaction subproblems;^{99,100} these methods cannot be used to address negative densities due to stochastic diffusive fluxes.

## IV. TEMPORAL INTEGRATORS

In this section, we develop temporal integrators for the spatially discretized FHD equation (19). Our goal is twofold. First, we construct numerical methods that allow for a large time step size even in the presence of fast diffusion. By treating diffusion implicitly, the severe restriction on time step size can be bypassed. Second, we construct methods that maintain accuracy even if the time step size is much larger than the diffusive hopping time. Since it is quite difficult to achieve second-order weak accuracy for general multiplicative noise,^{101} our goal here is to ensure second-order accuracy wherever possible. In the limit in which the number of molecules per cell is very large and one can replace random numbers by their means, our schemes reduce to standard second-order schemes for deterministic reaction-diffusion PDEs. For linearized FHD, our midpoint tau leaping-based schemes are second-order weakly accurate, and all midpoint schemes reproduce at least second-order accurate static correlations, i.e., structure factors.

We build on previous work by some of us in Refs. 63, 64, and 74 and propose two (semi-) implicit schemes as an alternative numerical method to conventional RDME methods. We mainly consider the case where diffusion is much faster than reaction and molecules on average diffuse more than a cell length per time step (i.e., $2dD\Delta t\u226b\Delta x2$). We focus here on *unsplit* schemes that do not rely on operator splitting. This is because we found that unsplit schemes give notably more accurate structure factors than corresponding split schemes in our case. In addition, including other transport processes (e.g., advection) and handling boundary conditions^{102} to second order is not straightforward for split schemes.

It is convenient to introduce dimensionless numbers, $\alpha $ and $\beta $, which measure how fast diffusion and reaction are relative to the given time step size $\Delta t$, respectively. For the single-species Equation (12), assuming $\Delta x1=\cdots =\Delta xd=\Delta x$, we define

where *r* is the chemical relaxation rate appearing in the linearized equation (12). Hence, we can express the well-mixed condition (i.e., the penetration depth $\u2113=D/r\u226b\Delta x$) as $\alpha \u226a\beta $. In addition, the numerical stability condition of a scheme can also be given in terms of $\alpha $ and $\beta $. That is, if reaction and/or diffusion are treated explicitly in a scheme, values of $\alpha $ and/or $\beta $ larger than a stability threshold cause numerical instability. For $\alpha \u226a\beta $, the stability limit is mainly determined by fast diffusion,

Note that the stability limit becomes severe for large diffusion coefficients and small grid spacing and worsens with increasing dimension.

In Section IV A, we present several numerical schemes for the FHD equation (19), including two implicit schemes, and analyze the temporal orders of accuracy for the structure factors using the linearized analysis described in Appendix C. In Section IV B, we analyze the stochastic accuracy of the numerical schemes for large $\Delta t$ by investigating the structure factor of the one-dimensional Schlögl model at different wavenumbers. Since analysis for the nonlinear equations is lacking at present, we numerically justify the handling of multiplicative noise in Section V A.

### A. Schemes

The simplest method for integrating Eq. (19) in time is the *Euler–Maruyama tau leaping* (EMTau) scheme,

where superscripts denote the point in time at which quantities are evaluated, e.g., $nsk=ns(k\Delta t)$ and $ark=ar(\bm{n}(k\Delta t))$, and we have used the compact notation for spatial discretization introduced in Section III B. Here, $\u222bk\Delta t(k+1)\Delta tW(t\u2032)dt\u2032$ has been replaced by $\Delta tWk$, where *W*^{k} denotes a collection of standard random Gaussian variables sampled independently for each species on each grid face at each time step. That is, the stochastic diffusive flux of species *s* on face ** f** at time step

*k*is proportional to $Ws,\bm{f}k$.

We also construct numerical schemes where reactions are treated by SSA, which is an exact (exponential) integrator for reactions. We denote by $\u211cs(\bm{n},\tau )$ the (random) change in the number density of species *s* for a cell with initial state ** n** obtained from SSA over at time interval $\tau $ (in the absence of diffusion). We can then write the

*Euler–Maruyama SSA*(EM-SSA) scheme as

The EMTau scheme is *explicit* in the sense that all terms on the right-hand side of Eq. (25) can be evaluated without knowing $nsk+1$. However, a by a stability condition (for derivation, see Eq. (C7))

which reduces to condition (24) for $\alpha \u226a\beta $. Since the EM-SSA scheme treats reactions using an exponential integrator, it is only subject to the stability limit (24) without a restriction on $\alpha $.

The stability limit imposed by fast diffusion can be overcome by using standard *implicit* methods such as the second-order implicit midpoint or Crank–Nicolson method, which gives a system of linear equations for $nsk+1$,

The linear system (28) can be solved efficiently iteratively using multigrid relaxation;^{103} for $\beta \u22721$, solving the linear system is not much more expensive than a step of a second-order explicit time stepping scheme. Note, however, that scheme (28) is only first order accurate overall for reaction-diffusion systems. Hence, in addition to improved *stability*, it is important to develop higher-order schemes to improve *accuracy*. Note that this is not as simple as replacing tau leaping in Eq. (28) with SSA as in Eq. (26); this would still be only first-order accurate even in the deterministic limit.

Here we construct numerical schemes based on the second-order temporal integrators for the linearized equations of FHD developed in Refs. 64 and 74. Those temporal integrators are second-order accurate in the weak sense for additive noise and are used here as the basis for handling diffusion. In order to add reactions into diffusion-only schemes, we consider two types of sampling methods: tau leaping and SSA. For tau leaping, we use the weakly second-order tau leaping method;^{65,66} a similar two-stage scheme has been originally proposed for the CLE^{104} to achieve second-order weak accuracy. Here we combine predictor-corrector midpoint schemes proposed in Ref. 64 (for diffusion) and the second-order tau leaping method (for reaction). Owing to the similar two-stage structures of those schemes, they fit together in a rather natural manner. In addition, since the resulting schemes still fit the framework of the implicit-explicit algorithms analyzed in Ref. 64, they are second-order weakly accurate for additive noise and an additional order of accuracy (i.e., third order) is gained for the structure factor.

We also develop midpoint schemes that use SSA instead of tau leaping for reactions. Unlike tau leaping-based schemes, the SSA-based schemes do not suffer from instability even in the presence of rapid reactions. The use of SSA may also help prevent the development of negative densities, which is one of the main numerical issues for large fluctuations. Hence, while SSA-based numerical schemes are computationally more expensive, they work better than tau leaping-based schemes when reactions are fast or when the number of molecules is small. The SSA-based schemes we propose here belong to a class of exponential Runge–Kutta schemes, and we construct them to ensure second-order deterministic accuracy, as well as second-order accuracy for the structure factor; a detailed analysis of their weak accuracy is at present missing even for linearized FHD.

#### 1. Explicit midpoint schemes

As a prelude to constructing two-stage implicit methods, we first consider improving the accuracy of the explicit EMTau scheme (25) by using an explicit two-stage Runge–Kutta (predictor-corrector) approach. By combining the explicit midpoint predictor-corrector scheme from Refs. 64 and 74 (for diffusion) and the midpoint tau leaping scheme from Refs. 65 and 66 (for reaction), we obtain the *explicit midpoint tau leaping* (ExMidTau) scheme,

where the superscripts (1) and (2) indicate that the terms correspond to the first and second half of the time step, respectively. That is, $P(1)$ (and similarly for *W*^{(1)} and other random increments) denotes the *same* random number in both predictor and corrector stages and is only sampled once per time step. Following Refs. 65 and 66, the mean reaction rate for the second half step is corrected to $(2ar\u22c6\u2212ark)+$, where $a+=max(a,0)$.

For the magnitude of the stochastic diffusive fluxes over the second half of the time step, we consider the following three options for the face average value $n\u223cs\u2219$:

While all options are consistent with the Ito interpretation, the effect of this choice on accuracy requires a nonlinear analysis that is not available at present. The option (30b) is used in Ref. 74 and shown to lead to second-order weak accuracy for FHD equations linearized around a time-dependent macroscopic state. The option (30c) is inspired by the midpoint tau leaping scheme.^{65,66} However, it does not actually lead to second-order weak accuracy for multiplicative noise because the fluctuating diffusion equation does not have the simple noise structure that the CLE has.^{104} For all our simulations, we use option (30c), as justified by numerical results in Section V A.

The reactions can also be treated using SSA, to give the *explicit midpoint SSA* (ExMidSSA) scheme,

Here the predictor stage (31a) + (31b) is a split reaction-diffusion step, but the corrector is not split. Note that two $\u211c(1)$ appearing in Eqs. (31b) and (31c) are the same random increment computed using SSA. In other words, the SSA algorithm is called once for each half of the time step; this has the same computational cost as calling SSA once to compute $\u211cs(\bm{n}k,\Delta t)$ in the EM-SSA scheme (26).

Since both ExMidTau and ExMidSSA schemes treat diffusion explicitly, they are subject to stability limits. The ExMidTau scheme has the same stability limit (27) as the EMTau scheme, whereas the ExMidSSA scheme is subject to the same limit (24) as the EM-SSA scheme. The ExMidTau scheme with the option (30b) is an instance of the explicit midpoint scheme analyzed in Ref. 74 for weak noise (i.e., linearized FHD) and therefore achieves second-order weak accuracy for linearized FHD and gives third-order accurate equilibrium structure factors. On the other hand, the ExMidSSA scheme gives only second-order accurate structure factors.

#### 2. Implicit midpoint schemes

Here we present two implicit midpoint schemes, where diffusion is treated implicitly based on the implicit midpoint predictor-corrector scheme.^{64,74} By treating reactions using the second-order midpoint tau leaping scheme,^{65,66} we obtain the *implicit midpoint tau leaping* (ImMidTau) scheme,

where the three options for $n\u223cs\u2219$ are given in Eqs. (30). When SSA is used for the reactions, we obtain the *implicit midpoint SSA* (ImMidSSA) scheme,

In the corrector stage, both schemes treat diffusion using the Crank–Nicolson method since this gives the most accurate structure factors for diffusion-only systems.^{64} For the predictor step to the midpoint, we have chosen to use backward Euler for diffusion because this was found to be optimal using the structure factor analysis discussed in more detail in Section IV B.

We point out again that the difference in how reactions are included in the ImMidTau and ImMidSSA schemes stems from the fact that SSA is an exponential integrator whereas the midpoint tau leaping method is only a second-order integrator. This difference must be taken into account when analyzing the accuracy of SSA-based schemes both in the deterministic limit and in structure factor analysis. Like the explicit midpoint schemes in Section IV A 1, for linearized FHD, the ImMidTau scheme is second-order weakly accurate and gives a third-order accurate structure factor, whereas the ImMidSSA scheme gives only a second-order accurate structure factor. Since diffusion is treated implicitly in both schemes, they are *not* subject to a stability limit depending on $\beta $. However, due to the explicit treatment of reactions, the ImMidTau scheme is subject to the stability condition $\alpha \u22642$. The ImMidSSA scheme is unconditionally linearly stable and has no stability restrictions on $\alpha $ and $\beta $ but can be considerably more expensive for systems with a large number of molecules.

A key element of this work that distinguishes it from our previous work based on the CLE^{61} is that here we replaced the GWN in the chemical noise with Poisson noise and used a weakly second-order tau leaping method^{65,66} to account for the non-Gaussian nature of the chemical fluctuations. It is important to note that Poisson noise does not have a continuous range limit, i.e., the Poisson distribution remains integer-valued even as the number of molecules per cell becomes very large. Although it is tempting to replace the Poisson distribution with a Gaussian distribution, this changes the large deviation functional and therefore we recommend using tau leaping even in the case of weak fluctuations; we note that sampling from a Poisson distribution can be done with a cost essentially independent of the mean using well-designed rejection Monte Carlo methods. Because of the use of Poisson variates, which cannot be split into a mean and a fluctuation like a Gaussian variate can, there is no strict “deterministic limit” for our FHD discretizations. While the handling of diffusion degenerates to a standard second-order deterministic scheme in the absence of the noise, the chemical noise is always present and increments or decrements the number of molecules by integer numbers.

### B. Structure factor analysis

Analyzing the accuracy of temporal integrators for stochastic differential equations is notably nontrivial, especially if driven by multiplicative noise. As mentioned above, because of the multiplicative noise, all of our midpoint schemes are formally only first-order weakly accurate. However, traditional weak-order accuracy is not the most important goal in FHD simulations. As first argued in Ref. 63 and then elaborated in Refs. 64 and 74, for FHD it is more important to attain *discrete fluctuation-dissipation balance* and higher-order accuracy for the spectrum of the equilibrium fluctuations.

Here, we analyze the accuracy of our numerical schemes by investigating the structure factor *S*(*k*) for the one-dimensional linearized FHD equation (12). The analytic expression for *S*(*k*) produced by a given scheme can be obtained as a function of $\Delta x$ and $\Delta t$ following the procedure described in Appendix C. Of specific interest to us is how accurately the implicit schemes reproduce *S*(*k*) at large wavenumbers corresponding to length scales comparable to $\Delta x$ (i.e., $k\u2273(D\Delta t)\u22121/2$) when diffusion is the fastest process, $\beta \u226b1\u226b\alpha $. This is because incorrect diffusive dynamics at grid scales for $\Delta t\u226b\Delta x2/D$ can lead to gross errors in the magnitude of the fluctuations at large wavenumbers.

Errors in the structure factor arise from two sources: spatial and temporal discretization. As explained in Appendix C, the predominant contribution of spatial discretization is to replace −*k*^{2} with the symbol of the standard discrete Laplacian. In one dimension, this simply amounts to replacing *k* in the continuum expressions with the modified wavenumber $k\u223c$ defined by

Note that exactly the same expression applies to the RDME, where diffusion is simulated by lattice hops. In order to focus our attention on temporal integration errors, we will plot discrete structure factors as a function of $k\u223c$ instead of *k*, which effectively removes the spatial errors.

Figure 2 illustrates how *S*(*k*) deviates from the exact result (14) at different wavenumbers for large $\Delta t$. We compare *S*(*k*) obtained from the four schemes by using values $\alpha =0.025$, $\beta =10\alpha =0.25$ for the explicit schemes and $\alpha =0.5$, $\beta =10\alpha =5$ for the implicit schemes. Note that these values correspond to a case where the time step size is chosen as half of the stability limit $\Delta tmax$ for the explicit schemes, and it is increased by a factor of 20 ($\Delta t=10\Delta tmax$) for the implicit schemes. As described below, the accuracy at diffusion-dominated scales $k\u2113\u226b1$ and reaction-dominated scales $k\u2113\u226a1$ largely depends on how diffusion (i.e., explicit or implicit) and reaction (i.e., tau leaping or SSA) are treated, respectively.

As the time step size approaches the diffusive stability limit, $\beta \u21921/2$ in one dimension, the explicit schemes become inaccurate and eventually numerically unstable at the largest wavenumbers, as seen in the figure for $\beta =1/4$. Due to the small values of $\alpha $, both schemes give accurate results for reaction-dominated scales. Note, however, that the SSA-based schemes give exact *S*(0) and they are in general more accurate at reaction-dominated scales than the corresponding tau leaping-based schemes. Hence, we see that, for $\beta \u226b\alpha $, the accuracy of the explicit schemes is largely affected by numerical instability arising from diffusion.

On the other hand, both implicit schemes give fairly good *S*(*k*) in the overall range of *k* even though $\Delta t$ is twenty times larger and $\beta =5$. As expected, the ImMidSSA scheme is more accurate for reaction-dominated scales $k\u2113\u226a1$, for which the ImMidTau scheme shows some errors because of the relatively large value of $\alpha =0.5$. However, for intermediate scales $k\u2113\u223c1$, the ImMidTau scheme is more accurate because it attains third-order accuracy for static covariances.

At diffusion-dominated scales $k\u2113\u226b1$, both implicit schemes give accurate results. This is not accidental, for we have selected these schemes from a family of schemes parameterized in Ref. 64 exactly for this reason. Specifically, the treatment of diffusion in both schemes is based on the implicit midpoint scheme (Crank–Nicolson), which gives the exact *S*(*k*) in the absence of reactions.^{64} In addition, the reaction is incorporated in such a way that maintains fluctuation-dissipation balance for $k\u2113\u226b1$ even for relatively large values of $\alpha $. For small $\alpha $, the time integration error of the ImMidTau scheme for the structure factor at the maximum wavenumber $kmax\u2113=\pi \beta /\alpha $ is estimated as

where $S0(k)=lim\Delta t\u21920S(k)$ is the structure factor in the absence of temporal integration errors (see Appendix C). Hence, for a given value of $\alpha $, *S*(*k*) gives accurate results at $k\u2113\u226b1$ for large $\beta $. For the ImMidSSA scheme, $\beta $ in the numerator is replaced by $\beta +(1+2\beta )(\Gamma \xafn\xafr\u22121)$ and a similar stable behavior for large $\beta $ is observed. By expanding *S*(*k*_{max}) − *S*_{0}(*k*_{max}) for small $\Delta t$, we also see that the error is $O(\alpha 2\beta )=O(\Delta t3)$ for the ImMidTau scheme, whereas it is $O(\Delta t2)$ for the ImMidSSA except at thermodynamic equilibrium (i.e., except when $\Gamma \xaf=n\xafr$), where it is third-order accurate.

## V. NUMERICAL RESULTS

We perform numerical simulations for the following three stochastic reaction-diffusion systems. The simulation code is available at https://github.com/BoxLib-Codes/FHD_ReactDiff.git. In Section V A, we use the equilibrium Schlögl model in one, two, and three dimensions to validate our numerical methods. The analysis in Section IV B assumed additive noise, reflecting a large number of molecules per cell. Here we present numerical results demonstrating that the methodology continues to work when there are a small number of molecules per cell and the effects of multiplicative noise are significant. In particular, we show that our numerical methods, including the modified arithmetic-mean averaging function discussed in Section III C, accurately reproduce the Poisson statistics that characterize the thermodynamic equilibrium distribution. In Sections V B and V C, we study the effects of fluctuations on chemical pattern formation. In Section V B, we test our numerical methods on a time-dependent problem: two-dimensional Turing-like pattern formation in the three-species BPM model.^{77,78} We investigate how accurately both time-transient and steady state behaviors are captured for the ImMidTau and ImMidSSA schemes when a large time step size is used. We consider the case where the populations of chemical species have different orders of magnitude, which is a frequently encountered situation where a conventional RDME-based method may not work efficiently. We demonstrate that the ImMidTau scheme scales very well with an increasing number of molecules per cell so that even the deterministic limit of vanishing fluctuations can be explored. In Section V C, we demonstrate the scalability to large systems and computational efficiency of our FHD approach by presenting a three-dimensional numerical simulation of chemical front propagation in a two-species model.^{7}

As a reference method for comparison, we use an RDME-based method, as proposed in Appendix A, which is constructed via a standard operator splitting technique by combining multinomial diffusion sampling^{71} and SSA. Such a split scheme is notably more efficient than ISSA when there are a large number of molecules per cell and becomes an exact sampling method for the RDME in the limit $\Delta t\u21920$. This RDME-based scheme works with non-negative integer populations and reproduces correct fluctuations at thermodynamic equilibrium. However, diffusion imposes the same restriction (24) on $\Delta t$, and the split scheme produces only a first-order accurate structure factor in general.

In order to set a desired magnitude of fluctuations without changing any parameters for the macroscopic limit (e.g., penetration depth), we introduce a factor *A*, which scales the cell volume $\Delta V=A\Delta x1\cdots \Delta xd$. It can be interpreted as the surface cross section in one dimension and the thickness of a system in two dimensions, and as a rescaling of the number density in three dimensions. Since the number of molecules in a cell is $ns,\bm{i}\Delta V$, the larger *A* is, the more molecules in a cell there are and the weaker the fluctuations become. However, the corresponding macroscopic system is unaffected by the value of *A*.

### A. Schlögl model at thermodynamic equilibrium

In this section, we test the numerical schemes constructed in Sections III and IV on the Schlögl reaction-diffusion model, which is first introduced in Section II C. Simulation parameters are chosen to correspond to a system in thermodynamic equilibrium, so that the equilibrium fluctuations are Poisson and the structure factor *S*(*k*) = *n*_{eq} is constant, both with and without (i.e., diffusion only) chemical reactions. Specifically, we set the rate constants as *k*_{1} = *k*_{2} = *k*_{3} = *k*_{4} = 0.1 (see Eq. (15)), which gives *n*_{eq} = 1 and $\alpha =0.2\Delta t$. We set the diffusion coefficient *D* = 1 and the grid spacing to unity for *d* = 1, 2, 3 and thus $\beta =\Delta t$. We consider the case where the mean number of molecules per cell is 10 by setting *A* = 10.

#### 1. Continuous-time FHD equation

Prior to evaluating the different temporal integration strategies, we first focus our attention on the continuous-time discrete-space FHD equation (19) to establish a baseline for comparison and to evaluate the effectiveness of the choice of averaging function (20). To eliminate temporal integration errors, we use the EMTau scheme (25) with a very small time step size $\Delta t=10\u22123$; results from the other FHD schemes are similar for sufficiently small $\Delta t$. The left panel of Fig. 3 shows the structure factors *S*(*k*) computed for *N*_{c} = 512 grid cells in one dimension for the arithmetic, geometric, and harmonic mean (AM, GM, HM) averaging functions $n\u223c(n1,n2)$. The correct flat spectrum is accurately reproduced by the modified AM averaging function (20). On the other hand, the GM and HM averaging functions give smaller *S*(*k*) at diffusion-dominated scales $k\u2113\u226b1$. This can be also observed from the diffusion-only system for all wavenumbers, as theoretically explained in Appendix B. Henceforth, we use the modified AM averaging function (20).

The right panel of Fig. 3 shows that using the AM averaging function, the correct Poisson distribution for the number *N* of molecules in a cell is accurately reproduced for both reaction-diffusion and diffusion-only systems at thermodynamic equilibrium. From the equilibrium number density distribution $\rho (n)$, we construct a discrete distribution for integer number of molecules *N* per cell,

and compare it with a Poisson distribution *P*_{Poisson}(*N*) with mean $neq\Delta V$, as well as a Gaussian distribution *P*_{Gauss}(*N*) having the same mean and variance as *P*_{Poisson}(*N*). The agreement between *P*(*N*) and *P*_{Poisson}(*N*) is remarkable in the sense that FHD was originally proposed to account for only second moments of (small) Gaussian fluctuations. Since *P*_{Poisson}(*N*) is significantly different from *P*_{Gauss}(*N*) for $neq\Delta V=10$, we confirm that our spatially discretized FHD equation (19) describes (large) Poisson fluctuations faithfully.

#### 2. Time integration errors

In order to investigate time integration errors of our numerical schemes, we compare the numerical equilibrium distribution for a given time step size $\Delta t$ with the target Poisson distribution *P*_{Poisson}(*N*) by using the following measures. First, we compute the Kullback–Leibler (KL) divergence (distance),

Second, we compute the probability of negative number densities,

Third, we compute the correlation coefficient $\zeta $ between neighboring cells,

Note that all three measures should be zero at thermodynamic equilibrium, as they would be for RDME.

We considered the three options for the stochastic flux amplitude $n\u223cs\u2219$ in Eqs. (30). For $\zeta $, the three options give similar values within standard errors of estimation. For *D*_{KL} and *P*_{neg}, option (30a) gives the largest values (i.e., least accurate) and option (30c) the smallest (not shown). Based on this result, we will adopt (30c) and use it for all of the simulations. Figure 4 shows how these measures deviate from zero as $\Delta t$ increases for *N*_{c} = 64 cells in one dimension for the different schemes. As expected, for small values of $\Delta t$, all schemes give similar values. *D*_{KL} converges to a small value, which is consistent with the good agreement between *P*(*N*) and *P*_{Poisson}(*N*) seen in the right panel of Fig. 3. *P*_{neg} is observed to converge to zero as $\Delta t\u21920$, which demonstrates the effectiveness of the approach described in Section III C and agrees with the analysis in Appendix B. Also, no correlation between neighboring cells is observed for small $\Delta t$ within statistical errors, which is consistent with the flat spectrum *S*(*k*) shown in the left panel of Fig. 3. As $\Delta t$ approaches the explicit stability limit $\Delta tmax$, rapid worsening is observed in both explicit schemes in all three measures. While *P*_{neg} and $\zeta $ behave similarly for both schemes, *D*_{KL} remains small for larger values of $\Delta t$ for the ExMidSSA scheme compared to the ExMidTau scheme.

For both implicit schemes, it can be clearly seen that not only is the diffusion instability bypassed but also the accuracy is well maintained for large $\Delta t$. For comparable accuracy, an order of magnitude larger time step size than the explicit schemes can be chosen. The ImMidTau scheme gives smaller *D*_{KL} than the ImMidSSA scheme if $\Delta t$ is smaller than a certain value. This is consistent with the observation that the former scheme has a higher temporal order of accuracy in *S*(*k*). However, due to inaccurate handling of reactions by tau leaping for large $\Delta t$, the ImMidTau scheme eventually gives larger *D*_{KL}. For *P*_{neg} and $\zeta $, a similar behavior is observed.

A similar behavior is observed for higher spatial dimensions *d* = 2 and *d* = 3 (not shown). However, for a given target accuracy tolerance, we find that a smaller time step size should be chosen, which is inversely proportional to the dimensionality *d*. This should not come as a surprise since the explicit stability limit $\Delta tmax\u223c1/d$, see Eq. (24). Therefore, we conclude that $\Delta t$ can be chosen as an order of magnitude larger for the implicit schemes than for the explicit schemes *independent* of the spatial dimension. As mentioned, the computational overhead for solving linear systems can be reduced by an efficient iterative solver. Using multigrid relaxation,^{103} the overall computational efficiency gain was roughly estimated to be a factor of 3. However, this factor largely depends on the problem as well as the implementation, especially on the linear solver used.

### B. Turing-like pattern formation

In this section, we investigate pattern formation in the three-species Baras–Pearson–Mansour (BPM) model,^{77,78}

We choose the rate constants so that the deterministic reaction-only system attains a limit cycle as its stable attractor, and we choose the diffusion coefficients so that a Turing-like pattern forms in the reaction-diffusion system.^{61} Specifically, we set $k1=k2=2\xd710\u22124$, *k*_{3} = 1, $k4=3.33\xd710\u22123$, *k*_{5} = 16.7, $k6=3.67\xd710\u22122$, *k*_{7} = 4.44, and *D*_{U} = 0.1, *D*_{V} = *D*_{W} = 0.01. We note that on the limit cycle, number densities of the three species oscillate in significantly different ranges: $nU\u2208(999,2024)$, $nV\u2208(302,645)$, and $nW\u2208(18.2,83.2)$. For a physical domain with side lengths *L*_{x} = *L*_{y} = 32, we use three spatial resolutions with grid sizes *N*_{c} = 64^{2}, 128^{2}, and 256^{2} cells. For the initial number densities, we choose a point on the limit cycle $(nU0,nV0,nW0)=(1686,534,56.4)$ and generate the initial number of molecules of each species *s* in each cell from a Poisson distribution with mean $ns0\Delta V=ns0A\Delta x\Delta y$. We use our implicit schemes in order to bypass the stiff stability limit imposed by the fast diffusion of U molecules. To obtain *reference FHD results* having minimal time integration errors, we use the ImMidTau scheme with $\Delta t=0.1$. To test the importance of fluctuations, a deterministic version of the ImMidTau scheme is used with $\Delta t=0.1$, with random initial conditions generated from a Poisson distribution corresponding to thickness *A* = 10.

Figure 5 shows snapshots of a final Turing pattern formed for *A* = 1 and 64^{2} cells. While the pattern is qualitatively correct, the quantitative behavior of our FHD formulation may be questioned since the mean number of W molecules in a cell can be as low as 4.5 at small *t* in this case. To confirm that the FHD description applies even for relatively small numbers of molecules per cell, we compare the FHD results to reference RDME results obtained using the SSA/2+MN+SSA/2 scheme (A3) with $\Delta t=0.01$. We find that the FHD reference simulations are qualitatively very similar to the RDME reference simulations over a wide range of thicknesses *A*, as we illustrate in Fig. 5. For our setup, after the initial formation of a disordered pattern of dots with low concentration of U molecules (blue dots in Fig. 5), the dots split and merge and diffuse to eventually form a stable regular pattern; note that the final patterns are nearly periodic lattice structures but their geometry is affected to some extent by the finite size of the domain. For *A* = 1, by $t=5\xd7104$, almost all samples had formed a steady pattern. Most samples formed a hexagonal (12 dots, see panels (a) and (c) in Fig. 5), and a few formed a monoclinic (11 dots, see panels (b) and (d) in Fig. 5) lattice of dots, for both FHD and RDME. Note that while FHD simulations using the ImMidTau scheme are equally efficient independent of *A*, RDME simulations become prohibitively expensive for large $A\u2273100$ due to the very large number of U molecules (as many as $2\xd7106A$) in the system. For weaker fluctuations, *A* = 1000, FHD simulations reveal that the annealing of the lattice defects takes much longer and we see several disordered or defective patterns even at $t=5\xd7104$ (not shown). Therefore, not only do fluctuations accelerate the formation of the initial (disordered) pattern, but they also appear to accelerate the annealing of the defects.

Since the formation of the pattern is driven by instability, it is itself a random process and a proper quantitative comparison between the different methods requires a careful statistical analysis of an ensemble of trajectories. In order to capture the time transient behavior of pattern formation, illustrated in the right panel of Fig. 6, we calculate the spatially averaged density $n\xafU(t)=1Nc\u2211\bm{i}nU,\bm{i}(t)$. In the left panel, we compare sample trajectories of $n\xafU(t)$ for *A* = 1 and *A* = 10 for RDME (similar results are obtained for FHD) and for deterministic reaction-diffusion. While $n\xafU(t)$ initially oscillates as in the limit cycle, as the Turing-like pattern begins to form, the oscillation amplitude decays and $n\xafU(t)$ eventually attains a steady value. By comparing *A* = 1 and 10, we see that larger fluctuations facilitate faster pattern formation, as observed in prior work^{61} by us and others. By comparing RDME results for thickness *A* = 10 with deterministic reaction-diffusion started from the same initial condition, we see that the effect of fluctuations on pattern formation is *not* just due to random initial conditions.

We generate 16 sample trajectories for each set of parameter values, fit each realization of $n\xafU(t)$ using seven fitting parameters $a1,\u2026,a7$ to

and compare the distributions of the fitting parameters. Note that *a*_{1} and *a*_{7} correspond to the decay onset time and the steady spatial average density, respectively. In the left panel of Fig. 7, we compare the empirical distributions of (*a*_{1}, *a*_{7}) from the RDME and FHD results for different values of the thickness *A* and spatial resolutions. For each value of *A*, we observe that distributions obtained from different methods and/or resolutions coincide. For *A* = 1 and *A* = 10, we reconfirm that the RDME approach produces statistically very similar results for three resolutions and the FHD results are statistically indistinguishable from the RDME results. It is quite remarkable that FHD works even for *A* = 1 and 256^{2} cells, which can have as low as 0.3 molecules per cell at small *t*, and this demonstrates the robustness of our treatment for a small number of molecules per cell. As the magnitude of the fluctuations increases (i.e., as *A* decreases), the pattern begins to form earlier (i.e., *a*_{1} decreases), as already seen in Fig. 6, while the steady spatial average density *a*_{7} becomes smaller. In addition, while the variance of *a*_{1} does not change significantly as *A* varies, the variance of *a*_{7} becomes larger for smaller *A*.

Finally, we investigate time integration errors of our implicit schemes for the Turing-like pattern formation. For the ImMidTau scheme, we increase $\Delta t$ up to the stability limit arising from reactions $\Delta tmax(R)\u22481.3$. The right panel of Fig. 7 shows the mean values $a\xaf7$ of the steady spatial average density *a*_{7} over 16 samples versus $\Delta t$ for 64^{2} cells. While both schemes give similar values to the RDME results for small $\Delta t$, they show a different behavior for large $\Delta t$, which also depends on *A*. As expected, in the ImMidTau scheme, the value of $a\xaf7$ rapidly deviates from the RDME result as $\Delta t$ approaches $\Delta tmax(R)$, especially if there are few molecules per cell, *A* = 0.1 or *A* = 1. On the other hand, in the ImMidSSA scheme, $\Delta t$ can be increased beyond $\Delta tmax(R)$, and deviations from the RDME results remain small even for the smallest value of *A*. Hence, handling reactions by SSA not only removes the reaction stability constraint but also improves the accuracy for a small number of molecules per cell. However, it should be noted that this improvement comes at a significant computational cost, since the SSA scheme is much more expensive than tau leaping especially as the number of molecules per cell increases. Therefore, the SSA-based schemes are impractical in the regime of weak fluctuations due to poor scaling. We discuss some alternatives to SSA that may significantly improve the computational cost for weaker fluctuations in the Conclusions.

### C. Front propagation

As a final example, we simulate three-dimensional front propagation in a two-species stochastic reaction-diffusion system having the following reaction network:

This model has been proposed to reproduce axial segmentation in Ref. 7, where ISSA simulations have been performed for the one-dimensional case. Following Ref. 7, we set *k*_{1} = 0.4, *k*_{2} = 0.137, *k*_{3} = 0.1, *k*_{4} = 1 and *D*_{A} = 1, *D*_{B} = 10. For a physical domain with side lengths *L*_{x} = *L*_{y} = *L*_{z} = 512, we use 256^{3} cells. To initiate front propagation, we generate initial number densities as follows. We first assign to each cell ** i** and species

*s*a mean number density,

where $r\bm{i}$ is the distance from the cell center to the center of the domain. This initializes a spherical region of radius *R* in the first uniform equilibrium state of the model, $(nA(1),nB(1))=(2.16,1.35)$, while the rest of the domain is initialized in the second uniform equilibrium state, $(nA(2),nB(2))=(0,10)$, with a smooth transition region between the two states of width $\u22482\xi $. Then, as in Section V B, we generate the initial number of molecules of each species in each cell from a Poisson distribution with mean $ns,\bm{i}0\Delta V=ns,\bm{i}0A\Delta x\Delta y\Delta z$. We simulate the system for parameters *A* = 1000, *R* = 16, and $\xi =4=2\Delta x$ using the ImMidTau scheme with $\Delta t=0.25$. For comparison, we also simulate the corresponding deterministic system using a deterministic version of the ImMidTau scheme with the same time step size and (random) initial conditions. Simulations are performed using a parallel implementation of the algorithm using the BoxLib software framework.^{105} We emphasize that a corresponding RDME system is too large to simulate with conventional RDME-based methods; while the total number of molecules in the system varies as the front propagates, it is of the order of 10^{12} for *A* = 1000.

Figure 8 shows the growth of the spherical region as the more stable phase propagates into the less stable phase via a spherical traveling wave. While the phase boundary having a peak population of species A propagates, a Turing-like pattern develops behind the wave front; in one dimension this pattern is periodic and more pronounced in the presence of fluctuations.^{7} In two and three dimensions, fluctuations not only enhance the pattern but they also make it disordered, as seen in the figure by comparing the stochastic and deterministic cases. In addition, the phase boundary becomes more irregular under fluctuations. Note that the numerical results for the deterministic case are not perfectly radially symmetric not only due to the noisy initial conditions but also due to grid artifacts introduced by the standard discrete Laplacian, which is not perfectly isotropic^{106} on length scales compared to the front width (i.e., the penetration depth); one would require an even finer grid to correct for this spatial discretization artifact.

## VI. CONCLUSIONS

In this work, we have formulated a fluctuating hydrodynamics (FHD) model for reaction-diffusion systems and developed numerical schemes to solve the resulting stochastic ordinary differential equations (SODEs) (19) for the number densities $ns,\bm{i}(t)$ of chemical species in each cell. We obtained the diffusion part of the SODEs from an FHD description of a microscopic system consisting of molecules undergoing independent Brownian motions and added reactions in an equivalent manner to the reaction-diffusion master equation (RDME). We presented two implicit predictor-corrector schemes, the ImMidTau (32) and ImMidSSA (33) schemes, that treat reactions using tau leaping and SSA, respectively. In these schemes, diffusion is treated implicitly so that the stability limit imposed by fast diffusion can be bypassed and the time step size can be chosen to be significantly larger than the hopping time scale of diffusing molecules. In addition, two-stage Runge–Kutta temporal integrators are employed to improve the accuracy. To confirm the validity of our FHD formulation and demonstrate the performance of our numerical schemes, we numerically investigated not only a system at steady state (Schlögl reaction-diffusion model) but also time-dependent two-dimensional Turing-like pattern formation and three-dimensional front propagation.

Based on our analytical and numerical investigation, we conclude that the ImMidTau scheme is an efficient and robust alternative numerical method for reaction-diffusion simulations. The reason is threefold. First, the cost of the scheme does not increase for increasing number of molecules per cell (weaker fluctuations). For small numbers of molecules per cell (large fluctuations), the integer-valued RDME description is more appropriate than the continuous-range FHD description. However, by using the approach proposed in Section III C, we ensured that the FHD description remains robust and gives accurate results even for a small number of molecules per cell, as shown in Section V A. Hence, as shown in Section V B, our numerical methods can efficiently simulate reaction-diffusion systems over a broad range of relative magnitude of the fluctuations. Second, the scheme allows a significantly larger time step size without degrading accuracy compared to existing RDME-based numerical methods,^{42–47} which use a fixed time step size for diffusion that is comparable to the hopping time scale. In particular, we found that the time step size could be chosen as an order of magnitude larger for the implicit schemes than for explicit methods, independent of the spatial dimension. Lastly, FHD can take advantage of the development of efficient parallel algorithms developed for computational fluid dynamics (CFD) that enable effective use of high-performance parallel architectures while providing the framework for treating more complex problems with additional physical phenomena. This enabled us to perform three-dimensional simulations of chemical front propagation involving as many as 10^{12} molecules using the BoxLib CFD software framework.^{105}

The explicit tau leaping methods used here are quite simple to implement but are subject to a stability limit for fast reactions, and can lead to negative densities when fluctuations are large. While some implicit tau leaping methods have been developed, as an alternative we developed methods that use SSA for reactions. The ExMidSSA and ImMidSSA methods, however, do not scale as the fluctuations become weaker. This deficiency can be corrected by replacing SSA by a recently developed hybrid algorithm termed asynchronous tau leaping^{107} that combines SSA and tau leaping in a dynamic manner by simulating multiple events with asynchronous time steps. Future work should develop FHD-based numerical schemes that are accurate and robust even for a small number of molecules per cell and also scale to the deterministic limit efficiently.

One of the advantages of the FHD approach for reaction-diffusion systems is its natural generalization to more complicated and realistic applications. Chemical reactions of interest usually occur in liquid solution, and often in a dense crowded environment such as the cytoplasm.^{108–110} It is well-known that Brownian motion of liquid molecules or suspended macromolecules in liquids is dominated by hydrodynamic effects related to viscous dissipation.^{81,110,111} This means that the diffusion model commonly used in reaction-diffusion models, including this work, which assumes that reactants are independent non-interacting Brownian walkers diffusing with a constant diffusion coefficient, does not apply in the majority of practical problems of interest. Notably, crowding or steric interactions affect the local hydrodynamic mobility of individual reactants, and hydrodynamic interactions (HI) among the diffusing reactants introduce strong correlations among the diffusive motions of the reacting particles (and also among reactants and passive crowding agents).^{110} Excluded volume due to steric repulsion introduces cross-diffusion effects, i.e., coupling between the diffusive fluxes for different species,^{112} as well as concentration-dependent diffusion coefficients. Furthermore, it has been observed that cross-diffusion may lead to qualitatively different Turing instabilities.^{113–115} Long-ranged contributions of hydrodynamic interactions can be captured by accounting for the advection of concentration fluctuations by the thermal velocity fluctuations, which follow a fluctuating Stokes equation.^{81,111} Additional thermodynamic contributions to the diffusive fluxes such as cross-diffusion, barodiffusion, and thermodiffusion do not seem straightforward in the RDME but are easily included in our FHD formulation.^{61} In future work, we will investigate these hydrodynamic effects on stochastic reaction-diffusion phenomena.

## ACKNOWLEDGMENTS

We would like to thank Jonathan Goodman for numerous informative discussions. This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under Contract No. DE-AC02-05CH11231 and Award No. DE-SC0008271. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

### APPENDIX A: RDME-BASED SPLIT SCHEME

As a reference algorithm for stochastic reaction-diffusion simulation, we construct a numerical scheme for the RDME through a standard operator splitting technique, as done in a number of prior works.^{25,44,45,47} This technique allows one to obtain a numerical scheme for the reaction-diffusion system by combining numerical methods for the diffusion-only and reaction-only systems. Here we combine multinomial diffusion sampling^{43,71} for diffusion and SSA for reaction via Strang splitting.^{70}

One distinguishing feature of the resulting scheme, compared to exact sampling methods such as ISSA, is the use of a fixed time step size $\Delta t$ for diffusion, which facilitates an efficient numerical simulation if diffusion is fast. As shown below, while there are several advantages of this scheme over ISSA, the choice of $\Delta t$ is constrained as it is for explicit FHD schemes.

In Appendix A 1, after reviewing the multinomial diffusion sampling method, we present an RDME-based split scheme and discuss its advantages and disadvantages. In Appendix A 2, we present a stochastic accuracy analysis for this scheme.

#### 1. Split scheme

In the multinomial diffusion (MN) sampling method, the number of molecules in each cell after time $\Delta t$ is calculated by sampling how many molecules have moved from a cell to a neighboring cell. Here we follow the simple algorithm described by Balter and Tartakovsky^{71} and only allow particles to move to nearest-neighbor cells. More complicated but also more accurate algorithms that allow particles to jump to further than nearest-neighbor cells are described by Lampoudi *et al.*^{43} We use the notation introduced in the body of the paper.

By denoting the number of molecules of species *s* diffusing from cell ** i** to cell $\bm{i}\u2032$ over the time interval $\Delta t$ as $Ns,\bm{i}\u2192\bm{i}\u2032$, the change in the number density $ns,\bm{i}$ can be expressed in terms of the sum of the inward and outward fluxes,

where $ns(t)={ns,\bm{i}(t)}$. For each cell ** i**, the outward fluxes $(Ns,\bm{i}\u2192\bm{i}+\bm{e}1,Ns,\bm{i}\u2192\bm{i}\u2212\bm{e}1,\u2026,Ns,\bm{i}\u2192\bm{i}+\bm{e}d,$$Ns,\bm{i}\u2192\bm{i}\u2212\bm{e}d,Ns,\bm{i}\u2192\bm{i})$ are random variables sampled from the multinomial distribution with $\u2211\bm{i}\u2032Ns,\bm{i}\u2192\bm{i}\u2032=ns,\bm{i}(t)\Delta V$ total trials and probabilities $(ps,ps,\u2026,ps,ps,1\u22122dps)$, where $ps=Ds\Delta t/\Delta x2$, where we have assumed $\Delta x1=\cdots =\Delta xd=\Delta x$.

For fast diffusion, this method becomes more efficient than treating hoping events one by one (as in ISSA). However, it is an approximate method since the actual distribution of the outward fluxes deviates from the multinomial distribution as $\Delta t$ increases. In fact, $\Delta t$ cannot be arbitrarily large and is limited by condition (24) because of the requirement $1\u22122dps\u22650$. We also note that the number of molecules in a cell never becomes negative due to the constraint $\u2211\bm{i}\u2032Ns,\bm{i}\u2192\bm{i}\u2032=ns,\bm{i}(t)\Delta V$. Hence, the fluxes on disjoint faces are correlated, which is different from the FHD description (18). In the deterministic limit, this scheme converges to a standard forward Euler scheme for the diffusion equation and is therefore only first-order accurate in time. In the stochastic setting, this scheme adds correlations between the fluxes through different faces of a given cell in such a way as to ensure discrete fluctuation-dissipation balance for *any* allowable time step size. In fact, for a system with diffusion only, this method ensures that the equilibrium fluctuations are strictly Poisson, as desired for independent Brownian walkers.

In order to handle chemical reactions, we use the SSA algorithm locally and independently in each cell, without any diffusive events. Let $\u211cs(\bm{n},\tau )$ denote the (random) change in the number density of species *s* when a cell with initial state ** n** is simulated using SSA over a time interval $\tau $. In the absence of diffusion, the SSA-based reaction scheme can be written as

If we combine the diffusion-only (A1) and reaction-only (A2) schemes using a Strang splitting approach,^{70} we obtain the SSA/2+MN+SSA/2 scheme

where superscripts denote time step or intermediate stage. (We note that a number of different splitting variants are possible; the version presented here gave the most accurate structure factor.) This split scheme has a number of advantages. It becomes an exact sampler (solver) for the RDME in the limit $\Delta t\u21920$, just like ISSA. It is notably more efficient than ISSA if there are many events per time interval $\Delta t$, and it can be parallelized in a straightforward manner using domain decomposition. Since the number of molecules is always a non-negative integer both in multinomial diffusion sampling and in SSA, this property is also preserved for this scheme. Moreover, since both sampling methods preserve the thermodynamic equilibrium distribution (i.e., the Poisson statistics), the split scheme also preserves it for any allowable time step size.

However, this scheme has some disadvantages. First, the time step size restriction (24) for diffusion becomes severe for fast diffusion. Second, SSA exhibits poor scalability with respect to the number of molecules in a cell. This can be resolved by replacing SSA with the tau leaping method, but the non-negativity is no longer guaranteed. Third, since the multinomial diffusion method used here is only first-order accurate, the accuracy of the scheme is first order even though Strang splitting is used. Constructing RDME-based diffusion methods that are more accurate is possible^{43,46} but nontrivial and is not the focus of our work.

#### 2. Structure factor analysis

In this section, we investigate the stochastic accuracy of the SSA/2+MN+SSA/2 scheme (A3). To this end, we consider the structure factor *S*(*k*) of an out-of-equilibrium monostable Schlögl model in one dimension (see Sections II B and II C). In the limit of many molecules per cell, an asymptotic expression of *S*(*k*) can be obtained for the scheme as a function of $\Delta x$ and $\Delta t$. The multinomial fluxes can be approximated by correlated Gaussian ones and the type of analysis summarized in Appendix C can be applied; we do not give the details here for brevity. Note that we present a similar structure factor analysis for our FHD-based schemes with some background and details in Section IV B.

Figure 9 illustrates how *S*(*k*) deviates from the exact result (14) as $\Delta t$ is increased to $\Delta tmax$ (equivalently, to $\beta =0.5$), for $\alpha =0.1\beta $, see Eq. (23). While *S*(*k*) is accurately reproduced at the reaction-dominated scales $k\u2113\u226a1$ for all values of $\beta $, it becomes inaccurate for smaller scales as $\Delta t$ approaches $\Delta tmax$. We recall that for a system at thermodynamic equilibrium, the split scheme *exactly* preserves the correct equilibrium distribution for any $\Delta t<\Delta tmax$. This property ensures that, for $\alpha \u226a\beta $, good structure factors are obtained even for systems outside of thermodynamic equilibrium, which exhibit a nonzero correlation length. For example, for $\beta =0.25$, the SSA/2+MN+SSA/2 scheme in Fig. 9 gives a notably more accurate *S*(*k*) than the FHD-based *explicit* schemes in Fig. 2. Hence, even though this split scheme is found to give only first-order accurate *S*(*k*), the resulting structure factor is relatively insusceptible to increasing $\Delta t$ until the stability limit is approached.

### APPENDIX B: AVERAGING FUNCTION $\mathit{n}\u223c$

In this appendix, we show that the arithmetic mean should be chosen as the averaging function $n\u223c(n1,n2)$. Here we consider the diffusion-only case for a single species and investigate the equilibrium distribution of the spatially discretized FHD equation (18). Since the *true* equilibrium distribution for a bulk (infinite) system is known to be a product Poisson distribution from the corresponding microscopic system consisting of molecules undergoing independent Brownian motions, we choose $n\u223c$ so that the resulting equilibrium distribution is as close as possible to the true distribution. In addition, since the prevention of negative cell number densities is one of the essential issues for the development of a robust FHD numerical scheme, special care is taken to modify the form of $n\u223c(n1,n2)$ for small values of *n*_{1} and *n*_{2}. Here we focus on the continuous-time case and do not assume any specific temporal integrator.

The corresponding microscopic system has $Ncn\xaf\Delta V$ molecules, where $n\xaf$ is the mean number density and *N*_{c} is the number of cells. The equilibrium distribution of the numbers of molecules in each cell follows the multinomial distribution with equal probabilities 1/*N*_{c}. Hence, it is straightforward to obtain the following second-order statistics of cell number density $n\bm{i}$:

If the FHD system (18) attains an equilibrium state, it can be shown that its second-order statistics are completely characterized by $\u27e8n\u223c\u27e9$, which is the equilibrium average of $n\u223c$ over all faces,

Comparing Eqs. (B4)–(B6) to the correct result Eqs. (B1)–(B3) suggests that one needs to choose $n\u223c$ so that $\u27e8n\u223c\u27e9$ is as close as possible to $n\xaf$. It is easy to see that the arithmetic mean (AM) *would* give the right answer: $\u27e8n\u223cAM\u27e9=\u27e812(n1+n2)\u27e9=n\xaf$. On the other hand, to calculate $\u27e8n\u223c\u27e9$ for the geometric mean $n\u223cGM=n1n2$ or the harmonic mean $n\u223cHM=2/(n1\u22121+n2\u22121)$, one needs to know the equilibrium distribution $\rho (n1,n2)$ of two neighboring cells. However, under the reasonable assumption that all three averaging functions give similar distributions $\rho (n1,n2)$ allowing only non-negative number densities, it can be easily shown that $\u27e8n\u223cHM\u27e9\u2264\u27e8n\u223cGM\u27e9\u2264\u27e8n\u223cAM\u27e9$ from the well-known inequalities among the Pythagorean means. In fact, in Fig. 3, this ordering is observed from the structure factor of the diffusion-only system, see Eq. (B6). Hence, we conclude that the arithmetic mean is the right choice for $n\u223c$.

However, if the arithmetic mean is employed without modification, Eq. (18) does not attain an equilibrium state. This is because almost surely at some point on some grid face we will have $n\u223c<0$ so that the stochastic diffusive flux becomes undefined. The non-negativity of cell number densities is guaranteed if the stochastic diffusive flux through a face is turned off when the number density of either cell sharing the face becomes zero. Specifically, under some technical assumptions, it can be proven that if $n\u223c(n1,n2)$ is a non-negative function that satisfies

then the number density of each cell never becomes negative. The non-negative arithmetic mean averaging function,

does not allow negative density (i.e., $\rho (n)=0$ for $n<0$). However, due to the discontinuity of $n\u223c$ at *n*_{1} = 0 or *n*_{2} = 0, $\rho (n)$ does not decrease to zero as *n* becomes zero (i.e., $limn\u21920+\rho (n)>0$) and a delta function is formed at *n* = 0, see Fig. 10.

To understand this behavior, we consider numerically integrating Eq. (18) with a small time step size $\Delta t>0$. Even if the density in a given cell *n*_{1}(*t*) has a small positive value, $n\u223c$ can be large if the density in the neighboring cell *n*_{2}(*t*) is large. In this case, $n1(t+\Delta t)$ can become negative due to the stochastic diffusive flux. But, once *n*_{1} has become negative, the stochastic diffusive flux is turned off, and due to the deterministic diffusion, *n*_{1} increases and becomes positive again. Hence, as shown in Fig. 10, $\rho (n)$ attains a peak in the negative density region near *n* = 0. The width and height of the peak are proportional to $\Delta t1/2$ and $\Delta t\u22121/2$, respectively, and the peak becomes a delta function in the limit $\Delta t\u21920$.

We note that the averaging in Eq. (B8) can be expressed as $n\u223c=12(n1+n2)H(n1\Delta V)H(n2\Delta V)$, where *H* is the Heaviside function. To avoid a discontinuity in $n\u223c$, we can use a smoothed Heaviside function to arrive at Eq. (20). Note that the smoothing is based on the number of molecules in a cell ($N=n\Delta V$) and the smoothing region $0\u2264N\u22641$ is chosen so that the stochastic diffusive flux is modified only when there is less than one molecule in a cell. In Fig. 10, we show the distribution $\rho (n)$ near *n* = 0 obtained by the averaging function (20), for a rather small mean number of molecules per cell, $n\xaf\Delta V=5$. With the use of a smoothed Heaviside function, the spurious delta function at *n* = 0 as $\Delta t\u21920$ is removed, and the probability of negative density is greatly reduced for small $\Delta t$.

### APPENDIX C: LINEARIZED EQUATION ANALYSIS

In this appendix, we summarize how the discrete structure factor is obtained as a function of $\Delta x$ and $\Delta t$ when a given spatiotemporal discretization is applied to the linearized FHD equation (12), following the Fourier-space analysis developed in Ref. 63. For simplicity, we consider here the one-dimensional case.

Applying the spatial discretization given in Section III to Eq. (12) and taking a discrete Fourier transform, we obtain an Ornstein–Uhlenbeck equation for the Fourier coefficient $\delta n^k(t)$,

where the modified wavenumber $k\u223c$ is defined in Eq. (34), and $Wk(D)(t)$ and $Wk(R)(t)$ are independent standard GWN processes. Compared to the continuous-space case (see Eq. (13)), *k* is replaced by $k\u223c$ due to the discrete Laplacian and divergence operators. Note that in this linearized analysis, Poisson processes have been replaced by Gaussian ones with the mean evaluated at the ensemble average (i.e., macroscopic) values of the density. For convenience, we have replaced complex-valued GWN processes by real-valued ones having the same noise intensities.

where *W*_{1} and *W*_{2} are independent standard normal random variables. We write this as

where $=d$ denotes being equal in distribution. For any given temporal integrator, we can straightforwardly obtain analytic expressions for *M*_{k} and $NkNk\u2217$. For example, for the EMTau scheme,

The covariance of the noise $NkNk\u2217$ for multinomial diffusion (A1) can most easily be obtained from Eq. (C8) from the observation that, in the absence of reactions, the exact structure factor $S(k)=n\xaf$ is obtained for any stable time step.

A similar procedure is applicable to numerical schemes having SSA for reactions. Since SSA is an exact integrator, the linearized reaction part in Eq. (C1) is exactly solved. For example, for the EM-SSA scheme (26), we have (cf. Eq. (C2))

While the expressions of *M*_{k} and $NkNk*$ become complicated for the predictor-corrector midpoint schemes, a theoretical analysis is still tractable with the help of symbolic algebra tools.

By calculating *M*_{k} and $NkNk*$ for a given numerical scheme, the stability condition and the structure factor can be obtained as follows. From the condition that the amplification factor *M*_{k} should satisfy

the stability condition is obtained. For the EMTau scheme, we obtain

The analytic expression for $S(k)=V\u27e8\delta n^k\delta n^k*\u27e9$ can be calculated from

which is obtained from the time invariance relation $\u27e8\delta n^k(t)\delta n^k*(t)\u27e9=\u27e8\delta n^k(t+\Delta t)\delta n^k*(t+\Delta t)\u27e9$ and Eq. (C3).^{63} From the analytical expressions of *S*(*k*) and $S0(k)=lim\Delta t\u21920$$S(k)$, [*S*(*k*_{max}) − *S*_{0}(*k*_{max})]/*S*_{0}(*k*_{max}) is easily obtained, the series expansion of which for small $\alpha $ (and fixed $\beta $) gives Eq. (35) for the ImMidTau scheme. Similarly, series expansions for small $\Delta t$ (i.e., fixed ratio $\alpha /\beta $) reveal the temporal order of accuracy of *S*(*k*).