# Adjustable surface tension independent of the collision operator for pseudo-potential lattice Boltzmann methods

Romain Noël, Feifei Qin, Linlin Fei, Jan Carmeliet

► To cite this version:

Romain Noël, Feifei Qin, Linlin Fei, Jan Carmeliet. Adjustable surface tension independent of the collision operator for pseudo-potential lattice Boltzmann methods. Physical Review Fluids, 2025, 10 (10), pp.104901. 10.1103/lrz9-ll2v. hal-05375568

HAL Id: hal-05375568

https://hal.science/hal-05375568v1

Submitted on 20 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Adjustable surface tension independently of the collision operator for
pseudo-potential lattice Boltzmann methods

Romain NOËL

Univ. Gustave Eiffel, Inria, Cosys/SII,
I4S, Bouguenais, F-44344, France and

Chair of Building Physics, Department of Mechanical and Process Engineering,
ETH Zürich (Swiss Federal Institute of Technology in Zürich), Zürich, 8093, Switzerland

Feifei QIN*

Institute of Extreme Mechanics, School of Aeronautics,
National Key Laboratory of Aircraft Configuration Design,
Key Laboratory for Extreme Mechanics of Aircraft
of Ministry of Industry and Information Technology,

Northwestern Polytechnical University,
Xi’an, 710072, Shannxi, China. and

Chair of Building Physics, Department of Mechanical and Process Engineering,
ETH Zürich (Swiss Federal Institute of Technology in Zürich), Zürich, 8093, Switzerland

Linlin FEI and Jan CARMELIET

Chair of Building Physics, Department of Mechanical and Process Engineering,
ETH Zürich (Swiss Federal Institute of Technology in Zürich), Zürich, 8093, Switzerland

(Dated: February 25, 2025)
Abstract

In this work, we propose an alternative surface tension adjustment approach in the pseudo-potential lattice BOLTZMANN (LB) model, which is not only decoupled from the density ratio but also independent of the collision operator. This is achieved by incorporating a generic source term obtained from the difference between a modified moment equilibrium and an original moment equilibrium distribution function in the LB equation. The explicit form of the source term is obtained through a third-order CHAPMAN–ENSKOG analysis of the LB equation, aiming to recover the targeted governing equations from a modified LANDAU free energy theory using a hybrid density-pseudo-potential form. The source term can be easily and straightforwardly incorporated into different widely used collision operators, such as single relaxation time (SRT or LBGK), multiple relaxation time (MRT) and entropic-MRT (KBC) operators. The proposed method is validated by three benchmarks: (i) the liquid-gas co-existence density, (ii) surface tension adjustment via the static droplet case and (iii) droplet deformation via oscillations. It is shown that thermodynamic consistency is restored in a large range of temperature ratio, and, moreover, a remarkable tunable surface tension range of 140 times can be achieved. Benefiting from the proposed surface tension adjustment method, we have successfully modeled the droplet impact and splashing dynamics with a WEBER number up to 10500, achieving one order of magnitude higher than LB simulations reported in the literature.

# I. INTRODUCTION

## A. Multiphase flows

Multiphase flows are widely observed in natural phenomena, daily life, and various engineering applications, ranging from raindrop dynamics, and water boiling process in cooking to inkjet printing, to name but a few [1, 2]. Thus, gaining a thorough understanding of the underlying physical mechanisms of complex multiphase interface dynamics is of key importance to the advancement in many fields of science and technology. However, multiphase flow is a typical non-linear, non-equilibrium and multi-scale process, involving complicated motions and deformations of thin interfaces (between phases or components) with complex shapes at varying time and space scales. For such a process, it is very challenging to experimentally quantify the key dynamic parameters at play, hampering the detailed analysis of the mechanisms behind. Therefore, there is an increasing need for developing robust and reliable direct numerical models to capture the complex interfacial dynamics in multiphase flow.

During the last three decades, the Lattice BOLTZMANN Method (LBM) [3–5], based on simpli-

<footnote>* fqin@nwpu.edu.cn</footnote>
fied lattice formulations of the kinetic BOLTZMANN equation, has obtained remarkable success for numerical simulations of various complex multiphase flows [6, 7]. The classical “collision-streaming” algorithm decouples the non-linearity and non-locality (which are often tangled in NAVIER-STOKES equations-based numerical models) [8]: the non-linear collision operator is entirely local and the non-local streaming is completely linear towards the discrete distribution, enabling it highly efficient in large-scale parallel computing [9, 10]. The mesoscale nature of LBM allows its easy incorporation of micro- and mesoscale physics, leading to the handy treatment of multiphase interface dynamics. Therefore, shortly after the birth of LBM, it had been extended to simulate multiphase/multicomponent flows [11]. Up to now, the existing multiphase Lattice BOLTZMANN (LB) models can be generally classified into four categories: the color-gradient model [11, 12], the pseudo-potential model [13, 14], the free-energy model [15, 16] and the mean-field model [17]. Among them, the pseudo-potential model (also known as SHAN-CHEN model), originally proposed by SHAN and CHEN [13, 14], is the simplest one in both concept and implementation and thus becomes particularly popular for LB simulations of multiphase flows.

In the pseudo-potential model, the non-ideal interactions among phases/components are mimicked by means of a force term, defined based on a local pseudo-potential function of the fluid density [13]. By proper choice of the function, the main features of multiphase flows, namely a non-ideal gas Equation of State (EoS) supporting the phase transition and a non-zero surface tension to minimize the interfacial energy, can be incorporated into the single-phase LB model with minimum programming effort [8]. With this treatment, the liquid-gas interface is automatically captured and dynamically evolves during the simulation, advantageous than other free-energy models that require extra ALLEN-CAHN [18] or CAHN-HILLARD [19] equations to track the interface. However, simplicity comes at the expense of some well-known drawbacks of the original pseudo-potential model. One drawback is that a spurious velocity field exists near the phase or component interfaces [20], which not only decreases the numerical stability but also contaminates the real velocity field. Another drawback is the thermodynamic inconsistency, in the sense that the liquid-vapor coexistence densities by the mechanical stability solution from the pseudo-potential force deviate from those of the MAXWELL's construction from the EoS [14, 21], limiting the model to multiphase systems to relatively small density ratios and unrealistic equation of state. By improving the original pseudo-potential model from many perspectives, for example, the calculation of the pseudo-potential force [10, 21–23], improvement in the forcing scheme [24–28], incorporation of different/modified EoS [29–32], and use of advanced collision operators [33–36], the above two drawbacks have been overcome significantly.
In addition, the pseudo-potential model suffers from the trouble of non-adjustable surface tension, which also obtained much attention. In the original pseudo-potential model, the surface tension, which plays an essential role in interfacial phenomena, depends on the density ratio (and the adopted EoS). In order to adjust the surface tension, SBRAGAGLIA et al. [22] developed a multi-range potential by combining the original nearest-neighbor interactions and the next-nearest-neighbor interactions. They claimed that the surface tension could be tuned independently of the density ratio by carefully tuning the weight coefficients in the multirange potential [22]. However, HUANG et al. [37] found that the density ratio changes when using the multi-range model, especially for relatively large density-ratio cases. By adding a source term into the Multiple Relaxation Time (MRT) LB equation to incorporate specific additional terms in the pressure tensor, LI et al. [38] proposed a nearest-neighbor-based method to tune the surface tension. It is demonstrated based on theoretical analysis and numerical results that the method can change the surface tension over a wide range while keeping the density ratio almost unchanged. The method has been extended to three-dimensional (3D) systems based on the classical MRT LBM by XU et al. [39] and non-orthogonal MRT LBM by FEI et al. [35]. It has also been further improved by considering the effects of higher-order terms based on third-order and even fourth-order analysis [27, 28, 33, 40]. Nevertheless, the methodology of LI et al. is mainly used in MRT-type collision operators, whose extension to other collision operators may need additional efforts. Theoretically, the surface tension could also be changed independently of the density ratio via the customized EoS [32], while its numerical verification, especially in unstable cases, is still missing. More recently, WANG et al. [10] proposed an extended combined pseudo-potential model based on the unified lattice BOLTZMANN framework [41, 42].

The motivation of the present work is to propose an alternative surface tension adjustment, which is not only decoupled from the density ratio but also independent of the collision operator. This is achieved based on a set of targeted governing equations that incorporate the thermodynamic pressure tensor from a modified free-energy functional. The rest of the paper is structured as follows: after this multiphase flows in LBM introduction, the following section I B is presenting the governing equations. Section II introduces the LB numerical for multiphase flows from a usual yet relatively general point of view. Section III details the new source term designed from a theoretical analysis and its numerical implementation. Section IV presents the assessment of the suggested models and conducts model validations against benchmark results. Finally, a brief conclusion is drawn in section V. Also, the CHAPMAN-ENSKOG expansion performed is detailed in section A, while section B gives extra results obtained for other collision operators and forcing schemes to
demonstrate the genericity of our work.

## B. Governing equations

As known, the motion of single-phase viscous fluid can be well described by the NAVIER-STOKES (NS) equations. To model multi-phase systems, the original NS equations need to be modified with the considerations of two key factors: (i) a correct relation that separates liquid from gas state, i.e., the EoS and (ii) the surface tension effect which determines the dynamics at the vicinity of the liquid-gas interface. Such an aim can be achieved by introducing a sophisticated pressure tensor into the NS equations.

To build a thermodynamic consistent pressure tensor, a common choice is using the LANDAU free energy theory. Following the classical capillarity theory of VAN DER WAALS, the free energy functional within a gradient-squared approximation is [43]:

$$ \mathcal{F}[\rho] = \int_V \Psi(\rho, \nabla \rho) d\boldsymbol{x} = \int_V \left( E_b(\rho, T) + \frac{\kappa}{2} (\nabla \rho)^2 \right) d\boldsymbol{x} \quad (1) $$

where $\Psi$ is the bulk free-energy of LANDAU, $\rho$ is the local density, $E_b$ is the bulk free-energy density at a given temperature $T$ and $\kappa$ is the capillarity coefficient. The square of the gradient term describes the free energy contribution from density gradients in an in-homogeneous system. The equilibrium state of the system is obtained by minimizing the free energy under the constraint of the constant total mass, leading to the non-local pressure expressed as [44]:

$$ p_{nl} = p_{EoS} - \kappa \rho \nabla^2 \rho - \frac{\kappa}{2} (\nabla \rho)^2 \quad (2) $$

where $p_{EoS}$ is the EoS given by

$$ p_{EoS} = \rho \frac{\partial E_b}{\partial \rho} - E_b \quad (3) $$

