Fluctuating hydrodynamics (FHD) provides a framework for modeling microscopic fluctuations in a manner consistent with statistical mechanics and nonequilibrium thermodynamics. This paper presents an FHD formulation for isothermal reactive incompressible liquid mixtures with stochastic chemistry. Fluctuating multispecies mass diffusion is formulated using a Maxwell–Stefan description without assuming a dilute solution, and momentum dynamics is described by a stochastic Navier–Stokes equation for the fluid velocity. We consider a thermodynamically consistent generalization for the law of mass action for non-dilute mixtures and use it in the chemical master equation (CME) to model reactions as a Poisson process. The FHD approach provides remarkable computational efficiency over traditional reaction-diffusion master equation methods when the number of reactive molecules is large, while also retaining accuracy even when there are as few as ten reactive molecules per hydrodynamic cell. We present a numerical algorithm to solve the coupled FHD and CME equations and validate it on both equilibrium and nonequilibrium problems. We simulate a diffusively driven gravitational instability in the presence of an acid-base neutralization reaction, starting from a perfectly flat interface. We demonstrate that the coupling between velocity and concentration fluctuations dominates the initial growth of the instability.

## I. INTRODUCTION

Thermal fluctuations in fluids arise from random molecular motions, driving both microscopic behavior and macroscopic behavior that deterministic models fail to predict. In diffusive mixing experiments, velocity fluctuations lead to giant fluctuations in concentration in the presence of concentration gradients.^{1} Buoyancy-driven instabilities can be triggered or affected by thermal fluctuations.^{2,3} In reaction-diffusion systems, thermal fluctuations can accelerate the formation of Turing patterns on a macroscopic time scale^{4} and induce long-time memory in the chemical kinetics of a diffusion-limited system.^{5}

In this paper, we develop a formulation and numerical methodology for the stochastic simulation of reactive microfluids. Here we incorporate a stochastic description of chemical reactions based on the chemical master equation (CME)^{6} into an isothermal fluctuating hydrodynamics (FHD)^{7,8} description of diffusive and advective mass transport. Hence, our proposed algorithm combines discrete processes (CME for reactive processes) and continuous processes (FHD for transport processes); a similar idea has been used for the simulation of the Boltzmann equation.^{9–11} The use of the CME enables us to correctly capture large fluctuations of composition, going beyond the Gaussian approximation inherent in the chemical Langevin equation (CLE) used in our prior work.^{12} While our previous work on reaction-diffusion systems^{4} also employed the CME, it was restricted to dilute solutions. Here we generalize the CME to non-dilute ideal mixtures with a complete Maxwell–Stefan formulation of diffusive transport in multispecies mixtures. This includes cross-diffusion coupling among distinct species and can account for deviations from ideality, unlike the standard reaction-diffusion master equation (RDME) approach.^{13–15} Finally, by including the fluctuating Navier–Stokes equations in the model, we account for advection by thermal velocity fluctuations, which is necessary to capture giant nonequilibrium composition fluctuations.^{10,11,16}

Our approach is related to, but also distinct from, prior work on fluctuating hydrodynamics for reactive liquid mixtures. An alternative Langevin-based approach proposed in Ref. 17, and extended to full hydrodynamics in Ref. 18, represents reactions as a diffusion process along an internal reaction coordinate, driven by the Gaussian noise. This description is fully consistent with nonequilibrium thermodynamics and fluctuating hydrodynamics, but is not easily extensible to multispecies mixtures, and, importantly, is expensive to use in numerical simulations because it requires introducing an additional reaction coordinate, thus effectively increasing the dimensionality of the problem. Instead, in our approach, we only consider the reactant and product states and consider reactions as a jump process between these two states, driven by the Poisson noise. The deterministic (macroscopic) as well as a linearized version of the FHD equations we consider here is the same as that obtained from a quasi-stationary approximation of the model developed in Ref. 18 [see Eq. (26) in Ref. 18 and the book by Keizer^{19}), as we have discussed in more detail in prior work.^{12} The key difference is that here we describe chemical fluctuations using a *nonlinear* FHD description based on a master equation, rather than a linearized Langevin description. This is common in stochastic reaction-diffusion models used in biochemical modeling,^{20,21} as we have discussed in more detail in prior work.^{4} However, traditional RDME descriptions have been restricted to dilute solutions and do not account for velocity (momentum) fluctuations. More broadly, biochemical reaction-diffusion models have largely been developed without input from the field of (non)equilibrium thermodynamics, and especially fluctuating hydrodynamics. Here we bridge this gap by combining features of the RDME with FHD, thus delivering on the promise made in Ref. 4 to “explore combining Langevin and CME approaches together, thus further bridging the apparent gap between the two.” Giant nonequilibrium fluctuations, which arise due to the coupling with velocity fluctuations, have been studied theoretically using linearized FHD for a dimerization reaction in Refs. 22 and 23. Here we study giant fluctuations in a liquid mixture undergoing a dimerization reaction numerically and show that a quantitatively accurate theoretical description is difficult due to the nonlinearity of the macroscopic steady state.

In this work, we simplify our previous variable-density low Mach FHD formulation by restricting it to miscible liquid mixtures,^{3} in which the density is essentially independent of composition at a fixed pressure and temperature. The resulting Boussinesq (incompressible) approximation of the momentum equation enables us to construct an efficient numerical method that accounts for inertial effects important in buoyancy-driven fluid flows yet remains robust for small Reynolds numbers and large Schmidt numbers. The spatio-temporal discretization of the FHD equations is based on our previous work^{3} but with some important improvements necessary for simulating complex reactive mixtures at small length scales. Notably, we extend our previous work on reaction-diffusion systems^{4} to general multispecies mixtures so that large deviations of composition are handled accurately and robustly, and negative densities are avoided.

We follow a general framework for the systematic construction of FHD numerical methods based on the stochastic version of the method of lines approach.^{24} Using this framework, we have previously developed stochastic simulation methods for gas mixtures^{25} and quasi-incompressible miscible liquid mixtures.^{3,26,27} For liquid mixtures, we have developed a computationally efficient low Mach number model that eliminates fast pressure waves while preserving the spatio-temporal spectrum of the slower diffusive fluctuations.^{27} To avoid severe restriction on the time step size when the Schmidt number is large, we have developed an implicit temporal discretization of viscous dissipation^{26} that relies on a variable-coefficient multigrid precondition to solve the coupled velocity-pressure Stokes system.^{28}

In this paper, we make three novel contributions to the numerical methodology developed in our prior work. First, by incorporating a second-order midpoint tau-leaping scheme^{29} into our prior algorithms for multispecies miscible liquid mixtures,^{3} we construct a numerical method that efficiently samples reactions at a cost no larger than that of integrating the chemical Langevin equation. Because our novel midpoint temporal integrator solves the CME by using tau leaping,^{30} it is robust for large composition fluctuations, while also being efficient for weak fluctuations. Second, the midpoint scheme is constructed to be robust for large Schmidt numbers, i.e., much faster momentum diffusion compared to mass diffusion, as is typical in liquid systems. In particular, the numerical method reproduces the correct spectrum of giant nonequilibrium fluctuations even for time step sizes much larger than the stability limit dictated by fast momentum diffusion, while also preserving the slow inertial momentum dynamics at large scales. Third, we take careful attention to handling vanishing species robustly both in the formulation of the multispecies diffusion model and in the numerical algorithm.

The rest of the paper is organized as follows. In Sec. II, we present the formulation of the FHD equations coupled with the CME formulation of reactions. In Sec. III, we present a numerical scheme that can solve these equations accurately and robustly even in the presence of large composition fluctuations and vanishing species. In Sec. IV, we present numerical results for four examples and discuss various aspects of our numerical method and the effects of thermal fluctuations. First, we verify that for dilute solutions our algorithm preserves the robustness and accuracy properties of our previous method for reaction-diffusion systems^{4} by modeling the hydrolysis of sucrose at micrometer scales. Second, to assess the fidelity of our approach in a non-dilute setting, we consider a binary mixture undergoing a dimerization reaction 2A *⇌* A_{2} at thermodynamic equilibrium with a small number of molecules per cell. Third, we also study such a mixture out of equilibrium in the presence of giant nonequilibrium fluctuations with a large number of molecules per cell. Fourth, we use our numerical algorithm to simulate a diffusively driven gravitational instability in the presence of an acid-base neutralization reaction recently studied experimentally^{2} and show that the coupling between velocity and concentration fluctuations triggers and drives the instability at early times. In Sec. V, we conclude the paper with a brief summary and a discussion of future directions.

## II. REACTIVE FLUCTUATING HYDRODYNAMICS

Our formulation relies on several approximations appropriate for many isothermal miscible liquid mixtures. First, we neglect the effects of thermodiffusion and barodiffusion on mass transport and assume constant temperature *T* and thermodynamic pressure *P*. Second, we assume that density variations due to composition are small enough that they have no effect on the flow field except through a buoyancy force. Hence, we formulate our FHD system as an isothermal Boussinesq simplification of the low Mach number multispecies model used in Ref. 3. While a numerical method can be potentially constructed without these approximations, the Boussinesq formulation greatly reduces the complexity of the numerical scheme without losing essential physics.

