Towards a consistent lattice Boltzmann model for two-phase fluid

We propose a kinetic framework for single-component non-ideal isothermal flows. Starting from a kinetic model for a non-ideal fluid, we show that under conventional scaling the Navier-Stokes equations with a non-ideal equation of state are recovered in the hydrodynamic limit. A scaling based on the smallness of velocity increments is then introduced, which recovers the full Navier-Stokes-Korteweg equations. The proposed model is realized on a standard lattice and validated on a variety of benchmarks. Through a detailed study of thermodynamic properties including co-existence densities, surface tension, Tolman length and sound speed, we show thermodynamic consistency, well-posedness and convergence of the proposed model. Furthermore, hydrodynamic consistency is demonstrated by verification of Galilean invariance of the dissipation rate of shear and normal modes and the study of visco-capillary coupling effects. Finally, the model is validated on dynamic test cases in three dimensions with complex geometries and large density ratios such as drop impact on textured surfaces and mercury drops coalescence.


Introduction
Multi-phase flows are omnipresent in science and technology. From micro-droplets coalescing in clouds, to solidification or melting of alloys and diesel droplets evaporation and subsequent combustion, all involve multiple interacting phases and moving interfaces. This ubiquity fueled wide efforts focused on the development of predictive mathematical models and numerical tools for multi-phase flows. While significant attention has been focused on sharp interface methods requiring efficient tracking of the evolving and deforming interfaces, and imposing jump conditions (Sethian & Smereka 2003;Scardovelli & Zaleski 1999;Popinet 2018;Prosperetti & Tryggvason 2009), the evergrowing range of temperatures and pressures involved in typical systems of interest is making thermodynamic consistency of the computational models essential. For instance, dramatically different thermodynamic regimes are encountered in diesel engines during the compression phase, in aeoronautical engines during take-off while most rocket engines operate in trans-and super-critical regimes, where the interface thickness becomes comparable to the flow scales. Nucleation and cavitation are yet another example, where the sharp interface limit does not hold and modifications to the classical nucleation theory (Debenedetti 1997), related to curvature-dependence of the surface tension, are required. In such cases, an accurate account of non-ideality of the fluid, including a finite † Email address for correspondence: ikarlin@ethz.ch arXiv:2112.01975v1 [physics.flu-dyn] 3 Dec 2021 interface thickness, is crucial for predictive numerical simulations of the flow physics. At a macorscopic level, a primer example for thermodynamics of non-ideal fluids is the second-gradient theory, first introduced by van der Waals for single-component fluids (van der Waals 1894), leading to the Navier-Stokes equations supplemented with the Korteweg stress tensor (Korteweg 1901), and is a starting point for numerical methods known as diffuse interface approach (Anderson et al. 1998). On the other hand, extension of the Boltzmann equation to dense gases within the Enskog hard-sphere collision model (Enskog 1921) and Vlasov mean-field approximation (Vlasov 1961) provides a kinetictheory basis for dynamics of non-ideal fluid (Chapman & Cowling 1939).
Since the pioneering work of Shan & Chen (1993), the lattice Boltzmann method (LBM) gained popularity as a viable numerical tool targeting the hydrodynamic regime of multi-phase flows. Despite their popularity and wide usage, most multi-phase models for LBM, apart from limited studies in He & Doolen (2002), lack a clear kinetictheory framework. For instance, the so-called pseudo-potential lattice Boltzmann models lack a clear and consistent continuous kinetic model and scaling law recovering the full target macroscopic system. Thorough analyses of the bulk thermodynamic and interface properties of the models (especially near the critical state) are also very scarce. Furthermore, ever since their inception, such models have continuously struggled with larger density ratio simulations achieving at best, via different strategies, ratios of the order of 10 3 as reflected by a number of recent reviews Li et al. 2016).
In this paper, we revisit the construction of the lattice Boltzmann model for isothermal two-phase flows. We propose a flexible kinetic framework for dense fluids with non-ideal equations of states. Using the lattice Boltzmann method discretization strategy, and under proper scaling, the model is shown to recover the full Navier-Stokes-Korteweg system of equations. Through a detailed study of thermodynamic properties, the model is shown to be well-posed and convergent to the capillary fluid thermodynamics. The well-posedness of the model and proper consideration of the proposed scaling is shown to guarantee recovery of the hydrodynamic-scale dynamics both at very large density ratio and near critical point.
The outline is as follows: We begin in section 2.1 with a summary of the second-gradient theory due to van der Waals (1894) and Korteweg (1901). In section 2.2, following a more microscopic approach, we consider a class of kinetic models suitable for a non-ideal fluid. We proceed in section 2.3 with a scaling assumption of small flow velocity increments which leads to a lattice Bhatnagar-Gross-Krook (LBGK) equation with a new realization of the nonlocal force that guarantees consistency with Korteweg's stress.
Thermodynamics of the LBGK model is validated in section 3. In section 3.1, we demonstrate convergence of vapour-liquid coexistence to the principle of corresponding states (Guggenheim 1945), independently of the equation of state and for liquid-vapour density ratio up to at least ∼ 10 11 . The remainder of the paper is based on the van der Waals equation of state. In section 3.2, we show that the surface tension in the present LBGK model obeys a temperature scaling in excellent agreement with the theory. After verifying that the proposed model allows for choosing surface tension independently of the density ratio in section 3.3, we show in section 3.4 that it is also consistent with Gibbs' theory of dividing surfaces (Gibbs 1874). Simulations presented in section 3.4 reveal a generalized Laplace law and uncover the effect of curvature on the surface tension, in agreement with the theory by Tolman (1949). Finally, in section 3.5, we show that the interface width scales with the temperature in accord with van der Waals theory.
We turn to probing hydrodynamic features of our model in section 4. In section 4.1, we demonstrate that it correctly implements the jump condition for the stresses at the liquid-vapour interface in the simulation of layered Poiseuille flow. In section 4.2, Galilean invariance is demonstrated by measuring dissipation of normal modes in a moving reference frame. The viscosity-capillarity coupling is probed in section 4.3 by measuring the frequency of higher-order capillary waves, in excellent agreement with Rayleigh's theory (Rayleigh 1879). We also demonstrate that damping rate of capillary wave agrees with analytical solution. Validation of bulk properties is concluded by measuring the isothermal speed of sound in section 4.4, where excellent comparison to the theoretical prediction is demonstrated for large density ratios. In section 4.5, the model is extended to the simulation of a fluid-solid interface, and is is validated by demonstrating the Young-Laplace law and a liquid column motion in a channel with non-uniform wettability. In section 5, the model is used to simulate water impact on textured superhydrophobic surfaces and mercury droplets coalescence to demonstrate its ability to handle simulations at extremely high density ratios. Conclusions are drawn in section 6.

Second-gradient theory: Korteweg's stress and capillary fluid equations
In the second-gradient theory as introduced by van der Waals (1894), free energy per unit volume is expressed as: where A is the bulk free energy per unit volume, ρ is the density and κ is the capillary coefficient. The second term represents the interface energy while the bulk free energy is solely a function of the local density and temperature (Giovangigli 2020). The equilibrium state of the corresponding system is obtained by minimizing free energy in a given volume under the constraint of constant total mass, leading to the stress tensor (Anderson et al. 1998), where I is unit tensor and L is the Lagrange function, and λ is the Lagrange multiplier for the mass constraint, or chemical potential, where ∇ 2 is the Laplace operator. This in turn leads to the following Korteweg's stress tensor (Korteweg 1901): is the thermodynamic pressure, or equation of state. From the local balance equations for mass and momentum one obtains the macroscopic governing laws for an iso-thermal capillary fluid: where u is the fluid velocity and the stress tensor T is The Navier-Stokes viscous stress tensor reads, where S is the trace-free rate-of-strain tensor, (2.11) and µ and η are the dynamic and the bulk viscosity, respectively. The momentum balance equation (2.8) can be recast in the following form, where Korteweg's force F K is the divergence of the Korteweg pressure tensor, The latter can be written in the following form, where we have introduced a reference pressure P 0 . Navier-Stokes momentum equations with Korteweg's force (2.14) shall be a target for reconstruction by a suitable kinetic model.

