Fluctuating hydrodynamics of reactive liquid mixtures

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.¹ Buoyancy-driven instabilities can be triggered or affected by thermal fluctuations.^2,3 In reaction-diffusion systems, thermal fluctuations can accelerate the formation of Turing patterns on a macroscopic time scale⁴ and induce long-time memory in the chemical kinetics of a diffusion-limited system.⁵ 

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)⁶ 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.¹² While our previous work on reaction-diffusion systems⁴ also employed the CME, it was restricted to dilute solutions. Here we generalize the CME to non-dilute ideal mixtures with a complete Maxwell–Stefan formulation of diffusive transport in multispecies mixtures. This includes cross-diffusion coupling among distinct species and can account for deviations from ideality, unlike the standard reaction-diffusion master equation (RDME) approach.^13–15 Finally, by including the fluctuating Navier–Stokes equations in the model, we account for advection by thermal velocity fluctuations, which is necessary to capture giant nonequilibrium composition fluctuations.^10,11,16

Our approach is related to, but also distinct from, prior work on fluctuating hydrodynamics for reactive liquid mixtures. An alternative Langevin-based approach proposed in Ref. 17, and extended to full hydrodynamics in Ref. 18, represents reactions as a diffusion process along an internal reaction coordinate, driven by the Gaussian noise. This description is fully consistent with nonequilibrium thermodynamics and fluctuating hydrodynamics, but is not easily extensible to multispecies mixtures, and, importantly, is expensive to use in numerical simulations because it requires introducing an additional reaction coordinate, thus effectively increasing the dimensionality of the problem. Instead, in our approach, we only consider the reactant and product states and consider reactions as a jump process between these two states, driven by the Poisson noise. The deterministic (macroscopic) as well as a linearized version of the FHD equations we consider here is the same as that obtained from a quasi-stationary approximation of the model developed in Ref. 18 [see Eq. (26) in Ref. 18 and the book by Keizer¹⁹), as we have discussed in more detail in prior work.¹² 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.⁴ 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,³ in which the density is essentially independent of composition at a fixed pressure and temperature. The resulting Boussinesq (incompressible) approximation of the momentum equation enables us to construct an efficient numerical method that accounts for inertial effects important in buoyancy-driven fluid flows yet remains robust for small Reynolds numbers and large Schmidt numbers. The spatio-temporal discretization of the FHD equations is based on our previous work³ but with some important improvements necessary for simulating complex reactive mixtures at small length scales. Notably, we extend our previous work on reaction-diffusion systems⁴ 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.²⁴ Using this framework, we have previously developed stochastic simulation methods for gas mixtures²⁵ 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.²⁷ To avoid severe restriction on the time step size when the Schmidt number is large, we have developed an implicit temporal discretization of viscous dissipation²⁶ that relies on a variable-coefficient multigrid precondition to solve the coupled velocity-pressure Stokes system.²⁸ 

In this paper, we make three novel contributions to the numerical methodology developed in our prior work. First, by incorporating a second-order midpoint tau-leaping scheme²⁹ into our prior algorithms for multispecies miscible liquid mixtures,³ 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,³⁰ it is robust for large composition fluctuations, while also being efficient for weak fluctuations. Second, the midpoint scheme is constructed to be robust for large Schmidt numbers, i.e., much faster momentum diffusion compared to mass diffusion, as is typical in liquid systems. In particular, the numerical method reproduces the correct spectrum of giant nonequilibrium fluctuations even for time step sizes much larger than the stability limit dictated by fast momentum diffusion, while also preserving the slow inertial momentum dynamics at large scales. Third, we take careful attention to handling vanishing species robustly both in the formulation of the multispecies diffusion model and in the numerical algorithm.

The rest of the paper is organized as follows. In Sec. II, we present the formulation of the FHD equations coupled with the CME formulation of reactions. In Sec. III, we present a numerical scheme that can solve these equations accurately and robustly even in the presence of large composition fluctuations and vanishing species. In Sec. IV, we present numerical results for four examples and discuss various aspects of our numerical method and the effects of thermal fluctuations. First, we verify that for dilute solutions our algorithm preserves the robustness and accuracy properties of our previous method for reaction-diffusion systems⁴ by modeling the hydrolysis of sucrose at micrometer scales. Second, to assess the fidelity of our approach in a non-dilute setting, we consider a binary mixture undergoing a dimerization reaction 2A ⇌ A₂ at thermodynamic equilibrium with a small number of molecules per cell. Third, we also study such a mixture out of equilibrium in the presence of giant nonequilibrium fluctuations with a large number of molecules per cell. Fourth, we use our numerical algorithm to simulate a diffusively driven gravitational instability in the presence of an acid-base neutralization reaction recently studied experimentally² 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 ρ₀,

