Generalized equilibria for color-gradient lattice Boltzmann model based on higher-order Hermite polynomials: A simplified implementation with central moments

We propose generalized equilibria of a three-dimensional color-gradient lattice Boltzmann model for two-component two-phase flows using higher-order Hermite polynomials. Although the resulting equilibrium distribution function, which includes a sixth-order term on the velocity, is computationally cumbersome, its equilibrium central moments (CMs) are velocity-independent and have a simplified form. Numerical experiments show that our approach, as in Wen et al. [Phys. Rev. E 100, 023301 (2019)] who consider terms up to third order, improves the Galilean invariance compared to that of the conventional approach. Dynamic problems can be solved with high accuracy at a density ratio of 10; however, the accuracy is still limited to a density ratio of $1\,000$. For lower density ratios, the generalized equilibria benefit from the CM-based multiple-relaxation-time model, especially at very high Reynolds numbers, significantly improving the numerical stability.


I. INTRODUCTION
Multiphase and multicomponent flows are ubiquitous phenomena observed in both industry and nature.The complexity and diversity of these flows render them a fascinating subject for numerical modeling, simulation, and theoretical research.The lattice Boltzmann (LB) method, initially proposed by McNamara and Zanetti [1] in 1988, has attracted attention as a powerful computational fluid dynamics tool for capturing multicomponent and multiphase flows [2][3][4].The multiphase LB model can be classified into the following categories based on the physical content of the fluid-fluid interface algorithm.
Among the abovementioned models, the CG model, which introduces virtually colored distribution functions, offers many advantages in simulating multiphase and multicomponent flows.First, the model features strict mass-conservation properties for each phase and high flexibility in setting the interfacial tension [20].In this model, the static drop test is not required to determine the interfacial tension, which can be obtained directly without any analysis or assumptions.Second, the CG model is known to have superior dissolution properties when simulating small droplets or bubbles compared to other multiphase LB models [15]; i.e., the small droplets or bubbles are less prone to disappear.In addition, the surface tension, density ratio, and viscosity ratio can be selected independently [21,22].Studies have been conducted comparing different multiphase LB models [23][24][25].Datadien et al. [25] quantitatively compared the CG and pseudopotential models in terms of accuracy and stability and found that the CG model exhibited better characteristics.Because of these outstanding features, the CG model has been applied to various problems related to multiphase flow, such as pseudo-boiling [26], water transport in membranes [27,28], surfactant transport [29,30], thermocapillary flows [31], and droplets on microstructure [32], etc.Because of the compatibility of the LB method with parallel computing, the CG model has also been used for GPU-accelerated computation [33] and to acquire teacher data for machine learning [22].In this study, we focus on the CG model.
The CG model originates from the two-component lattice-gas automata (LGA) of Rothman and Keller [34].Gunstensen et al. [5] proposed the first CG model by combining the single-phase LB method of McNamara and Zanetti [1] with the two-component LGA developed by Rothman and Keller.A perturbation step is introduced to recover Laplace's law at the interface by adding a binary fluid collision operator to the post-collision state at sites near the interface.This model is formulated in two dimensions, and a three-dimensional model was published shortly thereafter [35].Subsequently, Grunau et al. [6] introduced different densities and viscosities by incorporating the freedom of the rest particle equilibrium distribution.This model has been successfully applied to the two-dimensional Rayleigh-Taylor instability [36], which is known as dynamic interfacial instability.Reis and Phillips [21] extended the CG model to a twodimensional nine-velocity (D2Q9) lattice.They modi-fied the perturbation operator to correctly recover the Navier-Stokes equations and showed that it could simulate the dynamic process of merging two droplets at a density ratio of 18.5.Liu et al. [37] derived a generalized perturbation operator using the phase field (or order parameter) instead of a color gradient and described the CG model in three dimensions.The applicability of the CG model to dynamic problems, such as a droplet in a shear flow and a single bubble rising in a viscous fluid, was demonstrated.Leclaire et al. [38] generalized the CG model in two and three dimensions to include solid wall boundary conditions and implemented it using the open-source PALABOS library [39].Mora et al. [40] highlighted the importance of isotropy in the calculation of color gradients to accurately capture behaviors such as pore-scale phenomena.Recently, Subhedar [41] employed a velocity-based equilibrium function, initially proposed for the phase-field LB model [42], to handle high-density ratios within the CG LB framework.
A characteristic feature of the CG model is the recoloring operation, which plays an important role in maintaining immiscibility at the interface and mimics the separation mechanism.The original Gunstensen algorithm [5] is implemented by solving a maximization problem for the work performed by the color gradient against the color flux.To reduce velocity fluctuations, Tölke et al. [43] proposed a modified algorithm in which the phase separation was not as strong as that of the Gunstensen algorithm, but the numerical stability was enhanced.Instead of widening the interface width, Latva-Kokko-Rothman's recoloring algorithm [44], designed based on the work of D'Ortona et al. [45], addresses certain issues with the previous CG model, namely the lattice-pinning problem and spurious currents near the fluid-fluid interface.Subsequently, Halliday et al. [46] improved this algorithm slightly.Subhedar et al. [47] showed that Halliday's algorithm has smaller spurious currents than that of Latva-Kokko and Rothman [44].Leclaire et al. [48] integrated the Latva-Kokko-Rothmann algorithm into the Reis-Phillips CG model and demonstrated that integrating Latva-Kokko-Rothman's recoloring operator into Reis-Phillips' perturbation operator greatly improves the numerical stability and accuracy of solutions over a wide range of parameters.Using a higher-order isotropic gradient operator also enhances numerical stability and accuracy [49].Recently, the mathematical similarity between the recovered macroscopic equation from the recoloring algorithm and the conservative Allen-Cahn equation [50], which is an interface-capturing equation used in phase-field modeling, has been discussed in several studies [41,47,51].
The original CG model suffers from a lack of Galilean invariance because of the non-Navier-Stokes terms identified by Liu et al. [37].To restore Galilean invariance, correction terms should be added to the equilibrium distribution function [52,53].Note that a similar problem exists in the free-energy model [54][55][56].Following the analysis of Holdych et al. [57], a source term for improv-ing Galilean invariance was derived by Leclaire et al. [53] and incorporated into an equilibrium distribution function.This enhanced equilibrium distribution function improved the momentum discontinuity problem through numerical tests on a layered Couette flow and was then successfully adopted in many studies [58][59][60].Using a different approach, Ba et al. [20] modified an equilibrium distribution function based on the third-order Hermite expansion of the Maxwellian distribution, taking a cue from Li et al. [61].They also demonstrated that this modification improved the discontinuous velocity.The CG model in Ref. [20] was developed in two dimensions, whereas Wen et al. [62] later extended it to three dimensions using a D3Q19 lattice.
The collision operation in the LB method plays a pivotal role in ensuring the numerical stability and accuracy of computations [63], irrespective of whether the flow is single phase or multiphase.Coreixas et al. [64] have conducted an exhaustive review of collision models, presenting a systematic organization and defining the mathematical relationships among them.From the viewpoint of computational cost, all moment spaces should be approximately equivalent since the LB scheme is memorybound.This has been shown notably by Bauer et al. [65] on CPUs and by Latt et al. [66] on GPUs.The lattice Bhatnagar-Gross-Krook (BGK) equation [67,68] stands out as the simplest and most widely used collision operator in the LB method.Based on a single relaxation time approximation, the BGK equation [69] within the LB framework can be perceived as the propensity of the distribution function to gravitate toward its equilibrium state after a specified relaxation time [4].However, despite the success of the lattice BGK equation, it tends to become numerically unstable under the flow conditions of high Reynolds number (low-viscosity) [70].Several collision models have been proposed to counteract these shortcomings [63,64].Specifically, d'Humières [71] aimed to enhance the numerical stability of the LB method by introducing a collision step in the moment space coupled with an increase in the free parameters.This technique, which assigns different relaxation times to distinct moments, is called the multi-relaxation-time (MRT) collision model.The MRT model has garnered considerable attention owing to its enhanced numerical stability, particularly after the publication of its formulations in two-dimensional (D2Q9) [72] and three-dimensional (D3Q15 and D3Q19) spaces [73].In particular, determining higher-order relaxation parameters is notably challenging for three-dimensional configurations.As a solution, the two-relaxation-time [74] and entropic MRT models [75,76] have been proposed to offer closures for these elusive relaxation parameters.
Although the MRT model offers improved numerical stability compared to the lattice BGK model, it encounters instability at elevated Reynolds numbers.To address this issue and further bolster numerical stability, Geier et al. [77] introduced the cascaded LB method.This method executes the collision step in the moment space corresponding to the comoving reference frame.Notably, in contrast to that in the raw-moment-based MRT (hereafter referred to as RM-MRT) model mentioned earlier [71][72][73], collision operations are performed in the space of the central moments (CMs), and this is referred to as CM-MRT.Results of several numerical simulations using the CM-MRT model have shown that it has excellent numerical stability in simulations corresponding to high Reynolds numbers [77][78][79].Indeed, a critical point of CM-MRT is the relaxation of (1) high-order moments and (2) the moment related to bulk viscosity.By imposing a relaxation frequency of 1, the stability is increased through a higher hyperviscosity and a higher bulk viscosity [80,81].It is noted that using the same relaxation time for all moments, whether they be RMor CM-based, leads to the same stability domain [80].The relationship between the RM-and CM-MRT models can be described under the general MRT framework in Ref. [82,83].The recently proposed work of Luo et al. [84] aims at a unified framework, which seamlessly integrates the widely used existing collision models.However, it would only be a first step toward a unified framework for comparing different collision models.For example, the matrix approach used in their work cannot reproduce the cumulant collision model [85] since it is not possible to express the cumulant collision model in a linear matrix form, as pointed out in Ref. [64].
Because of its excellent numerical stability, CM-MRT has been successfully applied to multiphase flows [59,[86][87][88][89][90][91].Lycett-Brown and Luo [86] pioneered the pseudopotential model-based formulation.In the context of the CG LB model, Leclaire et al. [87] first introduced CM-MRT in a formulation with a unit density ratio.Subsequently, Saito et al. [59] developed a CG model with a density ratio using the nonorthogonal CM set proposed by De Rosis [78] and showed that the simulation of dynamic liquid jet flows with extremely high Reynolds numbers (up to O(10 6 )) is possible.To improve the potential Galilean invariance in the CG model, they used the equilibrium distribution function described by Leclaire et al. [53].Nevertheless, as highlighted by Wen et al. [62], this equilibrium distribution function leaves error terms in the recovered macroscopic momentum equation.Additionally, the functional form of the equilibrium CMs that transfers the equilibrium distribution function to the CM space is complex and cumbersome to implement.
In this study, we propose the generalization of an equilibrium distribution function to improve the Galilean invariance in the CG model using higher-order Hermite polynomials.Furthermore, we demonstrate that the equilibrium CMs of the proposed generalized equilibrium distribution function are extremely concise and easily executable.The remainder of this paper is organized as follows: In Sec.II, we outline the underlying D3Q27 CG LB model and describe its implementation within the framework of CM-MRT.In Sec.III, the equilibrium distribution function within the CG LB framework is expressed in its general form using Hermite polynomials, including a comparison with existing expressions.Furthermore, their characteristics and functional forms in the CM space are investigated.In Sec.IV, through several numerical experiments, the numerical properties of the generalized equilibrium distribution function derived in this paper are presented and compared with existing distribution functions.Sec.V concludes this paper.

A. Model description
In the present LB model, the distribution functions f i move on a D3Q27 lattice (i ∈ [0, • • • , 26]), which is a straightforward extension of the D2Q9 model [92], with the lattice velocity c i defined as: where c = δ x /δ t , δ x is the lattice spacing, and δ t is the time step.Hereafter, the formulation is given by δ x = δ t = 1, which is similar to the typical LB method.
In the CG LB model for two-phase flows, the two immiscible fluids are represented by introducing virtual red and blue fluids.Distribution functions f k i represent the fluids k, where k = r and b denote "red" and "blue," respectively.The total distribution function is expressed as follows: The time evolution equation of the distribution function is expressed as the following LB equation with a forcing term [93]: where x and t denote the position and time, respectively.The last term F i introduces the body force into the LB equation; several implementations of the forcing term in the phase space have been proposed, which are reviewed in Ref. [94].The LB equation ( 3) can be split into the collision step, described as: and the streaming step, described as: where f * i denotes the post-collision distribution functions.In the lattice BGK model [67,68], the collision term is obtained as: where τ is the relaxation time and f eq i is the equilibrium distribution function, which is discussed later in this paper.The following equilibrium distribution functions are often used in standard CG LB models [20,21,37,40,48,95,96]: where c s is the speed of sound defined by the LB method for standard single-phase flows (c 2 s = 1/3 for the D3Q27 lattice [92]) and φ i is the lattice-specific parameter, which can be expressed as: For the D3Q27 lattice, ξ = 9/19, and ξ values for other lattice models can be found in Ref. [38].The weight function is given by: The density of the k-phase fluid is given by: The total fluid density is given by ρ = k ρ k .The momentum is defined as follows: where F denotes the body force.Excluding a few recent CG models [41,51], the pressure p and density ρ are connected by an ideal gas equation of state (EOS).To obtain a continuous pressure at the interface, the CG model introduces different sound speeds for each phase.Using the k-phase speed of sound c k s , the pressure can be expressed as: where p k = ρ k (c k s ) 2 is the pressure of k-phase.This feature results in different sound speeds being defined for both phases.From the pressure balance p r = p b in the interface region, we obtain the following relationship between the density ratio and speed of sound [21,51,95,97]: where the superscript "0" indicates the initial density value of the simulation [53] and the relation (c k s ) 2 = ξρ k (1 − α k ) holds.Eq. ( 13) implies that the greater the density ratio, the greater the difference in the speed of sound between the two phases; the speed of sound in the liquid phase becomes progressively smaller than that in the gas phase.Note that to overcome such limitations, Lafarge et al. [51] reformulate the pressure definition under the assumption of mechanical equilibrium [98,99] and successfully introduced a mixture EOS for two-phase flow, inspired by the stiffened gas formulation [100], into the CG model.
The fluid interface is tracked using an order parameter that distinguishes the two components in a multicomponent flow, defined as [20,59]: The order parameter values ϕ = 1, −1, and 0 correspond to a purely red fluid, a purely blue fluid, and the interface between the two, respectively [43].The interfacial tension between the two fluids is introduced as a spatially varying body force F s based on the continuum surface force (CSF) [101], which is defined as: where σ is the interfacial tension coefficient and κ is the interface curvature [102], which is expressed as: where the surface gradient operator ∇ s = (I − nn) • ∇.
Using the order parameters in Eq. ( 14), the unit normal vector is defined as: A recoloring step is introduced to maintain immiscibility in the interfacial region, which is a feature of the CG LB model.Following the algorithm developed by Halliday et al. [46], the recoloring operation can be represented as: with where β is the parameter that controls the interface thickness.In Eq. ( 19), the pressure p comes from the equilibrium at rest as usually found in the literature (e.g.[38,59]).Unless otherwise stated, we set β = 0.7 to reproduce the correct interfacial behavior with as narrow an interface as possible [37,46,103].After summing f r, * i and f b, * i , the streaming operation is performed according to Eq. ( 5) and the boundary conditions are implemented as necessary.
The partial derivatives of the variable χ are evaluated using the second-order isotropic finite difference [37,104,105]: To ensure smoothness across the interface, the kinematic viscosity is interpolated using the following order parameter: where ν k denotes the k-phase kinematic viscosity.From the relationship between the relaxation time and viscosity [20], we obtain: where µ = ρν denotes dynamic viscosity.

B. CM-MRT collision operation
In this study, the collision step in Eq. ( 4) is implemented in the CM space.To achieve this, it is necessary to transform the distribution functions into CMs.This procedure first converts the distribution functions into RMs, which are further converted into CMs.The RMs and CMs are defined as [77]: where α, β, and γ range from zero to two, yielding 3 3 = 27 moments.Although the CMs can be computed directly from Eq. ( 24), fewer computations are required when using the RMs [Eq.( 23)].Specific conversion expressions between the distribution functions, RMs, and CMs can be obtained by executing get f m k relations.jl,as presented in the Supplementary Material.
In this paper, the collision step in the CM space is based on Appendix D of Geier et al. [85]; however, the correction terms added to the second-order moment collisions in Ref. [85] are not considered because they are based on a Taylor expansion that falls apart close to strong discontinuities such as walls and phase interfaces.It also differs from Ref. [85] in that it considers third-and fifth-order forcing terms, in accordance with Refs.[60,83,84].As a result, the collision step can be expressed as: where k * αβγ and k eq αβγ are the post-collision and equilibrium CMs, respectively.The equilibrium CMs are described in Sec.III.ω , ω 1 , • • • , ω 10 are the relaxation rates.ω 1 is related to the kinematic viscosity based on the relation ω 1 = 1/τ and Eq. ( 22), and ω 2 is associated with the bulk viscosity.The others are free and can be chosen from the range {0 • • • 2}.The correction terms Q x , Q y , and Q z appearing in the second-order collision stages are specified in Sec.III E.
Following the collision operation in the CM space, the CMs are transformed back into the equilibrium distribution functions via the RMs.The specific conversion expressions between them can be obtained by executing get f m k relations.jl,as described in the Supplementary Material.

III. EQUILIBRIA FOR CG LB
Before discussing equilibria in the CG model, we briefly review equilibria for single-phase flows.Within this framework, the most widely used equilibrium distribution function can be represented as [13,67]: This equation is obtained by Taylor expanding the Maxwell-Boltzmann distribution with macroscopic velocities assuming a low Mach number and considering terms up to O(u 2 ).In contrast, De Rosis and Luo [106] expanded the Maxwell-Boltzmann distribution to an arbitrary order O(u N ) using Hermite polynomials [107,108].The relaxation of the continuous Maxwellian distribution is equivalent to that of its discrete counterpart when the equilibrium state is constructed using sixth-order Hermite polynomials in three dimensions [106], described as follows: where H iαβγ denotes the Hermite polynomial presented in Appendix B. When Hermite polynomials of orders of three and higher are neglected, Eq. ( 27) can be reduced to Eq. (26).
We now return to the CG LB model.The equilibrium distribution function within the CG framework expressed by Eq. ( 7) is generally adopted when dealing with density contrasts.In this study, we reformulate the equilibrium distribution function for the CG model as: where E i is an isotropic operator that ensures interface isotropy [51], which can be expressed as: The correction operator Φ i is introduced to recover the Galilean invariance.The specific functional form of the correction operator Φ i is defined in this section.The second term of the RHS in Eq. ( 28) represents the deviation from the ideal gas; that is, when p = ρc 2 s for each phase (e.g., unit density ratio), the second term can be neglected and the formulation for the usual single-phase flow can be recovered.
The functional form of equilibrium CMs is also an important aspect in the following discussion.As depicted in Eq. ( 24), the equilibrium CMs can also be computed as follows: Here, we investigate the equilibrium CMs for several equilibrium distribution functions.Scripts with symbolic computations used to specifically compute the equilibrium CMs get equilibria.jlcan be found in the Supplementary Material.
A. Standard CG equilibria with g eq,2 i We begin with the most popular equilibrium distribution function in the CG model, which is expressed in Eq. (7).By setting N = 2 and Φ i = 0 in Eq. ( 28), Eq. ( 7) can be reformulated as: Substituting Eq. ( 31) into Eq.( 30), we compute the corresponding equilibrium CMs in each order as follows: Zeroth order: First order: Second order: Third order: Fourth order: Fifth order: Sixth order: Based on Eqs. ( 32)- (38), several properties of the equilibrium CMs in Eq. ( 31) can be identified.
• The underlined terms: products of (p − ρc 2 s ) and u α .
• The double underlined terms: products of ρ and u α .
• The remaining terms: independent of u α .
B. Standard CG equilibria with g eq,6 i In single-phase flow models, upon applying the correct set of Hermite polynomials to the discrete equilibrium, the velocity dependence in the derived equilibrium CMs vanishes, leading to Galilean invariance [106].To demonstrate the effects of similar procedures on the CG model, we set N = 6 and Φ i = 0 in Eq. ( 28): This procedure corresponds to the replacement of g eq,2 i in Eq. (31) using g eq,6 i [Eq.(27)].The equilibrium distribution function adopted in Ref. [60] is similar to Eq. (39).However, the previous study uses a different form of correction operator, unlike the method used in Refs.[53,59] [see Eq. (D5) in Ref. [60] for details].Substituting Eq. ( 39) into Eq.( 30), the corresponding equilibrium CMs are computed as: Zeroth order: First order: Second order: Third order: Fourth order: Fifth order: Sixth order: No change is observed in the equilibrium from the zeroth to second order; they are identical to those in Eqs. ( 32)- (34) in Sec.III A. This is because the equilibrium distribution functions in Eqs.(31) and (39) are identical from the zeroth to second order.Upon comparing with Eqs. ( 35)- (38), it can be seen that differences appear at the equilibrium CMs for the third and higher orders.
That is, the double-underlined terms in the form of the product of ρ and the velocity components, as seen in Eqs. ( 35)- (38), are now completely removed.Therefore, one of the third-order moments k eq 111 becomes zero, However, there are still underlined velocity-dependent terms proportional to (p − ρc 2 s ), as shown in Eqs. ( 43)- (46).With the standard equilibrium distribution function in Eq. ( 7) [or Eq. ( 31)], the recovered macroscopic equation includes an unwanted error term [37], as argued by Huang et al. [52].To reduce the effect of this error term, Li et al. [20,62] proposed an improved equilibrium distribution function based on their previous research [61].The original form of this improved equilibrium distribution function can be formulated as: where D denotes the spatial dimension.The last term of Eq. ( 47) is derived on the basis of a third-order Hermite expansion of the Maxwell-Boltzmann distribution [61,108].Upon rewriting Eq. ( 47) in a form equivalent to Eq. ( 28), the nonzero correction operator Φ i is obtained as: Here, we replace g eq,2 i in Eq. ( 48) using g eq,6 i , as described in Sec.III B and then substituting it into Eq.( 30) to obtain the equilibrium CMs, which yields the following results: Zeroth order: First order: Second order: Third order: Fourth order: Fifth order: Sixth order: The equilibrium CMs from the zeroth to second order remain unchanged; the higher-order terms above the fourth order change slightly compared to Eqs. ( 44)- (46).However, there still remain underlined velocity-dependent terms proportional to (p − ρc 2 s ).Notably, all of the thirdorder equilibrium moments in Eq. ( 43) become zero.By using the modified equilibrium distribution function in Eq. ( 48), based on the third-order Hermite expansion, the velocity-dependent terms proportional to (p − ρc 2 s ) in the third-order equilibrium CMs can be eliminated; i.e., it has been implied that the N th-order velocity-dependent terms proportional to (p − ρc 2 s ) in the equilibrium CMs can be eliminated by considering the appropriate form of Φ i in Eq. ( 28) with the N th-order Hermite polynomials.

D. Generalized equilibria
Based on the above results and considerations, we propose a new generalized equilibrium distribution function for the CG model.Similar to g eq,N i in Eq. ( 28), we adopt an equilibrium distribution function based on Hermite expansion up to the sixth order g eq,6 i [60,106].In addition, considering the fourth-, fifth-, and sixth-order effects in Φ i in Eq. ( 28), we define the following equilibrium distribution function: Naturally, if we neglect the terms involving the Hermite polynomials of the fourth, fifth, and sixth orders and replace g eq,6 i with g eq,2 i , Eq. ( 56) is reduced to the original equilibria, as proposed by Li et al. [20,62] and described by Eq. ( 47) or (48).
In the phase space, the new equilibrium distribution function proposed in Eq. ( 56) contains terms up to O(u 6 ), increasing its complexity compared to the conventional form in Eq. ( 7), and its implementation is cumbersome.However, computing the equilibrium moments using Eq. ( 30) yields the following results: Zeroth order: First order: Second order: k eq 110 = k eq 011 = k eq 101 = 0, k eq 200 = k eq 020 = k eq 002 = p, Third order: k eq 120 = k eq 102 = k eq 012 = k eq 210 = k eq 201 = k eq 021 = k eq 111 = 0, Fourth order: Fifth order: Sixth order: The equilibrium moments of the generalized equilibrium distribution function [Eq.( 56)] are no longer velocitydependent in the CM space.Nonzero equilibrium moments are limited to even orders, with only six such moments; their form is simple.All the odd moments are zero.For p = ρc 2 s as in the LB method for single-phase flows, this is attributed to the equilibrium CMs, as described in Refs.[77,83,85,106,109].The equilibrium in the CM space is extremely simple; thus, its execution has also advantages in terms of coding simplicity.

E. Correction term and implementation
Considering the third-order moments of the equilibrium distribution function given by Eq. ( 56), we obtain: The third-order moments of the proposed equilibrium distribution function are equivalent to those in Ref. [62] except for the third-order term ρu α u β u γ .As discussed in literature, even with corrections made to the equilibrium distribution function, the diagonal elements of the third-order moments contain the term ρc 2 s because of the low symmetry of the standard lattices.To address this issue, a correction term Q = [Q x , Q y , Q z ] is computed and considered in the collision process, as described below [20,62]: The divergence operation is approximated using Eq. ( 20).The correction term computed using Eq. ( 65) is considered in the collision step as the external force on the relevant second-order diagonal moments (k 200 , k 020 , and k 002 ), as shown in Eq. (25).
In this study, all relaxation coefficients except for ω 1 are set to unity as in the literature (e.g., Ref. [60]).Regarding the bulk viscosity, this setting is considered to enhance the numerical stability of the present simulations by over-dissipating acoustic waves, as discussed in Ref. [64].With such a choice of relaxation coefficients, the collision step in Eq. ( 25) can be reformulated as fol-lows: Here, since our results so far have shown that k eq 100 = k eq 010 = k eq 001 = k eq 110 = k eq 011 = k eq 101 = 0 and k eq 200 = k eq 020 = k eq 002 = p for all types of equilibria, we have substituted these values.When using the generalized equilibria proposed in this study [Eqs.( 60)-( 63)], the collision operations over the third order in Eq. ( 66) can be expressed more concisely: This implies that only six second-order CMs (k 110 , k 011 , k 101 , k 200 , k 020 , and k 002 ) must be computed immediately before the collision operation, and the other CMs need not be computed in this case.This not only simplifies the calculations but is also expected to reduce the computational cost.
The equilibria models considered in this study are summarized in Table I.

IV. NUMERICAL RESULTS AND DISCUSSION
To investigate the numerical properties of the generalized central equilibrium moments obtained in Sec.III, five numerical experiments are performed.Comparisons with existing equilibrium distribution functions are also presented.To ensure a fair comparison, all collision operations are performed in the CM space.All relaxation coefficients are set to 1, except those related to the kinematic viscosity.As wall boundaries in this paper, we implement the no-slip, slip, and moving wall boundary conditions using the fullway bounce-back scheme [4] for simplicity.The computation code was written in the Julia language [110], version 1.6.7,and parallelized using the MPI.jl package [111].

A. Stationary droplet
The first case involves the simulation of a stationary droplet.As a fundamental property of the numerical performance involving fluid-fluid interfaces, the validity of the surface tension obtained using the CG model described is checked, as described in Sec.II.By measuring the pressure difference between the inside and outside of the droplet at equilibrium, the surface tension is expressed, as shown in Eq. (15).It should be emphasized that, in principle, the CG model does not require prior computations to obtain the surface tension coefficient.However, this computation was performed for validation purposes.Simulations were performed using Model D (Table I).Because the simulation was performed in two dimensions, according to Laplace's law, the pressures inside and outside the droplet were theoretically predicted as: The computational setup was adapted from Ref. [20].A circular droplet of red fluid with a radius R = 25 was placed in a 100 × 100 discretized space.The surrounding area was filled with a blue fluid with ρ 0 b = 1.Periodic boundary conditions were imposed at all boundaries.The kinematic viscosities of each phase were set to ν r = ν b = 1/6.Other parameters were also set as in Ba et al. [20] and compared with the previous results.
Table II summarizes the simulation parameters and the error evaluations for certain density ratios and surface tension coefficients.For stationary cases, the simulations are stable for density ratios up to 1000.The errors in the table are calculated as E = |σ th −σ cal |/σ th ×100%, where the subscripts 'th' and 'cal' denote the theoretical and calculated surface tension coefficients, respectively.By comparing the results of the existing RM-MRT model [20] with those of the present CM-MRT model used in this study, we find that the evaluated errors are of the same order of magnitude.Therefore, we can state that our generalized equilibria [Eqs.( 57)-( 63)] can accurately predict the surface tension, which is pivotal for the simulation of two-phase flows.

B. Droplet in a moving tube
Modifications to the standard equilibrium distribution function cause a lack of Galilean invariance in some free-energy [54][55][56] and CG models [62].As reported previously, if the Galilean invariance is lacking, the initially circular droplet in the moving tube is largely deformed.In this section, we investigate the improvement in the equilibrium distribution function and the effect of the correction term described in Sec.III E for the recovery of the Galilean invariance.We investigate the behavior of the four types of equilibria mentioned in Table I.
We employ a computational setup similar to that in Wen et al. [62].The computational domain was discretized as 140 × 140 in a two-dimensional space.Initially, a droplet with radius R = 30 and density ρ 0 r = 3 was placed at rest.The droplet was surrounded by a stationary fluid with density ρ 0 b = 1.The top and bottom boundaries enforce a wall boundary condition moving at a constant velocity [96,112] with U = 0.02, whereas periodic boundary conditions are set for the left and right boundaries.The movements of the top and bottom walls drive the motion of the fluid inside, which in turn moves the droplets.The additional simulation parameters are as follows: Images of the density distribution obtained using the conditions listed in Table I are shown in Fig. 1.As shown in Figs.1(a) and (b), the droplets in both cases transition into elliptical shapes over time.This is because of the lack of Galilean invariance in the equilibria employed in these models.In the CG model, only considering g eq i up to the sixth order [106] does not inherently improve the Galilean invariance.In contrast, Figs.1(c) and (d) show nearly equivalent and better results, respectively, than those in Figs.1(a) and (b).In other words, the droplets maintain their circular shape even after a long period of time.Thus, we can see that Galilean invariance is restored in these models.This is attributed to the contribution of the improvement in the third-order moments.
To clarify whether Φ i or Q contributes to the improvement in Galilean invariance, an additional simulation based on Model D but without Q is also conducted [see Fig. 1(e)].Although the shape and motion of the droplet are different from those in Figs.1(a) and (b), the droplet is deformed over time.Therefore, we conclude that both the correction operator Φ i of third order or higher and the correction term Q are necessary to improve the Galilean invariance of the CG model.
However, looking at again the results in Figs.1(c) and (d), one can see that the droplet is not perfectly circular, e.g., for t = 54 000; therefore, it seems that the Galilean invariance is not perfectly restored.The gradient com-putation of the correction term Q [Eq.( 65)] by finite differences introduces numerical errors that distort the droplet.A higher-order lattice (e.g., D3Q39 lattice) and the corresponding third-order equilibrium [108] may be needed to solve it completely.Interestingly, no significant difference is observed between Figs. 1(c) and (d).This implies that at the Navier-Stokes level, the correction of the error term in the CG model is sufficient up to the third order in terms of accuracy.This was to be expected because moments up to the fourth order are sufficient to recover even the compressible Navier-Stokes-Fourier equations [108,113].
To observe the differences in more detail, the velocity distribution characteristics in Models A and D are illustrated in Fig. 2. In Fig. 2(a), an unphysical velocity distribution is observed around the deformed droplet.A similar conclusion was drawn in a previous study [20].However, as shown in Fig. 2(b), no significant velocity discontinuity exists near the interface, and the velocity distribution is smooth.This improvement in the unphysical flow is expected to result in a significant difference in the computational accuracy, especially for more complex computation targets.

C. Layered two-phase flow
To further investigate the effects of the four types of equilibria presented in Table I on the velocity profiles of the two-phase flow, immiscible layered flows between two parallel plates [3,52] were simulated.In the two-dimensional simulation, periodic boundary conditions were applied to the left and right boundaries, whereas no-slip boundaries were applied to the top and bottom boundaries.A constant body force F = [F x , 0] was applied to the entire domain as the driving force.In this flow problem, the vertical velocity component u y is assumed to be zero throughout the domain.
Assuming a Poiseuille-type flow in the channel, the analytical solution for the velocity profile is given by [3]: where the coefficients are defined as: where M = µ r /µ b denotes the dynamic shear-viscosity ratio.
In our simulations, the computational domain was discretized into a pseudo-one-dimensional setup with N x × N y = 10 × 100 lattices, where a = 25 and b = 50.Three conditions equivalent to those in Ref. [62] are considered, as summarized in Table III, where the maximum density ratio reaches 1 000.The surface tension coefficient and constant body force are set as σ = 0.002/4.5 and F x = 1.5 × 10 −8 , respectively, as in Wen et al. [62].
As described in Sec.IV B, the four equilibria models (Table I) are examined.Figure 3 shows the numerically obtained velocity profiles for each equilibrium presented in Table I together with the analytical solution of Eq. ( 69).Similar to the simulation results described in Sec.IV B, Models A and B and Models C and D show nearly equivalent results.Without correcting for the equilibria (Models A and B), the deviation between the numerical and analytical solutions increases, and a discontinuity near the interface  can be observed in all three cases.In contrast, when a correction term is applied (Models C and D), the velocity profiles obtained from the simulations agree well with the analytical solution.Again, the accuracy at the Navier-Stokes level of Models C and D is comparable, as described in the previous section.

D. Two-dimensional bubble rising
To verify the accuracy of the dynamic gas-liquid twophase system, we performed a well-known single-bubble rising benchmark simulation, based on a previous study Hysing et al. [114].This is an unsteady problem involving parameters such as the density ratio, viscosity ratio, surface tension, interface curvature, and gravity; it is more complex than the simulations described above.Because no exact solution to this problem can be obtained, a comparison with the available numerical solutions is made.The dimensionless numbers describing the problem are the density ratio ρ r /ρ b , viscosity ratio µ r /µ b , Reynolds number Re, and Eotvos number Eo.The Reynolds and Eotvos numbers are defined as follows: where D denotes the bubble diameter (characteristic length scale), U = (gD) 1/2 is the characteristic velocity, and σ is the surface tension coefficient.The time scale is characterized by D/U = (D/g) 1/2 .Here, we performed simulations for the two cases presented in Ref. [114]).
The cases and corresponding parameters in lattice units and dimensionless numbers are summarized in Table IV.
For details on the computational setup and parameters, refer to Ref. [114].
The simulation was performed in a rectangular computational domain with width W = 1 and height H = 2W in dimensionless units.A bubble of diameter D = 0.5 was placed at the center of the lower half of the domain in a domain filled with liquid.This single bubble rose with gravity, with F = −ρg as the driving force.No-slip boundary conditions were applied for the top and bottom edges, whereas slip boundary conditions were applied for the lateral sides.Throughout this simulation, the relation α b = 8/27 is used to set the speed of sound in the gas phase to c b s = 1/ √ 3, which is equivalent to that of the standard LB method using D3Q27 [58].In addition to the bubble shapes obtained from the simulation, according to [115], the centroid position y c and bubble velocity V c are measured using the following definitions: where the blue density ρ b (x) is used to track the bubble regions.In recent years, this benchmark simulation has also been addressed by many researchers in the framework of the LB method to verify the accuracy of interface tracking [90,[116][117][118][119]; it is noted that most of them are based on the phase-field LB model.Leclaire et al. [115] was the first to address this benchmark problem in the CG model; however, in their simulation, only Case 1 was addressed due to the limitation of the density ratio.
Here, simulations are performed using only Model D, the generalized equilibria (Table I).
Figure 4 shows the time evolution of the bubble shape, center of mass, and rise velocity obtained from the simulation results for Case 1.In this case, the domain was discretized into N x ×N y = 160×320, and we set β = 0.7.Using the parameter set in Table IV, the Mach number with respect to the speed of sound in the liquid phase was obtained as M a = U/c r s = 0.00548 and the number of iterations up to T = 3 was calculated to be 364 000.As shown in Fig. 4(a), the initially circular bubble rises under the effect of gravity and is deformed by the balance between gravity and surface tension.Similar to the benchmark solution [114], no bubble breakup occurs.In Fig. 4(b) and (c), the present results are compared with available numerical data: finite-element-method simulations with a level-set-based approach (with FreeLIFE solver) [114] and phase-field-based approach [120] and the phase-field LB model [119].The time histories of the center of mass and rise velocity are in good agreement with those of all existing studies.
Figure 5 presents the simulation results for Case 2. In this case, the domain was discretized into N x × N y = 320 × 640.After several test simulations, we set β = 0.2  III.The four types of equilibria (Table I) are examined.When a correction is applied (Models C and D), the velocity profiles obtained from the simulations agree well with the analytical solution.to widen the range of interfaces.Using the parameter set in Table IV, the Mach number with respect to the speed of sound in the liquid phase was obtained as M a = U/c r s = 0.0548 and the number of iterations up to T = 3 was calculated to be 728 000.From Fig. 5(a), it can be verified that the simulation is stable and does not break at a density ratio of 1 000.At T = 3, the ends of the bubble appear to break, similar to that observed in previous studies [114,119].From Fig. 5(b) and (c), it can be concluded that the time history of the center of mass agrees well with that observed in the existing studies; however, the discrepancy is large with respect to the rising velocity.There are two main reasons for this finding.The first corresponds to the definition of the bubble velocity given in Eq. (74).To manage high-density ratios, the interface was thickened by setting β = 0.2, and the velocity in the interfacial region was included in the evaluation within the definition of Eq. ( 74), which decreased the overall rise velocity.The second factor is the influence of the compressibility of the liquid phase.In a typical CG model, the speed of sound in the liquid phase decreases with an increasing density ratio because of the assumption of an ideal gas EOS for each phase.A reduction in the speed of sound in the liquid phase at such high-density ratios can be prevented by introducing a suitable EOS [51] to replace the ideal gas EOS.Nevertheless, to the best of our knowledge, this presents a novel approach to solving Case 2 using the CG LB model.

E. Three-dimensional Rayleigh-Taylor instability
The final numerical experiment involves the threedimensional Rayleigh-Taylor instability.The Rayleigh-Taylor instability is a fundamental interfacial instability induced when a heavy fluid is placed over a light fluid subjected to a slightly disturbed interface in a gravitational field [121].Under certain conditions, Kelvin-Helmholtz instability occurs owing to the velocity difference across the fluid-fluid interface.This problem involves complex interface deformation.
We refer to the computational setup adopted in a previous study He et al. [122].A schematic of the computational setup is shown in Fig. 8 of Ref. [58].The top and bottom boundaries are no-slip walls and the lateral boundaries are periodic.As described in Ref. [122], a single-mode initial perturbation is imposed as: in the mid-plane, where W is the computational domain width.The body force for this problem is incorporated as: where g = (0, 0, −g).The gravitational acceleration g is chosen to satisfy the relation (W g) 1/2 = 0.04 [122].IV.Results from finiteelement-method simulations with level-set-based approach (with FreeLIFE solver) [114] and phase-field-based approach [120] and the phase-field LB model [119] are also shown.IV.Results from finite-element-method simulations with level-set-based approach (with FreeLIFE solver) [114] and phase-field-based approach [120] and the phase-field LB model [119] are also shown.

The computational domain is set as
which is a dimensionless number that characterizes this problem, is fixed at 0.5 throughout the simulations.This corresponds to a density ratio of 3. Interface tension is neglected; thus, the perturbed interface is expected to always be unstable in the inviscid case.The kinematicviscosity ratio is set to unity.Another dimensionless number is the Reynolds number, which is defined as: In this problem, time is scaled by (W/g) 1/2 .Figure 6 shows the time evolution of the interface, that is, the isosurface with ϕ = 0, obtained from the simulations for Re = 1 024.The upper panels [Fig.6(a)] show the results obtained using the original equilibria (Model A in Table I), whereas the lower panels [Fig.6(b)] show the results obtained using the generalized equilibria (Model D in Table I).Surfaces regarded as interfaces are colored according to the velocity magnitude.In both cases, the overall trend of the interface shape is similar.For example, the spike tip descends over time to form a mushroom-like roll-up.In the initial stage, at T = 1, no significant difference is observed in the change in the shape of the interface because of the equilibrium employed.Focusing on the later stages, after T = 2, differences are observed in both the interface shape and velocity, depending on the equilibria used.This is because of the difference in accuracy owing to the applied equilibria, as discussed in Secs.IV B and IV C; the generalized equilibria proposed in this study accurately capture the velocity distribution near the interface, whereas the original equilibria result in an unphysical velocity discontinuity near the interface.Over time, these differences in behavior have a greater effect on the interface shape and velocity field.
To better quantify the time evolution of the Rayleigh-Taylor instability, the interfacial position was measured, and the results are shown in Fig. 7. Measurements were taken at three characteristic points: the bubble, saddle,  I), while the generalized equilibria proposed in this study are used for the bottom panels (Model D in Table I).
and spike (their locations are shown in Fig. 6).The result of the generalized equilibria is shown as a solid line, whereas the result of the original equilibria is shown as a broken line.The difference between the generalized equilibria (solid line) and original equilibria (broken line) is particularly noticeable in the time history of the spike.
The original equilibria exhibit a slower evolution of the tip position.The velocity discontinuity near the interface is considered to have prevented the development of the spike tip.In addition, the data obtained by the improved CG LB [62] and multiphase LB flux solver [123] are shown in Fig. 7 for comparison with previous numerical simulations.In regards to the saddle point and spike tip, the present generalized equilibria results are in good agreement with those in Refs.[62,123], with no noticeable differences.In regards to the bubble tip, the generalized equilibrium LB and multiphase LB flux solver results are in good agreement, although the improved CG model progresses slightly faster.
The following discussion focuses on whether any difference can be observed between the simulations of Models C and D (Table I).Based on the numerical experiments presented in Secs.IV B and IV C, it can be concluded that their computational accuracies are nearly equivalent.For more details, additional simulations of the Rayleigh-Taylor instability were performed.The same settings as those in Fig. 6 were used, and only the kinematic viscosity was changed to ν = 0.The simulations were performed at an infinite Reynolds number.Figure 8 shows the time history of the spike tip at an infinite Reynolds number.The results for Re = 1 024 are also presented for reference.Images of the interface shape at infinite Reynolds numbers are shown.We observe that the evolution of the spike tip is considerably larger as the Reynolds number is increased.Furthermore, the computation is stable.From the interface shape, it can be seen that the Kelvin-Helmholtz instability at the interface is particularly prominent in the later stages.However, in the case of Model C (red solid line), the computation was broken, even though it was computed within the same CM-MRT framework.This indicates that the velocity-dependent components of the equilibrium CMs in Eqs. ( 53)-(55) contribute to the numerical instability of the computation.The proposed generalized equilibrium CMs [Eqs.( 57)-( 63)] is effective in avoiding this numerical instability.It should be noted that the ability to stably simulate does not necessarily imply an accurate solution, as discussed in De Rosis [78].In simulations involving high Reynolds numbers, choosing appropriate values for both the mesh size and kinematic fluid viscosity is of paramount importance.In conclusion, the equilibria proposed in this study not only provide simplified velocity-independent forms in the CM space but also benefit from the high numerical stability of the CM-MRT and can stably solve multiphase flows with extremely high Reynolds numbers.

