# The multicomponent diffuse-interface model
and its application to water/air interfaces

E. S. Benilov111Email address for correspondence<EMAIL_ADDRESS>Department of Mathematics and Statistics, University of Limerick,
Limerick V94 T9PX, Ireland

###### Abstract

Fundamental properties of the multicomponent diffuse-interface model (DIM),
such as the maximum entropy principle and conservation laws, are used to
explore the basic interfacial dynamics and phase transitions in fluids. Flat
interfaces with monotonically-changing densities of the components are proved
to be stable. A liquid layer in contact with oversaturated but stable vapour
is shown to either fully evaporate or eternally expand (depending on the
initial perturbation), whereas a liquid in contact with saturated vapour
always evaporates. If vapour is bounded by a solid wall with a sufficiently
large contact angle, spontaneous condensation occurs in the vapour. The
external parameters of the multicomponent DIM – e.g., the Korteweg matrix
describing the long-range intermolecular forces – are determined for the
water–air combination. The Soret and Dufour effects are shown to be negligible
in this case, and the interfacial flow, close to isothermal.

##  1 Introduction

The diffuse-interface model (DIM) was proposed by Korteweg (1901) as an
attempt to avoid the abrupt change of parameters in the models of
liquid/vapour interfaces existing at the time. It is based on two assumption:

(i) The long-range attractive intermolecular force (the van der Waals force)
can be modelled by a _pair-wise_ potential, so that the force affecting a
molecule is the algebraic sum of those exerted on it by other molecules.

(ii) The characteristic distance over which the van der Waals force acts is
much smaller than the interfacial thickness.

The resulting representation of the molecular interaction has been
incorporated into models of various phenomena, such as phase transitions in
ferroelectric materials (Ginzburg, 1960), spinodal decomposition (Cahn, 1961;
Lowengrub & Truskinovsky, 1998), growth of, and oscillations in, crystals
(Collins & Levine, 1985; Tang et al., 2006; Heinonen et al., 2016; Philippe et
al., 2020), solidification of alloys (Stinner et al., 2004; Nestler et al.,
2005), phase separation in polymer blends (Thiele et al., 2007; Madruga &
Thiele, 2009), electrowetting (Lu et al., 2007), contact lines (Jacqmin, 2000;
Pismen & Pomeau, 2000; Ding & Spelt, 2007; Yue et al., 2010; Yue & Feng, 2011;
Sibley et al., 2014; Ding et al., 2017; Borcia et al., 2019), contact lines in
liquids with surfactants (Zhu et al., 2019, 2020), Faraday instability (Borcia
& Bestehorn, 2014; Bestehorn et al., 2021), Rayleigh–Taylor instability
(Zanella et al., 2020, 2021), cavitation (Petitpas et al., 2009), nucleation
and collapse of bubbles (Magaletti et al., 2015, 2016; Gallo et al., 2018,
2020), capillary condensation (Pomeau, 1986; Benilov, 2022a) liquid films
(Benilov, 2020c, 2022a), tumor growth (Frigeri et al., 2015; Rocca & Scala,
2016; Dai et al., 2017), classification of high-dimensional data (Bertozzi &
Flenner, 2012, 2016), _etc._

The present paper is mainly concerned with application of the DIM to
evaporation of drops and condensation of vapour on a solid. These phenomena
have been examined using the single-component version of the DIM (Benilov,
2022b, a), where the fluid/gas interface is modelled by that between the
liquid and vapour phases of the same fluid. The evaporation in this case was
shown to be due to a flow caused by a weak imbalance of chemical potentials of
the liquid and vapour phases.

It is unclear, however, how the results obtained via the single-component DIM
are modified by the effect of air, whose density (at, say,
$25^{\circ}\mathrm{C}$) exceeds that of water vapour by a factor of more than
$50$. Furthermore, there are three additional physical effects in
multicomponent fluids: diffusion (of water vapour in air), the Soret effect
(thermodiffusion), and the Dufour effect (heat flux due to density gradient).
Only one of these, the diffusion, has been examined before (e.g., Deegan et
al., 2000; Dunn et al., 2009; Eggers & Pismen, 2010; Colinet & Rednikov, 2011;
Rednikov & Colinet, 2013; Morris, 2014; Stauber et al., 2014, 2015; Janeček et
al., 2015; Saxton et al., 2016, 2017; Rednikov & Colinet, 2019; Wray et al.,
2019), but the models employed in these papers do not include the flow-induced
evaporative flux discovered via the DIM. Thus, the multicomponent diffuse-
interface model appears to be a tool describing all the mechanisms at work in
evaporation/condensation of liquids into/from air.

The same can hopefully be said about the dynamics of contact lines, as most of
the existing models work for some fluids (including water) only if the so-
called slip length – effectively, the interfacial thickness – is set to be
unrealistically small (Podgorski et al., 2001; Winkels et al., 2011;
Puthenveettil et al., 2013; Limat, 2014; Benilov & Benilov, 2015).

Before using a new model, one generally needs to examine its basic properties,
test it on problems with well-understood physics (to ensure that the
mathematics captures it), and parameterise this model for the intended
applications. These are the three aims of the present work in the context of
the multicomponent DIM.

The following results are reported:

  1. 1.

Multicomponent flat interfaces with monotonically-changing densities of the
species (components) are all stable. This conclusion follows from the entropy
principle and conservation laws.

  2. 2.

Several new physical effects are described, the most interesting of which is
evaporation of a flat liquid layer in contact with saturated vapour. This
phenomenon appears to be similar to evaporation of drops surrounded by
saturated vapour (Deegan et al., 2000; Eggers & Pismen, 2010; Benilov, 2020d,
2021, 2022b), but with one important difference: the drops evaporate because
the curvature of their boundary increases the effective saturation pressure
(the so-called Kelvin effect), making the saturated vapour effectively _under_
saturated and, thus, encouraging evaporation. This explanation is clearly
inapplicable to layers with flat boundaries – yet they evaporate anyway. The
actual mechanism is based on the long-range interaction of the liquid/vapour
interface with the substrate – which implies that, for macroscopic liquid
films, this effect is weak. It should be important, however, for nano-films
whose thickness is comparable to the interfacial thickness.

  3. 3.

The multicomponent DIM is parameterised for water/air interfaces at normal
conditions, which can be used in the future for modelling evaporation and
condensation of water in the Earth’s atmosphere or the dynamics of contact
lines of water drops.

The paper has the following structure. In section 2, the problem is formulated
mathematically. Section 3 examines the entropy principle and conservation
laws. In section 4, the governing equations are nondimensionalised and the
main nondimensional parameters are identified. Sections 5–7 examine basic
solutions of the DIM, and in section 8, the DIM is parameterised for water/air
interfaces. Section 9 summarises the results obtained, including the effects
partly mentioned in item 2 of the above list.

##  2 Formulation

###  2.1 Thermodynamics

When studying hydrodynamics of a compressible fluid, one has to deal with its
thermodynamic properties. In this section, they are described briefly and in a
self-contained form (for the benefit of readers specializing is incompressible
hydrodynamics).

Consider an $N$-component compressible non-ideal fluid, characterised by the
temperature $T$ and partial mass densities $\rho_{i}$ where $i=1...N$. The
fluid’s thermodynamic properties are fully described by two functions: the
internal energy $e(\rho_{1}...\rho_{N},T)$ and entropy
$s(\rho_{1}...\rho_{N},T)$ – both specific, i.e., per unit mass. The
dependence of the fluid pressure $p$ on $(\rho_{1}...\rho_{N},T)$, or the
equation of state, is defined by

$p=\rho\sum_{i}\rho_{i}\left(\frac{\partial
e}{\partial\rho_{i}}-T\frac{\partial s}{\partial\rho_{i}}\right),$ (1)

(e.g. Giovangigli & Matuszewski, 2013), where

$\rho=\sum_{i}\rho_{i}$ (2)

is the total density (here and hereinafter, the summation is implied to be
from $1$ to $N$ unless stated otherwise). The partial chemical potentials, in
turn, are given by

$G_{i}=\frac{\partial\left(\rho
e\right)}{\partial\rho_{i}}-T\frac{\partial\left(\rho
s\right)}{\partial\rho_{i}}.$ (3)

Note that $e(\rho_{1}...\rho_{N},T)$ and $s(\rho_{1}...\rho_{N},T)$ are not
fully arbitrary, but should satisfy the fundamental thermodynamic relation, or
the Gibbs relation, which can be written in the form

$\frac{\partial e}{\partial T}=T\frac{\partial s}{\partial T}.$ (4)

The equivalence of this equality and the standard form of the Gibbs relation
is shown in Appendix A.

Using equations (1)–(4), one can derive the following identities:

$\frac{\partial p}{\partial T}=\sum_{i}\rho_{i}\frac{\partial G_{i}}{\partial
T}+\rho s,$ (5) $\frac{\partial
p}{\partial\rho_{j}}=\sum_{i}\rho_{i}\frac{\partial G_{i}}{\partial\rho_{j}},$
(6)

then rewrite (6) in the form

$G_{i}=\frac{\partial}{\partial\rho_{i}}\left(\sum_{j}\rho_{j}G_{j}-p\right).$
(7)

In what follows, one also needs the specific heat capacity at constant volume,

$c=\frac{\partial e}{\partial T},$ (8)

where the traditional subscript V is omitted. In this paper, $c$ is assumed to
be positive, which is indeed the case for neutral fluids (Lynden-Bell, 1999).

Define also

$a_{i}=-\frac{\partial e}{\partial\rho_{i}},$ (9)

which will be referred to as the generalised van der Waals parameter (of the
$i$-th species), and

$B=-p-\rho\sum_{i}\rho_{i}a_{i}.$ (10)

$B$ is not one of the standard thermodynamic functions, but is convenient to
use when thermodynamics is coupled to hydrodynamics: as seen later, $B$
characterises the release (consumption) of heat due to the fluid’s mechanical
compression (expansion). Using equation of state (1) and definition (9) of the
van der Waals parameter, one can rearrange expression (10) in the form

$B=\rho T\sum_{i}\rho_{i}\frac{\partial s}{\partial\rho_{i}}.$ (11)

###  2.2 Examples: the Enskog–Vlasov and van der Waals fluids

In the low-density limit, the specific internal energy and entropy of any
fluid should tend to those of ideal gas, i.e.,

$e\sim\frac{T}{\rho}\sum_{i}c_{i}\rho_{i},\qquad s\sim\frac{\ln
T}{\rho}\sum_{i}c_{i}\rho_{i}-\frac{1}{\rho}\sum_{i}R_{i}\rho_{i}\ln\rho_{i}\qquad\text{as}\qquad\rho_{i}\rightarrow
0,$

where $c_{i}$ is the specific partial heat capacity, and $R_{i}=k_{B}/m_{i}$
is the specific gas constant ($m_{i}$ is the molecular mass of the $i$-th
species, and

$k_{B}=1.380649\times
10^{-23}\mathrm{m}^{2}\mathrm{s}^{-2}\mathrm{kg\,K}^{-1}$

is the Boltzmann constant). Note also that one can replace $\ln
T\rightarrow\ln T/\bar{T}$ and $\ln\rho\rightarrow\ln\rho/\bar{\rho}$, where
are $\bar{T}$ and $\bar{\rho}$ are suitable dimensional scales. This would
make the arguments of the logarithms nondimensional (with none of physically
measurable parameters depending on $\bar{T}$ and $\bar{\rho}$).

A simple description of non-ideal fluids is delivered by the Enskog–Vlasov
(EV) model, according to which

$e=\frac{T}{\rho}\sum_{i}c_{i}\rho_{i}-\frac{1}{\rho}\sum_{ij}a_{ij}\rho_{i}\rho_{j},$
(12) $s=\frac{\ln
T}{\rho}\sum_{i}c_{i}\rho_{i}-\frac{1}{\rho}\sum_{i}R_{i}\rho_{i}\ln\rho_{i}-\Theta(\rho_{1}...\rho_{N}),$
(13)

where $\Theta(\rho_{1}...\rho_{N})$ is an arbitrary analytic function, and
$a_{ij}$ characterises the long-range attraction between molecules of the
$i$-th and $j$-th species (with the implication that $a_{ij}=a_{ji}$). For a
pure fluid, $a_{11}$ coincides with the generalised van der Waals parameter
defined by (9) – hence, the notation.

One can readily verify that expressions (12)–(13) satisfy the Gibbs relation
(4). They can be viewed as two-term expansions in small $T$, under an extra
assumption that the heat capacity $c$ and the generalised van der Waals
parameter $a_{i}$ are independent of $T$ and are zeroth-degree homogeneous
functions of $\rho_{i}$ (for a pure fluid, $c$ and $a_{1}$ do not depend on
$\rho_{1}$ at all).

The EV fluid model originates from the Enskog–Vlasov kinetic theory (de
Sobrino, 1967; Grmela, 1971) and, as such, is naturally suited for the use
with the DIM which can be viewed as a hydrodynamic approximation of the EV
kinetic equation (Giovangigli, 2020, 2021). Equations (12)–(13) work very well
for inert gases (Benilov & Benilov, 2019), they will be shown to work
reasonably well for water, nitrogen, and oxygen (section 8 of this paper).

The EV involves too many parameters to be used as an illustration of general
results – so, in such cases, the simpler van der Waals model will be employed.
It is a particular case of equations (12)–(13) with

$\Theta=-\frac{1}{\rho}\sum_{i}R_{i}\rho_{i}\ln\left(1-{\displaystyle\sum\limits_{j}}b_{j}\rho_{j}\right),$
(14)

where $b_{i}$ is the reciprocal of the maximum density of the $i$-th species.
Physically, $1/b_{i}$ can be interpreted as density of the closest packing.

Given expressions (12)–(14), definition (1) of the pressure and definition (3)
of the chemical potential yield

$p=\frac{T}{1-{\displaystyle\sum\limits_{j}}b_{j}\rho_{j}}\sum_{i}R_{i}\rho_{i}-\sum_{ij}a_{ij}\rho_{i}\rho_{j},$
(15)
$G_{i}=TR_{i}\ln\frac{\rho_{i}}{1-{\displaystyle\sum\limits_{j}}b_{j}\rho_{j}}+\frac{Tb_{i}{\displaystyle\sum\limits_{j}}R_{j}\rho_{j}}{1-{\displaystyle\sum\limits_{j}}b_{j}\rho_{j}}-2\sum_{j}a_{ij}\rho_{j}+T\left(R_{i}+c_{i}-c_{i}\ln
T\right).$ (16)

For a pure fluid ($N=1$), equation (15) reduces to the classical van der Waals
equation of state (van der Waals, 1893).

###  2.3 Hydrodynamics

Consider a fluid flow characterised by the species densities
$\rho_{i}(\mathbf{r},t)$, mass-averaged velocity $\mathbf{v}(\mathbf{r},t)$,
and temperature $T(\mathbf{r},t)$, where $\mathbf{r}=\left[x,y,z\right]$ is
the position vector and $t$, the time. Let the species be affected by forces
$\mathbf{F}_{i}$ (which will be later identified with the van der Waals
forces), and the fluid as a whole, affected by viscosity. The shear viscosity
$\mu_{s}$ and bulk viscosity $\mu_{b}$ depend generally on $\rho_{i}$ and $T$.

Let the flow be governed by following equations:

$\frac{\partial\rho_{i}}{\partial
t}+\mathbf{\nabla}\cdot\left(\rho_{i}\mathbf{v}+\mathbf{J}_{i}\right)=0,$ (17)
$\frac{\partial\left(\rho\mathbf{v}\right)}{\partial
t}+\mathbf{\nabla}\cdot\left(\rho\mathbf{vv}\right)=\mathbf{\nabla}\cdot\boldsymbol{\Pi}+\sum_{i}\rho_{i}\mathbf{F}_{i}-\mathbf{\nabla}p,$
(18) $\frac{\partial}{\partial t}\left(\rho
e+\frac{1}{2}\rho\left|\mathbf{v}\right|^{2}\right)+\mathbf{\nabla}\cdot\left[\left(\rho
e+\frac{1}{2}\rho\left|\mathbf{v}\right|^{2}+p\right)\mathbf{v}-\boldsymbol{\Pi}\cdot\mathbf{v}+\mathbf{Q}\right]\\\
=\sum_{i}\mathbf{F}_{i}\cdot\left(\rho_{i}\mathbf{v}+\mathbf{J}_{i}\right).$
(19)

Here, the viscous stress tensor is

$\boldsymbol{\Pi}=\mu_{s}\left[\boldsymbol{\boldsymbol{\nabla}}\mathbf{v}+\left(\boldsymbol{\boldsymbol{\nabla}}\mathbf{v}\right)^{tr}-\frac{2}{3}\mathbf{I}\left(\boldsymbol{\boldsymbol{\nabla}}\cdot\mathbf{v}\right)\right]+\mu_{b}\,\mathbf{I}\left(\boldsymbol{\boldsymbol{\nabla}}\cdot\mathbf{v}\right),$
(20)

where the dotless product of two vectors (e.g.,
$\boldsymbol{\boldsymbol{\nabla}}\mathbf{v}$) produces a second-order tensor,
and the superscript tr denotes transposition. The diffusion fluxes
$\mathbf{J}_{i}$ and the heat flux $\mathbf{Q}$ are related to the forces
$\mathbf{F}_{i}$, temperature $T$, and chemical potentials $G_{i}$ by

$\mathbf{J}_{i}=\sum_{j}D_{ij}\left[\mathbf{F}_{j}-T\mathbf{\nabla}\left(\frac{G_{j}}{T}\right)\right]-\frac{\zeta_{i}}{T}\mathbf{\nabla}T,$
(21)
$\mathbf{Q}=\sum_{j}\zeta_{j}\left[\mathbf{F}_{j}-T\mathbf{\nabla}\left(\frac{G_{j}}{T}\right)\right]-\frac{\kappa}{T}\mathbf{\nabla}T,$
(22)

where $D_{ij}$, $\zeta_{i}$, and $\kappa$ are the transport coefficients (all
three depend generally on $\rho_{i}$ and $T$).

To understand the physical meaning of expressions (21)–(22), rewrite them in
the form

$\mathbf{J}_{i}=\sum_{j}D_{ij}\mathbf{F}_{j}-~{}\underset{\text{diffusion}}{\underbrace{\sum_{j}D_{ij}^{\prime}\mathbf{\nabla}\rho_{j}}}~{}-~{}\underset{\text{Soret
effect}}{\underbrace{\zeta_{i}^{\prime}\rho\mathbf{\nabla}T}},$ (23)
$\mathbf{Q}=\sum_{j}\zeta_{j}\mathbf{F}_{j}-~{}\underset{\text{Dufour
effect}}{\underbrace{\sum_{j}\zeta_{j}\sum_{k}\frac{\partial
G_{j}}{\partial\rho_{k}}\mathbf{\nabla}\rho_{k}}~{}}-~{}\underset{\text{heat
conduction}}{\underbrace{\kappa^{\prime}\mathbf{\nabla}T}},$ (24)