Here, $υ$ is the fluid velocity, π is the mechanical pressure (a Lagrange multiplier that ensures that the velocity remains divergence free³¹), η(⁠$w$⁠) is the viscosity, $∇¯=∇+∇T$ is a symmetric gradient, and Σ is the stochastic momentum flux. By denoting the number of species with N_spec, the vector of mass fractions (concentrations) is given by $w=(w1,…,wNspec)$⁠, where $ws$ is the mass fraction of species s and ∑_s$ws$ = 1. We compute the mass density of each species using ρ_s = ρ₀$ws$⁠, and thus the total mass density ∑_sρ_s = ρ₀ is strictly constant. The buoyancy force f(⁠$w$⁠) is a problem-specific function of $w$⁠. The total diffusive mass flux F_s of species s is decomposed into a dissipative flux $F¯s$ and fluctuating flux $F̃s$⁠,

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

where k_B 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,

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

where $m¯$ is the mixture-averaged molecular mass,

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

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

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

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

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

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 N_react elementary reversible reactions of the form

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

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.

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

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

where $κr±(T,P)=λr∏s⁡exp(ν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,

as required by statistical mechanics.

It is important to note that propensity density functions and equilibrium constants are expressed in terms of mole fractions x_s (for ideal mixtures) or activities x_sγ_s. This generalized LMA has a different form compared to the number density based LMA used for ideal gas mixtures in our prior work.¹² 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.

We believe that an accurate mesoscopic chemistry description should be based on a master equation approach, which leads to the CME²⁰ for well-mixed³³ 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,⁴ which is related to the computationally efficient tau leaping method.³⁰ The change in the number of molecules N_s of species s in a given cell during an infinitesimal time interval dt is expressed in terms of the number of occurrences $P(ar±ΔVdt)$ of each reaction r,

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

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

where $Zrreact(r,t)$ denotes a standard GWN field. The Langevin description is justified in the limit of small Gaussian fluctuations with respect to average concentrations.²⁰ However, the Langevin description predicts an unphysical equilibrium state with negative densities and does not correctly model large deviations of chemical fluctuations.¹² 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 A₂ molecules undergoing a dimerization reaction

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

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

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

and the entropy of the system is

Hence, we obtain the equilibrium distribution

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

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

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

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

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

where S^eq denotes the entropy at x^eq and N = ∑_sN_s. We can further approximate S_Stirling up to second order in $δxs=xs−xseq$ (s = 1, …, N_spec − 1), to get a Gaussian approximation to the Einstein distribution. This Gaussian approximation is described by linearized FHD, and we study it in more detail, including spatial dependence, next.

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

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

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

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

where ν = η/ρ₀, $D=−D12$⁠, and $μw$ is the second order derivative of Gibbs free energy with respect to concentration $w$⁠, given as $μw=kBT/[mw¯(1−w¯2)]$ for an ideal mixture. The linearized reaction term is denoted by $ρ0−1mΩ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,²⁷ 

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

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

where the linearized reaction rate is

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

One of the common assumptions in traditional reaction-diffusion modeling is that each chemical species is dilute and thus diffuses independently of other species. In this section, we explain how our formulation simplifies in the dilute limit. We consider a solution where all solute species are dilute (i.e., x_s ≪ 1), but the solvent is possibly a homogeneous mixture. We use index s here to denote only solute species. In the dilute limit, γ_s → 1 and the solute number densities are linearly proportional to their mole fractions, $ns≈(ρ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 n_s,

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

where $kr±$ denote the reaction rate constants.

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

where Ω_s represent the stochastic chemistry terms based on the LMA (35). Therefore, in the absence of fluid flow, our formulation is reduced to our previous reaction-diffusion model⁴ in the dilute limit. Note that jump processes and diffusion processes are combined in (36) and the time evolution of the probability distribution of n = {n_s} can be described by the differential Chapman–Kolmogorov equation.¹³ 

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;³⁵ 

• reduces to our previous method for reaction-diffusion systems⁴ 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.²⁴ 

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.²⁴ However, in order to prevent unphysical oscillations in mass densities in high Péclet number flows with sharp gradients, we can also use the Bell–Dawson–Shubin (BDS) second-order Godunov advection scheme.^37,38 We note that BDS advection adds artificial dissipation and does not obey a fluctuation-dissipation principle but is necessary for simulations where centered advection would fail due to insufficient spatial resolution. All simulations in this paper use centered advection unless otherwise noted. The discretization of the momentum equation is the same as our previous work.^3,26 We allow for periodic boundary conditions, impermeable walls, and no-flow reservoirs^27,36 held at fixed concentrations.

The first modification relative to our previous work³ 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,⁴ 

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

Specifically, we first convert cell-centered mass fractions to number densities, then apply ñ to obtain face-centered number densities, and finally convert these back to mass fractions that are used to compute $m¯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., n_iΔV ≤ 1), which prevents the occurrence of negative number densities. In most cases of interest, small numbers of molecules per cell correspond to dilute species. For dilute species [see (36)], the validity of using ñ has been justified in Ref. 4. In Sec. IV B, we numerically confirm that our approach is robust even when the total number of molecules in a cell is $O(10)$⁠.

Another key modification is the computation of the diffusion matrix χ for the deterministic and stochastic mass fluxes in the absence of some species or in the presence of vanishing species. In the vanishing limit, where one or more concentrations become zero, the diffusion matrix χ is not well conditioned since the corresponding diagonal component χ_ss diverges. This can cause numerical issues when one attempts to compute $F¯$ and $F̃$ since they depend on Wχ and $Wχ12$⁠, respectively. Unlike χ, however, the matrices Wχ and WχW are well defined in the vanishing limit, and we can construct Wχ 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 Wχ 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 Wχ in the vanishing limit, we introduce a mapping, $m(i)$⁠, used to expand/contract a subsystem matrix to/from a full matrix. For example, in a 6-species system having vanishing species $w2$ and $w4$⁠, we have $m(i)=(1,0,2,0,3,4)$⁠, i = (1, …, 6). As graphically illustrated using the 6-species system in Fig. 1, there are four cases to consider when one populates (Wχ)_ij,

Note that color names in the parentheses in (39) correspond to the colors in Fig. 1. A derivation of (39) and (40) is presented in  Appendix A. The full matrix $Wχ12$ can be obtained from the Cholesky decomposition of the symmetric matrix WχW. We note that if species s is vanishing, then (WχW)_is = (WχW)_sj = 0, so no stochastic mass flux is generated for species s.

We note that for each vanishing species only the diagonal element of Wχ remains nonzero. Hence, the diffusion of a dilute species s (⁠$ws$ ≪ 1) becomes decoupled from other species [see (A8)] and the effective diffusion coefficient D_s in (40) corresponds to the trace diffusion coefficient of s in the given fluid mixture. It is also important to note that the construction (39) guarantees that $F¯s=F̃s=0$ for vanishing species and ensures the mass conservation condition over all species, $∑s′F¯s′=∑s′F̃s′=0$⁠. Therefore, this procedure is robust to roundoff errors.

In our double-precision implementation, we treat any species s with $ws$ < 10⁻¹⁴ 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,

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.⁴ As explained below in detail, we treat viscous momentum dissipation implicitly and species diffusion explicitly. This is because in liquids the time step size is limited by the viscous Courant–Friedrichs–Lewy (CFL) condition (i.e., momentum diffusion) due to the large Schmidt number.

To combine a second-order midpoint tau-leap reaction sampling^29,39 with a predictor-corrector scheme for FHD, we adopt mass density updates from the ExMidTau (explicit midpoint tau leaping) scheme we previously developed for reaction-diffusion systems.⁴ 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 ν = η/ρ₀, unlike the trapezoidal scheme used in Refs. 3 and 26.

We advance the system from time t_n = nΔt to time t_n+1 = (n + 1)Δt in four steps:

1. Perform a predictor Stokes solve for the velocity $υn+1,*$ at t_n+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 t_n+1.

4. Perform a corrector Stokes solve for velocity $υn+1$ at t_n+1.

These steps are elaborated in detail in Algorithm 1. In the Algorithm 1 description, superscripts are used to denote the time level where a given quantity is evaluated, e.g., fⁿ = 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).

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,²⁸ as we have done previously.^3,26 The difference between the predictor and corrector Stokes solves is the temporal discretization of the advective term (explicit vs. trapezoidal) and the forcing term (explicit vs. midpoint); both Stokes solves are required for second-order deterministic accuracy.