V. CONCLUSIONS
In this study, we propose the generalized equilibria for a three-dimensional CG LB model and investigate its numerical properties within the framework of the CM-MRT model.
First, the equilibrium distribution function of the CG model is reformulated in a more general form using Hermite polynomials [Eq.( 28)], which can provide a prospec- tive indication of the deviation from the ideal gas scenario.An examination of the equilibrium distribution function and its CMs in the existing CG model shows that the velocity dependence of the equilibrium CM vanishes in the order of the single-phase equilibrium distribution function g eq,N i and the correction operator Φ i .Inspired by this fact, we formulate an equilibrium distribution function in Hermite polynomials up to the sixth order for both g eq,N i and Φ i .In the phase space, this equilibrium distribution function contains terms up to O(u 6 ), increasing its complexity compared to its original form, and its implementation becomes more cumbersome.However, most equilibrium CMs are zero, and the nonzero moments have simple functional forms with no velocity dependence.
Numerical experiments show that our approach improves the Galilean invariance, and the accuracy for static problems is comparable to that of the previous approach [20,62], considering up to the third order for Φ i .In a dynamic problem, the bubble rise benchmark set by Hysing et al. [114] was applied and was found to be as accurate as that of other methods for a density ratio of 10.
In contrast, at a density ratio of 1 000, the center of mass was comparable to that obtained with other studies, but the rise velocity was different.The numerical stability of this dynamic and high-density-ratio problem is noteworthy; however, further improvements are required to solve it more accurately within the framework of the CG model.
Finally, simulations of the three-dimensional Rayleigh-Taylor instability were performed under the condition of At = 0.5.In the simulations for Re = 1 024, a large difference is observed in the time evolution of the interface shape and tip position with and without the Galilean invariance correction.Furthermore, simulations with zero kinematic viscosity (infinite Reynolds number) show that the equilibrium CMs above the fourth order significantly contribute to the numerical stability.
The concept of the proposed generalized equilibria is applicable to the forcing-based free-energy LB model [56].Because the equilibrium CMs have simplified forms, those in the cumulant space are similarly simple.Therefore, the present CM-MRT-based model can be easily extended to a cumulant-based model [85].Furthermore, by following the strategy of De Rosis and Coreixas [88] starting from equilibrium distribution functions proposed in Ref. [64], the present framework could be implemented on more compact lattices (e.g., D3Q19) to improve the computational efficiency.

