We develop a low Mach number formulation of the hydrodynamic equations describing transport of mass and momentum in a multispecies mixture of incompressible miscible liquids at specified temperature and pressure, which generalizes our prior work on ideal mixtures of ideal gases [Balakrishnan *et al.*, “Fluctuating hydrodynamics of multispecies nonreactive mixtures,” Phys. Rev. E **89** 013017 (2014)] and binary liquid mixtures [Donev *et al.*, “Low mach number fluctuating hydrodynamics of diffusively mixing fluids,” Commun. Appl. Math. Comput. Sci. **9**(1), 47-105 (2014)]. In this formulation, we combine and extend a number of existing descriptions of multispecies transport available in the literature. The formulation applies to non-ideal mixtures of arbitrary number of species, without the need to single out a “solvent” species, and includes contributions to the diffusive mass flux due to gradients of composition, temperature, and pressure. Momentum transport and advective mass transport are handled using a low Mach number approach that eliminates fast sound waves (pressure fluctuations) from the full compressible system of equations and leads to a quasi-incompressible formulation. Thermal fluctuations are included in our fluctuating hydrodynamics description following the principles of nonequilibrium thermodynamics. We extend the semi-implicit staggered-grid finite-volume numerical method developed in our prior work on binary liquid mixtures [Nonaka *et al.*, “Low mach number fluctuating hydrodynamics of binary liquid mixtures,” arXiv:1410.2300 (2015)] and use it to study the development of giant nonequilibrium concentration fluctuations in a ternary mixture subjected to a steady concentration gradient. We also numerically study the development of diffusion-driven gravitational instabilities in a ternary mixture and compare our numerical results to recent experimental measurements [Carballido-Landeira *et al.*, “Mixed-mode instability of a miscible interface due to coupling between Rayleigh–Taylor and double-diffusive convective modes,” Phys. Fluids **25**, 024107 (2013)] in a Hele-Shaw cell. We find that giant nonequilibrium fluctuations can trigger the instability but are eventually dominated by the deterministic growth of the unstable mode, in both quasi-two-dimensional (Hele-Shaw) and fully three-dimensional geometries used in typical shadowgraph experiments.

## I. INTRODUCTION

The fluctuating hydrodynamic description of binary mixtures of miscible fluids is well-known,^{1,2} and has been used successfully to study long-ranged non-equilibrium correlations in the fluctuations of concentration and temperature.^{1} Much less is known about mixtures of more than two species (multicomponent mixtures), both theoretically and experimentally, despite their ubiquity in nature and technological processes. Part of the difficulty is in the increased complexity of the formulation of multispecies diffusion and the increased difficulty of obtaining analytical results, as well as the far greater complexity of experimentally measuring transport coefficients in multispecies mixtures. In fact, experimental efforts to characterize the thermo-physical properties of ternary mixtures are quite recent and rather incomplete.^{3}

At the same time, many interesting physical phenomena occur only in mixtures of more than two species. Examples include diffusion-driven gravitational instabilities that only occur when there are at least two distinct diffusion coefficients,^{4,5} as well as reverse diffusion, in which one of the species in a mixture of more than two species diffuses in the direction opposite to its concentration gradient. Another motivation for this work is to extend our models and numerical studies to chemically reactive liquid mixtures;^{6} interesting chemical reaction networks typically involve many more than two species. Giant nonequilibrium thermal fluctuations^{1,7,8} are expected to exhibit qualitatively new phenomena in multispecies mixtures due to their coupling with phenomena such as diffusion- and buoyancy-driven instabilities. Due to the difficulty in obtaining analytical results in multispecies mixtures, it is important to develop computational tools for modeling complex flows of multispecies mixtures. In previous work, we developed a fluctuating hydrodynamics finite-volume solver for ideal mixtures of ideal gases^{9} and studied giant fluctuations, diffusion-driven instabilities, and reverse diffusion in gas mixtures. In practice, however, such phenomena are much more commonly observed and measured experimentally in non-ideal mixtures of liquids. It is therefore important to develop fluctuating hydrodynamics codes that take into account the large speed of sound (small compressibility) of liquids, as well as the non-ideal nature of most liquid solutions and mixtures. While thermal transport does of course play a role in liquid mixtures as well, it is often the case that experimental measurements are made isothermally or in the presence of steady temperature gradients, or that temperature and concentration fluctuations essentially decouple from each other.^{10} In this work, we extend our work on isothermal low Mach number fluctuating hydrodynamics for non-ideal binary mixtures of liquids^{11,12} to multispecies mixtures.

Transport in multispecies fluid mixtures is a topic of great fundamental and engineering importance, and has been studied in the field of nonequilibrium thermodynamics^{13} and chemical engineering^{14} for a long time. While the case of a binary fluid mixture is well-understood and the complete hydrodynamic equations are well-known, including thermal fluctuations,^{1,2} there are some inconsistencies in the literature regarding the treatment of multispecies diffusion. Here, we combine several sources together to obtain a formulation that is, to our knowledge, new, although all of the required pieces are known. In particular, our focus here is on formulating the equations in a way that: (1) fits the framework of nonequilibrium thermodynamics, notably, the GENERIC framework;^{15} (2) is amenable to computer simulations for large numbers of species; (3) allows for straightforward inclusion of thermal fluctuations; (4) includes all standard mass transport processes (Fickian diffusion, thermodiffusion, and barodiffusion) and applies to non-ideal mixtures of non-ideal fluids; (5) expresses diffusive fluxes in a barycentric (center of mass) frame to allow seamless integration with the Navier-Stokes equations; and (6) is based on the Maxwell-Stefan (MS) (rather than Fickian) formulation of diffusion.^{14,16,17} These requirements inform our choice of the different elements of the formulation from different sources.

Our primary source is the recent monograph by Kuiken,^{16} which contains a nearly complete formulation except for some confusion between ideal and non-ideal mixtures that we clarify later on. Another primary source we rely on is the book of Krishna and Taylor.^{14} This book, like many other sources in chemical engineering, rely on separating one of the *N* species as a special “solvent” species in order to eliminate the redundancy in the description. While this simplifies the analytical formulation to some extent, it breaks the inherent symmetry of the problem by singling out a species. This complicates the numerical implementation and does not work well when the reference species vanishes. Following our prior work on mixtures of ideal gases,^{9} we rely on the monograph by Giovangigli^{18} to develop a formulation that treats all species equally and deals with the redundancy by using linear algebra techniques, see also a recent mathematical analysis and generalization to non-ideal gas mixtures.^{19} Here, we amend this formulation to account for non-ideality of the mixture, based in large part in Ref. 14 but now rewritten in terms of *N* rather than *N* − 1 variables. For thermodiffusion, we rely on the formulation in the book of Kjelstrup *et al.*,^{17} instead of that of Kuiken, in order to be in agreement to the standard definition of Soret and thermodiffusion coefficients for binary mixtures. Thermal fluctuations are formulated by fitting the formulation into the GENERIC framework, relying closely on the work of Öttinger.^{20}

In the remained of this section, we will define some notation and introduce the low Mach number formulation. In Sec. II, we discuss the formulation of both the deterministic and stochastic mass fluxes. By combining the proposed formulation of the mass fluxes with the low Mach number Navier-Stokes equations, we obtain a description of the fluctuating hydrodynamics of quasi-incompressible mixtures of incompressible miscible liquids. In Sec. III, we introduce and validate a discretization of the resulting system of equations using finite-volume methods, which exactly enforces the low Mach number quasi-incompressibility constraint.^{21} The method treats viscosity implicitly, allowing us to study flows over a broad range of Reynolds numbers, including steady Stokes flow. In Sec. IV, we use our algorithms to numerically study the development of diffusion-driven gravitational instabilities in a ternary solution of sugar and salt in water. We find a favorable comparison between our numerical results and recent experimental measurements in a Hele-Shaw cell.^{5} We propose that shadowgraph measurements in a different geometry may yield additional information about the possible coupling between the nonlinear instability and the giant concentration fluctuations that develop due to the presence of sharp gradients at the fluid-fluid interface.

### A. Basic notation

Our notation is based closely, though not entirely, on the work of Kuiken. We avoid the use of molar quantities and instead rely on “per molecule” equivalents, and also prefer the term “species” over “component.” Vectors (both in the geometrical and in the linear algebra sense), matrices (and tensors), and operators are denoted with bold letters. A diagonal matrix whose diagonal is given by a vector is denoted by the corresponding capital letter, for example, $X=Diag x $ implies *X _{ij}* =

*x*. The vector of

_{i}δ_{ij}*N*partial mass densities is $\rho w= \rho 1 , \u2026 , \rho N $, giving the total mass density $\rho = \u2211 i = 1 N \rho i $. The partial number densities are denoted with

*n*=

_{k}*ρ*/

_{k}*m*, where

_{k}*m*is the molecular mass of species

_{k}*k*, with total number density $n= \u2211 i = 1 N n i $.

The mass fractions are denoted with **w**, *w _{k}* =

*ρ*/

_{k}*ρ*, while the number or mole fractions are denoted with

**,**

*x**x*=

_{k}*n*/

_{k}*n*; both the mass and number fractions must sum to unity,

**1**

^{T}

**=**

*x***1**

^{T}

**w**= 1, where

**1**denotes a vector of ones. One can transform between mass and number fractions by

where the mixture-averaged molecular mass is

A useful formula that we will use later is the Jacobian of transforming from mass to number fractions

where $W=Diag w $.

### B. Low Mach number hydrodynamic equations

The hydrodynamics of miscible mixtures of incompressible liquids can be described using low Mach number equations, as explained in more detail in our prior works.^{11,12} The low Mach number equations can be obtained by making the ansatz that the thermodynamic behavior of the system is captured by a reference pressure $P r , t = P 0 r $, with the additional pressure contribution $\pi r , t =O Ma 2 $ capturing the mechanical behavior while not affecting the thermodynamics, where Ma is the Mach number. The reference pressure is determined from the condition of hydrostatic equilibrium in the absence of flow. In a gravitational field, the reference state is stratified and the reference pressure is in hydrostatic balance, ∇*P* = *ρ*** g**, where

**is the gravitational acceleration (see Ref. 22 for details of the construction of these types of models). In this work, we will assume that the reference pressure gradients are weak so that the thermodynamic properties of the system can be evaluated at a reference pressure**

*g**P*

_{0}that does not depend on space and time.

We will focus here on systems for which the temperature $T r , t = T 0 r $ is *specified* and not modeled explicitly. A constant-temperature model is appropriate, for example, if the system is in contact with an external heat reservoir and the thermal conductivity is sufficiently large to ensure a constant temperature (e.g., a constant temperature gradient) is maintained, and the Dufour effect is negligible. In the context of fluctuating hydrodynamics, one can often argue that temperature (more precisely, energy) fluctuations decouple from concentration fluctuations^{1} and one can model mass and energy transport separately. It is possible to extend our low Mach number formulation to include energy transport,^{23} as we will consider in future work. Here, we simply account for mass transport due to imposed fixed thermodynamic pressure and temperature gradients in the form of barodiffusion and Soret mass fluxes, but do not model the evolution of the thermodynamic pressure and temperature explicitly.