Kinetic model for non-ideal fluid
In order to introduce a kinetic model for non-ideal fluid, we begin with the first Bogolioubov-Born-Green-Kirkwood-Yvon (BBGKY) equation, where f (r, v, t) and f 2 (r, v, r 1 , v 1 , t) are the one-and the two-particle distribution functions, respectively, r, r 1 and v, v 1 are particles position and velocity, while V is a potential of pair interaction. The local equilibrium state is defined by the Maxwellian f eq at constant temperature T , parameterized by the local values of density ρ and flow velocity u, where R is gas constant. Furthermore, let us introduce a projector K onto local equilibrium at constant temperature, (2.17) The projector property, K 2 = K, can be verified by a direct computation. With the projector (2.17), the interaction term in (2.15) is split into two parts by writing an identity, J = (1 − K) J + KJ . (2.18) The first term, (2.19) satisfies the local conservation of both mass and momentum, It is conventional to model the locally conserving part of the interaction with a single relaxation time Bhatnagar-Gross-Krook (BGK) approximation, where the relaxation time τ is a free parameter. The second term in the identity (2.18), satisfies the local mass but not the local momentum conservation. Indeed, after integration by part in the velocity v and neglecting boundary integrals, we arrive at where the force F nloc reads, Collecting the BGK approximation together with the nonlocal contribution, a generic kinetic model may be written, Evaluation of the force (2.24) requires us to specify the particles interaction. It is customary to invoke the Enskog-Vlasov model (Enskog 1921;Vlasov 1961) where both hard-sphere collisions and a weak long-range attraction potential contribute to a non-local momentum transfer. For the hard-sphere Enskog part, a de-localization of the collision is responsible for a non-vanishing contribution of momentum transfer through the distance between the centres of the spheres upon their impact while the Vlasov approximation contributes non-locally to the momentum transfer from a distributed mean-field force. Evaluation of both the Enskog and Vlasov contributions to the force (2.24) proceeds along familiar lines (Chapman & Cowling 1939;He & Doolen 2002) and is reported in Appendices A and B for completeness, The first term is the lowest-order contribution of the collisional momentum transfer from the Enskog hard-sphere model, (2.27) Here b = 4v HS , with v HS = V HS /m the specific volume of hard-sphere of diameter d and mass m, while V HS = πd 3 /6 is the volume of the sphere. Moreover, χ is the equilibrium two-particle correlation function, evaluated at the local density reduced by the specific volume of hard sphere, χ = χ(bρ(r, t)); To the lowest order, χ = 1 + (5/8)bρ + O((bρ) 2 ), cf. Chapman & Cowling (1939). The second term in (2.26) is the contribution of a longrange attraction potential V in the mean field Vlasov approximation. To third order in the gradient of density, where parameters a and κ are, Thus, with the approximations specified, the non-local force (2.26) becomes, where the reference ideal gas pressure P 0 is provided by the local Maxwellian (2.16), while the equation of state of the Enskog-Vlasov gas is of van der Waals type, This allows us to extend the Enskog-Vlasov kinetic model and a phenomenological equation of state P (2.6) can be used instead of P EV (2.33). Moreover, the reference pressure P 0 can be made selective by rescaling the local equilibrium, where D is the space dimension. While the Enskog-Vlasov partition above corresponds to selecting P 0 = ρRT , an alternative is provided by Reyhanian et al. (2020), where P 0 = P is chosen. Using the rescaled equilibria (2.34), a family of kinetic models parameterized by the reference pressure reads, The kinetic equation (2.35) shall be considered as a semi-phenomenological model of nonideal fluid, with the relaxation time τ , the capillarity coefficient κ, pressure P and reference pressure P 0 as phenomenological input parameters, while the Enskog-Vlasov realization will serve as a representative example for estimates of various flow regimes. The analysis of the kinetic model (2.35) under the conventional scaling of a small deviation from a uniform state (Chapman & Cowling 1939), (2.36) is detailed in Appendix C. To second order in space derivatives, the resulting momentum balance equation reads, where the dynamic viscosity µ and the bulk viscosity η in the Navier-Stokes stress tensor (2.10) are defined by the reference pressure (D = 3), Thus, the momentum balance equation (2.37) is form-invariant with respect to the choice of reference pressure, provided P 0 satisfies a sub-isentropic condition, for some C > 0. With (2.40), the bulk viscosity (2.39) is positive and vanishes when the reference pressure follows an isentropic process for ideal monatomic gas, P 0 = Cρ 5/3 . For example, any polytropic process, P 0 = Aρ n , 1 n 5/3 satisfies the sub-isentropic condition and results in η = (5/3 − n)τ P 0 . Special case of isothermal process n = 1 returns η = (2/3)τ P 0 , and the viscous stress tensor becomes, On the other hand, when compared to the two-phase momentum equation (2.12), the macroscopic limit recovers only the nonideal gas component thereof while missing Korteweg's capillarity contribution. Indeed, the third-order term, ∼ 3 ρ∇∇ 2 ρ in (2.28) and (2.35), does not contribute to the momentum equation (2.37) under the scaling (2.36). This is consistent with the well-known results from kinetic theory (Chapman & Cowling 1939) and is not surprising: The scaling (2.36) is essentially based on the Knudsen number, which overrides the relative contribution of the capillarity term by two orders, cf. Appendix C. Thus, under weak non-uniformity assumption (2.36), the capillarity terms are seen as higher-order, Burnett-level contributions, and cannot appear in the main (first and second) orders in the momentum balance equation (2.37). In fact, condition (2.36) rules out situations at an interface between phases where gradients of density become large over a relatively short distance. Therefore, in order for the kinetic model (2.35) to recover in-full the momentum balance (2.12), a different scaling needs to be applied.

Time step and forcing
A rescaling of the kinetic model in this section shall be maintained by introducing a time step δt. As a preliminary consideration, we evaluate the contribution of the force term over the time step. To that end, as noted by Kupershtokh et al. (2009), for a generic force F , we can write the action of the force on the distribution function as a full derivative in a frame that moves with the local fluid velocity, 1 ρ ∂f eq ∂u · F = df eq dt . (2.42) Introducing the velocity increment, and integrating in time, leads to an approximation, This so-called exact difference method (EDM) becomes accurate for a gravity force, F /ρ = const, otherwise it often provides a reliable estimate for the force term and is widely used. In what follows, the scaling to be applied assumes smallness of the velocity increment (2.43) rather than smoothness of the spatial distribution of the force. Since the velocity increment is based on a time step, it is natural to proceed with a lattice Boltzmann realization of the kinetic equation.

Standard lattice and product-form
The lattice Boltzmann model shall be realized with the standard discrete velocity set D3Q27, where D = 3 stands for three dimensions and Q = 27 is the number of discrete velocities, We first define a triplet of functions in two variables, ξ α and ζ αα , and consider a product-form associated with the discrete velocities c i (2.45), All pertinent populations below shall be determined by specifying the parameters ξ α and ζ αα in the product-form (2.49). A two-dimensional version of the model on the D2Q9 lattice is obtained by omitting the z-component in all formulas below. Finally, we use notation δr i = c i δt for the lattice links, and denote the grid spacing in any direction α = x, y, z as δr = |c iα |δt, c iα = 0: For the D3Q27 discrete velocity set (2.45), the lattice spacing is same in all Cartesian directions.