Given these approximations, we recast the continuity equation for mass density as a divergence-free constraint on velocity and assume a constant density *ρ*_{0},

Here, $\upsilon $ is the fluid velocity, *π* is the mechanical pressure (a Lagrange multiplier that ensures that the velocity remains divergence free^{31}), *η*($w$) is the viscosity, $\u2207\xaf=\u2207+\u2207T$ is a symmetric gradient, and **Σ** is the stochastic momentum flux. By denoting the number of species with *N*_{spec}, the vector of mass fractions (concentrations) is given by $w=(w1,\u2026,wNspec)$, where $ws$ is the mass fraction of species *s* and *∑*_{s}$ws$ = 1. We compute the mass density of each species using *ρ*_{s} = *ρ*_{0}$ws$, and thus the total mass density *∑*_{s}*ρ*_{s} = *ρ*_{0} is strictly constant. The buoyancy force ** f**($w$) is a problem-specific function of $w$. The total diffusive mass flux

*F*_{s}of species

*s*is decomposed into a dissipative flux $F\xafs$ and fluctuating flux $F\u0303s$,

and *m*_{s}Ω_{s} represents a source term representing stochastic chemistry, where *m*_{s} is the molecular mass and Ω_{s} is the number density production rate for species *s*. Note that by summing up (3) over all species, we recover (2) since *∑*_{s}*F*_{s} = **0** and *∑*_{s}*m*_{s}Ω_{s} = 0. Based on the fluctuation-dissipation relation, the stochastic momentum flux **Σ** is modeled as

where *k*_{B} is Boltzmann’s constant, and $\mathcal{Z}mom(r,t)$ is a standard Gaussian white noise (GWN) tensor field with uncorrelated components having *δ*-function correlations in space and time.

We formulate multispecies diffusion in Sec. II A and chemistry in Sec. II B. It is important to note that both the diffusion and chemistry formulations are obtained from a general form of the specific chemical potential for each species,

where $\mu s0(T,P)$ is a reference chemical potential and *γ*_{s}(** x**,

*T*,

*P*) is the activity coefficient (for an ideal mixture,

*γ*

_{s}= 1). Here

**denotes mole fractions, which can be expressed in terms of $w$ as**

*x*where $m\xaf$ is the mixture-averaged molecular mass,

In Sec. II C, we confirm the thermodynamic consistency of our formulation by showing that thermodynamic equilibrium is determined by the chemical potentials and that transport processes and reactions do *not* change equilibrium statistics. In Sec. II D, we discuss the simplification of our model for dilute solutions.

### A. Multispecies diffusion

Here we summarize the FHD description of multispecies diffusion formulated in Ref. 3. Neglecting thermodiffusion and barodiffusion, the Maxwell–Stefan formulation of the diffusion driving force gives

where **Γ** is the matrix of thermodynamic factors that becomes the identity matrix for ideal mixtures and ** W** is a diagonal matrix with entries $w$. The symmetric matrix

**Λ**is defined via

where $\u2212Dss\u2032$ is the Maxwell–Stefan binary diffusion coefficient between species *s* and *s*′. Denoting a pseudo-inverse of **Λ** with ** χ**, we can rewrite (9) as

The stochastic mass fluxes $F\u0303$ are given by the fluctuation-dissipation relation

where $\chi 12$ is a “square root” of ** χ** satisfying $\chi 12(\chi 12)T=\chi $, and $\mathcal{Z}mass(r,t)$ is a standard GWN field with uncorrelated components. Modifications of this formulation in the presence of trace or vanishing species are discussed in Sec. III A.

### B. Chemical reactions

We consider a liquid mixture undergoing *N*_{react} elementary reversible reactions of the form

where $\nu sr\xb1$ are the molecule numbers and $Ms$ are chemical symbols. We define the stoichiometric coefficient of species *s* in the forward reaction *r* as $\Delta \nu sr+=\nu sr\u2212\u2212\nu sr+$ and the coefficient in the reverse reaction as $\Delta \nu sr\u2212=\nu sr+\u2212\nu sr\u2212$. We assume that mass conservation holds in each reaction *r*; i.e., $\u2211s\Delta \nu sr\xb1ms=0$ for all *r*. It is important to note that all reactions must be reversible for thermodynamic consistency.

To sample Ω_{s}, we need propensity density functions $ar\xb1$ for the forward/reverse (+/-) rates of reaction *r*. Specifically, the mean number of reaction occurrences in a locally well-mixed reactive cell of volume Δ*V* during an infinitesimal time interval *dt* is given as $ar\xb1\Delta Vdt$. Accordingly, the mean number density production rate of species *s* is given as

#### 1. Generalized law of mass action

Here we adopt the canonical form for the rate of chemical reactions.^{19,32} Propensity density functions are expressed as^{12}

where *λ*_{r}(*T*, *P*) ≥ 0 is a reaction rate parameter assumed to be independent of the composition and $\mu ^s=ms\mu s/kBT$ is the dimensionless chemical potential per particle. For the general form of chemical potential (6), we have

where $\kappa r\xb1(T,P)=\lambda r\u220fs\u2061exp(\nu sr\xb1\mu ^s0)$ denotes the forward/reverse reaction rate constant. From the condition $ar+=ar\u2212$ at chemical equilibrium, we can express the equilibrium constant as a purely thermodynamic quantity,

as required by statistical mechanics.

It is important to note that propensity density functions and equilibrium constants are expressed in terms of *mole fractions x*_{s} (for ideal mixtures) or activities *x*_{s}*γ*_{s}. This generalized LMA has a different form compared to the *number density* based LMA used for ideal gas mixtures in our prior work.^{12} However, this does not imply any incompatibility between the two forms of LMA. For isothermal gas mixtures, pressure changes significantly upon reaction due to changes in mole numbers and thus $\kappa r\xb1(T,P)$ cannot be assumed to be constant. On the other hand, in liquid mixtures, where pressure changes are not significant, $\kappa r\xb1(T,P)$ can be assumed to be constant.

#### 2. CME-based stochastic chemistry

We believe that an accurate mesoscopic chemistry description should be based on a master equation approach, which leads to the CME^{20} for well-mixed^{33} systems. As will be demonstrated in Sec. III A, both the CME description and the generalized LMA are crucial for achieving thermodynamic consistency. Note, however, that our CME-based description itself does not require reversible reactions. For modeling purposes, one can exclude some forward or reverse reactions by assuming that they have zero rates. However, we remind the reader that this is inconsistent with equilibrium thermodynamics.

For reactions in a closed well-mixed cell of volume Δ*V*, the CME describes the time evolution of the system in terms of the temporal change in the *probability* of the system to occupy each state (specified by the number of molecules of each species). We use an equivalent, but more direct, *trajectory-wise* representation,^{4} which is related to the computationally efficient tau leaping method.^{30} The change in the number of molecules *N*_{s} of species *s* in a given cell during an infinitesimal time interval *dt* is expressed in terms of the number of occurrences $P(ar\xb1\Delta Vdt)$ of each reaction *r*,

where $P(m)$ denotes the Poisson random variables with mean *m*. Note that the instantaneous rate of change is written as an Ito stochastic term. The tau leaping method discretizes (18) with a finite time step size Δ*t*. To faithfully model the discrete nature of reactions, we sample *integer-valued* reaction counts using Poisson random numbers as in the traditional tau leaping algorithm. However, it is important to note that we use *continuous-ranged* number densities for advection-diffusion, and therefore cells are not guaranteed to have an integer number of molecules.

We note that a Gaussian approximation of the Poisson random number $P(ar\xb1\Delta Vdt)$ in (18) leads to the chemical Langevin equation (CLE). In this Langevin (Gaussian noise) approximation,^{4,12}

where $Zrreact(r,t)$ denotes a standard GWN field. The Langevin description is justified in the limit of small Gaussian fluctuations with respect to average concentrations.^{20} However, the Langevin description predicts an unphysical equilibrium state with negative densities and does not correctly model large deviations of chemical fluctuations.^{12} By contrast, the tau leaping method correctly reproduces the large deviation functional of the CME, while still remaining as efficient as the CLE [see the discussion around Eq. (8.1) in Ref. 34].

### C. Thermodynamic consistency

We now demonstrate the thermodynamic consistency of our formulation for ideal mixtures at thermodynamic equilibrium. For the simplicity of exposition, we consider a binary liquid mixture of A atoms and A_{2} molecules undergoing a dimerization reaction

noting that this analysis also applies to multispecies ideal mixtures. In Sec. III A, we consider the single-cell (homogeneous) case. We obtain the thermodynamic equilibrium distribution of monomers and dimers to show that our chemistry model satisfies detailed balance with respect to the correct Einstein equilibrium distribution. In Sec. III B, we consider the spatially extended case. We show that the governing Boussinesq equations give flat structure factors at thermodynamic equilibrium in the Gaussian approximation, in agreement with statistical mechanics.

#### 1. Single-cell system