where

$D_{ik}^{\prime}=\sum_{j}D_{ij}\frac{\partial
G_{j}}{\partial\rho_{k}},\qquad\zeta_{i}^{\prime}=\frac{\zeta_{i}}{T\rho}+\sum_{j}\frac{D_{ij}}{\rho}\left(\frac{\partial
G_{j}}{\partial T}-\frac{G_{j}}{T}\right),$
$\kappa^{\prime}=\frac{\kappa}{T}+\sum_{j}\zeta_{j}\left(\frac{\partial
G_{j}}{\partial T}-\frac{G_{j}}{T}\right)$

are the standard diffusivities, thermodiffusivities, and thermal conductivity,
respectively. The second term in expression (23) corresponds to the classical
Fick Law (the fluxes depend linearly on the density gradients), and the last
term in expression (24) characterises heat conduction described by the Fourier
Law. The last term in (23) describes the Soret effect ($\mathbf{\nabla}T$
gives rise to diffusion) and the second term in (24), the Dufour effect
($\mathbf{\nabla}\rho_{j}$ gives rise to heat conduction).

The same four effects – diffusion, heat conduction, the Soret and Dufour
effects – are described, obviously, by the original expressions (21)–(22),
albeit in a form where the terms cannot be matched to a single effect each.

One might think that representing the fluxes in terms of
$\mathbf{\nabla}\rho_{j}$ would be more natural than using
$\mathbf{\nabla}\left(G_{j}/T\right)$. Observe, however, that the coefficient
of $\left(\mathbf{\nabla}T\right)/T$ in expression (21) coincides with the
coefficient of
$\left[\mathbf{F}_{j}-T\mathbf{\nabla}\left(G_{j}/T\right)\right]$ in (22).
This symmetry reflects the so-called Onsager reciprocal relations (Ferziger &
Kaper, 1972), which also imply

$D_{ij}=D_{ji},$

i.e., the diffusion of an $i$-th species in a $j$-th species occurs the same
way as that of the $j$-th species in the $i$-th species.

It should also be assumed that the extended transport matrix,

$D_{ij}^{(ext)}=\begin{bmatrix}&&&\zeta_{1}\\\ &D_{ij}&&\vdots\\\
&&&\zeta_{N}\\\ \zeta_{1}&\cdots&\zeta_{N}&\kappa\end{bmatrix}$

is positive semidefinite ($D_{ij}^{(ext)}\succeq 0$) – i.e.,

$\sum_{i=1}^{N+1}\sum_{j=1}^{N+1}d_{i}\,D_{ji}^{(ext)}d_{j}\geq 0,$

for all $\left(N+1\right)$-dimensional arrays $d_{i}$. As seen later, this
property is essential for the entropy principle to hold.

Furthermore, the transport coefficients should be such that

$\sum_{i}D_{ij}=0,\qquad\sum_{i}\zeta_{i}=0.$ (25)

As a result, the density equation (17) and expressions (21) for the diffusion
fluxes imply that

$\frac{\partial\rho}{\partial
t}+\mathbf{\nabla}\cdot\left(\rho\mathbf{v}\right)=0,$ (26)

where $\rho$ is the total density given by (2). Observe that, for a pure
fluid, restrictions (25) can be satisfied only if $D_{11}=0$ and
$\zeta_{1}=0$, which means that pure fluids neither diffuse nor thermodiffuse.

Equations (17)–(22) have been first derived from the thermodynamics of
irreversible processes (Meixner, 1941). They were also derived from
statistical mechanics (Bearman & Kirkwood, 1958; Mori, 1958) and non-
equilibrium statistical thermodynamics (Keizer, 1987) – for more references,
see (Giovangigli, 1999). A derivation of the small-density version of
(17)–(22) from the Boltzmann kinetic equation can be found in any textbook on
kinetic theory (e.g. Ferziger & Kaper, 1972).

In all these cases, the derived expressions for the transport coefficients
automatically satisfy the Onsager relations and the rest of the properties
listed above.

Note also that a reduction of the above equations for a binary fluid with no
Soret and Dufour effects ($N=2$, $\zeta_{i}=0$) was used by Liu et al. (2016)
to show that such a model is able to describe the phase equilibrium for a real
binary mixture of $\mathrm{CO}_{2}$ and ethanol.

###  2.4 Alternative forms of the momentum and energy equations

Since the transport fluxes (21)–(22) are expressed in terms of $G_{i}$ and $T$
(not $\rho_{i}$ and $T$), it is convenient to do the same for the pressure
gradient in the momentum equation (18). Recalling identities (6)–(5) which
imply

$\mathbf{\nabla}p=\sum_{i}\rho_{i}\mathbf{\nabla}G_{i}+\rho
s\mathbf{\nabla}T,$

and using equation (26) to simply the left-hand side of the momentum equation
(18), one reduces it to

$\rho\left[\frac{\partial\mathbf{v}}{\partial
t}+\left(\mathbf{v}\cdot\mathbf{\nabla}\right)\mathbf{v}\right]=\mathbf{\nabla}\cdot\boldsymbol{\Pi}+\sum_{i}\rho_{i}\left(\mathbf{F}_{i}-\mathbf{\nabla}G_{i}\right)-\rho
s\mathbf{\nabla}T.$ (27)

The energy equation, in turn, can be conveniently rewritten in terms of the
temperature (which is a measurable quantity, unlike the internal energy $e$).
Replacing, thus, in Eq. (19),

$\frac{\partial e}{\partial t}=\sum_{i}\frac{\partial
e}{\partial\rho_{i}}\frac{\partial\rho_{i}}{\partial t}+\frac{\partial
e}{\partial T}\frac{\partial T}{\partial t},$

one should use the density equation to eliminate $\partial\rho_{i}/\partial
t$. Using then equations (17) and (27), and recalling identities (8)–(9) and
(11), one obtains

$\rho c\left(\frac{\partial T}{\partial
t}+\mathbf{v}\cdot\mathbf{\nabla}T\right)+\mathbf{\nabla}\cdot\mathbf{Q}=\boldsymbol{\Pi}:\left(\mathbf{\nabla
v}\right)+B\mathbf{\nabla}\cdot\mathbf{v}+\sum_{i}\left(\mathbf{F}_{i}-\rho
a_{i}\mathbf{\nabla}\right)\cdot\mathbf{J}_{i},$ (28)

where the symbol “$:$” denotes double scalar product of two tensors.

The first term on the right-hand side of equation (28) describes heat
production by viscosity and the second term, that by fluid compression or
expansion.

###  2.5 The van der Waals force

Assume that a molecule of a $j$-th species exerts on a molecule of an $i$-th
species an isotropic force with a potential $\Phi_{ij}(r)$ where
$r=\left(x^{2}+y^{2}+z^{2}\right)^{1/2}$. Assuming for simplicity that the
fluid is unbounded, one can write the mass-averaged force affecting the $i$-th
species in the hydrodynamic equations (18)–(19) in the form

$\mathbf{F}_{i}(\mathbf{r},t)=\mathbf{\nabla}\sum_{j}\int\frac{\rho_{j}(\mathbf{r}^{\prime},t)}{m_{i}m_{j}}\Phi_{ij}(|\mathbf{r}^{\prime}-\mathbf{r|})\,\mathrm{d}^{3}\mathbf{r}^{\prime},$
(29)

where $m_{i}$ is the molecular mass an, the integration is implied to be over
the whole space. To guarantee the convergence of the integral in (29) and
those arising later, the potential $\Phi_{ij}(r)$ is assumed to decay
exponentially as $r\rightarrow\infty$.

Next, let the spatial scale of $\rho(\mathbf{r},t)$ be much larger than that
of $\Phi_{ij}(r)$, in which case expression (29) can be simplified
asymptotically. To do so, change in it
$\mathbf{r}^{\prime}\rightarrow\mathbf{r}^{\prime}+\mathbf{r}$ and then expand
$\rho_{j}(\mathbf{r}^{\prime}+\mathbf{r},t)$ about $\mathbf{r}^{\prime}$,
which yields

$\mathbf{F}_{i}(\mathbf{r},t)=\mathbf{\nabla}\sum_{j}\int\left[\rho_{j}(\mathbf{r},t)+\mathbf{r}^{\prime}\cdot\mathbf{\nabla}\rho_{j}(\mathbf{r},t)+\frac{1}{2}\mathbf{r}^{\prime}\mathbf{r}^{\prime}:\mathbf{\nabla\nabla}\rho_{j}(\mathbf{r},t)+\cdots\right]\frac{\Phi_{ij}(r^{\prime})}{m_{i}m_{j}}\mathrm{d}^{3}\mathbf{r}^{\prime}.$

Given the isotropy of $\Phi_{ij}(r^{\prime})$, the second integral in the
above expansion vanishes, and one obtains

$\mathbf{F}_{i}=\sum_{j}W_{ij}\mathbf{\nabla}\rho_{j}+\sum_{j}K_{ij}\mathbf{\nabla}\nabla^{2}\rho_{j}+\cdots,$
(30)

where

$W_{ij}=\int\frac{\Phi_{ij}(r^{\prime})}{m_{i}m_{j}}\mathrm{d}^{3}\mathbf{r}^{\prime},\qquad
K_{ij}=\int\frac{r^{\prime
2}}{2}\frac{\Phi_{ij}(r^{\prime})}{m_{i}m_{j}}\mathrm{d}^{3}\mathbf{r}^{\prime}.$

Since Newton’s Third Law implies that $\Phi_{ij}=\Phi_{ji}$, the matrices
$W_{ij}$ and $K_{ij}$ are symmetric.

Once expansion (30) is substituted into the hydrodynamic equations (18)–(19),
its first term can be absorbed into the internal energy – i.e., eliminated by
the change

$e\rightarrow e+\frac{1}{2\rho}\sum_{ij}W_{ij}\rho_{i}\rho_{j},\qquad
G_{i}\rightarrow G_{i}+\sum_{j}W_{ij}\rho_{j}.$

This does not come as a surprise, as the energy associated with potential
interactions of molecules can be viewed as a kind of internal energy – in
fact, the second term of expansion (30) could also be (and sometimes is)
absorbed into $e$. This is not done in this paper, however, as it would make
$e$ a functional (instead of a function), with the implication that all the
thermodynamic definitions and identities in Sect. 2.1 would need to be
rewritten in terms of functional derivatives.

Thus, without loss of generality, one can set in expression (30) $W_{ij}=0$.
Omitting also the small terms hidden in “$\cdots$”, one obtains

$\mathbf{F}_{i}=\sum_{j}K_{ij}\mathbf{\nabla}\nabla^{2}\rho_{j},$ (31)

which is the multicomponent extension of the standard DIM formula for the van
der Waals force (e.g., Pismen & Pomeau, 2000). The matrix $K_{ij}$ is the
extension of the so-called Korteweg parameter for pure fluids, and it will be
referred to as the Korteweg matrix. It should be positive definite,
$K_{ij}\succ 0$, as the van der Waals force should be attractive, not
repulsive.

Since the original representation (29) was a phenomenological model and, thus,
the pair-wise potential $\Phi_{ij}$ cannot be measured, the Korteweg matrix
should be viewed as a set of adjustable parameters. As seen later, they can be
deduced from the measurements of the equation of state and surface tension.
One should keep in mind, however, that the Korteweg matrix should not depend
on the temperature. Such a dependence would be physically unjustified, as the
intermolecular attraction (characterised by $K_{ij}$) should not depend on the
molecules’ velocities (characterised by $T$).

###  2.6 Boundary conditions at a solid wall

Let the fluid occupy a domain $\mathcal{D}$, bounded by a smooth surface
$\partial\mathcal{D}$. For simplicity, the so-called Navier slip is disallowed
in this work, so that the boundary condition for the velocity is

$\mathbf{v}=\mathbf{0}\qquad\text{at}\qquad\mathbf{r}\in\partial\mathcal{D}.$
(32)

To ensure mass conservation, one should require

$\mathbf{n\cdot
J}_{i}=0\qquad\text{at}\qquad\mathbf{r}\in\partial\mathcal{D},$ (33)

where $\mathbf{n}$ is the outward-pointing unit normal to
$\partial\mathcal{D}$.

The boundary condition for the temperature, in turn, depends on the problem at
hand. Since this paper is concerned, _inter alia_ , with the entropy principle
and energy conservation, it will be assumed that no heat escapes through the
boundary,

$\mathbf{n\cdot Q}=0\qquad\text{at}\qquad\mathbf{r}\in\partial\mathcal{D}.$
(34)

Boundary conditions (32)–(34) would be sufficient for the standard
compressible multi-component hydrodynamics, but the DIM requires an extra
condition [due to the presence of higher-order derivatives of the density
field in expression (31)].

The most common version of such condition – prescribing a linear combination
of the boundary-value of the density and its normal gradient – ascends to the
paper by Cahn & Hilliard (1958). In application to pure fluids, the Neumann
reduction of this condition was proposed by Seppecher (1996) and the Dirichlet
reduction, by Pismen & Pomeau (2000). As shown by Benilov (2020b), the latter
follows from the assumptions under which the whole DIM is derived (pair-wise
intermolecular interactions, slowly-varying density field), and so it is used
in the present paper. Thus, require that

$\rho_{i}=\rho_{0,i}\qquad\text{at}\qquad\mathbf{r}\in\partial\mathcal{D},$
(35)

where the constant $\rho_{0,i}$ is specific to the fluid/solid combination
under consideration. The general version of the boundary condition for
$\rho_{i}$ is discussed briefly in Appendix B.

To clarify the physical meaning of condition (35), consider the van der Waals
forces acting on a fluid molecule located infinitesimally close to the wall:
the solid attracts it _towards_ the wall, while the other fluid molecules pull
it _away_. The former force is fixed, whereas the latter grows with the near-
wall density – so that the balance is achieved when the density assumes a
certain value – namely, the parameter $\rho_{0,i}$ in condition (35). This
argument suggests that a smaller value of $\rho_{0,i}$ corresponds to a
hydrophobic wall (characterised by a large contact angle) and larger
$\rho_{0,i}$, to a hydrophilic one.

According to its physical meaning, $\rho_{0,i}$ does not depend on the
temperature. As shown later, its value can be deduced from a measurement of
the contact angle.

##  3 The entropy principle

###  3.1 Conservation laws and the H-theorem

It can be verified that the governing equations and boundary conditions
introduced above conserve the mass of each species

$M_{i}=\int_{\mathcal{D}}\rho_{i}\mathrm{d}^{3}\mathbf{r},$ (36)

and the total energy

$E=\int_{\mathcal{D}}\left[\rho
e+\frac{1}{2}\rho\left|\mathbf{v}\right|^{2}+\frac{1}{2}\sum_{ij}K_{ij}\left(\mathbf{\nabla}\rho_{i}\right)\cdot\left(\mathbf{\nabla}\rho_{j}\right)\right]\mathrm{d}^{3}\mathbf{r}.$
(37)

The three terms in expression (37) represent the internal, kinetic, and van
der Waals energies.

The governing equations and boundary conditions satisfy also an H-theorem,
reflecting the fact that the net entropy of a fluid in a thermally-insulated
container cannot decrease. To prove this, consider the following combination
of the governing equations:

$\left(\ref{2.28}\right)+T\rho\sum_{i}\frac{\partial
s}{\partial\rho_{i}}\times(\ref{2.17})+Ts\times\left(\ref{2.26}\right).$

After straightforward algebra involving the use of the thermodynamic
identities presented in subsection 2.1 and expressions (21)–(22) for the
fluxes, one obtains

$\frac{\partial\left(\rho s\right)}{\partial t}+\mathbf{\nabla}\cdot\left(\rho
s\mathbf{v}-\sum_{i}\frac{G_{i}}{T}\mathbf{J}_{i}+\frac{\mathbf{Q}}{T}\right)=\frac{\boldsymbol{\Pi}:\left(\mathbf{\nabla
v}\right)}{T}\\\
+\left\\{\sum_{ij}D_{ij}\left|\mathbf{F}_{j}-T\mathbf{\nabla}\left(\frac{G_{j}}{T}\right)\right|^{2}+2\sum_{i}\zeta_{i}\left[\mathbf{F}_{i}-T\mathbf{\nabla}\left(\frac{G_{i}}{T}\right)\right]\cdot\left(-\frac{\mathbf{\nabla}T}{T}\right)\phantom{\left|-\frac{\mathbf{\nabla}T}{T}\right|^{2}}\right.\\\
+\left.\vphantom{\sum_{ij}\left|\left(\frac{G_{j}}{T}\right)\right|^{2}}\kappa
T\left|-\frac{\mathbf{\nabla}T}{T}\right|^{2}\right\\}.$ (38)

The first term on the right-hand side of this equation is non-negative due to
the following (easily verifiable) identity,

$\boldsymbol{\Pi}:\left(\mathbf{\nabla
v}\right)=\frac{1}{2}\mu_{s}\left|\left(\mathbf{\nabla
v}\right)+\left(\mathbf{\nabla
v}\right)^{tr}-\frac{2}{3}\mathbf{I\nabla}\cdot\mathbf{v}\right|^{2}+\mu_{b}\left(\mathbf{\nabla}\cdot\mathbf{v}\right)^{2},$
(39)

whereas the expression in curly brackets is non-negative because the extended
transport matrix is positive semidefinite (see subsection 2.3). Thus,
integrating equation (38) over $\mathcal{D}$, and using boundary conditions
(32)–(34), one obtains

$\frac{\mathrm{d}S}{\mathrm{d}t}\geq 0,$ (40)

where

$S=\int_{\mathcal{D}}\rho s\,\mathrm{d}^{3}\mathbf{r}.$ (41)

Inequality (40) is the desired H-theorem.

It follows from (38)–(39) that the exact equality in (40) can only hold if the
velocity field is spatially uniform; together with the no-flow boundary
conditions, this requirement amounts to $\mathbf{v}=\mathbf{0}$ (i.e., the
fluid is static).

###  3.2 Stability via the entropy principle

The most common way to examine the stability of a steady solution of a set of
equations consist in linearizing these equations with respect to a small
perturbation, assuming the harmonic dependence of the perturbation on $t$, and
solving the resulting eigenvalue problem. In the problem at hand, however, it
is much simpler to examine stability using the entropy principle.

