## 1 Introduction
In order to address this disadvantage, a new flux reconstruction scheme is proposed in this work using particle distribution functions both before and after the collision. A local switch function related to the order parameter is also applied to adjust these particle distribution functions obtained before and after the collision. With this switch function, numerical diffusion can be controlled and the fluxes can be evaluated according to the characteristics of the local multiphase flows. Besides the flow field, fluid interfaces are solved with the DI method, which also needs improvement. In general, the mass conservation in the DI method cannot be exactly guaranteed, and much effort has been invested in solving this problem. [PERSON] et al. [30] found that a droplet in the shear flow would shrink gradually. [PERSON] et al. [31] analyzed a two-phase model and found that the static bubble or droplet showed some spontaneous shrinkage phenomena and could even disappear eventually. Some researchers [32; 33] opted for an improved mesh to overcome the issue. [PERSON] et al. [33] adopted a dual-resolution Cartesian grid to simulate multiphase problems, solving the flow field on a lower resolution grid while solving the phase field on a higher resolution grid. However, choosing the high-resolution grid often increases the computation cost, and the mass non-conservation issue still exists. Researchers have proposed improved models and methods [34; 35; 36]. [PERSON] et al. [34] proposed a method with the mass correction procedure based on the He-Chen-Zhang lattice Boltzmann equation. [PERSON] et al. [35] solved the mass conservation issue by introducing the mass-conserved term into the Cahn-Hilliard (C-H) equation and, on this basis, [PERSON] et al. [36] proposed a mass-conserved fractional step axisymmetric lattice Boltzmann flux solver. However, capturing interface inaccuracy in some local places has not been fully resolved yet. Specifically, the thickness of the fluid interface cannot be maintained due to convection and diffusion, so the calculation of surface tension force is not accurate and numerical stability is also affected. In order to overcome this, many methods [37; 38; 39] have been proposed. [PERSON] et al. [37] proposed an improved interfacial profile-corrected C-H model. [PERSON] et al. [38] proposed a phase-field method that involves correcting the flux. [PERSON] et al. [39] proposed two improved phase-field method formulations based on the profile-corrected and flux-corrected guidelines. [PERSON] et al. [40] combined the level set method with the DI method by introducing re-initialization to reconstruct the phase interface. [PERSON] et al. [41] proposed an interface-compressed method applying the gradient of the normal velocity of the local interface. Currently, methods for improving the interfacial condition are still being developed. In this study, we aimed to develop an interface-corrected DI method with an improved MLBFS for simulating incompressible multiphase flows with large density ratios. Two important improvements were made: one for the DI, allowing it to maintain the equilibrium shape and the thickness of the fluid interface; and the other for the MLBFS, allowing it to evaluate the fluxes according to the characteristics of the local multiphase flows. To enforce the shape and thickness of the fluid interface during the dynamic evolution of the flow field, an interface correction term is introduced into the DI model by applying the Hamilton-Jacobi equation. This term only takes into account effects around the fluid interfaces between different phases. The fluxes are reconstructed by applying the particle distribution functions both before and after collisions. A local switch function is also introduced to automatically adjust the contributions of these distribution functions to the fluxes. Moreover, a mass correction term is also introduced into the DI model.