We denote the number of monomers and dimers as *N*_{1} and *N*_{2}, respectively. By the constant density approximation, *N*_{1} and *N*_{2} satisfy *N*_{1} + 2*N*_{2} = *ρ*_{0}Δ*V*/*m* ≡ *N*_{0}, where *m* is the mass of a monomer and *N*_{0} is the total number of A atoms in a cell of volume Δ*V*. Hence, we denote the equilibrium distribution of the composition with *P*(*N*_{2}).

Statistical mechanics predicts that the equilibrium distribution is given by the Einstein distribution, $P\u223ceS/kB$, where *S* denotes the entropy of the system at a given state. Note that even though we consider isothermal systems, we can still use the Einstein distribution since the only contribution to the free energy that depends on composition is the entropy of mixing. For a binary ideal mixture, the entropy of mixing is given as

and the entropy of the system is

Hence, we obtain the equilibrium distribution

with $\u2211N2=0N0/2P(N2)=1$. Note that it is straightforward to obtain the ratio of occupation probabilities of adjacent states,

We now analyze when detailed balance is achieved for the dimerization reaction (20) with respect to the equilibrium distribution *P*(*N*_{2}). The detailed balance condition is given as

where *a*^{±}(*N*_{2}) denote the forward/reverse rates at the state with *N*_{2} dimers, which are to be determined. By using (17) and (24), one can show that the detailed balance condition (25) exactly holds for

with *N*_{1} = *N*_{0} − 2*N*_{2}. It is important to note that (26) reduces to $a+=\kappa +x12$ and *a*^{−} = *κ*^{−}*x*_{2} in the thermodynamic limit. Hence, (26a) can be considered as an integer-based correction to the generalized LMA (16); this correction makes sense because the probability of choosing a second monomer is (*N*_{1} − 1)/(*N*_{1} + *N*_{2} − 1). Such integer corrections are well known for low density solutions and used in most RDME models of reaction-diffusion systems, but to our knowledge they have not previously been formulated for non-dilute ideal mixtures.

In the thermodynamic limit, we can apply Stirling’s approximation to (21) and express chemical potentials in (22) in terms of equilibrium mole fractions $xseq$, to give

where *S*^{eq} denotes the entropy at *x*^{eq} and *N* = *∑*_{s}*N*_{s}. We can further approximate *S*_{Stirling} up to second order in $\delta xs=xs\u2212xseq$ (*s* = 1, …, *N*_{spec} − 1), to get a Gaussian approximation to the Einstein distribution. This Gaussian approximation is described by linearized FHD, and we study it in more detail, including spatial dependence, next.

#### 2. Spatially extended system

We can extend the dimerization results obtained for the single-cell case to the spatially extended case. For an ideal mixture, the total entropy of the system is additive over the individual cells,

where $N2(i)$ denotes the number of dimers in cell *i*. Therefore, the Einstein distribution for the spatially extended system is the product distribution

This means that the number of dimers in each cell is independent of those in the other cells and has the same distribution as the single-cell case.

We note that our FHD model of multispecies diffusion is constructed so that it reproduces the correct Einstein distribution under Stirling’s approximation; i.e., our model is consistent with (27) and (28). Hence, the combined chemistry and FHD model is expected to give the correct equilibrium distribution, as long as there are sufficiently many molecules of all species in each cell to justify the continuous approximation. In Sec. IV B, we numerically confirm that our method gives an accurate approximation to (29) even when there are significant fluctuations of composition, with as few as *N*_{2} ∼ 10 dimers per cell.

At the level of a Gaussian approximation, we can investigate the system analytically using the linearized FHD equations. We denote the mass fractions of monomers and dimers as $w$ and 1 − $w$, respectively. We assume that $w$ fluctuates around $w\xaf$. At equilibrium, our FHD equations (1)–(3) are linearized for $\upsilon $ = *δ*$\upsilon $ and $w=w\xaf+\delta w$ as follows:

where *ν* = *η*/*ρ*_{0}, $D=\u2212D12$, and $\mu w$ is the second order derivative of Gibbs free energy with respect to concentration $w$, given as $\mu w=kBT/[mw\xaf(1\u2212w\xaf2)]$ for an ideal mixture. The linearized reaction term is denoted by $\rho 0\u22121m\Omega 1lin$.

We denote the equilibrium structure factors (spectra) by $S\upsilon ,\upsilon eq(k)=\u27e8\delta \upsilon ^\delta \upsilon ^*\u27e9$ and $Sw,weq(k)=\u27e8\delta \u0175\delta \u0175*\u27e9$, where the hat denotes a Fourier transform and the asterisk denotes a conjugate transpose. Noting that the concentration equation is uncoupled from the momentum equation, these structure factors can be obtained separately. For the non-reactive case, they are independent of *k*,^{27}

It is easy to show that the spatial correlations of the composition fluctuations $Sw,weq$ are fully consistent with the Gaussian approximation of (29).

For the reactive case, (30c) contains the stochastic chemistry term

where the linearized reaction rate is

This is obtained by linearizing the Langevin expression (19). One can easily show that the inclusion of the reaction term does not change $Sw,weq$, consistent with thermodynamic equilibrium. This explicitly confirms that our formulation is consistent with equilibrium statistical mechanics at the level of a Gaussian approximation of the fluctuations.

### D. Dilute limit

One of the common assumptions in traditional reaction-diffusion modeling is that each chemical species is dilute and thus diffuses independently of other species. In this section, we explain how our formulation simplifies in the dilute limit. We consider a solution where all solute species are dilute (i.e., *x*_{s} ≪ 1), but the solvent is possibly a homogeneous mixture. We use index *s* here to denote only solute species. In the dilute limit, *γ*_{s} → 1 and the solute number densities are linearly proportional to their mole fractions, $ns\u2248(\rho 0/m\xafsol)xs$, where $m\xafsol$ is the mixture-averaged molecular mass among solvent species; see (8). Hence, $\mu ^s$ can be expressed in terms of *n*_{s},

and consequently, the generalized (mole fraction based) LMA can be cast into the form of the traditional (number density based) LMA,

where $kr\xb1$ denote the reaction rate constants.

In the dilute limit, multispecies diffusion also becomes simpler. In Appendix A, we consider the dilute limit of a single solute species dissolved in a solvent mixture and show that the diffusion of the solute species is decoupled from solvent species; see (A8). It is straightforward to extend this result to multiple solute species. The diffusion coefficient of each solute species *s* then becomes a constant, that is, decoupled from the other species, yielding

where Ω_{s} represent the stochastic chemistry terms based on the LMA (35). Therefore, in the absence of fluid flow, our formulation is reduced to our previous reaction-diffusion model^{4} in the dilute limit. Note that jump processes and diffusion processes are combined in (36) and the time evolution of the probability distribution of ** n** = {

*n*

_{s}} can be described by the differential Chapman–Kolmogorov equation.

^{13}

## III. NUMERICAL METHOD

In developing a numerical method to solve (1)–(3), we seek an approach that

exhibits second-order accuracy in space and time deterministically and second-order weak accuracy in time for the linearized FHD equations;

^{35}reduces to our previous method for reaction-diffusion systems

^{4}in the dilute limit, in the absence of fluid flow;generates accurate structure factors for both equilibrium and giant fluctuations, even for large Schmidt numbers;

is robust in the presence of trace or vanishing species.

We explain below how our design decisions satisfy these requirements. In Sec. III A, we review our spatial discretization scheme and discuss robust numerical approaches for avoiding negative densities and treating vanishing species. In Sec. III B, we present our temporal integration scheme. In Sec. III C, we analyze the weak accuracy of our temporal integrator.

### A. Spatial discretization

Our spatial discretization is identical to the one used in our previous work on non-reactive FHD,^{3,26,27,36} with a few modifications noted below. The numerical framework is a structured-grid finite-volume approach with cell-averaged densities and pressure and face-averaged (staggered) velocities. We use standard second-order stencils for the gradient, divergence, and spatial averaging in order to satisfy discrete fluctuation-dissipation balance.^{24}

For the densities, we construct all mass fluxes on faces and employ the standard conservative divergence. For the advective mass fluxes, we implement two options. Centered advection uses two-point averaging of densities to faces, is nondissipative, and thus preserves the spectrum of fluctuations.^{24} However, in order to prevent unphysical oscillations in mass densities in high Péclet number flows with sharp gradients, we can also use the Bell–Dawson–Shubin (BDS) second-order Godunov advection scheme.^{37,38} We note that BDS advection adds artificial dissipation and does not obey a fluctuation-dissipation principle but is necessary for simulations where centered advection would fail due to insufficient spatial resolution. All simulations in this paper use centered advection unless otherwise noted. The discretization of the momentum equation is the same as our previous work.^{3,26} We allow for periodic boundary conditions, impermeable walls, and no-flow reservoirs^{27,36} held at fixed concentrations.

The first modification relative to our previous work^{3} is that for the stochastic mass fluxes $F\u0303$, we compute the matrix $2m\xaf\rho 0W\chi 12$ directly on the face using spatially averaged densities rather than computing this matrix at cell centers and averaging to faces. To compute the spatial averages, we use a *modified* arithmetic averaging function,^{4}

where *n*_{1} and *n*_{2} denote the number densities at the cell centers of two neighboring cells and *H* is a smoothed Heaviside function, defined as