This leads to the following thermodynamic pressure tensor (also known as the KORTEWEG's stress tensor) [45]:

$$ p_{th} = p_{nl} I + \frac{\partial \Psi}{\partial \nabla(\rho)} \nabla(\rho) = \left( p_{EoS} - \kappa \rho \nabla^2 \rho - \frac{\kappa}{2} (\nabla \rho)^2 \right) I + \kappa \nabla(\rho) \otimes \nabla(\rho), \quad (4) $$

where $I$ is the identity matrix, and $\cdot \otimes \cdot$ is the tensor product. Incorporating the above pressure tensor into the NS equations leads to the target equations to be reconstructed by the free-energy-based kinetic models. It may be noted that in the original free-energy lattice BOLTZMANN models [15, 16], the pressure tensor is directly incorporated into the Equilibrium Distribution Function (EDF), which
destroys its Maxwellian shape, leading to worse control on the numerical stability. To improve the model, one can add the pressure tensor by means of an interaction force, either in a pressure form [17] or a chemical potential form [46].

Alternatively, one can address the problem by using a different free-energy density. The classical formulation $\Psi(\rho, \nabla \rho)$ in eq. (1), despite its manifest link to underlying molecular interactions [47], is yet to be checked by first-principle derivation. Therefore, it leaves us with some freedom to design a modified free-energy functional based on a bottom-up strategy. Among others, a hybrid density-pseudo-potential form $(\rho - \psi)$ can be considered [48]:

$$ \Psi(\rho, \nabla \psi) = E_b(\rho, T) + \frac{\kappa}{2} |\nabla \psi|^2. \qquad (5) $$

It has been confirmed that the above form is compliant with the free energy approach, supporting minimum principle [48]. Following the same procedure from eqs. (1) to (4), the modified KORTEWEG's stress tensor using the density-pseudo-potential form can be expressed as:

$$ p_{th} = \left( p_{EoS} - \kappa \rho \frac{d\psi}{d\rho} \nabla^2 \psi - \frac{\kappa}{2} |\nabla \psi|^2 \right) I + \kappa \nabla(\psi) \otimes \nabla(\psi) \qquad (6) $$

Re-injecting the expression of the stress tensor from eq. (6) in the NAVIER-STOKES equations, and using $(p_{EoS} = \rho T - \frac{G}{2}\psi^2)$ one obtains:

$$ \partial_t(\rho) + \partial_x(\rho v^{eq}) = 0, \qquad (7) $$

$$ \begin{aligned} & \partial_t (\rho v^{eq}) + \nabla_x \cdot \left( \rho v^{eq} \otimes v^{eq} - \nu \left[ \partial_x (v^{eq}) + \partial_x (v^{eq})^T \right] \right) \\ & + \nabla_x \cdot \left( \left\{ \rho T - \frac{G}{2}\psi^2 - \kappa \rho \frac{d\psi}{d\rho} \nabla^2 \psi - \frac{\kappa}{2} |\nabla \psi|^2 \right\} I + \kappa \nabla(\psi) \otimes \nabla(\psi) \right) = \rho g. \end{aligned} \qquad (8) $$

The above NS equations with the modified KORTEWEG stress tensor are the targeted physical governing equations for kinetic modeling in the present work.

## II. NUMERICAL METHOD: LATTICE BOLTZMANN

### A. Classical pseudo-potential LBM for multiphase flows

To realize the NAVIER-STOKES equations with additional force term, the following LB equation is generally applied

$$ f_i(\boldsymbol{x} + \boldsymbol{\xi}_i \Delta t, t + \Delta t) = f_i(\boldsymbol{x}, t) + \Omega_i + F_i, \qquad (9) $$

where $f_i$ is the discrete density distribution function, $\boldsymbol{\xi}_i = (\xi_{ix}, \xi_{iy})$ is the discrete velocity in the $i^{th}$ direction, $\Delta t$ is the time step, $F_i$ represents the forcing term in the velocity space. $\Omega_i$
is the collision operator depending on different models. For instance, in the widely-used Lattice BOLTZMANN, BHATNAGAR, GROSS and KROOK (LBGK) model [4],

$$ \Omega_{i,\text{BGK}} = -\frac{1}{\tau}(f_i - f_i^{eq}), \qquad (10) $$

where $\tau$ represents the relaxation time determined by the fluid viscosity, and $f_i^{eq}$ is the discrete equilibrium density distribution function. In the MRT model [49],

$$ \Omega_{i,\text{MRT}} = -(M_M^{-1}\Lambda_M M_M)_{ij}(f_j - f_j^{eq}), \qquad (11) $$

where $\Lambda_M = (\tau_\rho^{-1}, \tau_e^{-1}, \tau_\zeta^{-1}, \tau_j^{-1}, \tau_q^{-1}, \tau_j^{-1}, \tau_q^{-1}, \tau_v^{-1}, \tau_v^{-1})$ represents the diagonal matrix, and its elements represent the relaxation in moment space. $M_M$ is the orthogonal transformation matrix to transfer variables from discrete space to moment space. In the KARLIN-BOSCH-CHIKATAMARLA (KBC) model [36, 50],

$$ \Omega_{i,\text{KBC}} = -\beta(f_i - f_i^{\text{mirr}}), \qquad (12) $$

where $\beta$ is the parameter determining the fluid viscosity and $f_i^{\text{mirr}}$ the mirror state of the $f_i$. When $f_i^{\text{mirr}}$ is given as $2f_i^{eq} - f_i$, the KBC model degenerates into the LBGK model.

For the implementation of the force term, there are two widely-used schemes in the LB community, i.e., the Exact Difference Method (EDM) [24, 51] and GUO's scheme [52]. Both schemes are applicable to LBGK and KBC collision schemes. However, the liquid-vapor density ratio is found limited to a small value when using the original GUO's scheme, due to the discretization of LBM. Li et al. [25] improved GUO's scheme to realize a high density ratio up to 1000, and further extended it to the MRT collision operator [26] for lower fluid viscosity. For the LBM users, the generally applied models for the study of multiphase flows are, what we denoted in this work, “LBGK/KBC+EDM” and “MRT+GUO” combinations, where the former represents the LBGK or KBC collision operator coupled with the EDM force scheme and the latter represents the MRT collision operator coupled with GUO's force scheme with Li's improvement, respectively. In the LBGK/KBC+EDM framework, the force terms is written as

$$ F_{i,\text{EDM}} = f_i^{eq}(\rho, \mathbf{v} + \Delta\mathbf{v}) - f_i(\rho, \mathbf{v}), \qquad (13) $$

where $\mathbf{v} = \sum_i f_i \xi_i / \sum_i f_i$ and $\Delta\mathbf{v} = \mathbf{F} \Delta t / \rho$. $\mathbf{F}$ is the force term accounting for all external forces such as fluid-fluid interaction, fluid-solid interaction and body force. In the MRT+GUO framework, the forces are normally represented in the moment space. In current work, the simulations are conducted using the MRT+GUO framework, to realize large ranges of viscosity and surface tension. Therefore, we introduce the framework in detail as follows.
By the multiplication of the transformation matrix, the right-hand side of eq. (9) with MRT collision operator is written as

$$ \mathbf{m}^c = \mathbf{m} - \Lambda_M (\mathbf{m} - \mathbf{m}^{eq}) + \Delta t (I - \frac{\Lambda_M}{2}) \hat{\mathbf{F}}, \quad (14) $$

where $\mathbf{m} = M_M \mathbf{m}$, $I$ is the unit tensor, $\hat{\mathbf{F}}$ is the forcing term in the moment space with $(I - \frac{\Lambda_M}{2}) \hat{\mathbf{F}} = M_M F_i$, and the equilibria $\mathbf{m}^{eq}$ are given by:

$$ \mathbf{m}^{eq} = \rho(1, -2 + 3|\mathbf{v}|^2, 1 - 3|\mathbf{v}|^2, v_x, -v_x, v_y, -v_y, v_x^2 - v_y^2, v_x v_y)^T. \quad (15) $$

In the diagonal matrix, the parameters are selected as $\tau_{\rho}^{-1} = \tau_{j}^{-1} = 1.0$, $\tau_{e}^{-1} = \tau_{\zeta}^{-1} = \tau_{q}^{-1} = 1.1$ to achieve good stability [26]. The relaxation time $\tau_v$ is determined using the fluid kinematic viscosity with $\nu = c_s^2(\tau_v - 0.5)\Delta t$, where $c_s = c_l/\sqrt{3}$ is the speed of sound, and $c_l = 1$ is the lattice speed. The dynamic viscosity is represented as $\eta = \rho\nu$. The propagation process of the MRT LBM is written as:

$$ f_i(\mathbf{x} + \boldsymbol{\xi}_i \Delta t, t + \Delta t) = f_i^c(\mathbf{x}, t), \quad (16) $$

with the post-collision distribution $\mathbf{f}^c = M_M^{-1} \mathbf{m}^c$. The macroscopic forcing terms $\hat{\mathbf{F}}$ from [53] are used

$$ \begin{aligned} \hat{F}_0 &= 0, & \hat{F}_1 &= 6(v_x F_x + v_y F_y), & \hat{F}_2 &= -\hat{F}_1, & \hat{F}_3 &= F_x, & \hat{F}_4 &= -F_x, & \hat{F}_5 &= F_y, & \hat{F}_6 &= -F_y, & \hat{F}_7 &= 2(v_x F_x - v_y F_y), & \hat{F}_8 &= v_x F_y + v_y F_x. \end{aligned} \quad (17) $$

In this work, the force $\mathbf{F}$ is fluid-fluid interaction $\mathbf{F}_f$ given by the local mean gradient approximation [24]

$$ \mathbf{F}_f = (1 - 2B)\mathbf{F}^{\text{local}} + B\mathbf{F}^{\text{mean}}, \quad (18) $$

where $\mathbf{F}^{\text{local}}$ is the local force term

$$ \mathbf{F}^{\text{local}} = \psi(\mathbf{x}) G_1 \sum_{i=0}^{q} W(|\boldsymbol{\xi}_i|^2) \psi(\mathbf{x} + \boldsymbol{\xi}_i \Delta t) \boldsymbol{\xi}_i, \quad (19) $$

and $\mathbf{F}^{\text{mean}}$ is the mean force term

$$ \mathbf{F}^{\text{mean}} = G_1 \sum_{i=0}^{q} W(|\boldsymbol{\xi}_i|^2) \psi^2(\mathbf{x} + \boldsymbol{\xi}_i \Delta t) \boldsymbol{\xi}_i. \quad (20) $$

The coefficient $B$ is a tuning parameter to achieve thermodynamic consistency. $G_1 = -1$ is the interaction strength and $W(|\boldsymbol{\xi}_i|^2)$ are the force weights. The interaction potential $\psi$ is given by
incorporating the non-ideal EoS $p_{EoS}$ as $\psi = \sqrt{2(p_{EoS} - \rho c_s^2)/G_1c^2}$ [29]. In this paper, we use the CARNAHAN-STARLING EoS:

$$p_{EoS} = \rho RT \frac{1 + b\rho/4 + (b\rho/4)^2 - (b\rho/4)^3}{(1 - b\rho/4)^3} - a\rho^2 \quad (21)$$

where $a = 0.4963R^2T_c^2/p_c$, $b = 0.18727RT_c/p_c$, $R$ is the gas constant, $T_c$ and $p_c$ are critical temperature and pressure, respectively. The parameters $a, b$ and $R$ are used to determine $T_c$ and $p_c$, and they may affect the numerical performance. In their work, HUANG et al. [37] have studied the influences of these parameters on interface thickness, density ratio and numerical stability. In the current simulations, $a = 1, b = 4, R = 1$ are used following [29, 37], to realize an acceptable interface thickness of around 5 lattices and small spurious current of $\sim 10^{-3}$, while ensuring numerical stability.

By CHAPMAN-ENSKOG expansion, the NAVIER-STOKES equations can be derived from the above lattice BOLTZMANN modeling. However, this model is not able to independently tune the surface tension (i.e., without varying the temperature ratio). In the following, we first introduce several approaches to independently tune the surface tension in literature, and then introduce our own method to realize it. We note that, the scheme that we propose here to tune the surface tension is independent of the collision operator, i.e., it can be applied to LBGK/KBC+EDM and MRT+GUO frameworks.

## B. Pseudo-potential phase change with surface tension models

LI & LUO [38] suggest, in a MRT framework, to introduce an additional source term $C_{LL}$ which reads in the moment space:

$$C_{LL} = \begin{bmatrix} 0 \\ 1.5\tau_e^{-1}(K_{xx,LL} + K_{yy,LL}) \\ -1.5\tau_\zeta^{-1}(K_{xx,LL} + K_{yy,LL}) \\ 0 \\ 0 \\ 0 \\ -\tau_v^{-1}(K_{xx,LL} - K_{yy,LL}) \\ -\tau_v^{-1}K_{xy,LL} \end{bmatrix}, \quad (22)$$
in which, the tensor $K_{LL}$ is given by a thermodynamic expression and reads:

$$K_{\alpha\beta,LL} = \kappa \frac{G_1}{2} \psi(\mathbf{x}) \sum_{i=1} W(|\xi_i|^2) (\psi(\mathbf{x} + \xi_i \Delta t) - \psi(\mathbf{x})) \xi_{i\alpha} \xi_{i\beta} \quad (23)$$

where $\kappa$ is a parameter to control the surface tension. LYNCETT-BROWN et al. suggest a modification of the local source tensor [27], as follows:

$$K_{\alpha\beta,\text{LyB}} = \frac{G_1 c_s^2}{2} \psi(\mathbf{x}) \left[ -\kappa \sum_{i=1} W(|\xi_i|^2) (\psi(\mathbf{x} + \xi_i \Delta t) - \psi(\mathbf{x})) \xi_{i\alpha} \xi_{i\beta} \right. \\ \left. + \frac{\kappa+1}{6} \sum_{i=1} W(|\xi_i|^2) (\psi(\mathbf{x} + \xi_i \Delta t) - 2\psi(\mathbf{x}) + \psi(\mathbf{x} - \xi_i \Delta t)) \delta_{\alpha\beta} \right]. \quad (24)$$

In both cases, the control of surface tension is obtained in an MRT framework, and the modification of $\kappa$ is affecting the density ratio non-linearly in eq. (24), which could appear as a limitation since it leads to the deviation from Maxwell construction, i.e., the thermodynamic consistency is violated.

### III. SOURCE TERM OF SURFACE TENSION TUNING

#### A. A new generic source term

Inspired by the EDM concept of generating a source term from the difference between a modified moment equilibrium and an original moment equilibrium distribution function, and combined with the fact that extra source terms from [27, 38] leads to higher order moment modification, we decide that building a new source term based on equilibrium and enabling the control of higher moments can be a benefit. This operation requires a “equilibrium” distribution that has free moments, and previous works [54] have introduced such concept up to the N-th order and gave an example at third order. Such distribution can be written in its simple form as:

$$\tilde{f}(\xi, \rho, \mu, \Lambda, \alpha) = \frac{\rho \left( 1 + \operatorname{erf} \left( \frac{\alpha^T \lambda^{-1}}{\sqrt{2}} (\xi - \mu) \right) \right)}{(2\pi)^{D/2} (\det \Lambda)^{1/2}} \exp \left( -\frac{1}{2} (\xi - \mu)^T \Lambda^{-1} (\xi - \mu) \right), \quad (25)$$

where $\Lambda$ is a definite positive matrix, $\lambda$ is the diagonal standard deviation matrix of $\Lambda$ given by $\lambda_{ij} = \delta_{ij} \sqrt{\Lambda_{ij}}$ ($\delta_{ij}$ is the KRONECKER symbol), and $\alpha$ is a vector that drives the skewness of the distribution. However, in LBM implementation, only the projection of this distribution over the
HERMITE'S polynomials is at stake. We define the following macroscopic quantities:

$$ \int \tilde{f} \, d\xi = \rho, \qquad (26) $$

$$ \int \xi \tilde{f} \, d\xi = \rho v, \qquad (27) $$

$$ \int (\xi - v) \otimes (\xi - v) \tilde{f} \, d\xi = \rho\sigma. \qquad (28) $$

The generalized out-equilibrium distribution function from eq. (25) projected on the HERMITE'S polynomials and truncated at the second order becomes:

$$ \tilde{f}_i(\mathbf{x}, t, \rho, \mathbf{v}, \sigma) = w_i \rho \left( 1 + \frac{\mathbf{v} \cdot \boldsymbol{\xi}_i}{c_s^2} + \frac{(\mathbf{v} \cdot \boldsymbol{\xi}_i)^2}{2c_s^4} + \frac{\sigma : \boldsymbol{\xi}_i \otimes \boldsymbol{\xi}_i}{2c_s^4} - \frac{\boldsymbol{\xi}_i^2}{2c_s^2} - \frac{\mathbf{v}^2}{2c_s^2} - \frac{\mathrm{Tr}(\sigma)}{2c_s^2} + \frac{D}{2} \right), \quad (29) $$

where $c_s$ is commonly chosen as $c_s = \frac{1}{\sqrt{3}}$ in the lattice units.

Therefore, to generalize the EDM, we introduce a new source term that will allow tuning second order moment (e.g. the surface tension):

$$ S_i = F_i + \tilde{C}_i = [f_i^{eq}(\mathbf{x}, t, \rho, \mathbf{v} + g\Delta t) - f_i^{eq}(\mathbf{x}, t, \rho, \mathbf{v})] + [\tilde{f}_i(\mathbf{x}, t, \rho, \mathbf{v}, p + \tilde{K}) - f_i^{eq}(\mathbf{x}, t, \rho, \mathbf{v})], \quad (30) $$

where $F_i$ is a generic source term relative to the body force (e.g. $F_i = f_i^{eq}(\rho, \mathbf{v} + g\Delta t) - f_i^{eq}(\rho, \mathbf{v})$ when using EDM scheme) and $\tilde{K}$ is a second order tensor that can be adjusted to tune the surface tension. Note that $f_i^{eq}(\mathbf{x}, t, \rho, \mathbf{v}) = \tilde{f}(\mathbf{x}, t, \rho, \mathbf{v}, p)$ when $p = \rho c_s^2 I$ is the ideal pressure of the fluid. Using the previous definitions and expressions, the source term to tune the surface tension is given by

$$ \tilde{C}_i = w_i \rho \left( \frac{\tilde{K} : \boldsymbol{\xi}_i \otimes \boldsymbol{\xi}_i}{2c_s^4} - \frac{\mathrm{Tr}(\tilde{K})}{2c_s^2} \right). \qquad (31) $$

The expression of the tensor $\tilde{K}$ remains to be determined and it will be objective of the next section.

Building such source term from simple BHATNAGAR, GROSS and KROOK (BGK) expressions has the advantage to enable its adaptation to many posterior frameworks. Thus, the latter equation can be, after adjustments, included in the different collision operator and force scheme. As an example and a link with previous works, in MRT framework, the use of $M_M \cdot \tilde{C}_i$ as additional source term can be exactly equal to $C_{LL}$ if the tensor $\tilde{K}$ is taken such that $\tilde{K} = -K + \frac{3}{4} \mathrm{Tr}(K)I$, where K can be $K_{LyB}$ or $K_{LL}$.

## B. CHAPMAN-ENSKOG expansion

In order to determine the expression of the tensor $\tilde{K}$ to use in our source term, we performed a third order CHAPMAN-ENSKOG expansion for a generic source term. Both the equilibrium distribution in $\Omega_i$ and the force term $F_i$ should be expanded to the third order, since they both appear
at the right-hand side of eq. (9) in the collision step. This allows us to better fit with eq. (8) and adapt the expression of $\tilde{K}$ to the force scheme used. This genericity also offers the possibility to examine the effects of previous surface tension source terms. Further detail about the algebra used is provided in section A.

The TAYLOR-MAC LAURIN expansion of the Lattice BOLTZMANN Equation (LBE), followed by a small perturbation parameter, $\epsilon$, leads to the following equations up to the third order of perturbation:

$$ \epsilon^0: \quad f_i^{eq} = \hat{f}_i^{(0)} \qquad (32) $$

$$ \epsilon^1: \quad D_t^{(1)} (\hat{f}_i^{(0)}) - \Delta t S_i^{(1)} = -\hat{\omega} \Delta t \hat{f}_i^{(1)}, \qquad (33) $$

$$ \epsilon^2: \quad \tau_2 D_t^{(1)^2} \hat{f}_i^{(0)} + \partial_t^{(2)} \hat{f}_i^{(0)} + \frac{1}{\hat{\omega}} D_t^{(1)} S_i^{(1)} - \Delta t S_i^{(2)} = -\hat{\omega} \Delta t \hat{f}_i^{(2)}, \qquad (34) $$

$$ \epsilon^3: \quad \tau_3 D_t^{(1)^3} \hat{f}_i^{(0)} + \tau_2 \left[ \partial_t^{(2)} D_t^{(1)} + D_t^{(1)} \partial_t^{(2)} \right] \hat{f}_i^{(0)} + \partial_t^{(3)} \hat{f}_i^{(0)} \\ + \frac{1}{\hat{\omega}} \left\{ \left[ \partial_t^{(2)} + \tau_2 D_t^{(1)^2} \right] S_i^{(1)} + D_t^{(1)} S_i^{(2)} \right\} = -\hat{\omega} \Delta t \hat{f}_i^{(3)}, \qquad (35) $$

in which commonly, the superscript $(j)$ indicates the decomposition per order of perturbation, $\hat{\cdot}$ is the change of variable in lattice BOLTZMANN equation, $D_t^{(1)}$ is the material derivative. While the variables $\tau_2$ and $\tau_3$ are defined such as $\tau_2 = \frac{1}{2} - \frac{1}{\Delta t \hat{\omega}}$ and $\tau_3 = \frac{1}{6} + \frac{1}{\Delta t^2 \hat{\omega}^2} - \frac{1}{\Delta t \hat{\omega}}$.

Projecting the previous equations on conserved moments, one can obtain, for any forcing scheme, the zero-th moment macroscopical equation:

$$ \begin{aligned} \partial_t (\rho) + \partial_x \cdot (\rho \mathbf{v}^{eq}) - \frac{1}{\tau_2} \partial_x \cdot \partial_x \cdot (\Pi^{(1)}) & \\ - \frac{1}{12} A D_t^{*2} (M_0 [S_i^{(1)}]) + \left(1 - \frac{1}{\hat{\omega}}\right) D_t^* (M_0 [S_i^{(2)}]) = A M_0 [S_i^{(1)} + S_i^{(2)}], & \end{aligned} \qquad (36) $$

and first moment macroscopical equation:

$$ \begin{aligned} \partial_t (\rho \mathbf{v}^{eq}) + \partial_x \cdot (\rho \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} + p + \Pi^{(1)}) - \frac{1}{\tau_2} \partial_t \partial_x \cdot (\Pi^{(1)}) - \frac{1}{\tau_2} \partial_x \cdot D_t^* (M_2 [f_i^{(1)}]) \\ + \partial_x \cdot (M_2 [f_i^{(2)}]) - \frac{1}{12} A D_t^{*2} (M_1 [S_i^{(1)}]) + \left(1 - \frac{1}{\hat{\omega}}\right) D_t^* (M_1 [S_i^{(2)}]) = A M_1 [S_i^{(1)} + S_i^{(2)}], & \end{aligned} \qquad (37) $$

