We formulate and study computationally the fluctuating compressible Navier-Stokes equations for reactive multi-species fluid mixtures. We contrast two different expressions for the covariance of the stochastic chemical production rate in the Langevin formulation of stochastic chemistry, and compare both of them to predictions of the chemical master equation for homogeneous well-mixed systems close to and far from thermodynamic equilibrium. We develop a numerical scheme for inhomogeneous reactive flows, based on our previous methods for non-reactive mixtures [Balakrishnan , Phys. Rev. E **89**, 013017 (2014)]. We study the suppression of non-equilibrium long-ranged correlations of concentration fluctuations by chemical reactions, as well as the enhancement of pattern formation by spontaneous fluctuations. Good agreement with available theory demonstrates that the formulation is robust and a useful tool in the study of fluctuations in reactive multi-species fluids. At the same time, several problems with Langevin formulations of stochastic chemistry are identified, suggesting that future work should examine combining Langevin and master equation descriptions of hydrodynamic and chemical fluctuations.

## I. INTRODUCTION

Chemical reactions are of central importance in both natural and industrial processes spanning the range of length scales from the microscopic, through the mesoscopic, and up to macroscopic scales. It is the rule, rather than the exception, that chemical reactions are strongly coupled to hydrodynamic transport processes, such as advection, diffusion, and thermal conduction. Prominent examples include diffusion-limited aggregation, pattern and chemical wave formation in reactive solutions, reaction-driven convective instabilities, heterogeneous catalysis, combustion, complex biological processes, and others. Even in a homogeneous system with only slightly exothermic reactions, the chemistry is coupled to the hydrodynamics, leading to non-trivial effects such as fluctuation induced transitions.^{1}

Fluctuations affect reactive systems in multiple ways. In stochastic biochemical systems, such as reactions inside the cytoplasm, or in catalytic processes, some of the reacting molecules are present in very small numbers and therefore, discrete stochastic models are necessary to describe the system. In diffusion-limited reactive systems, such as simple coagulation 2*A* → *A*_{2} or annihilation *A* + *B* → *C*, spatial fluctuations in the concentration of the reactants grow as the reaction progresses and must be accounted for to accurately model the correct macroscopic behavior.^{2,3} In unstable systems, such as diffusion-driven Turing instabilities,^{4–8} detonation,^{1} or buoyancy-driven convective instabilities,^{9} fluctuations are responsible for initiating the instability and may profoundly affect its subsequent temporal and spatial development. In systems with a marginally stable manifold, fluctuations lead to a drift along this manifold that cannot be described by the traditional law of mass action, and has been suggested as being an important mechanism in the emergence of life.^{10–12}

Much of the work on modeling stochastic chemistry has been for homogeneous, “well-mixed” systems, such as continuously stirred tank reactors (CSTRs), but there is increasing interest in spatial models.^{13} When hydrodynamic transport is included, the focus has almost exclusively been on species diffusion, and there is a large body of literature on stochastic reaction-diffusion models. A master equation approach, notably, the Chemical Master Equation (CME), is widely accepted for modelling well-mixed systems. The Reaction-Diffusion Master Equation (RDME) extends this type of approach to spatially varying systems.^{14–16} In the RDME, the system is subdivided into reactive subvolumes (cells) and diffusion is modeled as a discrete random walk by particles moving between cells, while reactions are modeled using local CMEs.^{17,18} A large number of efficient and elaborate event-driven kinetic Monte Carlo algorithms for solving the CME and RDME, exactly or approximately, have been developed with many tracing their origins to the Stochastic Simulation Algorithm (SSA) of Gillespie.^{19,20} The issue of convergence as the RDME grid is refined is delicate.^{21–23} Variants of the RDME have been proposed that improve or eliminate the sensitivity of the results to the grid resolution, such as the convergent RDME (CRDME) of Isaacson^{23} in which reactions can happen between molecules in neighboring cells as well. Particle-based spatial methods for stochastic chemistry include reactive Brownian dynamics,^{24,25} Green’s function reaction dynamics,^{26,27} first-passage kinetic Monte Carlo,^{28–31} the small-voxel tracking algorithm,^{32} and others.

In large part, the development of stochastic reaction-diffusion models has been divorced from work in the fluid dynamics community. Full hydrodynamic transport including advection, sound waves, viscous stress, and thermal conduction, as well as nonequilibrium thermodynamics and chemistry, are fairly common in the reacting flow community. See, for example, textbooks by Kuo^{33} and Law.^{34} However, work in this area focuses on macroscopic modeling; spontaneous thermal fluctuations, either chemical or hydrodynamic, are typically not considered.

Within the field of nonequilibrium thermodynamics, the fluctuation-dissipation theorem provides the connection between hydrodynamic transport and spontaneous fluctuations. In particular, as an extension of conventional hydrodynamic theory, fluctuating hydrodynamics incorporates mesoscopic fluctuations in a fluid by adding stochastic flux terms to the deterministic fluid equations.^{35} These noise terms are white in space and time and are formulated using fluctuation-dissipation relations to yield equilibrium covariances of the fluctuations in agreement with equilibrium statistical mechanics. Linearized fluctuating hydrodynamics was first introduced by Landau and Lifshitz^{36} and has since been used to study simple and binary fluid systems in and out of equilibrium.^{35}

A number of numerical algorithms for solving the equations of fluctuating hydrodynamics have also been developed.^{37–44} These algorithms draw from a wealth of deterministic computational fluid dynamics (CFD) techniques and handle transport such as diffusion in a much more sophisticated fashion than random hopping between cells. For example, they include effects such as cross-diffusion, barodiffusion, thermodiffusion (i.e., Soret effect) as well as advection by fluid motion. Furthermore, semi-implicit temporal discretizations and higher-order spatial discretizations can be used, even as fluctuations due to the discrete nature of the fluid are accounted for. In this spirit, Koh and Blackwell^{45} propose using a more traditional gradient-driven diffusive flux formulation consistent with CFD practice within a CME-based description; however, their treatment of fluctuations is rather *ad hoc* and not consistent with the formulation of stochastic mass fluxes in fluctuating hydrodynamics. Very recently, a spatial chemical Langevin formulation (SCLE) was proposed^{46} in which the chemical Langevin approximation^{47} is applied to the RDME treating diffusive hops as another reaction in a very large reaction network. While this leads to a formulation similar to fluctuating hydrodynamics, it has several shortcomings, notably, it does not allow one to treat diffusion using advanced CFD algorithms.

Investigations utilizing fluctuating hydrodynamics have revealed the crucial importance of hydrodynamic fluctuations in transport mechanisms, especially mass and heat diffusion. Notably, it is now well-known that all nonequilibrium diffusive mixing processes are accompanied by long-range correlations of fluctuations. In certain scenarios, these nonequilibrium fluctuations grow in physical extent well beyond molecular scales with magnitudes far greater than those of equilibrium fluctuations. These so-called “giant fluctuations” are observed in laboratory experiments,^{48–50} and arise because of the coupling between thermal velocity fluctuations and concentration or temperature fluctuations. In fact, it has recently been shown using nonlinear fluctuating hydrodynamics that mass diffusion in liquids is dominated by advection by thermal velocity fluctuations.^{51} Therefore, modeling diffusion using collections of independent random walkers, as done in the RDME, is fundamentally inappropriate for describing the nature of hydrodynamic fluctuations at microscopic and mesoscopic scales; instead, hydrodynamic coupling (correlations) between the diffusing particles must be taken into account.^{52} Including fluctuations within the continuum description has also been shown to be important in particle-continuum hybrids,^{53} and should also benefit hybrid models for reaction-diffusion systems.^{54}

In the hydrodynamic equations, chemical reactions may be treated as a white noise source term,^{16,55} in a fashion analogous to the stochastic transport fluxes. The study of fluctuating hydrodynamic models that include chemical reactions is relatively recent and there are few computational studies in the literature. Stochastic reaction-diffusion equations are considered by Atzberger in Ref. 56, but only within the reaction-diffusion framework and not accounting for fluctuations in the chemical production rates. A thorough discussion of stochastic formulations of chemical reactions within the framework of statistical mechanics can be found in the monograph by Keizer;^{57} Keizer does not, however, consider hydrodynamic transport in spatially extended systems in depth. In a sequence of important papers,^{58–61} chemical reactions have been incorporated in a *nonlinear* nonequilibrium thermodynamic formalism, making it possible to combine realistic nonlinear deterministic models based on the traditional law of mass action (LMA) kinetics with fluctuating hydrodynamics. When considering fluctuations, however, a linearized approximation was used by the authors, limiting the range of applicability to modeling small Gaussian fluctuations around a macroscopic state that evolves in a manner unaffected by the fluctuations. A more phenomenological approach was followed to fit the LMA into the nonequilibrium thermodynamics GENERIC formalism by Grmela and Öttinger,^{62} but fluctuations were not considered. Here, we formulate a complete set of fluctuating hydrodynamic equations for a reactive multispecies mixture of ideal gases. We account for mass, momentum and energy transport, and chemical reactions, and consider a nonlinear formalism for describing the thermal fluctuations.

Hydrodynamics is a macroscopic coarse-grained description, and fluctuating hydrodynamics is a mesoscopic coarse-grained description. As such, both descriptions rely on the approximation that the length and time scales under consideration are much larger than molecular, i.e., that each coarse-grained degree of freedom involves an average over many molecules. In fact, although formally written as a continuum model, fluctuating hydrodynamics is, in truth, a discrete model that only makes sense when seen as a coarse-grained description for the evolution of a collection of spatially discrete hydrodynamic variables involving averages over many nearby molecules.^{63,64} The fact that many molecules are involved in the reactions allows for a Langevin-like continuum description (i.e., diffusion processes) of the fluctuations instead of discrete models such as master equations (i.e., jump processes). The accuracy of Langevin formulations for chemically reacting systems has long been a topic of debate.^{65,66} In this work, we take the first step in combining realistic fluid dynamics with a stochastic chemical description and adopt a Langevin approach to describing fluctuations. In future work, we will explore combining Langevin and CME approaches together, thus further bridging the apparent gap between the two.

Here, we first formulate the fluctuating reactive Navier-Stokes-Fourier equations, discuss their physical validity, and develop numerical methods for solving the stochastic partial differential equations. The methodology is a direct extension of our previous work on fluctuating hydrodynamics for non-reactive multispecies gas mixtures^{42} to include a Langevin model of chemical reactions. We consider two distinct Langevin models, which are identical when very close to chemical equilibrium but differ far from thermodynamic equilibrium. The first model, which we term the Log-Mean equation (LME), is based on the GENERIC formulation of Grmela and Öttinger,^{62} but can be traced to older work on the subject as well.^{57,67,68} The second Langevin model is the more familiar chemical Langevin equation (CLE).^{47,57,69}