If, at a certain steady state, the total entropy $S$ has a local maximum
constrained by the conditions of fixed energy $E$ and mass $M_{i}$, this state
is stable. The inverse is also true: if $S$ does _not_ have a maximum, the
corresponding state is _un_ stable – because a perturbed solution with a
higher entropy cannot evolve ‘back’. Neutrally stable oscillations are also
prohibited by the entropy principle – hence, the system can only evolve
further away, towards a steady state where the entropy does have a maximum.

Let a fluid be enclosed in a container (which can be later assumed to be
infinitely large, if need be) and seek a maximum of $S$, constrained by the
conditions of fixed $M_{i}$ and $E$. This problem amounts to finding the
stationary points of the functional

$H[\rho_{1}...\rho_{N},T,\mathbf{v}]=S+\sum_{i}\eta_{i}M_{i}+\lambda E,$ (42)

where $\lambda$ and $\mu_{i}$ are the Lagrange multipliers, and $S$, $E$, and
$M_{i}$ are given by (41), (37), and (36), respectively. Varying $H$ with
respect to $\mathbf{v}$ and equating the variation to zero, one obtains

$\mathbf{v}=\mathbf{0},$ (43)

i.e., a steady state must be (unsurprisingly) static. Next, varying $H$ with
respect to $T$, one obtains

$\frac{\partial s}{\partial T}+\lambda\frac{\partial e}{\partial T}=0.$

Comparison of this equality with the Gibbs relation (4) yields

$\lambda=-\frac{1}{T}.$ (44)

Since $\lambda$ is a constant, equation (44) implies that $T$ is spatially
uniform – i.e., a steady state must be isothermal.

Finally, varying $H$ with respect to $\rho_{i}$, recalling expressions
(43)–(44) for $\mathbf{v}$ and $\lambda$, and keeping in mind definition (3)
of the chemical potential $G_{i}$, one obtains

$\frac{1}{T}\left(G_{i}-\sum_{j}K_{ij}\nabla^{2}\rho_{j}\right)+\eta_{i}=0.$
(45)

This equation describes all steady states of the governing equations, and it
will be extensively used in the rest of this paper. The temperature $T$ in
equation (45) should be treated as a known parameter, whereas the Lagrange
multiplier $\eta_{i}$ is to be deduced from the boundary conditions. The
latter will be illustrated for the case of an infinite domain, plus an
assumption that the fluid at infinity be spatially uniform and characterised
by a coordinate-independent chemical potential $G_{i}=G_{\infty,i}$. In this
case equation (45) yields

$\eta_{i}=-\frac{G_{\infty,i}}{T},$ (46)

as required.

Note also that equation (45) can also be recovered by adapting the governing
equations for the state of equilibrium. To do so, one should set
$\mathbf{v}=\mathbf{0}$, $T=\operatorname{const}$, and
$\partial\rho_{i}/\partial t=0$ in equations (17)–(18) and, recalling
expression (21) for the diffusion fluxes, obtain (45), as required.

To examine a solution of equation (45) for stability, one needs to calculate
the second variation of $H$. Omitting the algebra (involving the use of
identities (3)–(4), definition (8) of the heat capacity $c$, and expression
(44) for $\lambda$), one obtains

$\delta^{2}H=\frac{1}{T}\int_{\mathcal{D}}\left\\{\sum_{ij}\left[-\frac{\partial
G_{i}}{\partial\rho_{j}}\left(\delta\rho_{i}\right)\left(\delta\rho_{j}\right)-K_{ij}\left(\mathbf{\nabla}\delta\rho_{i}\right)\cdot\left(\mathbf{\nabla}\delta\rho_{j}\right)\right]-\frac{\rho
c}{T}\left(\delta
T\right)^{2}-\rho\left|\delta\mathbf{v}\right|^{2}\right\\}\mathrm{d}^{3}\mathbf{r}.$

Evidently, perturbations of the temperature and velocity are negative and can
only lower the total entropy – hence, the type of the stationary point
(maximum vs saddle) is fully determined by the variation of the density.
Setting, thus, $\delta T=0$ and $\delta\mathbf{v}=\mathbf{0}$, one obtains

$\delta^{2}H=\frac{1}{T}\int_{\mathcal{D}}\sum_{ij}\left[-\frac{\partial
G_{i}}{\partial\rho_{j}}\left(\delta\rho_{i}\right)\left(\delta\rho_{j}\right)-K_{ij}\left(\mathbf{\nabla}\delta\rho_{i}\right)\cdot\left(\mathbf{\nabla}\delta\rho_{j}\right)\right]\mathrm{d}^{3}\mathbf{r}.$
(47)

Expression (47) is the main tool used in this paper for studying the stability
properties of steady states described by equation (45).

##  4 Nondimensionalization and the governing parameters

Introduce a characteristic density scale $\varrho$, a pressure scale $P$, a
temperature scale $T_{0}$, and a typical value $R$ of the specific gas
constant $R_{i}$ introduced in subsection 2.2. These scales allow one to
nondimensionalise all thermodynamics variables and functions introduced in
subsection 2.1,

$\left(\rho_{i}\right)_{nd}=\frac{\rho_{i}}{\varrho},\qquad\rho_{nd}=\frac{\rho}{\varrho},\qquad
T_{nd}=\frac{T}{T_{0}},\qquad e_{nd}=\frac{\varrho}{P}e,\qquad
s_{nd}=\frac{s}{R},$
$p_{nd}=\frac{p}{P},\qquad\left(G_{i}\right)_{nd}=\frac{\rho}{P}G_{i},\qquad
c_{nd}=\frac{c}{R},\qquad a_{nd}=\frac{\varrho^{2}}{P}a,\qquad
B_{nd}=\frac{B}{P}.$

It is also convenient to nondimensionalise the coefficients which appear in
the governing equations. Using their respective scales, one obtains

$\left(K_{ij}\right)_{nd}=\frac{K_{ij}}{K},\qquad\left(\mu_{s}\right)_{nd}=\frac{\mu_{s}}{\mu},\qquad\left(\mu_{b}\right)_{nd}=\frac{\mu_{b}}{\mu},$
$\left(R_{i}\right)_{nd}=\frac{R_{i}}{R},\qquad\kappa_{nd}=\frac{\kappa}{\varkappa},\qquad\left(D_{ij}\right)_{nd}=\frac{D_{ij}}{D},\qquad\left(\zeta_{i}\right)_{nd}=\frac{\zeta_{i}}{\zeta}.$

Finally, introduce

$\mathbf{r}_{nd}=\frac{\mathbf{r}}{l},\qquad
t_{nd}=\frac{Vt}{l},\qquad\mathbf{v}_{nd}=\frac{\mathbf{v}}{V},$
$\boldsymbol{\Pi}_{nd}=\frac{l}{\mu
V}\boldsymbol{\Pi},\qquad\left(\mathbf{F}_{i}\right)_{nd}=\frac{l^{3}}{K\varrho}\mathbf{F}_{i},\qquad\left(\mathbf{J}_{i}\right)_{nd}=\frac{\mathbf{J}_{i}}{\varrho
V},\qquad\mathbf{Q}_{nd}=\frac{\mathbf{Q}}{PV},$

where

$l=\sqrt{\frac{\varrho^{2}K}{P}},$

is the characteristic interfacial thickness, and

$V=\frac{Pl}{\mu}$

is the velocity scale reflecting the three-way balance of the pressure
gradient, viscous stress, and van der Waals force (Benilov, 2020a). Note also
that the characteristic interfacial thickness $l$ is on a nano scale
(Magaletti et al., 2016; Gallo et al., 2020; Benilov, 2020b).

Rewriting equations (17), (27)–(28), (21)–(22), and (31), and omitting the
subscript nd, one obtains

$\frac{\partial\rho_{i}}{\partial
t}+\mathbf{\nabla}\cdot\left(\rho_{i}\mathbf{v}+\mathbf{J}_{i}\right)=0,$ (48)
$\framebox{$\alpha$}\rho\left[\frac{\partial\mathbf{v}}{\partial
t}+\left(\mathbf{v}\cdot\mathbf{\nabla}\right)\mathbf{v}\right]=\mathbf{\nabla}\cdot\boldsymbol{\Pi}+\sum_{i}\rho_{i}\left(\mathbf{F}_{i}-\mathbf{\nabla}G_{i}\right)-\rho
s\mathbf{\nabla}T,$ (49) $\framebox{$\tau$}\rho c\left(\frac{\partial
T}{\partial
t}+\mathbf{v}\cdot\mathbf{\nabla}T\right)+\mathbf{\nabla}\cdot\mathbf{Q}=\boldsymbol{\Pi}:\left(\mathbf{\nabla
v}\right)-B\mathbf{\nabla}\cdot\mathbf{v}+\sum_{i}\left(\mathbf{F}_{i}+\rho
a_{i}\mathbf{\nabla}\right)\cdot\mathbf{J}_{i},$ (50)
$\boldsymbol{\Pi}=\mu_{s}\left[\mathbf{\nabla v}+\left(\mathbf{\nabla
v}\right)^{tr}-\frac{2}{3}\left(\mathbf{\nabla}\cdot\mathbf{v}\right)\mathbf{I}\right]+\mu_{b}\left(\mathbf{\nabla}\cdot\mathbf{v}\right)\mathbf{I},$
(51)
$\framebox{$\delta$}\mathbf{J}_{i}=\sum_{j}D_{ij}\left[\mathbf{F}_{j}-T\mathbf{\nabla}\left(\frac{G_{j}}{T}\right)\right]-\framebox{$\dfrac{\nu\delta}{\beta}$}\frac{\zeta_{i}}{T}\mathbf{\nabla}T,$
(52)
$\framebox{$\beta$}\mathbf{Q}=\framebox{$\nu$}\sum_{j}\zeta_{j}\left[\mathbf{F}_{j}-T\mathbf{\nabla}\left(\frac{G_{j}}{T}\right)\right]-\kappa\mathbf{\nabla}T,$
(53) $\mathbf{F}_{i}=\mathbf{\nabla}\sum_{j}K_{ij}\nabla^{2}\rho_{j},$ (54)

where the nondimensional parameters

$\alpha=\frac{K\varrho^{3}}{\mu^{2}},\qquad\tau=\frac{\varrho
RT_{0}}{P},\qquad\beta=\frac{KP\varrho^{2}}{\mu\varkappa
T_{0}},\qquad\nu=\frac{\zeta P}{\varkappa
T_{0}\varrho},\qquad\delta=\frac{K\varrho^{4}}{D\mu P},$ (55)

are placed for better ‘visibility’ in boxes.

Since $\alpha$ appears in front of the material derivative in equation (48),
it should be interpreted as the microscopic Reynolds number (associated with
the flow near the interface, not the global flow). $\tau$ is the
nondimensional temperature. $\beta$ is the isothermality parameter (Benilov,
2020a): if it is small, the temperature field is close to being spatially
uniform (isothermal). The Nusselt number $\nu$ characterises the importance of
heat diffusion relative to heat conduction – see expression (53) for the heat
flux. Finally, the position of $\delta$ in equation (52) suggests that this
parameter characterises advection relative to diffusion (recall that the flux
$\mathbf{J}_{i}$ was nondimensionalised using the advection scale $\varrho
V$).

One can also introduce the Prandtl and Schmidt numbers,

$Pr=\frac{\beta\tau}{\alpha}=\frac{\mu R}{\varkappa},\qquad
Sc=\frac{\delta}{\alpha}=\frac{\mu\varrho}{DP},$

characterizing viscosity relative to heat conduction and diffusion,
respectively.

The nondimensional versions of the boundary conditions look exactly the same
as their dimensional counterparts, (32)–(35) – as do the nondimensional
versions of the thermodynamic identities of subsection 2.1 except definition
(8) of the heat capacity, which becomes

$c=\frac{1}{\tau}\frac{\partial e}{\partial T}.$

This paper does not aim to present a comprehensive classification of
asymptotic regimes of the multicomponent DIM (for the pure-fluid DIM, see
Benilov, 2020a). Only the simplest regime will be described and used later as
a qualitative illustration of theoretically-predicted behaviors. It
corresponds to the following assumptions:

$\alpha\ll 1,\qquad\beta\ll 1,\qquad\nu\ll 1.$

The smallness of $\alpha$ allows one to take advantage of the slow-flow
approximation – whereas the other two assumptions and equations (50) and (53)
imply that

$T=1+\mathcal{O}(\beta),$

i.e., the fluid is almost isothermal. Setting, thus, $T=\operatorname{const}$
in expression (52) for the diffusion flux and substituting it into the density
equation (48), one obtains

$\frac{\partial\rho_{i}}{\partial
t}+\mathbf{\nabla}\cdot\left[\rho_{i}\mathbf{v}+\sum_{j}D_{ij}\mathbf{\nabla}\left(\sum_{n}K_{jn}\nabla^{2}\rho_{n}-G_{j}\right)\right]=0,$
(56)

where it was assumed, without loss of generality, that $\delta=1$. Simplifying
similarly equation (49) and substituting into it expression (54) for
$\mathbf{F}_{i}$, one obtains

$\mathbf{\nabla}\cdot\boldsymbol{\Pi}+\sum_{i}\rho_{i}\mathbf{\nabla}\left(\sum_{j}K_{ij}\nabla^{2}\rho_{j}-G_{i}\right)=0.$
(57)

Observe that equations (56)–(57) do not involve the (small) temperature
variations – hence, the temperature equation (50) can be simply omitted.

Equations (56)–(57) and expression (51) for the viscous stress form a closed
set of approximate equations for the unknowns $\rho_{i}$ and $\mathbf{v}$. The
chemical potential $G_{i}$ in these equations should be treated as a known
function of $\rho_{1}...\rho_{N}$, and the temperature $T$, as a known
parameter.

Unlike the exact set – which describes fast acoustic waves and slow
interfacial flow – the approximate equations describe only the latter. This is
a clear advantage: in a numerical simulation, for example, waves necessitate a
small timestep and, thus, dramatically slow down the computation. At the same
time, the two sets of equations have very similar properties: they share the
same steady solutions, both conserve mass and energy, and satisfy the
H-theorem.

##  5 Basic solutions and their stability

###  5.1 Spatially uniform states

Consider a uniform fluid where there is no flow and all species are in vapour
phase. If the temperature drops, one of the species may become overcooled,
giving rise to rapid condensation. A similar instability may occur when all or
some of the species are in _liquid_ phase, and the temperature _increases_ –
giving rise to rapid _evaporation_.

To determine which exactly states are unstable (thus, do not occur in real
world), one could perform the usual linear analysis. For a pure fluid ($N=1$),
this is a straightforward task yielding the following stability criterion:

$\frac{\partial p}{\partial\rho}>0.$ (58)

Thus, instability occurs if an increase in density _lowers_ the pressure, so
that the flow generated by the pressure gradient brings even more fluid to
this region.

For an arbitrary $N$, however, the analysis of linearised equations is
extremely cumbersome. Instead, it will be examined via the entropy principle –
i.e., using expression (47). Taking the limit
$\mathcal{D}\rightarrow\mathbb{R}^{3}$ (unbounded fluid), assuming that
$\rho_{i}$ is spatially uniform, and recalling that $K_{ij}\succ 0$ and $c>0$,
one can deduce from (47) that $H$ has a maximum if

$\frac{\partial G_{i}}{\partial\rho_{j}}\succ 0.$ (59)

This is the standard stability criterion for a spatially uniform state of a
multicomponent fluid (Glansdorff & Prigogine, 1971).

Four comments are in order:

  * •

For a physically meaningful $G_{i}(\rho_{1}...\rho_{N},T)$, condition (59)
holds for a sufficiently rarefied vapour or a sufficiently dense liquid.

  * •

The states with marginal stability are sometimes referred to as “spinodal
points” and stable vapour, as “subspinodal vapour”.

  * •

Interestingly, the viscosity and transport coefficients do not appear in
criterion (59). The corresponding effects can only slow the instability down,
but not eliminate it.

  * •

To reconcile criterion (59) with its pure-fluid counterpart (58), note that
for $N=1$, $\partial G/\partial\rho$ and $\partial p/\partial\rho$ have the
same sign [as implied by identity (6)].

###  5.2 Two-phase saturated states

Consider an interface separating liquid and vapour of the same pure fluid. If
they are in equilibrium, their temperatures are equal, and the rest of the
parameters satisfy the so-called Maxwell construction (Maxwell, 1875),

$G(\rho^{(l)},T)=G(\rho^{(v)},T),\qquad p(\rho^{(l)},T)=p(\rho^{(v)},T),$ (60)

where the superscripts (l) and (v) mark the parameters of the liquid and
vapour, respectively. The former and latter equalities in (60) guarantee
thermodynamic and mechanical equilibria of the interface, respectively.

Assume also that both phases are stable,

$\left(\frac{\partial
G}{\partial\rho}\right)_{\rho=\rho^{(l)}}>0,\qquad\left(\frac{\partial
G}{\partial\rho}\right)_{\rho=\rho^{(v)}}>0,$

and the density of the liquid exceeds that of the vapour,

$\rho^{(l)}>\rho^{(v)}.$

Subject to these conditions, the Maxwell construction (60) uniquely determines
the _saturated_ densities $\rho^{(l)}$ and $\rho^{(v)}$ as functions of $T$.

In what follows, the multicomponent version of the Maxwell construction will
be shown to follow from the DIM’s entropy principle.

Consider an insulated container with fluid subdivided between two states,
liquid and vapour. If the liquid/vapour and fluid/wall interfaces are
sufficiently thin, the corresponding full entropy, mass, and energy can be
approximately written in the form

$S=V^{(l)}\rho^{(l)}\,s(\rho_{1}^{(l)}...\rho_{N}^{(l)},T^{(l)})+V^{(v)}\rho^{(v)}\,s(\rho_{1}^{(v)}...\rho_{N}^{(v)},T^{(v)}),$
$M_{i}=V^{(l)}\rho_{i}^{(l)}+V^{(v)}\rho_{i}^{(v)},$ (61)
$E=V^{(l)}\rho^{(l)}\,e(\rho_{1}^{(l)}...\rho_{N}^{(l)},T^{(l)})+V^{(v)}\rho^{(v)}\,e(\rho_{1}^{(v)}...\rho_{N}^{(v)},T^{(v)}),$
(62)

where $V^{(l)}$ and $V^{(v)}$ are the volumes of the liquid and vapour phases,
respectively. Introduce also the full volume of the container,

$V=V^{(l)}+V^{(v)}.$ (63)

The Maxwell construction can be derived by maximizing $S$ subject to the
constraints of fixed $M_{i}$, $E$, and $V$ (which are now _functions_ , as
opposed to being _functionals_ in the previous subsection). Straightforward
calculations show that the maximum of entropy is achieved if $T^{(l)}=T^{(v)}$
(isothermality) and