where $M_n(\cdot)$ is the n-th moment of a quantity, $D_t^*(M_n(\cdot)) = \partial_t^{(1)}(M_n(\cdot)) + \partial_x^{(1)}(M_{n+1}(\cdot)), \Pi^{(1)}$ is the second order moment of the $f_i^{(1)}$, and $A = \frac{2\Delta t}{2-\hat{\omega}\Delta t}$.
It is worth noting that previous equations can serve to analyze the consistency and precision of various sourcing scheme at third order. In the further considerations, the model choice has been made that forcing term is affecting the first perturbation order and the surface tension source term is related to second order perturbation. In the following, we focus specifically on EDM and GUO's forcing schemes, both including the local mean gradient approximation introduced in [24, 51]. Inspired from [28], the previous equations are simplified adopting steady-state hypotheses. Then, the expression of $\tilde{K}$ is obtained by equalizing the simplified equations with eq. (8). Within these conditions, the expression of source tensor $\tilde{K}$ has been exhibited for EDM and GUO's forcing schemes thus, using $P_n(\cdot)$ the n-cyclic permutation and the discretization constants $\lambda_1 = 1$, $\lambda_2 = \frac{1}{6}$, we obtained the two following expressions:

EDM:

$$ \tilde{K}_{\alpha\beta,\text{EDM}} = \left( \left[ \frac{G_1\lambda_1}{2} - \frac{G}{2} \right] \psi^2 + \left[ \psi G_1\lambda_2 - \kappa\rho \frac{d\psi}{d\rho} \right] \partial_{x_\gamma}^2(\psi) + \left[ (1+2B)\frac{G_1\lambda_2}{2} - \frac{\kappa}{2} \right] \left|\partial_{x_\gamma}(\psi)\right|^2 I_{\alpha\beta} \right. \\ \left. + \left[ \kappa - (1+2B)G_1^2\lambda_1^2 \frac{\Delta t^2}{\hat{\omega}} \frac{\psi^2}{\rho} - (1-2B)G_1\lambda_2 \right] \partial_{x_\alpha}(\psi) \otimes \partial_{x_\beta}(\psi) \right. \\ \left. + \frac{A}{12} \partial_{x_\gamma} \cdot \left\{ \frac{G_1\lambda_1\Delta t}{\rho} P_3(\psi \nabla_{x_\alpha}(\psi) \otimes p_{\beta\gamma}) \right\} - \frac{\tau_3}{\tau_2} \partial_{x_\gamma} \cdot \left( M_{3,\alpha\beta\gamma} [f_i^{(1)}] \right) \right) \quad (38) $$

Guo:

$$ \tilde{K}_{\alpha\beta,\text{Guo}} = \left( \left[ \frac{G_1\lambda_1}{2} - \frac{G}{2} \right] \psi^2 + \left[ G_1\lambda_2\psi + \frac{G_1\lambda_1Bc_s^2}{6\Delta t}\psi - \kappa\rho\frac{d\psi}{d\rho} \right] \partial_{x_\gamma}^2(\psi) + \left[ (1+2B)\frac{G_1\lambda_2}{2} + \frac{G_1\lambda_1Bc_s^2}{6\Delta t} - \frac{\kappa}{2} \right] \partial_{x_\gamma}(\psi) \right. \\ \left. + \left[ \kappa - (1-2B)G_1\lambda_2 + \frac{G_1\lambda_1Bc_s^2}{3\Delta t} \right] \partial_{x_\alpha}(\psi) \otimes \partial_{x_\beta}(\psi) + \left[ \frac{4G_1\lambda_1Bc_s^2}{12\Delta t} \right] \psi \partial_{x_\alpha}(\partial_{x_\beta}(\psi)) \right. \\ \left. - \frac{\tau_3}{\tau_2} \partial_{x_\gamma} \cdot \left( M_{3,\alpha\beta\gamma} [f_i^{(1)}] \right) \right) \quad (39) $$

These expressions of source tensors are both, up to the third perturbation order, preserving thermodynamic consistency from LANDAU free energy theory by using a hybrid density-pseudo-potential form while allowing a surface tension tuning.

### C. Independent of the collision operator implementation

In section II, we have introduced two frameworks of LBM. Here we introduce how the surface tension term is incorporated in the two frameworks. For the LBGK/KBC+EDM framework, the implementation is quite straightforward. We can directly add eq. (31) to the right-hand side of eq. (13). We note that, for the implementation of pressure tensor $\tilde{K}_{\alpha\beta}$, we neglect the last
term representing the high order terms, in order to achieve better numerical stability for low fluid viscosity. With $G_1\lambda_1 = G$, $\lambda_1 = 1$, $\lambda_2 = \frac{1}{6}$, $\tilde{K}_{\alpha\beta}$ is simplified as

$$ \tilde{K}_{\alpha\beta,\text{EDM}} = \left( 2B \left[ \frac{1}{6} - \kappa \rho \frac{d\psi}{d\rho} \right] \partial_{x_\gamma}^2 (\psi) + \left[ (1+2B) \frac{G_1}{12} - \kappa \rho \frac{d\psi}{d\rho} \right] \left| \partial_{x_\gamma} (\psi) \right|^2 I_{\alpha\beta} \right. \\ \left. + \left[ \kappa - (1+2B)G_1^2 \frac{\Delta t^2}{\hat{\omega}} \frac{\psi^2}{\rho} - \frac{1}{6}(1-2B)G_1 \right] \partial_{x_\alpha} (\psi) \otimes \partial_{x_\beta} (\psi) \right. , \quad (40) \\ \left. + \frac{A}{12} \partial_{x_\gamma} \cdot \left\{ \frac{G_1 \Delta t}{\rho} P_3 (\psi \nabla_{x_\alpha} (\psi) \otimes p_{\beta\gamma}) \right\} \right) $$

where $\hat{\omega} = \frac{2}{1+2\nu c_s^2}$ and $A = \frac{2}{2-\hat{\omega}}$.

For the MRT+GUO framework, an additional term $M_M \tilde{C}_i$ is added to the right-hand side of eq. (14), which is written as

$$ M_M \cdot \tilde{C}_i = \begin{bmatrix} 0 \\ 3\rho(\tilde{K}_{xx} + \tilde{K}_{yy}) \\ -3\rho(\tilde{K}_{xx} + \tilde{K}_{yy}) \\ 0 \\ 0 \\ 0 \\ \rho(\tilde{K}_{xx} - \tilde{K}_{yy}) \\ \rho \tilde{K}_{xy} \end{bmatrix} . \quad (41) $$

In eq. (41), the tensor $\tilde{K}_{\alpha\beta}$ is also implemented without the last term representing the higher order term, i.e.,

$$ \tilde{K}_{\alpha\beta,\text{Guo}} = \left( \left[ \frac{G_1}{6}\psi + \frac{G_1 B}{18\Delta t}\psi - \kappa\rho\frac{d\psi}{d\rho} \right] \partial_{x_\gamma}^2 (\psi) + \left[ (1+2B)\frac{G_1}{12} + \frac{G_1 B}{18\Delta t} - \frac{\kappa}{2} \right] \left| \partial_{x_\gamma} (\psi) \right|^2 I_{\alpha\beta} \right. \\ \left. + \left[ \kappa - \frac{1}{6}(1-2B)G_1 + \frac{G_1 B}{9\Delta t} \right] \partial_{x_\alpha} (\psi) \otimes \partial_{x_\beta} (\psi) + \frac{G_1 B}{9\Delta t} \psi \partial_{x_\alpha} \left( \partial_{x_\beta} (\psi) \right) \right) \quad (42) $$

## IV. NUMERICAL EXPERIMENTS

As we explained, the source term that we proposed to tune the surface tension is applicable to various collision operator. We have tested both the LBGK/KBC+EDM and MRT+GUO frameworks, and found the source term works better for the latter framework, in terms of achieving large ranges of fluid viscosity and surface tension. Therefore, in the numerical simulations, we use the
<figure><img src="image_2.png" /><figcaption>Figure 1. Comparison of liquid-vapor density ratios at various temperature ratios between current LBM simulations and Maxwell equal area rule.</figcaption></figure>

MRT+GUO framework. The results of LBGK/KBC+EDM framework are given in section B for interested readers.

In the following, we first try to achieve the liquid-vapor density ratio curve under different temperature ratios. In eq. (42) of section III, the parameter *B* can be adjusted to reach thermo-dynamic consistency. By selecting the value of *B* = -0.35, good agreement between our simulation results and MAXWELL equal area rule is achieved in the large range of temperature ratio of *T/T*<sub>c</sub> ∈ (0.55, 0.95), as shown in fig. 1.

Then we conduct 3 series of numerical experiments to evaluate the capability of the developed current model. The first series of simulations is to model the stationary liquid-vapor coexistence of a droplet, and check the magnitude of spurious current, as well as a tunable range of surface tension. The second and third series of numerical experiments are to check the model performance of modeling slight and intense dynamic behaviors. The 3 subsections are given as follows in details.

A. Suspended droplet

In this subsection, we model the stationary liquid droplet suspended in a vapor environment. We use a 2D simulation domain of 240 × 240 *lattices*² with droplet initial radius of 48 lattices, while periodic boundary is applied in both *x* and *y* directions. Liquid and vapor kinematic viscosity are set as $ν_l = 0.01$ and $ν_v = 1/6$, respectively. The temperature ratio is selected at $T/T_c = 0.75$, with
<figure><img src="image_3.png" /><figcaption>Figure 2. Modeling of liquid droplet suspended in vapor. (a) Density contour with zoom-in of the interfacial area. (b) Distribution of spurious current.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. The variation of (a) surface tension ratio, (b) maximum spurious current with parameter $κ$.</figcaption></figure>

the obtained liquid-vapor density ratio of $ρ_l/ρ_v = 22.3$. The reason we select this density ratio is to maintain a relatively low level of spurious current, in order not to affect the simulation accuracy. In eq. (42), the coefficient $κ$ is used to adjust the surface tension. We first apply $κ = 0$ in the current simulation, and the obtained result is shown in fig. 2. fig. 2a shows the density contour with a diffusive interface, while fig. 2b shows the spurious current with the maximum value less than $4.0 \times 10^{-3}$. A dashed black circle is plotted upon the droplet profile in fig. 2a. No deviation of the droplet profile from the circular shape can be directly observed. Moreover, a zoom-in of a quarter of the droplet is inserted in fig. 2a. Still, no variation of the density landscape is observed. The result shows that the spurious current is small enough to ensure accurate modeling of two-phase dynamics. In the pseudo-potential model, the surface tension is normally measured by $σ = (p_{in} - p_{out}) \cdot r$, where $p_{in}$ ($p_{out}$ respectively) and $r$ are the pressure inside (outside respectively) the droplet and
the droplet radius. The obtained surface tension is $\sigma = 0.0091$ at $\kappa = 0$.