The lattice Boltzmann equation
The local density ρ and flow velocity u are defined using the populations f i , With a reference pressure P 0 , and by setting the parameters, the local equilibrium populations are represented with the product-form (2.49), (2.54) The LBGK equation, supplemented with a forcing term, is written, where ω is a dimensionless relaxation parameter, the equilibrium populations are provided by (2.54) while the extended equilibrium populations f * i are defined by the productform (2.49) with the following assignment for the parameters, where Φ αα /ρ is a correction term for the diagonals of the non-equilibrium momentum flux tensor, where ς = δr/ √ 3δt is the so-called lattice speed of sound. Thus, the extended equilibrium is explicitly written as, (2.59) Comments are in order: • If the correction term (2.58) is omitted in (2.57), the population (2.59) becomes the equilibrium with the increment δu (2.43) due to force action added to the flow velocity u. This corresponds to the EDM forcing maintained by the second bracket in the right hand side of the LBGK equation (2.55).
• The correction term (2.56) is added following a proposal by Saadat et al. (2021). Its purpose is to compensate for the bias of the D3Q27 velocities (2.45), c 3 iα = c iα , and to restore the Galilean invariance of the normal components of the viscous stress tensor.
The LBGK model (2.55) is generic with respect to the choice of reference pressure P 0 and the force F . We now proceed with the special case of Korteweg's force in order to establish a representation thereof matched to the lattice Boltzmann system.

Pseudo-potential and capillarity
Following the representation (2.14), Korteweg's force includes two distinct parts, the term supplying the nonideal gas equation of state and the capillarity term responsible for the surface tension. Introducing a pseudo-potential ψ, Korteweg's force is written, In the lattice Boltzmann setting, the pseudo-potential part is represented as a linear combination of the first-and second-neighbours contributions, δtψ(r)∇ψ(r) = ψ(r) where the weights w i are defined by the product-form (2.49) at ξ α = 0, ζ αα = ς 2 , (2.63) and where the parameters G 1 and G 2 satisfy the conditions, (2.61) is represented in a similar way but using the density instead of the pseudo-potential, (2.66) whereκ = κδr 2 and the parameters G 3 and G 4 satisfy the conditions, Condition (2.67) cancels the first-order derivative, while condition (2.68) maintains the capillarity contribution. Combining both the pseudo-potential (2.62) and the capillarity (2.66) contributions, we obtain the lattice Boltzmann form of Korteweg's force (2.61) as, where the sign convention follows (2.61).
In the context of the lattice Boltzmann method, various pseudo-potential representations have long been in use, and comments are in order to make a distinction to the present formulation.
• By setting G 2 = 0 in (2.62) and ignoring Korteweg's capillarity term (2.66) by choosing G 3 = G 4 = 0, one recovers a model first proposed by Shan & Chen (1993) for special equations of state (see Eqs. (G 1), (G 2), (G 3) in Appendix G) and extended to a general equation of state by Yuan & Schaefer (2006). Unlike the above condition (2.65) which eliminates the third-order error in the pseudo-potential part of Korteweg's force (2.14), the model of Shan & Chen (1993) requires the third-order error to be retained in order that it mimics surface tension effects. Consequently, the force in the models of Shan & Chen (1993); Yuan & Schaefer (2006) becomes, The force (2.70) neither conforms with the van der Waals second-gradient theory and Korteweg's stress of section 2.1 (unless ψ = Aρ which, however, does not lead to phase separation), nor can it be derived from the BBGKY equation with a central-force particles interaction of section 2.2. Another, relatively minor issue is the fixed parameter 1/3, which is mimicking the capillarity coefficient and is the result of the third-order error retained in the expansion.
• By imposing condition (2.64) while still discarding Korteweg's capillarity contribution (2.66), one arrives at a dual-range force model of Sbragaglia et al. (2007), The model (2.71) is an improvement on the first-neighbour model (2.70) in that it allows for a variable capillarity-like coefficient. At the same time, it does not resolve the inconsistency with Korteweg's stress tensor.
• Various other modifications of the original models of Shan & Chen (1993) were proposed to improve on the main inconsistency by introducing and tuning ad hoc numerical coefficients tailored to a selected equation of state (Luo et al. 2021;Li et al. 2013;Huang et al. 2019;Kupershtokh et al. 2009). However, to the best of our knowledge, none of these proposals came to recognize that the problem lies in the fact that it is impossible to represent both parts of Korteweg's force, the non-ideal equation of state and the capillarity term while using a pseudo-potential alone. This fact follows both from the phenomenological thermodynamics reviewed in section 2.1, as well as from a more microscopic approach of section 2.2. The pseudo-potential in both cases represents only the equation of state while the capillarity term requires the density field to be used, as featured by our Eq. (2.66).
• Pseudo-potential is a convenient form of representing the equation of state contribution to Korteweg's force, tailored to the lattice Boltzmann setting. If other numerical methods are used to evaluate the force (2.61), such as higher-order finite difference, the model is usually renamed to a free energy approach.
With the generic LBGK model (2.55) and Korteweg's force (2.69) both specified, we proceed to the analysis of the hydrodynamic limit under a suitable scaling.

Hydrodynamic limit under small velocity increment scaling
Chapman-Enskog analysis of the LBGK equation (2.55) was performed under the following scaling: With the characteristic values of the flow velocity U, the flow scale L, the density ρ, the force F and the velocity increment δu, the following conditions apply: (2.73) The first scaling condition (2.72) refers to a smallness of velocity increment, that is, to the smallness of the force action over time δt. The second scaling condition (2.73) is a resolution requirement. Both conditions are assumed to hold simultaneously. Details of the analysis are provided in Appendix D while the result is summarized below. Let us introduce a transformed velocity U by shifting the local velocity u by half of the velocity increment δu (2.43), (2.74) Then the following mass and momentum balance equations to second order in ε are recovered when the force (2.69) is used under the scaling (2.72) and (2.73), where the dynamic and the bulk viscosity in the Navier-Stokes stress tensor (2.10) are related to the relaxation parameter ω and the reference pressure P 0 as follows: Unlike the previous result (2.37), the momentum balance (2.76) includes not only the nonideal gas pressure but also the capillarity term, and is thus consistent with Korteweg's force in the momentum balance. It should be pointed out that the scaling (2.72) refers to smallness of the increment of the flow velocity rather that to smallness of either the time step or of the force. Thus, rescaling the kinetic model (2.35) based on the smallness of flow velocity increments results in both the non-ideal gas equation of state and the capillarity revealed at the Euler level O(ε) of the momentum balance (2.76). This is in a contrast to the conventional scaling, which is tight to the nonuniformity and surface tension would appear only at a Burnett level O( 3 ).
In the remainder of this paper, and without loss of generality, we set the reference pressure P 0 = ς 2 ρ in order to minimize the correction term (2.58). In the next section, the model is scrutinized by a set of numerical tests probing various aspects of thermoand hydrodynamic consistency.