Specifically, we first convert cell-centered mass fractions to number densities, then apply *ñ* to obtain face-centered number densities, and finally convert these back to mass fractions that are used to compute $m\xafW\chi 12$. We note that *ñ* drives the average (and thus stochastic flux) to zero if the number of molecules in either neighboring cell is sufficiently small (i.e., *n*_{i}Δ*V* ≤ 1), which prevents the occurrence of negative number densities. In most cases of interest, small numbers of molecules per cell correspond to dilute species. For dilute species [see (36)], the validity of using *ñ* has been justified in Ref. 4. In Sec. IV B, we numerically confirm that our approach is robust even when the total number of molecules in a cell is $O(10)$.

Another key modification is the computation of the diffusion matrix ** χ** for the deterministic and stochastic mass fluxes in the absence of some species or in the presence of vanishing species. In the vanishing limit, where one or more concentrations become zero, the diffusion matrix

**is not well conditioned since the corresponding diagonal component**

*χ**χ*

_{ss}diverges. This can cause numerical issues when one attempts to compute $F\xaf$ and $F\u0303$ since they depend on

**and $W\chi 12$, respectively. Unlike**

*Wχ***, however, the matrices**

*χ***and**

*Wχ***are well defined in the vanishing limit, and we can construct**

*WχW***using a special procedure. The basic idea is that we first compute a diffusion sub-matrix**

*Wχ*

*χ*^{sub}of

**with the rows and columns corresponding to each vanishing species omitted. Then we expand this sub-matrix into the full matrix**

*χ***and approximate the remaining components using the mathematical limit of vanishing species, $ws$ → 0**

*Wχ*^{+}, for all vanishing species

*s*.

To formally describe the procedure for computing ** Wχ** in the vanishing limit, we introduce a mapping, $m(i)$, used to expand/contract a subsystem matrix to/from a full matrix. For example, in a 6-species system having vanishing species $w2$ and $w4$, we have $m(i)=(1,0,2,0,3,4)$,

**= (1, …, 6). As graphically illustrated using the 6-species system in Fig. 1, there are four cases to consider when one populates (**

*i***)**

*Wχ*_{ij},

Note that color names in the parentheses in (39) correspond to the colors in Fig. 1. A derivation of (39) and (40) is presented in Appendix A. The full matrix $W\chi 12$ can be obtained from the Cholesky decomposition of the symmetric matrix ** WχW**. We note that if species

*s*is vanishing, then (

**)**

*WχW*_{is}= (

**)**

*WχW*_{sj}= 0, so no stochastic mass flux is generated for species

*s*.

We note that for each vanishing species only the diagonal element of ** Wχ** remains nonzero. Hence, the diffusion of a dilute species

*s*($ws$ ≪ 1) becomes decoupled from other species [see (A8)] and the effective diffusion coefficient

*D*

_{s}in (40) corresponds to the trace diffusion coefficient of

*s*in the given fluid mixture. It is also important to note that the construction (39) guarantees that $F\xafs=F\u0303s=0$ for vanishing species and ensures the mass conservation condition over all species, $\u2211s\u2032F\xafs\u2032=\u2211s\u2032F\u0303s\u2032=0$. Therefore, this procedure is robust to roundoff errors.

In our double-precision implementation, we treat any species *s* with $ws$ < 10^{−14} as a vanishing species.

### B. Temporal integration scheme

The spatial discretization of the non-reactive FHD equations for the mass densities yields a set of stochastic ordinary differential equations. It is straightforward to incorporate our CME-based chemistry model from Sec. II B via additional Poisson-noise terms,

where $ws,i$ denotes the mass fraction of species *s* in cell ** i**.

Our overall temporal integration strategy is a predictor-corrector approach for both species and velocity. Our goal is to develop a scheme that is second-order accurate in space and time deterministically, exhibits second-order weak accuracy in time for the linearized FHD equations, and treats reactions in a manner consistent with the CME.^{4} As explained below in detail, we treat viscous momentum dissipation implicitly and species diffusion explicitly. This is because in liquids the time step size is limited by the viscous Courant–Friedrichs–Lewy (CFL) condition (i.e., momentum diffusion) due to the large Schmidt number.

To combine a second-order midpoint tau-leap reaction sampling^{29,39} with a predictor-corrector scheme for FHD, we adopt mass density updates from the ExMidTau (explicit midpoint tau leaping) scheme we previously developed for reaction-diffusion systems.^{4} Hence, our new scheme uses a *midpoint* predictor for mass densities, which differs from our earlier *trapezoidal* scheme for non-reactive FHD systems.^{3,26} We have compared several combinations of mass and momentum updates to identify the variant of the scheme that gives the most accurate spectrum of the fluctuations (structure factors) in both equilibrium and giant fluctuation settings. In Sec. III C, we provide an analysis of the structure factors, which both guides and verifies the design decisions we made to ultimately choose this particular temporal integration scheme, and demonstrate the advantages of the midpoint scheme. One important observation is that our temporal integrator is robust in the large Schmidt number limit, Sc = *ν*/*D* → *∞*, where *ν* = *η*/*ρ*_{0}, unlike the trapezoidal scheme used in Refs. 3 and 26.

We advance the system from time *t*_{n} = *n*Δ*t* to time *t*_{n+1} = (*n* + 1)Δ*t* in four steps:

Perform a predictor Stokes solve for the velocity $\upsilon n+1,*$ at

*t*_{n+1}.Calculate predictor mass densities $\rho sn+12,*$ at the midpoint time $t=tn+12\Delta t$.

Calculate corrector mass densities $\rho sn+1$ at time

*t*_{n+1}.Perform a corrector Stokes solve for velocity $\upsilon n+1$ at

*t*_{n+1}.

These steps are elaborated in detail in Algorithm 1. In the Algorithm 1 description, superscripts are used to denote the time level where a given quantity is evaluated, e.g., *f*^{n} = ** f**($wn$). Also, $\mathcal{W}momn$ and $\mathcal{W}(i)massn$ (

*i*= 1, 2) denote the collections of i.i.d. (independent and identically distributed) standard normal random variables generated on control volume faces independently at each time step, and $\mathcal{W}\xafmom\u2261\mathcal{W}mom+\mathcal{W}momT$. We denote collections of independent Poisson random variables generated at cell centers independently at each time step with $P(i)$ (

*i*= 1, 2) and denote [•]

^{+}≡ max(•, 0).

In our time-advancement scheme, each Stokes problem couples a Crank–Nicolson discretization of viscous dissipation to the divergence-free constraint on velocity, to simultaneously solve for the velocity and mechanical pressure. To solve the Stokes system, we use a variable-coefficient (tensor) multigrid-preconditioned GMRES (generalized minimal residual) solver,^{28} as we have done previously.^{3,26} The difference between the predictor and corrector Stokes solves is the temporal discretization of the advective term (explicit vs. trapezoidal) and the forcing term (explicit vs. midpoint); both Stokes solves are required for second-order deterministic accuracy.

As mentioned above, Steps 2 and 3 of the present scheme become essentially the same as the ExMidTau scheme in the dilute limit in the absence of advection. The only difference is a Stratonovich-type update of the stochastic mass flux in (47). While our previous analysis for RDME systems^{4} adopted the Ito interpretation, we choose the Stratonovich-type update here since a general analysis for weak fluctuations (linearized FHD)^{40} guarantees second-order weak accuracy of the overall scheme for this choice. It can be shown that the Stratonovich and Ito interpretations become identical in the dilute limit. Hence, our numerical method not only achieves second-order weak accuracy for weak fluctuations but also inherits nice features of the ExMidTau scheme carefully designed for strong fluctuations.

### C. Structure factor analysis

We analyze our new temporal integrator by investigating time integration errors in the spectrum of giant concentration fluctuations for a binary mixture undergoing a dimerization reaction. We assume that a weak uniform concentration gradient is applied along the *y*-axis with gravity pointing in the positive *y*-direction. The Fourier-transformed linearized equations for *δ*$v\u2225$ ≡ *δ*$vy$ and *δ*$w$ (see Appendix C in Ref. 12) take the form

Here ** k** ≡

*k*_{⊥}is a wavevector in the plane perpendicular to the gradient,

*g*is the gravitational acceleration,

*ζ*=

*ρ*

^{−1}(

*∂ρ*/

*∂*$w$) is the solutal expansion coefficient, and

*h*is the concentration gradient,

**$w$ =**

*∇**h*

*e*_{y}. Using the method developed in Ref. 24, we analytically compute the resulting structure factors when our temporal integrator is used to solve (50a). For the non-reactive case (

*r*= 0), we also compute structure factors obtained from two schemes developed in our previous work.

^{3,26}The overdamped scheme (see Algorithm 2 in Ref. 26) uses the steady Stokes equation, i.e., it eliminates the inertial term

*∂*

_{t}$\upsilon $ =

**0**by taking an overdamped limit. We refer to the previous scheme for solving the inertial equation as the inertial trapezoidal scheme (see Algorithm 1 in Ref. 26) and to our new scheme as the inertial midpoint scheme (see Algorithm 1).

