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.

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 scale4 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 systems4 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 Keizer19), 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 work3 but with some important improvements necessary for simulating complex reactive mixtures at small length scales. Notably, we extend our previous work on reaction-diffusion systems4 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 mixtures25 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 dissipation26 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 scheme29 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 systems4 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 A2 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 experimentally2 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.

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,

ρ0υt+π=ρ0(υυT)+(η¯υ+Σ)+f,
(1)
υ=0,
(2)
ρ0wst=ρ0(wsυ)Fs+msΩs.
(3)

Here, υ is the fluid velocity, π is the mechanical pressure (a Lagrange multiplier that ensures that the velocity remains divergence free31), η(w) is the viscosity, ¯=+T is a symmetric gradient, and Σ is the stochastic momentum flux. By denoting the number of species with Nspec, the vector of mass fractions (concentrations) is given by w=(w1,,wNspec), where ws is the mass fraction of species s and sws = 1. We compute the mass density of each species using ρs = ρ0ws, 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 Fs of species s is decomposed into a dissipative flux F¯s and fluctuating flux F̃s,

Fs=F¯s+F̃s,
(4)

and msΩs represents a source term representing stochastic chemistry, where ms 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 sFs = 0 and smsΩs = 0. Based on the fluctuation-dissipation relation, the stochastic momentum flux Σ is modeled as

Σ=ηkBT𝒵mom+(𝒵mom)T,
(5)

where kB is Boltzmann’s constant, and 𝒵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,

μs(x,T,P)=μs0(T,P)+kBTmslog(xsγs),
(6)

where μs0(T,P) is a reference chemical potential and γs(x, T, P) is the activity coefficient (for an ideal mixture, γs = 1). Here x denotes mole fractions, which can be expressed in terms of w as

x=m¯w1m1,,wNspecmNspec,
(7)

where m¯ is the mixture-averaged molecular mass,

m¯=swsms1.
(8)

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.

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

Γx=ρ01ΛW1F¯,
(9)

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

Λss=xsxsDssifss and Λss=ssΛss,
(10)

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

F¯=ρ0WχΓx.
(11)

The stochastic mass fluxes F̃ are given by the fluctuation-dissipation relation

F̃=2m¯ρ0Wχ12𝒵mass,
(12)

where χ12 is a “square root” of χ satisfying χ12(χ12)T=χ, and 𝒵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.

We consider a liquid mixture undergoing Nreact elementary reversible reactions of the form

s=1Nspecνsr+Mss=1NspecνsrMs(r=1,,Nreact),
(13)

where νsr± are the molecule numbers and Ms are chemical symbols. We define the stoichiometric coefficient of species s in the forward reaction r as Δνsr+=νsrνsr+ and the coefficient in the reverse reaction as Δνsr=νsr+νsr. We assume that mass conservation holds in each reaction r; i.e., sΔνsr±ms=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± 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±ΔVdt. Accordingly, the mean number density production rate of species s is given as

Ω¯s=rα=±Δνsrαarα.
(14)

In Sec. II B 1, we give a generalized law of mass action (LMA) based on thermodynamically consistent ar±, and in Sec. II B 2, we present a CME-based stochastic formulation of chemical reactions.

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 as12 

ar±=λrseνsr±μ^s,
(15)

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

ar±=κr±s(xsγs)νsr±,
(16)

where κr±(T,P)=λrsexp(νsr±μ^s0) denotes the forward/reverse reaction rate constant. From the condition ar+=ar at chemical equilibrium, we can express the equilibrium constant as a purely thermodynamic quantity,

Kr(T,P)=κ+κ=expsΔνsr+μ^s0,
(17)

as required by statistical mechanics.

It is important to note that propensity density functions and equilibrium constants are expressed in terms of mole fractions xs (for ideal mixtures) or activities xsγ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 κr±(T,P) cannot be assumed to be constant. On the other hand, in liquid mixtures, where pressure changes are not significant, κr±(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 CME20 for well-mixed33 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 Ns of species s in a given cell during an infinitesimal time interval dt is expressed in terms of the number of occurrences P(ar±ΔVdt) of each reaction r,

dNs=ΩsΔVdt=rα=±ΔνsrαP(arαΔVdt),
(18)

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±ΔVdt) in (18) leads to the chemical Langevin equation (CLE). In this Langevin (Gaussian noise) approximation,4,12

ΩsCLE=rα=±Δνsrαarα+arαZrreact,
(19)

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].

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 A2 molecules undergoing a dimerization reaction

2AA2,
(20)

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 N1 and N2, respectively. By the constant density approximation, N1 and N2 satisfy N1 + 2N2 = ρ0ΔV/mN0, where m is the mass of a monomer and N0 is the total number of A atoms in a cell of volume ΔV. Hence, we denote the equilibrium distribution of the composition with P(N2).

Statistical mechanics predicts that the equilibrium distribution is given by the Einstein distribution, PeS/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

Smix(N1,N2)=kBlog(N1+N2)!N1!N2!,
(21)

and the entropy of the system is

S(N1,N2)=Smix(N1,N2)kB(N1μ^10+N2μ^20).
(22)

Hence, we obtain the equilibrium distribution

P(N2)eS(N02N2,N2)/kB
(23)

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