Liquid-vapour coexistence: The principle of corresponding states
We begin with the validation of liquid-vapour coexistence. Two-dimensional flat interface simulations were performed on a grid 800 × 10, filled with the vapour phase of a fluid with a specified equation of state and periodic boundary conditions. A column of the liquid phase over 400 grid-points was placed at the centre of the domain. Simulation were ran until steady-state was reached. Steady-state was monitored via a L 0 norm convergence criterion based on the liquid density ρ l at the centre of the drop and the vapour density ρ v at a location outside the drop. A theoretical prediction for the coexistence density ratio ρ l /ρ v is readily obtained via the equilibrium condition leading to the Maxwell equal-area construction, where P sat (T ) is the saturation pressure at which the liquid and vapor phases coexist at a given temperature T below the critical point. Initially four generally adopted equations of state (EoS) were considered: the van der Waals EoS (van der Waals 1873), where parameters a and b are related to critical temperature T c and pressure P c as, a = 27 64 the Peng-Robinson EoS (Peng & Robinson 1976), where ω the acentric factor (ω = 0.344 for water), and the Riedlich-Kwong-Soave EoS (Redlich & Kwong 1949;Soave 1972), and the Carnahan-Starling EoS (Carnahan & Starling 1969), (3.11) Fig. 1 demonstrates that stationary density ratios ρ l /ρ v obtained from the simulation are in excellent agreement with the theoretical coexisting liquid-vapour density ratios that are defined by Maxwell's equal-area rule (3.1), for all four equations of state and for ratios as high as at least ρ l /ρ v ∼ 10 11 . It is noted that high coexisting density ratios were obtained without any tuning parameters, universally for all equations of state considered. Therefore, and without loss of generality, in the remainder of this article we only consider the van der Waals equation of state (3.2).
A discussion on the principle of corresponding states and the necessity of adherence to it in the simulation of realistic systems at large density ratios is in order. According to Guggenheim (1945), the principle of corresponding states is "the most useful byproduct of van der Waals' equation of state". The principle maintains that all properties that depend on inter-molecular forces are related to the critical properties of the substance in a universal way, regardless of the molecular compound of interest. For the equation of state, the principle of corresponding states implies that the reduced pressure P r = P/P c is a universal function of the reduced temperature T r = T /T c and of the reduced density The universality of the reduced pressure (3.12) can be used to write Maxwell's equal-area rule in reduced form, The coexistence density ratio ρ l /ρ v at a given reduced temperature T r is therefore also universal. In order to probe the consistency with the principle of corresponding states in our model, we can be explained as follows: The local value of the scaling parameter ε (2.72) is estimated as ε l,v ∼ δtF/ρ l,v U on the liquid and on the vapour sides of the interface, respectively. Hence, their relative magnitude scales as the inverse of the density ratio, Thus, even if the scaling condition (2.72) is satisfied on the liquid side, ε l 1, it is still prone to violation on the vapour side, if the density ratio ρ l /ρ v becomes sufficiently large.
Furthermore, due to the scaling relation between interface width and temperature detailed in section 3.5, the interface gets thinner at lower temperatures. This in turn means that, for a given δr, fewer grid-points resolve the interface with decreasing temperature. Whenever the parameters δr/W and δtF /ρU increase, contributions of higher orders in ε become significant and lead to deviations from the analytical predictions. The characteristic interface width scales as W ∝ 1/ √ a (Jamet et al. 2001). As such lower values of a lead to larger W , in parallel with smaller force increments over δt, which restores dominance of order ε terms over Burnett and super-Burnett level contributions, and therefore the corresponding states principle. This last point is demonstrated by the convergence of the coexistence density ratio to the analytical predictions with decreasing a. For well-resolved simulations, as shown in Fig. 1, the model correctly recovers the coexistence densities and thus complies with the principle of corresponding states. Below, we refer to the convergence of the scheme to the principle of corresponding states as the thermodynamic convergence.

Temperature dependence of the surface tension near the critical point
Surface tension at liquid-vapour interface decreases with increasing temperature and vanishes at the critical point (Guggenheim 1945). For the van der Waals equation of state, the surface tension coefficient σ follows a scaling law as T r → 1 (van der Waals 1894; Blokhuis & Kuipers 2006), (3.14) In order to probe the consistency of the proposed lattice Boltzmann model, the temperature dependence of the surface tension was numerically measured in two different configurations, the flat interface and the circular drop, in a temperature interval T r ∈ [0.85, 1].
In the first configuration, the surface tension coefficient was evaluated using its definition for the flat interface (Kirkwood & Buff 1949), where the interface is normal to the x-axis in a two-dimensional simulation setup. The normal P xx and the tangential P yy components of the discrete pressure tensor were computed using a formalism developed in Appendix E, following a proposal by Shan (2008).
In the second configuration, simulations of circular liquid drops surrounded with vapour at the center of a square domain were conducted. At each temperature, four different initial drop radii were considered, R 0 ∈ {45, 55, 65, 75}, chosen in such a way that the interface width is sufficiently small compared to initial drop radius. In the simulation, we used W/R 0 0.1, where W = (ρ l − ρ v )/ max|∇ρ| is the interface width, in order to minimize curvature-dependence of surface tension. The corresponding surface tension coefficient was evaluated using the Laplace law (D = 2) in a form, where the radius R e corresponds to the equimolar dividing surface (Gibbs 1874). A brief reminder of Gibbs' theory of dividing surfaces is in order. The total mass in both the diffuse and sharp interface pictures can be written as: where ρ l V l and ρ v V v are the masses in the bulk liquid and vapor phases in the sharp interface picture, while Γ is the excess mass on a dividing surface Σ, or mass adsorbance (Gibbs 1874). By requiring that no mass be stored on the dividing surface we get the definition of the equimolar surface: The family of dividing surfaces in the case of drop or bubble are concentric spheres (D = 3) or concentric circles (D = 2) parameterized by their radius R. In particular, for a two-dimensional drop, the mass adsorbance can be written as a function of the radius of the dividing circle, while the zero-adsorbtion condition (3.18), Γ (R e ) = 0, implies the equimolar radius R e , (3.20) The drop configuration along with the scaling of the pressure difference across the interface with drop radius are shown in Fig. 3. The results as obtained from both the flat interface and the drop configurations are shown in Fig. 4. It is clearly observed that the surface tensions as obtained from the proposed formulation (using either one of the considered configurations) agree very well with Eq. (3.14), provided that W R e . Discussion of curvature-dependence of surface tension shall be continued in sec. 3.4.

Control of surface tension
The present formulation allows us to select the surface tension in the simulation via the capillarity parameter κ, independently from the density ratio and temperature. In order to demonstrate this feature, flat interface simulations were performed at three reduced temperatures, T r = {0.99, 0.98, 0.97}, for different values ofκ. The surface tension was evaluated using Eq. (3.15). Results in Fig. 5 demonstrate that the surface tension can be effectively tuned usingκ and that changes in surface tension do not affect the equilibrium density ratios. This is further detailed in Table 1 where the values of surface tension and deviations of the vapor and liquid densities from theoretical values are given as a function  ofκ for T r = 0.97. It is clearly observed that while the surface tension covers two orders of magnitude the deviation from the predicted density of vapour is at most 0.041 percent. Furthermore, in the limit of vanishing κ the surface tension also vanishes, as expected from the theory, Eq. (3.14).