Using the formal definition of the surface tension with the modified hybrid KORTEWEG's stress tensor [24], one can obtain:

$$ \sigma = \int_{-\infty}^{+\infty} (p_{th,xx} - p_{th,yy}) dx = \int_{-\infty}^{+\infty} \kappa \left[ \partial_x (\psi)^2 - \partial_y (\psi)^2 \right] dx. \quad (43) $$

For the 1D case, it is simplified as $\sigma = \int_{-\infty}^{+\infty} \kappa \left(\frac{d\psi}{dx}\right)^2 dx$. In the case of using traditional density-based KORTEWEG's stress tensor [45], the pseudo-potential $\psi$ is replaced by density $\rho$, i.e., $\sigma = \int_{-\infty}^{+\infty} \kappa \left(\frac{d\rho}{dx}\right)^2 dx$. With both forms of stress tensors, a linear relation between the surface tension $\sigma$ and the capillary coefficient $\kappa$ is represented. We further vary the values of $\kappa$ to obtain the tunable range of the surface tension. As shown in fig. 3a, when $\kappa$ ranges from -3.2 to 1.0, the surface tension ratio $\sigma/\sigma_{\kappa=0}$ varies from 0.01 to 1.4, i.e., a remarkable tunable range of surface tension of 140 times is achieved. Moreover, an approximately linear relation of $\sigma/\sigma_{\kappa=0} = 1.022 + 0.327\kappa$ is obtained between $\kappa$ and $\sigma/\sigma_{\kappa=0}$. In the work of a widely applied pseudo-potential LBM [26], the tunable range is around 76 times from $\sigma = 0.0215$ to $\sigma = 1.64$, with $\kappa$ decreasing from 0.99 to 0. For the tuning approach by multi-range pseudo-potential [22], LI & LUO [26] checked the performance and found the surface tension variation is only 2.4 times, with the tuning parameter varying 100 times. In the work of WANG et al. [55], they also proposed a method to tune the surface tension. Their plot shows a surface tension variation range of 20 times, with the tuning parameter varying from 0 to 2.4. In terms of other LB method, in the work of MAZLOOMI et al. [56], they applied a free energy model to model multiphase flows. In their plot, with the tuning parameter increasing from $\kappa = 0.00468$ to $\kappa = 0.027$ (about 6 times), the surface tension increases from $\sigma = 0.353$ to $\sigma = 0.501$ with the variation of only 1.4 times. Compared with these tuning approaches, our proposed method is able to adjust the surface tension with a much wider range. In terms of spurious current, as shown in fig. 3b, the maximum values corresponding to $\kappa \in (-3.2, 1.0)$ are in the range of $(3.1 \times 10^{-3}, 1.0 \times 10^{-2})$, which are normally small enough to be neglected. We note that, the maximum spurious current does not vary monotonically with the tuning parameter $\kappa$ (the surface tension). The reason is that, the variation of parameter $\kappa$ leads to the change of complex non-linear terms such as $\kappa\rho\frac{d\psi}{d\rho}\partial_{x\gamma}^2(\psi)$, $\kappa\partial_{x\alpha}(\psi)\otimes\partial_{x\beta}(\psi)$ in the pressure tensor, which show non-monotonic influence on the fluid-fluid interactions, thus leading to non-monotonic behaviors of the maximum spurious current. On the other hand, if $\kappa$ is fixed and the surface tension is tuned by varying the temperature ratio $(T/T_c)$, the maximum spurious current shows a monotonic negative correlation, being consistent with the results in literature [18, 57, 58]. A small drawback of the current model is the liquid vapor density ratio variation to some extent,
<figure><img src="image_5.png" /><figcaption>Figure 4. Modeling results of droplet oscillation. (a) The oscillation process from $t_0$ to $t_4$ of Case1. (b) The comparison of normalized droplet radius along $x$ direction between 3 cases.</figcaption></figure>

i.e., decreasing from $\rho_l/\rho_v = 32.8$ at $\kappa = -3.2$ to $\rho_l/\rho_v = 18.3$ at $\kappa = 1.0$ with the variation range of around 45%. The main reason of the density ratio variation is also because $\kappa$ affects the nonlinear terms of $\kappa \rho \frac{d\psi}{d\rho} \partial_{x_\gamma}^2(\psi)$, $\kappa \partial_{x_\alpha}(\psi) \otimes \partial_{x_\beta}(\psi)$, etc., which changes the fluid-fluid interaction strength to some extent. We note that, for the surface tension tunable range of 140 times, the density ratio variation range of 45% is acceptable in modeling multiphase flows. From the above results we can observe that, our model can realize high surface tension range with negligible spurious current and acceptable density ratio variation.

## B. Oscillating droplet

In this subsection, we investigate the droplet oscillation under different liquid properties. The simulation setups are the same as those in the previous subsection, except the initial droplet is elliptical with radii in $x$ and $y$ directions of $R_{x,0} = 53$ and $R_{y,0} = 47$ lattices, respectively. The initial profile of the droplet is illustrated in fig. 4a at $t_0$. The liquid-vapor density ratio is $\rho_l/\rho_v = 22.3$, the liquid kinematic viscosity and surface tension are $\nu_l = 0.01$ and $\sigma = 0.0091$, respectively. Since the droplet is non-circular initially (at $t_0$ in fig. 4a), the surface tension makes it oscillate until a final circular shape is obtained.

The oscillation process of one period from $t_0$ to $t_4$ is also given in fig. 4a. We can observe that, the radius in $x$ direction $R_x$ first shrinks ($t_1$) to an quasi-equilibrium state, and then reaches its minimum ($t_2$). Afterwards, it expands gradually from another quasi-equilibrium state $t_3$ to maximum $R_x$ at $t_4$, which is the starting shape for the next period. The variation of $R_y$ is the
opposite to $R_x$. We recorded the evolution of $R_x$, normalized by the final radius of the circular shape at equilibrium $R_0 = \sqrt{R_{x,0}R_{y,0}}$, as shown in fig. 4b. The green dashed curve notated “Case1” corresponds to the results in fig. 4a, while the curves notated “Case2” and “Case3” correspond to the results with varied kinematic viscosity of $\nu_{l,2} = 0.167$ and varied surface tension of $\sigma_3 = 0.0059$, respectively. Comparing Case1 and Case2, we find the droplet oscillation periods are the same, while the oscillation radius evolves faster to equilibrium in Case2. Reason for this behavior is explained by the equal liquid surface tensions, while the viscosity is higher in Case2, leading to higher viscous dissipation. By comparing Case1 and Case3, we find the droplet oscillation period is longer in Case3, while the normalized radius evolution is similar. Reason for this observation is linked the same liquid viscosity, while the surface tension is smaller in Case3, leading to less confinement of liquid-vapor interface. Moreover, the oscillation period of a 2D elliptical droplet can be given as follows [26, 59]:

$$T_o = 2\pi \left[ n(n^2 - 1) \frac{\sigma}{\rho_l R_0^3} \right]^{-1/2} \quad (44)$$

where $n$ denotes the model of oscillation, which is given by $n = 2$ for an initial elliptic shape. From eq. (44), we obtain the calculated oscillation periods for the three cases as $T_{o,1} = T_{o,2} = 5445$, $T_{o,3} = 6871$. In current simulations, the modeled oscillation periods are $T_{o,1} = T_{o,2} = 5400$, $T_{o,3} = 6800$. The good agreement between theoretical and modeling results (relative error less than 1%) shows our model is accurate in capturing the simple droplet dynamics such as oscillations. A small drawback of the simulation result lies in the slightly higher steady state droplet radius in Case3, compared with Case1 and Case2. The reason is attributed to the density ratio variation induced by surface tension change. Compared to Case1 and Case2 with a surface tension of $\sigma_1 = 0.0091$, the surface tension in Case3 with $\sigma_3 = 0.0059$ is lower. The corresponding density ratios in Case1/Case2 and Case3 are $(\rho_l/\rho_v)_{Case1,2} = 0.334/0.015 \approx 22.3$ and $(\rho_l/\rho_v)_{Case3} = 0.334/0.0128 \approx 26.1$, respectively. We can see the vapor density decreases in Case3, while the liquid density remains almost constant. Since the total mass of liquid and vapor is conserved, at a lower vapor density in Case3, more mass is present in the liquid form. With a constant liquid density, the liquid volume increases, leading to a slightly higher (less than 2%) droplet radius.

### C. Splashing droplet

In this subsection, we apply the model to investigate the droplet impact on a thin liquid film and resulting droplet splashing. To capture the dynamics of thin liquid sheet, we follow the setups
<figure><img src="image_6.png" /><figcaption>Figure 5. Droplet impact and splashing dynamics on thin liquid film at different conditions. (a) Simulation setup. (b) Case1: $Re = 100$, $We = 220$. (c) Case2: $Re = 2000$, $We = 220$. (d) Case3: $Re = 2000$, $We = 4500$. $t^*$ denotes dimensionless time since the droplet touches the liquid film. (e) The zoom-in of the liquid sheet and secondary droplets with mesh.</figcaption></figure>

of [36], as shown in fig. 5a. The domain size is 2400 × 1200 lattices², with the left and right sides set as periodic boundaries and the top and bottom as no-slip walls. The droplet diameter is $D_0 = 200$
<figure><img src="image_7.png" /><figcaption>Figure 6. Normalized droplet spreading radius against dimensionless time. (a) The illustration of droplet spreading radius. (b) The scaling relation. Symbols: simulations with various REYNOLDS Re and WEBER We numbers, solid line: square root of time law.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Droplet impact and splashing dynamics on thin liquid film at maximum reachable Reynolds and Weber numbers, respectively. (a) Maximum reachable REYNOLDS number of Re = 24000 at We = 103. (b) Maximum reachable WEBER number of We = 10500 at Re = 1000.</figcaption></figure>

lattices, the thin film thickness is $h_2 = 100$ lattices and the central height of the droplet from the film is $h_1 = 220$ lattices. The droplet initial falling velocity is given as $U = 0.12$ neglecting gravity. The liquid kinematic viscosity is set at a value of $\nu_l = 0.012$, while the vapor kinematic viscosity $\nu_v$ is set about 20 times of $\nu_l$ following [26]. We note that, the collision operator we apply in this work is MRT, which is more advantageous than the single-relaxation-time (SRT) collision operator in realizing a wider range of viscosity. Results in literature [60] showed that, a variation of
several hundred times in fluid viscosity can be realized when using MRT LBM. The liquid viscosity and surface tension are varied to realize a large range of REYNOLDS number $Re \in (100, 2000)$ and WEBER number $We \in (220, 4500)$. These dimensionless number are defined in such configuration by:

$$Re = UD_0/v_l, \quad We = \rho_l U^2 D_0 / \sigma. \qquad (45)$$

We have conducted a series of simulations with different REYNOLDS and WEBER numbers. Here we compare the results of three cases, i.e., Case1 with $Re = 100$, $We = 220$, Case2 with $Re = 2000$, $We = 220$ and $Re = 2000$, $We = 4500$, respectively. The snapshots of droplet dynamics are shown in fig. 5, where the dimensionless time is defined as $t^* = Ut/D_0$. We observe that, after the droplet impacts the film at $t^* = 0.6$, two thin liquid sheets are generated at the left and right sides of the droplet. As time goes on, the droplet continuously merges with the liquid film, the liquid sheets become thicker/higher and the droplet becomes thinner (even lower than the liquid film). As shown at $t^* = 3.0$ and afterwards, a “crown” shape is formed. In Case1 with low REYNOLDS and WEBER numbers, the inertial force is not high compared to viscous and capillary forces, and the liquid sheets are relatively thick. By increasing the REYNOLDS from $Re = 100$ in Case1 to $Re = 2000$ in Case2 while maintaining the same WEBER number, we observe the liquid sheets are thinner and higher in general, indicating the inertial force is much stronger compared to the viscous force. However, the liquid sheets still hold without breaking, indicating the inertial force is not strong enough compared to the capillary force. When we increase the WEBER number from $We = 220$ in Case2 to $We = 4500$ in Case3, we find the thin liquid sheets break up into multiple secondary droplets at the tip, i.e., the droplet splashes as shown in the 4th subfigure of fig. 5c ($t^* = 6.6$).

A not expected observation in fig. 5c is that, the secondary droplets vanish with time, which is not typically physical. Therefore, we have compared the mesh resolution with secondary droplets size. As shown in fig. 5e, the diameters of secondary droplets are much larger than the film thickness. Since the liquid sheets evolve physically instead of disappearing, we think the mesh resolution is not the reason. We also conducted the simulation using the same setups as fig. 5c but using the widely-applied improved model of LI et al. [26], and the secondary droplets are also found to disappear with time. These behaviors remind us of shrinkage effects observed in [61], in the sense that the droplets or bubbles below a critical size tend to shrink spontaneously. The shrinkage effect is an inherent feature of diffused interface numerical models, which are driven by the redistribution of the interface and bulk energies to minimize the system's total energy. It is noted the critical size increases with the system size and interface thickness for suspending droplets/bubbles [61, 62]
and depends on the contact angle for sessile droplets [63]. Such dependencies have been discussed extensively in the phase-field LB model [61–63] and recently also validated in the pseudo-potential model [64]. Accordingly, in our cases, the large system size leads to a relatively large critical size, so that the secondary droplets below this size are unstable and finally disappear. A possible solution of this vanishing issue might be to apply the sharp-interface volume of fluid (VoF) approach. For instance, in the work of KULKRA et al. [65], they are capable of capturing very small droplets during the atomizing pulsed jet. However, how the droplet would behave if their radii are smaller than the mesh size is still a problem required to be explored.

From fig. 5, we can observe that the radius of the “crown” continuously increases with time. This phenomenon is normally quantified by the droplet spreading radius $r$, as shown in fig. 6a. As theoretically explained in [66], the spreading radius $r$ follows a square root of time law

$$r/D_0 \approx C\sqrt{Ut/D_0}, \qquad (46)$$

where the coefficient $C \approx 1.1$ for 3D or axisymmetric modeling. In the current 2D simulations, we use $C = 1.06$. fig. 6b shows the comparison of simulation results with the empirical relation, where good agreement is achieved for all cases.