P(N2+1)P(N2)=(N02N2)(N02N21)(N0N2)(N2+1)exp(2μ^10μ^20).
(24)

We now analyze when detailed balance is achieved for the dimerization reaction (20) with respect to the equilibrium distribution P(N2). The detailed balance condition is given as

P(N2)a+(N2)=P(N2+1)a(N2+1),
(25)

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

a+(N2)κ+N1N1+N2N11N1+N21,
(26a)
a(N2)κN2N1+N2,
(26b)

with N1 = N0 − 2N2. It is important to note that (26) reduces to a+=κ+x12 and a = κx2 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 (N1 − 1)/(N1 + N2 − 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

SStirling=SeqkBNsxslog(xs/xseq),
(27)

where Seq denotes the entropy at xeq and N = sNs. We can further approximate SStirling up to second order in δxs=xsxseq (s = 1, …, Nspec − 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,

Stot=iS(N02N2(i),N2(i)),
(28)

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

Ptot=iP(N2(i)).
(29)

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 N2 ∼ 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¯. At equilibrium, our FHD equations (1)–(3) are linearized for υ = δυ and w=w¯+δw as follows:

t(δυ)=ρ01π+ν2(δυ)+νkBTρ01𝒵mom+(𝒵mom)T,
(30a)
(δυ)=0,
(30b)
t(δw)=D2(δw)+2DkBTρ01μw1𝒵mass+ρ01mΩ1lin,
(30c)

where ν = η/ρ0, D=D12, and μw is the second order derivative of Gibbs free energy with respect to concentration w, given as μw=kBT/[mw¯(1w¯2)] for an ideal mixture. The linearized reaction term is denoted by ρ01mΩ1lin.

We denote the equilibrium structure factors (spectra) by Sυ,υeq(k)=δυ^δυ^* and Sw,weq(k)=δŵδŵ*, 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 

Sυ,υeq=kBTρ0I,Sw,weq=kBTρ0μw=mρ0w¯(1w¯2).
(31)

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

ρ01mΩ1lin=r(δw)+8m2(1w¯)ρ02(1+w¯)κ𝒵react,
(32)

where the linearized reaction rate is

r=4mρ0w¯(1+w¯)2κ.
(33)

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.

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., xs ≪ 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(ρ0/m¯sol)xs, where m¯sol is the mixture-averaged molecular mass among solvent species; see (8). Hence, μ^s can be expressed in terms of ns,

μ^s=μ^s0+logm¯solρ0+logns,
(34)

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

ar±=kr±snsνsr±,
(35)

where kr± 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

tns=Ds2ns+2Dsns𝒵s+Ωs,
(36)

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 model4 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 = {ns} can be described by the differential Chapman–Kolmogorov equation.13 

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 systems4 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.

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 reservoirs27,36 held at fixed concentrations.

The first modification relative to our previous work3 is that for the stochastic mass fluxes F̃, we compute the matrix 2m¯ρ0Wχ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 

ñ(n1,n2)=n1+n22H(n1ΔV)H(n2ΔV),
(37)

where n1 and n2 denote the number densities at the cell centers of two neighboring cells and H is a smoothed Heaviside function, defined as

H(x)=0for x0,xfor 0x1,1for x1.
(38)

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¯Wχ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., niΔ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¯ and F̃ since they depend on and Wχ12, respectively. Unlike χ, however, the matrices and WχW are well defined in the vanishing limit, and we can construct using a special procedure. The basic idea is that we first compute a diffusion sub-matrix χ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+, for all vanishing species s.

To formally describe the procedure for computing 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), i = (1, …, 6). As graphically illustrated using the 6-species system in Fig. 1, there are four cases to consider when one populates ()ij,

(Wχ)ij=(39a)wiχm(i)m(j)sub,                m(i)0,m(j)0(yellow),(39b)miDim¯,                  m(i)=0,  j=i(red),(39c)0,                        m(i)=0,  ji,(39d)wiDjkm(k)0xkDkjχm(i)m(k)submjm¯,          m(i)0,m(j)=0(blue),
where

Dj=km(k)0xkDkj1.
(40)

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χ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.

FIG. 1.

Graphical depiction of the expansion of sub-matrix χsub into a full matrix for a 6-species system having vanishing species w2 and w4. Depicted is the full matrix , where the colors correspond to the cases in (39).

FIG. 1.

Graphical depiction of the expansion of sub-matrix χsub into a full matrix for a 6-species system having vanishing species w2 and w4. Depicted is the full matrix , where the colors correspond to the cases in (39).

Close modal

We note that for each vanishing species only the diagonal element of remains nonzero. Hence, the diffusion of a dilute species s (ws ≪ 1) becomes decoupled from other species [see (A8)] and the effective diffusion coefficient Ds 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¯s=F̃s=0 for vanishing species and ensures the mass conservation condition over all species, sF¯s=sF̃s=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.

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,

ρ0dws,idt=ρ0(wsυ)Fsi+msrα=±ΔνsrαP(arαΔVdt)ΔVdti,
(41)

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 sampling29,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 tn = nΔt to time tn+1 = (n + 1)Δt in four steps:

  1. Perform a predictor Stokes solve for the velocity υn+1,* at tn+1.

  2. Calculate predictor mass densities ρsn+12,* at the midpoint time t=tn+12Δt.

  3. Calculate corrector mass densities ρsn+1 at time tn+1.

  4. Perform a corrector Stokes solve for velocity υn+1 at tn+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., fn = f(wn). Also, 𝒲momn and 𝒲(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 𝒲¯mom𝒲mom+𝒲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).

Algorithm 1.

Advancing the mass densities ρsn=ρ0wsn and velocity υn from time tn to tn+1 = tn + Δt. We list both centered and BDS options for advection; the BDS notation is defined in Sec. III B in Ref. 26.

Advancing the mass densities ρsn=ρ0wsn and velocity υn from time tn to tn+1 = tn + Δt. We list both centered and BDS options for advection; the BDS notation is defined in Sec. III B in Ref. 26.
Advancing the mass densities ρsn=ρ0wsn and velocity υn from time tn to tn+1 = tn + Δt. We list both centered and BDS options for advection; the BDS notation is defined in Sec. III B in Ref. 26.

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 systems4 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.

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δvy and δw (see Appendix C in Ref. 12) take the form

t(δv^)=νk2(δv^)+2νkBTρ01ik𝒵^mom+gζ(δŵ),
(50a)
t(δŵ)=h(δv^)Dk2(δŵ)+2DkBTρ01μw1ik𝒵^massr(δŵ)+2rkBTρ01μw1Z^react.
(50b)

Here kk 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 = hey. 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υ = 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 kBT/ρ0 = 1. To denote how fast momentum diffusion, species diffusion, and reaction are, we define the following dimensionless Courant numbers:

α=νΔtΔx2,β=DΔtΔx2,γ=rΔt.
(51)

To consider the case of a relatively large Δt with a large Schmidt number Sc = 103 (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 ξ=D/r=10Δx and ξ=1210Δx, respectively. Other parameters are chosen so that μwh2Δt2 = 100, and gζhΔt2 = 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,

Sw,wneq=kBTρ0D(D+ν)k4h2.
(52)

However, the power law is suppressed at small k by gravity27 or reaction.12 

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,wk−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).

FIG. 2.

Structure factors for giant concentration fluctuations. Panel (a) shows Sw,w for the non-reactive case with gravity. Results from our numerical scheme (inertial midpoint) and two earlier schemes are compared with the exact result. Panel (b) shows Sw,wneq=Sw,wSw,weq for the reactive case with no gravity. For two rate constants, results from our numerical scheme are compared with the exact result. The x-axis is the dimensionless wavenumber K = kΔx.

FIG. 2.

Structure factors for giant concentration fluctuations. Panel (a) shows Sw,w for the non-reactive case with gravity. Results from our numerical scheme (inertial midpoint) and two earlier schemes are compared with the exact result. Panel (b) shows Sw,wneq=Sw,wSw,weq for the reactive case with no gravity. For two rate constants, results from our numerical scheme are compared with the exact result. The x-axis is the dimensionless wavenumber K = kΔx.

Close modal

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],