In low Mach number models, the total mass density $\rho w ; T 0 , P 0 $ is a specified function of the local composition at the given reference pressure and temperature. Here, we consider mixtures of incompressible liquids that do not change density upon mixing. A straightforward multispecies generalization of the binary formulation we proposed in Refs. 11 and 12 is given in Eq. (2.1) in Ref. 24, and takes the form of an equation of state (EOS) constraint

where $ \rho \u0304 i $ are the pure-component densities, which we will assume to be specified constants. We note that even if the specific EOS (2) is not a very good approximation over the entire range of concentrations, it may nevertheless be a very good approximation over the range of concentrations of interest if the densities $ \rho \u0304 i $ are adjusted accordingly. In this case, $ \rho \u0304 i $ are to be interpreted not necessarily as pure component densities, since some of the components may not even exist as fluid phases at the reference pressure and temperature, but rather, as numerical parameters describing locally the dependence of the mass density on the composition at the specified reference pressure and temperature.

The *low Mach number* equations for the center of mass velocity $v r , t $, the mechanical component of the pressure $\pi r , t $, and the partial densities $ \rho 1 r , t , \u2026 , \rho N r , t $ of a multispecies mixture of *N* fluids can be written in conservation form as^{11}

where $\eta w ; T 0 , P 0 $ is the viscosity, $ \u2207 \u0304 =\u2207+ \u2207 T $ is a symmetric gradient, and ** g** is the gravitational acceleration. Note that the bulk viscosity terms have been absorbed into the pressure

*π*in the low Mach formulation.

^{11}Thermal fluctuations are accounted for through the stochastic momentum flux

**Σ**, formally modeled as

^{1,25}

where *k _{B}* is Boltzmann’s constant and

**𝒲**(

**,**

*r**t*) is a standard white noise Gaussian tensor field with uncorrelated components,

Here, $F= F 1 , \u2026 , F N $ is a composite vector of diffusive deterministic, $ F \xaf $, and stochastic, $ F \u02dc $, fluxes in the barycentric (center of mass) frame, where $ F i = F \xaf i + F \u02dc i $ is the flux for species *i*. We will use a compact matrix notation in which we can write the mass conservation equations without subscripts,

The diffusive fluxes preserve mass conservation because they sum to zero (**1**^{T}** F** = 0 in matrix notation),

which ensures that the total mass density obeys the usual continuity equation

Differentiating EOS constraint (2) in time, we get

## II. DIFFUSION

In this section, we develop a formulation of the diffusive fluxes in the barycentric frame for a non-ideal mixture, in a manner suitable for numerical modeling. Nonequilibrium thermodynamics expresses the deterministic diffusive fluxes in terms of the thermodynamic driving force (gradients of the chemical potential); in short-hand matrix notation,^{9}

where $ \mu k x , T , P $ is the chemical potential of species *k* and **∇**_{T} denotes a gradient at constant temperature. It is important to note that we use chemical potential per unit mass,^{13} which differs from more commonly used chemical potential per mole^{14} by a factor of *m _{k}N_{A}*, where

*N*is Avogadro’s number. The matrix of Onsager coefficients

_{A}**is symmetric (by Onsager’s reciprocity) and positive-semidefinite (to ensure dissipation, i.e., positive entropy production), and has zero row and column sums (to ensure mass conservation (7)), see (26) for the explicit form. The vector of thermal diffusion ratios**

*L***also sums to zero to ensure mass conservation (7), see (27) for the explicit form.**

*ξ*In order to make (9) suitable for computation, we need to express the gradients of the chemical potential and the Onsager matrix in terms of more readily computable quantities. The gradient of chemical potential at constant temperature can be expressed in terms of the gradient of composition and pressure using the chain rule,

We now explain how to relate ∂** μ**/∂

**, ∂**

*x***/∂**

*μ**P*, and

**to more familiar thermodynamic and transport quantities. This will allow us to express the deterministic component of**

*L***as a function of the local gradients of composition, temperature, and pressure (see (25) for the final result), and will provide us with a model for the stochastic mass fluxes (see (30) for the final result).**

*F*### A. Chemical potentials

We use the specific (per mass) Gibbs density $g w , T , P =u\u2212Ts+Pv$ as the thermodynamic potential, where *u* is the specific internal energy density, *s* the specific entropy density, and *v* = *ρ*^{−1} is the specific volume. The chemical potentials per unit mass are ** μ** = ∂

*g*/∂

**w**. For non-ideal mixtures, we can express chemical potentials as a sum of ideal and excess contributions,

where $ \mu k 0 T , P $ is a reference chemical potential (e.g., pure liquid state at standard conditions) and $ \gamma k x , T , P $ is the activity coefficient of species *k*; for an ideal mixture, *γ _{k}* = 1. In the low Mach number setting we consider here, the chemical potentials depends on pressure

*only*through the reference state. This is always true for an ideal mixture, but may be assumed more generally so long as the activities only depend on composition and

*not*on pressure.

#### 1. Thermodynamic factors

Note that all thermodynamic functions are in principle only defined for valid compositions, **1**^{T}**w** = 1; however, any analytical extension of these functions can be used to work with unconstrained derivatives instead of the more traditional constrained derivatives. The matrix of thermodynamic factors **Γ** is defined via the *unconstrained* derivatives

which we can write in component form as

We note that the thermodynamic factors are incorrectly defined in the book by Kuiken^{16} to have pressure *P* in the denominator instead of *nk _{B}T* (the two are of course equal for ideal gases). Using the matrix of thermodynamic factors, we can express the contribution to the gradient of the chemical potentials in terms of gradients of composition,

Note that the constraint $ \u2211 i = 1 N x i =1$ is automatically taken into account since $ \u2211 i = 1 N \u2207 x i =0$ and the component of the unconstrained derivatives normal to the constraint does not actually matter.

For nonideal (dense) gas mixtures, it is possible to relate **Γ** to the equation of state, see Ref. 19 for example, calculations for a dense-gas EOS. In order to model the thermodynamic factors as a function of composition in liquid mixtures, several different models have been defined and experimentally parameterized, such as the Wilson, NTLR, or UNIQUAC models, as described in more detail in Appendix D of the book by Krishna and Taylor.^{14} These models are all based on the normalized excess Gibbs energy density per particle $ g \u0303 ex x , T , P $. Converting this to specific excess Gibbs energy density, we can write

giving the thermodynamic factors in the form

By converting the second-order derivatives with respect to **w** to the more traditional derivatives with respect to ** x** by using Jacobian (1), we obtain the final relation

where the symmetric matrix ** H** = ∂

^{2}

*g*

_{ex}/∂

*x*^{2}is the Hessian of the excess free energy per particle.

In the neighborhood of a stable point (far from phase separation), the total Gibbs energy density is locally a convex function of composition. By also including the ideal contribution to the free energy, it is not hard to show that this implies the stability condition

where one of the eigenvalues is always zero and the rest must be non-negative. The physically key quantity required to model the thermodynamics of non-ideal mixtures is ** H**, rather than the more traditional

**Γ**. To avoid a large departure from the literature, we continue to use

**Γ**but we note that in our numerical codes, the input is

**and**

*H***Γ**is calculated from (13). If one tries to model

**Γ**directly, it is difficult to ensure the correct symmetry structure, which is obscured in

**Γ**but directly evident in

**. We therefore disagree with statements in the literature that it is more accurate to use models for activities (first derivatives) than to use models for the excess free energy and then take second derivatives. While the former may indeed be more accurate, it may also lead to inconsistent thermodynamics; for thermodynamic consistency, one must model the excess free energy as a function of composition.**

*H*#### 2. Partial volumes

In order to compute (10), we express the partial volumes ** θ** = ∂

**/∂**

*μ**P*by using a Maxwell relation,

where $v w , T , P = \rho \u2212 1 $ is the specific volume. For a mixture of incompressible liquids given by EOS (2), the above relates $ \theta k = \rho \u0304 k \u2212 1 $ to the pure-component densities. Instead of partial volumes, we will use the volume fractions *φ _{k}* =

*ρ*. In ideal gas mixtures,

_{k}θ_{k}*φ*=

_{i}*x*and in the low Mach number setting, $ \phi k = \rho k / \rho \u0304 k $; note that $ \u2211 i = 1 N \phi i =1$.

_{i}### B. Diffusion driving force

The thermodynamic driving forces for diffusion are the chemical potential gradients, **∇**_{T}** μ**. Note, however, that the definition of the thermodynamic force is not unique, as becomes evident when we consider the average local entropy production rate due to mass diffusion,

where *a*_{i} is the acceleration of the particles of species *i* due to external fields (e.g., gravity or electric fields). Because the fluxes add to zero, we can add an arbitrary vector ** α** to all of the chemical potential gradients without changing the entropy production rate. Let us therefore write the entropy production rate as

where the above defines the thermodynamic driving force *d*_{k} for the diffusion of the *k*th species. Thermodynamic equilibrium corresponding to a vanishing of the entropy production rate, more specifically, to a vanishing of both the driving forces and the fluxes, *d*_{eq} = 0 and $ F \xaf eq =0$.

The fluxes above are defined in the barycentric frame. In order to determine the appropriate value of ** α**, let us consider transforming to a frame of reference that is moving with velocity

**v**

_{ref}relative to the center of mass of the mixture. This changes the fluxes to

*F*_{k}↦

*F*_{k}−

*ρw*

_{k}**v**

_{ref}and changes the entropy production rate by $ \rho / T v ref \u22c5 \u2211 i = 1 N d i $. This implies that if we want to have Galilean invariance of the entropy production rate, we should ensure that the driving forces sum to zero,

**1**

^{T}

**= 0.**

*d*If we take a gradient at constant temperature of both sides of the Gibbs-Duhem relation,

we get the relation

which shows that **1**^{T}** d** = 0 implies

In this work, we will only consider gravity, for which all species accelerations are equal to the gravitational acceleration, $ a k = \u2211 i = 1 N w i a i =g$.

This leads us to define the diffusion driving force as^{14,16,18}

Note that Kuiken^{16} puts pressure *P* in the denominator of the barodiffusion term instead of *nk _{B}T*, which leads to an inconsistency with the majority of the literature and the standard definition of the MS diffusion coefficients.

^{17}Note that each of the two terms in the driving force separately sums to zero, since

**1**

^{T}

**=**

*ϕ***1**

^{T}

**w**= 1 and

### C. Maxwell-Stefan description of diffusion

In order to compute the diffusive fluxes using (9), we need to relate the Onsager matrix ** L** to the more familiar MS diffusion coefficients. The Maxwell-Stefan relations are obtained by equating the driving force to the frictional force on a species due to the difference in its velocity relative to other species,

where