Finally, we would like to discuss about the limits of Re and We in our proposed model. For the lower-bound, similar to other models, the minimum values can be approximately zero since the droplet falling velocity $U$ can be close to zero. For the upper-bound, with the set of $U = 0.12$, the maximum Re and We can be realized are $Re = 24000$ at $We = 103$ and $We = 10500$ at $Re = 1000$, respectively. The results of droplet impact and splashing at these two conditions are shown in fig. 7. For the maximum reachable Re and We at the same time in one simulation, they are approximately $Re = 2000$ and $We = 4500$, as discussed in fig. 5d. To summarize this subsection, we can accurately model the droplet impact and splashing dynamics with large ranges of REYNOLDS and WEBER numbers. Especially for the WEBER number, benefiting from the proposed model for tuning surface tension, we can reach WEBER number as high as $We = 10500$. In LI et al. [26], the authors applied the standard MRT collision operator without the tuning method of surface tension. The Weber number in their simulations is $We = 103$. In QIN et al. [36], the authors applied an entropic-multiple-relaxation-time (EMRT) collision operator with the multi-range pseudo-potential to tune the surface tension. The Weber number in their simulations is $We = 220$. Compared to those studies, we can reach the WEBER number at a higher order of magnitude.
V. CONCLUSIONS

Pseudo-potential multiphase lattice BOLTZMANN models are widely applied in various scientific and engineering fields. However, a major drawback that restricts the applications of this class of model is the difficulty in tuning surface tension, as well as the limited tunable range. In this work, we have offered an alternative approach by proposing an additional adjustment term originated from a third-order CHAPMAN-ENSKOG analysis. This term allows us to tune the surface tension independently from the temperature ratio. More importantly, our simulation results show the proposed term is valid in various collision operators and force treatments, including the widely applied LBGK/KBC+EDM and MRT+GUO frameworks. A few benchmarks are given to validate the model. The first one is the achievement of various liquid-vapor density ratios spanning a wide range of temperature ratios. The second one lies in the linear adjustment of surface tension range over remarkably 140 times for a suspended droplet, while maintaining small spurious current with its maximum value less than 0.01. Droplet deformation by oscillation is also modeled under different viscosity and surface tension, with the obtained oscillation period agreeing well with the theoretical model. The proposed term is also validated as practical in intensive interfacial flows when a droplet impacts on a liquid film, at various WEBER numbers up to 10500. Detailed information of liquid sheet growth, breakup and splashing are all well-captured. The crown radii at different WEBER numbers also agree with the empirical equation, indicating the high accuracy of the proposed approach. With a solid mathematical analysis and detailed numerical validation, we hope the proposed adjustment approach of surface tension will benefit the LB community to broaden the application fields of pseudo-potential models.

ACKNOWLEDGMENTS

Swiss National Science Foundation (SNSF Grant No. 175793), National Natural Science Foun-
dation of China (grant No. 12302337, No. 12388101), ETH Research Grant (project no. ETH-08
21-2), the 111 Project of China through Grant No. B17037, are acknowledged for their financial
support. The Computing Center in Xi'an and the Swiss National Super Computing Center (project
Nos. s1081 and go09) are acknowledged for the computational support.
Appendix A: CHAPMAN-ENSKOG expansion details

The LBE is solved numerically by LBM algorithm and is given by:

$$ \hat{f}_i(\mathbf{x} + \boldsymbol{\xi}_i \Delta t, t + \Delta t) = \hat{f}_i(\mathbf{x}, t) - \Delta t \hat{\omega} (\hat{f}_i(\mathbf{x}, t) - f_i^{eq}(\mathbf{x}, t)) + \Delta t S_i(\mathbf{x}, t) + \mathcal{O}(\Delta t^3), \quad (A1) $$

where $\hat{f}_i$ is the distribution function (from the BOLTZMANN equation) after the change of variable $\hat{f}_i = f_i + \frac{\Delta t \hat{\omega}}{2-\hat{\omega}\Delta t} (f_i - f_i^{eq}) - \frac{\Delta t}{2-\hat{\omega}\Delta t} S_i$, and $S_i$ is a generic source term (which can be from EDM or GUO's forcing scheme plus LYNCETT-BROWN surface tensor corrector, for example). Using the TAYLOR-MAC LAURIN expansion, combine with the time scale decomposition $\Delta t \partial_t = \sum_{k=1}^{\infty} \epsilon^k \partial_t^{(k)}$, $\Delta t \nabla_x = \epsilon^1 \nabla_x^{(1)}$, $\hat{f}_i = \sum_{k=0}^{\infty} \epsilon^k \hat{f}_i^{(k)}$, $S_i = \epsilon^1 S_i^{(1)} + \epsilon^2 S_i^{(2)}$, one can get the following expansion:

$$ \sum_{m=1}^{\infty} \frac{1}{m!} \left[ \sum_{k=2}^{\infty} \epsilon^k \partial_t^{(k)} + \epsilon^1 (\partial_t^{(1)} + \boldsymbol{\xi}_i \cdot \nabla_x^{(1)}) \right]^m \sum_{r=0}^{\infty} \epsilon^r \hat{f}_i^{(r)} = -\Delta t \hat{\omega} \left( \sum_{r=0}^{\infty} \epsilon^r \hat{f}_i^{(r)} - f_i^{eq} \right) + \Delta t \epsilon^1 S_i^{(1)} + \Delta t \epsilon^2 S_i^{(2)}. \quad (A2) $$

One can also introduce the material (or particle) derivative of any quantity $a$ at first order of perturbation ($D_t^{(1)}(a) = \partial_t^{(1)}(a) + \boldsymbol{\xi}_i \cdot \nabla_x^{(1)}a$); then, the previous equation is gathered per perturbation order:

$$ \epsilon^0: \quad f_i^{eq} = \hat{f}_i^{(0)} \quad (A3) $$

$$ \epsilon^1: \quad D_t^{(1)} (\hat{f}_i^{(0)}) - \Delta t S_i^{(1)} = -\hat{\omega} \Delta t \hat{f}_i^{(1)}, \quad (A4) $$

$$ \epsilon^2: \quad \tau_2 D_t^{(1)^2} \hat{f}_i^{(0)} + \partial_t^{(2)} \hat{f}_i^{(0)} + \frac{1}{\hat{\omega}} D_t^{(1)} S_i^{(1)} - \Delta t S_i^{(2)} = -\hat{\omega} \Delta t \hat{f}_i^{(2)}, \quad (A5) $$

$$ \epsilon^3: \quad \tau_3 D_t^{(1)^3} \hat{f}_i^{(0)} + \tau_2 [\partial_t^{(2)} D_t^{(1)} + D_t^{(1)} \partial_t^{(2)}] \hat{f}_i^{(0)} + \partial_t^{(3)} \hat{f}_i^{(0)} \\ + \frac{1}{\hat{\omega}} \left\{ [\partial_t^{(2)} + \tau_2 D_t^{(1)^2}] S_i^{(1)} + D_t^{(1)} S_i^{(2)} \right\} = -\hat{\omega} \Delta t \hat{f}_i^{(3)}, \quad (A6) $$

with the following introduced variables $\tau_2 = \frac{1}{2} - \frac{1}{\Delta t \hat{\omega}}$ and $\tau_3 = \frac{1}{6} + \frac{1}{\Delta t^2 \hat{\omega}^2} - \frac{1}{\Delta t \hat{\omega}}$. Noting $M_n(\cdot)$ the n-th moment operator, the solvability conditions can be written:

$$ \widehat{M}_n^{(k)} = M_n(\hat{f}_i^{(k)}) = \sum_i \xi_i^n \hat{f}_i = \frac{2}{2 - \hat{\omega} \Delta t} \left[ M_n(f_i^{(k)}) - \frac{\Delta t}{2} M_n(S_i^{(k)}) \right] \quad \forall k \ge 1, \forall n \ge 0, \quad (A7a) $$

$$ \text{with } M_0(f_i^{(k)}) = 0 \text{ and } M_1(f_i^{(k)}) = 0 \quad \forall k \ge 1. \quad (A7b) $$

Then, extracting moments from eqs. (A3) to (A6), one can get moments of non-equilibrium distributions.
a. At $\epsilon^1$ scale: Taking the zeroth order moment and first moment of eq. (A4) and noting $A = \frac{2\Delta t}{2-\hat{\omega}\Delta t}$ one can get the following equation:

$$ \partial_t^{(1)}(\rho) + \partial_x^{(1)} \cdot (\rho \mathbf{v}^{eq}) = AM_0(S_i^{(1)}), \qquad (A8) $$

$$ \partial_t^{(1)} (\rho \mathbf{v}^{eq}) + \partial_x^{(1)} \cdot (\rho \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} + p) = \frac{2\Delta t}{2-\hat{\omega}\Delta t} \sum_i \xi_i S_i^{(1)} = AM_1(S_i^{(1)}). \qquad (A9) $$

The second moment is given by:

$$ \partial_t^{(1)} (\rho \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} + p) + \partial_x^{(1)} \cdot (\rho \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} + P_3 (\mathbf{v}^{eq} \otimes p)) = \frac{2\Delta t}{(2-\hat{\omega}\Delta t)} (M_2(S_i^{(1)}) - \hat{\omega}\Pi^{(1)}). \qquad (A10) $$

From the latter equation, the expression of $\Pi^{(1)}$ can be obtained:

$$ \Pi^{(1)} = - \frac{1}{A\hat{\omega}} \left\{ p \cdot \partial_x^{(1)} (\mathbf{v}^{eq}) + \left( p \cdot \partial_x^{(1)} (\mathbf{v}^{eq}) \right)^T + \partial_x^{(1)} \cdot (\mathbf{v}^{eq}) p - \partial_x^{(1)} \cdot (\mathbf{v}^{eq}) \rho \partial_\rho (p) \right\} \\ + \frac{1}{\hat{\omega}} \left[ M_2(S_i^{(1)}) - (M_1(S_i^{(1)}) \otimes \mathbf{v}^{eq} + \mathbf{v}^{eq} \otimes M_1(S_i^{(1)})) + M_0(S_i^{(1)}) (\mathbf{v}^{eq} \otimes \mathbf{v}^{eq} - \partial_\rho (p)) \right]. \qquad (A11) $$

Continuing to extract higher non-conserved moments, the third moment reads:

$$ \partial_t^{(1)} (\rho \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} + P_3 (\mathbf{v}^{eq} \otimes p)) + \partial_x^{(1)} \cdot (\rho \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} + P_6 (\mathbf{v}^{eq} \otimes \mathbf{v}^{eq} \otimes p) + p \otimes p) \\ = AM_3 (S_i^{(1)}) - A\hat{\omega}M_3 (f_i^{(1)}) \qquad (A12) $$

The latter is leading to the third moment of the first perturbation order distribution, given by:

$$ M_{3,\alpha\beta\gamma} (f_i^{(1)}) = [AM_{2,\alpha\beta} (S_i^{(1)}) - A\hat{\omega}\Pi_{\alpha\beta}^{(1)}] \mathbf{v}_{\gamma}^{eq} + \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot (p_{\alpha\beta} p_{\gamma\delta} + \mathbf{v}_{\alpha}^{eq} \mathbf{v}_{\beta}^{eq} p_{\gamma\delta} + \mathbf{v}_{\alpha}^{eq} \mathbf{v}_{\delta}^{eq} p_{\gamma\beta} + \mathbf{v}_{\beta}^{eq} \mathbf{v}_{\delta}^{eq} p_{\gamma\alpha}) \\ + [AM_{1,\alpha} (S_i^{(1)}) - \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot (\rho \mathbf{v}_{\alpha}^{eq} \mathbf{v}_{\delta}^{eq} + p_{\alpha\delta})] (\mathbf{v}_{\beta}^{eq} \mathbf{v}_{\gamma}^{eq} + \frac{p_{\beta\gamma}}{\rho}) + [\rho \mathbf{v}_{\alpha}^{eq} \mathbf{v}_{\beta}^{eq} \mathbf{v}_{\delta}^{eq} + P_3 (\mathbf{v}_{\alpha}^{eq} p_{\beta\delta})] \cdot \partial_{\mathbf{x}_{\delta}}^{(1)} (\mathbf{v}_{\gamma}^{eq}) \\ + [AM_{1,\beta} (S_i^{(1)}) - \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot (\rho \mathbf{v}_{\beta}^{eq} \mathbf{v}_{\delta}^{eq} + p_{\beta\delta})] \frac{p_{\alpha\gamma}}{\rho} + [AM_{1,\gamma} (S_i^{(1)}) - \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot (\rho \mathbf{v}_{\gamma}^{eq} \mathbf{v}_{\delta}^{eq} + p_{\gamma\delta})] \frac{p_{\alpha\beta}}{\rho} \\ + [AM_0 (S_i^{(1)}) - \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot (\rho \mathbf{v}_{\delta}^{eq})] [\mathbf{v}_{\alpha}^{eq} \left(\partial_{\rho} (p_{\beta\gamma}) - \frac{p_{\beta\gamma}}{\rho}\right) + \mathbf{v}_{\beta}^{eq} \left(\partial_{\rho} (p_{\alpha\gamma}) - \frac{p_{\alpha\gamma}}{\rho}\right) - \mathbf{v}_{\gamma}^{eq} \frac{p_{\alpha\beta}}{\rho} - \mathbf{v}_{\alpha}^{eq} \mathbf{v}_{\beta}^{eq} \mathbf{v}_{\gamma}^{eq}] \qquad (A13) $$

b. At $\epsilon^2$ scale: One can obtain the second order perturbation of the zero-th moment equation:

$$ \partial_t^{(2)} (\rho) + 0 = \frac{2\Delta t}{2-\hat{\omega}\Delta t} \sum_i S_i^{(2)}, \qquad (A14) $$

and the second order perturbation of the first moment equation:

$$ \partial_t^{(2)} (\rho \mathbf{v}^{eq}) + \partial_x^{(1)} \cdot (\Pi^{(1)}) = \frac{2\Delta t}{2-\hat{\omega}\Delta t} M_1 (S_i^{(2)}). \qquad (A15) $$
With the notation $D_t^{*(1)}(M_n(\hat{f}_i^{(k)})) = \partial_t^{(1)}(M_n(\hat{f}_i^{(k)})) + \partial_x^{(1)}(M_{n+1}(\hat{f}_i^{(k)}))$, the second moment of the second perturbation order equation yields:

$$ \partial_t^{(2)} (\rho \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} + p) + D_t^{*(1)} (\Pi^{(1)}) = A M_2 (S_i^{(2)}) - A \hat{\omega} M_2 (f_i^{(2)}), \quad (A16) $$

which itself leads to expression of the second moment of the second perturbation order distribution:

$$ M_{2,\alpha\beta} (f_i^{(2)}) = \frac{1}{A\hat{\omega}} \left\{ \partial_{\mathbf{x}_{\gamma}}^{(1)} \cdot (\Pi_{\alpha\gamma}^{(1)}) \otimes \mathbf{v}_{\beta}^{eq} + \mathbf{v}_{\alpha}^{eq} \otimes \partial_{\mathbf{x}_{\gamma}}^{(1)} \cdot (\Pi_{\beta\gamma}^{(1)}) - D_t^{*(1)} [M_{2,\alpha\beta} (f_i^{(1)})] \right\} \\ + \frac{1}{\hat{\omega}} \left\{ M_{2,\alpha\beta} (S_i^{(2)}) \right\} \\ + \frac{1}{\hat{\omega}} \left\{ M_0 (S_i^{(2)}) [\mathbf{v}_{\alpha}^{eq} \otimes \mathbf{v}_{\beta}^{eq} - \partial_{\rho} (p_{\alpha\beta}(\rho))] - M_{1,\alpha} (S_i^{(2)}) \otimes \mathbf{v}_{\beta}^{eq} - \mathbf{v}_{\alpha}^{eq} \otimes M_{1,\beta} (S_i^{(2)}) \right\}. \quad (A17) $$

c. At $\epsilon^3$ scale: The third perturbation order, which depends on previous computed moments, yields the following zero-th moment:

$$ \partial_t^{(3)} (\rho) - \frac{1/12}{\tau_2} \partial_x^{(1)} \cdot \partial_x^{(1)} \cdot (\Pi^{(1)}) - \frac{1}{12} A D_t^{*(1)2} M_0 (S_i^{(1)}) + \left(1 - \frac{1}{\hat{\omega}}\right) D_t^{*(1)} M_0 (S_i^{(2)}) = 0, \quad (A18) $$

and the following first moment:

$$ \partial_t^{(3)} (\rho \mathbf{v}^{eq}) - \frac{1/12}{\tau_2} \partial_t^{(1)} \partial_x^{(1)} \cdot (\Pi^{(1)}) - \frac{1/12}{\tau_2} \partial_x^{(1)} \cdot D_t^{*(1)} (M_2 [f_i^{(1)}]) + \partial_x^{(1)} \cdot (M_2 [f_i^{(2)}]) \\ - \frac{1}{12} A D_t^{*(1)2} (M_1 [S_i^{(1)}]) + \left(1 - \frac{1}{\hat{\omega}}\right) D_t^{*(1)} (M_1 [S_i^{(2)}]) = 0. \quad (A19) $$

d. Recomposition ($\epsilon^1+\epsilon^2+\epsilon^3$) scale: To obtain the recovered equations, one should recompose the time/perturbation scales, which leads to the mass conservation up to the third perturbation order:

$$ \partial_t (\rho) + \partial_x \cdot (\rho \mathbf{v}^{eq}) - \frac{1/12}{\tau_2} \partial_x \cdot \partial_x \cdot (\Pi^{(1)}) \\ - \frac{1}{12} A D_t^{*2} (M_0 [S_i^{(1)}]) + \left(1 - \frac{1}{\hat{\omega}}\right) D_t^* (M_0 [S_i^{(2)}]) = A M_0 [S_i^{(1)} + S_i^{(2)}], \quad (A20) $$

and the associated momentum equation up to the third perturbation order:

$$ \partial_t (\rho \mathbf{v}^{eq}) + \partial_x \cdot (\rho \mathbf{v}^{eq} \otimes \mathbf{v}^{eq} + p + \Pi^{(1)}) - \frac{1/12}{\tau_2} \partial_t \partial_x \cdot (\Pi^{(1)}) - \frac{1/12}{\tau_2} \partial_x \cdot D_t^* (M_2 [f_i^{(1)}]) \\ + \partial_x \cdot (M_2 [f_i^{(2)}]) - \frac{1}{12} A D_t^{*2} (M_1 [S_i^{(1)}]) + \left(1 - \frac{1}{\hat{\omega}}\right) D_t^* (M_1 [S_i^{(2)}]) = A M_1 [S_i^{(1)} + S_i^{(2)}]. \quad (A21) $$
e. Using GUO forcing scheme Once the forcing scheme is specified, the previous equations
can be adapted to such scheme. When using GUO's forcing scheme, the different moments are:

(A22)

$$
S_i^{(1)} = \left(1 - \frac{\hat{\omega} \Delta t}{2}\right) w_i \left( \frac{\boldsymbol{\xi}_i - \mathbf{v}^{\mathrm{eq}}}{c_s^2} + \frac{(\mathbf{v}^{\mathrm{eq}} \cdot \boldsymbol{\xi}_i) \boldsymbol{\xi}_i}{c_s^4} \right) \cdot \mathbf{F} \quad (\mathrm{A}23)
$$

$$
S_i^{(2)} = \tilde{f}(\rho, \mathbf{v}, p + \tilde{K}) - f_i^{eq}(\rho, \mathbf{v}, p) \quad (A24)
$$

$$
M_0 \left[ S_i^{(1)} \right] = 0 \qquad (A25)
$$

$$
M_1 \left[ S_i^{(1)} \right] = \left( 1 - \frac{\hat{\omega} \Delta t}{2} \right) \mathbf{F} = \frac{1}{A \Delta t} \mathbf{F} \qquad (\text{A26})
$$

$$
M_2 \left[ S_i^{(1)} \right] = \left( 1 - \frac{\hat{\omega} \Delta t}{2} \right) \left( \mathbf{v} \otimes \mathbf{F} + \mathbf{F} \otimes \mathbf{v} \right) = \frac{1}{A \Delta t} \left( \mathbf{v} \otimes \mathbf{F} + \mathbf{F} \otimes \mathbf{v} \right) \quad (\text{A27})
$$

$$
M_3 \left[ S_i^{(1)} \right] = \frac{c_s^2}{A \Delta t} P_3 \left( \mathbf{F}_\alpha \otimes I_{\beta\gamma} \right) = \frac{c_s^2}{A \Delta t} \left( \mathbf{F}_\alpha \otimes I_{\beta\gamma} + \mathbf{F}_\beta \otimes I_{\alpha\gamma} + \mathbf{F}_\gamma \otimes I_{\alpha\beta} \right) \quad (A28)
$$

$$
M_0 \left[ S_i^{(2)} \right] = 0 \qquad (A29)
$$

$$
M_1 \left[ S_i^{(2)} \right] = 0 \qquad (A30)
$$

$$
M_2 \left[ S_i^{(2)} \right] = \tilde{K} \quad (\text{free parameter by construction}). \qquad (A31)
$$

These moments from the source terms lead to non equilibrium moments variables:

$$
M_2 [f_i^{(1)}] = \Pi^{(1)} = -\frac{1}{A\hat{\omega}} \left\{ p \cdot \partial_x^{(1)} (\mathbf{v}^{\mathrm{eq}}) + \left(p \cdot \partial_x^{(1)} (\mathbf{v}^{\mathrm{eq}})\right)^T + \partial_x^{(1)} (\mathbf{v}^{\mathrm{eq}}) [p - \rho \partial_\rho (p)] \right\} + 0, \quad (1)
$$

(A32)

$$
\begin{equation}
\begin{aligned}
M_{3,\alpha\beta\gamma} [f_i^{(1)}] ={}& \partial_{\mathbf{x}_\delta}^{(1)} \cdot \left( p_{\alpha\beta} p_{\gamma\delta} + \mathbf{v}_\alpha^{\text{eq}} \mathbf{v}_\beta^{\text{eq}} p_{\gamma\delta} + \mathbf{v}_\alpha^{\text{eq}} \mathbf{v}_\delta^{\text{eq}} p_{\gamma\beta} + \mathbf{v}_\beta^{\text{eq}} \mathbf{v}_\delta^{\text{eq}} p_{\gamma\alpha} \right) \\
& + \left[ \frac{1}{\Delta t} (\mathbf{v}_\alpha \mathbf{F}_\beta + \mathbf{F}_\alpha \mathbf{v}_\beta) + p_{\alpha\delta} \cdot \partial_{\mathbf{x}_\delta}^{(1)} (\mathbf{v}_\beta^{\text{eq}}) + (p_{\beta\delta} \cdot \partial_{\mathbf{x}_\delta}^{(1)} (\mathbf{v}_\alpha^{\text{eq}})) + \partial_{\mathbf{x}_\delta}^{(1)} \cdot (\mathbf{v}_\delta^{\text{eq}}) [p_{\alpha\beta} - \rho \partial_\rho (p_{\alpha\beta})] \right] \mathbf{v}_\gamma^{\text{eq}} \\
& + \left[ \frac{1}{\Delta t} \mathbf{F}_\alpha - \partial_{\mathbf{x}_\delta}^{(1)} \cdot (\rho \mathbf{v}_\alpha^{\text{eq}} \mathbf{v}_\delta^{\text{eq}} + p_{\alpha\delta}) \right] \left( \mathbf{v}_\beta^{\text{eq}} \mathbf{v}_\gamma^{\text{eq}} + \frac{p_{\beta\gamma}}{\rho} \right) + \left[ \rho \mathbf{v}_\alpha^{\text{eq}} \mathbf{v}_\beta^{\text{eq}} \mathbf{v}_\delta^{\text{eq}} + P_3 (\mathbf{v}_\alpha^{\text{eq}} p_{\beta\delta}) \right] \cdot \partial_{\mathbf{x}_\delta}^{(1)} (\mathbf{v}_\gamma^{\text{eq}}) \\
& + \left[ \frac{1}{\Delta t} \mathbf{F}_\beta - \partial_{\mathbf{x}_\delta}^{(1)} \cdot (\rho \mathbf{v}_\beta^{\text{eq}} \mathbf{v}_\delta^{\text{eq}} + p_{\beta\delta}) \right] \frac{p_{\alpha\gamma}}{\rho} + \left[ \frac{1}{\Delta t} \mathbf{F}_\gamma - \partial_{\mathbf{x}_\delta}^{(1)} \cdot (\rho \mathbf{v}_\gamma^{\text{eq}} \mathbf{v}_\delta^{\text{eq}} + p_{\gamma\delta}) \right] \frac{p_{\alpha\beta}}{\rho} \\
& - \partial_{\mathbf{x}_\delta}^{(1)} \cdot (\rho \mathbf{v}_\delta^{\text{eq}}) \left[ \mathbf{v}_\alpha^{\text{eq}} \left( \partial_\rho (p_{\beta\gamma}) - \frac{p_{\beta\gamma}}{\rho} \right) + \mathbf{v}_\beta^{\text{eq}} \left( \partial_\rho (p_{\alpha\gamma}) - \frac{p_{\alpha\gamma}}{\rho} \right) - \mathbf{v}_\gamma^{\text{eq}} \frac{p_{\alpha\beta}}{\rho} - \mathbf{v}_\alpha^{\text{eq}} \mathbf{v}_\beta^{\text{eq}} \mathbf{v}_\gamma^{\text{eq}} \right],
\end{aligned}
\tag{A33}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
M_{2,\alpha\beta} [f_i^{(2)}] ={}& \frac{1}{A\hat{\omega}} \left\{ \partial_{\mathbf{x}_\gamma}^{(1)} \cdot \left( \Pi_{\alpha\gamma}^{(1)} \right) \otimes \mathbf{v}_\beta^{\text{eq}} + \mathbf{v}_\alpha^{\text{eq}} \otimes \partial_{\mathbf{x}_\gamma}^{(1)} \cdot \left( \Pi_{\beta\gamma}^{(1)} \right) \right\} \\
& - \frac{1}{A\hat{\omega}} \left\{ \partial_t^{(1)} \left( \Pi_{\alpha\beta}^{(1)} \right) + \partial_{\mathbf{x}_\gamma}^{(1)} \cdot \left( M_{3,\alpha\beta\gamma} (f_i^{(1)}) \right) \right\} + \frac{\tilde{K}_{\alpha\beta}}{\hat{\omega}}.
\end{aligned}
\tag{A34}
\end{equation}
$$
f. Using EDM forcing scheme Similarly for EDM's forcing scheme, the source terms moments are:

$$S_i^{(1)} = f_i^{eq}(\rho, \mathbf{v} + \frac{\Delta t \mathbf{F}}{\rho}) - f_i^{eq}(\rho, \mathbf{v}) \qquad (A35)$$

$$S_i^{(2)} = \tilde{f}(\rho, \mathbf{v}, p + \tilde{K}) - f_i^{eq}(\rho, \mathbf{v}, p) \qquad (A36)$$

$$M_0 [S_i^{(1)}] = 0 \qquad (A37)$$

$$M_{1,\alpha} [S_i^{(1)}] = \mathbf{F}_{\alpha} \Delta t \qquad (A38)$$

$$M_{2,\alpha\beta} [S_i^{(1)}] = \mathbf{F}_{\alpha} \Delta t \otimes \frac{\mathbf{F}_{\beta} \Delta t}{\rho} + \mathbf{v}_{\alpha} \otimes \mathbf{F}_{\beta} \Delta t + \mathbf{F}_{\alpha} \Delta t \otimes \mathbf{v}_{\beta} \qquad (A39)$$

$$M_{3,\alpha\beta\gamma} [S_i^{(1)}] = \rho \left( \mathbf{v}_{\alpha} + \frac{\mathbf{F}_{\alpha} \Delta t}{\rho} \right) \otimes \left( \mathbf{v}_{\beta} + \frac{\mathbf{F}_{\beta} \Delta t}{\rho} \right) \otimes \left( \mathbf{v}_{\gamma} + \frac{\mathbf{F}_{\gamma} \Delta t}{\rho} \right) - \rho \mathbf{v}_{\alpha} \otimes \mathbf{v}_{\beta} \otimes \mathbf{v}_{\gamma} + P_3 \left( \frac{\mathbf{F}_{\alpha} \Delta t}{\rho} \otimes \mathbf{p}_{\beta\gamma} \right) \qquad (A40)$$

$$M_0 [S_i^{(2)}] = 0 \qquad (A41)$$

$$M_1 [S_i^{(2)}] = 0 \qquad (A42)$$

$$M_{2,\alpha\beta} [S_i^{(2)}] = \tilde{K}_{\alpha\beta} \text{ (free parameter by construction).} \qquad (A43)$$

Therefore, the resulting non-equilibrium moments are:

$$M_{2,\alpha\beta} [f_i^{(1)}] = \Pi_{\alpha\beta}^{(1)} = -\frac{1}{A\hat{\omega}} \{V_{\alpha\beta}\} + \frac{1}{\hat{\omega}} \mathbf{F}_{\alpha} \Delta t \otimes \frac{\mathbf{F}_{\beta} \Delta t}{\rho}, \qquad (A44)$$