Sw,wneq=kBTρ0(Dk2+r)[(D+ν)k2+r]h2.
(53)

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

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.

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

sucrose+H2Oglucose+fructose.
(54)

Sucrose is particularly dilute, with only ∼10 molecules per computational cell, whereas there are ∼107 glucose and fructose molecules and ∼1010 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, Ds (s = 1, 2, 3),41,42 and the self-diffusion coefficient of water Dwater,43 the Maxwell–Stefan binary diffusion coefficients are assigned as in Ref. 44,

Đs4=Ds,Đss=DsDsDwater(s,s=1,2,3).
(55)

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+n3, a reverse rate a = kn1n2, 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ΔV=10 and n1eq=n2eq are subsequently determined from equilibrium. The resulting equilibrium mass fractions are w1eq=w2eq=4.9×1003, w3eq=5.7×1009, 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 Ds,maxΔtx2 = 0.07 for species diffusion and νΔtx2 = 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 tT/10.

We recall that the number of sucrose molecules in a cell N = n3ΔV has a continuous range in FHD simulations. We define its discrete distribution as P(N)=N12N+12ρ(N)dN, 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 PPoisson(N), and ρ(N) is essentially zero for negative values of N. We note that PPoisson(N) is significantly different from its Gaussian approximation PGauss(N). For the larger Δt, while the remarkable agreement with the Poisson distribution is still observed, negative values of N start to appear, yielding 0ρ(N)dN3×105; 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 theory3), indicating that there are no spurious correlations between cells.

FIG. 3.

Equilibrium distribution for a dilute sugar solution undergoing a hydrolysis reaction. Numerical results for the distribution (histogram) P(N) of the number of sucrose molecules in a cell N are compared with the physically correct Poisson distribution PPoisson(N) and its Gaussian approximation PGauss(N). In the inset, numerical results for the continuous distribution ρ(N) are shown near N = 0. Results from a smaller time step size Δt = 10−5 are plotted with error bars corresponding to two standard deviations, whereas those from a larger time step size Δt = 10−4 are plotted without errorbars for clarity.

FIG. 3.