Effect of curvature on surface tension: Tolman length
The drop simulation in section 3.2 made use of the Laplace law, relying on the equimolar dividing surface of radius R e (3.20). Further discussion on the non-uniqueness of the choice of dividing surface and curvature-dependence of surface tension is in order. Following Gibbs (1874), the free energy of a drop or bubble separated from the surrounding vapour or liquid by a dividing circle (D = 2) or sphere (D = 3) of length or area Σ is, A = U − T S + σΣ, where U and S are the internal energy and entropy of bulk phases while the last term is the adsorbance of free energy. The equilibrium condition requires vanishing of the variation δA; for the isothermal case we have, where P l,v and P v,l are the pressures inside and outside the liquid drop or vapour bubble. Using δV l,v = −δV v,l = 2(D − 1)πR D−1 δR and δΣ = 2(D − 1) 2 πR D−2 δR leads to a generalized Laplace law, The derivative of surface tension dσ/dR is termed a notional derivative by some authors (Blokhuis & Bedeaux 1992) in order to stress that it refers to arbitrariness of the dividing surface. Apart from the equimolar surface (3.20), the surface of tension is another possible choice to lift the ambiguity of the dividing surface. The notional derivative vanishes at the surface of tension, dσ dR R=Rs = 0, (3.23) thereby reducing the generalized Laplace law (3.22) to a standard form, Integrating (3.22) from R s to R, and eliminating ∆P using (3.24), one obtains analytic expression for the notional surface tension σ(R) relative to its minimum σ s at the surface of tension R s , Eq. (3.25) provides for a simple way of identifying the surface of tension and the corresponding surface tension. Two-dimensional simulations at T r = 0.98 were conducted with different initial drop and bubble sizes, R 0 ∈ {30, 40, 50, 60, 70, 80, 100, 120, 140}. For each dividing surface of radius R, the corresponding surface tension was evaluated as (Blokhuis & Bedeaux 1992), where P (r) is the tangential component of the pressure tensor, computed via the discrete pressure tensor detailed in Appendix E, and is the normal pressure component in the sharp interface system, where H is the Heaviside step function. For each drop or bubble, surface tension σ(R) (3.26) was probed at seventeen equidistant dividing circles between R min = 5δr and R max = 165δr. The discrete values σ(R) obtained in these simulations were fitted with Eq. (3.25), with σ s and R s as the free fitting parameters. Fig. 6 shows that the data for all drops and bubbles collapsed on a single master curve, in excellent agreement with the theory (3.25). Thus, the proposed model correctly identifies the surface of tension for the van der Waals fluid. While the notional derivative dσ(R)/dR vanishes at the surface of tension, the same does not hold for the surface tension at the surface of tension, dσ s /dR s = 0. In other words, surface tension σ s depends on the curvature of the surface of tension. In a seminal paper, Tolman (1949) characterized the curvature-dependence of surface tension by the Tolman length: For sufficiently large R s , the leading-order curvature-dependence of the surface tension may be written (Tolman 1949), where δ T is the Tolman length and σ 0 is the flat interface surface tension coefficient.
Here, the negative (positive) sign corresponds to drops (bubbles), respectively. We first compare the leading-order Tolman model (3.28) with the simulation. The values of R s obtained in the previously described drops and bubbles simulations of Fig.  6 are plotted against the pressure difference for different drop and bubble sizes in Fig. 7. It is clear that, for smaller drops and bubbles, the pressure difference deviates from the Laplace law with constant σ s = σ 0 , indicating a curvature-dependent surface tension.
Fitting the data points with Eq. (3.28), the Tolman length can be extracted from the simulation, here δ T = 9δr for both drops and bubbles, at the reduced temperature T r = 0.98. While the leading-order Tolman correction (3.28) improves agreement with data at moderate R s , deviation persists for smaller drops and bubbles at δr/R s > 0.03. Higherorder terms in the inverse powers of R s , important for droplets or bubbles of small radius, are neglected by (3.28) and were addressed by Helfrich (1978): where k andk are the bending and Gaussian rigidities; note that the latter vanishes for D = 2. Taking the second-order term (3.29) into account, the best fit in Fig. 7 results in bending rigidity k = 1.049 × 10 5 σ 0 δr 2 . Finally we consider the limit of flat interface where the Tolman length can be derived independently from the above considerations. In this case, location of the surface of tension X s can be found as the normalized first-order moment of the normal stress difference (Rao & Berne 1979), (3.30) With the dividing surface as a vertical straight line at X in two dimensions, the mass adsorbance Γ (X) is defined as (see example in Fig. 8), (3.31) Similar to the case of cylindrical symmetry considered in sec. 3.2, the equimolar surface is found by annihilating the mass adsorbance, Γ (X e ) = 0. The Tolman length in the limit of flat interface is the distance between the surface of tension and the equimolar surface (Blokhuis & Bedeaux 1992;Blokhuis & Kuipers 2006), (3.32) An example from the simulation is presented in Fig. 9. As shown by Blokhuis & Bedeaux (1992), the Tolman length scales with the reduced temperature as, In order to validate the scaling (3.33) in our model, flat interface simulations were conducted over a range of temperatures in the vicinity of the critical state. The Tolman length (3.32) for various temperatures is shown in Fig. 9. The results obtained from simulations are in excellent agreement with the theoretical scaling (3.33). A similar study using the Shan-Chen equations of states and the pseudo-potential model was presented in (Lulli et al. 2021). Furthermore, the flat interface simulations lead a to δ T = 9.2δr at T r = 0.98, in agreement with the value obtained from drop and bubble simulations, δ T = 9δr at same temperature. The relatively small discrepancy can be attributed to higher-order terms in the curvature, neglected in the fitting process.

Interface width: Temperature scaling and control
In the present work we use a definition for the interface width bearing numerical information as to how well the stiff gradients are resolved on a given mesh, making it directly related to the velocity increment per time-step: where ρ l and ρ v are densities of saturated liquid and vapor, respectively. It can readily be observed that in the limit of a sharp interface, i.e. resolved with δr, δr/W → 1. As previously noted for the co-existence densities and surface tension, the model recovers thermodynamical properties/scalings of the second-gradient fluid in the limit ε → 0, akin to δr/W → 0. As such in the limit of δr/W → 1 one expects numerical effects to dominate and observe deviations from thermodynamics of the van der Waals fluid. Surface tension vanishes as the temperature approaches the critical, cf. section 3.2, Eq. (3.14) and Fig. 4, while the interface diverges as T → T c . As noted by Widom (1965), the van der Waals theory predicts the temperature scaling of the interface width as, In order to validate consistency of the proposed model, simulations of flat interfaces were carried out in a range of reduced temperatures T r near the critical point, and corresponding interface widths W (T r ) (3.34) were measured. Given the effect of the choice of a on interface thickness simulations were also carried out with different a.
The results obtained from simulations near the critical point are shown in Fig. 10. As noted by many authors in the literature (Jamet et al. 2001), the parameter a can be used to control the interface thickness, as done in the present work, at a given density ratio, leaving the ratios and Maxwell construction unaffected. In agreement with the equivalent states theory, the right hand side plot in Fig. 10  a. In addition it is interesting to note that for a fixed grid-size δr, as (1 − T r ) → 1, δr/W → 1 (equivalent to the scaling parameter ε introduced in the multi-scale analysis) indicating deviation from the thermodymically converged state. This is illustrated by the deviation of the numerical interface thickness, starting at T r ≈ 0.98 from the theoretical predictions. Lowering the value of a, i.e. rescaling the interface by a factor 1/ √ a and therefor lowering ε, it is observed that interface is again well-resolved and the scaling (3.35) restored. In agreement with previous sections, in the limit of ε → 0, the model is shown to be thermodynamically converged and recovers the properties of the secondgradient fluid.