As mentioned above, Steps 2 and 3 of the present scheme become essentially the same as the ExMidTau scheme in the dilute limit in the absence of advection. The only difference is a Stratonovich-type update of the stochastic mass flux in (47). While our previous analysis for RDME systems⁴ adopted the Ito interpretation, we choose the Stratonovich-type update here since a general analysis for weak fluctuations (linearized FHD)⁴⁰ 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

Here k ≡ k_⊥ is a wavevector in the plane perpendicular to the gradient, g is the gravitational acceleration, ζ = ρ⁻¹(∂ρ/∂$w$⁠) is the solutal expansion coefficient, and h is the concentration gradient, ∇$w$ = he_y. Using the method developed in Ref. 24, we analytically compute the resulting structure factors when our temporal integrator is used to solve (50a). For the non-reactive case (r = 0), we also compute structure factors obtained from two schemes developed in our previous work.^3,26 The overdamped scheme (see Algorithm 2 in Ref. 26) uses the steady Stokes equation, i.e., it eliminates the inertial term ∂_t$υ$ = 0 by taking an overdamped limit. We refer to the previous scheme for solving the inertial equation as the inertial trapezoidal scheme (see Algorithm 1 in Ref. 26) and to our new scheme as the inertial midpoint scheme (see Algorithm 1).