Equilibrium distribution for a dilute sugar solution undergoing a hydrolysis reaction. Numerical results for the distribution (histogram) P(N) of the number of sucrose molecules in a cell N are compared with the physically correct Poisson distribution PPoisson(N) and its Gaussian approximation PGauss(N). In the inset, numerical results for the continuous distribution ρ(N) are shown near N = 0. Results from a smaller time step size Δt = 10−5 are plotted with error bars corresponding to two standard deviations, whereas those from a larger time step size Δt = 10−4 are plotted without errorbars for clarity.

Close modal

We next consider a liquid mixture undergoing the dimerization reaction (20). This binary system contains monomers A (s = 1) and dimers A2 (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 (N1 and N2) do not vary independently due to the constant density assumption N1 + 2N2 = ρ0ΔV/m, where m is the mass of a monomer. Therefore, we investigate the equilibrium distribution of N2, P(N2), for N1 + 2N2 = 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, D12=1, and m = 1 with kB = 1. We set ρ0 = 40 with w1eq=w2eq=0.5 so that N1 + 2N2 = 40 with N1eq=20 and N2eq=10. We set the reaction rates a± as in (26) with a modification

a+=κ+N1+N1++N2+(N11)+N1++N2+1,
(56)

where N+ = max(N, 0). Note that (56) turns off unphysical reactions when 0 < N1 < 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 N2=N2N2P(N2)=N2eq. The magnitude of κ± is determined so that the linearized reaction rate r = 0.1 [see (33)] gives a penetration depth ξD12/r=10Δx. We set η = 103 and T=103. We use a small Δt = 10−2 to minimize temporal integration errors. For 16 independent samples with 105 time steps, we collect data every 102 time steps, discarding the first 104 time steps.

In Fig. 4, we compare the simulation result for the equilibrium distribution P(N2) with theoretical results obtained in Sec. III A. We denote the exact Einstein distribution obtained from the entropy expression (22) by Pexact, Stirling’s approximation result obtained from (27) by PStirling ∼ exp(SStirling/kB), and the Gaussian approximation of PStirling by PGauss. Note that Pexact is a discrete distribution, whereas PStirling and PGauss have continuous ranges, 0 < N2 < 20 and − < N2 < , respectively. Significant deviations of PGauss from Pexact 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 Pexact. Even beyond this range, it gives sensible values with accuracy better than or comparable to PStirling. Measurable deviations are observed only for larger values N2 = 19 and 20, for which the occupation probabilities are already very small (Pexact(N2) < 10−6).

FIG. 4.

Equilibrium distribution for a binary ideal mixture undergoing a dimerization reaction. The distribution P(N2) of the number of dimers in each cell N2 is computed using our numerical method and compared with various theoretical results (see text). Error bars correspond to two standard deviations.

FIG. 4.

Equilibrium distribution for a binary ideal mixture undergoing a dimerization reaction. The distribution P(N2) of the number of dimers in each cell N2 is computed using our numerical method and compared with various theoretical results (see text). Error bars correspond to two standard deviations.

Close modal

It is important to note that statistically identical results for P(N2) 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 limit4 and analyzed only for this case, the validity of its extension to non-dilute solutions cannot be taken for granted.

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,w1Sw1,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 = 103, 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 Lx = Ly = 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 = 106 so that there are sufficiently many monomers and dimers in a cell, N1 + 2N2 = 2.5 × 105 for ρ0 = 1 and m = 1. We set D12=1, η = 103, and kBT = 103. 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 ξ=10Δ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 D12Δt/Δx2=0.1 and νΔtx2 = 100. We ran 105 steps discarding the first 104 steps and computed the steady-state monomer concentration profile w¯1(y) and Sw1,w1neq(kx).

To impose a concentration gradient in the y-direction, no-slip reservoir conditions27 are imposed with (w1, w2) = (0.49, 0.51) at y = 0 and (w1, w2) = (0.51, 0.49) at y = Ly. 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¯1(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¯1(y) becomes more evident. This is because a larger region around y = Ly/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).

FIG. 5.

Giant fluctuations with a dimerization reaction. (a) Steady-state monomer concentration profile w¯1(y). Results from two linearized reaction rates r and the non-reactive case are compared. (b) Nonequilibrium enhancement Sw1,w1neq in the structure factor of the monomer mass fraction. Numerical results (depicted by symbols) are compared with theoretical predictions obtained under a linear gradient setting (depicted by a solid line for r = 0.4, dashed line for r = 0.04, and dotted line for r = 0). A slope marker for the k−4 decay is drawn and arrows denoting kx=10r/D12 are depicted for r = 0.4 (solid) and r = 0.04 (dashed). Note that nonlinear gradients develop in reactive cases, which explains discrepancies between numerical and theoretical results at small to intermediate wavenumbers.

FIG. 5.

Giant fluctuations with a dimerization reaction. (a) Steady-state monomer concentration profile w¯1(y). Results from two linearized reaction rates r and the non-reactive case are compared. (b) Nonequilibrium enhancement Sw1,w1neq in the structure factor of the monomer mass fraction. Numerical results (depicted by symbols) are compared with theoretical predictions obtained under a linear gradient setting (depicted by a solid line for r = 0.4, dashed line for r = 0.04, and dotted line for r = 0). A slope marker for the k−4 decay is drawn and arrows denoting kx=10r/D12 are depicted for r = 0.4 (solid) and r = 0.04 (dashed). Note that nonlinear gradients develop in reactive cases, which explains discrepancies between numerical and theoretical results at small to intermediate wavenumbers.

Close modal

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 wavenumber36 

k̃x=sin(kxΔx/2)Δx/2
(57)

is used instead of kx. For the non-reactive case (r = 0), a clear k−4 power law is observed until the confinement effect becomes significant for small kxLy1. For the reactive cases, the k−4 power law is only observed at large kxr/D12. For larger r, the k−4 power law appears in a narrower range of kx 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 factor47 

1+41cosh(kxLy)kxLykxLy+sinh(kxLy),
(58)

the theoretical prediction is further improved at small kx, 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 h2 is not obvious. Considering that Sw1,w1neq is an averaged structure factor for different values of y, we estimate the value of h2 from the profile of w¯1(y) using a spatial average,

h2=1Ly0Lydw¯1dy2dy.
(59)

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 kx 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.

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

HCl+NaOHNaCl+H2O
(60)

plays a negligible role in this problem.48 We consider a Hele-Shaw cell with side lengths Lx = Ly = 1.6 and Lz = 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

f(w)=ρ0s=13αsMswsgey,
(61)

where αs is the solutal expansion coefficient and Ms 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+ = kn1n2. However, we note that neutralization is a diffusion-limited reaction. In other words, reaction occurs extremely fast (with rate λ ∼ 1011 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 cm3 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],

δw=1ρ0ΔV(Iw01T)(W0M)12zmass,
(62)

where W0=diag(ws0), M = diag(ms), and zmass is a vector with i.i.d. standard normal random variables. The initial momentum fluctuations are generated in a similar manner,

δυ=kBTρ0ΔVzmom,
(63)

where zmom is a vector with i.i.d. standard normal random variables, followed by an L2 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.

TABLE I.

Four simulations performed for the buoyancy-driven instability with an acid-base neutralization reaction.

Initial fluctuationsStochastic fluxes
ChemistryMassMomentumMassMomentum
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 fluctuationsStochastic fluxes
ChemistryMassMomentumMassMomentum
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

FIG. 6.

Asymmetric growth of convective chemo-hydrodynamic fingering patterns in a Hele-Shaw cell, induced by a double-diffusive instability in the presence of a neutralization reaction. The left and middle columns, (a) and (b), depict the mass density profiles of chemical species at t = 25 and t = 40 (simulation A), whereas the right column (c) displays the non-reactive case at t = 40 (simulation B). The top, middle, and bottom rows show ρHCl, ρNaOH, and ρNaCl, respectively. For each species, two-dimensional slices of the three-dimensional field ρs(x, y, z) are shown. The square images show ρs(x, y, z = Lz/2) (halfway between the glass plates), and the thin vertical side bars show the slice ρs(x = 0, y, z) corresponding to the left edge of the square images. Both simulations were initiated with natural mass and momentum fluctuations without any artificial perturbation.

FIG. 6.

Asymmetric growth of convective chemo-hydrodynamic fingering patterns in a Hele-Shaw cell, induced by a double-diffusive instability in the presence of a neutralization reaction. The left and middle columns, (a) and (b), depict the mass density profiles of chemical species at t = 25 and t = 40 (simulation A), whereas the right column (c) displays the non-reactive case at t = 40 (simulation B). The top, middle, and bottom rows show ρHCl, ρNaOH, and ρNaCl, respectively. For each species, two-dimensional slices of the three-dimensional field ρs(x, y, z) are shown. The square images show ρs(x, y, z = Lz/2) (halfway between the glass plates), and the thin vertical side bars show the slice ρs(x = 0, y, z) corresponding to the left edge of the square images. Both simulations were initiated with natural mass and momentum fluctuations without any artificial perturbation.

Close modal

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.

FIG. 7.

Influence of different types of thermal fluctuations on the formation of fingering patterns. We compare the mass density profiles ρs(x, y, z = Lz/2) of NaCl at t = 30 halfway between the glass plates for three simulations. Simulation A (left) corresponds to the full fluctuating hydrodynamics equations. Compared with simulation A, all stochastic mass components (stochastic mass flux and stochastic chemistry) are omitted in simulation C (middle). Simulation D (right) is similar to simulation C but with initial velocity fluctuations also removed. Red circles indicate the areas with the biggest differences.

FIG. 7.

Influence of different types of thermal fluctuations on the formation of fingering patterns. We compare the mass density profiles ρs(x, y, z = Lz/2) of NaCl at t = 30 halfway between the glass plates for three simulations. Simulation A (left) corresponds to the full fluctuating hydrodynamics equations. Compared with simulation A, all stochastic mass components (stochastic mass flux and stochastic chemistry) are omitted in simulation C (middle). Simulation D (right) is similar to simulation C but with initial velocity fluctuations also removed. Red circles indicate the areas with the biggest differences.

Close modal

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.

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 micropumps50 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

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.

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

        wN0+,
(A1a)
wiwi0>0(i=1,,N1) with i=1N1wi0=1.
(A1b)

We show next that each component of converges to

(Wχ)ijwi0χijsub  (i,j=1,,N1),
(A2a)
(Wχ)NNmNDNm¯0,
(A2b)
(Wχ)Ni0    (i=1,,N1),
(A2c)
(Wχ)iNwi0DNk=1N1χiksubxk0ĐkNmNm¯0(i=1,,N1),
(A2d)

where χsub is the diffusion matrix of the subsystem consisting of non-vanishing species with w0=(w10,,wN10), and x0 and m¯0 are computed from (7) and (8) using w0. We also show that the trace diffusion coefficient DN of species N in the fluid mixture with composition w0 is expressed as

DN=k=1N1xk0ĐkN1.
(A3)

By rearranging the definition of χ,3 

χ=(Λ+αwwT)11α11T,
(A4)

where α ≠ 0, and using 1TΛ = 0, χw = 0, and 1Tw = 1, we obtain

χΛ+1wT=I.
(A5)

Looking at the (N, N)-component of (A5) and using (10), we have

xNχNNi=1N1xiĐiN1χNiχNN+wN=1.
(A6)

Noting that χNN=O(wN1) and

χNiχNN0,
(A7)

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

F¯N=ρ0mNDNm¯0xN.
(A8)

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

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

ws0=aUws0,lower+(1aU)ws0,upper,
(B1)

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 ρ0w3(r, t)dr, as well as the L1-norm of the vy field ∥vy∥ = |vy(r, t)|dr.

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−22k ≲ 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Δ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).