The resulting algorithms are used to assess the importance of thermal fluctuations in several simple but relevant examples. The first example is a simple dimerization reaction, which has been studied theoretically in prior work by others.^{59–61} Our second example is the Baras-Pearson-Mansour (BPM) reaction network,^{70,71} which exhibits a rich behavior ranging from bistability to limit cycles. We study these examples in both well-mixed small-scale systems, comparing with the chemical master equation, and in spatially extended systems, comparing with fluctuating hydrodynamic theory and previous numerical work. For the latter, rather than imposing the non-equilibrium constraint by fixing concentrations in the bulk, the constraints are applied as boundary conditions, thus maintaining *strict* consistency with equilibrium thermodynamics, including microscopic reversibility (detailed balance), in all of the models we study. These examples illustrate how thermal fluctuations drive giant concentration fluctuations and how they affect the rate of pattern formation in an inhomogeneous system.

## II. THEORY

In this section, we summarize the mathematical formulation of the complete fluctuating Navier-Stokes (FNS) equations for compressible reactive multispecies fluid mixtures. The details for non-reactive fluid mixtures are presented in Ref. 42; here, we focus on the chemistry. The formulation is first presented in its general form; the specific case of reactions in ideal gas mixtures is treated in Sec. II C.

The species density, momentum, and energy equations of hydrodynamics for a mixture of *N*_{S} species are given by

where *ρ _{s}*,

*m*, and Ω

_{s}_{s}are the mass density, molecular mass, and number density production rate for species

*s*. The variables

**v**,

*p*, and

*E*denote, respectively, fluid velocity, pressure, and specific total energy for the mixture. The total density is $\rho = \u2211 s = 1 N S \rho s $, where

**g**is gravitational acceleration and superscript

*T*denotes transpose.

We consider a system with *N*_{R} elementary reactions with reaction *r* written in the form,

Here, 𝔐_{s} are the chemical symbols and $ \nu s i \xb1 $ are the molecule numbers for the forward and reverse reaction *r*. The stoichiometric coefficients are $ \nu s r = \nu s r \u2212 \u2212 \nu s r + $ and mass conservation requires that $ \u2211 s \nu s r m r =0$.^{57} For simplicity of notation, when there is no ambiguity we omit the range of the sums and write $ \u2211 s $ for sums over all species, and write $ \u2211 r $ for sums over all reactions. Note that chemistry does not appear explicitly in energy equation (3) since the species heat of formation is included in the specific total energy.

Transport properties are given in terms of the species diffusion flux, **ℱ**, viscous tensor, **Π**, and heat flux, **𝒬**. Mass conservation requires that the species diffusion flux and the production rate due to chemical reactions satisfy the constraints

so that summing the species equations gives the continuity equation

The detailed form of the transport terms is summarized in Appendix A, see Ref. 42 for details. It is important to note that we neglect any possible effect of the chemical reactions on the transport coefficients of the mixture.

We write the chemical production rate as the sum of a deterministic and a stochastic contribution, $ \Omega s = \Omega \xaf s + \Omega \u02dc s $, with the stochastic rate going to zero in the deterministic limit.^{72,73} To formulate these production rates, we define the dimensionless chemical affinity as

where *μ _{s}* is the specific chemical potential (i.e., per unit mass) and $ \mu \u02c6 s = m s \mu s / k B T$ is the dimensionless chemical potential per particle;

*T*and

*k*are temperature and Boltzmann’s constant, respectively. Summing over reactions gives the deterministic production rate for species

_{B}*s*as

^{62}

where

and *τ _{r}* is a time scale characterizing the reaction rate. This form of the deterministic equations, while at first sight appearing different from the more familiar LMA, is fully consistent with it. The production rate, as given by (7) and (8), is also consistent with nonequilibrium thermodynamics;

^{59}this way of expressing the production rate in terms of a thermodynamic driving force (difference of exponentials of chemical potentials) can be seen as a generalization of the LMA to non-ideal systems, as elaborated in Sec. II C.

For a binary mixture undergoing a dimerization reaction, the deterministic part of the complete set of hydrodynamic equations including chemical reactions has been fit into a *nonlinear* nonequilibrium thermodynamics formalism in Ref. 59 by introducing an additional reaction coordinate, as inspired by earlier work of Pagonabarraga *et al.*^{58} This extends earlier considerations of dimerization reactions in a strictly linear fluctuating chemistry formalism.^{74} In the limit of high reaction barrier, the equations written in Ref. 59 are equivalent to the ones we employ here even though our notation is different. However, fluctuating contributions in Ref. 59 are only considered in a linearized approximation, severely limiting the range of applicability to describing small Gaussian fluctuations around a deterministic average flow.

In Secs. II A and II B, we develop two *nonlinear* forms for the stochastic contribution to the reactive production rates, one coming from irreversible thermodynamics cast in the GENERIC formalism^{62} and the other being a generalization of the more familiar form associated with the CLE.^{47,57,69}

### A. The log-mean equation

Grmela and Öttinger^{62} cast phenomenological LMA (7) in the GENERIC formalism and obtain a nonlinear form for the dissipative matrix, under the *assumption* of a quadratic dissipative potential. Note that the entropy production rate can *uniquely* be written as a quadratic function of the thermodynamic driving force *only* for a single reaction; the resulting peculiar form of the mobility (dissipative) matrix (see Eq. (113) in Ref. 62) involving a logarithmic mean has recently been justified from a model reduction perspective.^{75} Here, we *assume* that there is no cross-coupling between different reactions and thus, associate an independent stochastic production rate with each reaction. Coupling between distinct reactions has been considered within a nonequilibrium thermodynamic framework only in some very specific cases^{76,77} and a general formulation requires more detailed knowledge about the coupling mechanism than is available in practice.

Following the general principles for including fluctuations^{78} in the GENERIC formalism,^{79} it is straightforward to write a Gaussian stochastic production rate assuming independence among the different reaction channels,

where

and where $ Z r \Omega $ are independent white-noise random scalar fields with covariance

with each $ Z r \Omega $ driving the stochastic production rate of a single chemical reaction *r*. We refer to this formulation for the stochastic chemistry as the “log-mean” form; the reasoning behind this name will become evident when presented in Sec. II C for ideal mixtures. Note that (9) uses the kinetic or Klimontovich interpretation^{80,81} of the stochastic integral, formally denoted as a kinetic stochastic product with a ♢ symbol in (9). The variance of the stochastic forcing $ D r LM $ can be seen to be positive because $ A \u02c6 $ and $A$ always have the same sign.^{62} Note that $ \u2211 s m s \Omega \xaf s = \u2211 s m s \Omega \u02dc s LM =0$, as required by mass conservation.

For the purposes of exposition, it is useful to consider a homogeneous “well-mixed” system of volume Δ*V*, which will correspond to a single hydrodynamic cell after spatial discretization of (1). The dynamics of the (intensive) number density *n _{s}* =

*N*/Δ

_{s}*V*, where

*N*= Δ

_{s}*Vρ*/

_{s}*m*is the number of molecules of species

_{s}*s*in the cell, is given by

which can also be written in Ito form as

where $ W r \Omega ( t ) $ are independent scalar white noise processes with covariance $\u3008 W r \Omega ( t ) W r \u2032 \Omega ( t \u2032 ) \u3009= \delta r , r \u2032 \delta ( t \u2212 t \u2032 ) $. We call this system of stochastic ordinary differential equations (SODEs) the LME.

The derivative in the last stochastic drift term in (12) is the directional derivative of $ D r LM $ along the reaction coordinate. Unlike the more familiar Ito or Stratonovich interpretations of the noise, the kinetic form of the noise ensures that the corresponding Fokker-Planck equation has the traditional form^{82}

This ensures that the LME is in detailed balance with respect to the Einstein distribution ∼*S*/(*k _{B}T*) for a closed system at thermodynamic equilibrium, where

*S*is the total entropy of the system.

^{79}We note that it is not possible to obtain the LME from the CME with a systematic procedure; one must invoke some guiding principles about the structure of coarse-grained Fokker-Planck equations to “derive” this form of the noise.

^{67,68,82}

### B. The chemical Langevin equation

Since both $A$ and $ A \u02c6 $ are equal to zero at chemical equilibrium, near chemical equilibrium we can linearize (10) to first order in the affinity $A$, and approximate the amplitude of the stochastic production rate in terms of a sum over each forward and reverse reaction, that is,

Since this sum of products of exponentials is evidently positive, we can potentially use it even far from chemical equilibrium, and write the stochastic production rate as

where^{83}

Here, $ Z r , + \Omega $ are independent white-noise scalar random fields that give the stochastic contribution from the forward reaction, while $ Z r , \u2212 \Omega $ correspond to the reverse reactions; the forward and reverse reactions are taken to be independent. In Sec. II C, the production rate factors, $ D r LM $ and $ D r CL $, are further simplified for the case of ideal gas mixtures.

Form (15) for the amplitude of the stochastic production rate is found in most work on the subject.^{47,59–61,84} For example, though not written in this form, Eq. (8f) in Ref. 61 contains a sum of two exponential terms and is equivalent to (14) for the specific reaction considered there.^{85} For a well-mixed homogeneous system of volume Δ*V*, the number densities of molecules of the different species follow the system of SODEs,

Stochastic equation (17) is commonly referred to as the CLE following Gillespie,^{47} and can be obtained from the CME by an uncontrolled truncation of the Kramers-Moyal expansion to second order. It is traditional to assume an Ito interpretation of the noise in the CLE, even though no precise justification for this can be made within the accuracy to which the CLE approximates the CME.^{69} Mathematically, the nonlinear CLE contains similar information to the central limit theorem (i.e., linearized fluctuating hydrodynamics) corresponding to the CME in the limit of weak noise (large number of reactant molecules).

As seen from (14), the two stochastic differential equations for the number densities, LME using kinetic noise (12) and CLE using Ito noise (17), are equivalent near chemical equilibrium. They are, however, different far from chemical equilibrium, as we illustrate in more detail in Sec. IV A. Notably, the forward and reverse reactions are treated together in the LME, consistent with the fact that, due to microscopic reversibility, there is only one independent rate coefficient for each reaction.^{57} The ratio of the forward and reverse reaction rates is related to the equilibrium reaction constant, which is a *thermodynamic* and not a kinetic quantity. In fact, the LME is closely related to the notion of the existence of a state of thermodynamic equilibrium in which each pair of forward and reverse reactions is in detailed balance with each other; one cannot write a LME for a system with irreversible reactions, which fundamentally violate detailed balance.