$G_{i}(\rho_{1}^{(l)}...\rho_{N}^{(l)},T)=G_{i}(\rho_{1}^{(v)}...\rho_{N}^{(v)},T),$
(64)
$p(\rho_{1}^{(l)}...\rho_{N}^{(l)},T)=p(\rho_{1}^{(v)}...\rho_{N}^{(v)},T).$
(65)

Four comments are in order:

  * •

Since the solution describing coexistence of two phases satisfies the entropy
principle, it is automatically stable.

  * •

The Maxwell construction can sometimes yield a solution with negative
$V^{(l)}$or $V^{(v)}$. In such cases, the two-phase equilibrium is irrelevant,
and the fluid evolves towards the one-phase state with the same masses of the
species and total energy.

  * •

The multicomponent Maxwell construction (64)–(65) comprises $N+1$ equations
for $2N$ unknowns – hence, does not uniquely fix the liquid and vapour
densities.
To close set (64)–(65), one should assume the masses $M_{i}$ and the total
energy $E$ to be known and view equalities (61)–(63) as additional equations.
They bring the total number of equations to $2N+3$ (making the ‘expanded’
Maxwell construction appear overdetermined) – but in this formulation,
$V^{(l)}$, $V^{(v)}$, and $T$ should also be viewed as unknowns.
Physically, if a certain amount of fluid, with a certain amount of energy, is
placed in a box, the entropy principle uniquely determines the final
temperature and the proportion in which the box is subdivided between the
liquid and vapour phases.

  * •

The closure of the Maxwell construction described in the previous bullet is,
obviously, inapplicable to containers of infinite volume. To understand how
conditions (64)–(65) should be closed in this case, consider the interface
between the Earth’s atmosphere and ocean. For this setting, one should
prescribe the (atmospheric) pressure and composition of dry air above the
ocean’s surface. With these parameters given, (64)–(65) yield the saturated
moisture content of the air and the saturated amounts of nitrogen, oxygen,
_etc._ dissolved in the oceanic water.

To illustrate the use of the Maxwell construction, consider a van der Waals
fluid, whose pressure and chemical potential are described by expressions
(15)–(16). Assume for simplicity that the fluid is pure ($N=1$) and monatomic
($c_{i}=3R_{i}/2$) – and let the scales $\varrho$, $P$, and $R$ used for
nondimensionalization be such that

$a_{11}=1,\qquad b_{1}=1,\qquad R_{1}=1.$ (66)

If $T<8/27$ (which is a _subcritical_ temperature of the van der Waals fluid),
two states exist representing the liquid and vapour phases, with some (but not
all) of the states in between being unstable (see figure 1). If the
temperature is _supercritical_ , only one phase exists and interfaces do not.

Figure 1: The pressure [panel (a)] and chemical potential [panel (b)] of a
pure van der Waals fluid as functions of the density, for three values of the
temperature: (1) $T=0.26$ (subcritical); (2) $T=8/27$ (critical), (3) $T=0.33$
(supercritical). The critical point on curve (2) is marked with a filled
circle; the saturation and spinodal points on curve (1) are shown by filled
and empty squares, respectively. The labels “undersat[urated] vapour”,
“oversat[urated] v[apour]”, _etc._ apply only to curve (1).

The simplest model describing a water/air interface is that with $N=2$. The
Maxwell construction in this case should be complemented with one extra
condition, setting the pressure above the interface equal to its atmospheric
value,

$p(\rho_{1}^{(v)},\rho_{2}^{(v)},T)=p_{A}.$

The difference between a pure and a multicomponent fluids is illustrated in
figure 2 for parameters (66) and

$a_{22}=0.2,\qquad a_{12}=0,\qquad b_{2}=1,\qquad R_{2}=0.6,\qquad
p_{A}=0.03.$ (67)

These values reflect a compromise between simplicity, illustrative purposes
(the curves with $N=1$ and $N=2$ should be visibly different), and an attempt
to loosely match the parameters of water and air (loosely, because the van der
Waals model includes few adjustable constants). In particular, parameters (67)
make the critical temperature of the second species noticeably smaller than
that of the first species. As a result, a temperature range exists where the
former is definitely vapour, whereas the latter can be either vapour or liquid
(as is indeed the case with water under normal conditions).

Figure 2: The saturated densities [panel (a)] and pressure [panel (b)] vs the
temperature, for a pure (solid line) and two-component (dotted line) van der
Waals fluids, with the parameters (66)–(67). The empty circles mark the
boiling point, the filled circle marks the critical point. The saturated
densities of the second species (of the two-component fluid) are not shown.

The most important difference between a multicomponent fluid at fixed pressure
and a pure fluid is that the former exhibits a boiling point. It occurs when
the saturated pressure of the liquid-phase species matches the applied
(atmospheric) pressure: as a result, the second species is _completely_
replaced by the vapour of the first species – and so the fluid becomes pure.
If it is heated beyond the boiling point, the pressure can no longer be kept
fixed, but ought to increase (to match the saturated pressure of the first
species, which grows with $T$).

One should keep in mind that the small-$T$ part of figure 2 is physically
meaningless due to freezing (which is not described by the DIM).

##  6 One-dimensional steady states

###  6.1 Flat liquid/vapour interfaces

Consider an unbounded fluid involving liquid and vapour phases in equilibrium,
separated by a static flat interface. Its spatial structure is described by
equation (45) derived from the entropy principle. To adapt it specifically for
a flat interface, let $\rho_{i}$ depend on a single coordinate – say, $z$ – so
that the interface is horizontal. Let the fluid below the interface be liquid
and above, vapour,

$\rho_{i}\rightarrow\rho_{i}^{(l)}\qquad\,\text{as}\qquad
z\rightarrow-\infty,$ (68)
$\rho_{i}\rightarrow\rho_{i}^{(v)}\qquad\text{as}\qquad z\rightarrow+\infty.$
(69)

Taking in equation (45) the limit $z\rightarrow+\infty$, one can determine the
constant $\eta_{i}$ and rewrite (45) in the form

$\sum_{j}K_{ij}\frac{\mathrm{d}^{2}\rho_{j}}{\mathrm{d}z^{2}}=G_{i}-G_{i}(\rho_{1}^{(v)}...\rho_{N}^{(v)},T).$
(70)

Boundary-value problem (68)–(70) is invariant with respect to the change
$z\rightarrow z+\operatorname{const}$ – hence, its solution is not unique. To
make it such, an extra condition is needed – say,

$\rho_{1}=\frac{1}{2}\left(\rho_{1}^{(l)}+\rho_{1}^{(v)}\right)\qquad\text{at}\qquad
z=0.$ (71)

As shown before, the parameters $\rho_{i}^{(l)}$ and $\rho_{i}^{(v)}$ are not
arbitrary but have to satisfy the Maxwell construction, i.e., conditions
(64)–(65). Thus, they must be intrinsic to boundary-value problem (68)–(71)
and can be derived directly from it, as a condition for existence of solution.

Indeed, consider equation (70) in the limit $z\rightarrow-\infty$ and recall
(68) – which immediately yields the first half of the Maxwell construction,
equation (64). To derive the second half, consider

$\sum_{i}\int_{-\infty}^{\infty}\left(\ref{6.3}\right)\times\frac{\mathrm{d}\rho_{i}}{\mathrm{d}z}\mathrm{d}z.$

The integral in this equation can be evaluated via identity (7), and the
constant of integration can be fixed via boundary condition (69). After
straightforward algebra, one obtains

$\frac{1}{2}\sum_{ij}K_{ij}\frac{\mathrm{d}\rho_{i}}{\mathrm{d}z}\frac{\mathrm{d}\rho_{j}}{\mathrm{d}z}=\sum_{i}\rho_{i}\left[G_{i}-G_{i}(\rho_{1}^{(v)}...\rho_{N}^{(v)},T)\right]-p+p(\rho_{1}^{(v)}...\rho_{N}^{(v)},T).$
(72)

Taking in this equation the limit $z\rightarrow-\infty$ and recalling boundary
condition (68), one recovers the second part of the Maxwell construction,
equation (65), as required.

Boundary-value problem (68)–(71) was solved numerically for a two-component
van der Waals fluid, using the MATLAB function BVP4c based on the three-stage
Lobatto IIIa formula (Kierzenka & Shampine, 2001). Examples of numerical
solutions with parameters (66)–(67) and

$K_{11}=K_{22}=1,\qquad K_{12}=0,$ (73)

are shown in figure 3. Observe that an increase in temperature makes the
interface thicker, and it becomes infinitely thick when the critical
temperature is reached.

Figure 3: Comparison between liquid/vapour interfaces in pure (solid line) and
two-component (dotted line) van der Waals fluids with parameters (66)–(67) and
(73). Curves (1)–(3) correspond to $T=0.17,~{}0.23,~{}0.28132$. The last of
these values is the (approximate) boiling point, so that the solid and doted
curves virtually coincide and $\rho_{2}$ is virtually zero.

As shown in Appendix C.1, all solutions describing liquid/vapour interfaces,
are stable as long as $\mathrm{d}\rho_{i}/\mathrm{d}z$ never vanishes
($\rho_{i}(z)$ is strictly monotonic) for all $i$.

###  6.2 One-dimensional drops and bubbles

Apart from liquid/vapour interfaces, equation (45) admits spatially localised
solutions, describing drops or bubbles floating in vapour or liquid,
respectively. In this subsection, the simplest – one-dimensional – solutions
of this kind are discussed, describing a flat layer of increased or decreased
density. They both correspond to the following boundary condition

$\rho_{i}\rightarrow\rho_{\infty,i}\qquad\text{as}\qquad
z\rightarrow\pm\infty,$ (74)

where $\rho_{\infty,i}$ is the density of the vapour (liquid). Unlike the
interface solutions examined above, the fluid outside the drop (bubble) does
not have to be saturated vapour (liquid); it can actually be any stable or
metastable state.

The 1D reduction of the steady-state equation (45) and expression (46) for
$\eta_{i}$ can be written in the form

$\sum_{j}K_{ij}\frac{\mathrm{d}^{2}\rho_{j}}{\mathrm{d}z^{2}}=G_{i}-G_{\infty,i},$
(75)

where $G_{\infty,i}=G_{i}(\rho_{\infty,1}...\rho_{\infty,N},T)$ is the
chemical potential at infinity.

Figure 4 illustrates typical solutions of boundary-value problem (74)–(75)
computed for a pure fluid. The following features can be observed in panel (a)
depicting some 1D drop solutions:

Figure 4: Solutions of boundary-value problem (74)–(75) for a pure van der
Waals fluid with $T=0.2$. (a) Drop solutions with
$\rho_{\infty,1}=0.1,~{}0.05,~{}0.03587,~{}0.0358485$ (curves (1)–(4),
respectively). (b) Bubble solutions with
$\rho_{\infty,1}=0.63,~{}0.68,~{}0.7286,~{}0.728671$ (curves (1)–(4),
respectively). The labels “u[ndersaturated] v[apour]”, “o[versaturated]
v[apour]”, _etc._ characterise $\rho_{\infty,1}$.

  * •

Drop solutions exist only if $\rho_{\infty,1}$ corresponds to oversaturated
vapour.

  * •

As $\rho_{\infty,1}\rightarrow\rho_{1}^{(v)}+0$, the drop becomes increasingly
thick.

  * •

Once $\rho_{\infty,1}$ passes $\rho_{1}^{(v)}$, drop solutions cease to exist.
Physically, this is because liquid drops surrounded by undersaturated vapour
evaporate.

Similar tendencies have been observed for 1D bubble solutions, illustrated in
figure 4(b).

Most importantly, all drops (bubbles) in oversaturated vapour (undersaturated
liquid) are likely to be unstable. For pure fluids ($N=1$), the instability
can be proven rigorously: as shown in Appendix C.2, the entropy has a saddle
(not maximum) on these solutions.

To understand how slightly perturbed drops and bubbles evolve, they were
simulated numerically using the simplest asymptotic version of the DIM – the
one based on equations (56)–(57) and (51). It is adapted for $N=1$ and a
single spatial coordinate in Appendix D, which also outlines the numerical
method used.

Various initial conditions have been simulated for the van der Waals pure
fluid, with only two patterns of evolution observed. If the initial condition
includes a steady-drop solution $\rho_{1}^{(sd)}(z)$ _plus some extra fluid_ ,
spontaneous condensation is typically triggered off, giving rise to two shock
waves propagating away from the drop’s center. This behavior is illustrated in
figure 5(a) for

$T=0.2,\qquad\rho_{\infty,1}=0.05,$ (76)

and the initial condition

$\rho_{1}=\rho_{1}^{(sd)}(1.01\,z),\qquad w=0\qquad\text{at}\qquad t=0,$ (77)

where $w$ is the vertical velocity (the other two velocity components are zero
due to the 1D nature of the flow).

Figure 5: Evolution of a perturbed drop, described by equations (134)–(135)
and boundary conditions (136)–(137), for a pure van der Waals fluid and
parameters (76). (a) Initial condition (77) (the steady solution plus extra
fluid); the time $t$ of a snapshot and the corresponding curve number $n$ are
inter-related via $t=5000\,n$. (b) Initial condition (78) (the steady solution
minus extra fluid); $t=500\,n$.

If, however, the initial condition includes _less_ fluid than
$\rho_{1}^{(nm)}(z)$, the drop would typically evaporate. This pattern is
illustrated in figure 5(b) for parameters (76) and the following initial
condition,

$\rho_{1}=\rho_{1}^{(sd)}(0.99\,z),\qquad w=0\qquad\text{at}\qquad t=0.$ (78)

The two patterns can be understood physically on the basis of the fact that
the steady-drop solutions exist only for drops surrounded by oversaturated
metastable vapour – which is stable with respect to infinitesimally small __
perturbations, but may be unstable with respect to finite ones. One can thus
assume that the drop solution provides a lower bound for the mass of
perturbations capable of triggering off instability of the surrounding vapour.

It is also interesting to see how a drop would evolve if it is surrounded by
_saturated_ (not _over_ saturated) vapour.

The fact that the thickness of a steady drop becomes infinite as
$\rho_{\infty,1}\rightarrow\rho_{1}^{(v)}$ suggests that a _finite_ -size
liquid layer, surrounded by saturated vapour, evaporates. The mechanism of
this evaporation, however, is not immediately clear. 2D and 3D drops, for
example, evaporate because the curvature of their boundary gives rise to a
liquid-to-vapour mass flux (Benilov, 2020d, 2021, 2022b) – but the boundaries
of 1D drops are flat. One can only assume that they evaporate due to a long-
range interaction of the drop’s upper and lower interfaces.

To explore this hypothesis, 1D drops floating in saturated vapour were
simulated numerically (using the same model and numerical method as before).
It has turned out that these drops do evaporate, and their evolution can be
subdivided into three distinct stages:

  1. (i)

The boundaries of the drop rapidly assume the profile of a steady
liquid/vapour interface – of the kind examined in the previous subsection.

  2. (ii)

The drop begins to get thinner, but the density of the drop’s core remains
close to $\rho_{1}^{(l)}$.

  3. (iii)

Once the thickness of the drop becomes comparable to the interfacial
thickness, the density at the drop’s center begins to rapidly decrease, and
the drop disappears.

The main characteristic of such behavior is the evaporation time $t_{e}$,
which can be defined as the interval over which the density at the drop’s
center falls by a factor of, say, $10$ – i.e.,

$\frac{\rho_{1}(0,t_{e})}{\rho_{1}(0,0)}=0.1.$

$t_{e}$ was computed as a function of the drop’s initial size $W$, for the
following initial condition:

$\rho_{1}=\rho_{1}^{(v)}+\frac{\rho_{1}^{(l)}-\rho_{1}^{(v)}}{1+\left(\dfrac{2z}{W}\right)^{6}},\qquad
w=0\qquad\text{at}\qquad t=0.$

It has turned out that, for $W$ changing from $0.2$ to $10$, the evaporation
time $t_{e}$ grows from $48.5$ to $10^{6}$ – i.e., exponentially. This agrees
with the hypothesis that the drop evaporation is caused by long-range
interaction of the drop’s boundaries and the fact that the density in
liquid/vapour interfaces tends to $\rho_{1}^{(v)}$ and $\rho_{1}^{(l)}$ (as
$z\rightarrow\pm\infty$) _exponentially_ quickly.

###  6.3 Solid/fluid interfaces

Let the fluid be bounded below by a flat horizontal wall (substrate) located
at $z=0$, so that boundary condition (35) reduces to

$\rho_{i}=\rho_{0,i}\qquad\text{at}\qquad z=0.$ (79)

Far above the substrate, the fluid is homogeneous,

$\rho_{i}\rightarrow\rho_{\infty,i}\qquad\text{as}\qquad z\rightarrow+\infty,$
(80)

where $\rho_{\infty,i}$ corresponds to a stable or metastable state. Finally,
equation (75) used in the previous subsection for 1D drops and bubbles applies
to the present case as well.

Boundary-value problem (75), (79)–(80) was solved numerically for a van der
Waals pure fluid with

$T=0.2,\qquad\rho_{0,1}=0.3.$ (81)

Its typical solutions are illustrated in figure 6 ($\rho_{0,i}$ varies,
$\rho_{\infty,1}$ is fixed) and in figure 7 ($\rho_{0,i}$ is fixed,
$\rho_{\infty,1}$ varies). The following features can be observed.

Figure 6: Examples of solid/fluid interfaces for a pure van der Waals fluid
with parameters (81), for a given $\rho_{\infty,1}$: (a)
$\rho_{\infty,1}=\rho_{1}^{(v)}$, $\rho_{0,1}$ varies from $0$ to
$\rho_{1}^{(l)}$; (b) $\rho_{\infty,1}=\rho_{1}^{(l)}$, $\rho_{0,1}$ varies
from $1$ to $\rho_{1}^{(v)}$. The arrows show the direction towards the
boundary of the pool of allowable values of $\rho_{0,1}$. Figure 7: Examples
of solid/fluid interfaces for a pure van der Waals fluid with parameters (81).
The fluid at infinity is: (1) undersaturated vapour, (2) oversaturated vapour
with two possible solutions (the unstable one is shown in thin line), (3)
oversaturated vapour with a single solution, (4) undersaturated liquid with a
single solution, _etc._ Observe that there are no solutions between curves (3)
and (4).

(i) For a given $\rho_{\infty,1}$, there exists a certain pool of $\rho_{0,i}$
which can be ‘connected’ to this $\rho_{\infty,1}$. If, for example,
$\rho_{\infty,1}$ equals the saturated vapour density (the case illustrated in
figure 6a), than $\rho_{0,1}$ must be smaller than the saturated liquid
density (and _vice versa_ : if $\rho_{\infty,1}=\rho_{1}^{(l)}$, then
$\rho_{0,i}>\rho_{1}^{(v)}$ – see figure 6b). When $\rho_{\infty,1}$
approaches the pool’s boundary, a clearly visible liquid/vapour interface
emerges in the solution and is gradually moving away from the substrate.