We set Δ*x* = 1 and *k*_{B}*T*/*ρ*_{0} = 1. To denote how fast momentum diffusion, species diffusion, and reaction are, we define the following dimensionless Courant numbers:

To consider the case of a relatively large Δ*t* with a large Schmidt number Sc = 10^{3} (as is typical of liquid mixtures), we set *α* = 250 and *β* = 0.25. For the reactive case, we consider two reaction rates *γ* = 0.025 and *γ* = 0.1, corresponding to penetration depths $\xi =D/r=10\Delta x$ and $\xi =1210\Delta x$, respectively. Other parameters are chosen so that $\mu w$*h*^{2}Δ*t*^{2} = 100, and *gζh*Δ*t*^{2} = 0.025 if gravity is present.

The structure factor can be decomposed into the sum $Sw,w=Sw,weq+Sw,wneq$, where $Sw,weq$ is the equilibrium structure factor (31) and $Sw,wneq$ is the nonequilibrium enhancement. In the non-reactive case with no gravity, the nonequilibrium enhancement exhibits a *k*^{−4} power law in the entire range of wavenumbers *k*,

For the non-reactive case with gravity, we compare $Sw,w$ obtained from the three schemes with the exact result in Fig. 2(a). A power-law spectrum $Sw,w$ ∼ *k*^{−4} develops for intermediate wavenumbers *k*. At small wavenumbers, $Sw,w$ becomes constant due to gravity. At large wavenumbers, the *k*^{−4} decay in the nonequilibrium part is hidden due to the flat equilibrium structure factor $Sw,weq$. Our numerical scheme reproduces $Sw,w$ accurately for all but the largest *k* values, whereas both earlier schemes exhibit significant deviations at either large or small *k* values. Significant deviations of the previous inertial scheme at large *k* are due to temporal integration errors in the nonequilibrium part $Sw,wneq$, as can be seen more clearly by examining the cross correlation between fluctuations of $w$ and $v$ (not shown). The divergence of $Sw,w$ for the overdamped scheme at small *k* demonstrates that the overdamped limit does not apply for sufficiently small *k* with gravity. Thus, our new scheme combines the favorable features of our previous trapezoidal inertial scheme (correct behavior for small *k* with gravity) and the overdamped scheme (correct behavior for large *k*).

In Fig. 2(b), we show the nonequilibrium enhancement $Sw,wneq$ in the structure factor for the reactive case with no gravity. We obtain the exact $Sw,wneq$ by analyzing (50a) without the stochastic mass fluxes and with deterministic reaction [see also Eq. (58) in Ref. 22 or Eq. (44) in Ref. 23],

Our midpoint scheme reproduces $Sw,wneq$ accurately for both rate constants. We emphasize that these results are remarkable given that *α* = O(10^{2}). Our new scheme remains accurate for *α* = *β*Sc ≫ 1 because the relative error in $Sw,wneq$ for our midpoint scheme has the form $[1+O(Sc\u22121)]O(\Delta t2)$, indicating robust behavior for large Schmidt numbers. On the other hand, the relative error for the trapezoidal scheme has an $O(Sc)O(\Delta t4)$ term, which results in significant deviations at large *k*, as observed in Fig. 2(a).

## IV. NUMERICAL EXAMPLES

In this section, we consider four examples that demonstrate the capabilities of our numerical methodology. In Sec. IV A, we model the hydrolysis of sucrose in an aqueous solution with very dilute solutes. In Sec. IV B, we investigate a binary liquid mixture undergoing a dimerization reaction at thermodynamic equilibrium. In Sec. IV C, to verify the correct coupling of mass and momentum fluctuations, we study nonequilibrium giant fluctuations in a mixture undergoing a dimerization reaction. In Sec. IV D, to demonstrate the scalability and practical utility of our method, we investigate the effects of fluctuations for a reactive fingering instability.

### A. Hydrolysis of sucrose

We consider a dilute solution of sugar in water at equilibrium, undergoing the reversible hydrolysis reaction

Sucrose is particularly dilute, with only ∼10 molecules per computational cell, whereas there are ∼10^{7} glucose and fructose molecules and ∼10^{10} water molecules per cell. We investigate the equilibrium distribution of the number of sucrose molecules in a cell to demonstrate that our approach correctly models the dilute limit.

We use cgs units and choose physical parameters assuming *T* = 293, atmospheric pressure, *ρ*_{0} = 1, and *η* = 0.01. The four species are glucose (*s* = 1), fructose (*s* = 2), sucrose (*s* = 3), and water (*s* = 4). Using the trace diffusion coefficients of the solutes, *D*_{s} (*s* = 1, 2, 3),^{41,42} and the self-diffusion coefficient of water *D*_{water},^{43} the Maxwell–Stefan binary diffusion coefficients are assigned as in Ref. 44,

Since we consider the dilute limit, we assume that the system is an ideal mixture and obeys the traditional LMA, with a forward rate *a*^{+} = *k*^{+}*n*_{3}, a reverse rate *a*^{−} = *k*^{−}*n*_{1}*n*_{2}, and an equilibrium constant $K=n1eqn2eq/n3eq$. The reaction equilibrium lies almost completely in the direction of the formation of glucose and fructose,^{45} but uncatalyzed sucrose hydrolysis is extremely slow (with a half-life of 500 years).^{46} While we use an experimental value of *K*,^{45} we artificially increase the reaction rates to *k*^{+} = 10 and *k*^{−} = *K*/*k*^{+} so that the forward reaction occurs about 100 times per cell per simulation.

We set up a two-dimensional system consisting of 32 × 32 cells with dimensions Δ*x* = Δ*y* = 10^{−4} and periodic boundary conditions. The thickness of the system is Δ*z* = 10^{−4} and the cell volume Δ*V* = Δ*x*Δ*y*Δ*z*. We consider the case where there are ten sucrose molecules per cell. Hence, $n3eq$ is determined from $n3eq\Delta V=10$ and $n1eq=n2eq$ are subsequently determined from equilibrium. The resulting equilibrium mass fractions are $w1eq=w2eq=4.9\u2009\xd7\u200910\u221203$, $w3eq=5.7\u2009\xd7\u200910\u221209$, and $w4eq=0.990$. We use two time step sizes, Δ*t* = 10^{−5} and 10^{−4}, to check the continuous-time limit and quantify time integration errors. Note that the larger Δ*t* corresponds to diffusive CFL numbers *D*_{s,max}Δ*t*/Δ*x*^{2} = 0.07 for species diffusion and *νΔt*/Δ*x*^{2} = 100 for momentum diffusion. For each value of Δ*t*, we ran 16 independent samples up to time $T=1$, collecting data every *t* = 10^{−4} for $t\u2265T/10$.

We recall that the number of sucrose molecules in a cell *N* = *n*_{3}Δ*V* has a continuous range in FHD simulations. We define its discrete distribution as $P(N)=\u222bN\u221212N+12\rho (N\u2032)dN\u2032$, where *ρ*(*N*) is the continuous distribution of *N*. Figure 3 shows that for the smaller Δ*t*, *P*(*N*) is remarkably close to the physically correct Poisson distribution *P*_{Poisson}(*N*), and *ρ*(*N*) is essentially zero for negative values of *N*. We note that *P*_{Poisson}(*N*) is significantly different from its Gaussian approximation *P*_{Gauss}(*N*). For the larger Δ*t*, while the remarkable agreement with the Poisson distribution is still observed, negative values of *N* start to appear, yielding $\u222b\u2212\u221e0\rho (N)dN\u22483\u2009\xd7\u200910\u22125$; see the inset of Fig. 3. The same results were obtained in our previous reaction-diffusion model of dilute solutions,^{4} confirming that our treatment of the stochastic mass flux coefficients (see Sec. III A) is consistent with the dilute limit, even in the presence of random advection. We have also confirmed that the equilibrium structure factor of each species (not shown) has a flat spectrum (as predicted by theory^{3}), indicating that there are no spurious correlations between cells.

### B. Dimerization: Equilibrium distribution