By contrast, the forward and reverse reactions are treated completely independently in the CLE and there is no difficulty in writing a CLE for a system with irreversible reactions. The CLE is evidently inconsistent with the notion of detailed balance and is, in fact, inconsistent with equilibrium thermodynamics. Although written in a different form, Keizer’s (4.8.37) is the CLE, and Keizer’s (4.8.36) is the LME;^{57} Keizer argues that the CLE is the correct equation and concludes: “Although the theoretical description of nonequilibrium ensembles would be greatly simplified if the phenomenological choice [LME] were correct, this appears not to be the case.” We will compare and contrast these two equations on some specific examples in Sec. IV A.

### C. The law of mass action and ideal gas mixtures

In the formulation of hydrodynamic transport one normally works with the specific chemical potential, which has the general form^{86}

where $ \mu s o $ is the chemical potential at a reference state, $ X s = N s / \u2211 s \u2032 N s \u2032 $ is the mole fraction, and *γ _{s}*(

*p*,

*T*,

**) is the activity coefficient of species**

*X**s*. For chemistry, it is more convenient to work with a dimensionless chemical potential per particle,

where $ \mu \u02c6 s o = ( m s \mu s o ) / ( k B T ) $. Note that *X _{s}γ_{s}* is the activity (i.e., effective concentration) and for an ideal mixture,

*γ*= 1.

_{s}^{87}This gives

which leads to a generalized LMA of the form

where $ \kappa r \xb1 ( T , p , X ) $ are the more familiar forward/reverse reaction rates (per unit time and per unit volume). Since there is only one independent time scale parameter, *τ _{r}*, the forward and reverse rates are not independent and the LMA gives the ratio to be the equilibrium constant,

which is a purely thermodynamic quantity (closely related to the dimensionless reference Gibbs energy for the reaction at a unit reference pressure) that can be calculated from pure component data.^{84,88} Note that in chemistry texts, the equilibrium constant is typically defined in terms of concentrations rather than activities as we have done here.

For ideal gas mixtures, we can further simplify generalized LMA (18) to the more familiar form using number densities instead of mole fractions. From classical statistical mechanics, for an ideal gas mixture we can write^{89}

where $ \Lambda s =h/ 2 \pi m s k B T $ is the thermal wavelength of a structure-less particle and *j*(*T*) is the partition function for the internal degrees of freedom.^{90} In general, *j*(*T*) is a complicated function depending on the quantized energy levels of a molecule but in the classical approximation $j ( T ) = ( T / T o ) 1 2 z $, where *z* is the number of classical internal degrees of freedom and *T _{o}* is a reference temperature.

For ideal gas mixtures, chemical production rate (18) can be written in the familiar power-law form,

where $ k r \xb1 $ are the forward/reverse reaction rates for the LMA formulated in terms of number density instead of activity. For unimolecular reactions (e.g., 𝔐_{1} → ⋯), the “decay time” for a particle is usually assumed to be constant, in which case the corresponding reaction rate (e.g., $ k r + $) is a constant. For bimolecular reactions (e.g., 𝔐_{1} + 𝔐_{2} → ⋯) the production rate is usually assumed to be proportional to the collision frequency times an Arrhenius factor, in which case the corresponding reaction rate is only a function of temperature.^{91}

For the stochastic production rate in an ideal gas mixture, using

gives

where logmean and arthmean are the logarithmic and arithmetic mean, respectively. Note that $ D r LM $ is zero if either reaction rate equals zero while $ D r CL $ is non-zero if either reaction rate is non-zero. The logarithmic mean form of the noise in the ideal mixture case has appeared, phenomenologically, in several early papers^{67,68} and in more recent work^{75,92,93} by other authors.

## III. NUMERICAL SCHEME

The numerical integration of (1)-(3) is based on a method of lines approach in which we discretize the equations in space and then use a SODE integration algorithm to advance the solution using the basic overall approach described in Ref. 42. The spatial discretization uses a finite volume representation with cell volume Δ*V*, where $ U i , j , k n $ denotes the average value of **U** = (*ρ _{s}*,

*ρ*

**v**,

*ρE*) in cell-(

*i*,

*j*,

*k*) at time step

*n*. To ensure that the algorithm satisfies discrete fluctuation-dissipation balance, the spatial discretizations for the hydrodynamic fluxes are done using centered discretizations; see Refs. 40 and 42 for details.

Discretization of the system in space results in a system of SODEs driven by a collection of independent white-noise processes **𝒲**(*t*) that represent a spatial discretization of the random Gaussian fields **𝒵**(** r**,

*t*) used to construct the noise. After temporal discretization, these white noise processes are represented by a collection of independent and identically distributed (i.i.d.) standard normal variates

**, which can be thought of as a spatio-temporal discretization of**

*Z***𝒵**; the discretization is reflected in the presence of a prefactor $1/ \Delta V \Delta t $ in the expressions for $ \Omega \u02dc $ given below.

^{94}

For temporal integration, we use the low-storage third-order Runge-Kutta (RK3) scheme previously used to solve the single and two-component FNS equations,^{40} using the weighting of the stochastic forcing proposed by Delong *et al.*^{94} With this choice of weights, the temporal integration is weakly second-order accurate for additive noise (e.g., the linearized equations of fluctuating hydrodynamics^{95}). As discussed at length in Ref. 95, the hydrodynamic stochastic fluxes should be considered as additive noise in a linearized approximation.

The implementation of the methodology supports three boundary conditions in addition to periodicity. The first is a specular, adiabatic wall which is impermeable to momentum, mass, or heat transport (i.e., all fluxes are zero at the wall). A second type of boundary condition is a no slip, reservoir wall at which the normal velocity vanishes (i.e., the total mass flux at the boundary vanishes) and the other velocity components, mole fractions and temperature satisfy inhomogeneous Dirichlet boundary conditions; this mimics a permeable membrane connected to a reservoir on the other side of the boundary. The third boundary condition is a variant of the no slip condition for which the wall is impermeable to mass but conducts heat. When a Dirichlet condition is specified for a given quantity, the corresponding diffusive flux is computed as a difference of the cell-center value and the value on the boundary. In such cases, the corresponding stochastic flux is multiplied by $ 2 $ to ensure discrete fluctuation-dissipation balance, as explained in detail.^{43,96}

Because the noise arising from the chemical reactions is multiplicative, special care must be taken to capture the correct stochastic drift terms arising from the kinetic interpretation of the noise in the LME.^{97,98} We have chosen to write the equations in Ito form, which leads to an additional stochastic drift term in LME (12). To integrate the Ito form in time, we evaluate the amplitude of the noise at the beginning of the time step and reuse the same random increments in all three stages of the RK3 scheme. The stochastic drift term arising in the LME is treated as a deterministic term but is also only evaluated at the beginning of the time step. The resulting scheme is only first-order weakly accurate. It is possible to construct second-order weak integrators by using the special one-dimensional nature (i.e., there is only a single reaction coordinate for each reaction even if there are many species involved) of the chemical noise.^{99} However, in our simulations the time step is typically limited by stability considerations for advective and diffusive hydrodynamic processes, and therefore chemistry is accurately resolved even by a first-order scheme. Alternative temporal integration strategies will be discussed in the Conclusions.

The chemical Langevin form of the noise in (17) is discretized as

where $ Z r \Omega i , j , k $ are zero-mean normal Gaussian variates generated independently in each cell at the beginning of each time step, $\u3008 ( Z r \Omega ) i , j , k n ( Z r \u2032 \Omega ) i \u2032 , j \u2032 , k \u2032 n \u2032 \u3009= \delta i , i \u2032 \delta j , j \u2032 \delta k , k \u2032 \delta n , n \u2032 \delta r , r \u2032 $. Note that the two terms in (15) have been combined into a single white-noise process with amplitude $ D r CL = D r , + CL + D r , \u2212 CL $. As discussed above, LME noise (9) leads to an Ito correction in (12),

The directional derivative of $ D r LM $ in the last term can be evaluated analytically, or, for simplicity of implementation, it can be efficiently approximated numerically using a finite difference along the reaction coordinate.

## IV. NUMERICAL RESULTS

In this section we describe several test problems that demonstrate the capabilities of the numerical methodology. We consider two reaction systems, the first being simple dimerization,

where 𝔐 = (*A*, *A*_{2}), that is, species 1 is the monomer and species 2 is the dimer. In Sec. IV A 1 we investigate simple dimerization in a homogeneous system; in Sec. IV B we investigate the “giant fluctuation” phenomenon in the presence of dimerization for a system with an applied concentration gradient.

The second model we consider is based on the Gray-Scott (GS) model,^{100,101} which is known to exhibit a rich morphology of stationary and time-dependent patterns.^{102} This model, as formulated by Pearson,^{102} consists of the reactions

where the concentrations of the “feed species,” *U _{f}* and

*V*, are held fixed and species

_{f}*S*is inert. Since elementary reactions are rarely trimolecular in nature, we consider a variation of the GS model developed by Baras

*et al.*

^{70,71}The BPM model

^{103}is

with 𝔐 = (*U*, *V*, *W*, *S*, *U _{f}*,

*V*). This model was developed as a more realistic variant of the Gray-Scott model suitable for particle simulations of dilute gases.

_{f}^{104}Note that the first and third reactions are irreversible in the original BPM model, which is not consistent with detailed balance. The BPM model has not been studied as extensively as the GS model but its dynamics are expected to be qualitatively similar. Note that the BPM model replaces the trimolecular reaction in the GS model with a pair of bimolecular reactions and introduces

*W*as an intermediary species. In the standard BPM model,

^{70}the number densities of

*S*,

*U*, and

_{f}*V*are held fixed so, being an open system, total mass is not conserved and detailed balance is not satisfied.

_{f}^{57}In Sec. IV A 2 we investigate this standard BPM model in a homogeneous “well-mixed” system; in Sec. IV C we simulate a two dimensional domain with full hydrodynamic transport with species

*S*,

*U*, and

_{f}*V*held fixed only at the boundaries.

_{f}### A. Homogeneous systems

We first consider homogeneous “well-mixed”^{105} systems of volume Δ*V* with only chemistry (i.e., no hydrodynamics). In this section we compare the results obtained using the LME form, (22), and the CLE form, (23), with results from CME simulations performed using the SSA, also known as the Gillespie algorithm.^{19} The CME is widely accepted as an accurate model for well-mixed chemical systems and SSA is a popular scheme for simulating the stochastic process described by the CME.^{20} As we will see, the two forms for the Langevin noise have their advantages and disadvantages and both forms are only approximations of the CME with limited ranges of validity.