$$M_{3,\alpha\beta\gamma} [f_i^{(1)}] = [V_{\alpha\beta} + A\mathbf{v}_{\alpha} \otimes \mathbf{F}_{\beta} \Delta t + A\mathbf{F}_{\alpha} \Delta t \otimes \mathbf{v}_{\beta}] \mathbf{v}_{\gamma}^{eq} \\ + \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot \left( p_{\alpha\beta}p_{\gamma\delta} + \mathbf{v}_{\alpha}^{eq}\mathbf{v}_{\beta}^{eq}p_{\gamma\delta} + \mathbf{v}_{\alpha}^{eq}\mathbf{v}_{\delta}^{eq}p_{\gamma\beta} + \mathbf{v}_{\beta}^{eq}\mathbf{v}_{\delta}^{eq}p_{\gamma\alpha} \right) \\ + \left[ A\mathbf{F}_{\alpha} \Delta t - \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot (\rho\mathbf{v}_{\alpha}^{eq}\mathbf{v}_{\delta}^{eq} + p_{\alpha\delta}) \right] \left( \mathbf{v}_{\beta}^{eq}\mathbf{v}_{\gamma}^{eq} + \frac{p_{\beta\gamma}}{\rho} \right) + \left[ \rho\mathbf{v}_{\alpha}^{eq}\mathbf{v}_{\beta}^{eq}\mathbf{v}_{\delta}^{eq} + P_3(\mathbf{v}_{\alpha}^{eq}p_{\beta\delta}) \right] \cdot \partial_{\mathbf{x}_{\delta}}^{(1)} (\mathbf{v}_{\gamma}^{eq}) \\ + \left[ A\mathbf{F}_{\beta} \Delta t - \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot (\rho\mathbf{v}_{\beta}^{eq}\mathbf{v}_{\delta}^{eq} + p_{\beta\delta}) \right] \frac{p_{\alpha\gamma}}{\rho} + \left[ A\mathbf{F}_{\gamma} \Delta t - \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot (\rho\mathbf{v}_{\gamma}^{eq}\mathbf{v}_{\delta}^{eq} + p_{\gamma\delta}) \right] \frac{p_{\alpha\beta}}{\rho} \\ - \partial_{\mathbf{x}_{\delta}}^{(1)} \cdot (\rho\mathbf{v}_{\delta}^{eq}) \left[ \mathbf{v}_{\alpha}^{eq} \left( \partial_{\rho}(p_{\beta\gamma}) - \frac{p_{\beta\gamma}}{\rho} \right) + \mathbf{v}_{\beta}^{eq} \left( \partial_{\rho}(p_{\alpha\gamma}) - \frac{p_{\alpha\gamma}}{\rho} \right) - \mathbf{v}_{\gamma}^{eq} \frac{p_{\alpha\beta}}{\rho} - \mathbf{v}_{\alpha}^{eq}\mathbf{v}_{\beta}^{eq}\mathbf{v}_{\gamma}^{eq} \right], \qquad (A45)$$

$$M_{2,\alpha\beta} [f_i^{(2)}] = -\frac{1}{(A\hat{\omega})^2} \left\{ \partial_{\mathbf{x}_{\gamma}}^{(1)} \cdot (V_{\alpha\gamma}) \otimes \mathbf{v}_{\beta}^{eq} + \mathbf{v}_{\alpha}^{eq} \otimes \partial_{\mathbf{x}_{\gamma}}^{(1)} \cdot (V_{\beta\gamma}) \right\} \\ + \frac{1}{A(\hat{\omega})^2} \left\{ \partial_{\mathbf{x}_{\gamma}}^{(1)} \cdot \left( \mathbf{F}_{\alpha} \Delta t \otimes \frac{\mathbf{F}_{\gamma} \Delta t}{\rho} \right) \otimes \mathbf{v}_{\beta}^{eq} + \mathbf{v}_{\alpha}^{eq} \otimes \partial_{\mathbf{x}_{\gamma}}^{(1)} \cdot \left( \mathbf{F}_{\beta} \Delta t \otimes \frac{\mathbf{F}_{\gamma} \Delta t}{\rho} \right) \right\} \\ + \frac{1}{(A\hat{\omega})^2} \partial_t^{(1)} (V_{\alpha\beta}) - \frac{1}{A(\hat{\omega})^2} \partial_t^{(1)} \left( \mathbf{F}_{\alpha} \Delta t \otimes \frac{\mathbf{F}_{\beta} \Delta t}{\rho} \right) - \frac{1}{A\hat{\omega}} \partial_{\mathbf{x}_{\gamma}}^{(1)} \cdot \left( M_{3,\alpha\beta\gamma} [f_i^{(1)}] \right) \\ + \frac{\tilde{K}_{\alpha\beta}}{\hat{\omega}} \qquad (A46)$$
g. Simplification at steady-state: Following the methodology from [28], we can simplify the problem by introducing the following steady-state hypothesis:

$$ \bullet \mathbf{v}^{eq} = 0, $$

$$ \bullet \partial_t (*) = 0. $$

With these assumptions eq. (A21), for Guo's scheme, becomes:

$$ \begin{aligned} \partial_t (\rho \mathbf{v}_\alpha^{eq}) + \partial_{\mathbf{x}_\beta} \cdot \left( \rho \mathbf{v}_\alpha^{eq} \otimes \mathbf{v}_\beta^{eq} + p_{\alpha\beta} - \frac{1}{A\hat{\omega}} (\rho c_s^2 \partial_{\mathbf{x}_\alpha} (\mathbf{v}_\beta^{eq}) + \rho c_s^2 \partial_{\mathbf{x}_\alpha} (\mathbf{v}_\beta^{eq})) \right) \\ + \frac{\tau_3}{\tau_2} \partial_{\mathbf{x}_\beta} \cdot \partial_{\mathbf{x}_\gamma} \cdot \left( M_{3,\alpha\beta\gamma} [f_i^{(1)}] \right) \\ - \frac{c_s^2}{12\Delta t} \partial_{\mathbf{x}_\beta} \cdot \partial_{\mathbf{x}_\gamma} \cdot (\mathbf{F}_\alpha \otimes I_{\beta\gamma} + \mathbf{F}_\beta \otimes I_{\alpha\gamma} + \mathbf{F}_\gamma \otimes I_{\alpha\beta}) + (\partial_{\mathbf{x}_\beta} \cdot \tilde{K}_{\alpha\beta}) = \frac{1}{\Delta t} \mathbf{F}. \end{aligned} \quad (A47) $$

and similarly, for EDM, the same equation becomes:

$$ \begin{aligned} \partial_t (\rho \mathbf{v}_\alpha^{eq}) + \partial_{\mathbf{x}_\beta} \cdot \left( \rho \mathbf{v}_\alpha^{eq} \otimes \mathbf{v}_\beta^{eq} + p_{\alpha\beta} - \frac{1}{A\hat{\omega}} (\rho c_s^2 \partial_{\mathbf{x}_\alpha} (\mathbf{v}_\beta^{eq}) + \rho c_s^2 \partial_{\mathbf{x}_\alpha} (\mathbf{v}_\beta^{eq})) \right) \\ + \frac{\tau_3}{\tau_2} \partial_{\mathbf{x}_\beta} \cdot \partial_{\mathbf{x}_\gamma} \cdot \left( M_{3,\alpha\beta\gamma} [f_i^{(1)}] \right) \\ + \frac{1}{\hat{\omega}} \partial_{\mathbf{x}_\beta} \cdot \left( \mathbf{F}_\alpha \Delta t \otimes \frac{\mathbf{F}_\beta \Delta t}{\rho} \right) \\ - \frac{A}{12} \partial_{\mathbf{x}_\beta} \cdot \partial_{\mathbf{x}_\gamma} \cdot \left\{ \rho \left( \frac{\mathbf{F}_\alpha \Delta t}{\rho} \right) \otimes \left( \frac{\mathbf{F}_\beta \Delta t}{\rho} \right) \otimes \left( \frac{\mathbf{F}_\gamma \Delta t}{\rho} \right) + P_3 \left( \frac{\mathbf{F}_\alpha \Delta t}{\rho} \otimes p_{\beta\gamma} \right) \right\} \\ + \partial_{\mathbf{x}_\beta} \cdot \tilde{K}_{\alpha\beta} = A \mathbf{F}_\alpha \Delta t \end{aligned} \quad (A48) $$

h. Expanding pseudo-potential forces with Local Mean Gradient Approximation Using Kupershtokh Local Mean Gradient Approximation (LMGA) [24] force expression:

$$ \mathbf{F}_{\alpha} = (1 - 2B) \mathbf{F}_{\alpha}^{\text{local}} + B \mathbf{F}_{\alpha}^{\text{mean}}, \quad (A49) $$

with

$$ \mathbf{F}_{\alpha}^{\text{local}} = \psi(\mathbf{x}) G_1 \sum_{i=0}^{q} W(|\xi_i|^2) \psi(\mathbf{x} + \xi_i \Delta t) \xi_i, \quad (A50) $$

and

$$ \mathbf{F}_{\alpha}^{\text{mean}} = G_1 \sum_{i=0}^{q} W(|\xi_i|^2) \psi^2(\mathbf{x} + \xi_i \Delta t) \xi_i . \quad (A51) $$

This forces can be rewritten using a third order Taylor expansion:

$$ \mathbf{F}_{\alpha}^{\text{local}} = G_1 \lambda_1 \psi \nabla_{\mathbf{x}_{\alpha}} (\psi) + G_1 \lambda_2 \psi \nabla_{\mathbf{x}_{\alpha}} (\nabla^2_{\mathbf{x}_{\beta}} \psi) + \mathcal{O}(\Delta \mathbf{x}^3), $$
$$ F_{\alpha}^{\text{local}} = -\partial_{x_{\beta}} \cdot \left\{ \left( -\frac{G_1 \lambda_1}{2} \psi^2 - G_1 \lambda_2 \psi \partial_{x_{\gamma}}^2 (\psi) - \frac{G_1 \lambda_2}{2} |\partial_{x_{\gamma}} (\psi)|^2 \right) I_{\alpha\beta} + G_1 \lambda_2 \partial_{x_{\alpha}} (\psi) \otimes \partial_{x_{\beta}} (\psi) \right\} + \mathcal{O}(\Delta x^3) . \quad (A52) $$

$$ F_{\alpha}^{\text{mean}} = G_1 \lambda_1 \partial_{x_{\alpha}} (\psi^2) + G_1 \lambda_2 \partial_{x_{\alpha}} (\partial_{x_{\gamma}}^2 (\psi^2)) + \mathcal{O}(\Delta x^3) , $$

$$ F_{\alpha}^{\text{mean}} = -\partial_{x_{\beta}} \cdot \left\{ \left( -G_1 \lambda_1 \psi^2 - 2G_1 \lambda_2 \psi \partial_{x_{\gamma}}^2 (\psi) - 2G_1 \lambda_2 |\partial_{x_{\gamma}} (\psi)|^2 \right) I_{\alpha\beta} \right\} + \mathcal{O}(\Delta x^3) . \quad (A53) $$

Incorporating these forcing expression eqs. (A52) to (A53) into eq. (A47) (and eq. (A48) respectively), one can obtain the expression of stress tensor source $\tilde{K}$ to verify eq. (8) given in eq. (A54) (and eq. (A55) respectively).

i. Local Mean Gradient Approximation pseudo-potential forces with GUO's scheme:

$$ \tilde{K}_{\alpha\beta} = \left( \left[ \frac{G_1 \lambda_1}{2} - \frac{G}{2} \right] \psi^2 + \left[ G_1 \lambda_2 \psi + \frac{G_1 \lambda_1 B c_s^2}{6\Delta t} \psi - \kappa \rho \frac{d\psi}{d\rho} \right] \partial_{x_{\gamma}}^2 (\psi) + \left[ (1+2B) \frac{G_1 \lambda_2}{2} + \frac{G_1 \lambda_1 B c_s^2}{6\Delta t} - \frac{\kappa}{2} \right] |\partial_{x_{\gamma}} (\psi)|^2 \right. \\ \left. + \left[ \kappa - (1-2B)G_1 \lambda_2 + \frac{G_1 \lambda_1 B c_s^2}{3\Delta t} \right] \partial_{x_{\alpha}} (\psi) \otimes \partial_{x_{\beta}} (\psi) + \left[ \frac{4G_1 \lambda_1 B c_s^2}{12\Delta t} \right] \psi \partial_{x_{\alpha}} (\partial_{x_{\beta}} (\psi)) \right. \\ \left. - \frac{\tau_3}{\tau_2} \partial_{x_{\gamma}} \cdot \left( M_{3,\alpha\beta\gamma} [f_i^{(1)}] \right) \right) . \quad (A54) $$

j. Local Mean Gradient Approximation pseudo-potential forces with EDM scheme:

$$ \tilde{K}_{\alpha\beta} = \left( \left[ \frac{G_1 \lambda_1}{2} - \frac{G}{2} \right] \psi^2 + \left[ \psi G_1 \lambda_2 - \kappa \rho \frac{d\psi}{d\rho} \right] \partial_{x_{\gamma}}^2 (\psi) + \left[ (1+2B) \frac{G_1 \lambda_2}{2} - \frac{\kappa}{2} \right] |\partial_{x_{\gamma}} (\psi)|^2 \right) I_{\alpha\beta} \\ + \left[ \kappa - (1+2B)G_1^2 \lambda_1^2 \frac{\Delta t^2}{\hat{\omega}} \frac{\psi^2}{\rho} - (1-2B)G_1 \lambda_2 \right] \partial_{x_{\alpha}} (\psi) \otimes \partial_{x_{\beta}} (\psi) \\ + \frac{A}{12} \partial_{x_{\gamma}} \cdot \left\{ \frac{G_1 \lambda_1 \Delta t}{\rho} P_3 (\psi \nabla_{x_{\alpha}} (\psi) \otimes p_{\beta\gamma}) \right\} - \frac{\tau_3}{\tau_2} \partial_{x_{\gamma}} \cdot \left( M_{3,\alpha\beta\gamma} [f_i^{(1)}] \right) . \quad (A55) $$

## Appendix B: The results of LBGK/KBC+EDM framework

As we described in section II, the LBGK model is a simplification of KBC model with the simple form of $f_i^{\text{mirr}} = 2f_i^{\text{eq}} - f_i$. Therefore, to make the model more general, we implement the KBC+EDM framework referring to [36, 67]. The tuning coefficient to achieve the thermodynamic consistency is eq. (40) tested as $B = -0.17$. Similar to section IVA, we model the stationary liquid droplet suspended in a vapor environment. We use a 2D simulation domain of 240 × 240 lattices² with droplet radius of 48 lattices, while periodic boundary is applied in both x and y directions. Liquid and vapor kinematic viscosity are set as $\nu_l = 0.01$ and $\nu_v = 1/6$, respectively.
<figure><img src="image_9.png" /><figcaption>Figure 8. Modeling of liquid droplet suspended in vapor with (a) density contour and (b) spurious current using the KBC+EDM framework.</figcaption></figure>