We next consider a liquid mixture undergoing the dimerization reaction (20). This binary system contains monomers A (*s* = 1) and dimers A_{2} (*s* = 2) and is representative of cyclic dimer formation in pure liquid acetic acid. We demonstrate here our ability to model a system with strong fluctuations in the absence of a dominant solvent by considering a small number of molecules (∼10) of each species per cell. As in the sugar solution example, we investigate the equilibrium distribution of monomers and dimers; however, since the system is not dilute, the distribution of each species is not Poisson. The numbers of monomers and dimers in a cell (*N*_{1} and *N*_{2}) do not vary independently due to the constant density assumption *N*_{1} + 2*N*_{2} = *ρ*_{0}Δ*V*/*m*, where *m* is the mass of a monomer. Therefore, we investigate the equilibrium distribution of *N*_{2}, *P*(*N*_{2}), for *N*_{1} + 2*N*_{2} = 40.

We simulate a two-dimensional system consisting of 32 × 32 cells under periodic boundary conditions. Here we use arbitrary units that give Δ*x* = Δ*y* = Δ*z* = 1, $\u2212D12=1$, and *m* = 1 with *k*_{B} = 1. We set *ρ*_{0} = 40 with $w1eq=w2eq=0.5$ so that *N*_{1} + 2*N*_{2} = 40 with $N1eq=20$ and $N2eq=10$. We set the reaction rates *a*^{±} as in (26) with a modification

where *N*^{+} = max(*N*, 0). Note that (56) turns off unphysical reactions when 0 < *N*_{1} < 1. The rate constants *κ*^{+} = 0.8724 and *κ*^{−} = 1.125 are chosen as follows. The ratio *K* = *κ*^{+}/*κ*^{−} = 0.7755 is determined so that the resulting theoretical distribution gives $\u27e8N2\u27e9=\u2211N2N2P(N2)=N2eq$. The magnitude of *κ*^{±} is determined so that the linearized reaction rate *r* = 0.1 [see (33)] gives a penetration depth $\xi \u2261\u2212D12/r=10\Delta x$. We set *η* = 10^{3} and $T=103$. We use a small Δ*t* = 10^{−2} to minimize temporal integration errors. For 16 independent samples with 10^{5} time steps, we collect data every 10^{2} time steps, discarding the first 10^{4} time steps.

In Fig. 4, we compare the simulation result for the equilibrium distribution *P*(*N*_{2}) with theoretical results obtained in Sec. III A. We denote the exact Einstein distribution obtained from the entropy expression (22) by *P*_{exact}, Stirling’s approximation result obtained from (27) by *P*_{Stirling} ∼ exp(*S*_{Stirling}/*k*_{B}), and the Gaussian approximation of *P*_{Stirling} by *P*_{Gauss}. Note that *P*_{exact} is a discrete distribution, whereas *P*_{Stirling} and *P*_{Gauss} have continuous ranges, 0 < *N*_{2} < 20 and −*∞* < *N*_{2} < *∞*, respectively. Significant deviations of *P*_{Gauss} from *P*_{exact} indicate that the system is subject to strong fluctuations, as expected from the small number $N2eq=10$. Over a remarkably wide range, our numerical method accurately matches *P*_{exact}. Even beyond this range, it gives sensible values with accuracy better than or comparable to *P*_{Stirling}. Measurable deviations are observed only for larger values *N*_{2} = 19 and 20, for which the occupation probabilities are already very small (*P*_{exact}(*N*_{2}) < 10^{−6}).

It is important to note that statistically identical results for *P*(*N*_{2}) are obtained from the corresponding non-reactive system with *κ*^{±} = 0 (not shown). This confirms thermodynamic consistency of our overall formulation. In addition, this also confirms the validity of our overall numerical treatment for diffusion with strong fluctuations. In particular, considering that our multiplicative GWN modeling for strong fluctuations was developed in the dilute limit^{4} and analyzed only for this case, the validity of its extension to non-dilute solutions cannot be taken for granted.

### C. Dimerization: Giant fluctuations

We now investigate a system where velocity fluctuations are coupled to diffusion. We consider the same dimerization reaction but examine giant fluctuations in the presence of a weak concentration gradient with no gravity. We focus on the nonequilibrium contribution to the structure factor, $Sw1,w1neq=Sw1,w1\u2212Sw1,w1eq$, for wavevectors perpendicular to the concentration gradient. We neglect stochastic mass fluxes and use deterministic chemistry so that we eliminate the equilibrium fluctuations and thus obtain $Sw1,w1neq$ with greater statistical accuracy. We have previously considered a gas mixture in a similar setting;^{12} here we consider a liquid mixture with a large Schmidt number Sc = 10^{3}, which quantitatively changes the spectrum of giant fluctuations.

A detailed theoretical analysis of giant fluctuations using linearized FHD first appeared in Ref. 22 assuming that the system is near chemical equilibrium everywhere. It was later extended in Ref. 23 to account for the nonlinearity caused by the fact that the system is not everywhere in chemical equilibrium; this theoretical analysis assumes a liquid mixture, so it was not applicable for the gas mixture we considered in Ref. 12. In these theoretical studies, the concentration gradient was imposed via the Soret effect by applying a temperature gradient, unlike the case we consider here where the concentration is fixed at the *y*-walls using reservoir boundary conditions. Furthermore, the theoretical studies in Refs. 12 and 22 did not account for the boundary conditions for the fluctuating fields.

We consider a two-dimensional square domain of side length *L*_{x} = *L*_{y} = 64 (in arbitrary units) and periodic boundary conditions in the *x*-direction. The system is divided into 128 × 128 grid cells with grid spacing Δ*x* = Δ*y* = 0.5. To remain in the linearized FHD regime, we increase the cell depth to Δ*z* = 10^{6} so that there are sufficiently many monomers and dimers in a cell, *N*_{1} + 2*N*_{2} = 2.5 × 10^{5} for *ρ*_{0} = 1 and *m* = 1. We set $\u2212D12=1$, *η* = 10^{3}, and *k*_{B}*T* = 10^{3}. For the dimerization reaction, the equilibrium constant *K* = *κ*^{+}/*κ*^{−} = 0.75 is chosen to give a reference equilibrium state with $w1eq=w2eq=0.5$. Two sets of reaction constants are considered: (*κ*^{+}, *κ*^{−}) = (8.438 × 10^{−2}, 0.1125), corresponding to linearized reaction rate *r* = 0.4 and penetration depth $\xi =10\Delta x$, and (*κ*^{+}, *κ*^{−}) = (8.438 × 10^{−3}, 1.125 × 10^{−2}), corresponding to *r* = 0.04 and *ξ* = 10Δ*x*. The time step size is set to Δ*t* = 0.025, which gives Courant numbers $\u2212D12\Delta t/\Delta x2=0.1$ and *νΔt*/Δ*x*^{2} = 100. We ran 10^{5} steps discarding the first 10^{4} steps and computed the steady-state monomer concentration profile $w\xaf1(y)$ and $Sw1,w1neq(kx)$.

To impose a concentration gradient in the *y*-direction, no-slip reservoir conditions^{27} are imposed with ($w1$, $w2$) = (0.49, 0.51) at *y* = 0 and ($w1$, $w2$) = (0.51, 0.49) at *y* = *L*_{y}. While a linear concentration profile is formed in the steady state for the non-reactive case, a nonlinear profile is generated by the dimerization reaction. In Fig. 5(a), the profiles of $w\xaf1(y)$ for reaction rates *r* = 0.4 and 0.04 are compared with the one for the non-reactive case. As *r* increases, the nonlinearity in $w\xaf1(y)$ becomes more evident. This is because a larger region around *y* = *L*_{y}/2 is constrained to be in chemical equilibrium due to faster reactions, resulting in larger concentration gradients at the boundaries. Identical concentration profiles are obtained from the corresponding deterministic reaction-diffusion systems (not shown).

In Fig. 5(b), we show numerical results of $Sw1,w1neq$. To account for errors in the discrete approximation to the continuum Laplacian, the modified wavenumber^{36}

is used instead of *k*_{x}. For the non-reactive case (*r* = 0), a clear *k*^{−4} power law is observed until the confinement effect becomes significant for small $kx\u226aLy\u22121$. For the reactive cases, the *k*^{−4} power law is only observed at large $kx\u226br/\u2212D12$. For larger *r*, the *k*^{−4} power law appears in a narrower range of *k*_{x} values and the prefactor of the power law becomes larger.

For the non-reactive case, the prefactor of the power law is accurately predicted by the theoretical prediction (53). By multiplying (53) by the confinement factor^{47}

the theoretical prediction is further improved at small *k*_{x}, as shown in Fig. 5(b). We note, however, that this factor is obtained for impermeable walls and the resulting correction is not exact for our reservoir boundaries. For the reactive cases, the validity of (53) is questionable due to the nonlinear concentration gradients. In fact, how to estimate the value of *h*^{2} is not obvious. Considering that $Sw1,w1neq$ is an averaged structure factor for different values of *y*, we estimate the value of *h*^{2} from the profile of $w\xaf1(y)$ using a spatial average,

Theoretical results obtained from (53), (58), and (59) are shown in Fig. 5(b). Remarkably, the *k*^{−4} power law region is accurately predicted. However, the theoretical prediction overestimates $Sw1,w1neq$ at small *k*_{x} by several orders of magnitude for the reactive cases. This is expected since the local linear gradient approximation eventually fails at large length scales. The FHD equations linearized around a nonlinear stationary profile were studied in Ref. 23; however, an explicit result for the static structure factor that we could compare with our numerical result was not obtained.

### D. Fingering instability with a neutralization reaction

In this section, we examine the development of asymmetric fingering patterns arising from a diffusion-driven gravitational instability in the presence of a neutralization reaction. We perform three-dimensional large-scale simulations of a double-diffusive instability occurring during the mixing of HCl and NaOH solution layers in a vertical Hele-Shaw cell (two parallel glass plates separated by a narrow gap). This system has been studied experimentally and theoretically using a two-dimensional Darcy advection-diffusion-reaction model.^{2,48} Thermal fluctuations may play a key role in triggering the instability. To the best of our knowledge, our simulations are the first ones to use a three-dimensional model and the first to include thermal fluctuations. We investigate the effects of each stochastic component (mass flux, momentum flux, and chemistry) on the evolution of the system. We initialize our simulations with natural mass and momentum fluctuations without any artificial perturbation, and therefore our simulation can be regarded as an *ideal* experiment.

#### 1. Model setup

For the model setup and physical parameters, we follow the experiment of Lemaigre *et al.*^{2} We use cgs units unless otherwise specified and assume *T* = 293 and atmospheric pressure. The isothermal approximation has been justified by a linear stability analysis showing that the heat generated by the neutralization reaction