Shear stress: Layered Poiseuille flow
The setup consists of a rectangular domain filled with the liquid phase at the bottom and the vapor phase on top. The flow is driven by a body force. Top and bottom are subject to no-slip boundary conditions while the inlet and outlet are fixed by periodicity. The momentum balance at steady state reduces to a well-posed system of ordinary differential equations, ∂ y µ l ∂ y u + ρ l g = 0,∀y : 0 y h l , (4.1a) ∂ y µ v ∂ y u + ρ v g = 0,∀y : h l y H, (4.1b) closed by the following boundary conditions: u| y=H = 0, (4.2b) (4.2d ) Here H is the height of the channel, h l the height of the section filled with the liquid phase, µ l and µ v are the dynamic viscosities of the liquid and vapor, respectively, and g is the acceleration due to the body force. Analytical solution is given in Appendix F. This configuration is of particular interest as it probes ability of a two-phase model to capture jump conditions at the interface which happens only if deviatoric components of the viscous stress tensor (2.11) are correctly recovered. Two sets of parameters were considered: (a) µ l /µ v = 10.1, ρ l /ρ v = 10.1 (T r = 0.77) and (b) µ l /µ v = 11.3, ρ l /ρ v = 1030 (T r = 0.36). Simulations were conducted on a grid of size 5 × 200, with a constant acceleration g = 10 −8 δr/δt 2 . Furthermore, the same simulations were performed with the conventional second-order equilibrium (Succi 2002) instead of the product-form (2.54). The steady-state results are compared to analytical solution in Fig. 11. Contrary to the conventional second-order equilibrium, the use of the product form (2.54) in our model recovers the correct jump conditions and therefore results in continuous velocity profiles, also for a large density ratio, in excellent agreement with analytical solution.

Normal stress: Dissipation of acoustic waves
To assess Galilean invariance of the diagonal components of the viscous stress tensor and the effect of the correction term, dissipation of acoustic waves was measured numerically. Acoustic waves were initialized as a small perturbation of density with an amplitude δρ and subject to uniform background velocity U 0 representing a moving reference frame, where ρ 0 is the density of saturated liquid or vapor at a given temperature and λ is the wavelength of the perturbation. The maximum velocity in the domain was monitored over time and fitted to an exponential function of the form, max(u) − U 0 = exp (−Λt), where the coefficient Λ is tied to the dissipation rate of the normal modes as, Simulations were performed on a periodic domain of the size L = 128δr × δr with λ = 64δr at the temperature T r = 0.36, corresponding to density ratio ρ l /ρ v ≈ 10 3 . Perturbation amplitude was set to δρ = 10 −4 ρ 0 while the reference frame velocity varied in the range U 0 δt/δr ∈ [0, 0.3]. The dissipation rate (4.4) with and without the correction term are shown in Fig. 12. It is observed that the correction term for the third-order diagonal moments of the equilibrium populations (2.58) restores the Galilean invariance (independence from the reference frame velocity) of the dissipation rate of normal modes.

Viscosity-capillarity coupling: Oscillating drop
Rayleigh's classical study of a freely oscillating drop of non-viscous liquid under small deformations laid ground for the understanding of capillary waves (Rayleigh 1879). Assuming purely axisymmetric oscillation modes, Rayleigh's oscillation frequency of nth mode, n = 2, 3, . . . , for large density ratios, ρ l /ρ v 1, is given by, (4.5) Two-dimensional simulations were performed for a drop of initial radius R 0 = 40δr, for a relatively low kinematic viscosity, ν l δt/δr 2 = 0.01. Modes of order n = 2, 3, 4, 5, 6 were initiated with a monochromatic perturbation, R(θ, 0) = R 0 1 + a n cos(nθ) − 1 4 na 2 n , with the amplitude a n = 0.1. Reduced temperature of van der Waals fluid was set to T r = 0.36, corresponding to density ratio ρ l /ρ v ≈ 10 3 . The initial and mid-cycle shapes of Rayleigh's modes are shown in Fig. 13. Simulations were performed over eight oscillation cycles and corresponding oscillation periods were identified. Fig. 13 demonstrates excellent comparison of oscillation frequencies measured in the simulation with corresponding theoretical values (4.5), even for higher-order modes. While Rayleigh's analysis was restricted to non-viscous liquids, a number of analytical solutions for viscous oscillating drops have been obtained over the years. In particular, Aalilija et al. (2020) derived a time-dependent solution of the following form, where n = a n exp(−λ n t) cos (2πf n t), while λ n is the damping rate of the nth mode, In order to validate the damping rate of capillary waves, the previously described twodimensional drop was simulated in a range of kinematic viscosities, ν l δt/δr 2 ∈ [0.01, 0.1], subject to a monochromatic perturbation with the lowest mode n = 2 and the initial perturbation amplitude was set to a 2 = 0.2. The envelope of oscillations was fitted using the exponential function exp(−Λt), see Fig. 14. The effective damping coefficient Λ was used along with Eq. (4.7) to evaluate the apparent kinematic viscosity. As shown in Fig. 14, the dissipation rates obtained by simulation agree well with the analytical solution (4.7).

Isothermal speed of sound
Finally, we validate the speed of sound at different temperatures in both the liquid and vapor phases. For the van der Waals fluid, the isothermal sound speed is, The sound speed was measured by monitoring the position of a pressure front over time in a quasi-one-dimensional simulation at different temperatures. The obtained results are compared to theory (4.8) in Fig. 15. It is observed that the simulations accurately capture the speed of sound in both liquid and vapor phases, for entire range of density ratios on the coexistence diagram Fig. 1, up to at least ρ l /ρ v ∼ 10 11 . Results for other equations of state are compared in Appendix G.

Interaction with solid boundaries: Equilibrium contact angle
A comprehensive review of various implementations of the contact angle on solid boundaries in the lattice Boltzmann setting can be found in Li et al. (2014). Here we follow a proposal by Benzi et al. (2006) where a virtual density, and therefor corresponding virtual pseudo-potential, is attributed to the solid nodes. The calculations of the force (2.69) is then carried at all fluid nodes. The no-slip condition is imposed via the modified bounce-back scheme of Bouzidi et al. (2001) for curved boundaries. The virtual density attributed to the solid nodes is the free parameters controlling the contact angle. The ability of this approach to correctly reproduce the dynamics of the contact line in a high-density ratio regime will be validated below.

Young-Laplace equation
As a first validation of the solid/fluid interaction the case of a two-dimensional channel of height H and length L is considered. Initially a rectangular column of liquid of length L/3 and height H surrounded by vapor was placed at the center of the channel. Simulations were performed until the system reaches steady-state, and the angle between the liquid/vapor interface and the solid wall at the triple-point was measured. The static contact angle measured at the triple-point should verify the Young-Laplace law for this configuration: where θ is the equilibrium contact angle and σ the liquid/vapor surface tension.
Here to validate the static angle, the reduced temperature of van der Waals fluid was set to T r = 0.36 which corresponds to density ratio ρ l /ρ v = 10 3 . The pressure difference ∆P is computed as the difference in pressure between two monitoring points located respectively inside and outside the liquid column. The convergence criterion used to assess steady-state is based on the time-variations of pressure at these two monitoring points. The contact angles were measured using the low-bond axisymmetric drop shape analysis module (Stalder et al. 2010) in ImageJ. The results obtained for the entire contact angle range are shown in Fig. 16. The measured contact angles show excellent agreement with the Young-Laplace equation. Furthermore, the results show that, for ρ w → ρ v , the contact angle tends to θ → 180 • , while as ρ w → ρ l , the contact angles tends to θ → 0 • therefore covering the entire range of contact angles.