(ii) It follows from the above that, for some pairs
$\left(\rho_{0,1},~{}\rho_{\infty,1}\right)$, no solution exists. As
illustrated in figure 7, some of such pairs involve unstable values of
$\rho_{\infty,1}$ (and, thus, are unimportant), but there are also ones with a
metastable $\rho_{\infty,1}$.

To clarify what happens in such cases, the evolution was simulated numerically
(using the model and method used in the previous two subsections). It has
turned out that, if boundary-value problem (75), (79)–(80) does not have a
solution, the substrate triggers off a spontaneous phase change. This result
can be readily interpreted physically: if a sufficiently _hydrophilic_
substrate (with a sufficiently small $\rho_{0,1}$) touches _oversaturated
vapour_ , it triggers off spontaneous condensation, whereas a sufficiently
_hydrophobic_ substrate touching _undersaturated liquid_ triggers off
evaporation. In either case, a shock wave of phase change propagates away from
the substrate, and the solid/fluid interface cannot be stationary.

(iii) For some pairs $\left(\rho_{0,1},~{}\rho_{\infty,1}\right)$, there are
two different solutions – a monotonic and non-monotonic ones (see figure 7).
It is unclear which of the two occurs in reality.

It has turned out that only monotonic solutions of boundary-value problem
(79)–(80), (75) maximise the entropy, whereas non-monotonic solutions do not
(the former claim is proved in the general case, but the latter, only for
$N=1$ – see Appendix C.3). That is, non-monotonic solutions cannot be ruled
out with certainty for $N\geq 2$ – yet the mere fact that the general
stability proof that works for monotonic solutions cannot be extended to non-
monotonic ones seems to resolve the dilemma in favour of the former.

(iv) Another physically important conclusion follows from the fact that
boundary-value problem (75), (79)–(80) has no more than one stable solution –
and, thus, does _not_ admit solutions describing a liquid layer of a finite
thickness on a substrate, with saturated vapour above it. All such layers,
regardless of their thickness, evaporate – and the Kelvin effect cannot be
responsible for this effect (because the liquid/vapour interface is flat). The
evaporation in this case can only be caused by long-range interaction between
the interface and the substrate.

(v) It is interesting to compare the interfacial profiles shown in figures 6–7
to those computed by Evans et al. (2017) via the density functional theory and
Monte Carlo method for a Lennard-Jones fluid bounded by a single wall or
contained between two parallel walls. The single-wall profiles of Evans et al.
(2017) (see their figure 3) are qualitatively similar to those computed in
this paper, but their two-wall profiles are riddled with short-scale
oscillations (see figures 6, 14, and 20). It is not clear whether the
oscillations are caused by ‘interference’ of the walls: on the one hand, the
distance between the walls exceeds the spatial scale of wall–molecule
interaction by a factor of $30$ (hence, the ‘interference’ should be weak) –
but, on the other hand, no oscillations occur when this distance is infinite.

Whatever the nature of the oscillations is, one should not expect the DIM to
describe them due to the omission of the higher-order derivatives of
$\rho_{i}$ when obtaining expression (31) for the van der Waals force.

##  7 Surface tension and contact angle

Consider a fluid bounded below by a horizontal substrate, and an oblique
liquid/vapour interface intersecting the substrate – see figure 8. Note that
this figure depicts to a hydrophilic substrate, such that the contact angle
$\theta$ is smaller than $\pi/2$.

Figure 8: A schematic illustrating boundary conditions (85)–(86) for a static
contact line. $\theta$ is the contact angle, $\phi^{(s/v)}(z)$ describes the
_s_ olid/_v_ apour interface, _etc._

If in equilibrium, the setting outlined is described by the two-dimensional
reduction of equation (45). Given expression (46) for $\eta_{i}$, one obtains

$G_{i}-\sum_{j}K_{ij}\left(\frac{\partial^{2}\rho_{j}}{\partial
x^{2}}+\frac{\partial^{2}\rho_{j}}{\partial z^{2}}\right)+G_{\infty,i}=0.$
(82)

Impose also the boundary condition

$\rho_{i}\rightarrow\rho_{i}^{(v)}\qquad\text{as}\qquad z\rightarrow+\infty,$
(83)

which implies that the constant in (82) is
$G_{\infty,i}=G_{i}(\rho_{1}^{(v)}...\rho_{N}^{(v)},T)$. At the substrate, the
standard boundary condition,

$\rho_{i}=\rho_{0,i}\qquad\text{at}\qquad z=0,$ (84)

is assumed.

To close boundary-value problem (82)–(84), one should set boundary conditions
as $x\rightarrow\pm\infty$. The setting depicted in figure 8 corresponds to

$\displaystyle\rho_{i}$ $\displaystyle\rightarrow\rho_{i}^{(s/v)}(z)\hskip
66.01059pt\text{as}\qquad x\rightarrow-\infty,$ (85) $\displaystyle\rho_{i}$
$\displaystyle\rightarrow\rho_{i}^{(s/l)}(z)+\rho_{i}^{(l/v)}(\xi)\qquad\text{as}\qquad
x\rightarrow+\infty.$ (86)

Here, the function $\rho_{i}^{(s/v)}(z)$ describes a solid/vapour interface
(i.e., satisfies boundary-value problem (75), (79)–(80) with
$\rho_{\infty,i}=\rho_{i}^{(v)}$); the function $\rho_{i}^{(s/l)}(z)$
describes a solid/liquid interface; and $\rho_{i}^{(l/v)}(\xi)$ (where
$\xi=z\cos\theta-x\sin\theta$) describes a liquid/vapour interface tilted at
an angle $\theta$ (it satisfies boundary-value problem (68)–(70) with $z$
changed to $\xi$).

In the framework of the DIM, the contact angle $\theta$ is not arbitrary, but
is fully determined by the fluid’s thermodynamic properties, the Korteweg
matrix $K_{ij}$, and the near-wall density $\rho_{0,i}$. To derive an
expression for $\theta$ (similar to that derived by Pismen & Pomeau (2000) for
pure fluids), consider

$\sum_{i}\int_{0}^{\infty}\left(\ref{7.1}\right)\times\frac{\partial\rho_{i}}{\partial
x}\mathrm{d}z.$

After straightforward algebra involving integration by parts and the use of
boundary conditions (83)–(84) and identity (7), one obtains

$\frac{\partial}{\partial
x}\int_{0}^{\infty}\left[\sum_{j}\rho_{j}\left(G_{j}-G_{j,\infty}\right)-p+p_{\infty}+\frac{1}{2}\sum_{ij}K_{ij}\left(\frac{\partial\rho_{i}}{\partial
z}\frac{\partial\rho_{j}}{\partial z}-\frac{\partial\rho_{i}}{\partial
x}\frac{\partial\rho_{j}}{\partial x}\right)\right]\mathrm{d}z=0,$

where $p_{\infty}=p(\rho_{1}^{(v)}...\rho_{N}^{(v)},T)$. Finally, integrating
the above expression from $x=-\infty$ to $x=+\infty$ and recalling boundary
conditions (85)–(86), one obtains

$\int_{0}^{\infty}\left\\{\left[\sum_{i}\rho_{i}\left(G_{i}-G_{\infty,i}\right)-p+p_{\infty}\right]_{\rho_{i}=\rho_{i}^{(s/v)}}+\frac{1}{2}\sum_{ij}K_{ij}\frac{\mathrm{d}\rho_{i}^{(s/v)}}{\mathrm{d}z}\frac{\mathrm{d}\rho_{j}^{(s/v)}}{\mathrm{d}z}\right\\}\mathrm{d}z\\\
=\int_{0}^{\infty}\left\\{\left[\sum_{i}\rho_{i}\left(G_{i}-G_{\infty,i}\right)-p+p_{\infty}\right]_{\rho_{i}=\rho_{i}^{(s/l)}}+\frac{1}{2}\sum_{ij}K_{ij}\frac{\mathrm{d}\rho_{i}^{(s/l)}}{\mathrm{d}z}\frac{\mathrm{d}\rho_{j}^{(s/l)}}{\mathrm{d}z}\right\\}\mathrm{d}z\\\
+\int_{0}^{\infty}\left\\{\left[\sum_{i}\rho_{i}\left(G_{i}-G_{\infty,i}\right)-p+p_{\infty}\right]_{\rho_{i}=\rho_{i}^{(l/v)}}\vphantom{\left[\sum_{i}\right]_{\rho_{i}=\rho_{i}^{(l/v)}}\frac{\mathrm{d}\rho_{i}^{(l/v)}}{\mathrm{d}\xi}}\right.\\\
+\left.\frac{1}{2}\left(\cos^{2}\theta-\sin^{2}\theta\right)\sum_{ij}K_{ij}\frac{\mathrm{d}\rho_{i}^{(l/v)}}{\mathrm{d}\xi}\frac{\mathrm{d}\rho_{i}^{(l/v)}}{\mathrm{d}\xi}\vphantom{\left[\sum_{i}\right]_{\rho_{i}=\rho_{i}^{(l/v)}}\frac{\mathrm{d}\rho_{i}^{(l/v)}}{\mathrm{d}\xi}}\right\\}\frac{\mathrm{d}\xi}{\cos\theta}.$

This equality can be simplified using identity (72) (with $\rho_{i}$ changed
to$\rho_{i}^{(l/v)}$), and similar identities for $\rho_{i}^{(s/v)}$ and
$\rho_{i}^{(s/l)}$. After straightforward algebra involving re-denoting
$\xi\rightarrow z$, one obtains the DIM version of Young’s formula,

$\cos\theta=\frac{\sigma^{(s/v)}-\sigma^{(s/l)}}{\sigma^{(l/v)}},$ (87)

where

$\sigma^{(s/v)}=\sum_{ij}K_{ij}{\displaystyle\int_{0}^{\infty}}\dfrac{\mathrm{d}\rho_{i}^{(s/v)}}{\mathrm{d}z}\dfrac{\mathrm{d}\rho_{j}^{(s/v)}}{\mathrm{d}z}\mathrm{d}z,\qquad\sigma^{(s/l)}=\sum_{ij}K_{ij}{\displaystyle\int_{0}^{\infty}}\dfrac{\mathrm{d}\rho_{i}^{(s/l)}}{\mathrm{d}z}\dfrac{\mathrm{d}\rho_{j}^{(s/l)}}{\mathrm{d}z}\mathrm{d}z,$
$\sigma^{(l/v)}=\sum_{ij}K_{ij}{\displaystyle\int_{-\infty}^{\infty}}\dfrac{\mathrm{d}\rho_{i}^{(l/v)}}{\mathrm{d}z}\dfrac{\mathrm{d}\rho_{j}^{(l/v)}}{\mathrm{d}z}\mathrm{d}z$
(88)

are the surface tension coefficients of the solid/vapour, solid/liquid, and
liquid/vapour interfaces, respectively. Treating hydrophobic substrates
($\theta>\frac{1}{2}\pi$) in a similar fashion, one can show that formula (87)
applies to that case as well.

To calculate the surface tension – say, $\sigma^{(s/l)}$ – one first needs to
solve boundary-value problem (68)–(71) which determines the function
$\rho_{i}=\rho_{i}^{(l/v)}(z)$. For a pure fluid, however, the expression for
$\sigma^{(s/l)}$ can be rewritten as a closed-form integral. To do so, let
$N=1$ in equation (72), which yields

$\frac{\mathrm{d}\rho_{1}}{\mathrm{d}z}=-\sqrt{\frac{2}{K_{11}}\left\\{\rho_{1}\left[G_{1}(\rho_{1},T)-G_{1}(\rho_{1}^{(v)},T)\right]-p(\rho_{1},T)+p(\rho_{1}^{(v)},T)\right\\}}.$

Using this equation to change the variable of integration in the expression
for $\sigma^{(l/v)}$ in (88), one obtains

$\sigma^{(l/v)}=\sqrt{2K_{11}}{\displaystyle\int_{\rho_{1}^{(v)}}^{\rho_{1}^{(l)}}}\sqrt{\rho_{1}\left[G_{1}(\rho_{1},T)-G_{1}(\rho_{1}^{(v)},T)\right]-p(\rho_{1},T)+p(\rho_{1}^{(v)},T)}\mathrm{d}\rho.$
(89)

Thus, to calculate $\sigma^{(l/v)}$, one should first use the given
$p(\rho_{1},T)$ and $G_{1}(\rho_{1},T)$ to find $\rho_{1}^{(l)}$ and
$\rho_{1}^{(v)}$ (via the Maxwell construction) and then evaluate integral
(89).

##  8 Parameterizing the diffuse–interface model for water/air interfaces

To use the DIM in applications, one needs the following external parameters:
the fluid’s thermodynamic properties, the dependence of the viscosity and
transport coefficients on $\left(\rho_{1}...\rho_{N},T\right)$, and the
Korteweg matrix. In this section, an approach to specifying these parameters
is described and applied to water/air interfaces.

For reasons described in subsection 2.2, the fluid’s thermodynamic properties
will be approximated by the Enskog–Vlasov (EV) model. It involves an
undetermined matrix $a_{ij}$ and an undetermined function
$\Theta(\rho_{1}...\rho_{N})$, which should be fixed as the best fits of the
empiric characteristics of the fluid under consideration.

In subsections 8.1–8.2, it will be explained how the fitting should be carried
out for a pure fluid (in application to water, nitrogen, and oxygen). An
approach to determining $K_{11}$ for a pure fluid will be outlined in
subsection 8.3. Subsection 8.4 describes how $a_{ij}$, $\Theta$, and $K_{ij}$
can be determined for water–air mixture, and its viscosity and transport
coefficients are dealt with in subsection 8.5.

All these tasks will be carried out in terms of the original (dimensional)
variables.

###  8.1 The van der Waals parameter of pure water, nitrogen, and oxygen

For a pure fluid, the EV expression (12) for the internal energy yields (the
subscript 1 is omitted)

$e=cT-a\rho,$

which suggests that the van der Waals parameter $a$ can be determined as the
_linear_ fit of the empiric dependence of $cT-e$ on $\rho$. The heat capacity
$c$ in this expression should be the same as that in the Enskog–Vlasov kinetic
theory – i.e., $3R$ for water and $5R/2$ for nitrogen and oxygen. For
simplicity, the fitting was carried out using only the data on the critical
isobar $p=p_{cr}$ [as done by Benilov & Benilov (2019)], but the resulting
straight line fits the isobars $p=p_{cr}/2$ and $p=2p_{cr}$ reasonably well
too (see the top panels of figure 9). The parameters of the critical points
for the fluids under consideration, as well as the other parameters needed
here and hereinafter, can be found in Tables 1–2.

Figure 9: Thermodynamic properties of $\mathrm{H}_{2}\mathrm{O}$ (left-hand panels) and $\mathrm{N}_{2}$ (right-hand panels). The non-connected symbols show the empiric data from Lindstrom & Mallard (1997) presented in isobaric form, for three values of the pressure $p$ relative to its critical value, $p_{cr}$ (see the legend). The gap between $\rho^{(v)}$ and $\rho^{(l)}$ in the empiric data for $p=p_{cr}/2$ reflects the impossibility (difficulty) of measurements in an unstable (metastable) fluid. The dotted lines show the parameters calculated via the EV fluid model. Fluid |  $T_{tp}~{}(\mathrm{K})$ |  $\rho_{tp}^{(l)}~{}(\mathrm{kg\,m}^{-3})$ |  $\rho_{tr}^{(v)}~{}(\mathrm{kg\,m}^{-3})$ |  $T_{cr}~{}(\mathrm{K})$ |  $\rho_{cr}~{}(\mathrm{kg\,m}^{-3})$ |  $p_{cr}~{}(\mathrm{bar})$
---|---|---|---|---|---|---
$\mathrm{H}_{2}\mathrm{O}$ |  $273.160$ |  $~{}999.79$ |  $0.0048546$ |  $647.096$ |  $322.00$ |  $220.64$
$\mathrm{N}_{2}$ |  $63.151$ |  $~{}867.22$ |  $0.6742700$ |  $126.192$ |  $313.30$ |  $33.958$
$\mathrm{O}_{2}$ |  $54.361$ |  $1306.10$ |  $0.0103580$ |  $154.581$ |  $436.14$ |  $50.430$
Table 1: The parameters of the triple and critical points (subscripts tr and cr, respectively) of $\mathrm{H}_{2}\mathrm{O}$, $\mathrm{N}_{2}$, and $\mathrm{O}_{2}$ Lindstrom & Mallard (1997). Fluid |  $m~{}(10^{-26}\mathrm{kg})$ |  $R~{}(\mathrm{m}^{2}\mathrm{s}^{-2}\mathrm{K}^{-1})$ |  $a~{}(\mathrm{m}^{5}\mathrm{s}^{-2}\mathrm{kg}^{-1})$ |  $K~{}(10^{-17}\mathrm{m}^{7}\mathrm{s}^{-2}\mathrm{kg}^{-1})$
---|---|---|---|---
$\mathrm{H}_{2}\mathrm{O}$ |  $2.9915$ |  $461.53$ |  $2112.1$ |  $1.87810$
$\mathrm{N}_{2}$ |  $4.6516$ |  $296.81$ |  $222.23$ |  $1.50780$
$\mathrm{O}_{2}$ |  $5.3135$ |  $259.84$ |  $172.73$ |  $0.84587$
air |  $4.7706$ |  $289.41$ |  $211.84$ |  $1.36880$
Table 2: The parameters of $\mathrm{H}_{2}\mathrm{O}$, $\mathrm{N}_{2}$,
$\mathrm{O}_{2}$, and air: $m$ is the molecular mass, $R$ is the specific gas
constant, $a$ is the van der Waals parameter, $K$ is the Korteweg constant.
The parameters of air are calculated as the 79-21 weighted averages of the
corresponding parameters of nitrogen and oxygen, respectively.

The resulting values of the van der Waals parameter $a$ for
$\mathrm{H}_{2}\mathrm{O}$, $\mathrm{N}_{2}$, and $\mathrm{O}_{2}$ are
presented in Table 2. Interestingly, $a$ of water exceeds significantly those
of nitrogen and oxygen. This is likely to be caused by the difference in the
molecular structure of these fluids: the water molecule has a non-zero dipolar
moment, whereas nitrogen and oxygen molecules are symmetric – hence, do not.
Since the van der Waals force is of electric nature, one can assume that
dipolar molecules interact stronger than neutral ones.