is the mass-averaged velocity of species *k* augmented by the thermodiffusion “slip.” Here, the symmetric matrix of MS *binary* diffusion coefficients ** D** has zero diagonal,

*D*= 0, and the off-diagonal elements are positive (there may be exceptions to this rule for ionic solutions

_{kk}^{26}) diffusion coefficients that have a physical interpretation of suitably dimensionalized inverse friction coefficients between

*pairs*of species. This positivity of

**ensures a positive entropy production, and thus consistency with the second law. It is observed that the MS diffusion coefficients show less variation with changes in composition than alternatives such as Fickian diffusion coefficients.**

*D*^{16,17}The off-diagonal elements of

**can therefore be interpolated as a function of composition relatively easily.**

*D*^{27–29}The thermodiffusive fluxes are expressed in terms of the thermodiffusion coefficients

*D*^{(T)}. Since only differences of $ D k ( T ) $’s appear, there are only

*N*− 1 thermodiffusion coefficients; in order to ensure that the mass fluxes sum to zero, $ \u2211 i = 1 N \rho i v i =0$, we require that $ \u2211 i = 1 N \rho i D i ( T ) =0$. This gives the constraint

which removes the redundancy in the specification of the thermodiffusion coefficients. (We thank an anonymous reviewer for pointing out this relation.)

We can write (18) in matrix form as

where the symmetric matrix **Λ** is defined via

It is relatively straightforward to show that **Λ** is positive semidefinite if *D _{ij}* > 0 for

*i*≠

*j*. Here, we introduced the vector of thermal diffusion ratios

where $ \u2211 i = 1 N \zeta i =0$ by construction. In our algorithm, the primary inputs are the MS diffusion coefficients ** D** and the thermodiffusion coefficients

*D*^{(T)};

**Λ**and

**are calculated from them.**

*ζ* which now relates the deterministic diffusive fluxes $ F \xaf $ with the gradients in composition, pressure, and temperature. By solving the above linear system for $ F \xaf $ subject to the condition $ \u2211 i = 1 N F \xaf i =0$, we can obtain a formula for the fluxes in terms of gradients of ** x**,

*P*, and

*T*. In order to carry out this computation, we follow Giovangigli,

^{18}where linear algebra tools are developed to solve linear system (23).

### D. Fick’s law

Let us introduce (Giovangigli^{18} attributes the introduction of **χ** to Ref. 30) the symmetric positive-semidefinite *diffusion matrix* **χ** as a pseudo-inverse of **Λ**,^{18} see Appendix A for more details,

where *α* > 0 is an arbitrary constant, for example, the choice $\alpha =Trace \Lambda $ guards against roundoff errors. One can directly compute **χ** using the above formula, but we will discuss numerical alternatives in Appendix A. While the MS diffusion coefficients are binary friction coefficients, the matrix **χ** is a multispecies construct that takes into account the composition of the mixture. One can in fact start the formulation from the matrix **χ**; however, we prefer to use the more-standard MS coefficients as input and compute **χ** from them. The reason behind this choice is the belief that the MS diffusion coefficients change more slowly with composition and thus are easier to tabulate and interpolate, than would be **χ**.

It can be shown^{18} that the solution of the linear system of equations (23) can be written in the Fickian form

This expression will be used in our numerical codes to compute the fluxes from the gradients in composition, pressure, and temperature. We use gradients of number fractions (composition) rather than gradients of chemical potential since the later is numerically ill-behaved due to the logarithmic divergence of the chemical potentials for nearly vanishing species. It is, however, also possible to isolate the singularity of the chemical potentials and to use the gradient of the non-singular part of the chemical potentials directly, instead of using **Γ** to convert to gradients of composition, as done in Ref. 31. Observe that the fluxes automatically add up to zero, $ 1 T F \xaf =0$, since $ 1 T W\chi = w T \chi = \chi w T =0$. Note that the expression inside the brackets in (25) adds to zero over all species, since **1**^{T}** ζ** = 0 and

**1**

^{T}

**= 0. It is important to preserve these “sum to zero” properties in spatial discretizations of Fick’s law, as we discuss in more detail in Sec. III B.**

*d*### E. Thermal fluctuations

To formulate the stochastic mass flux, let us first relate the diffusion matrix **χ** to the more familiar Onsager matrix. By comparing the Onsager and the Maxwell-Stefan expressions for the fluxes,

we can directly identify (see (2.17) in Ref. 19)

which makes it clear that the Onsager matrix is symmetric positive semidefinite (SPD) since **χ** is SPD. For the Soret effect, we can identify ** ξ** and

**as rescaled versions of each other**

*ζ*The fact ** L** is SPD by construction is crucial for adding thermal fluctuations (stochastic mass fluxes), since that requires the “square root” of the Onsager matrix, notably, a matrix $ L 1 2 $ that satisfies $ L 1 2 L 1 2 \u22c6 =L$, where star denotes an adjoint (transpose for real matrices or conjugate transpose for complex matrices).

^{1,9,15}It is easy to see that

meets this criterion, where $ \chi 1 2 \chi 1 2 \u22c6 =\chi $; for example, $ \chi 1 2 $ can be taken to the lower-triangular Cholesky factor of **χ**. In fluctuating hydrodynamics, we simply add a stochastic contribution to the mass flux of the form

where **𝒵** denotes a collection of *N* spatio-temporal white noise random fields (note that only *N* − 1 random fields are actually required since one of the eigenvalues of ** L** is zero), i.e., a random Gaussian field with correlations,

Observe that the stochastic fluxes sum to zero, $ 1 T F \u02dc =0$ because $ \chi 1 2 w=0$ follows from **χw** = 0.

We are finally in a position to write the complete equation for mass fractions (4),

In Appendix C, we demonstrate that stochastic mass fluxes (29) can also be derived by following the Maxwell-Stefan construction and augmenting the dissipative frictional forces between pairs of species by corresponding (Langevin) fluctuating forces. That formulation gives another physical interpretation to the stochastic mass fluxes, but is not as useful for computational purposes because the number of pairs of species (and thus stochastic forces) is much larger than the number of species, so in computations, we use (30).

Important quantities that can be derived from fluctuating equation (30) are the spectrum of the time correlation functions and the amplitude of the fluctuations at thermodynamic equilibrium, referred to as the dynamic and static structure factors, respectively. The matrix of equilibrium structure factors can be expressed either in terms of mass or mole fractions. Here, we define the matrix of static covariances in terms of the fluctuations in the mass fractions *δ***w** around the equilibrium concentrations **w**. The dynamic structure factor matrix $ S w k , t $ is defined as

where *i* and *j* are two species (including *i* = *j*), ** k** is the wavevector, hat denotes a Fourier transform, and star denotes a complex conjugate. The equal time covariance in Fourier space is the static structure factor $ S w k = \delta w \u0302 \delta w \u0302 \u22c6 $,

The equilibrium static factors were computed for a ternary mixture in Ref. 32. In Appendix D, we use (30) to obtain the equilibrium static structure factor for a mixture with an arbitrary number of species,

If stability condition (14) is satisfied then ** S_{w}**⪰

**0**will be symmetric positive semidefinite, as it must be since it is a covariance matrix. If the mixture is unstable then the above calculation is invalid because the fluctuations around the mean will not be small and linearized fluctuating hydrodynamics will not apply. In the low Mach number setting, the structure factor for density is

## III. NUMERICAL ALGORITHM

In this section, we give some details about our numerical algorithms, and then present some validation studies that verify the deterministic and stochastic order of accuracy of our schemes. In particular, we confirm that we can accurately model equilibrium and non-equilibrium concentration fluctuations in multispecies ternary mixtures. In Sec. IV, we use the algorithms described here to model the development of instabilities during diffusive mixing in ternary mixtures.

### A. Low Mach integrator

The numerical algorithms we use to solve multispecies low Mach number Eqs. (3), (5), and (30) are closely based on the binary mixture algorithms described in detail in Ref. 12. In particular, the spatial discretization of quasi-incompressible flow EOS constraint (5) and velocity equation (3), as well as the temporal integration algorithms, is identical to the binary case.^{12} Some of the key features of the algorithms developed in Refs. 11 and 12 are as follows.

We employ a uniform staggered-grid finite-volume (flux-based) spatial discretization because of the ease of enforcing the constraint on the velocity divergence (note that our compressible algorithm

^{9}uses a collocated grid) and incorporating thermal fluctuations.^{33}Our spatial discretization strictly preserves mass and momentum conservation, as well as the EOS constraint

^{11}(but see Sec. III A 1).By using the high-resolution Bell-Dawson-Shubin (BDS) scheme

^{34}for mass advection, we can robustly handle the case of no mass diffusion (no dissipation in (30)).Our temporal discretization uses a predictor-corrector integrator that treats all terms except momentum diffusion (viscosity) explicitly. We have developed two different temporal integrators one for

*inertial*momentum equation (3), and one for the viscous-dominated or*overdamped*limit^{35}in which the velocity equation becomes the steady Stokes equation.^{12}We treat viscosity implicitly without splitting the pressure update, relying on a recently developed variable-coefficient multigrid-preconditioned Stokes solver.

^{36}This makes our algorithms efficient and accurate over a broad range of Reynolds number, including the zero Reynolds number limit, even in the presence of nontrivial boundary conditions.