#### 1. Dimerization reaction

We start by considering dimerization reaction (25) in a closed system for which the deterministic production rate for species 1 (monomers) is

and by mass conservation, for dimers $ \Omega 2 =\u2212 1 2 \Omega 1 $ since *m*_{2} = 2*m*_{1}. The constraint of mass conservation can also be expressed as *n*_{1} + 2*n*_{2} = *n*_{0}, where *n*_{0} is the initial number density of A particles (in either monomer or dimer form). We may then write

and limit our attention to the monomer species. For simplicity, we take the ratio of the rate constants to be *k*^{+}/*k*^{−} = 1/*n*_{0} so the equilibrium mass fraction $ ( Y 1 ) eq = 1 2 $ (i.e., $ \Omega \xaf 1 =0$ for $ Y 1 = 1 2 $).

The log-mean stochastic production rate for *n*_{1} is

with

From corresponding Fokker-Planck equation (FPE) (13) one finds the LME is in detailed balance with respect to the equilibrium distribution

where *Z* is a normalization constant. This Einstein distribution $ P eq LM ( Y 1 ) \u223cexp ( S ( Y 1 ) / k B ) $ is in agreement with the correct thermodynamic entropy *S* in the limit of Stirling’s approximation, as we demonstrate in Appendix B.

For the chemical Langevin equation, the stochastic production rate for *n*_{1} is

with

The equilibrium distribution can also be found from the stationary solution of the FPE corresponding to the CLE, which we do not write here for brevity.^{106} We do note that, for this example, $ P eq CL ( Y 1 ) $ is quite close to a Gaussian. We further observe that, unlike the LME, no thermodynamic interpretation can be given to $ln\u2009 P eq CL $. In fact, the tails of $ P eq CL $ are quite different from those of $ P eq LM $ and, being nearly Gaussian, the former includes unphysical values of the concentration (i.e., *Y*_{1} is not constrained between 0 and 1).

Figure 1 shows numerical results for the equilibrium distribution of the number of monomers *N*_{1} = *n*_{1}Δ*V* = *ρY*_{1}Δ*V*/*m*_{1}. At thermal equilibrium, the simulation results using the LME form for the noise are in excellent agreement with equilibrium statistical mechanics (see Appendix B) and with CME/SSA simulation results. Other work has also shown that, when detailed balance is obeyed, the LME correctly reproduces the equilibrium transition rates for rare jumps between stable minima in bistable systems.^{68} On the other hand, the CLE result has the noticeable flaw that, being a Gaussian, the distribution extends to unphysical negative values of *N*_{1}.

However, the LME does not compare favorably with the numerical solution of the CME for time-dependent situations, such as when a system is relaxing toward equilibrium. To illustrate this, we simulate an ensemble of systems prepared with an initial condition far from equilibrium, specifically with *Y*_{1}(*t* = 0) ≈ 1, and measure the time-dependent probability distribution *P*(*Y*_{1}, *t*) as the system relaxes toward chemical equilibrium ($ ( Y 1 ) eq = 1 2 $). As expected, for the ensemble mean value of the number density $ n \u0304 1 ( t ) =\u3008 n 1 ( t ) \u3009$, we find close agreement among LME, CLE, and CME results (not shown), even when fluctuations are quite large. However, if we consider the standard deviation of the number of monomers, the left panel in Fig. 2 clearly demonstrates that the CLE is in much better agreement with the CME (as shown by the SSA results) for describing relaxation toward equilibrium. Also shown on this graph is the theoretical solution for the standard deviation obtained by first linearizing the CLE^{107} around the solution of the deterministic law of mass action (which is the law of large numbers corresponding to the CME^{69}), and then writing a system of ODEs for the mean and variance of *n*_{1}(*t*). Specifically, we have that $d n \u0304 1 ( t ) /dt= \Omega 1 $ and, using Ito’s formula, we get the central limit theorem corresponding to the CME,^{69}

where $ C 1 ( t ) =\u3008 n 1 ( t ) \u2212 n \u0304 1 ( t ) 2 \u3009$.

The agreement between CLE and CME in the left panel of Fig. 2 is not surprising since it is well-known that the central limit theorem for the CME is a linear Langevin equation with the noise covariance given by the CLE form rather than the LME form.^{69,108} The agreement between the CLE and SSA results is less impressive when we look more closely at the probability distribution during the relaxation toward equilibrium. The right panel of Fig. 2 shows histograms of the probability distribution for the monomers, *P*(*N*_{1}, *t*), at an early time in the relaxation. This distribution is *not* close to Gaussian for the CME/SSA solution and we see that, in this regard, the CLE result is no better than the LME result. In fact, for the probability distribution of the dimers (not shown) the CLE results have the unphysical feature of non-zero probability for negative values of dimer concentration. Of course, one can argue that the number of molecules in the system is too small for a Langevin approximation to apply. If the fluctuations are decreased, the probability distribution *P*(*N*_{1}, *t*) will become closer to Gaussian and then the CLE will provide a better description; note, however, that the tails of the distribution will always be incorrect for the CLE, even at thermodynamic equilibrium.

#### 2. Bistable BPM model

As a less trivial homogeneous example, we consider BPM model (27) for an open system held at a non-equilibrium steady state for which the probability distribution function is bimodal.^{71} For the chemistry-only study in this section, the number of molecules of species 1, 2, and 3 (*U*, *V*, and *W*) is allowed to vary while the number of molecules of all other species is fixed. The relevant parameters are given in Table I. Note that a similar system (with all reactions being bimolecular) was studied by Baras *et al.*, who found good agreement between SSA and molecular simulations using the Direct Simulation Monte Carlo (DSMC) algorithm.^{71} Baras *et al.* also examined the accuracy of the CLE *linearized* around the solution of the deterministic equations, and, not surprisingly, found it to be a very poor approximation of the CME for the parameters they chose. Gillespie^{47} suggests that “A repetition of the study of Baras and co-workers using the Langevin equation [CLE] instead of the [linearized CLE] should show the chemical Langevin equation in a fairer light.” This section presents such a study using both the CLE and the LME. The parameters were selected such that the number of particles is large (roughly *O*(10^{2} − 10^{3}) for each species) but not so plentiful as to prevent the SSA simulation from accurately sampling the bimodal distribution in a reasonable amount of computation time.

Species . | S
. | U
. _{f} | V
. _{f} |
---|---|---|---|

n (fixed) _{r} | 7.0 × 10^{2} | 5 × 10^{11} | 5.65 × 10^{10} |

Species . | S
. | U
. _{f} | V
. _{f} |
---|---|---|---|

n (fixed) _{r} | 7.0 × 10^{2} | 5 × 10^{11} | 5.65 × 10^{10} |

… . | ℜ_{1}
. | ℜ_{2}
. | ℜ_{3}
. | ℜ_{4}
. | ℜ_{5}
. |
---|---|---|---|---|---|

k^{+} | 2 × 10^{−3} | 10^{−3} | 0.020 093 6 | 0.28 | 0.28 |

k^{−} | 2 × 10^{−9} | 2.198 × 10^{−4} | 2.009 × 10^{−8} | 2.8 × 10^{−10} | 2.8 ×10^{−10} |

… . | ℜ_{1}
. | ℜ_{2}
. | ℜ_{3}
. | ℜ_{4}
. | ℜ_{5}
. |
---|---|---|---|---|---|

k^{+} | 2 × 10^{−3} | 10^{−3} | 0.020 093 6 | 0.28 | 0.28 |

k^{−} | 2 × 10^{−9} | 2.198 × 10^{−4} | 2.009 × 10^{−8} | 2.8 × 10^{−10} | 2.8 ×10^{−10} |

A phase-space picture of a typical trajectory is shown in the left panel of Fig. 3. The trajectory moves between two basins centered around the two stable deterministic steady states, which are labeled state *A* (corresponding to *N*_{1} ≈ 1740, *N*_{2} ≈ 448, *N*_{3} ≈ 328 molecules) and state *B* (corresponds to *N*_{1} ≈ 1224, *N*_{2} ≈ 936, *N*_{3} ≈ 1424 molecules). Based on this picture, we chose to define a collective coordinate *x*(*t*) which is the projection of the state (*n*_{1}, *n*_{2}, *n*_{3}) onto the line connecting the two stable points (red line in the figure). This simple linear collective coordinate has the property that *x* = 0 at state *A* and *x* = 1 at state *B*; note that *x*(*t*) is *not* bounded between zero and one. The inset in Fig. 3 (left panel) shows *x*(*t*) for a typical trajectory as the system moves between the basins.

The right panel in Fig. 3 shows the steady-state probability distribution *P*(*x*) for the collective coordinate *x*, clearly illustrating its bimodal form; similar results are found for the probability distributions for *n*_{1}, *n*_{2}, and *n*_{3}. The results for the LME and CLE Langevin approximations are qualitatively similar to those from the CME/SSA but quantitatively different; the LME result is in better agreement with the CME for this specific example. To also examine the long-time dynamics of the well-mixed bistable BPM system, we assign each (discrete) point in time to either state *A* or state *B* (see inset in left panel of Fig. 3). The assignment is performed by defining two sets $A= x < 0 . 3 $ and $B= x > 0 . 6 $ and assigning each point in the trajectory to the *last* set that was visited. The distribution of waiting times spent in the two states before transitioning to another state is related to the transition rate, and the ratio of the average waiting times gives the ratio of the probabilities to be found in each of the two states. For large Δ*V* (weak noise), the transitions are rare events and the distribution of waiting times should be approximately exponential (recall that for an exponential distribution the variance is the square of the mean). Numerical results for the mean and variance of the time spent in state B before transiting to state A, and vice versa, are shown in Table II.

State A . | Mean . | Variance . | State B . | Mean . | Variance . |
---|---|---|---|---|---|

SSA | 1.13 × 10^{5} | 1.53 × 10^{10} | SSA | 2.86 × 10^{5} | 8.74 × 10^{10} |

LME | 3.23 × 10^{5} | 9.03 × 10^{10} | LME | 8.45 × 10^{5} | 1.91 × 10^{11} |

CLE | 1.38 × 10^{5} | 1.70 × 10^{10} | CLE | 5.01 × 10^{5} | 7.91 × 10^{10} |

State A . | Mean . | Variance . | State B . | Mean . | Variance . |
---|---|---|---|---|---|