FIG. 8.

Effects of the reaction rate constant k on the fingering instability observed when a layer of NaOH is placed on top of HCl solution, for deterministic simulations with a randomly perturbed initial interface. Panel (a) shows the time profiles of the norm of the vy field for various values of k. Panel (b) shows the gross NaCl production up to time t as a function of k. Solid lines denote the results for Δt = 10−3, whereas dotted lines in the same colors depict the results for Δt = 10−2. Arrows indicate k=4/(nHCl0Δt) for Δt = 10−2 (dotted line) and Δt = 10−3 (solid line), where nHCl0 is the initial number density of HCl in the lower layer.

FIG. 8.

Effects of the reaction rate constant k on the fingering instability observed when a layer of NaOH is placed on top of HCl solution, for deterministic simulations with a randomly perturbed initial interface. Panel (a) shows the time profiles of the norm of the vy field for various values of k. Panel (b) shows the gross NaCl production up to time t as a function of k. Solid lines denote the results for Δt = 10−3, whereas dotted lines in the same colors depict the results for Δt = 10−2. Arrows indicate k=4/(nHCl0Δt) for Δt = 10−2 (dotted line) and Δt = 10−3 (solid line), where nHCl0 is the initial number density of HCl in the lower layer.

Close modal

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 νΔtx2 = 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.