The key difference between binary mixtures^{11,12} and multispecies mixtures is the handling of density equation (8) and the computation of the diffusive and stochastic mass fluxes. In the binary case, the conserved variables we use are *ρ* and *ρ*_{1}, with the corresponding primitive variables being *ρ* and the mass fraction *c* ≡ *w*_{1} = *ρ*_{1}/*ρ*. In the multispecies case, our conserved variables are the partial densities *ρ _{k}*; the total density $\rho = \u2211 i = 1 N \rho i $ is computed from those as needed. The corresponding primitive variables are

*ρ*and

*w*. In the binary case, we expressed all of the diffusive fluxes in terms of gradients of mass fractions, but in the multispecies case, we rely on the more traditional formulation in terms of gradients of number (mole) fractions, and we also include

_{k}*x*as primitive variables. Further details on the computation of the multispecies diffusive and stochastic mass fluxes are given in Sec. III B.

_{k}#### 1. EOS drift

Our low Mach number algorithms are specifically designed to ensure that the evolution remains on the EOS constraint, i.e., that the partial densities or equivalently the density and the composition in each grid cell strictly satisfy (2).^{11} Nevertheless, due to roundoff error and finite solver tolerance in the fluid Stokes solver, a slow drift off the EOS constraint occurs over multiple time steps. To correct this, we occasionally need to project the state (partial densities) back onto the constraint.^{11} A similar projection onto the EOS is required in the BDS advection scheme for average states extrapolated to the faces of the grid.^{12}

For binary mixtures, we used a simple *L*_{2} projection onto the EOS. For mixtures of many species, some of the species may be trace species or not present at all, and in this case, it seems more appropriate to use a mass-fraction-weighted *L*_{2} projection step. Given a state $ \rho , w $ that does not necessarily obey (2), the weighted *L*_{2} projection consists of correcting *ρ _{k}* as follows:

where the correction is

which vanishes for species not present (*w _{k}* = 0). When performing a global projection onto the EOS, one should additionally re-distribute the total change in the mass of species

*k*over all of the grid cells to ensure that the projection step does not change the total mass of any species.

^{11}

### B. Diffusive and stochastic mass fluxes

The computation of the diffusive deterministic and stochastic mass fluxes for binary mixtures is described in detail in Ref. 11. We follow a similar but slightly different procedure for multispecies mixtures, primarily guided by the desire to make the algorithm efficient for mixtures of many species.

#### 1. Mixture model

The user input to our fluid dynamics code, i.e., the *mixture model*, is a specification of the required thermodynamic (e.g., non-ideality factors) and transport properties (e.g., shear viscosity) of the mixture as a function of state. The state of the mixture is described by the variables $ w , P , T $, or, equivalently, $ x , P , T $, where we recall that in our low Mach number model, the pressure and temperature are specified and not modeled explicitly, and the density is not an independent variable since it is determined from EOS constraint (2). Therefore, the mixture model in our low Mach code consists of specifying the thermodynamic and transport properties as a function of the composition **w**.

In the multispecies case, the mixture model requires specifying binary Maxwell-Stefan diffusion coefficients for each pair of species, i.e., the lower triangle of the matrix ** D**. Additional input is the vector of thermodiffusion coefficients

*D*^{(T)}(recall that only

*N*− 1 of these are independent since an arbitrary constant can be added to this vector), and the Hessian of the excess free energy per particle

**. MS diffusion coefficients can be interpolated as a function of composition using Vignes- or Darken-type formulas,**

*H*^{27–29}based on data obtained experimentally

^{3}or from molecular dynamics simulations.

^{37,38}The thermodynamics can be parametrized using Wilson, NTLR, or UNIQUAC models, and Hessian matrices

**can be computed from the formulas presented in Appendix D of the book by Krishna and Taylor,**

*H*^{14}based on experimental or molecular dynamics data.

^{39}We are not aware of any models for parameterizing the thermodiffusion coefficients as a function of composition in liquid mixtures. We note, however, that despite the availability of various mixture models, experimental efforts to obtain the parameters required in these models and compare various models are very recent. We are not aware of any mixture of more than two species for which there is reliable and reproducible data for the mass and thermal diffusion and thermodynamic coefficients even in the vicinity of a reference state, yet alone over a broad range of compositions.

From the mixture model input, i.e., *η*, ** D**,

*D*^{(T)}, and

**, we compute the following quantities. First, we obtain the matrix**

*H***Λ**using (21), and then from

**Λ**, we compute the diffusion matrix

**χ**using (A2), as discussed in more detail in Appendix A. We also compute the matrix of thermodynamic factors

**Γ**using (13), as well as the vector of thermal diffusion ratios

**using (22). These computations provide all of the matrices and vectors required to compute the non-advective mass fluxes in (30). We remind the reader that the species volume fractions**

*ζ***are easily computable for our model of a mixture of incompressible components, $ \phi k = \rho k \theta k = \rho k / \rho \u0304 k $.**

*ϕ*#### 2. Spatial discretization

The basic spatial discretization of the fluid equations and mass advection is unchanged from our previous work on binary mixtures^{11,12} and we do not discuss it further here. Here, we explain how we handle the diffusive and stochastic mass fluxes in the multispecies setting. The deterministic and stochastic mass fluxes are computed on the faces of the grid, and the divergence of the flux is computing using a conservative difference, in two dimensions

where $F= F \xaf + F \u02dc $.

Our spatial discretization of deterministic diffusive fluxes (25) closely mimics the one described in Section IV.A of Ref. 11, and is based on centered differences and centered averaging. In order to avoid division by zero in the absence of certain species in some parts of the domain, in each cell (i.e., for each cell center), we modify the densities $ \rho k \u2190max \u03f5 , \u2009 \rho k $ to be no-smaller than a small constant *ϵ* on the order of the roundoff tolerance; this modification is only done for the purpose of the diffusive flux calculation. In each cell, we compute the primitive variables *ρ*, **w**, and ** x** and then use the user-provided mixture model to compute

**Γ**(if the mixture is non-ideal),

**χ**,

**, and**

*ζ***. Next, in each cell (**

*ϕ**i*,

*j*), we compute the matrices

*ρ*

**and**

*W*χ**Γ**and the vectors

**/**

*ζ**T*, and $ \varphi \u2212 w / n k B T $. Then, we average (interpolate) these matrices and vectors to the faces of the grid using arithmetic averaging, for example,

and compute gradients of composition, pressure, and temperature using centered differences, for example,

Note that the key property $ \u2211 i = 0 N \u2207 x i =0$ is preserved, for example,

Finally, we compute the deterministic fluxes using (25) by a matrix-vector product, for example,

Note that the important properties of (25) discussed in Sec. II D are maintained by this discretization. To ensure mass conservation, it is crucial that the mass fluxes for different species add up to zero, for example, it must be that $ 1 T F i + 1 2 , j ( x ) =0$ on every *x* face of the grid. In the continuum formulation, this is true because **1**^{T}** Wχ** =

**w**

^{T}

**χ**= 0; it is not hard to show that the arithmetic averaging procedure used above preserves this property,

since **χw** = 0 in each cell. Similarly, the continuum properties **1**^{T}** ζ** = 0,

**1**

^{T}

**Γ∇**= 0, and $ 1 T \varphi \u2212 w / n k B T =0$ are preserved discretely due to their linearity and the linearity of the averaging process. This shows the importance of the linearity of the interpolation from the cell centers to the cell faces.

*x*Upon spatial discretization, the stochastic fluxes acquire a prefactor of $\Delta V \u2212 1 2 $ due to the delta function correlation of white-noise, where Δ*V* is the volume of a grid cell.^{40} This converts the spatio-temporal white-noise process $Z r , t $ into a collection of independent temporal white-noise processes $Y t $, one process for each face of the grid, for example,

In our code, we compute the Onsager matrix ** L** in every cell and then compute $ L 1 2 $ by Cholesky factorization; an equally good alternative is to compute $ \chi 1 2 $ by Cholesky factorization.

^{56}We then use arithmetic averaging to compute face-centered Cholesky factors $ L 1 2 $. An alternative procedure, which is likely better at maintaining discrete fluctuation-dissipation balance

^{41}but is the number of dimensions times more expensive, is to average

**to the faces, and then perform a Cholesky factorization on each face of the grid. Note, however, that achieving strict discrete fluctuation-dissipation balance requires expressing the fluxes in terms of the discrete gradients of the chemical potential, which is rather inconvenient and not numerically well-behaved. In this work, we chose to work with gradients of number fractions and thus only achieve discrete fluctuation-dissipation balance approximately.**

*L*### C. Numerical tests

In the deterministic setting, we have confirmed the second-order accuracy of our numerical method by repeating the lid-driven cavity test used in our previous work on binary mixtures.^{12} The essential difference is that the bubble being advected through a pure liquid of a first species in the lid-driven cavity is now composed of a mixture of two other species, making this a ternary mixture test. Our numerical results show little to no difference between the ternary and binary mixture cases, and show second-order pointwise deterministic convergence for our low Mach number scheme.

In this section, we focus on tests in the context of fluctuating hydrodynamics, in particular, we examine the matrix of dynamic structure factors $ S w k , t $ defined in (31), for a ternary mixture. We use the computer algebra system Maple to evaluate (33) and (D3) and obtain explicit formulas to which we compare our numerical results below. We also examine *S*_{ρ}, since, according to (34), by examining the fluctuations in density, we are examining the correlations among all pairs of species.

#### 1. Equilibrium fluctuations

One of the key quantities used to characterize the intensity of *equilibrium* thermal fluctuations is the static structure factor or static spectrum of the fluctuations. We perform these tests in the steady Stokes regime since the velocity fluctuations decouple from density fluctuations at equilibrium; the only purpose of the fluid solver at uniform equilibrium is to ensure that the density remains consistent with the composition.

In this equilibrium test, we use a ternary mixture with Stefan-Maxwell diffusion matrix and non-ideality matrix

in some arbitrary units in which *k _{B}* = 1. The molecular masses for the ternary mixture are

*m*

_{1}= 1.0,

*m*

_{2}= 2.0,

*m*

_{3}= 3.0 and the pure component densities are $ \rho \u0304 1 =2.0,\u2009 \rho \u0304 2 =3.0,\u2009 \rho \u0304 3 =3.857$. The system is a two-dimensional periodic system at equilibrium with equilibrium densities

*ρ*

_{1}= 0.6,

*ρ*

_{2}= 1.05,

*ρ*

_{3}= 1.35. At these conditions, the equilibrium density variance is $\Delta V\u2009 \delta \rho 2 = S \rho =0.3,$ where Δ

*V*is the volume of a grid cell. We employ a square grid of 32 × 32 cells with grid spacing Δ

*x*= Δ

*y*= 1 for these investigations, with the thickness in the third direction set to give a large Δ

*V*= 10

^{6}and thus ensure consistency with linearized fluctuating hydrodynamics. A total of 2 × 10

^{4}time steps are skipped in the beginning to allow equilibration of the system, and statistics are then collected for an additional 10

^{6}steps.

At equilibrium, the static structure factors are independent of the wavenumber due to the local nature of the correlations, $ S w i , j k = S eq , \u2009 w i , j =const.$ Since we include mass diffusion using an explicit temporal integrator, for finite time step sizes Δ*t*, we expect to see some deviation from a flat spectrum at the largest wavenumbers (i.e., for *k* ∼ Δ*x*^{−1}).^{40,41} In Fig. 1, we show the spectrum of density fluctuations at equilibrium for two different time step sizes, a large time step size Δ*t* = 0.1 (left panel), and a smaller time step size Δ*t* = 0.05 (right panel). Since the largest eigenvalue of the diffusion matrix is around χ ≈ 2, the largest stable time step size is Δ*t*_{max} ≈ 0.12. As seen in the figure, for Δ*t* = 0.1, which is close to the stability limit, we see a significant enlargement of the fluctuations at the corners of the Fourier grid; when we reduce the time step by a factor of 2 and we reduce the error by a factor of around 8, consistent with the fact that the explicit midpoint method used in our overdamped algorithm^{12} is third-order accurate for static covariances.^{40} Therefore, in the limit of sufficiently small time steps, we will recover the correct flat spectrum, demonstrating that our equations and our numerical scheme obey a fluctuation-dissipation principle.

In the left panel of Fig. 2, we show numerical results for the dynamic structure factors $ S w ( i , j ) k , t $ for several *i* ≠ *j* for $k= 4 , 4 \u22c52\pi /L$, where *L* = 32 is the length of the square domain. Note that the factors for *i* = *j* are not statistically independent due to the constraint that mass fractions sum to unity, and are thus not shown. In the right panel of Fig. 2, we show numerical results for $ S \rho k , t $, given by (34), for several different wavenumbers $k= \kappa x , \kappa y \u22c52\pi /L$. We compare the numerical results to the theoretical prediction (D3), which is a sum of two exponentially decaying functions. Excellent agreement is seen between simulation and theory, demonstrating that our numerical method correctly reproduces both the statics and dynamics of the compositional fluctuations.