SSA | 1.13 × 10^{5} | 1.53 × 10^{10} | SSA | 2.86 × 10^{5} | 8.74 × 10^{10} |

LME | 3.23 × 10^{5} | 9.03 × 10^{10} | LME | 8.45 × 10^{5} | 1.91 × 10^{11} |

CLE | 1.38 × 10^{5} | 1.70 × 10^{10} | CLE | 5.01 × 10^{5} | 7.91 × 10^{10} |

Our results indicate that for the BPM model the CLE and LME provide a reasonably good qualitative description of the long-time dynamics and rare-event statistics for the parameters studied here. However, both approximations are in general uncontrolled and the quantitative match between the CME and either CLE or LME will not improve even if the cell volume increases and the fluctuations decrease in amplitude. In Ref. 68 it is observed that the LME correctly reproduces the very long-time dynamics (more precisely, the large deviation theory) of the CME for the bistable Schlogl model. This conclusion is, however, specific to this simple one-dimensional model because the system obeys detailed balance even though it is not in thermodynamic equilibrium; the BPM model studied here is not in detailed balance and there is no *a priori* reason to expect the LME to be more accurate than the CLE. The fact that the CLE is not able to describe rare events is well-known, see, for example, the discussion by Gillespie in Ref. 109 and after Eq. (9b) (which is the CLE) in Ref. 66, or recent numerical studies of noise-induced multistability.^{110} It can, in fact, easily be proven that this problem is shared by *all* diffusion process (SODE) approximations of the CME,^{111} and fundamentally stems from the difference between the rare-event statistics of Gaussian and Poisson noise.^{112} A promising alternative is to use tau-leaping to approximately integrate the CME in time^{66} since it uses Poisson noise, and thus has the potential to correctly approximate the long-time behavior of the CME. This point is discussed further in the Conclusions.

### B. Giant fluctuations

We now consider a system for which concentration fluctuations are affected by both chemistry and hydrodynamics in an interesting fashion. In the absence of chemistry, a gradient of concentration induces a long-ranged correlation of concentration fluctuations.^{35,113,114} These correlations are closely tied to the experimentally observed “giant fluctuation” phenomenon.^{48–50} In an isothermal, nonreacting binary mixture the static structure factor for fluctuations in the mass fraction of the first species contains two contributions,

where “hat” denotes a Fourier component; the equilibrium part is

The non-equilibrium enhancement of the static structure factor due to a concentration gradient is $ S neq k \u223c \u2207 Y 1 2 / k 4 $, where the wavevector ** k** is perpendicular to the imposed concentration gradient.

The nature of these long-ranged correlations is modified in the presence of chemical reactions, as predicted by linearized fluctuating hydrodynamics.^{60,61,115} Some preliminary numerical studies of fluctuations in the presence of chemistry have been done in Ref. 116 using a RDME-based description. However, these simulations are for a simpler isomerization *A*⇄*B* in one dimension and, furthermore, they are concerned with reaction-diffusion only and do not account for the hydrodynamic velocity fluctuations that are responsible for the giant concentration fluctuation phenomenon.

We consider here dimerization reaction (25) in a spatially inhomogeneous system. A rather detailed linearized fluctuating hydrodynamic theory for this example has been developed by Bedeaux *et al.* in Ref. 61, for a system in which a concentration gradient is imposed via a temperature gradient through the Soret effect. However, this analysis assumes a liquid mixture (large Schmidt and Lewis numbers) and thus does not apply to gas mixtures. Therefore, a simplified theoretical analysis of giant fluctuations in binary gas mixtures in the presence of an imposed constant concentration gradient and reactions is developed in Appendix C.

Our simplified theory decouples the temperature equation and uses a concentration equation (specifically the mass fraction of the first species) coupled to an incompressible fluctuating velocity equation. For the case of a *liquid* mixture with very large Schmidt number, which is the case considered in Ref. 61, the calculation predicts that the nonequilibrium enhancement of the static structure factor of concentration fluctuations for the monomer species is (see Eq. (C1))

where *D* is the diffusion coefficient and *η* is the viscosity. The last term on the r.h.s. depends on the penetration depth *d*,^{117}

We see that for large wavenumbers (*k* ≫ 1/*d*) the spectrum is ∼*k*^{−4}, as in the absence of the chemical reaction. However, for small wavenumbers (*k* ≪ 1/*d*) there is a transition to a *k*^{−2} spectrum. For gas mixtures, however, a more detailed model is required that takes into account the finite value of the Schmidt number. The result of this calculation is Eq. (C2), which predicts a further transition to a flat (constant) spectrum at small wavenumbers, with a finite $ S neq ( k = 0 ) = ( k B T ) \u2207 Y 1 2 / ( 9 \rho ( k \u2212 ) 2 ) $. The calculation in Appendix C indicates that this effect is important even in liquids and the more refined theory ought to be used if quantitative agreement with experiments or simulations is sought.

We performed a series of simulations to investigate these predictions using the full hydrodynamic equations plus chemistry. It is important to note that even though we use the full nonlinear equations, nonlinearities in the fluctuations are negligible for the simulations reported here.^{118} In fact, the noise is very weak (since the domain is quite thick in the *z* direction) and the numerical method is effectively performing a computational linearization of the fluctuating equations around the solution of the (nonlinear) deterministic equations;^{95} this is similar to what is done analytically in Refs. 60 and 61 but does not require any approximations. In the small Gaussian noise regime the linearized CLE equation applies, and therefore in these simulations we use the CLE form for the stochastic chemical production rate in agreement with the theory in Refs. 60 and 61. Identical results (to within statistical errors) are obtained by using the LME form of the noise (not shown); this is not unexpected since the important noise here is the stochastic momentum tensor driving the velocity fluctuations; the stochastic mass flux and production rates only affect the reaction-diffusion part of the spectrum, which is much smaller than the nonequilibrium enhancement we study here.

Here, we assume that traditional number-density based LMA (21) holds with constant rates *k*^{+} and *k*^{−}. This requires that the time scale for the reaction is proportional to the number density, that is, *τ* ∼ *p*/*k _{B}T* =

*n*. From (19) and (20), for the dimerization of an ideal gas the ratio of these rates is

which is a complicated function that depends on the form of the internal degrees of excitation. These details determine the number fraction (*X*_{1})_{eq} or, equivalently, the mass fraction (*Y*_{1})_{eq} at chemical equilibrium; here we set the ratio of the forward and reverse rates to ensure (*Y*_{1})_{eq} = 1/2, assuming that Λ_{1}, *j*_{1}, and *j*_{2} are consistent with this choice. Since the reaction here changes the number density and thus the pressure, the reaction is strongly coupled to the momentum and energy transport equations. In order to minimize this coupling, we adjust the number of internal degrees of freedom of the dimer particles. Specifically, we set the heat capacities to $ c p , 1 = 5 2 k B / m 1 $ (corresponding to three translational and zero internal degrees of freedom) and *c*_{p,2} = 5*k _{B}*/

*m*

_{2}(corresponding to three translational and five internal degrees of freedom). This choice ensures that

*c*of the mixture is independent of composition so that at constant pressure the reaction is isothermal.

_{p}The fluid was taken to be a dilute binary mixture of hard-sphere gases, using kinetic theory formulae for the transport coefficients.^{39} In CGS units, the species diameters are *σ*_{1} = 2.58 × 10^{−8} and *σ*_{2} = 3.23 × 10^{−8}, and *m*_{1} = 6.64 × 10^{−23}. At equilibrium the density *ρ*_{eq} = 1.78 × 10^{−3}, the temperature *T*_{eq} = 300, and the concentration (*Y*_{1})_{eq} = 0.5. For these parameters, at the equilibrium conditions the mass diffusion coefficient is *D* = 0.2698 while the momentum diffusion coefficient (kinematic viscosity) is *ν* = 0.2374. The ratio of the reverse and forward reaction rates is fixed at *k*^{−}/*k*^{+} = *ρ*/*m*_{1} = 2.679 85 × 10^{19}. We vary the penetration depth *d* by changing the value of the reaction rates.

The simulations used a 128^{2} grid and time step size Δ*t* = 2.5 × 10^{−8}, grid spacing Δ*x* = Δ*y* = 10^{−3}, and thickness in the *z* direction of Δ*z* = 10^{−3}. The first 6 × 10^{4} time steps were skipped and then statistics collected for 2 × 10^{6} time steps. A steady concentration gradient was imposed by using Dirichlet boundary conditions at top and bottom boundaries (*y* = *L* and *y* = 0). Specifically, we take $ Y 1 y = 0 , t =0.3$ and $ Y 1 y = L y , t =0.7$, with temperature fixed at *T* = 300 and no-slip boundary conditions for the velocity. Periodic boundary conditions were used in the other direction. Concentration profiles for various values of penetration depth, *d*, are shown in Fig. 4. As expected, when the chemistry is slow (*d* ≫ *L*) the concentration profile is nearly linear; when the chemistry is fast the concentration is nearly constant (at its chemical equilibrium value of (*Y*_{1})_{eq} = 1/2) except near the boundaries. Note that we set the thermal diffusion ratio to zero (i.e., no Soret effect) so that the system is isothermal and the simple theory presented in Appendix C applies.

For comparison, we also performed simulations in which we turn off all hydrodynamics except Fickian diffusion, giving us a reference reaction-diffusion structure factor *S _{rd}*(

*k*). For the case of a binary mixture

^{39}with

*k*

^{−}/

*k*

^{+}=

*ρ*/

*m*

_{1}, the reaction-diffusion CLE reduces to

where *n*_{0} = *ρ*_{1}/*m*_{1} + 2*ρ*_{2}/*m*_{2} is the total number density of A particles contained in both monomers and dimers, which is spatially constant in this reaction-diffusion approximation. It is important to note that reaction-diffusion model (33) is thermodynamically inconsistent because it ignores the coupling of the chemistry to the energy and momentum transport. It is well-known that adding chemistry should not change the fluctuations at thermodynamic equilibrium,^{60} and this is indeed the case for the complete set of hydrodynamic equations that we study in this work. By contrast, for reaction-diffusion (33) the equilibrium structure factor for (*Y*_{1})_{eq} = 1/2 is given by

which only approaches thermodynamically correct answer (31) for *kd* ≫ 1. Note that the inconsistency between full hydrodynamics and reaction-diffusion is not evident in Eqs. (27a) and (28) in Ref. 60, because the authors of that work “neglect the dependence of the specific Gibbs energy difference on pressure.” This inconsistency is not of any importance in our study because we only use the reaction-diffusion simulations to obtain a baseline to subtract from the full hydrodynamic runs at large wavenumbers; at small wavenumbers the nonequilibrium enhancement is many orders of magnitude larger than the difference between (31) and (34).

The reaction-diffusion runs are not limited by the Courant condition so we increased the time step to Δ*t* = 2.5 × 10^{−7}; a total of 6 × 10^{4} steps were skipped initially and then statistics were collected for 3 × 10^{6} steps. As seen in Fig. 4, the average concentration profiles are nearly the same whether the simulations used the full hydrodynamic equations or simply species diffusion. For the structure factor, however, we find that the reaction-diffusion simulations do *not* reproduce giant fluctuation result (31), rather, they follow (34), which does not show a power-law enhancement, as seen in the top panel of Fig. 5. This is expected since the giant fluctuation effect arises from the coupling of concentration fluctuations with the velocity fluctuations.

Because equilibrium structure factor (30) is derived assuming a uniform bulk state, which is not actually the case here, we define *S*_{neq}(*k*) = *S*(*k*) − *S _{rd}*(

*k*) as a measure of the “giant” nonequilibrium fluctuations coming from the coupling with the velocity equation. Results for

*S*

_{neq}(

*k*) for several penetration depths are shown in the bottom panel of Fig. 5, comparing simulation results with the simple theory, Eq. (C2). Since chemistry should have minimal effects for large

*k*according to the theory, we compute an effective concentration gradient by approximately matching the tail of the numerical result to the tail of the theory. We see that the theory correctly reproduces the qualitative trends, namely, that the giant fluctuations level off to a constant value at a wavenumber of order

*d*

^{−1}. However, except for the case of no reaction (

*d*→ ∞),

^{119}the theory is not in quantitative agreement with the simulations. To confirm that the issue is not under-resolution of the penetration depth by the grid, we perform runs with a finer grid of 256

^{2}cells,

^{120}and we get the same result over the common range of wavenumbers, showing these runs are sufficiently resolved for the purpose of computing

*S*(

*k*). Note that in the plots the numerical wavenumber

*k*is corrected to account for discretization artifacts in the standard 5-point Laplacian,

_{x}*k*

^{2}= sin

^{2}(

*k*Δ

_{x}*x*/2)/(

*k*Δ

_{x}*x*/2)