1.
A.
Vailati
and
M.
Giglio
, “
Giant fluctuations in a free diffusion process
,”
Nature
390
,
262
265
(
1997
).
2.
L.
Lemaigre
,
M. A.
Budroni
,
L. A.
Riolfo
,
P.
Grosfils
, and
A.
De Wit
, “
Asymmetric Rayleigh–Taylor and double-diffusive fingers in reactive systems
,”
Phys. Fluids
25
,
014103
(
2013
).
3.
A.
Donev
,
A.
Nonaka
,
A. K.
Bhattacharjee
,
A. L.
Garcia
, and
J. B.
Bell
, “
Low Mach number fluctuating hydrodynamics of multispecies liquid mixtures
,”
Phys. Fluids
27
,
037103
(
2015
).
4.
C.
Kim
,
A.
Nonaka
,
J. B.
Bell
,
A. L.
Garcia
, and
A.
Donev
, “
Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach
,”
J. Chem. Phys.
146
,
124110
(
2017
).
5.
A.
Donev
,
C.-Y.
Yang
, and
C.
Kim
, “
Efficient reactive Brownian dynamics
,”
J. Chem. Phys.
148
,
034103
(
2018
).
6.
N. G.
Van Kampen
,
Stochastic Processes in Physics and Chemistry
(
North-Holland
,
Amsterdam
,
1983
).
7.
L. D.
Landau
and
E. M.
Lifshitz
, in
Fluid Mechanics
, Volume 6 of Course of Theoretical Physics, 1st English ed. (
Pergamon Press
,
1959
), Chap. XVII, pp.
523
529
.
8.
J. M.
Ortiz de Zárate
and
J. V.
Sengers
,
Hydrodynamic Fluctuations in Fluids and Fluid Mixtures
(
Elsevier
,
2006
).
9.
G. A.
Bird
,
Molecular Gas Dynamics and the Direct Simulation of Gas Flows
, 2nd ed. (
Clarendon Press
,
1994
).
10.
M.
Malek Mansour
and
F.
Baras
, “
Microscopic simulation of chemical systems
,”
Physica A
188
,
253
276
(
1992
).
11.
F.
Baras
and
M.
Malek Mansour
, “
Microscopic simulation of chemical instabilities
,”
Adv. Chem. Phys.
100
,
393
474
(
1997
).
12.
A. K.
Bhattacharjee
,
K.
Balakrishnan
,
A. L.
Garcia
,
J. B.
Bell
, and
A.
Donev
, “
Fluctuating hydrodynamics of multi-species reactive mixtures
,”
J. Chem. Phys.
142
,
224107
(
2015
).
13.
C. W.
Gardiner
,
Handbook of Stochastic Methods
, 2nd ed. (
Springer-Verlag
,
1985
).
14.
F.
Baras
and
M.
Malek Mansour
, “
Reaction-diffusion master equation: A comparison with microscopic simulations
,”
Phys. Rev. E
54
,
6139
6148
(
1996
).
15.
G.
Nicolis
and
I.
Prigogine
,
Self-Organization in Nonequilibrium Systems
(
Wiley-Interscience
,
1977
).
16.
P.
Dziekan
,
L.
Signon
,
B.
Nowakowski
, and
A.
Lemarchand
, “
Reaction-diffusion approach to prevertebrae formation: Effect of a local source of morphogen
,”
J. Chem. Phys.
139
,
114107
(
2013
).
17.
I.
Pagonabarraga
,
A.
Pérez-Madrid
, and
J. M.
Rubí
, “
Fluctuating hydrodynamics approach to chemical reactions
,”
Physica A
237
,
205
219
(
1997
).
18.
D.
Bedeaux
,
I.
Pagonabarraga
,
J. M.
Ortiz de Zárate
,
J. V.
Sengers
, and
S.
Kjelstrup
, “
Mesoscopic non-equilibrium thermodynamics of non-isothermal reaction-diffusion
,”
Phys. Chem. Chem. Phys.
12
,
12780
12793
(
2010
).
19.
J.
Keizer
,
Statistical Thermodynamics of Nonequilibrium Processes
(
Springer-Verlag
,
New York
,
1987
).
20.
D. T.
Gillespie
,
A.
Hellander
, and
L. R.
Petzold
, “
Perspective: Stochastic algorithms for chemical kinetics
,”
J. Chem. Phys.
138
,
170901
(
2013
).
21.
R.
Erban
,
J.
Chapman
, and
P.
Maini
, “
A practical guide to stochastic simulations of reaction-diffusion processes
,” e-print arXiv: 0704.1908 (
2007
).
22.
J. M.
Ortiz de Zárate
,
J. V.
Sengers
,
D.
Bedeaux
, and
S.
Kjelstrup
, “
Concentration fluctuations in nonisothermal reaction-diffusion systems
,”
J. Chem. Phys.
127
,
034501
(
2007
).
23.
D.
Bedeaux
,
J. M.
Ortiz de Zárate
,
I.
Pagonabarraga
,
J. V.
Sengers
, and
S.
Kjelstrup
, “
Concentration fluctuations in non-isothermal reaction-diffusion systems. II. The nonlinear case
,”
J. Chem. Phys.
135
,
124516
(
2011
).
24.
A.
Donev
,
E.
Vanden-Eijnden
,
A.
Garcia
, and
J.
Bell
, “
On the accuracy of finite-volume schemes for fluctuating hydrodynamics
,”
Commun. Appl. Math. Comput. Sci.
5
,
149
197
(
2010
).
25.
K.
Balakrishnan
,
A. L.
Garcia
,
A.
Donev
, and
J. B.
Bell
, “
Fluctuating hydrodynamics of multispecies nonreactive mixtures
,”
Phys. Rev. E
89
,
013017
(
2014
).
26.
A.
Nonaka
,
Y.
Sun
,
J. B.
Bell
, and
A.
Donev
, “
Low Mach number fluctuating hydrodynamics of binary liquid mixtures
,”
Commun. Appl. Math. Comput. Sci.
10
,
163
204
(
2015
).
27.
A.
Donev
,
A.
Nonaka
,
Y.
Sun
,
T. G.
Fai
,
A. L.
Garcia
, and
J. B.
Bell
, “
Low Mach number fluctuating hydrodynamics of diffusively mixing fluids
,”
Commun. Appl. Math. Comput. Sci.
9
,
47
105
(
2014
).
28.
M.
Cai
,
A.
Nonaka
,
J. B.
Bell
,
B. E.
Griffith
, and
A.
Donev
, “
Efficient variable-coefficient finite-volume Stokes solvers
,”
Commun. Comput. Phys.
16
,
1263
1297
(
2014
).
29.
Y.
Hu
,
T.
Li
, and
B.
Min
, “
A weak second order tau-leaping method for chemical kinetic systems
,”
J. Chem. Phys.
135
,
024113
(
2011
).
30.
D. T.
Gillespie
, “
Approximate accelerated stochastic simulation of chemically reacting systems
,”
J. Chem. Phys.
115
,
1716
1733
(
2001
).
31.
A. J.
Chorin
and
J. E.
Marsden
,
A Mathematical Introduction to Fluid Mechanics
(
Springer
,
1990
), Vol. 3.
32.
M. R.
Marcelin
, “
Sur la mécanique des phénomènes irréversibles
,”
Comptes rendus hebdomadaires des séances de l’Académie des sciences, tome
151
,
1052
1055
(
1910
), http://visualiseur.bnf.fr/CadresFenetre?O=NUMM-3104&M=tdm.
33.

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 ξ=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 