The temperature ratio is selected at $T/T_c = 0.75$, with the obtained liquid-vapor density ratio of $\rho_l/\rho_v = 22.0$. We first apply $\kappa = 0$ in the current simulation, and the obtained result is shown in fig. 8. Figure 8a shows the density contour with a diffusive interface, while fig. 8b shows the spurious current with the maximum value around $8.0 \times 10^{-3}$, which is approximately two times of that in the MRT+GUO framework (fig. 2b).

We further vary the values of $\kappa$ to obtain the tunable range of the surface tension. As shown in fig. 9a, when $\kappa$ ranges from -2.5 to 2.0, the surface tension ratio $\sigma/\sigma_{\kappa=0}$ varies from 0.26 to 1.64. The tunable surface tension range is around 6 times, which is much lower than 140 times in the MRT+GUO framework. But still, an approximately linear relation of $\sigma/\sigma_{\kappa=0} = 1.013 + 0.306\kappa$ is obtained between $\kappa$ and $\sigma/\sigma_{\kappa=0}$. In terms of the spurious current, as shown in fig. 9b, the maximum values corresponding to $\kappa \in (-2.5, 2.0)$ are in the range of $(4.9 \times 10^{-3}, 2.9 \times 10^{-2})$, which is also higher than $(3.1 \times 10^{-3}, 1.0 \times 10^{-2})$ in the MRT+GUO framework. For the liquid-vapor density ratio, it decreases from $\rho_l/\rho_v = 31.5$ at $\kappa = -2.5$ to $\rho_l/\rho_v = 14.7$ at $\kappa = 2.0$ with the variation range of around 55%, which is a little higher than 45% in the MRT+GUO framework. From the above results we can observe that, the proposed source term to tune the surface tension works better in the MRT+GUO framework, rather than the LBGK/KBC+EDM framework. On the other hand, the results of both MRT+GUO and LBGK/KBC+EDM frameworks show that the proposed model is applicable independently of collision operators and force schemes.
<figure><img src="image_10.png" /><figcaption>Figure 9. The variation of (a) surface tension ratio, (b) maximum spurious current with parameter $κ$ using the KBC+EDM framework.</figcaption></figure>

[1] D. Richard, C. Clanet, and D. Quéré, Contact time of a bouncing drop, Nature **417**, 811 (2002).

[2] Q. Li, K. H. Luo, Q. J. Kang, Y. L. He, Q. Chen, and Q. Liu, Lattice Boltzmann methods for multiphase flow and phase-change heat transfer, Prog. Energy Combust. Sci. **52**, 62 (2016).

[3] G. R. McNamara and G. Zanetti, Use of the boltzmann equation to simulate lattice-gas automata, Phys. Rev. Lett. **61**, 2332 (1988).

[4] Y. H. Qian, D. D'Humières, and P. Lallemand, Lattice BGK Models for Navier-Stokes Equation, Europhysics Letters (EPL) **17**, 479 (1992).

[5] S. Chen and G. D. Doolen, Lattice Boltzmann Method for Fluid Flows, Annu. Rev. Fluid Mech. **30**, 329 (1998).

[6] S. Succi, *The lattice Boltzmann equation: for complex states of flowing matter* (Oxford University Press, 2018).

[7] Z. Huang, Y. Zhu, P. Jiang, and Y. Xiong, Investigation of a porous transpiration-cooled strut injector, Journal of Propulsion and Power **31**, 278 (2015).

[8] S. Succi, Lattice boltzmann 2038, EPL (Europhysics Letters) **109**, 50001 (2015).

[9] L. Fei, J. Yang, Y. Chen, H. Mo, and K. H. Luo, Mesoscopic simulation of three-dimensional pool boiling based on a phase-change cascaded lattice boltzmann method, Physics of Fluids **32**, 103312 (2020).

[10] G. Wang, L. Fei, and K. H. Luo, Unified lattice boltzmann method with improved schemes for multi-phase flow simulation: Application to droplet dynamics under realistic conditions, Physical Review E **105**, 045314 (2022).

[11] A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti, Lattice Boltzmann model of immiscible
fluids, Phys. Rev. A **43**, 4320 (1991).

[12] D. Grunau, S. Chen, and K. Eggert, A lattice Boltzmann model for multiphase fluid flows, Phys. Fluids A Fluid Dyn. **5**, 2557 (1993).

[13] X. Shan and H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E **47**, 1815 (1993).

[14] X. Shan and H. Chen, Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E **49**, 2941 (1994).

[15] M. R. Swift, W. R. Osborn, and J. M. Yeomans, Lattice Boltzmann Simulation of Nonideal Fluids, Phys. Rev. Lett. **75**, 830 (1995).

[16] M. Swift, E. Orlandini, W. Osborn, and J. Yeomans, Lattice Boltzmann simulations of liquid-gas and binary fluid systems, Phys. Rev. E **54**, 5041 (1996).

[17] X. He, S. Chen, and R. Zhang, A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh-Taylor Instability, J. Comput. Phys. **152**, 642 (1999).

[18] H. Liang, J. Xu, J. Chen, H. Wang, Z. Chai, and B. Shi, Phase-field-based lattice boltzmann modeling of large-density-ratio two-phase flows, Physical Review E **97**, 033309 (2018).

[19] D. Zhang, Y. Li, Y. Wang, and C. Shu, An improved multiphase lattice boltzmann flux solver with a modified cahn-hilliard equation for multiphase flow with super large density ratio, Physics of Fluids **36**, 10.1063/5.0189032 (2024).

[20] X. Shan, Analysis and reduction of the spurious current in a class of multiphase lattice Boltzmann models, Phys. Rev. E **73**, 047701 (2006).

[21] X. Shan, Pressure tensor calculation in a class of nonideal gas lattice Boltzmann models, Phys. Rev. E **77**, 066702 (2008).

[22] M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama, and F. Toschi, Generalized lattice Boltzmann method with multirange pseudopotential, Phys. Rev. E **75**, 026702 (2007).

[23] S. F. Kharmiani, H. Niazmand, and M. Passandideh-Fard, An alternative high-density ratio pseudopotential lattice boltzmann model with surface tension adjustment capability, Journal of Statistical Physics **175**, 47 (2019).

[24] A. L. Kupershtokh, D. A. Medvedev, and D. I. Karpov, On equations of state in a lattice Boltzmann method, Computers and Mathematics with Applications **58**, 965 (2009).

[25] Q. Li, K. H. Luo, and X. J. Li, Forcing scheme in pseudopotential lattice Boltzmann model for multiphase flows, Physical Review E **86**, 016709 (2012).

[26] Q. Li, K. H. Luo, and X. J. Li, Lattice Boltzmann modeling of multiphase flows at large density ratio with an improved pseudopotential model, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics **87**, 10.1103/PhysRevE.87.053301 (2013), arXiv:1211.6932.

[27] D. Lycett-Brown and K. H. Luo, Improved forcing scheme in pseudopotential lattice Boltzmann methods for multiphase flow at arbitrarily high density ratios, Physical Review E **91**, 023305 (2015).

[28] R. Huang and H. Wu, Third-order analysis of pseudopotential lattice Boltzmann model for multiphase
flow, Journal of Computational Physics **327**, 121 (2016).

[29] P. Yuan and L. Schaefer, Equations of state in a lattice Boltzmann model, Physics of Fluids **18**, 10.1063/1.2187070 (2006).

[30] C. E. Colosqui, G. Falcucci, S. Ubertini, and S. Succi, Mesoscopic simulation of non-ideal fluids with self-tuning of the equation of state, Soft matter **8**, 3798 (2012).

[31] Q. Li and K. Luo, Thermodynamic consistency of the pseudopotential lattice boltzmann model for simulating liquid-vapor flows, Applied Thermal Engineering **72**, 56 (2014).

[32] C. Peng, L. F. Ayala, Z. Wang, and O. M. Ayala, Attainment of rigorous thermodynamic consistency and surface tension in single-component pseudopotential lattice boltzmann models via a customized equation of state, Physical Review E **101**, 063309 (2020).

[33] D. Lycett-Brown and K. H. Luo, Cascaded lattice Boltzmann method with improved forcing scheme for large-density-ratio multiphase flow at high Reynolds and Weber numbers, Phys. Rev. E **94**, 053313 (2016).

[34] Q. Li, D. Du, L. Fei, and K. H. Luo, Three-dimensional non-orthogonal MRT pseudopotential lattice Boltzmann model for multiphase flows, Comput. & Fluids **186**, 128 (2019).

[35] L. Fei, J. Du, K. H. Luo, S. Succi, M. Lauricella, A. Montessori, and Q. Wang, Modeling realistic multiphase flows using a non-orthogonal multiple-relaxation-time lattice Boltzmann method, Phys. Fluids **31** (2019).

[36] F. Qin, A. Mazloomi Moqaddam, Q. Kang, D. Derome, and J. Carmeliet, Entropic multiple-relaxation-time multirange pseudopotential lattice Boltzmann model for two-phase flow, Physics of Fluids **30**, 10.1063/1.5016965 (2018).

[37] H. Huang, M. Krafczyk, and X. Lu, Forcing term in single-phase and Shan-Chen-type multiphase lattice Boltzmann models, Physical Review E **84**, 046710 (2011).

[38] Q. Li and K. H. Luo, Achieving tunable surface tension in the pseudopotential lattice Boltzmann modeling of multiphase flows, Physical Review E **88**, 053307 (2013).

[39] A. Xu, T. Zhao, L. An, and L. Shi, A three-dimensional pseudo-potential-based lattice boltzmann model for multiphase flows with large density ratio and variable surface tension, International Journal of Heat and Fluid Flow **56**, 261 (2015).

[40] Y. Wu, N. Gui, X. Yang, J. Tu, and S. Jiang, Fourth-order analysis of force terms in multiphase pseudopotential lattice boltzmann model, Computers & Mathematics with Applications **76**, 1699 (2018).

[41] L. Fei, K. H. Luo, and Q. Li, Three-dimensional cascaded lattice Boltzmann method: Improved implementation and consistent forcing scheme, Phys. Rev. E **97**, 053309 (2018).

[42] K. H. Luo, L. Fei, and G. Wang, A unified lattice boltzmann model and application to multiphase flows, Philosophical Transactions of the Royal Society A **379**, 20200397 (2021).

[43] J. D. Van der Waals, *Over de Continuiteit van den Gas-en Vloeistoftoestand*, Vol. 1 (Sijthoff, 1873).

[44] D. M. Anderson, G. B. McFadden, and A. A. Wheeler, Diffuse-interface methods in fluid mechanics, Annual review of fluid mechanics **30**, 139 (1998).
[45] D. J. Korteweg, On the form taken by the equations of fluid motion if we take into account the capillary forces caused by considerable but known variations in density and on the theory of capillarity under the assumption of a continuous variation in density, Dutch Archives of Natural Sciences **6**, 1 (1901).

[46] T. Lee and P. F. Fischer, Eliminating parasitic currents in the lattice boltzmann equation method for nonideal gases, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics **74**, 046709 (2006).

[47] Y. Rocard and V. Thermodynamique, Théorie cinétique (1967).

[48] M. Sbragaglia, H. Chen, X. Shan, and S. Succi, Continuum free-energy formulation for a class of lattice Boltzmann multiphase models, EPL (Europhysics Letters) **86**, 24005 (2009).

[49] P. Lallemand and L.-s. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Physical Review E **61**, 6546 (2000), arXiv:0803.2826.

[50] F. Bösch, S. S. Chikatamarla, and I. V. Karlin, Entropic multirelaxation lattice Boltzmann models for turbulent flows, Physical Review E **92**, 043309 (2015).

[51] AL. Kupershtokh, New method of incorporating a body force term into the lattice Boltzmann equation, in *Proceeding of the 5th International EHD Workshop* (Poitiers, France, 2004) pp. 241–246.

[52] Z. Guo, C. Zheng, and B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E **65**, 046308 (2002).

[53] M. E. McCracken and J. Abraham, Multiple-relaxation-time lattice-boltzmann model for multiphase flow, Phys. Rev. E **71**, 036701 (2005).

[54] R. Noël, *The Lattice Boltzmann Method for Numerical Simulation of Continuum Medium Aiming Image-Based Diagnostics*, Ph.D. thesis, Université de Lyon, Saint-Étienne, France (2019).

[55] G. Wang, L. Fei, and K. H. Luo, Unified lattice boltzmann method with improved schemes for multi-phase flow simulation: Application to droplet dynamics under realistic conditions, Physical Review E **105**, 10.1103/PhysRevE.105.045314 (2022).

[56] A. Mazloomi M, S. S. Chikatamarla, and I. V. Karlin, Entropic Lattice Boltzmann Method for Multi-phase Flows, Phys. Rev. Lett. **114**, 174502 (2015).

[57] B. Wen, L. Zhao, W. Qiu, Y. Ye, and X. Shan, Chemical-potential multiphase lattice boltzmann method with superlarge density ratios, Physical Review E **102**, 013303 (2020).

[58] A. Montessori, L. A. Hegele, and M. Lauricella, A thread-safe lattice boltzmann model for multicomponent turbulent jet simulations, AIAA Journal **63**, 1005 (2025).

[59] Horace Lamb, *Hydrodynamics*, 6th ed., 0 (Cambridge University Press, 1993) pp. 1–23.

[60] L. Chen, Q. Kang, Y. Mu, Y.-L. He, and W.-Q. Tao, A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications, Int. J. Heat Mass Transfer **76**, 210 (2014).

[61] P. Yue, C. Zhou, and J. J. Feng, Spontaneous shrinkage of drops and mass conservation in phase-field simulations, Journal of Computational Physics **223**, 1 (2007).

[62] L. Zheng, T. Lee, Z. Guo, and D. Rumschitzki, Shrinkage of bubbles and drops in the lattice boltzmann equation method for nonideal gases, Physical Review E **89**, 033302 (2014).

[63] C. Zhang and Z. Guo, Spontaneous shrinkage of droplet on a wetting surface in the phase-field model,
Physical Review E **100**, 061302 (2019).

[64] L. Fei, F. Qin, G. Wang, K. H. Luo, D. Derome, and J. Carmeliet, Droplet evaporation in finite-size systems: Theoretical analysis and mesoscopic modeling, Physical Review E **105**, 025101 (2022).

[65] Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet, and S. Zaleski, The atomising pulsed jet, Journal of Fluid Mechanics **1009**, A35 (2025).

[66] C. Josserand and S. Zaleski, Droplet splashing on a thin liquid film, Physics of Fluids **15**, 1650 (2003).

[67] F. Bösch, B. Dorschner, and I. Karlin, Entropic multi-relaxation free-energy lattice Boltzmann model for two-phase flows, EPL (Europhysics Letters) **122**, 14002 (2018).