^{2}.

The mismatch between theory and simulation is not so surprising since the theory is for a weak gradient applied to a system that is essentially near equilibrium; this is not true in this setup. The only way to get this isothermal system out of equilibrium is via the boundaries, so the system is actually far from chemical equilibrium near the boundaries and then goes to chemical equilibrium in the middle of the domain, but in the middle the gradient disappears. A new more sophisticated theory is required that linearizes not around a constant state but rather around a spatially dependent state (this is automatically done in our code). Also, boundaries (i.e., confinement effects) may need to be included, especially for penetration depths comparable to the system size.

### C. Pattern formation

Since the seminal work of Turing,^{121} pattern formation in deterministic reaction-diffusion systems has been investigated extensively, mostly in theoretical studies but also by laboratory experiments.^{122,123} The study of stochastic systems is more recent and, as described in the Introduction, has primarily focused on models based on a RDME. Such a model was introduced by Mansour and Houard^{124} as a practical numerical scheme for the study of correlations in spatially distributed reactive systems. Subsequently, RDME-based models have been used to study the influence of fluctuations on pattern formation for a variety of reaction-diffusion systems.^{4–7} Recently, a SCLE was proposed^{46} and the RDME, SCLE, and deterministic equations were compared for the development of patterns in the Gray-Scott model; it was observed that the SCLE is qualitatively similar to the RDME for the majority of examined sets of parameters, but not always. All of the RDME-based models usually use simplified descriptions of diffusion, but recently it has been observed that accounting for cross-diffusion effects (which are included in complete generality in our formulation) may lead to qualitatively different behavior for Turing instabilities.^{125–128} Particle simulations including full hydrodynamics have also been performed using the DSMC method^{129} and molecular dynamics;^{8} these are, however, limited to small systems in (quasi) one dimension because of the high computational cost of particle simulations.

In our final example, we consider pattern formation in BPM model (27) for a dilute gas mixture with full hydrodynamic transport. The system is initialized in a uniform constant “reference” state in which the number densities of the different species are as specified in Table III. These number densities and the reaction rates are set so that the reference state is similar to that investigated in Ref. 70. Specifically, under the assumption that the number densities of the reservoir or “solvent” species *S*, *U _{f}*, and

*V*are fixed, the deterministic dynamics of the three reactive species

_{f}*U*,

*V*, and

*W*starts close to the single unstable fixed point; the stable attractor of the dynamics is a limit cycle around this unstable point. Of course, when the number densities of the solvent species are not fixed the dynamics is six-dimensional and much more complex. Since we consider a time-dependent non-equilibrium scenario, the chemical Langevin form of noise (24) was used for the stochastic chemistry.

Reaction . | $ R 1 $ . | $ R 2 $ . | $ R 3 $ . | $ R 4 $ . | $ R 5 $ . |
---|---|---|---|---|---|

k^{+} | 2.0 × 10^{−11} | 2.0 × 10^{−11} | 3.3333 × 10^{9} | 3.3333 × 10^{12} | 3.3333 × 10^{12} |

k^{−} | 2.0 × 10^{−17} | 2.0 × 10^{−12} | 3.3333 × 10^{−1} | 3.3333 × 10^{5} | 3.3333 × 10^{5} |

Reaction . | $ R 1 $ . | $ R 2 $ . | $ R 3 $ . | $ R 4 $ . | $ R 5 $ . |
---|---|---|---|---|---|

k^{+} | 2.0 × 10^{−11} | 2.0 × 10^{−11} | 3.3333 × 10^{9} | 3.3333 × 10^{12} | 3.3333 × 10^{12} |

k^{−} | 2.0 × 10^{−17} | 2.0 × 10^{−12} | 3.3333 × 10^{−1} | 3.3333 × 10^{5} | 3.3333 × 10^{5} |

Species . | U
. | V
. | W
. | S
. | U
. _{f} | V
. _{f} |
---|---|---|---|---|---|---|

Boundary n _{r} | 1.513 × 10^{19} | 4.38 × 10^{18} | 3.8 × 10^{17} | 5.0 × 10^{20} | 1.33 × 10^{20} | 5.0 × 10^{20} |

Species . | U
. | V
. | W
. | S
. | U
. _{f} | V
. _{f} |
---|---|---|---|---|---|---|

Boundary n _{r} | 1.513 × 10^{19} | 4.38 × 10^{18} | 3.8 × 10^{17} | 5.0 × 10^{20} | 1.33 × 10^{20} | 5.0 × 10^{20} |

In the standard BPM model, the solvent species (*S*, *U _{f}*, and

*V*) have fixed concentrations but in our hydrodynamic simulations they were fixed only at the boundaries. These three species are also made abundant to buffer them from having rapid variations in concentration (see Table III). In the limit of infinite concentrations of the solvent species (

_{f}*S*,

*U*, and

_{f}*V*), the dynamics approaches a reaction-diffusion model in which advection as well as momentum and heat transport become negligible. The reference (initial) values for mole fraction are used to prescribe Dirichlet boundary conditions for species on each side of the domain. This setup mimics open reservoirs in the form of permeable membranes;

_{f}^{43}note, however, that implementing boundaries that are also open for advective mass transport (i.e., inflow and outflow) is quite challenging

^{130}and not presently supported in our implementation. At the boundaries, the temperature is fixed at

*T*= 300 K and the fluid velocity is set to zero (no-slip) so species transport is primarily due to mass diffusion.

In our hydrodynamic simulations, the fluid is modeled as a hard sphere dilute gas so the transport coefficients depend on the masses and diameters of the particles of each species. The particles for all species in the BPM model have equal mass (*m* = 6.64 × 10^{−26} g) so as to ensure that the reactions conserve mass. For all species, the particles have only translational energy and no internal degrees of freedom (i.e., *z* = 0) so pressure and enthalpy are unaffected by reactions. In the BPM model, species *U* plays the role of the “inhibitor” while species *V* is the “activator.”^{123} Typically, pattern formation occurs when the inhibitor diffuses faster than the activator so we set the diameter of species *U* particles to be smaller than that of *V* particles, specifically *d*_{1} = 0.125 nm and *d*_{2} = 0.5 nm. Since we take $ k 1 + \u226b k 1 \u2212 $ (see Table III), the intermediary species, *W*, supports the activation of *V* and thus we set the diameters of *V* and *W* to be equal (this makes these two species hydrodynamically indistinguishable). The diameters of the other species (*S*, *U _{f}*, and

*V*) are small,

_{f}*d*

_{4}=

*d*

_{5}=

*d*

_{6}= 0.025 nm, so they diffuse rapidly from the boundaries and within the system. These specific values of the diameters are chosen such that the self-diffusion coefficient of

*S*(and

*U*,

_{f}*V*) is roughly an order of magnitude larger than the diffusion coefficient of

_{f}*U*, which is itself an order of magnitude larger than the diffusion coefficient of

*V*(and

*W*).

The system is simulated in a rectangular domain that is divided into 128 × 128 × 1 cells with Δ*x* = Δ*y* = 100 nm. The magnitude of the noise is varied by varying the domain thickness, which was either Δ*z* = 100 nm (low noise) or 10 nm (high noise). The reference value for the number of molecules per cell for the species of interest, *U*, is *O*(10^{4}) in the former case and *O*(10^{3}) in the latter. The total number of solvent molecules per cell is *O*(10^{6}) for weak noise and *O*(10^{5}) for strong noise. The time step is Δ*t* = 10 ps, as determined from stability requirements for the explicit temporal integrator.

Figure 6 illustrates the pattern formation observed in the system for low noise (top row), high noise (middle row), and deterministic evolution (bottom row) started from a perturbed initial condition generated by the high noise simulation (see figure caption). The boundaries take some time to influence the center of the domain, so in the center the reservoir species are depleted and the system moves toward chemical equilibrium. However, the boundary continuously forces the system so eventually spotted patterns develop, starting near the boundary, eventually filling the system. The resulting patterns are qualitatively similar to the “*λ* pattern” observed by Pearson^{102} in the GS model. In simulations with only species diffusion (i.e., setting all other transport to zero) we find similar patterning, indicating that this system is well-approximated by reaction-diffusion due to very large solvent concentrations.