34.
D.
Kelly
and
E.
Vanden-Eijnden
, “
Fluctuations in the heterogeneous multiscale methods for fast-slow systems
,”
Res. Math. Sci.
4
,
23
(
2017
).
35.
S.
Delong
,
Y.
Sun
,
B. E.
Griffith
,
E.
Vanden-Eijnden
, and
A.
Donev
, “
Multiscale temporal integrators for fluctuating hydrodynamics
,”
Phys. Rev. E
90
,
063312
(
2014
).
36.
F. B.
Usabiaga
,
J. B.
Bell
,
R.
Delgado-Buscalioni
,
A.
Donev
,
T. G.
Fai
,
B. E.
Griffith
, and
C. S.
Peskin
, “
Staggered schemes for fluctuating hydrodynamics
,”
Multiscale Model. Simul.
10
,
1369
1408
(
2012
).
37.
J. B.
Bell
,
C. N.
Dawson
, and
G. R.
Shubin
, “
An unsplit, higher order Godunov method for scalar conservation laws in multiple dimensions
,”
J. Comput. Phys.
74
,
1
24
(
1988
).
38.
A.
Nonaka
,
S.
May
,
A. S.
Almgren
, and
J. B.
Bell
, “
A three-dimensional, unsplit Godunov method for scalar conservation law
,”
SIAM J. Sci. Comput.
33
,
2039
2062
(
2011
).
39.
D. F.
Anderson
and
J. C.
Mattingly
, “
A weak trapezoidal method for a class of stochastic differential equations
,”
Commun. Math. Sci.
9
,
301
318
(
2011
).
40.
S.
Delong
,
B. E.
Griffith
,
E.
Vanden-Eijnden
, and
A.
Donev
, “
Temporal integrators for fluctuating hydrodynamics
,”
Phys. Rev. E
87
,
033302
(
2013
).
41.
A.
Venâncio
and
J. A.
Teixeira
, “
Characterization of sugar diffusion coefficients in alginate membranes
,”
Biotechnol. Tech.
11
,
183
185
(
1997
).
42.
J. F.
Tilley
and
R.
Mills
, “
Relations between mutual and intradiffusion coefficients in aqueous sucrose solutions
,”
J. Phys. Chem.
71
,
2756
2757
(
1967
).
43.
K.
Krynicki
,
C. D.
Green
, and
D. W.
Sawyer
, “
Pressure and temperature dependence of self-diffusion in water
,”
Faraday Discuss. Chem. Soc.
66
,
199
208
(
1978
).
44.
X.
Liu
,
A.
Bardow
, and
T. J. H.
Vlugt
, “
Multicomponent Maxwell–Stefan diffusivities at infinite dilution
,”
Ind. Eng. Chem. Res.
50
,
4776
4782
(
2011
).
45.
R. N.
Goldberg
,
Y. B.
Tewari
, and
J. C.
Ahluwalia
, “
Thermodynamics of the hydrolysis of sucrose
,”
J. Biol. Chem.
264
,
9901
9904
(
1989
).
46.
R.
Wolfenden
and
Y.
Yuan
, “
Rates of spontaneous cleavage of glucose, fructose, sucrose, and trehalose in water, and the catalytic proficiencies of invertase and trehalas
,”
J. Am. Chem. Soc.
130
,
7548
7549
(
2008
).
47.
J. M.
Ortiz de Zárate
,
T. R.
Kirkpatrick
, and
J. V.
Sengers
, “
Non-equilibrium concentration fluctuations in binary liquids with realistic boundary conditions
,”
Eur. Phys. J. E
38
,
99
(
2015
).
48.
C.
Almarcha
,
P. M. J.
Trevelyan
,
P.
Grosfils
, and
A.
De Wit
, “
Thermal effects on the diffusive layer convection instability of an exothermic acid-base reaction front
,”
Phys. Rev. E
88
,
033009
(
2013
).
49.
R.
Erban
and
S. J.
Chapman
, “
Stochastic modeling of reaction-diffusion processes: Algorithms for bimolecular reactions
,”
Phys. Biol.
6
,
046001
(
2009
).
50.
M. J.
Esplandiu
,
A. A.
Farniya
, and
D.
Reguera
, “
Key parameters controlling the performance of catalytic motors
,”
J. Chem. Phys.
144
,
124702
(
2016
).
51.
M. B.
Andersen
,
K. M.
Wang
,
J.
Schiffbauer
, and
A.
Mani
, “
Confinement effects on electroconvective instability
,”
Electrophoresis
38
,
702
711
(
2017
).
52.
J.-P.
Péraud
,
A.
Nonaka
,
A.
Chaudhri
,
J. B.
Bell
,
A.
Donev
, and
A. L.
Garica
, “
Low Mach number fluctuating hydrodynamics for electrolytes
,”
Phys. Rev. Fluids
1
,
074103
(
2016
).
53.
B. E.
Griffith
and
C. S.
Peskin
, “
Electrophysiology
,”
Commun. Pure Appl. Math.
66
,
1837
1913
(
2013
).
54.
J. L.
Moran
and
J. D.
Posner
, “
Electrokinetic locomotion due to reaction-induced charge auto-electrophoresis
,”
J. Fluid Mech.
680
,
31
66
(
2011
).
55.
J. L.
Moran
and
J. D.
Posner
, “
Phoretic self-propulsion
,”
Annu. Rev. Fluid Mech.
49
,
511
540
(
2017
).