###  8.2 The equation of state of pure $\mathrm{H}_{2}\mathrm{O}$,
$\mathrm{N}_{2}$, and $\mathrm{O}_{2}$

Pure fluids will be described by the EV model (12)–(13) with

$\Theta=R\left[-q^{(0)}\ln\left(1-0.99\frac{\rho}{\rho_{tp}}\right)+\sum_{n=1}^{4}q^{(n)}\left(\frac{\rho}{\rho_{tp}}\right)^{n}\right],$
(90)

where $q^{(0)}$… $q^{(4)}$ are undetermined coefficients and $\rho_{tp}$ is
the fluid’s density at the triple point ($\rho_{tp}$ is just a convenient
density scale; the fact that, at the triple point, the liquid and vapour are
in equilibrium with the solid phase is irrelevant, as solids are not described
by the DIM). Note that the first term in expression (90) sets the maximum
density at $\rho_{tp}/0.99$.

The coefficients $q^{(n)}$ were determined for $\mathrm{H}_{2}\mathrm{O}$,
$\mathrm{N}_{2}$, and $\mathrm{O}_{2}$ by ensuring that the expressions for
$p(\rho,T)$ and $G(\rho,T)$ corresponding to (90) yield the ‘correct’ – i.e.,
empiric – values for the critical density, temperature, and pressure, as well
as the liquid and vapour densities at the triple-point temperature $T=T_{tp}$
(five equations for the five unknown coefficients). All the necessary empiric
data can be found in Tables 1–2, and the computed values of $q^{(n)}$ are
listed in Table 3. Such an approach to calibrating the EV fluid model is a
slight modification of that of Benilov (2020b) – which is, in turn, a
modification of the approach of Benilov & Benilov (2018).

Fluid |  $q^{(0)}$ |  $q^{(1)}$ |  $q^{(2)}$ |  $q^{(3)}$ |  $q^{(4)}$
---|---|---|---|---|---
$\mathrm{H}_{2}\mathrm{O}$ |  $0.071894$ |  $1.4139$ |  $8.1126$ |  $-8.3669$ |  $4.0238$
$\mathrm{N}_{2}$ |  $-0.0013920$ |  $0.72934$ |  $6.4799$ |  $-8.1143$ |  $5.0186$
$\mathrm{O}_{2}$ |  $0.010770$ |  $0.58901$ |  $8.5357$ |  $-12.034$ |  $8.0872$
Table 3: The coefficients of the equation of state (90) for
$\mathrm{H}_{2}\mathrm{O}$, $\mathrm{N}_{2}$, and $\mathrm{O}_{2}$.

To illustrate how well the EV model describes real fluids, the isobaric (with
$p$ fixed) dependence of the temperature on the density is plotted for
$\mathrm{H}_{2}\mathrm{O}$ and $\mathrm{N}_{2}$ in the lower panels of figure
9. The results for $\mathrm{O}_{2}$ (not presented) are similar to those for
$\mathrm{N}_{2}$.

Note that the DIM has been coupled with realistic equations of state before
(e.g. Caupin, 2005; Magaletti et al., 2021); for pure water, this was
typically done using the IAPWS-95 equation (Wagner & Pruß, 2002). The EV model
used here is undoubtedly less accurate than the IAPWS-95 model – but it allows
one to describe consistently all of the fluids involved
($\mathrm{H}_{2}\mathrm{O}$, $\mathrm{N}_{2}$, and $\mathrm{O}_{2}$).

###  8.3 The Korteweg parameter of pure $\mathrm{H}_{2}\mathrm{O}$,
$\mathrm{N}_{2}$, and $\mathrm{O}_{2}$

(Benilov, 2020b) proposed to deduce $K$ from the requirement that the DIM
predict the correct value of the surface tension $\sigma^{(l/v)}$ of
liquid/vapour interface at the triple point. The same was done in the present
work: $\sigma^{(l/v)}$ was calculated via the DIM formula (89) with $G$ and
$p$ of the EV fluid, and the value of $K$ was chosen for which (89) yields the
same result as the empiric formula of Somayajulu (1988). The resulting $K$ for
$\mathrm{H}_{2}\mathrm{O}$, $\mathrm{N}_{2}$, and $\mathrm{O}_{2}$ can be
found in Table 2.

It should be emphasised that the values of $K$ computed via the above approach
depend on the chosen fluid model. This explains the difference between the
result for $K$ in this paper and the one computed by Benilov (2020b): the
former is based on the EV fluid model and the latter, on the van der Waals
model. The resulting 30% difference in $\sigma^{(l/v)}$ reflects the fact that
the latter model is much less accurate.

As a test of the whole approach based on the DIM and EV models, the
theoretical values for the saturation characteristics and $\sigma^{(l/v)}$
have been compared to their empiric counterparts for the whole temperature
range where liquid/vapour interfaces exist, i.e., between the triple and
critical points. The results are shown in figure 10: evidently, the
theoretical prediction are reasonably accurate.

Figure 10: Comparison of the results obtained via the DIM for an EV fluid
(solid line) with the corresponding empiric data (dashed line): (a) the
densities of saturated vapour and liquid, (b) saturated pressure, (c) surface
tension of liquid/vapour interface. The empiric data in panels (a) and (b) are
from (Lindstrom & Mallard, 1997), and those in panel (c) are from (Somayajulu,
1988).

###  8.4 Parameters of water–air interaction

Generally, the parameters of a fluid should be deduced from measurements of
its equation of state, surface tension, _etc._ – but these are rarely known
for _multicomponent_ fluids. The water–air mixture _at normal conditions_ is
an exception in this respect: the density of air is small in this case, and
its equation of state is close to that of an ideal gas. In addition, air will
be treated as a pure fluid with parameters equal to the weighted averages of
those of nitrogen and oxygen (see Table 2).

Thus, let the function $\Theta$ (the non-ideal part of the entropy of an EV
fluid) be independent of the air density $\rho_{2}$, so that expressions
(15)–(16) with $N=2$ yield

$p=T\left[R_{1}\rho_{1}+R_{2}\rho_{2}+\left(\rho_{1}+\rho_{2}\right)\frac{\mathrm{d}\Theta(\rho_{1})}{\mathrm{d}\rho_{1}}\right]-a_{11}\rho_{1}^{2}-2a_{12}\rho_{1}\rho_{2},$
(91)
$G_{1}=T\left[R_{1}\ln\rho_{1}+\Theta(\rho_{1})+\left(\rho_{1}+\rho_{2}\right)\frac{\mathrm{d}\Theta(\rho_{1})}{\mathrm{d}\rho_{1}}\right]-2\left(a_{11}\rho_{1}+a_{12}\rho_{2}\right)+T\left[R_{1}+c_{1}\left(1-\ln
T\right)\right],$ (92)
$G_{2}=T\left[R_{2}\ln\rho_{2}+\Theta(\rho_{1})\right]-2a_{12}\rho_{1}+T\left[R_{2}+c_{2}\left(1-\ln
T\right)\right].$ (93)

Here, the function $\Theta(\rho_{1})$ is the one given by formula (90), with
$\rho$ changed to $\rho_{1}$, and with its coefficients corresponding to
water.

Before using expressions (91)–(93), one needs to fix the non-diagonal term
$a_{12}$ of the matrix $a_{ij}$, responsible for the interaction of water and
air molecules. It can be deduced from $\rho_{2}^{(l)}$ (the amount of air
dissolved in water): one should choose such $a_{12}$ that the Maxwell
construction based on (91)–(93) predicts the correct value of
$\rho_{2}^{(l)}$. Since $a_{12}$ is supposed to not depend on $T$, such a
calculation should be done for a single temperature and the atmospheric
pressure. For $T=25^{\circ}\mathrm{C}$ and $p=1\,\mathrm{atm}$, for example,
$\rho_{2}^{(l)}=0.0227\,\mathrm{kg\,m}^{-3}$ (Engineering ToolBox, 2004). It
can then be deduced that the Maxwell construction based on (91)–(93) yields
the correct value of $\rho_{2}^{(l)}$ if

$a_{12}=208\,\mathrm{m}^{5}\mathrm{s}^{-2}\mathrm{kg}^{-1},$

i.e., approximately equal to $a_{22}$ and ten times smaller than $a_{11}$.

The accuracy of expressions (91)–(93) can be illustrated by the corresponding
value of the boiling point ($T\approx 109^{\circ}\mathrm{C}$, as opposed to
the exact value of $T=100^{\circ}\mathrm{C}$). Furthermore, at the “room
temperature” (say, $T=25^{\circ}\mathrm{C}$), the error should be even
smaller, because the room temperature is close to the triple point of water
($T\approx 0^{\circ}\mathrm{C}$) where (91)–(93) were calibrated.

It still remains to determine the non-diagonal term $K_{12}$ of the Korteweg
matrix $K_{ij}$.

Since $K_{ij}$ is supposed to be positive definite, $K_{12}$ should satisfy

$-1<\frac{K_{12}}{\left(K_{11}K_{22}\right)^{1/2}}<1.$ (94)

One would think that $K_{12}$ could be determined by comparing the surface
tension of liquid-water/air interface to that of liquid-water/vapour _-water_
interface. It turns out, however, that the difference between these is so
small that the existing experimental techniques cannot detect it (see
Vargaftik et al., 1983, and references therein). This can be due to smallness
of $K_{12}$ – but more likely, because the density of air is small. Either
way, the determination of $K_{12}$ via measurements of surface tension is
impossible.

To at least confirm that air cannot affect significantly the surface tension
of water, boundary-value problem (68)–(71) was solved for the parameters of
water and air at $T=25^{\circ}\mathrm{C}$, and $K_{12}$ varying through range
(94). It has turned out that the resulting variation of the surface tension is
only $7\%$.

Even though the dependence of the surface tension (which is a _global_
characteristic) on $K_{12}$ is weak, this parameter does influence the _local_
structure of the interface.

Figure 11 depicts the density profiles $\rho_{1}(z)$ and $\rho_{2}(z)$ for the
water/air interface at room temperature [i.e., the above-mentioned solution of
boundary-value problem (68)–(71)]. One can see that, as $K_{12}$ approaches
the right endpoint of range (94), the profile of water density is getting
steeper, and a layer of high air density is developing inside the interface.
For curve 5 of figure 11(b) – which is half-way through range (94) – the
maximum of $\rho_{2}(z)$ exceeds the atmospheric air density by a factor of
approximately $20$. Since the amount of air drawn into the interface grows
with $K_{12}$, one concludes that the air is ‘pulled in’ by the van der Waals
attraction exerted by the bulk of the liquid.

Figure 11: The spatial structure of water/air interfaces, for
$T=25^{\circ}\mathrm{C}$ and $K_{12}=\left(K_{11}K_{22}\right)^{1/2}\times
n/10$, where $n$ is the curve number: (a) density of water; (b) density of
air. The horizontal scales of panel (a) and (b) differ by a factor of 5; the
latter and the horisontal scale of the shaded inset in panel (b) differ by a
factor of 1500. The inset shows that the concentration of the gas dissolved in
liquid is non-zero.

This suggests a possibility of deducing $K_{12}$ from empiric data on
evaporation of drops (which depends on both global and local characteristics
of the interface (Benilov, 2022b)). Such an approach, however, requires an
extension of the evaporation model of Benilov (2022b) to multicomponent
fluids, which is currently in progress.

Observe that the interfaces depicted in figure 11 are such that $\rho_{2}(z)$
is not monotonic – hence, the sufficient stability criterion from subsection
6.1 and Appendix C.1 does not apply. This does not, however, mean that this
interface is unstable. It is, in fact, difficult to imagine that a microscopic
non-monotonicity of the density of air dissolved in water could destabilise
the whole interface – but this issue is still of interest theoretically and,
thus, deserves further investigation.

###  8.5 The viscosity and thermal conductivity of water–air mixture

There is a large body of work on the viscosity of multicomponent fluids, with
even the simplest theoretical models – e.g., that of Enskog–Chapman for dilute
gases (e.g., Ferziger & Kaper, 1972) – yielding a fairly complicated
dependence of $\mu_{s}$ and $\mu_{b}$ on $\rho_{i}$. Phenomenological results,
on the other hand, tend to include many _ad hoc_ parameters specific to the
fluid under consideration (e.g. Davidson, 1993, and references therein).
Overall, the simplest option seems to be the expression for the shear
viscosity proposed on a phenomenological basis by Hind et al. (1960) and
justified, under certain approximations, via statistical mechanics by Bearman
& Jones (1960),

$\mu_{s}=\mu_{s,1}f_{1}^{2}+\left(\mu_{s,1}+\mu_{s,2}\right)f_{1}f_{2}+\mu_{s,2}f_{2}^{2},$
(95)

where $\mu_{s,i}$ is the shear viscosity of the $i$-th species and

$f_{i}=\frac{\rho_{i}/m_{i}}{\rho_{1}/m_{1}+\rho_{2}/m_{2}}$

is its mole fraction. Expression (95) does not include any fluid-specific
parameters, but is capable of predicting the shear viscosity with a reasonable
accuracy. If, for example, air is treated as a mixture of $\mathrm{N}_{2}$ and
$\mathrm{O}_{2}$, the error of (95) for the range $T=0-100^{\circ}\mathrm{C}$
is less than 1%. A similar “mixture rule” can be assumed for the bulk
viscosity and thermal conductivity,

$\mu_{b}=\mu_{b,1}f_{1}^{2}+\left(\mu_{b,1}+\mu_{b,2}\right)f_{1}f_{2}+\mu_{b,2}f_{2}^{2},$
(96)
$\kappa=\kappa_{1}f_{1}^{2}+\left(\kappa_{1}+\kappa_{2}\right)f_{1}f_{2}+\kappa_{2}f_{2}^{2}.$
(97)

Generally, various mechanical properties of a multicomponent fluid are often
described by the same mixture rule, in which case it is referred to as
“generalised”.

It still remains to fix the dependence of $\mu_{s,i}$, $\mu_{b,i}$, and
$\kappa_{i}$ on $\left(\rho_{1},\rho_{2},T\right)$. In application to air –
which can be treated as a dilute gas – one can assume these parameter to
depend only on $T$ (as predicted by the kinetic theory of dilute gases), i.e.,

$\mu_{s,2}(\rho_{2},T)=\mu_{s,2}(0,T),\qquad\mu_{b,2}(\rho_{2},T)=\mu_{b,2}(0,T),\qquad\kappa_{2}(\rho_{2},T)=\kappa_{2}(0,T),$
(98)

where $\mu_{s,2}(0,T)$, $\mu_{b,2}(0,T)$, and $\kappa_{2}(0,T)$ are the small-
density limiting values of the corresponding parameters.

For water, whose liquid phase cannot be treated as a dilute gas, such an
approximation is inapplicable. Aiming again for simplicity, one can
approximate both viscosities by a parabola passing through two reference
points – the zero-density limit and the saturated liquid state,

$\mu_{s,1}(\rho_{1},T)=\mu_{s,1}(0,T)+\left[\mu_{s,1}(\rho_{1}^{(l)},T)-\mu_{s,1}(0,T)\right]\left[\frac{\rho_{1}}{\rho_{1}^{(l)}(T)}\right]^{2},$
(99)
$\mu_{b,1}(\rho_{1},T)=\mu_{b,1}(0,T)+\left[\mu_{b,1}(\rho_{1}^{(l)},T)-\mu_{b,1}(0,T)\right]\left[\frac{\rho_{1}}{\rho_{1}^{(l)}(T)}\right]^{2}.$
(100)

A similar approximation, but by a linear function, will be adopted for the
thermal conductivity,

$\kappa_{1}(\rho_{1},T)=\kappa_{1}(0,T)+\left[\kappa_{1}(\rho_{1}^{(l)},T)-\kappa_{1}(0,T)\right]\frac{\rho_{1}}{\rho_{1}^{(l)}(T)}.$
(101)

The accuracy of the above approximations for $\mu_{s,1}$ and $\kappa_{1}$ is
illustrated in figure 12.

Figure 12: The parameters of pure water at $T=327^{\circ}\mathrm{C}$: (a)
shear viscosity, (b) thermal conductivity. The empiric data (Lindstrom &
Mallard, 1997) are shown in solid line and approximations (99) and (101), in
dashed line. The dots mark the reference points. The gap in the empiric data
reflects the impossibility (difficulty) of measurements in an unstable
(metastable) fluid.

To use formulae (95)–(101), one needs to know how the viscosities and thermal
conductivities of water and air depend on $T$. For $\mu_{s}$ and $\kappa$,
such data are widely available (e.g., Lindstrom & Mallard, 1997; White, 2005).
Measurements of $\mu_{b}$, on the other hand, are scarce, but can still be
found in (Holmes et al., 2011), (Cramer, 2012), and (Shang et al., 2019) for
liquid water, vapour water, and air, respectively.

The empiric formulae proposed in these papers will not be discussed in detail.
The characteristic values they yield for normal conditions are listed in Table
4.

Fluid |  $\mu_{s}~{}(\mathrm{Pa\,s})$ |  $\mu_{b}~{}(\mathrm{Pa\,s})$ |  $\kappa~{}(\mathrm{W\,K}^{-1}\mathrm{m}^{-1})$
---|---|---|---
liquid water |  $0.890\times 10^{-3}$(Lindstrom & Mallard, 1997) |  $2.47\times 10^{-3}$(Holmes et al., 2011) |  $0.606460$(Lindstrom & Mallard, 1997)
vapour water |  $0.971\times 10^{-5}$(Lindstrom & Mallard, 1997) |  $7.20\times 10^{-5}$(Cramer, 2012) |  $0.018433$(Lindstrom & Mallard, 1997)
air |  $1.840\times 10^{-5}$(White, 2005) |  $1.75\times 10^{-5}$(Shang et al., 2019) |  $0.026089$ (White, 2005)
Table 4: The empiric viscosity and thermal conductivity of liquid water, water
vapour, and air at $25^{\circ}\mathrm{C}$, and the corresponding references.

###  8.6 The transport coefficients

Recalling restrictions (25), one can deduce that the extended transport matrix
for a two-component fluid is

$D_{ij}^{(ext)}=\begin{bmatrix}D_{11}&-D_{11}&\zeta_{1}\\\
-D_{11}&D_{11}&-\zeta_{1}\\\ \zeta_{1}&-\zeta_{1}&\kappa\end{bmatrix}.$

Evidently, it involves only three independent coefficients, one of which (the
thermal conductivity $\kappa$) has already been discussed in subsection 8.5.
The other two, $D_{11}$ and $\zeta_{1}$, are discussed below.

It can be safely assumed that diffusion is of importance only where the water
density is comparable to that of air. In the region where the former is high,
the latter is miniscule, and so is its influence on the global dynamics. This
effectively means that the diffusivities can be represented using the
Enskog–Chapman method: even though it does not apply to liquid water, the
error due to using it anyway is negligible.

