## 1 Introduction  
### Background  
The Lattice Boltzmann Method (LBM) has proven to be an efficient and reliable method for Computational Fluid Dynamics (CFD) for several decades now. Historically, it follows the seminal work of [PERSON] [1] on the use of Lattice Gas Cellular Automata (LGCA) for simulating the Navier-Stokes equations. This work was a remarkable breakthrough, as it proposed simulating fluids at the particle level on an Eulerian grid using interaction rules. Unfortunately, this came at the cost of computation time to overcome the numerical statistical noise due to the Boolean formulation of particle interactions. Decades before, [PERSON], [PERSON] and [PERSON] achieved a linearization of the collision operator of the Boltzmann Equation (BE) [2] that allows for the reduction in the complexity of the equation by avoiding the implementation of an integral collision model. [PERSON] [3] and [PERSON] [4] proposed using discretized versions of the Boltzmann equation to overcome the LGCA' limitations, and gave birth to the so-called Lattice Boltzmann Method. Following this line of thought, [PERSON] et al. proposed a Lattice Boltzmann-BGK Equation (LBGKE) with an optimized local equilibrium and proved the H-theorem for it [5].  
A large amount of work was since done in the LBM field, and recent years saw some remarkable improvements for the numerical simulation of high Reynolds numbersfluids, with for example, the lattice-kinetic theory [6]. One can also mention the Multiple Relaxation Times (MRT) LBM which also allows better stability of the method for high Reynolds numbers, using a multiple-relaxation time term [7]. Such matrix collision operator was first introduced by [8]. All these new insights are relevant for simulating high Reynolds number fluids, but also to open up on the idea of modifying the collision operator. Indeed, the collision operator is key in the method as it embeds the majority of the physics to be simulated.  
The simulation of multi-component flows is another important research interest in the LBM field. This represents one of the major interests of the method: the LBM is conceptually relatively straightforward as it describes particle physics on a mesoscopic scale. However, the main challenge to overcome is the interface definition which can be sorted roughly into four categories. The first category, called color-gradient method [9], has been developed for two-phase flows. It consists in adding a perturbation to the linearized collision operator to make the pressure tensor locally anisotropic near a fluid-fluid interface. This results in surface tension at interfaces while retaining the compatibility with the Navier-Stokes Equations (NSE) in non-interface regions. The second category is based on the introduction of an inter-particle potential [10]. It allows the simulation of multi-phase and multi-component immiscible fluids with different masses at constant temperature, with a high efficiency. The third category can simulate hydrodynamics of phase separation and two-phase flow [11]. The principle is to use a non-ideal pressure tensor in the collision operator to ensure thermodynamic consistency. The goal of this approach is to improve physical consistency. The fourth category takes into account molecular interactions [12; 13] for the simulation of incompressible two-phase flows. A review of multi-phase LBM is proposed in [14]. It gives a comprehensive overview of the methods and the associated algorithmic aspects.  
The idea of combining specifically adapted collision operators and the multi-phase approach leads to the question of the possibility of simulating solid matter with the LBM. In fact, the multi-phase LBM is of great interest, as it may simplify the notion of contact. As an example, the work of [PERSON] et al. [15] applied the multi-phase method to the simulation of ligament break-up and gives insight into its potential in bio-mechanical soft tissue simulations.