#### 2. Non-equilibrium fluctuations

Fluctuations in systems out of equilibrium are known to be long-range correlated and significantly enhanced compared to equilibrium. In particular, in the presence of an imposed (macroscopic) concentration gradient, concentration fluctuations exhibit a characteristic power-law static structure factor ∼*k*^{−4}.^{1} In Section IV.C in Ref. 9, we studied the long-ranged (giant) concentration fluctuations in a ternary mixture in the presence of a gradient imposed via the boundary conditions and confirmed that our multispecies compressible algorithm correctly reproduced theoretical predictions; here, we repeat this test but for a mixture of three incompressible liquids.

In order to simplify the theoretical calculations, see Appendix B in Ref. 9, we take the first two of the three species to be dynamically identical (indistinguishable), and take the molecular masses to be equal, *m*_{1} = *m*_{2} = *m*_{3} = 1.0 (this makes mass and mole fractions identical). The Stefan-Maxwell diffusion matrix is taken to be

and the mixture is assumed to be ideal, ** H** = 0, and isothermal,

**∇**

*T*= 0. In order to focus our attention on the nonequilibrium fluctuations, we set the stochastic mass flux to zero, $ F \u02dc =0$; this ensures that all concentration fluctuations come from the coupling to the velocity fluctuations via the gradient and eliminate the statistical errors coming from a finite background spectrum. The pure component densities are $ \rho \u0304 1 = \rho \u0304 2 = \rho \u0304 3 =1.0$, giving an incompressible fluid,

**∇**⋅

**v**= 0, consistent with the theoretical calculations. A weak concentration gradient is imposed by enforcing Dirichlet (reservoir

^{11}) boundary conditions for the mass fractions at the top and bottom boundaries, $w y = 0 , \u2009 t = 0 . 2493 , 0 . 245 , 0 . 5057 $ and $w y = L , \u2009 t = 0 . 250 \u2009 729 , 0 . 255 , 0 . 494 \u2009 271 $. Thesevalues are chosen so that the deterministic diffusive flux of the first species vanishes at

*y*=

*L*/2,$ F \u0304 1 y = L / 2 , t =0$, leading to a diffusion barrier for the first species, as in Ref. 9.

The computational grid has 128 × 64 grid cells with grid spacing Δ*x* = Δ*y* = 1, with the thickness in the third direction set to give Δ*V* = 10^{6}, and time step Δ*t* = 0.1. In order to study the spectrum of the giant concentration fluctuations, we compute the Fourier spectrum of the mass fractions averaged along the direction of the gradient; this corresponds to *k _{y}* = 0 and thus

*k*=

*k*. A total of 1.8 × 10

_{x}^{5}time steps are preformed at the beginning of the simulation to allow the system to equilibrate. Statistics are then collected for 5 × 10

^{5}steps.

In this nonequilibrium example, the coupling to the velocity equation is crucial and is the cause of the giant fluctuations. The normal component of the velocity at the two physical boundaries follows from the EOS and the diffusive fluxes through the boundary, see Eq. (15) in Ref. 11. For the tangential (*x*) component of velocity, we use either no-slip (zero velocity) or free-slip (zero shear stress) boundary conditions. In the limit of infinite Schmidt number, *ν* = *η*/*ρ* ≫ χ, where χ is a typical mass diffusion coefficient, the overdamped equations apply and $\nu \u2009 S w ( i , j ) k $ approaches a limit independent of the actual value of the Schmidt number. For finite Schmidt numbers, however, the actual value of the Schmidt number affects the spectrum, see Appendix B in Ref. 9 for the explicit formulas.

For the overdamped integrator, the actual value of the viscosity does not matter beyond simply rescaling the amplitude of the fluctuations, since the velocity equation is a time-independent (steady) Stokes equation in the viscous-dominated limit. In the left panel of Fig. 3, we show numerical results for $\nu \u2009 S w ( i , j ) k $ for *i* ≠ *j*, obtained using our *overdamped* algorithm.^{12} Excellent agreement with the theoretical prediction in Appendix B in Ref. 9 is seen for wavenumbers larger than *L*^{−1}; for small wavenumbers, the confinement suppresses the giant fluctuations in a manner that depends on the specific boundary conditions imposed.^{10} In the right panel of Fig. 3, we show $\nu \u2009 S w ( 2 , 2 ) k $ for several values of the kinematic viscosity, as obtained using our *inertial* algorithm.^{12} The implicit-midpoint (Crank-Nicolson) scheme used to treat viscosity in the inertial algorithm is unconditionally stable and allows an arbitrary time step size to be used. It is, however, well-known that this kind of scheme can produce unphysical results for very large viscous Courant numbers due to fact it is not *L*-stable (see discussion in Appendix B in Ref. 40). This is seen in the results in the right panel of Fig. 3 for the largest viscosity *ν* = 1000 (corresponding to viscous Courant number *ν*Δ*t*/Δ*x*^{2} = 100) at larger wavenumbers. It is actually quite remarkable that we can use the inertial integrator with rather large time step sizes and get very good results over most of the wavenumbers of interest; this is a property that stems from a specific fluctuation-dissipation balance in the implicit midpoint scheme.^{40} These results demonstrate that both our overdamped and inertial methods are able to reproduce the correct spectrum of the nonequilibrium concentration fluctuations.

#### 3. Thermodiffusion and barodiffusion

In the giant fluctuation example shown in Fig. 3, the system was kept out of equilibrium by imposed concentrations on the boundaries, which is difficult to realize in experiments. Instead, experiments that measure giant fluctuations in liquid mixtures typically rely on the Soret effect to induce a concentration gradient via an imposed temperature gradient.^{42} A concentration gradient can also be induced via barodiffusion in the presence of large gravitational accelerations, as used in ultracentrifuges for the purposes of separation of macromolecules and isotopes.^{14} Barodiffusion and thermodiffusion enter in density equations (30) in the same manner; however, the key difference is that barodiffusion requires gravity which also enters via the buoyancy term in velocity Equations (3) and (5). Furthermore, the steady state gradient induced by barodiffusion is determined by equilibrium thermodynamics only and does not involve any kinetic transport coefficients.

It is well known that there is no nonequilibrium enhancement of the fluctuations at steady state for a system in a gravitational field^{43} in the absence of external forcing (see Eq. (28) in Ref. 44). (Note, however, that the dynamics of the fluctuations is affected by gravity and by barodiffusion.^{44}) This is because the system is still in thermodynamic equilibrium, despite the presence of spatial nonuniformity (sedimentation). In particular, without doing any calculations, we know that the equilibrium distribution of the fluctuations is the Gibbs-Boltzmann distribution, with a local free-energy functional that now includes a gravitational energy contribution. In this section, we demonstrate that our low Mach number approach captures this important distinction between (ordinary) equilibrium fluctuations in the presence of barodiffusion and (giant) nonequilibrium fluctuations in the presence of thermodiffusion.

We consider a solution of potassium salt and sucrose in water (see Sec. IV for more details) in an ultracentrifuge. The physical parameters of this ternary mixture are given in Sec. IV A and a brief theoretical analysis is given in Appendix B. We perform two dimensional simulations of a system of physical dimensions 0.8 × 0.8 × 0.1 cm divided into 64 × 64 × 1 finite-volume cells. Periodic boundary conditions are used in the *x* direction and impermeable no-slip boundaries are used in the *y* direction. The average mass fractions over the domain are set to $ w av = 0 . 0492 , \u2009 0 . 0229 , \u2009 0 . 9279 $. A total of 0.5 ⋅ 10^{6} time steps are performed at the beginning of the simulation to allow the system to equilibrate before statistics are collected for 10^{6} steps.

In order to induce a strong sedimentation in this mixture, we need to increase the ratio $ m 2 g/ k B T $ by six orders of magnitude relative to its reference value on earth (see (B2)). In actual experiments, this would be accomplished by increasing the effective gravity (i.e., centrifugal acceleration) in an ultracentrifuge; however, increasing gravity by such a large factor makes the system of equations (3), (5), and (30) numerically too stiff for our semi-implicit temporal integrator. This is because buoyancy changes the time scale for relaxation of large-scale (small wavenumber) concentration fluctuations from the usual slow diffusive relaxation to a very fast non-diffusive relaxation.^{45} Therefore, instead of increasing *g*, we artificially decrease *k _{B}* by six orders of magnitude and apply earth gravity

*g*= − 981 along the negative

*y*direction. With these parameters, our inertial temporal integrator is stable with time step size up to about Δ

*t*= 0.5 s; the results reported below are for Δ

*t*= 0.25 s.

For comparison, we use the same parameters but turn gravity off and induce a concentration gradient via thermodiffusion. Specifically, we set the temperature at the bottom wall to 293 K and 300 K at the top wall, and set the thermodiffusion constants to the artificial values $ D ( T ) = \u2212 5 \u22c5 1 0 \u2212 4 , \u2212 2 \u22c5 1 0 \u2212 4 , 7 \u22c5 1 0 \u2212 4 $. These values ensure that the steady state vertical profiles of the mass fraction of salt and sugar are very similar between the barodiffusion and thermodiffusion simulations (see (B3)), as shown in the left panel of Fig. 4. In the right panel of the figure, we show the spectrum $ S \rho k $ of the vertically averaged density. For comparison, in addition to the cases of gradients induced by baro and thermodiffusion, we also show the spectrum for a spatially uniform system at thermodynamic equilibrium, in the absence of gravity and at a constant temperature of 293 K. We see that the spectrum for barodiffusion is very similar to that for uniform thermodynamic equilibrium, while the spectrum for thermodiffusion shows the *k*^{−4} power-law behavior as in Fig. 3. This demonstrates that our numerical code correctly reproduces equilibrium fluctuations even in the presence of strong sedimentation.

## IV. DIFFUSION-DRIVEN GRAVITATIONAL INSTABILITIES

In this section, we use our numerical methods to study the development of diffusion-driven gravitational instabilities in ternary mixtures. In Section IV.D in Ref. 9, we studied the development of a diffusion-driven Rayleigh-Taylor (RT) instability in a ternary gas mixture. Here, we simulate similar instabilities in a ternary mixture of incompressible liquids, using realistic parameters corresponding to recent experimental measurements, and study the effect of (nonequilibrium) thermal fluctuations on the development of the instabilities. Our investigations are inspired by the body of work by Anne De Wit and collaborators on diffusion- and buoyancy-driven instabilities.^{4–6} In particular, in Ref. 5, a classification of these instabilities in a ternary mixture are proposed, and several of the instabilities are investigated experimentally. In the first part of this section, we perform simulations of the experimental measurements of a mixed-mode instability (MMI). In the second part, we investigate diffusive layer convection (DLC) (just as we did in Section IV.D in Ref. 9 for gases), in a hypothetical shadowgraphy or light scattering experiment that could, in principle, be performed in the laboratory.