According to the leading-order Enskog–Chapman formula (e.g., Ferziger & Kaper,
1972),

$D_{11}=\frac{\rho_{1}\rho_{2}\left(\rho_{1}m_{2}+\rho_{2}m_{1}\right)}{\left(\rho_{1}+\rho_{2}\right)^{2}k_{B}T}\mathfrak{D}(T),$
(102)

where $\mathfrak{D}(T)$ does not depend on $\rho_{1}$ and $\rho_{2}$. Its
dependence on $T$ can be found in (Engineering ToolBox, 2018), from which the
following characteristic value can be deduced,

$\mathfrak{D}(25^{\circ}\mathrm{C})=2.49\times
10^{-5}\mathrm{m}^{2}\mathrm{s}^{-1}.$ (103)

The Enskog–Chapman expression for the thermodiffusivity is fairly bulky – as a
result, one often uses the simpler formula of de Groot & Mazur (1962),

$\zeta_{1}=\left(\rho_{1}+\rho_{2}\right)\frac{\left(\rho_{1}/m_{1}\right)\left(\rho_{2}/m_{2}\right)}{\left(\rho_{1}/m_{1}+\rho_{2}/m_{2}\right)^{2}}\mathfrak{U}(T).$
(104)

Unfortunately, no data on $\mathfrak{U}(T)$ for water–air mixture are
available in the literature; the author of this paper was able to find only an
estimate for a single temperature value (Lidon et al., 2021),

$\mathfrak{U}(21^{\circ}\mathrm{C})=-4.98\times
10^{-6}\,\mathrm{m}^{2}\mathrm{s}^{-1}.$ (105)

It is shown later, however, that thermodiffusion in water-air interfaces is
negligible, so the lack of data on $\mathfrak{U}(T)$ is unimportant.

In the next subsection, characteristic values of the coefficients $D_{11}$ and
$\zeta_{1}$ will be needed. These will be estimated for the particular case
$\rho_{1}=\rho_{2}=1.17\,\mathrm{kg\,m}^{-3}$ – i.e., when the water density
matches that of air at normal conditions. Substituting these values and the
molecular masses of water and air from Table 2 into formulae (102)–(105), one
obtains

$D_{11}\approx 1.37\times
10^{-10}\mathrm{m}^{3}\mathrm{s\,kg},\qquad\zeta_{1}\approx-2.76\times
10^{-6}\mathrm{m}^{-1}\mathrm{s}^{-1}\mathrm{kg}.$ (106)

Even though these values apply to different temperatures
($T=25^{\circ}\mathrm{C}$ and $T=21^{\circ}\mathrm{C}$), they will be used
simultaneously in the same qualitative estimate (under the implied assumption
that the four-degree difference does not matter).

###  8.7 The nondimensional parameters

To understand which effects are important for water/air interfaces under
normal conditions, one should estimate the nondimensional parameters (55).
They will be calculated using the characteristics of water from Tables 2–4 and
estimates (106). The characteristic pressure scale will be assumed to be
$P=a\varrho^{2}$ where $\varrho=997\,\mathrm{kg\,m}^{-3}$ is the density of
liquid water at $25^{\circ}\mathrm{C}$ and $a$ equals the van der Waals
parameter of water from Table 2. The following expression for the viscosity
scale will be used:

$\mu=\frac{4}{3}\mu_{s}+\mu_{b},$

which arises naturally in problems involving liquid films (Benilov, 2020c,
2022a) and drops (Benilov, 2022b). In the estimates for this paper, the
viscosities $\mu_{s}$ and $\mu_{b}$ are those for water at
$25^{\circ}\mathrm{C}$.

The following values of parameters (55) have been obtained:

$\alpha\approx 0.00139,\qquad\tau\approx 0.0653,\qquad\beta\approx
0.0593,\qquad\nu\approx 0.0321,\qquad\delta\approx 0.0176.$ (107)

The smallness of the microscopic Reynolds number $\alpha$ suggests that
inertia plays little role in the dynamics of water/air interfaces (hence, one
can use the Stokes approximation). The smallness of the nondimensional
temperature $\tau$ does not have physical implications, but enables one to use
asymptotic tools when calculating, say, the profile of the equilibrium
interface (Benilov, 2020c). The smallness of the Nusselt number $\nu$ implies
that the Soret and Dufour effects are negligible, and so the lack of empiric
data on them for water and air is unimportant. The smallness of $\beta$, in
turn, implies that the flow is almost isothermal, whereas the small value of
$\delta$ suggests that diffusion dominates advection.

Similar estimates have also been carried for the parameters of air at
$p=1\,\mathrm{atm}$ and $T=25^{\circ}\mathrm{C}$. It has turned out that
$\alpha$, $\beta$, $\nu$, and $\delta$ are even smaller than those for water,
but the nondimensional temperature $\tau$ for air is order-one. The latter is
not surprising, as the room temperature is much higher than the freezing
temperatures of nitrogen and oxygen, but is close to that for water.

All these observations should be helpful when using the DIM to examine
asymptotically the problems of evaporation of drops and moving contact lines.

One should keep in mind, however, that parameters (55) are sensitive to the
choice of the viscosity scale $\mu$. Benilov (2020a), for example, chose $\mu$
equal the half-sum of the shear viscosities of liquid water and vapour water –
as a result, $\alpha$ and $\beta$ were noticeably larger than those calculated
above. More generally, one should choose the viscosity scale, as well as the
other characteristic scales, to reflect the essentially physics of the setting
at hand.

Note also that estimates (107) apply to water under _normal_ conditions – but
in industrial applications, the governing nondimensional parameters can be
very different. In steam turbines, for example, the temperature can be as high
as $540^{\circ}\mathrm{C}$, and the pressure can exceed $230\,\mathrm{atm}$
(Vasserman & Shutenko, 2017).

##  9 Concluding remarks

The results obtained in this paper can be briefly summarised as follows:

  1. 1.

The entropy principle and conservation laws of the multicomponent diffuse-
interface model have been used to examine the stability of liquid/vapour
interfaces. Several physically important results are reported:

     * •

Flat liquid/vapour interfaces in an unbounded fluid are stable if the density
profiles of all species are monotonic.

     * •

There can exist up to two different solutions describing a solid/vapour
interface, one monotonic and one non-monotonic. The former is stable and the
latter is likely to be unstable (_definitely_ unstable for pure fluids).
Similar conclusions apply to solid/liquid interfaces.

     * •

For certain values of the near-substrate density (which is an external
parameter in DIM, linked to the contact angle of the substrate), no steady
solution exists describing a solid/vapour interfaces. Physically, such
substrates are too hydrophilic and, thus, trigger off spontaneous condensation
of adjacent vapour. Similarly, some substrates are too hydrophobic and trigger
off spontaneous evaporation of adjacent liquid.

     * •

A liquid layer between a flat substrate and a semi-space filled with vapour
can exist as a steady state only if the vapour is _oversaturated_. All such
layers are unstable, however: depending on the perturbation, they either fully
evaporate or grow indefinitely due to condensation of vapour on its upper
boundary.
If the vapour above the layer is _saturated_ or _undersaturated_ , the liquid
evaporates and no steady solution exists.
Similar conclusions apply to 1D vapour layers between a flat substrate and a
semi-space filled with liquid.

     * •

Similar conclusions to those in the previous bullet apply to a liquid layer
with vapour both above and below it, and a vapour layer with liquid below and
above it. Such solutions can be viewed as 1D drops and bubbles, respectively.

  2. 2.

The multicomponent DIM has been fully parameterised for water/air interfaces
under normal conditions. It is shown that the Soret and Dufour effects are
weak in this case, which agrees with the results of Jiang et al. (2020). It is
also argued that the interfacial flow in this case is isothermal.
These are of importance when studying evaporation of water drops or the
dynamics of their contact lines.

## Declaration of interests

The author reports no conflict of interests.

##  Appendix A The Gibbs relation

Introduce the volumetric densities of energy
$\mathcal{U}(\rho_{1}...\rho_{N},T)$, entropy
$\mathcal{S}(\rho_{1}...\rho_{N},T)$, and partial chemical potentials
$\mathcal{G}_{i}(\rho_{1}...\rho_{N},T)$ – related to the corresponding
specific quantities by

$\mathcal{U}=e\rho,\qquad\mathcal{S}=s\rho,\qquad\mathcal{G}_{i}=G_{i}.$ (108)

Now, definition (3) of $G_{i}$ can be rewritten in the form

$\mathcal{G}_{i}=\frac{\partial\mathcal{U}}{\partial\rho_{i}}-T\frac{\partial\mathcal{S}}{\partial\rho_{i}}.$
(109)

The standard volumetric version of the Gibbs relation (e.g., Eq. (2.4) of
Giovangigli & Matuszewski, 2013) is

$\mathrm{d}\mathcal{U}=T\,\mathrm{d}\mathcal{S}+\sum_{i}\mathcal{G}_{i}\,\mathrm{d}\rho_{i}.$

Substituting (108)–(109) into this equality and rewriting it in terms of
$\mathrm{d}T$ and $\mathrm{d}\rho_{i}$ (instead of $\mathrm{d}\mathcal{S}$ and
$\mathrm{d}\rho_{i}$), one obtains (4) as required.

##  Appendix B The general boundary condition at a solid wall

Consider the following generalization of boundary condition (35):

$\sum_{j}C_{ij}\left(\mathbf{n}\cdot\mathbf{\nabla}\rho_{j}\right)+\rho_{i}=\rho_{0,i}\qquad\text{at}\qquad\mathbf{r}\in\partial\mathcal{D},$
(110)

where $C_{ij}$ is a symmetric matrix, depending generally on
$\rho_{1}...\rho_{N}$. (110) is a multicomponent extension of a boundary
condition often used for pure fluids (e.g. Madruga & Thiele, 2009; Gallo et
al., 2021).

To understand how the energy conservation law is affected by switching to a
different boundary condition, observe that the governing equations (17)–(22)
and the other boundary conditions, (32)–(34), imply that

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathcal{D}}\left[\rho
e+\frac{1}{2}\rho\left|\mathbf{v}\right|^{2}+\frac{1}{2}\sum_{ij}K_{ij}\left(\mathbf{\nabla}\rho_{i}\right)\cdot\left(\mathbf{\nabla}\rho_{j}\right)\right]\mathrm{d}^{3}\mathbf{r}\\\
+\int_{\partial\mathcal{D}}\sum_{i}\left(\mathbf{n}\cdot\mathbf{\nabla}\rho_{i}\right)\sum_{j}K_{ij}\frac{\partial\rho_{j}}{\partial
t}\mathrm{d}A=0,$ (111)

where $\mathrm{d}A$ is the elemental aria on $\partial\mathcal{D}$. The ‘old’
boundary condition (35) implies that the second integral in this equality
vanishes, making the integral in the first term be a conserved quantity (the
energy) in this case.

Next, assume that the ‘new’ boundary condition (110) is imposed. Differentiate
it with respect to $t$, change the indices – first $j\rightarrow k$, then
$i\rightarrow j$ – and use the resulting equality to rearrange (111) in the
form

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathcal{D}}\left[\rho
e+\frac{1}{2}\rho\left|\mathbf{v}\right|^{2}+\frac{1}{2}\sum_{ij}K_{ij}\left(\mathbf{\nabla}\rho_{i}\right)\cdot\left(\mathbf{\nabla}\rho_{j}\right)\right]\mathrm{d}^{3}\mathbf{r}\\\
-\int_{\partial\mathcal{D}}\sum_{ij}K_{ij}^{\prime}\left(\mathbf{n}\cdot\mathbf{\nabla}\rho_{i}\right)\frac{\partial\left(\mathbf{n}\cdot\mathbf{\nabla}\rho_{j}\right)}{\partial
t}\mathrm{d}A=0,$ (112)

where

$K_{ij}^{\prime}=\sum_{k}K_{ik}C_{kj}.$

Equation (112) implies that the following quantity is conserved:

$E=\int_{\mathcal{D}}\left[\rho
e+\frac{1}{2}\rho\left|\mathbf{v}\right|^{2}+\frac{1}{2}\sum_{ij}K_{ij}\left(\mathbf{\nabla}\rho_{i}\right)\cdot\left(\mathbf{\nabla}\rho_{j}\right)\right]\mathrm{d}^{3}\mathbf{r}\\\
-\frac{1}{2}\int_{\partial\mathcal{D}}\sum_{ij}K_{ij}^{\prime}\,\left(\mathbf{n}\cdot\mathbf{\nabla}\rho_{i}\right)\left(\mathbf{n}\cdot\mathbf{\nabla}\rho_{j}\right)\mathrm{d}A.$

The second term in this expression is the surface energy corresponding to the
new boundary condition (110).

By comparison to the Dirichlet boundary condition (35), the mixed condition
(110) includes additional $N\left(N-1\right)/2$ adjustable parameters. Even
though the extra parameters may enable the DIM to be more accurate
quantitatively, one should think that a physics-based model does not need many
adjustable parameters to capture the qualitative nature of the phenomenon it
is applied to.

##  Appendix C Stability of one-dimensional steady states

In this appendix, the stability of three families of 1D solutions are
examined:

  1. 1.

liquid/vapour interfaces,

  2. 2.

1D drops and bubbles,

  3. 3.

solid/fluid interfaces.

###  C.1 Liquid/vapour interfaces

As mentioned before, equation (45) describes all steady states in the problem
at hand, and expression (47) describes the corresponding second variation of
the entropy. To adapt the latter for flat liquid/vapour interfaces, one needs
to replace the 3D integral over the domain $\mathcal{D}$ with a 1D integral
over the domain $-\infty<z<\infty$,

$\delta^{2}H=\frac{1}{T}\int_{-\infty}^{\infty}\sum_{ij}\left[-\frac{\partial
G_{i}}{\partial\rho_{j}}\left(\delta\rho_{i}\right)\left(\delta\rho_{j}\right)-K_{ij}\frac{\mathrm{d}\left(\delta\rho_{i}\right)}{\mathrm{d}z}\frac{\mathrm{d}\left(\delta\rho_{j}\right)}{\mathrm{d}z}\right]\mathrm{d}z.$

This expression can be rewritten in the form

$\delta^{2}H=\frac{1}{T}\int_{-\infty}^{\infty}\sum_{ij}\delta\rho_{i}\,\mathrm{\hat{O}}_{ij}\,\delta\rho_{j}\,\mathrm{d}z,$
(113)

where

$\mathrm{\hat{O}}_{ij}=-\frac{\partial
G_{i}}{\partial\rho_{j}}+K_{ij}\frac{\mathrm{d}^{2}}{\mathrm{d}z^{2}}$ (114)

is a second-order differential operator. Since the matrices $\partial
G_{i}/\partial\rho_{j}$ and $K_{ij}$ are symmetric (see definition (3) of
$G_{i}$ and subsection 2.5, respectively), $\mathrm{\hat{O}}_{ij}$ is self-
adjoint – hence, its spectrum is real. Note also that $\mathrm{\hat{O}}_{ij}$
depends on the steady state $\rho_{i}(z)$ via the coefficient $\partial
G_{i}/\partial\rho_{j}$.

It follows from equation (113) that the functional $H$ has a maximum at
$\rho_{i}(z)$ if and only if the corresponding operator
$\mathrm{\hat{O}}_{ij}$ is negative definite – or, equivalently, its spectrum
(the set of all discrete and continuous eigenvalues) is negative.

Theorem 1. If the solution $\rho_{i}(z)$ of boundary-value problem (68)–(70)
is monotonic for all $i$, the spectrum of the corresponding operator
$\mathrm{\hat{O}}_{ij}$ is negative.

Proof. Let $\Lambda$ be an eigenvalue of the _discrete_ spectrum (if any) and
$\psi_{i}$, the corresponding eigenfunction,

$\sum_{j}\left(-\frac{\partial
G_{i}}{\partial\rho_{j}}\psi_{j}+K_{ij}\frac{\mathrm{d}^{2}\psi_{j}}{\mathrm{d}z^{2}}\right)=\Lambda\psi_{i},$
(115) $\psi_{i}\rightarrow 0\qquad\text{as}\qquad z\rightarrow\pm\infty.$
(116)

Introduce

$\phi_{i}=\left(\frac{\mathrm{d}\rho_{i}}{\mathrm{d}z}\right)^{-1}\psi_{i},$
(117)

and observe that, since $\rho_{i}(z)$ is monotonic,
$\mathrm{d}\rho_{i}/\mathrm{d}z$ never vanishes and $\phi_{i}(z)$ is non-
singular.

Rewriting boundary-value problem (115)–(116) in terms of $\phi_{i}$ and
keeping in mind that $\rho_{i}(z)$ satisfies equation (70), one obtains

$\sum_{j}K_{ij}\left(2\frac{\mathrm{d}^{2}\rho_{j}}{\mathrm{d}z^{2}}\frac{\mathrm{d}\phi_{j}}{\mathrm{d}z}+\frac{\mathrm{d}\rho_{j}}{\mathrm{d}z}\frac{\mathrm{d}^{2}\phi_{j}}{\mathrm{d}z^{2}}\right)=\Lambda\frac{\mathrm{d}\rho_{i}}{\mathrm{d}z}\phi_{i},$
(118) $\frac{\mathrm{d}\rho_{i}}{\mathrm{d}z}\phi_{i}\rightarrow
0\qquad\text{as}\qquad z\rightarrow\pm\infty.$ (119)

Let $\left(K_{ij}\right)^{-1}$ be the inverse matrix to $K_{ij}$ (the latter
is positive definite – hence, invertible) and rewrite equation (118) in the
form

$\frac{\mathrm{d}}{\mathrm{d}z}\left[\left(\frac{\mathrm{d}\rho_{j}}{\mathrm{d}z}\right)^{2}\frac{\mathrm{d}\phi_{j}}{\mathrm{d}z}\right]=\Lambda\sum_{i}\left(K_{ij}\right)^{-1}\frac{\mathrm{d}\rho_{j}}{\mathrm{d}z}\frac{\mathrm{d}\rho_{i}}{\mathrm{d}z}\phi_{i}.$

Multiplying this equation by $\phi_{j}$, summing it with respect to $j$,
integrating from $z=0$ to $z=\infty$, integrating the left-hand side by parts,
and recalling boundary conditions (119), one obtains

$-\int_{-\infty}^{\infty}\sum_{j}\left(\frac{\mathrm{d}\rho_{j}}{\mathrm{d}z}\right)^{2}\left(\frac{\mathrm{d}\phi_{j}}{\mathrm{d}z}\right)^{2}\mathrm{d}z=\Lambda\int_{-\infty}^{\infty}\sum_{ij}\left(K_{ij}\right)^{-1}\left(\phi_{i}\frac{\mathrm{d}\rho_{i}}{\mathrm{d}z}\right)\left(\phi_{j}\frac{\mathrm{d}\rho_{j}}{\mathrm{d}z}\right)\mathrm{d}z.$
(120)