C.
Improved equilibria by Li et al.

FIG. 1 .
FIG. 1. Density field of droplet in a moving tube: (a) Model A, (b) Model B, (c) Model C, (d) Model D.Initially (t = 0), a circular droplet is placed in a domain.Owing to the movement of the top and bottom walls, the droplet also begins to move.In Models A and B, the droplets transition to an elliptic shape over time owing to the lack of Galilean invariance.In Models C and D, the droplets remain nearly circular over time owing to improved Galilean invariance.To clarify whether Φi or Q contributes to the improvement in Galilean invariance, (e) an additional simulation based on Model D but without Q, is also conducted.In this case, the droplet still deforms over time.

FIG. 2 .
FIG. 2. Comparison of the velocity field at t = 80 000δt: (a) Model A and (b) Model D. The solid line in the figures represents the interface position.The correction of Galilean invariance significantly improves the discontinuity in the velocity distribution.

FIG. 3 .
FIG. 3. Velocity profiles of the simulated layered two-phase flow in a channel together with the analytical solution: (a) Case A, (b) Case B, and (c) Case C, as presented in TableIII.The four types of equilibria (TableI) are examined.When a correction is applied (Models C and D), the velocity profiles obtained from the simulations agree well with the analytical solution.

FIG. 6 .
FIG. 6. Interface evolution of three-dimensional Rayleigh-Taylor instability for Re = 1024 and At = 0.5.The original equilibria for the CG model are used for the top panels (Model A in TableI), while the generalized equilibria proposed in this study are used for the bottom panels (Model D in TableI).

FIG. 7 .
FIG.7.Time histories of interface position for At = 0.5 and Re = 1024 with generalized (solid line) and original equilibria (broken line).Results from the improved CG RM-MRT model[62] and multiphase LB flux solver[123] are also shown.

FIG. 8 .
FIG. 8. Time history of spike position during Rayleigh-Taylor instability with images of interfaces at infinite Reynolds number with zero viscosity.

TABLE I .
Summary of equilibria models considered in this study.

TABLE III .
[62]meters used in the simulation of two-phase layered flow[62].

TABLE IV .
Physical parameters in lattice units and dimensionless numbers defining the test cases.