We begin with a brief summary of the experimental setup of Carballido-Landeira *et al.*^{5} for getting a diffusion-driven gravitational instability in a simple ternary mixture: a solution of salt in water on top of a solution of sugar in water. The concentrations of salt and sugar are small so that even though this is a ternary mixture in this very dilute limit, one can think of salt and sugar diffusing in water without significant interaction. The key here is the difference in the diffusion coefficient between sugar (a larger organic molecule diffusing slower) and salt (a smaller ion diffusing faster) in water. Both sugar and salt solutions have a density that grows with the concentration of the solute.

In the experiments, one starts with an almost (to within experimental controls) flat and almost sharp interface between the two solutions. Even if one starts in a stable configuration, with the denser solution on the bottom, the differential diffusion effects can create a local minimum in density below the contact line and a local maximum above the contact line. This leads to an unstable configuration and the development of DLC at symmetric distances above and below the contact line. If one starts with an unstable configuration of the denser solution on top, before the RT instability has time to develop and perturb the interface, differential diffusion effects can lead to the development of local extrema in the density above and below the contact line, which are outside the range of the initial densities. The dynamics is then a combination of RT and DLC giving rise to a MMI. The DLC leads to characteristic “Y shaped” convective structures developing around the interface at the locations of the local adverse density gradients, which evolve around an interface that is slowly perturbed by the RT growth to a finite amplitude modulation. See Sec. III in Ref. 5 for more details, and the bottom row of panels in Fig. 1 in Ref. 5, as well as our numerical results in Fig. 6, for an illustration of the development of the instability.

### A. Physical parameters

We use CGS units in what follows (centimeters for length, seconds for time, grams for mass). Following the experiments of Carballido-Landeira *et al.*, we consider a ternary mixture of potassium salt (KCl, species 1, molar mass *M*_{1} = 74.55, denoted by A in Ref. 5), sugar (sucrose, species 2, molar mass *M*_{2} = 342.3, denoted by B in Ref. 5), and water (species 3, molar mass *M*_{3} = 18.02), giving molecular masses ** m** = (1.238 ⋅ 10

^{−22}, 5.684 ⋅ 10

^{−22}, 2.99 ⋅ 10

^{−23}). The initial configuration is salt solution on top of sugar solution. In Ref. 5, it is assumed that the density dependence on the concentration can be captured by (this is a good approximation for very dilute solutions)

where $ \rho 0 = \rho \u0304 3 =1.0$ is the density of water, *α*_{1} = 48 for KCl, and *α*_{2} = 122 for sucrose and *Z _{k}* is the molar density of each component, related to the partial density via

*ρ*=

_{k}*Z*, where

_{k}M_{k}*M*is the molar mass. Noting that we can write our EOS (2) in the form

_{k}which gives us the EOS parameters

Note that here these should not be thought of as pure component densities since the solubility of the solvents in water is finite, rather, they are simply parameters that enable us to match our model EOS (2) to the empirical density dependence in the dilute regime.

For the dilute solutions we consider here, it is sufficient to assume that ** D** is constant for the range of compositions of interest and the mixtures are essentially ideal,

**=**

*H***0**. We also assume isothermal conditions with the ambient temperature to

*T*= 293 K, and assume constant viscosity

*η*= 0.01. Since the ternary mixture under consideration is what can be considered “infinite dilution,” we rely on the approximation proposed in Ref. 27,

where from Table I in Ref. 5, we read the diffusion coefficient of low-dilution KCl in water as *D*_{1} = 1.91 ⋅ 10^{−5} and the diffusion coefficient of dilute sucrose in water as *D*_{2} = 0.52 ⋅ 10^{−5}. Here, *D*_{3} is the self-diffusion coefficient of pure water, *D*_{3} = 2.3 ⋅ 10^{−5}, giving *D*_{12} ≈ 4.32 ⋅ 10^{−6}.

In the simulations reported below, we have neglected barodiffusion and assumed that the pressure is equal to the atmospheric pressure throughout the system. If barodiffusion were included, we would obtain partial sedimentation in the final state of the mixing experiment (where the system has completely mixed), as discussed in more detail in Sec. III C 3 and Appendix B. In the final steady state, (B2) suggests that the ratio of the mass fraction of sugar between the top and bottom of a cell of height *H* = 1 will be

which is indeed negligible and likely not experimentally measurable.

### B. Mixed-mode instability in a Hele-Shaw cell

In this section, we examine the mixed-mode instability illustrated in the bottom row of panels in Fig. 1 in Ref. 5. The geometry of the system is a Hele-Shaw cell, i.e., two parallel glass plates separated by a narrow gap, as illustrated in the left panel of Fig. 5. In our simulation, we model a domain of length 0.8 × 0.8 × 0.025, with gravity *g* = − 981 along the negative *y* direction. It is well-known that the flow averaged along the *z* axes in a Hele-Shaw setup can be *approximated* with a two-dimensional Darcy law; this is used in the simulations reported in Ref. 5. Here, we do not rely on any approximations but rather model the actual three-dimensional structure of the flow in all directions. We divide our domain into 256 × 256 × 8 cells and impose periodic boundary conditions in the *x* direction and no-slip walls in the *z* direction. In the y direction, we impose a reservoir (Dirichlet) boundary conditions for concentration to match the initial concentrations in the top and bottom half of the domain and impose a free-slip boundary condition on velocity to model an open reservoir of salt solution on the top and sugar solution at the bottom of the domain.

Note that the momentum diffusion time across a domain of length *L* ∼ 1 is *τ _{L}* ∼

*L*

^{2}/(2

*ν*) ≈ 50. The experimental snapshots in Fig. 1 of Ref. 5 indicate that this is comparable to the time it takes for the instability to fully develop. It is therefore not safe to rely on the overdamped approximation, so we use the inertial integrator described in Ref. 12 in our simulations. Nevertheless, the overdamped limit is a relatively good approximation in practice since the characteristic length scale of the

*Y*-shaped DLC fingers observed in the experiments is

*L*∼ 0.013, corresponding to viscous diffusion time

*τ*∼ 8.5 ⋅ 10

_{L}^{−3}. We have compared numerical simulations using the overdamped and inertial integrators and found little difference. We simulate the development of the instability to

*t*≈ 63.9 with a fixed time step size of Δ

*t*= 6.39 ⋅ 10

^{−2}, which corresponds to 75% of the stability limit dictated by our explicit treatment of mass diffusion. We use a centered discretization of advection; in these well-resolved simulations, little difference is observed between centered advection and the more sophisticated BDS advection scheme

^{34}summarized in Ref. 12.

The initial concentrations, denoted with superscript zero in what follows, on the top and bottom are determined from the dimensionless ratio reported in Ref. 5,

Specifically, we set the initial mass fractions of salt and sugar to $ w 1 0 =0.0864$ and $ w 2 0 =0.1368$, respectively, more precisely,

which, from the EOS (2) gives a density difference of about 0.8%. Similar results are observed for other values of the concentrations if the dimensionless ratio *R* in (38) is kept fixed, with the main difference being that lower concentrations lead to a slower development and growth of the instability.

In the experiments, the initial interface between the two solutions is not, of course, perfectly flat. To model this effect, it is common in the literature to add a small random perturbation to the initial conditions. Assuming that the growth of the unstable modes is exponential in time, the time it will take for the instability to reach a certain point in its development (e.g., to first split the Y-shaped fingers) will depend on the amplitude of the initial perturbation. Since the initial condition is not known to us and is impossible to measure experimentally to high accuracy, it is not possible to directly compare snapshots of the instability at “the same point in time” between simulations and experiments, or between simulations that use different initial perturbations. Instead, we compare a simulation in which thermal fluctuations are accounted for with one in which thermal fluctuations are not accounted for, both starting from the *same* randomly perturbed initial interface.

Specifically, we randomly perturb the concentrations in the layer of cells just above the interface, setting the mass fractions to $w=r w bottom 0 + ( 1 \u2212 r ) w top 0 $ next to the interface, with *r* being a uniformly distributed random number between 0 and 0.1. As we explained above, a direct comparison between simulations and experiments is not possible. Nevertheless, our numerical results shown in Fig. 6 are very similar, in both time development and visual appearance, to the experimental images shown in the bottom row of panels in Fig. 1 in Ref. 5. Here, we use density *ρ* as an indicator of the instability even though in actual light scattering or shadowgraph experiments it is the index of refraction rather than density that is observed; both density and index of refraction are (approximately) linear combinations of the concentrations and should behave similarly. More detailed comparisons between simulations and experiments require a careful coordination of the two and will not be attempted here. Instead, we focus our attention on examining the role (if any) played by the thermal fluctuations in the triggering and evolution of the instability.

Comparing the left and right column of panels in Fig. 6 shows only minor differences between the deterministic and the fluctuating hydrodynamics calculation. This indicates that the dynamics is dominated by the unstable growth of the initial perturbation, with thermal fluctuations adding a weak perturbation, which, though weak, does to some extent affect the details of the patterns formed at late times but not their generic features. By tuning the strength of the initial perturbation, one can also tune the impact of the fluctuations on the dynamics; after all, if one starts with a perfectly flat interface, the instability will be triggered by thermal fluctuations only. Nevertheless, we can conclude that once the instability develops sufficiently (e.g., the Y-shaped fingers form), the dynamics becomes dominated by the deterministic growth of the unstable modes. This can be seen from the fact that stochastic simulations (not shown) starting from a perfectly flat interface develop the same features at long times as deterministic simulations with a random initial perturbation.

To get a more quantitative understanding of the development of the instability, in Fig. 7, we show the Fourier spectrum of the density $ \rho \u0304 x $ averaged along the *y* and *z* directions, which is a measure of the fluctuations of the diffusive “interface.” In these investigations, we used a simulation box of size 1.6 × 0.2 × 0.05 (grid size 512 × 64 × 16 cells), and set the initial mass fractions to be four times smaller (note that this does not change the dimensionless number *R* in (38)), $ w 1 0 =0.0216$ and $ w 2 0 =0.0342$, in order to slow down the development of the instability and allow (giant) nonequilibrium thermal fluctuations to develop. The remaining parameters were identical to those reported above.

For a gravitationally stable configuration, the spectrum of $ \rho \u0304 x $, called the static structure factor *S*(*k*), exhibits a characteristic giant fluctuation *k*^{−4} power-law decay at large wavenumbers *k* at long times,^{1} as illustrated in Fig. 3 by showing *S*(*k*) for a simulation in which gravity is switched off. At small wavenumbers, gravity and confinement damp the amplitude and affect the dynamics of the giant fluctuations, as is well-understood for stable steady states.^{1,10,44–46} We emphasize that the thermal fluctuations in these investigations are completely dominated by the nonequilibrium (giant) fluctuations induced by the presence of large concentration gradients at the interface; the equilibrium concentration fluctuations are many orders of magnitude weaker in comparison. This is seen by the fact that we observe identical results if we completely turn off the stochastic mass fluxes and keep *only* the stochastic momentum flux (velocity fluctuations). Similarly, turning off the stochastic momentum flux but keeping the stochastic mass flux, in this particular setup, leads to very small fluctuations that are below the tolerance of our linear solvers.