It can be shown that the integrands of both integrals in equation (120) decay
exponentially as $z\rightarrow+\infty$ – hence, the integrals converge.
Keeping also in mind that $K_{ij}$ is positive definite (hence,
$\left(K_{ij}\right)^{-1}$ is such), one concludes that equality (120) implies
that $\Lambda<0$, as required.

Next, let $\Lambda$ be an eigenvalue of the _continuous_ spectrum and
$\psi(z)$, the corresponding eigenfunction. They satisfy equation (115) and
the boundary conditions

$\psi_{i}\sim A_{i\pm}\operatorname{e}^{ik_{\pm}z}\qquad\text{as}\qquad
z\rightarrow\pm\infty,$ (121)

where $A_{i\pm}$ and $k_{\pm}$ are undetermined constants (the latter is
real). Substituting asymptotics (121) into equation (115), one obtains

$\sum_{j}\left[-\left(\frac{\partial
G_{i}}{\partial\rho_{j}}\right)_{\rho_{1}=\rho_{1}^{(l)}...\rho_{N}=\rho_{N}^{(l)}}-k_{-}^{2}K_{ij}\right]A_{j-}=\Lambda
A_{-},$ $\sum_{j}\left[-\left(\frac{\partial
G_{i}}{\partial\rho_{j}}\right)_{\rho_{1}=\rho_{1}^{(v)}...\rho_{N}=\rho_{N}^{(v)}}-k_{+}^{2}K_{ij}\right]A_{j+}=\Lambda
A_{+},$

Evidently, $\Lambda$ is a common eigenvalue of the matrices

$-G_{ij-}^{\prime}-k_{-}^{2}K_{ij}\qquad\text{and}\qquad-
G_{ij+}^{\prime}-k_{+}^{2}K_{ij}$

where

$G_{ij-}^{\prime}=\left(\frac{\partial
G_{i}}{\partial\rho_{j}}\right)_{\rho_{1}=\rho_{1}^{(l)}...\rho_{N}=\rho_{N}^{(l)}},\qquad
G_{ij+}^{\prime}=\left(\frac{\partial
G_{i}}{\partial\rho_{j}}\right)_{\rho_{1}=\rho_{1}^{(v)}...\rho_{N}=\rho_{N}^{(v)}}.$

Recall that, at infinity, the liquid and vapour are supposed to be stable,
which implies that $G_{ij\pm}^{\prime}$ are both positive definite. The
Korteweg matrix $K_{ij}$ is also positive definite – hence, all eigenvalues of
$-G_{ij\pm}^{\prime}-k_{\pm}^{2}K_{ij}$ are negative, as required.

###  C.2 1D drops and bubbles

The entropy maximization problem in this case can again be reduced to the
analysis of the operator $\mathrm{\hat{O}}_{ij}$. Furthermore, the second part
of the proof of Theorem 1 can be applied to 1D drops and bubbles without
modifications, and so the _continuous_ spectrum of $\mathrm{\hat{O}}_{ij}$ is,
again, negative.

The _discrete_ -spectrum part of Theorem 1, however, cannot be generalised for
non-monotonic $\rho_{i}(z)$. Indeed, if $\mathrm{d}\rho_{i}/\mathrm{d}z$
vanishes somewhere, the function $\phi_{i}(z)$ defined by (117) is singular,
and the integral on the left-hand side of (120) diverges. Thus, for 1D drops
and bubbles, $\mathrm{\hat{O}}_{ij}$ may have positive discrete eigenvalues,
but proving that it definitely does is not easy. In what follows, such a proof
is presented only for $N=1$.

Theorem 2. The operator $\mathrm{\hat{O}}_{11}$ with $\rho_{1}(z)$ satisfying
boundary-value problem (74)–(75) with $N=1$, has at least one positive
discrete eigenvalue (corresponding to an even eigenfunction).

Proof. For $N=1$, boundary-value problem (74)–(75) becomes

$K_{11}\frac{\mathrm{d}^{2}\rho_{1}}{\mathrm{d}z^{2}}=G_{1}-G_{\infty,1},$
(122) $\rho_{1}\rightarrow\rho_{\infty,1}\qquad\text{as}\qquad
z\rightarrow\pm\infty.$ (123)

It can be readily shown that

$\rho_{1}(z)\sim\rho_{\infty,1}+A\operatorname{e}^{\sqrt{G^{\prime}}z}\qquad\text{as}\qquad
z\rightarrow-\infty,$ (124)

where $G^{\prime}=\left(\partial
G_{11}/\partial\rho_{1}\right)_{\rho_{1}=\rho_{\infty,1}}$ and $A$ is a real
constant (which can be expressed via a certain integral, but its exact value
is unimportant). Note that $A$ is positive for drops [i.e., solutions of the
kind illustrated in figure 4(a)] and negative for bubbles [solutions
illustrated in figure 4(b)].

Next, eigenvalue problem (115)–(116) with $N=1$ has the form

$-\frac{\partial
G_{1}}{\partial\rho_{1}}\psi_{1}+K_{11}\frac{\mathrm{d}^{2}\psi_{1}}{\mathrm{d}z^{2}}=\Lambda\psi_{1},$
(125) $\psi_{1}\rightarrow 0\qquad\text{as}\qquad z\rightarrow-\infty,$ (126)
$\psi_{1}\rightarrow 0\qquad\text{as}\qquad z\rightarrow+\infty.$ (127)

Note that, for 1D drops an bubbles, $\rho_{1}(z)$ is even – hence, so is
$\partial G_{1}/\partial\rho_{1}$ which appears in equation (125). Since all
other coefficients in this (linear second-order) equation are constants, one
concludes that $\psi_{1}$ is either even or odd. Assuming the latter, one can
reduce the interval where equation (125) is to be solved to
$\left(0,\infty\right)$ and replace boundary condition (127) with

$\frac{\mathrm{d}\psi_{1}}{\mathrm{d}z}=0\qquad\text{at}\qquad z=0.$ (128)

Next, define an auxiliary function $\chi(z)$ which satisfies equation (125)
and boundary condition (126) – but not necessarily condition (128). To still
have two boundary conditions, require that

$\chi(z)\sim
A\sqrt{G^{\prime}+\Lambda}\operatorname{e}^{\sqrt{G^{\prime}+\Lambda}z}\qquad\text{as}\qquad
z\rightarrow-\infty,$ (129)

where $A$ is the constant from asymptotics (122). Note that $\chi(z)$ exists
for _any_ value $\Lambda$ (depends on it as a parameter). In particular,

$\chi(z)=\frac{\mathrm{d}\rho_{1}}{\mathrm{d}z}\qquad\text{for}\qquad\Lambda=0,$
(130)

which can be verified by comparing equation (125) with the derivative of
equation (122) and asymptotics (129) with the derivative of asymptotics (124).
For a large $\Lambda$, on the other hand, it can be deduced from equation
(125) and boundary condition (129) that

$\chi(z)\sim
A\sqrt{\Lambda}\operatorname{e}^{\sqrt{\Lambda}z}\qquad\text{as}\qquad\Lambda\rightarrow+\infty.$
(131)

This asymptotics applies to all finite $z$.

Now, consider a drop solution [e.g., one of those illustrated in figure 4(a).
It is evident from the figure (and common sense) that
$\left(\mathrm{d}^{2}\rho_{1}/\mathrm{d}z^{2}\right)_{z=0}<0$ and $A>0$ –
hence, (130)–(131) imply

$\left(\frac{\mathrm{d}\chi}{\mathrm{d}z}\right)_{z=0}<0\qquad\text{for}\qquad\Lambda=0,$
$\left(\frac{\mathrm{d}\chi}{\mathrm{d}z}\right)_{z=0}>0\qquad\text{as}\qquad\Lambda\rightarrow+\infty.$

These two inequalities indicate that there exists a
$\Lambda\in\left(0,+\infty\right)$ such that
$\left(\mathrm{d}\chi/\mathrm{d}z\right)_{z=0}=0$ – hence, the corresponding
function $\chi(z)$ satisfies boundary condition (128) [as well as condition
(126) and equation (125)]. The corresponding (positive) $\Lambda$ is,
obviously, an eigenvalue, as required.

The case of bubble solutions is similar to that of drops (examined above) and
yields the same conclusion.

###  C.3 Solid/fluid interfaces

In this case, $\rho_{i}(z)$ satisfies equation (75) and boundary conditions
(79)–(80).

Theorem 3. The spectrum of the operator $\mathrm{\hat{O}}_{ij}$ corresponding
to a _monotonic_ $\rho_{i}(z)$ is real and negative.

The proof of this theorem is similar to that of Theorem 1. The only difference
is in the boundary conditions: since the density is constrained to a fixed
value at the substrate, its variation is zero – and the eigenfunctions of the
operator $\mathrm{\hat{O}}_{ij}$ should satisfy

$\psi_{i}=0\qquad\text{at}\qquad z=0.$

This boundary condition should hold for both discrete and continuous spectra.

Theorem 4. If $N=1$, the operator $\mathrm{\hat{O}}_{11}$ with non-monotonic
$\rho_{1}(z)$ has at least one positive discrete eigenvalue.

The proof of this theorem is similar to that of Theorem 2

##  Appendix D 1D pure-fluid reduction of equations (56)–(57) and (51)

To adapt asymptotic equation (56)–(57) for a pure fluid, recall that the
transport coefficients $D_{11}$ and $\zeta_{1}$ are zero in this case (because
this is the only possibility to satisfy restrictions (25) for $N=1$). Thus,
one obtains

$\frac{\partial\rho_{1}}{\partial
t}+\mathbf{\nabla}\cdot\left(\rho_{1}\mathbf{v}\right)=0,$ (132)
$\mathbf{\nabla}\cdot\boldsymbol{\Pi}+\rho_{1}\mathbf{\nabla}\left(K_{11}\nabla^{2}\rho_{1}-G_{1}\right)=0.$
(133)

Let the flow be 1D, so that the unknowns depend only on $z$ and $t$, and the
velocity has only one component $\mathbf{v}=\left[0,0,w\right]$. Then,
equations (132)–(133) and expression (51) for the viscous stress yield

$\frac{\partial\rho_{1}}{\partial
t}+\frac{\partial\left(\rho_{1}w\right)}{\partial z}=0,$ (134)
$\frac{\partial}{\partial
z}\left[\left(\frac{4}{3}\mu_{s}+\mu_{b}\right)\frac{\partial w}{\partial
z}\right]+\rho_{1}\frac{\partial}{\partial
z}\left(K_{11}\frac{\partial^{2}\rho_{1}}{\partial z^{2}}-G_{1}\right)=0.$
(135)

Let the substrate be located at $z=0$, so that boundary conditions (32) and
(35) yield

$w=0,\qquad\rho_{1}\rightarrow\rho_{0,1}\qquad\text{at}\qquad z=0.$ (136)

At infinity, the density tends to a fixed value and there should be no stress
– hence,

$\rho_{1}\rightarrow\rho_{\infty,1},\qquad\frac{\partial w}{\partial
z}\rightarrow 0\qquad\text{as}\qquad z\rightarrow+\infty.$ (137)

Using identity (6), one can replace in equation (135) $\rho_{1}\partial
G_{1}/\partial z\rightarrow\partial p/\partial z$ and then integrate this
equation with respect to $z$. Fixing the constant of integration via boundary
conditions (137)–(138), one obtains

$\left(\frac{4}{3}\mu_{s}+\mu_{b}\right)\frac{\partial w}{\partial
z}+K_{11}\left[\rho_{1}\frac{\partial^{2}\rho_{1}}{\partial
z^{2}}-\frac{1}{2}\left(\frac{\partial\rho_{1}}{\partial
z}\right)^{2}\right]=0.$ (138)

Given a suitable initial condition, equations (134) and (138), with boundary
conditions (136) fully determine $\rho_{1}(z,t)$ and $w(z,t)$.

For numerical simulations, equations (134) and (138) were discretised in $z$
using central differences, and the resulting set of ordinary differential
equations were solved using the MATLAB function ODE23tb (which can handle
stiff problems). This approach is usually referred to as the “method of lines”
(Schiesser, 1978).

## References

  * Bearman & Jones (1960) Bearman, R. J. & Jones, P. F. 1960 Statistical mechanical theory of the viscosity coefficients of binary liquid solutions. J. Chem. Phys. 33, 1432–1438.
  * Bearman & Kirkwood (1958) Bearman, R. J. & Kirkwood, J. G. 1958 Statistical mechanics of transport processes. XI. equations of transport in multicomponent systems. J. Chem. Phys. 28, 136–145.
  * Benilov (2020a) Benilov, E. S. 2020a Asymptotic reductions of the diffuse-interface model, with applications to contact lines in fluids. Phys. Rev. Fluids 5, 084003.
  * Benilov (2020b) Benilov, E. S. 2020b The dependence of the surface tension and contact angle on the temperature, as described by the diffuse-interface model. Phys. Rev. E 101, 042803.
  * Benilov (2020c) Benilov, E. S. 2020c Dynamics of liquid films, as described by the diffuse-interface model. Phys. Fluids 32, 112103\.
  * Benilov (2020d) Benilov, E. S. 2020d Nonexistence of two-dimensional sessile drops in the diffuse-interface model. Phys. Rev. E 102, 022802\.
  * Benilov (2021) Benilov, E. S. 2021 Can a liquid drop on a substrate be in equilibrium with saturated vapor? Phys. Rev. E 104, L032103.
  * Benilov (2022a) Benilov, E. S. 2022a Capillary condensation of saturated vapor in a corner formed by two intersecting walls. Phys. Fluids 34, 062103.
  * Benilov (2022b) Benilov, E. S. 2022b Dynamics of a drop floating in vapor of the same fluid. Phys. Fluids 34, 042104.
  * Benilov & Benilov (2015) Benilov, E. S. & Benilov, M. S. 2015 A thin drop sliding down an inclined plate. J. Fluid Mech. 773, 75–102.
  * Benilov & Benilov (2018) Benilov, E. S. & Benilov, M. S. 2018 Energy conservation and ${H}$ theorem for the Enskog-Vlasov equation. Phys. Rev. E 97, 062115\.
  * Benilov & Benilov (2019) Benilov, E. S. & Benilov, M. S. 2019 Peculiar property of noble gases and its explanation through the Enskog-Vlasov model. Phys. Rev. E 99, 012144.
  * Bertozzi & Flenner (2012) Bertozzi, A. L. & Flenner, A. 2012 Diffuse interface models on graphs for classification of high dimensional data. Multiscale Model. Simul. 10, 1090–1118.
  * Bertozzi & Flenner (2016) Bertozzi, A. L. & Flenner, A. 2016 Diffuse interface models on graphs for classification of high dimensional data. SIAM Review 58 (2), 293–328.
  * Bestehorn et al. (2021) Bestehorn, M., Sharma, D., Borcia, R. & Amiroudine, S. 2021 Faraday instability of binary miscible/immiscible fluids with phase field approach. Phys. Rev. Fluids 6, 064002.
  * Borcia & Bestehorn (2014) Borcia, R. & Bestehorn, M. 2014 Phase field modeling of nonequilibrium patterns on the surface of a liquid film under lateral oscillations at the substrate. Intl. J. Bifurc. Chaos 24, 1450110.
  * Borcia et al. (2019) Borcia, R., Borcia, I. D., Bestehorn, M., Varlamova, O., Hoefner, K. & Reif, J. 2019 Drop behavior influenced by the correlation length on noisy surfaces. Langmuir 35, 928–934.
  * Cahn (1961) Cahn, J. W. 1961 On spinodal decomposition. Acta Metallurg. 9, 795–801.
  * Cahn & Hilliard (1958) Cahn, J. W. & Hilliard, J. E. 1958 Free energy of a nonuniform system. i. Interfacial free energy. J. Chem. Phys. 28, 258–267.
  * Caupin (2005) Caupin, F. 2005 Liquid-vapor interface, cavitation, and the phase diagram of water. Phys. Rev. E 71, 051605.
  * Colinet & Rednikov (2011) Colinet, P. & Rednikov, A. 2011 On integrable singularities and apparent contact angles within a classical paradigm. Eur. Phys. J. Spec. Top. 197, 89–113.
  * Collins & Levine (1985) Collins, J. B. & Levine, H. 1985 Diffuse interface model of diffusion-limited crystal growth. Phys. Rev. B 31, 6119–6122.
  * Cramer (2012) Cramer, M. S. 2012 Numerical estimates for the bulk viscosity of ideal gases. Phys. Fluids 24, 066102.
  * Dai et al. (2017) Dai, M., Feireisl, E., Rocca, E., Schimperna, G. & Schonbek, M. E. 2017 Analysis of a diffuse interface model of multispecies tumor growth. Nonlinearity 30, 1639–1658.
  * Davidson (1993) Davidson, T. A. 1993 A simple and accurate method for calculating viscosity of gaseous mixtures. Tech. Rep.. U.S. Department of the Interior, Bureau of Mines.
  * Deegan et al. (2000) Deegan, R. D., Bakajin, O., Dupont, T. F., Huber, G., Nagel, S. R. & Witten, T. A. 2000 Contact line deposits in an evaporating drop. Phys. Rev. E 62, 756–765.
  * Ding & Spelt (2007) Ding, H. & Spelt, P. D. M. 2007 Wetting condition in diffuse interface simulations of contact line motion. Phys. Rev. E 75, 046708.
  * Ding et al. (2017) Ding, H., Zhu, X., Gao, P. & Lu, X.-Y. 2017 Ratchet mechanism of drops climbing a vibrated oblique plate. J. Fluid Mech. 835, R1.
  * Dunn et al. (2009) Dunn, G. J., Wilson, S. K., Duffy, B. R., David, S. & Sefiane, K. 2009 The strong influence of substrate conductivity on droplet evaporation. J. Fluid Mech. 623, 329–351.
  * Eggers & Pismen (2010) Eggers, J. & Pismen, L. M. 2010 Nonlocal description of evaporating drops. Phys. Fluids 22, 112101.
  * Engineering ToolBox (2004) Engineering ToolBox 2004 Solubility of air in water. [online] Available at: https://www.engineeringtoolbox.com/air-solubility-water-d_639.html.
  * Engineering ToolBox (2018) Engineering ToolBox 2018 Air - Diffusion coefficients of gases in excess of air. [online] Available at: https://www.engineeringtoolbox.com/air-diffusion-coefficient-gas-mixture-temperature-d_2010.html.