plays a negligible role in this problem.^{48} We consider a Hele-Shaw cell with side lengths *L*_{x} = *L*_{y} = 1.6 and *L*_{z} = 0.05, with the *y*-axis pointing in the vertical direction, and the *z*-axis being perpendicular to the glass plates. The domain is divided into grid cells with grid spacing Δ*x* = Δ*y* = Δ*z* = 6.25 × 10^{−3}, so there are 256 × 256 × 8 cells. We impose periodic boundary conditions in the *x*-direction and no-slip walls in the *z*-direction. In the *y*-direction, we impose free-slip reservoir boundary conditions with concentrations that match the initial conditions of each layer.

We start with a gravitationally stable initial configuration, where an aqueous solution of NaOH with a molarity of 0.4 mol/l is placed on top of a denser aqueous solution of HCl with a molarity of 1 mol/l. Each reactant and product is treated as a single charge-neutral species, giving the four species HCl (*s* = 1), NaOH (*s* = 2), NaCl (*s* = 3), and water (*s* = 4). Under the approximation that the solution density *ρ* is linearly dependent on the solute concentrations,^{2} the buoyancy force is expressed as

where *α*_{s} is the solutal expansion coefficient and *M*_{s} is the molecular weight (in g/mol) of the solute *s*. We set *g* = 981, *ρ*_{0} = 1, and *η* = 0.01. The initial density difference between the two layers is approximately 4 × 10^{−4}. The Maxwell–Stefan binary diffusion coefficients are determined using (55) from the known trace diffusion coefficients of the solutes and the self-diffusion coefficient of water. The values of *α*_{s} and the trace diffusion coefficients are obtained from Table II of Ref. 2.

Since the neutralization equilibrium lies far to the product side, we only consider the forward reaction. We assume that the rate is given by the traditional LMA for a dilute solution, *a*^{+} = *kn*_{1}*n*_{2}. However, we note that neutralization is a diffusion-limited reaction. In other words, reaction occurs extremely fast (with rate *λ* ∼ 10^{11} s^{−1}), as soon as reactants encounter each other. Because of this, the validity of the local LMA is questionable.^{5} The estimated value of *k* ∼ 10^{−11} cm^{3} s^{−1} is impractically large [converted using (31) in Ref. 49] and would require an unreasonably small Δ*t* for numerical stability. For our simulations, we choose a smaller value *k* = 10^{−18} and Δ*t* = 10^{−3} based on a deterministic numerical study presented in Appendix B.

The initial mass fractions in each cell are generated as the sum of mean values $w0$ and natural fluctuations *δ*$w$. The mean mass fractions $w0$ are set to $w0,upper$ = (0, 0.0157, 0, 0.9843) in the upper half-domain and $w0,lower$ = (0.0358, 0, 0, 0.9642) in the lower half-domain. Assuming that natural fluctuations are Gaussian, we sample them using the known equilibrium structure factor at the mean state [Eq. (D4) in Ref. 3],

where $W0=diag(ws0)$, ** M** = diag(

*m*

_{s}), and

*z*^{mass}is a vector with i.i.d. standard normal random variables. The initial momentum fluctuations are generated in a similar manner,

where *z*^{mom} is a vector with i.i.d. standard normal random variables, followed by an *L*^{2} projection onto the space of divergence-free vector fields.

We use the Langevin chemistry description given in (19) and the BDS advection scheme. We can justify the use of the CLE by noting that the system is near the macroscopic limit because fluctuations are weak. For centered advection, we observe oscillations around the interface of fingers (not shown) for the chosen grid spacing as expected due to the high cell Péclet number. When the grid spacing is reduced to half, oscillations become less pronounced without changing the results significantly (not shown).

#### 2. Effects of thermal fluctuations

We perform four FHD simulations changing how chemistry is treated, whether natural mass/momentum fluctuations are initially imposed, and whether subsequently stochastic mass/momentum fluxes are included, as summarized in Table I.

. | . | Initial fluctuations . | Stochastic fluxes . | ||
---|---|---|---|---|---|

. | Chemistry . | Mass . | Momentum . | Mass . | Momentum . |

Simulation A | Stochastic | On | On | On | On |

Simulation B | No reaction | On | On | On | On |

Simulation C | Deterministic | Off | On | Off | On |

Simulation D | Deterministic | Off | Off | Off | On |

. | . | Initial fluctuations . | Stochastic fluxes . | ||
---|---|---|---|---|---|

. | Chemistry . | Mass . | Momentum . | Mass . | Momentum . |

Simulation A | Stochastic | On | On | On | On |

Simulation B | No reaction | On | On | On | On |

Simulation C | Deterministic | Off | On | Off | On |

Simulation D | Deterministic | Off | Off | Off | On |

By comparing the results of simulations A, B, C, and D, we can assess the effects of chemo-hydrodynamic coupling and thermal fluctuations on the fingering pattern formation. For a perfectly flat initial interface, thermal fluctuations play an essential role in perturbing the interface at early times, but once an uneven interface appears, the dynamic instability dominates and thermal fluctuations are expected to play a secondary role in subsequent pattern formation, as we previously confirmed in the absence of reactions.^{3}

We compare the reactive case (simulation A) with the non-reactive case (simulation B) in Fig. 6. As also seen in the experiment,^{2} an asymmetric growth of fingers is observed in the reactive case. In addition, the growth of fingers is much faster when the reaction is present. This is due to the coupling of the fast neutralization reaction and the fast diffusion of the acid species. Disparities between the acid and base species can also be seen in the concentration profiles around the fingers; such disparities are not observed in the non-reactive case. We point out that the concentration develops three-dimensional profiles that are not constant across the thickness of the Hele-Shaw cell, as can be seen from the side bars (*y*-*z* cross sections) in the figure. Such a structure would not be captured by the two-dimensional Darcy approximation used in prior computational studies.^{2,48}

In Fig. 7, we compare simulations C and D with simulation A to investigate the contribution of each stochastic component. Compared with the full fluctuating hydrodynamics (simulation A), all stochastic *mass* components are omitted in simulation C. However, the resulting fingering patterns are essentially the same as in simulation A. This indicates that contributions of stochastic mass fluxes and stochastic chemistry are negligible in this example. Instead, concentration fluctuations driven by the stochastic momentum flux dominate the formation of an uneven interface. This can be confirmed by the comparison of simulation A with simulation D, where initial velocity fluctuations are turned off compared with simulation C, and only stochastic momentum fluxes are included. The resulting fingering patterns are virtually the same with slight differences caused by differences in the initial velocity conditions. This is consistent with the fact that any initial momentum conditions are quickly damped out in a liquid with a high Schmidt number. In fact, in a simulation similar to simulation C but *without* subsequent stochastic momentum fluxes, it takes more time (∼5 s longer) for fingering patterns to start to grow. Hence, we confirm that velocity fluctuations driving giant composition fluctuations dominate the triggering of the instability starting from a perfectly flat interface.

It is important to note that our simulation results show that the thermal fluctuations are sufficiently large to kick off the instability on a time scale comparable to that when a macroscopic initial perturbation is imposed. The fingering patterns observed in simulation A at *t* = 40 s are quite comparable to the experimental result shown in Fig. 1(e) of Ref. 2 at *t* = 30 s. Of course, in the actual experiments, the initial interface is not perfectly flat due to imperfections in the preparation of the initial configuration.

## V. SUMMARY AND DISCUSSION

We have developed a fluctuating hydrodynamics (FHD) formulation and numerical methodology for stochastic simulation of reactive liquid mixtures. Our approach robustly models a wide range of microliquids, including dilute solutions as well as mixtures with no single dominant solvent. Our multispecies transport model is based on the Maxwell–Stefan cross-diffusion, incorporates a stochastic chemistry description based on the chemical master equation (CME), and couples reaction-diffusion with a stochastic Navier–Stokes equation for the fluid velocity. Our numerical method is based on several techniques that helped us gain computational efficiency without compromising accuracy. Specifically, the implicit treatment of momentum dissipation allowed us to avoid the severe restriction on time step size imposed by the small Reynolds number and large Schmidt number. The use of the tau leaping method enabled us to sample CME-based chemistry efficiently while correctly sampling large deviations from chemical equilibrium.^{34} For a binary liquid mixture undergoing a dimerization reaction, we demonstrated the thermodynamic consistency of our overall formulation beyond the Gaussian approximation and accurately reproduced the equilibrium Einstein distribution for both dilute solutions and liquid mixtures. Owing to a careful treatment of strong concentration fluctuations and vanishing species, our numerical method remained robust even for cells with as few as ten molecules; coarse-graining at such small scales is at the limits of fluctuating hydrodynamics.

Our numerical results for the spectrum of giant nonequilibrium fluctuations in a binary mixture undergoing a dimerization reaction were not in good agreement with theoretical predictions for smaller wavenumbers. We believe that this mismatch is due to the fact that we used a very simple theory which assumes that the gradient is constant and weak. A more accurate theory requires linearizing the FHD equations around the nonlinear steady-state solution of the macroscopic equations and proper treatment of the boundary conditions. Such a linearization was carried out in Ref. 23 without accounting for the boundary conditions [see in particular Eq. (15) of Ref. 23]. Nevertheless, analytical computation of the structure factor proved too difficult and the authors used a perturbative analysis for which the zeroth-order approximation is the simple approximation that the applied gradient is constant and weak and the system is everywhere near chemical equilibrium. An explicit formula for the next-order correction was not obtained for the static structure factor. Our computations showed that the simple zeroth-order theory, while giving a qualitatively correct picture of how reactions affect giant fluctuations, overestimates the fluctuations by orders of magnitude at small wavenumbers.