To our knowledge, nonequilibrium fluctuations have not been examined in an unstable situation like the one we study here; some studies have been carried out near the onset of the Rayleigh-Bernard thermal convection instability.^{47,48} It is important to emphasize that the very validity of fluctuating hydrodynamics has *not* been established for unstable flows; one can justify linearized fluctuating hydrodynamics,^{1,49,50} but in a linearization, the fluctuations do *not* affect the deterministic flow and thus cannot trigger the instability. The results in Fig. 7 show that in the absence of thermal fluctuations, there is a clear band of wavenumbers *k* ≲ 64 whose amplitude grows in time due to the instability, while the remaining modes decay to zero rapidly. It is important to note that the most unstable wavelength of λ ≈ 2*π*/45 ≈ 0.14 cm is very close to the experimentally reported value of 0.13 cm.^{5} When thermal fluctuations are present, the same band of wavenumbers grows in a similar manner, but at the same time, larger wavenumbers show a nontrivial power-law spectrum. Eventually, however, the unstable modes completely dominate the dynamics and there is essentially no difference between the simulations with and without thermal fluctuations; to make a more precise quantitative statement, multiple Monte Carlo simulations are required to perform ensemble averages and reduce the statistical error, as well as to eliminate the effects of the boundary conditions along the horizontal and vertical directions. Note that at later times, the fluid flow becomes chaotic and further information may be gained by also examining the spectrum of the fluid velocity and not just the vertical projection of the density; we leave such detailed investigations for future studies.

### C. DLC instability

In this section, we use our algorithm to model a hypothetical experiment in which the gravity points in the direction perpendicular to the fluid-fluid interface and to the confining glass slips, as illustrated in the right panel of Fig. 5. This kind of geometry has already been used to experimentally investigate isothermal free diffusive mixing in a *binary* liquid mixture^{45} in a *stable* configuration. (Note that for a binary mixture, the only alternative to a stable configuration is to have a RT-unstable configuration.) To our knowledge, similar experiments have not been performed for ternary mixtures, although it is feasible they could be if the initial configuration can be prepared with sufficient control.

To be specific, we select the initial concentrations of the sugar and salt solutions to get a DLC instability; all other parameters are the same as for the MMI setup summarized above. The dimensionless parameter used for Fig. 1(c) in Ref. 5 is

giving the required ratio of mass fractions,

We set $ w 1 0 =0.022$ and $ w 2 0 =0.05$, which give an initial density of 1.014 375 on the top and 1.018 062 on the bottom, for a density difference of about −0.4% (recall that the sign of the difference is opposite for DLC and MMI configurations).

The dimensions of the domain in our simulations are 1 × 0.5 × 1, with periodic boundary conditions in *x* and *z*, and reservoir (Dirichlet) boundary condition for concentration together with free-slip boundary condition on velocity in the *y* direction. The grid is 256 × 128 × 256 grid cells and the time step size is fixed at Δ*t* = 0.025, which corresponds to a maximum advective Courant number of ∼0.5. We treat advection using the unlimited bilinear BDS high-resolution scheme^{51} summarized in Ref. 12. The initial interface is now perfectly flat, so that the instability is triggered by the nonequilibrium (giant) thermal fluctuations. We focus our investigation here on the differences in the spectrum of the thermal fluctuations with (unstable) and without (marginally stable) gravity. Note that experimental measurements of giant concentration fluctuations in microgravity are feasible and have been performed for binary liquid mixtures;^{8} microgravity measurements on ternary mixtures are ongoing or in the planning stage.^{3}

In this geometry, the interface between the two fluids cannot be visualized, rather, one sees the average index of refraction (a linear combination of the average concentrations) along the thickness of the sample (direction of the gradient and of gravity). In Fig. 8, we show color plots for the vertically averaged density at two points in time, one as the instability is just beginning to dominate the dynamics and the other as the instability has fully developed. In the presence of gravity, at early times, there are pronounced giant thermal fluctuations in addition to the fluctuations coming from the growth of the unstable modes, which dominate the dynamics at late times similarly to the MMI case in Fig. 6. The giant nonequilibrium fluctuations are more clearly visualized by considering stable free diffusive mixing in microgravity (*g* = 0, bottom panels).

In Fig. 9, we show the radially averaged power spectrum $S ( k = k ) $ corresponding to the vertically averaged density $ \rho \u0304 x , z $. Behavior similar to the MMI instability shown in Fig. 6 is observed. The band of wavenumbers between 16 ≲ *k* ≲ 128 is seen to grow in time, and the characteristic *k*^{−4} giant fluctuation spectrum is seen at large wavenumbers, especially clear in the results for microgravity. Deterministic simulations (not shown) starting from a randomly perturbed interface show similar unstable growth and are eventually indistinguishable from simulations in which the instability is triggered and enhanced by thermal fluctuations.

We defer detailed studies of these phenomena for future work, in the hope that our work will stimulate careful experimental investigations that can directly be compared to our computer simulations. By adjusting the initial concentrations and the choice of the two solutes, one can change the range of unstable wavenumbers and the time scale for the development of the instability; this is best done by performing a linear stability analysis.^{4} It is feasible that for some choice of parameters, a more substantial interaction between the spectrum at large wavenumbers, dominated by giant fluctuations, and the spectrum at small wavenumbers, dominated by the instability, will occur. If this is the case, one may be able to observe a measurable influence of nonequilibrium concentration fluctuations on the development and grown of the instability, allowing us, for the first time, to experimentally access the validity of nonequilibrium *nonlinear* fluctuating hydrodynamics.

## V. CONCLUSIONS

We have developed a low Mach number fluctuating hydrodynamics formulation of momentum and mass transport in non-ideal mixtures of incompressible liquids with given pure-component densities. In the present low Mach number model, energy transport is not modeled explicitly and the temperature is assumed constant in time. Momentum transport is accounted for in a quasi-incompressible framework in which pressure is in mechanical (hydrostatic) equilibrium and fast sound waves are eliminated from the system of equations because density instantaneously follows the local composition. The thermodynamics of the mixture is described by the Hessian of the normalized excess Gibbs energy per particle. The transport properties are given by the shear viscosity, the Maxwell-Stefan diffusion coefficients, and the thermal diffusion coefficients or ratios, as a function of the composition of the mixture.

We used our low Mach number algorithm to study the development of gravitational instabilities during diffusive mixing in ternary mixtures. These mixed mode and diffusive convection instabilities are specific to mixtures of more than two species and occur because of an interaction between the familiar buoyancy-driven Rayleigh-Taylor instability with differential diffusion effects. Our simulations of the mixed-mode instability closely mimic recent experiments^{5} performed in a Hele-Shaw setup, and we found good qualitative agreement between experimental and computational results. A more quantitative comparison is not possible at this stage and should be the subject of future work. A particularly challenging aspect is to experimentally control or measure the initial conditions with sufficient accuracy to be able to directly compare simulations to experiments.

Because of the presence of sharp gradients at the interface between the two diffusively mixing solutions, giant concentration fluctuations develop in the form of power-law tails in the spectrum of the concentration fluctuations. These nonequilibrium fluctuations are much larger than equilibrium ones and have the potential to trigger and feed the instability and affect the growth of the unstable structures. It is important to observe that such a coupling of the fluctuations back to the macroscopic dynamics requires nonlinearity, and is not possible in linearized fluctuating hydrodynamics. While some nonlinear effects in fluctuating hydrodynamics have been verified to occur in real liquids, it is not clear whether our nonlinear fluctuating hydrodynamic simulations can account for the effect of the fluctuations on the evolution of the mean flow. In particular, the nonlinear fluctuating hydrodynamic equations are ill-posed and some regularization, implicit in our finite-volume discretization, is required to even give meaning to the equations.^{52} The conclusions of the simulations performed here is that fluctuations, despite being “giant,” are quickly overwhelmed by the deterministic growth of the unstable modes. This suggests that in actual experiments, the development of the instability is primarily triggered by the imperfections in the initial condition or external fluctuations (e.g., vibrations).

In the Hele-Shaw geometry studied experimentally in Ref. 5, the no-slip boundaries strongly damp the velocity fluctuations and reduce the giant fluctuations. Furthermore, the quasi-two dimensional geometry limits the possibilities for interactions between the fluctuations and the instability. For this reason, we also reported simulations of diffusive layer convection for a three-dimensional geometry, as used in existing experiments in binary mixtures. Because giant fluctuations develop on a slow diffusive time scale for large and intermediate wavenumbers, it takes some time after the initial contact between the two fluids for the power-law spectrum to develop. At the same time, the unstable smaller wavenumbers have an exponentially growing amplitude, so that the instability can develop much faster than the giant fluctuations. In our simulations, we reduced the concentrations of solutes to make the density difference very small and thus slow down the unstable growth; it remains to be seen what is possible and can be observed experimentally.

An important challenge for experimental physicists is to devise ways to measure the thermophysical properties in multispecies mixtures.^{3} While our formulation can handle mixtures of arbitrary numbers of species with essentially complete physical fidelity, it is not possible to use our codes for realistic non-dilute mixtures because many of the parameters, notably the thermodynamic factors and the Maxwell-Stefan diffusion coefficients, are missing. We believe that it will be necessary to use computer simulations using the types of methods described here, together with Monte Carlo sampling of parameters, and potentially also molecular dynamics calculations,^{37,53} *in addition* to experimental measurement of observation functions, in order to obtain robust estimates of thermophysical parameters with a quantified uncertainty. The traditional approach of tabulating values with error bars, that has worked for binary mixtures, fails for multispecies mixtures due to proliferation of a larger number of parameters that are not all independent but are constrained by thermodynamic symmetries and stability conditions.

In future work, we will consider the inclusion of chemical reactions in our low Mach number models. It is also possible to include thermal effects in our formulation by accounting for the effects of thermal expansion in the quasi-incompressibility constraint. Two key difficulties are constructing a spatial discretization that ensures preservation of an appropriately generalized EOS, as well as developing suitable temporal integrators that can handle the multitude of time scales that appear due the presence of slow mass, fast momentum, and intermediate energy diffusion, and, potentially, ultrafast chemical reactions. Another challenge for future work on low Mach number fluctuating hydrodynamics is to account for the effects of surface tension in mixtures of immiscible or partially miscible liquids. While some thermodynamically consistent constructions of diffuse-interface models exist for ternary mixtures,^{24} we believe that much more work is required to formulate a complete multispecies system of equations in the presence of square-gradient terms in the free energy functional, even in the isothermal setting.

## Acknowledgments