Wettability-driven drop motion in a channel
Difference in contact angle on opposite sides of a liquid drop placed on a surface with a non-uniform wettability may initiate a motion of the drop. As a further validation of the solid/liquid interaction of the proposed model, the case of a liquid column placed in a channel and subjected to a wettability step function is considered. The configuration, illustrated in Fig. 17, consists of a channel of length L and width/height of H. Initially a liquid column of length L l is placed at the channel center. Once the liquid column has reached equilibrium, the wettability on the right-half of the channel is changed, which results in unequal contact angles on the different sides of the liquid column. An approximate analytical solution for the center-of-mass position X and centroid velocity U was found by Esmaili et al. (2012), (4.10) where the saturated centroid velocity U ∞ and transition time τ ∞ are computed respectively as: , (4.12) . (4.13) Here H = 70δr while L = 1260δr, L l = 420δr, ν l = ν v = 0.075δr 2 /δt and density ratio is set to 10 3 as in previous cases. Initially the contact angles on both sides of the water column were set to θ − = 59.3 • . Once a steady-state was reached, the contact angle on the right-hand side of the column was changed to θ + = 63.4 • , which resulted in a net force acting on the liquid. The liquid column centroid velocity and center-of-mass position are compared to analytical solutions in Fig. 18. Both quantities closely match analytical predictions.

Thermodynamic convergence and resolution requirements
Two-phase lattice Boltzmann models are routinely used to simulate impact of liquid drops on solids. While experiments are typically conducted with water and other liquids in air, the applicability of the two-phase LBM is justified whenever the dynamic effects are concerned. However, majority of LB models do not reach the experimentally relevant density ratios of, say, water/air density. The proposed LBM was demonstrated to be valid and thermodynamically well-posed at high density ratios. However, thermodynamic convergence, especially at higher density ratios comes at a cost in terms of interface resolution requirements. For instance, at T r = 0.36 resulting in a density ratio of 10 3 , to guarantee deviation below one percent from theory in the vapor phase one must have W > 8δr. In large-scale hydrodynamic-dominated cases targeted in the present section full thermodynamic convergence is not necessary. The only manifestations of deviation from the thermodynamically converged state that must be kept under control are spurious currents. Here, for density ratios 10 3 the interface thickness, (3.34), is fixed at W = 5δr resulting in a seven percent deviation of the vapor phase density from theory and spurious currents below 10 −3 δr/δt. The Tolman length for this choice of parameter is δ T /δr = 2.4 which for the resolutions considered below (R 0 = 75δr) results in δ T /R 0 = 0.032 therefor minimizing the impact of curvature on effective surface tension. Below, we consider benchmark simulations of dynamic effects at realistic density ratios of increasing complexity to probe the accuracy and numerical stability of our model in the dynamic setting.

Contact time on flat non-wetting surfaces
Extensive studies of drop impact dynamics have shown that the contact time on socalled super-hydrophobic surfaces is independent of the Weber number, We = ρ l D 0 U 2 0 /σ, and only scales with the inertio-capillary time, τ i = ρ l D 3 0 /8σ, i.e. for a given drop initial diameter D 0 the contact time is not affected by impact velocity (Gauthier et al. 2015). The different stages of the impact process are illustrated in Fig. 19 through water drops impacting two surfaces with different contact angles. Simulations carried out under the same conditions are shown and already point to very good agreement between simulations and experiments. To quantify the ability of the proposed model to capture the dynamics of drop impact and to investigate the Weber-independence of the contact time on nonwetting surfaces, simulations were performed for a wide range of Weber numbers We ∈ [1, 40], resulting in Reynolds numbers in the range of Re ∈ [400, 1400], with Re=U 0 D 0 /ν l . Simulations have been performed in boxes of size 4D 0 × 4D 0 × 4D 0 with D 0 = 150δr. The temperature was set to T r = 0.36 resulting in a density ratio of 10 3 . The obtained data show very good agreement with experimental contact times measured by Gauthier et al. (2015) and confirm the Weber-independence of the contact times. Both numerical and experimental data are shown in Fig. 20.

Reducing the contact time: Pancake bouncing
To further reduce the drop contact times, a number of different strategies have been devised during the past decades. Recently, Liu et al. (2014) proposed to use macroscopic structures, in the form of tapered posts, to reduce the contact time. It has been shown that above a certain threshold Weber number these structures can decrease the contact time by approximately 75 percent. This mechanism is also known as pancake bouncing, due to the pancake-like shape of the drop at take off. A detailed numerical study of pancake bouncing using LBM has been is presented in (Mazloomi Moqaddam et al. 2017), for a limited range of densities. Here, to further demonstrate the versatility of our model we consider a realistic density ratio of 10 3 . The geometry consists of a flat substrate populated by tapered posts of base and tip radii R B and R b and height h with a center-to-center distance w in a simple square arrangement. The simulations were run following the the geometrical configurations considered in experiments (Liu et al. 2014). Curved wall boundaries were implemented via the modified bounce-back scheme introduced by Bouzidi et al. (2001).The geometry is illustrated in Fig. 21. As for the flat substrate simulations, the domain size was set to 4D 0 × 4D 0 × 4D 0 with D 0 = 150δr. Snapshots from the different stages of impact for two different Weber numbers, one below the pancake bouncing threshold and one above, from both simulations and experiments are shown in Fig. 22. The results show excellent agreement between simulations and experiments. The contact times for different Weber numbers as obtained from simulations are compared to experimental results from (Liu et al. 2014) in Fig. 23. It is shown that the numerical simulations not only accurately capture the contact time reduction due to pancake bouncing and the shape of the drop at take off, they also correctly predict the onset of pancake bouncing. The shape of the drop at take off is assessed via the pancake quality parameter Q defined as the ratio of the drop radius at take off to its radius at maximum spreading (Liu et al. 2014).

Extreme density ratios: Inertia-dominated coalescence of mercury drops
The sudden and pronounced topological changes involved in the coalescence of drops, especially the formation and subsequent evolution of the neck between them has been of interest and subject to study for many years. The dynamics of the dimensionless neck radius r neck /R 0 (with R 0 the drops initial radius and r neck the neck radius) have been shown to belong to one of two regimes (Eggers et al. 1999) characterized by the Ohnesorge number, Oh = µ/ √ 2ρ l σR 0 : highly viscous Stokes or inertial. The neck radius evolution over time has been shown to scale either as r neck /R 0 ∝ t/τ i (in the viscous regime) or r neck /R 0 ∝ t/τ i (in the inertial regime). Here, to better illustrate the ability of the proposed scheme to model flows with extreme density ratios, we will focus on the first regime, more specifically the coalescence of mercury droplets, with Ohnesorge numbers of the order of 10 −4 (for 1 g drops). We consider a case following the experimental configuration of Menchaca-Rocha et al. (2001) with two 1 g mercury droplets. To match the proper density ratio of a mercury/air system (ρ Hg = 13600 kg/m 3 ) the temperature is set to T r = 0.267, resulting in a density ratio of ρ l /ρ v ∼ 11500. The two drops, resolved with 150 points on their radii are placed in a rectangular domain of size 400 × 400 × 800 with a center-to-center distance of 304 points. The drops are connected initially via a neck of 4 points radius. The evolution of the drops shape over time are compared to experiments from Menchaca-Rocha et al. (2001) in Fig. 24. To further illustrate the agreement of numerical simulations with experiments, the evolution of the neck radii r neck over time are shown in Fig. 25. Both qualitative comparison of the drops shape and quantitative comparison of the neck radius show that the presented solver is able to correctly model multi-phase flow dynamics with extreme density ratio. To the authors knowledge this is the first lattice Boltzmann simulation at such high density ratios, regardless of the formulation, i.e. pseudo-potential, free energy, phase-field etc.