In Ref. 121 Turing writes, “Another implicit assumption concerns random disturbing influences. Strictly speaking one should consider such influences to be continuously at work. This would make the mathematical treatment considerably more difficult without substantially altering the conclusions.” However, we see from Fig. 6 that the evolution is qualitatively different when large spontaneous fluctuations are present (middle row), as compared to when they are absent (bottom row). Specifically, the speed at which the patterns develop and propagate is noticeably accelerated by the spontaneous fluctuations (top and middle row), though the patterns themselves are qualitatively unchanged in this particular case. Other studies using reaction-diffusion models and particle simulations have reached similar conclusions.^{4–7} In Ref. 131 the authors investigated the Gray-Scott model by RDME simulations and concluded, “Complex spatiotemporal patterns, including spiral waves, Turing patterns, self-replicating spots and others, which are not captured or correctly predicted by the deterministic reaction-diffusion equations, are induced by internal reaction fluctuations.”

## V. CONCLUSIONS AND FUTURE WORK

In this work, we have formulated a fluctuating hydrodynamics model for chemically reactive ideal gas mixtures, and developed a numerical algorithm to solve the resulting system of stochastic partial differential equations. In our Langevin formalism, the stochastic mass, momentum, and heat flux as well as the stochastic chemical production rate, are modeled using uncorrelated white noise processes, and the local number densities are real variables. This is in contrast to the more traditional CME description of reactions that accounts for every individual reaction as a small jump of the (potentially very large) integer number of reactant molecules. We formulated the thermodynamic driving force for chemical reactions in agreement with nonlinear nonequilibrium thermodynamics^{59} and considered two different Langevin formulations of the stochastic chemical production rate. The first formulation is based on the law of mass action cast in the GENERIC framework^{62} and leads to a noise covariance that is the logarithmic mean (LME) of the forward and reverse production rates. This formulation is fully consistent with *equilibrium* statistical mechanics, more specifically, the resulting dynamics is time reversible (i.e., satisfies detailed balance) with respect to the Einstein distribution for a closed system. The second formulation is based on the CLE,^{47,57,69} and while it is not consistent with equilibrium statistical mechanics this form has its own merits.

We compared the two formulations on two chemical reaction networks for a well-mixed system, for both a simple dimerization reaction and a more complex network exhibiting bistability. We confirmed that at thermodynamic equilibrium the LME is more appropriate than the CLE, however, this is reversed for systems away from equilibrium, when compared with the CME. Not unexpectedly, neither is found to be entirely appropriate for describing rare events or large deviations from equilibrium. These examples remind us that a stochastic differential equation, which is a diffusion process, cannot be a uniformly accurate approximation for the CME, which is a Poisson process; the large-deviation statistics for Poisson noise is different than that of Gaussian noise. To further complicate the picture, there are known examples in which the discrete (integer-valued) nature of molecular populations plays a key role, implying that descriptions using real-valued concentrations such as fluctuating hydrodynamics must fail. For example, in wave fronts of the Fisher, Kolmogorov, Petrovski, Piskunov (FKPP) type it has been shown that the discreteness of the ME induces a logarithmic correction to the wave speed,^{132} similar to that observed when introducing a small cutoff in the leading edge of a FKPP front.^{133}

Another alternative coarse-grained description of stochastic chemistry, which we did not consider in this work, is tau leaping.^{66,134,135} Tau leaping is usually seen as a numerical method for approximately solving the CME, and the CLE can be seen as an approximation to tau leaping in a specific (central) limit in which a Poisson and a Gaussian variable become indistinguishable;^{66} observe, however, that the two kinds of distributions always have different tails. A different and more interesting characterization of tau leaping is to see it instead as an *alternative* to the CLE that maintains the Poisson nature of the noise rather than replacing it with Gaussian noise. In the limit that the time step Δ*t* → 0, for Gaussian noise one gets the CLE, and for Poisson noise one gets, in principle, the CME. In this limit, using the SSA algorithm is an efficient method to solve the CME exactly, and tau leaping is only useful as a numerical scheme for large Δ*t*. As mentioned earlier, computational fluctuating hydrodynamics should only be considered useful (or even valid) when each update of the coarse-grained degrees of freedom involves an average over many molecular events, such as many molecular collisions for momentum and energy transport, or many reactive collisions for chemistry. In other words, to distinguish it from the CME and the associated SSA algorithm, for tau leaping one should choose the time step size in a way that ensures that many reactions occur in each reactive channel. In this sense, tau leaping can be seen as a coarse-graining in time of the CME jump process, and, when combined with spatial coarse-graining, has the potential to be a useful coarse-grained model that bypasses the need for Langevin models of chemistry in our numerical schemes for reactive fluctuating hydrodynamics. Additional studies are needed to access the accuracy of tau leaping in situations where Langevin descriptions do poorly and such investigations are in progress.

As is well-known, the failure of Langevin approximations to describe large deviations is in fact closely connected to the fact that traditional linear nonequilibrium thermodynamics fails to describe chemical reactions because the entropy production rate is generally a non-quadratic function of the thermodynamic driving force (affinity). The mesoscopic Kramers picture of chemical reactions, as developed for isomerization in Ref. 59, is an interesting approach, which, however, remains mostly of theoretical utility; numerical simulations of this model would need to handle additional dimensions, as well as very slow diffusion across the reaction barrier. Furthermore, it is not obvious how to extend this description to general multispecies fluids with complex reaction networks.

While the CME is a well-agreed upon and well-justified description of statistically homogeneous systems, the story is much less clear for systems with spatial inhomogeneity; in fact, the precise mathematical meaning of “well-mixed” and the range of validity of the RDME remains obscure.^{21–23} A law of large numbers has been rigorously proven for several particle models and takes the expected form of a deterministic reaction-diffusion partial differential equation.^{72,73} Regarding fluctuations, however, there are very few particle models for which central limit theorems^{136} or large deviations functionals^{137} are known explicitly. It has been demonstrated that inhomogeneity leads to qualitative changes in the nature of phase transitions in bistable systems.^{138} It is also known that fluctuations can effectively renormalize the macroscopic transport and lead to non-analytic corrections to the law of mass action.^{2,3} It remains to be seen whether *nonlinear* spatially extended fluctuating hydrodynamics models can correctly reproduce this effect compared to particle simulations.^{139}

Furthermore, most particle models used for reaction-diffusion problems assume independent random Brownian walkers that react when coming near each other, and thus completely neglect transport mechanisms such as advection, sound waves, and thermal conduction. Furthermore, the mechanism of diffusion used in these models implicitly neglects the long-ranged hydrodynamic correlations present among particles diffusing in a viscous solvent.^{52} Notable exceptions are variants of the DSMC, which use a dilute^{140,141} or dense^{142} gas kinetic theory description of momentum and energy transport fully consistent with fluctuating hydrodynamics.^{142} While chemical reactions are commonly included in DSMC schemes,^{143,144} further studies of fluctuations and their consistency with nonequilibrium thermodynamics are needed. While some investigations of spatially distributed reactive systems have been performed using DSMC,^{129} a careful comparison to coarse-grained mesoscopic descriptions such as fluctuating hydrodynamics is needed. To model realistic chemistry, modern DSMC codes use either the Total Collision Energy (TCE) model^{140} or the more recent Quantum-Kinetic (QK) model,^{145} both of which we plan to compare with our formulations of reactive fluctuating hydrodynamics.

In Sec. IV B we studied the coupling of velocity fluctuations and chemistry in a system kept in a non-equilibrium steady state via boundary conditions. We found that, in agreement with existing theoretical computations, the chemical reactions have a strong effect on the giant long-range correlated concentration fluctuations. In this work we focused on gas mixtures, for which the Schmidt number is of order unity. We found that reactions profoundly change the nature of the giant fluctuations, whose spectrum switches from the well-known ∼*k*^{−4} at large wavenumbers to a flat plateau (with a value controlled by the reaction rate) at small wavenumbers. This could be useful, for example, for experimentally measuring reaction rates. However, we found that the simple quasi-periodic near-equilibrium theory we constructed was not in quantitative agreement with the simulations, indicating that a more precise theory is needed.

Reactive *liquid* mixtures are common in practice and exhibit interesting coupling of hydrodynamics with transport that has been studied both theoretically^{146} and experimentally.^{9} The Schmidt number in liquids is, however, very large, and the compressibility is very small (i.e., the speed of sound is very large), making the method presented here computationally infeasible. In the future, we will consider extending low Mach number (quasi-incompressible) methods that treat momentum diffusion implicitly to include stochastic chemistry, by combining Langevin or tau-leaping based descriptions of chemistry with the formulation and numerical methods developed in Refs. 44 and 147 for general non-ideal liquid mixtures.

Finally, the influence of hydrodynamic fluctuations on reactions will likely be very important for surface chemistry.^{148,149} An important application is heterogeneous catalysis in which a highly reactive catalytic surface facilitates bond breaking and bond rearrangement of adsorbed molecules. In this context, mesoscale simulations are particularly useful for the study of nano-catalytic systems.^{150} Microscopic catalytic particles have many advantages, such as higher activity, increased selectivity, and longer lifetime. However, to operate effectively these particles must have electrical contact with a substrate (typically a flat surface) and they must not be so small as to not have enough electrons to catalyze a reaction. For example, nanowire catalysts are typically 10-100 nm in size so the hydrodynamic environment in which the chemistry and transport occur is of mesoscopic scale (a few micrometers).

As in the case with chemistry in bulk flow, particle-based simulations will be useful benchmarks for comparison with reactive fluctuating hydrodynamics modeling surface chemistry. Furthermore, molecular simulations can be embedded within a fluctuating hydrodynamic code to create an Algorithm Refinement (AR) hybrid.^{39,151} The idea is to use a particle-based simulation for the domain near the surface to capture the physics at the molecular scale, such as adsorption of reactants onto the surface, surface diffusion, surface reactions, and desorption of products. This particle-based simulation would then be coupled to a fluctuating hydrodynamic simulation that treats the bulk fluid. Previous work for non-reactive fluids has already proven the utility of AR hybrids for simple interfaces^{53} and this numerical framework should prove useful in the study of surface reactions and active membrane transport.

## Acknowledgments

We would like to thank M. Malek-Mansour, Jonathan Goodman, Eric Vanden-Eijnden, Samuel Isaacson, Hans Christian Öttinger, Dick Bedeaux, Annie Lemarchand, Florence Baras, John Pearson, Sorin Tanase Nicola, and Signe Kjelstrup for informative discussions. This material is based upon work supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award No. DE-SC0008271 and under Contract No. DE-AC02-05CH11231.

### APPENDIX A: DETERMINISTIC AND STOCHASTIC TRANSPORT IN IDEAL GAS MIXTURES