We would like to thank Anne De Wit and Jorge Carballido Landeira for numerous discussions regarding their experiments on gravitational instabilities in ternary mixtures. This material is based upon work supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award No. DE-SC0008271, and under Contract No. DE-AC02-05CH11231.

### APPENDIX A: COMPUTING THE DIFFUSION MATRIX

An essential complication with treating all species equally is that the *N* mass or mole fractions are not all independent but must sum to unity. For example, Eq. (23) are not independent, and one must supplement them with the condition that the total mass flux is zero. Similarly, the rows and columns of **Λ** sum to zero, **Λ1** = **0**, where **1** denotes a vector of ones, and therefore, **Λ** is not invertible; the range of **Λ** are vectors that sum to zero. To deal with these complications, let us now introduce the projection matrix^{18}

which satisfies *Q*^{T}**w** = ** Q1** =

**0**; therefore, pre-multiplying a matrix by

*Q*^{T}ensures that the matrix has a range consisting of vectors that sum to zero, and post-multiplying by

*Q*^{T}ensures that the matrix has

**w**in its null-space. Now, let us define a diffusion matrix as

where *α* ≠ 0 is an arbitrary constant, for example, $\alpha =Trace \Lambda $; an alternative equivalent formula is (24).

Here is a quick summary of some useful relations derived in Ref. 54.

**Λ1**=**0**and**Λ**is SPD on**w**^{⊥}.**χw**=**0**and**χ**is SPD on**1**^{⊥}.**χ**and**Λ**are generalized inverses of each other,**χΛχ**=**χ**and**ΛχΛ**=**Λ**, more specifically,**Λχ**=*Q*^{T}and**χΛ**=.*Q*- The SPD matricesare inverses of each other, and $ \Lambda \u0303 \u2261\Lambda $ on the subspace$ \Lambda \u0303 = \Lambda + \alpha w w T = \chi \u0303 \u2212 1 = \chi + \alpha \u2212 1 1 1 T \u2212 1 $
**1**^{⊥}and $ \chi \u0303 =\chi $ on the subspace**w**^{⊥}.

In order to numerically compute **χ** from **Λ**, one could use a matrix inverse, as in (24). This fails, however, when some of the species vanish since the matrix **Λ** + *α***ww**^{T} is not invertible; this can be fixed by using a pseudoinverse, computed via the singular value decomposition (SVD) of **Λ**. However, both matrix inversion and SVD involve *O*(*N*^{3}) operations and can be expensive for large numbers of species. An alternative is to use an iterative numerical procedure^{54}

where ** N** =

**−**

*I*

*M*^{−1}

**Λ**and

**is a diagonal matrix with entries**

*M* The sum converges rapidly so only a few terms in the sum are needed to compute **χ**_{M} without having to do matrix inverses, but the speed of convergence is hard to access *a priori* and in practice, we set *M* to a fixed integer such as *M* = 5 or *M* = 10. It is important to note that, as proven in Ref. 18, the truncated sum to *M* terms gives an approximation **χ**_{M} ≈ **χ** that is symmetric positive semi-definite, and satisfies the properties **χ**_{M}**w** = **0** as it must. It is therefore perfectly consistent with thermodynamics to just approximate **χ** with **χ**_{M}. After all, since the Maxwell-Stefan diffusion coefficients are only known to at most two decimal places in practice, it makes little sense to invert **Λ** exactly or compute its (expensive) SVD. Instead, **χ**_{M} for small *M* is in practice an equally good approximation to the true diffusion matrix.

### APPENDIX B: STATIC EQUILIBRIUM OF A DILUTE SOLUTION IN GRAVITY

Barodiffusion is often neglected for liquid mixtures since its contribution to the diffusive fluxes is typically negligible compared to those due to compositional or temperature gradients, unless there are very large pressure gradients due to large accelerations, as in ultracentrifuges. Nevertheless, including barodiffusion is crucial in order to obtain the correct equilibrium steady state of a mixture in the presence of gravity. This is a direct consequence of the fact that barodiffusion has thermodynamic origin. If barodiffusion is neglected, at thermodynamic equilibrium, the mixture would mix uniformly even in the presence of gravity, reaching the state of maximal entropy. However, physically, we know that heavier species will migrate toward the bottom, minimizing the free energy by balancing the gain in potential energy with the loss of entropy. Here, we show that once barodiffusion is accounted for the hydrodynamic equations correctly reproduce the statistical mechanics of systems in a gravitational field.

As a simple case in which we know the “correct” answer, let us focus on the case when the mixture is a dilute suspension of a number of solutes such as colloids or a macromolecule (e.g., solution of sugar in water). In the dilute limit, the different solutes do not affect each other, and we can, in fact, focus on one of the solutes only and consider a binary solution. Let us take the first species to be the solute and the second species to be the solvent. For a dilute solution, *w*_{1} ≪ 1, **Γ** ≈ ** I**, $ m \u0304 \u2248 m 2 $, $\rho \u2248 \rho 2 \u2248 \rho \u0304 2 $, giving $ x 1 \u2248 m 2 / m 1 w 1 $ and $ \varphi 1 =w\rho / \rho \u0304 1 \u2248w \rho \u0304 2 / \rho \u0304 1 $. Note that hydrostatic equilibrium is built into the low Mach number equations through the reference pressure, which satisfies

*dP*/

*dh*=

*ρg*, where

*h*is the height. Thermodynamic equilibrium, i.e., a vanishing of the chemical potential gradients, corresponds to a vanishing of the diffusion driving force for the solute,

which directly gives the gravitational sedimentation profile of an ideal gas with molecular mass *m _{e}*,

where *m _{e}* is the excess mass of the colloids over that of the fluid,

The result (B2) is in agreement with the statistical-mechanical notion that the solute subsystem can be thought of as an ideal gas of particles subject to the Archimedean gravity force *m _{e}g*.

A similar calculation for the case of thermodiffusion gives

which is similar in form but does not have an origin in equilibrium statistical mechanics.

### APPENDIX C: FLUCTUATING MAXWELL-STEFAN DESCRIPTION

In this Appendix, we show that thermal fluctuations can be directly added to the Maxwell-Stefan formulation in a manner that is intuitive and consistent with the one we presented here, derived from the principles of nonequilibrium thermodynamics.

In the Maxwell-Stefan description of diffusion, one considers a frictional force between species *i* and *j* proportional to the difference in the particular velocities of the two species, with friction coefficient *γ _{ij}* =

*x*/

_{i}x_{j}*D*. It is natural to include thermal fluctuations in this description by including a fluctuating component to the frictional force with covariance proportional to

_{ij}*γ*, following the traditional Langevin approach. This leads us to the MS equation with fluctuations,

_{ij}where $ Z \u02dc ij r , t =\u2212 Z \u02dc j i r , t $ are space-time white-noise processes, one process per pair of species, for a total of $N N \u2212 1 /2$ random forces. From the above, we may augment (20) to account for fluctuations as

where **𝒵** is a vector of $N N \u2212 1 /2$ independent white noise processes and ** K** is a matrix that can directly be read from (C1). Following the same linear algebra steps as before, we can solve for the fluxes to obtain the additional fluctuating contribution

The covariance of this stochastic flux is

It can be shown that $\chi K K \u22c6 \chi =\chi $, which shows that the covariance

is proportional to Onsager matrix (26), and therefore, the stochastic mass fluxes obtained from the fluctuating MS description are statistically identical to (29). This justifies the prefactor $ 2 / n $ in the noise in (C1); one can also justify this factor based on kinetic theory considerations.

In this fluctuating MS description, the stochastic mass flux is constructed explicitly without using Cholesky factorization of the Onsager matrix. That is, this construction gives a “square root” of the Onsager matrix as an explicit construct,

which is different from (28). Note, however, that this formulation requires using $N N \u2212 1 /2$ random processes instead of only *N* − 1 random processes required when (28) is used; this increase in the number of random numbers is typically not justified by the savings of a Cholesky factorization.

### APPENDIX D: STATIC AND DYNAMIC STRUCTURE FACTORS

In this section, we compute equilibrium static (32) and dynamic (31) structure factors for a mixture at equilibrium. In linearized fluctuating hydrodynamics, the coupling between the mass (concentration) equations and the velocity equation disappears at equilibrium, and we can focus our attention on mass equation (30). When (30) is linearized around a uniform state at thermodynamic equilibrium, we can omit the advective term $\u2207\u22c5 \rho w v $, and we can treat density as constant.

Since we want to compute fluctuations in the mass fractions, we first convert (30) to use gradients of mass fractions instead of gradients of mole fractions, by using the chain rule

The symmetric Hessian matrix

We can therefore write the mass flux due to composition gradients in terms of mass fraction gradients as

where

At thermodynamic equilibrium, the fluctuating diffusion equation (30) can be expressed entirely in terms of the fluctuations of mass fractions by using (D1),

where $N= 2 / n \u2009W \chi 1 2 $. After taking a Fourier transform (D2) becomes a multi-dimensional Ornstein-Uhlenbeck (OU) equation, one system of *N* equations for each wavenumber. The dynamic structure factor is

A standard result for OU processes^{41,55} states that at steady state, the covariances satisfy the linear system of equations

This equation needs to be supplemented by the conditions that the static structure factor $ S w k $ is symmetric, and that the row and column sums of ** S_{w}** are zero because the mass fractions sum to one,

This system of two equations for ** S_{w}** has a unique solution; after some algebraic manipulations, the solution can be written in evidently symmetric form (33). For ideal mixtures,

**= 0 and (33) can be simplified to**

*H* which matches the corresponding expression for a mixture of ideal gases.^{2}

In actual experiments, what can be measured through dynamic light scattering or shadowgraphy is the spectrum of refractive index $ n r w ; P 0 , T 0 $. A related quantity is the density structure factor $ S \rho = \delta \rho \u0302 \delta \rho \u0302 \u22c6 $, which can easily be obtained from ** S_{w}** in the low Mach setting as follows. First, observe that in the low Mach number limit, density fluctuations do not include the contribution from pressure fluctuations (sound peaks in the static structure factor), rather, they only include a contribution due to fluctuations in mass fractions (central peak).

^{11}Physically, this corresponds to observing density fluctuations at a longer time scale, i.e., averaging over the fast pressure fluctuations on the sonic time scale. If we expand EOS constraint (2) to account for small thermal fluctuations,

*ρ*←

*ρ*+

*δρ*and

*w*←

_{i}*w*+

_{i}*δ*

**w**

_{i}, we get

giving the density fluctuations

From this relation, we can derive all properties of the density fluctuations, such as static and dynamic structure factors, from the corresponding values for the mass fractions. For example, the dynamic or static structure factor of density can be obtained from the corresponding result for the mass fractions via (34).

## REFERENCES

_{2}− O

_{2}− N

_{2}flames

Note that it is straightfoward to modify the standard Cholesky factorization algorithm to work for semi-definite matrices by simply avoiding division by zero pivot entries; the factorization process remains numerically stable and works even when some of the species vanish.