We performed the first three-dimensional simulations of a buoyancy-driven instability in the presence of an acid-base neutralization reaction. Our results demonstrate that velocity fluctuations generate giant concentration fluctuations that are sufficiently large to drive the initial growth of the instability, even when the initial interface is perfectly flat. In particular, we found that thermal fluctuations can trigger the instability on a time scale comparable to that observed in recent experiments, although a direct comparison is not possible because the exact initial condition in experiments is hard to control or measure.

In our prior work on reaction-diffusion systems,^{4} we treated diffusion implicitly. This allowed us to use time step sizes an order of magnitude larger than the stability limit imposed by species diffusion. In this work, we treated diffusion explicitly because momentum diffusion is much faster than mass diffusion in liquids, and thus the time step size was primarily limited by the requirement to accurately resolve momentum dynamics at small scales. Nevertheless, in a number of problems, such as catalytic micropumps^{50} or electroconvective instabilities at large applied voltages,^{51} the time scales of interest are those at which diffusion reaches a quasi-steady state in at least one direction. In this case, one must treat diffusion implicitly. This is straightforward in principle but requires the development of several nontrivial components. First, because the diffusion of all species is coupled in generic mixtures, one must develop either temporal integrators that treat only the diagonal part of the diffusion matrix implicity or develop a multispecies multigrid solver for coupled implicit discretizations. Second, the semi-implicit temporal integrators developed in Ref. 4 must be modified to integrate the momentum equation in a way that is robust for large Schmidt numbers. Third, boundary conditions need to be handled, both in the diffusion solver and in the coupling between diffusion and advection for reservoir boundaries.

In this work, we assumed the validity of a Boussinesq approximation, neglecting the change in density with composition at a given pressure and temperature, as dictated by the equation of state (EOS) of the mixture. This is a limiting approximation in practice, especially for reactive mixtures in which reactions can rapidly change density locally. In prior work,^{3} we accounted for the density dependence on composition using low Mach asymptotics. It is important to observe that the multispecies low Mach model proposed in Ref. 3 applies even to ideal gas mixtures, not just liquid mixtures. There are several difficulties with extending the formulation and algorithms we developed in prior work to reactive low Mach number models. First, reactions can lead to local changes in pressure which, in the low Mach limit, must get instantaneously distributed throughout the system as a global adjustment of the background thermodynamic pressure. It is anticipated that barodiffusion will have to be accounted for to achieve thermodynamic consistency when the chemical potentials depend nontrivially on pressures. Second, enforcing the EOS will require a nonlinear iteration of a coupled mass-momentum diffusion system, unlike the simpler case considered in Ref. 3 where we could enforce a linear EOS with only a decoupled linear Stokes solve. Both of these difficulties are compounded if one wishes to treat diffusion implicitly or to account for energy transport in a non-isothermal generalization.

In future work, we will account for electrochemistry by incorporating charged species into our model, similar to the developments in Ref. 52 for the non-reactive case. By using electroneutral asymptotics,^{53} we will be able to model the species (HCl, NaOH, and NaCl) in the acid-base fingering instability as separate ions (H^{+}, OH^{−}, Na^{+}, and Cl^{−}), which is physically correct given that HCl and NaOH are both strong electrolytes. Resolving the diffusion of each ion individually is required to correctly model electrodiffusion in mixtures with more than two ions. Incorporating charged species will also allow us to simulate weak electrolyte solutions (in which molecules do not fully disassociate into ions), catalytic motors,^{50} and electrokinetic locomotion.^{54,55}

## ACKNOWLEDGMENTS

We would like to thank Anne De Wit for helpful discussions regarding gravitational instabilities in the presence of neutralization reactions and Eric Vanden-Eijnden for discussions regarding tau leaping and large deviation theory. This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under Contract No. DE-AC02-05CH11231 and Award No. DE-SC0008271. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

### APPENDIX A: DIFFUSION MATRIX WITH VANISHING SPECIES

In this appendix, we derive the analytic expressions (39) for the matrix ** Wχ** in the limit of vanishing species. For simplicity, we consider a case where the last species among

*N*species vanishes,

We show next that each component of ** Wχ** converges to

where *χ*^{sub} is the diffusion matrix of the subsystem consisting of non-vanishing species with $w0=(w10,\u2026,wN\u221210)$, and *x*^{0} and $m\xaf0$ are computed from (7) and (8) using $w0$. We also show that the trace diffusion coefficient *D*_{N} of species *N* in the fluid mixture with composition $w0$ is expressed as

By rearranging the definition of ** χ**,

^{3}

where *α* ≠ 0, and using **1**^{T}**Λ** = **0**, ** χ**$w$ =

**0**, and

**1**

^{T}$w$ = 1, we obtain

Noting that $\chi NN=O(wN\u22121)$ and

we obtain (A2b) and (A3) by taking the limit of (A6) using (7). Then (A2b) and (A7) imply (A2c). Applying the same technique to the (*i*, *N*)-component of (A5) for *i* = 1, …, *N* − 1, we obtain (A2d).

We also observe that (A2) gives

Hence, the diffusion of the dilute species becomes decoupled from those of other species and its trace diffusion coefficient is given by (A3).

### APPENDIX B: RATE CONSTANT OF NEUTRALIZATION REACTION

In this appendix, we determine an appropriate value for the rate constant *k* of the neutralization reaction (60) for the simulations of the fingering example reported in Sec. IV D. As mentioned in the main text, the estimated value of *k* ∼ 10^{−11} is too large, as it requires impractically small time step sizes. By performing deterministic simulations, we investigate a range of values for *k* to determine at what point increasing *k* stops changing the results. We also examine the convergence of the results using different time step sizes.

For these deterministic simulations, we use a smaller domain (half the length in the *x*- and *y*-directions) with the same grid spacing. To generate an initial configuration with an uneven interface, we introduce random perturbations of composition in each cell immediately above the interface and set

where *a* = 0.1 and *U* is a standard normal random number generated independently in each cell. We compute fingering patterns for several values of *k* from 10^{−23} to 10^{−15}, with several values of Δ*t* ranging from 10^{−3} to 10^{−2}, using the *same* random initial configuration. To assess the similarity of two simulation results, we compute the gross NaCl production *ρ*_{0}*∫*$w3$(** r**,

*t*)

*d*

**, as well as the**

*r**L*

^{1}-norm of the $vy$ field ∥$vy$∥ =

*∫*|$vy$(

**,**

*r**t*)|

*d*

**.**

*r*Figure 8(a) shows the time evolution of ∥$vy$∥ for various values of *k* for Δ*t* = 10^{−3}. As *k* increases, ∥$vy$∥ grows faster, indicating that fingers grow faster. For 10^{−22} ≲ *k* ≲ 10^{−19}, time profiles change significantly depending on the value of *k*. On the other hand, for *k* ≳ 10^{−19}, the change becomes less significant. Also, time profiles for *k* ≲ 10^{−22} coincide with that of the non-reactive case. This suggests that there are three different regimes for *k*: slow, intermediate, and fast reaction regimes. The gross NaCl production shown in Fig. 8(b) exhibits similar behavior. While more NaCl is produced as *k* increases, the growth slows down around *k* ≈ 10^{−19} and a plateau is observed beyond this value. Hence, from a modeling point of view, one can simulate the neutralization reaction using a value of *k* from the plateau region. It is important to note, however, that one cannot choose an arbitrarily large value of *k* due to the stability limit imposed by our explicit tau-leaping treatment of reactions. In fact, fingering patterns obtained using Δ*t* = 10^{−2} and 10^{−3} (not shown) are essentially the same for *k* ≲ 10^{−18}. However, both results start to show unphysical behaviors for $k\Delta t>4/nHCl0$, where $nHCl0$ is the initial number density of HCl species in the lower layer, as can be seen from the abrupt increase of the gross NaCl production in Fig. 8(b).

Based on these observations, we choose *k* = 10^{−18} and Δ*t* = 10^{−3}. The value of Δ*t* is much smaller than the mass diffusion stability limit. As shown in Fig. 8(b), Δ*t* ≲ 10^{−2} is required to guarantee stability when the reaction is stiff and *k* ≈ 10^{−18}. It is noted, however, that Δ*t* ≲ 10^{−3} is required to give a reasonable CFL number for momentum diffusion *νΔt*/Δ*x*^{2} = 0.256. This is because small time-integration errors in the velocity field at early times can cause significant perturbations at later times because of the growing instability. If the exact time evolution at early times is not important, one can safely use Δ*t* = 10^{−2} without sacrificing physical fidelity.

## REFERENCES

In a reaction-diffusion setting, this means that diffusion dominates on the length scale Δ*x* of a reactive cell. Equivalently, for typical diffusion coefficient *D* and linearized reaction rate *r*, the penetration depth $\xi =D/r$ is significantly larger than Δ*x*, which is itself much larger than molecular scales. In this reaction-limited case, the validity of the mesoscopic description of reactions is guaranteed. On the other hand, for a diffusion-limited system, where *ξ* becomes comparable to molecular scales, the validity of the CME remains to be investigated.^{5}