This appendix summarizes fluctuating hydrodynamics for ideal gas mixtures; for a more general and detailed exposition see Ref. 42. Each of the hydrodynamic transport terms in (1)-(3) contains a deterministic term, denoted with an overbar and a stochastic term denoted by a tilde (e.g., $F= F \xaf + F \u02dc $). The deterministic viscous tensor is

where *η* and *κ* are the shear and bulk viscosity, respectively. We neglect any possible effect of the chemical reactions on the transport coefficients of the mixture. For example, we neglect any coupling between bulk viscosity and chemical reactions, which, in principle is allowed by the Curie principle since both are scalar processes.^{74} The corresponding stochastic viscous flux tensor can be written as^{36,152}

The symmetric Gaussian random tensor field, $ Z \u0302 \Pi $, is formulated as $ Z \u0302 \Pi = Z \Pi + ( Z \Pi ) T / 2 $, where $ Z \Pi $ is a white-noise random Gaussian tensor field, that is, $\u3008 Z \Pi \u3009=0$, and

with *α*, *β*, *γ*, *δ* = {*x*, *y*, *z*} being spatial components.

The deterministic mass flux and heat flux depend on the gradients of concentration, pressure, and temperature. For an ideal gas, we can write the deterministic fluxes as

and

where *h* is a vector of specific enthalpies. Here, the diffusion driving force^{86} is

and $X$, $Y$, and $M$ are diagonal matrices of mole fractions *X*, mass fractions *Y*, and molecular masses *M*. The matrix of multicomponent flux diffusion coefficients, **D**, the vector of rescaled thermal diffusion ratios, $ \chi \u0303 $, and the thermal conductivity, *λ*, can be obtained from standard software libraries, such as EGLIB,^{153} or from standard references, such as Hirschfelder *et al.*^{154}

For the ideal gas transport coefficients described above, we define

where $ m \xaf = ( \u2211 s Y s / m s ) \u2212 1 $ is the mixture-averaged molecular mass. The stochastic terms for the combined species equations and energy equation are determined from the phenomenological equations of nonequilibrium thermodynamics that relate fluxes to thermodynamic driving forces through the Onsager matrix, and the fluctuation-dissipation balance principle.^{42} Specifically, the stochastic mass flux is

where **BB**^{T} = 2*k _{b}*

**L**, and

**𝒵**

^{ℱ}is a white-noise random Gaussian vector field with uncorrelated components, that is, 〈

**𝒵**

^{ℱ}〉 = 0 and

with *α*, *β* = {*x*, *y*, *z*} being spatial components. Note that the matrix **B** can be obtained using Cholesky decomposition; the constraint $ \u2211 s F \u02dc s =0$ is ensured by construction.^{42} The stochastic heat flux is

where $\u3008 Z \alpha Q ( r , t ) Z \beta Q ( r \u2032 , t \u2032 ) \u3009= \delta \alpha \beta \u2009\delta ( r \u2212 r \u2032 ) \u2009\delta ( t \u2212 t \u2032 ) $.

### APPENDIX B: EQUILIBRIUM DISTRIBUTION FOR A DIMERIZATION REACTION

The entropy of mixing for a system undergoing dimerization is *S* = *k _{B}*ln

*N*, where

_{p}*N*is the number of distinct ways of forming

_{p}*M*dimers out of a total of

*N*monomers. This is straightforward to compute. The number of ways to choose 2

*M*out of the

*N*atoms to be in pairs is

Next, we need to group the 2*M* atoms into pairs; the number of distinct ways of pairing 2*M* objects is (2*M*)!/2^{M}*M*! This gives the entropy

where we have included an additional reference chemical potential *μ* to set the equilibrium concentration.

By expanding the logarithm of the right-hand side of (B1) using Stirling’s leading-order approximation, we obtain the thermodynamic limit of the entropy as a function of the monomer mass fraction $ Y 1 = N \u2212 2 M /N$. After fixing the chemical potential from the requirement that the most probable mass fraction is (*Y*_{1})_{eq} = 1/2, we get

This gives an Einstein distribution *P* ∼ *e*^{S/kB} exactly matching (28).

It is useful to compare this equilibrium distribution of the LME with that of the CME for a well-mixed system of volume Δ*V*. It is not hard to show that the CME for a dimerization reaction is in detailed balance with respect to the Einstein distribution with entropy (B1), with the reference chemical potential set to

We can use this to construct a rather accurate continuum approximation to this exact microscopic result by including the next order term in the Stirling formula,

### APPENDIX C: GIANT FLUCTUATIONS IN THE PRESENCE OF REACTIONS

This appendix outlines the fluctuating hydrodynamics theory for the long-range correlations of concentration fluctuations in a binary mixture undergoing a dimerization reaction, as studied numerically in Sec. IV B. We neglect the Dufour effect and assume the system to be isothermal, taking contributions from temperature fluctuations to be of higher order. Furthermore, we neglect gravity, assume the system is incompressible, and take the density and transport coefficients to be constant. We consider a “bulk” system,^{35} i.e., we neglect the influence of the boundaries. This gives an accurate approximation for wavenumbers that are large compared to the inverse height of the domain; for smaller wavenumbers the boundaries are expected to suppress the giant fluctuations.^{35,155} We also neglect the stochastic mass flux in the concentration equation since we are concerned with the nonequilibrium contribution due to the forcing by the velocity fluctuations.

We assume that all of the concentration gradients are in the same direction (say, the *y* axis). The incompressibility constraint is most easily handled by applying a **∇** × **∇** × operator to the momentum equation^{35} to obtain a system involving only the component of the velocity parallel to the gradient, *v*_{∥} ≡ *v _{y}*. The same calculation can easily be generalized to a multispecies mixture as well, see Appendix B in Ref. 42. This system of equations can be most easily solved in the Fourier domain, by considering wavevectors

**in the plane perpendicular to the gradient,**

*k***=**

*k*

*k*_{⊥}.

It is very straightforward to derive (31) in the limit of large Schmidt number by considering the fluctuating concentration equation forced by an *overdamped* (steady Stokes) fluctuating velocity. The steady Stokes equation in Fourier space has the form

where *T*_{0} is the constant temperature. Here $ W \u0302 ( t ) $ is a white-noise process (one per wavenumber), giving the white-in-time velocity $ v \u02c6 \u2225 ( t ) = k \u2212 1 2 k B T 0 / \eta i W \u0302 ( t ) $. A general form of the linearized equation for the concentration fluctuations, $c=\delta Y 1 = Y 1 \u2212 Y 1 $, in Fourier space will have the form

where *ψ* is the reaction rate at the equilibrium point (around which we are linearizing), and $f=d Y 1 /dy$ is the applied concentration gradient. For our specific reaction and the choice of equilibrium point (*Y*_{1})_{eq} = 1/2, we have that *ψ* = 3*k*^{−} (to see this simply linearize (33) around (*Y*_{1})_{eq} = 1/2). The resulting nonequilibrium static structure factor (recall that here we neglect the contribution due to the stochastic mass flux) is straightforward to calculate (see, for example, Appendix B in Ref. 42)

which is exactly (31) with the penetration depth

Since the Schmidt number is not very large for gases, we should improve the theory by not taking the overdamped limit but rather adding a velocity equation and considering the linearized *inertial* equations,

where *ρ*_{0} is the equilibrium density and *ν* = *η*/*ρ*_{0} is the kinematic viscosity. Solving this system of linear SODEs gives the improved structure factor

which shows that the finite value of the Schmidt number *S _{c}* =

*ν*/

*D*has an important effect on the giant fluctuations. In particular, in the more complete theory (C2),

is finite and not infinite as in (C1), which assumes infinite Schmidt number.

In the more complete theory, we see the following three regimes.

- For large wavenumbers we haveas if there were no reaction.$ S neq k \u226b D \psi \u2248 k B T 0 \eta D k 4 f 2 $
- If the Schmidt number is small,
*D*∼*ν*, as for gases, then for small wavenumbers we instantly switch to a flat spectrum$ S neq k \u226a D \psi \u2248 k B T 0 \rho 0 \psi 2 f 2 . $ - If the Schmidt number
*S*is large,_{c}*D*≪*ν*, as for liquids, then for intermediate wavenumberswe observe a change in the power law to$ S neq D \psi \u226a k \u226a S c 1 2 D \psi \u2248 k B T 0 \eta \psi k 2 f 2 $*S*_{neq}(*k*) ∼*k*^{−2}. Note, however, that this range spans only $ S c $ orders of magnitude in*k*, so even for*S*∼ 10_{c}^{4}(which applies to macromolecular solutions), the*k*^{−2}power-law only extends over at most two decades. Therefore, even for liquids with large Schmidt numbers the inertial equations should be used to model giant fluctuations in reactive mixtures.

## REFERENCES

Note that the two terms in (15) can be combined together to give a single white-noise process with amplitude $DrCL=Dr,+CL+Dr,\u2212CL$; this leaves the covariance of the stochastic forcing unchanged.

This is in the classical regime, that is, when the molecules’ mean spacing is much larger than their de Broglie wavelength.

To implement the BPM model in a molecular simulation with binary collisions, Baras *et al.* modify the model by making all the reactions bimolecular and by introducing an auxiliary species, *A*.

The precise mathematical definition of a “well-mixed” system is delicate and will not be discussed here at length. Roughly speaking, it means that diffusion is sufficiently fast compared to reactions; see Ref. 73 for a theorem on the limit of infinite diffusion.

In the literature, this equation is called the CFPE (chemical FPE), see (4) in Ref. 65.

The linearized CLE is known as the second-order van Kampen expansion^{156} in the physics literature and has the same formal order of accuracy as the CLE^{69} but is much simpler to solve analytically due to its linearity. Analysis in Ref. 65 suggests that the nonlinear CLE may be more accurate than the linearized CLE for large noise but we are not aware of rigorous mathematical estimates.

While Langevin approximations cannot approximate atypical statistics of the CME, the CLE is likely appropriate for describing the typical behavior of the CME.^{65}

The nonlinearity of the deterministic (macroscopic) equations is fully taken into account.

Note that in this case the mismatch between the theory and simulations at very small *k* is due to the effect of boundaries.^{155}

For the more resolved runs Δ*x* = Δ*y* = 5 × 10^{−4}, Δ*t* = 6.25 × 10^{−8} for the problem without hydrodynamics and Δ*t* = 6.25 × 10^{−9} with hydrodynamics.

Fluctuating hydrodynamics simulations have been shown to correctly describe the renormalization of diffusion coefficients in binary fluid mixtures,^{114} but we are not aware of any work in this direction for reactive systems.