Conclusion
Multi-phase flows are a well-established area in kinetic theory and more specifically in the context of the LBM. While dynamic 2-D/3-D simulations are routinely carried out with the different LB-based formulations, which are subject to continuous numerical improvements, a clear and concise kinetic framework along with a characterisation of the resulting fluid thermodynamic properties is lacking. The aim of the present work was to develop a general framework for iso-thermal singlecomponent multi-phase flow simulations consistent with the capillary fluid thermodynamics. Starting from the first-order BBGKY equation and using a projector operator in phase-space, a flexible model, in terms of pressure contribution partition, is proposed. A specific realization of the partition minimizing deviations from the first-neighbor stencil reference (optimal) state was then discretized using the LBM. The discrete solver was shown to recover the full Navier-Stokes-Korteweg under the proposed scaling. The resulting discrete model was then probed for thermodynamic consistency and shown to abide by all the scaling laws of the corresponding meanfield, i.e. surface tension, interface thickness and Tolman length. Finally the model was shown to allow for not only static, but dynamic large density ratios simulations with complex geometries, illustrated best by the extreme case of mercury drops coalescence, therefor effectively removing the wellknown limitations on maximum density ratio.
The detailed study of the properties of the proposed model and accompanying theoretical analyses lead us to believe that the present framework fills important gaps in the multi-phase literature, more specifically in the context of discrete kinetic models such as the LBM. It provides a consistent framework for the simulation of challenging physics as demonstrated by the curvature-dependence of the surface tension. Furthermore, it paves the way for extension to fully compressible non-ideal fluids, which will be topic of future publications.
while J (1) E is the non-local contribution to the lowest order, Since the Boltzmann collision integral conserves the mass and momentum locally, it is annulled by the projector, Furthermore, J 1 SR is evaluated at the local equilibrium f eq (2.16) to get (Chapman & Cowling 1939), which, after integration in v 1 and k, for the isothermal flow results in, where b = 2πd 3 /3m. Finally, applying the iso-thermal projector K, we obtain While the phenomenological Enskog's collision integral (Enskog 1921) was used above, the lowest-order approximation (A 9) is identical in other versions of hard-sphere kinetic equations such as the revised Enskog theory (RET) (Van Beijeren & Ernst 1973) or kinetic variational theory (Karkheck & Stell 1981).

Appendix B. Long-range interaction
Vlasov's mean-field approximation for a long-range interaction is reviewed next. Assuming absence of correlations, the two-particle distribution function is approximated as while the long-range interaction integral can be simplified, With a Taylor expansion around r, and neglecting higher-order terms, Eq. (B 2) leads to, where parameters a and κ are, after integration over a unit sphere, Applying the projector, we obtain, Finally, we can estimate the relative magnitude of the parameters a and κ by introducing a range of the attraction potential δ. Assuming d δ, we have a ∼V δ 3 and κ ∼V δ 5 , whereV is a characteristic value of the potential, thus κ/a ∼ δ.
Appendix C. Hydrodynamic limit of the Enskog-Vlasov-BGK kinetic model

C.1. Rescaled kinetic equation
For the Enskog-Vlasov-BGK kinetic model, let us introduce the following parameters: characteristic flow velocity U, characteristic flow scale L, characteristic flow time T = L/U, characteristic densityρ, isothermal speed of sound of ideal gas c s = √ RT , and kinematic viscosity of the BGK model of ideal gas ν = τ c 2 s . With the above, the variables are reduced as follows (primes denote nondimensional variables): time t = T t , space r = Lr , flow velocity u = Uu , particle velocity v = c s v , density ρ=ρρ , distribution function f =ρc −3 s f . Furthermore, the following non-dimensional groups are introduced: Viscosity-based Knudsen number Kn = τ c s /L, Mach number Ma = U/c s , Enskog number En = bρKn/Ma and Vlasov number Vs = a/bRT . With this, the Enskog-Vlasov-BGK kinetic model is rescaled as follows: where δ is the range of the attraction potential, estimated according to Eq. (B 8). The following scaling assumptions are applied: Acoustic scaling, Ma ∼ 1; Hydrodynamic scaling, Kn ∼ En ∼ δ/L ∼ ; Enskog-Vlasov parity, Vs ∼ 1. In other words, the conventional hydrodynamic limit treats all non-dimensional groups that are inversely proportional to the flow scale L (Kn, En and δ/L) as a small parameter while the Enskog-Vlasov parity ensures that both the short-and long-range contributions to the pressure are treated on equal footing. Returning to dimensional variables, we may write, where, Finally, taking into account the reference equilibrium (2.34) and the reference pressure P 0 , the rescaled kinetic equation (2.35) takes the form (C 3) with where we have only retained terms up to order two. Then introducing characteristic flow size L and velocity U the equation is made non-dimensional as: δr where primed variables denote non-dimensional form and where c = δr/δt. Assuming acoustic, i.e. U c ∼ 1 and hydrodynamic, i.e. δr L ∼ δu U ∼ ε, scaling and dropping the primes for the sake of readability: Then introducing multi-scale expansions: the following equations are recovered at scales ε and ε 2 : The moments of the non-local contributions (including both non-ideal contributions to the thermodynamic pressure, surface tension and the correction for the diagonals of the third-order moments tensor) are: Taking the moments of the Chapman-Enskog-expanded equation at order ε:

∂
(1) t ρu + ∇ · ρu ⊗ u + ∇ · P 0 I + F = 0, (D 11) while at order ε 2 the continuity equation is: Summing up Eqs. D 10 and D 12 we recover the continuity equation as: Furthermore, using U = u + δt 2ρ F and: (D 23) in combination with the Euler-level equation, and keeping in mind that errors of the form ∇ · δt 2 F ⊗F 4ρ in the convective term and δt∇µ ∇ F ρ + ∇ F ρ † in the viscous stress are of order ε 3 one recovers: Appendix E. Discrete non-local contributions to pressure tensor Following the analysis presented by Shan (2008), we write the force contribution for the present model as, The pressure tensor contributions from forces F A and F C can be readily written as: while F B and F D contribute to the pressure tensor as follows: w i ς 2 c i ⊗ c i ρ(r − c i δt)ρ(r + c i δt) .

(E 9)
These expressions allow to compute the discrete pressure tensor with high accuracy. computed from both the discrete and continuous pressure tensors, While the discrete evaluation method correctly results in a uniform pressure distribution throughout the domain, also across the interface, the continuous approximation evaluated using a finite differences approximation fails to do so, indicating errors due to higherorder terms. This points to the necessity of using the discrete pressure tensor instead of Eq. (E 10) for evaluation of quantities such as surface tension and Tolman length in sec Starling (3.10), along with two equations of state proposed by Shan & Chen (1993), with: while P 0 = ς 2 ρ. Corresponding critical densities ρ c and pseudo-temperatures G c are readily computed by solving the conditions at critical point, ∂P SC ∂ρ Gc,ρc = 0, ∂ 2 P SC ∂ρ 2 Gc,ρc = 0, (G 4) which results in (G c = −4, ρ c = 0.6931) and (G c = −7.389, ρ c = 1) for the pseudopotentials (G 2) and (G 3), respectively. Fig. 27 shows excellent agreement between simulations using present model and theory for all equations of state. It is interesting to note that, contrary to the conventional equations of state of the van der Waals type, equations of state (G 1), (G 2) and (G 3) demonstrate higher speed of sound in the vapor phase than in the liquid. While this, by itself, does not contradict thermodynamics, it remains unclear which substances may feature such an unusual behaviour.