## 1 Introduction  
As an important class of flow problems involving two or more fluids, multiphase flows exist widely in nature, as well as in industrial applications, such as droplets hitting a fluid or solid surface, bubble deformation in liquids, and gas-liquid flows in pipes in oil exploitation. In order to simulate multiphase flows, many methods, including the volume of fluid (VOF) method [1, 2], the level set method [3, 4], the front tracking method [5, 6, 7], and the diffuse interface (DI) method [8, 9], have been proposed. The VOF method is able to capture the phase interface accurately [1, 2]. The level set method applies a signed distance function to indicate different phases [3, 4]. The front tracking method has the advantage of accurately capturing phase interfaces with small or moderate changes in a very simple way [6, 7]. However, maintaining mass conservation is still challenging in multiphase flow simulations. DI methods have been widely applied in recent years due to their firm theoretical basis in the phase field method and advantages in dealing with the complicated evolution of the interface among different fluids, and they have been successfully applied to the Rayleigh-Taylor instability in binary fluids [10], bubbles rising or merging in liquids [11], droplet deformation on impact [12], and so on. Like other numerical methods, the DI methods have two important aspects for simulating multiphase flows: the solution of the flow field and the tracking of the fluid interfaces. For the solution of the flow field, current methods are usually classified into two types; i.e., the continuum method and the mesoscopic method, based on the gas kinetic theory. The continuum method is carried out by solving the Navier-Stokes (N-S) equations with the finite difference method, finite volume method, and/or finite element method. The continuum methods include those initially developed by [PERSON] et al. [13] and [PERSON] [14]. Different from the continuum method, the mesoscopic kinetic method, which is also called the lattice Boltzmann method (LBM), is realized by solving the lattice Boltzmann equations. The LBM originated from lattice gas automata [15], and now it has developed into a very useful approach to solve complex multiphase flow problems. The LBM also has the advantage of being able to compute multiphase flow with complicated geometric shapes. Many LB models [16; 17; 18] have been proposed, such as the color model [19]; the interparticle-potential model [20]; the thermodynamically consistent model [21]; the double-distribution-function model (DDF), which includes pressure and density [22]; the projection-like pseudo-potential model [23]; the stable-discretization DDF model [24]; the three-distribution-function model, which includes phase, velocity, and temperature [25]; and the improved single-relaxation-time model [26]. Recently, multiphase lattice Boltzmann flux solvers (MLBFSs) [10; 27; 28; 29] have also been developed and applied to many complex multiphase flow problems, such as the adaptive mesh refinement multiphase lattice Boltzmann flux solver (AMR-MLBFS) developed by [PERSON] et al. [27] and the improved multiphase lattice Boltzmann flux solver (IMLBFS) developed by [PERSON] et al. [28]. In the MLBFS, the flow properties are directly updated at each cell center and the fluxes are calculated locally at cell surfaces through local reconstruction of lattice Boltzmann solutions. The fluxes are obtained through the operations of equilibrium and non-equilibrium particle distribution functions conforming to the lattice Boltzmann equations. This means that the particle distribution functions are obtained after the streaming and collision process and numerical diffusion can be brought into the MLBFS in the process of flux computation. However, numerical diffusion is not enough, and this affects the computation of the MLBFS, resulting in a larger value for the characteristic velocity for simulating multiphase flows at large density ratios. The efficiency of the computation is thus limited.