## 1 Introduction  
Multi-component flows are common in many industrial and technological processes. Extraction and transport of hydrocarbons, liquid jet machining, coating materials, fuel sprays, pollutants treatment, droplets dispersion, and porous media are just a few notable examples. Over the years, many numerical strategies have been developed to deal with such flow problems. The first examples were given by the Lagrangian [1] and arbitrary Lagrangian-Eulerian methods [2]. The accuracies of these methods are very high but limited by the level of interface deformation, which may lead to severe mesh distortions. To reduce such distortions, front-tracking methods [3] were developed. They consider a fixed computational grid for the bulk flow and track the interface using Lagrangian tracers. Despite being less sensible to interface deformations, they become geometrically challenging as the deformation enlarges. To circumvent the aforementioned issues, other strategies, based on the Eulerian description of the flow field, were devised. In particular, the level-set method [4] and the volume of fluid method [5] are among the most successful ones. The former exploits the advection of a suitable function whose zero-level set denotes the interface. The latter tracks the free surface using a volume fraction variable and a geometrical reconstruction technique, e.g., the piecewise linear interface calculation [6]. All methods mentioned above treat the interface as a \(d\)-dimensional surface, with \(d\) being the problem dimension. A different interface-capturing approach is provided by diffuse interface methods [7; 8; 9; 10], which allow the different fluid components to develop a thin mixing region, usually described as a smooth (but possibly sharp) variation of a field variable. Since no reconstruction or tracking of the interface is required, complex three-dimensional interface topologies with possible coalescence/break-up of droplets can easily be captured with no extra cost. In this context, the variable density incompressible (VDI) flow model [11] exhibits interesting features for the simulation of incompressible flow problems involving immiscible fluids of different densities. The interface between components is managed in a diffuse fashion as smooth variations of the density field. Due to the coexistence of the divergence constraint and the mass conservation, the density is treated as a purely advected property. Thus, it acts as a field variable, allowing it to distinguish different components of different densities. As is the case for any diffuse interface method, the VDI flow model needs careful control of interface spreading to avoid meaningless descriptions of the flow problem under investigation. Numerical methods able to deliver high-order spatial accuracy are very attractive in this regard as they account for steep gradients of the field variable, thereby allowing to represent a sharp diffuse interface. Among the many high-order methods proposed in the literature, we focused on discontinuous Galerkin (DG) formulations, which are particularly well suited in the computational fluid dynamics context [12; 13; 14; 15; 16; 17]. Indeed, they provide favorable dispersion and dissipation properties and the ability to achieve very high accuracy even on arbitrarily shaped polytopic mesh elements [18]. In this work, we investigated some numerical challenges involved in implicit high-order accurate DG discretizations of the VDI model introduced in [19] for the simulation of multi-component incompressible flow problems. This discretization follows the idea of coupling the artificial compressibility method to the VDI flow model [19; 20; 21; 22; 23; 24] in order to introduce suitable Godunov fluxes at the inter-element boundaries. In particular, the formulation relies on the exact solution of local Riemann problems based on an artificial compressibility perturbation of the incompressible Navier-Stokes equations.