We set Δx = 1 and k_BT/ρ₀ = 1. To denote how fast momentum diffusion, species diffusion, and reaction are, we define the following dimensionless Courant numbers:

To consider the case of a relatively large Δt with a large Schmidt number Sc = 10³ (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 $μw$h²Δt² = 100, and gζhΔt² = 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⁻⁴ power law in the entire range of wavenumbers k,

However, the power law is suppressed at small k by gravity²⁷ or reaction.¹² 

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

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

Our midpoint scheme reproduces $Sw,wneq$ accurately for both rate constants. We emphasize that these results are remarkable given that α = O(10²). Our new scheme remains accurate for α = βSc ≫ 1 because the relative error in $Sw,wneq$ for our midpoint scheme has the form $[1+O(Sc−1)]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 is particularly dilute, with only ∼10 molecules per computational cell, whereas there are ∼10⁷ glucose and fructose molecules and ∼10¹⁰ water molecules per cell. We investigate the equilibrium distribution of the number of sucrose molecules in a cell to demonstrate that our approach correctly models the dilute limit.

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

Since we consider the dilute limit, we assume that the system is an ideal mixture and obeys the traditional LMA, with a forward rate a⁺ = k⁺n₃, a reverse rate a⁻ = k⁻n₁n₂, and an equilibrium constant $K=n1eqn2eq/n3eq$⁠. The reaction equilibrium lies almost completely in the direction of the formation of glucose and fructose,⁴⁵ but uncatalyzed sucrose hydrolysis is extremely slow (with a half-life of 500 years).⁴⁶ While we use an experimental value of K,⁴⁵ 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⁻⁴ and periodic boundary conditions. The thickness of the system is Δz = 10⁻⁴ 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×10−03$⁠, $w3eq=5.7×10−09$⁠, and $w4eq=0.990$⁠. We use two time step sizes, Δt = 10⁻⁵ and 10⁻⁴, to check the continuous-time limit and quantify time integration errors. Note that the larger Δt corresponds to diffusive CFL numbers D_s,maxΔt/Δx² = 0.07 for species diffusion and νΔt/Δx² = 100 for momentum diffusion. For each value of Δt, we ran 16 independent samples up to time $T=1$⁠, collecting data every t = 10⁻⁴ for $t≥T/10$⁠.

We recall that the number of sucrose molecules in a cell N = n₃ΔV has a continuous range in FHD simulations. We define its discrete distribution as $P(N)=∫N−12N+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 P_Poisson(N), and ρ(N) is essentially zero for negative values of N. We note that P_Poisson(N) is significantly different from its Gaussian approximation P_Gauss(N). For the larger Δt, while the remarkable agreement with the Poisson distribution is still observed, negative values of N start to appear, yielding $∫−∞0ρ(N)dN≈3×10−5$⁠; see the inset of Fig. 3. The same results were obtained in our previous reaction-diffusion model of dilute solutions,⁴ confirming that our treatment of the stochastic mass flux coefficients (see Sec. III A) is consistent with the dilute limit, even in the presence of random advection. We have also confirmed that the equilibrium structure factor of each species (not shown) has a flat spectrum (as predicted by theory³), indicating that there are no spurious correlations between cells.

We next consider a liquid mixture undergoing the dimerization reaction (20). This binary system contains monomers A (s = 1) and dimers A₂ (s = 2) and is representative of cyclic dimer formation in pure liquid acetic acid. We demonstrate here our ability to model a system with strong fluctuations in the absence of a dominant solvent by considering a small number of molecules (∼10) of each species per cell. As in the sugar solution example, we investigate the equilibrium distribution of monomers and dimers; however, since the system is not dilute, the distribution of each species is not Poisson. The numbers of monomers and dimers in a cell (N₁ and N₂) do not vary independently due to the constant density assumption N₁ + 2N₂ = ρ₀ΔV/m, where m is the mass of a monomer. Therefore, we investigate the equilibrium distribution of N₂, P(N₂), for N₁ + 2N₂ = 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 k_B = 1. We set ρ₀ = 40 with $w1eq=w2eq=0.5$ so that N₁ + 2N₂ = 40 with $N1eq=20$ and $N2eq=10$⁠. We set the reaction rates a^± as in (26) with a modification

where N⁺ = max(N, 0). Note that (56) turns off unphysical reactions when 0 < N₁ < 1. 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 η = 10³ and $T=103$⁠. We use a small Δt = 10⁻² to minimize temporal integration errors. For 16 independent samples with 10⁵ time steps, we collect data every 10² time steps, discarding the first 10⁴ time steps.

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

It is important to note that statistically identical results for P(N₂) are obtained from the corresponding non-reactive system with κ^± = 0 (not shown). This confirms thermodynamic consistency of our overall formulation. In addition, this also confirms the validity of our overall numerical treatment for diffusion with strong fluctuations. In particular, considering that our multiplicative GWN modeling for strong fluctuations was developed in the dilute limit⁴ 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,w1−Sw1,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;¹² here we consider a liquid mixture with a large Schmidt number Sc = 10³, which quantitatively changes the spectrum of giant fluctuations.

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

We consider a two-dimensional square domain of side length L_x = L_y = 64 (in arbitrary units) and periodic boundary conditions in the x-direction. The system is divided into 128 × 128 grid cells with grid spacing Δx = Δy = 0.5. To remain in the linearized FHD regime, we increase the cell depth to Δz = 10⁶ so that there are sufficiently many monomers and dimers in a cell, N₁ + 2N₂ = 2.5 × 10⁵ for ρ₀ = 1 and m = 1. We set $−D12=1$⁠, η = 10³, and k_BT = 10³. 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⁻², 0.1125), corresponding to linearized reaction rate r = 0.4 and penetration depth $ξ=10Δx$⁠, and (κ⁺, κ⁻) = (8.438 × 10⁻³, 1.125 × 10⁻²), 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 νΔt/Δx² = 100. We ran 10⁵ steps discarding the first 10⁴ 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 conditions²⁷ are imposed with (⁠$w1$⁠, $w2$⁠) = (0.49, 0.51) at y = 0 and (⁠$w1$⁠, $w2$⁠) = (0.51, 0.49) at y = L_y. While a linear concentration profile is formed in the steady state for the non-reactive case, a nonlinear profile is generated by the dimerization reaction. In Fig. 5(a), the profiles of $w¯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 = L_y/2 is constrained to be in chemical equilibrium due to faster reactions, resulting in larger concentration gradients at the boundaries. Identical concentration profiles are obtained from the corresponding deterministic reaction-diffusion systems (not shown).

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

is used instead of k_x. For the non-reactive case (r = 0), a clear k⁻⁴ power law is observed until the confinement effect becomes significant for small $kx≪Ly−1$⁠. For the reactive cases, the k⁻⁴ power law is only observed at large $kx≫r/−D12$⁠. For larger r, the k⁻⁴ power law appears in a narrower range of k_x values and the prefactor of the power law becomes larger.

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