## 3 Numerical Framework
### Set of Working Variables  
We consider the following set of working variables  
\[\mathbf{w}=\{p,u_{i},\widetilde{\rho}\}, \tag{10}\]  
which comprises the pressure \(p\), the velocity components \(u_{i=1,\ldots,d}\) and the _working density_\(\widetilde{\rho}\), such that \(\rho=f(\widetilde{\rho})\). Hence, in two space dimensions (\(d=2\)) the change of the variable matrix reads  
\[\mathbf{P}(\mathbf{w})=\frac{\partial\mathbf{q}}{\partial\mathbf{w}}= \begin{bmatrix}1&0&0&0\\ 0&\rho&0&\rho_{\widetilde{\rho}}u_{1}\\ 0&0&\rho&\rho_{\widetilde{\rho}}u_{2}\\ 0&0&0&\rho_{\widetilde{\rho}}\end{bmatrix} ,\]  
where \(\rho_{\widetilde{\rho}}=\mathrm{d}\rho/\mathrm{d}\widetilde{\rho}\). In this work, we consider three alternatives for the definition of \(\widetilde{\rho}\), exploiting the change of variable to enforce particular properties of the density field at the discrete level:  
* linear \[\rho=f_{1}(\widetilde{\rho})=\widetilde{\rho};\] (11)
* exponential \[\rho=f_{2}(\widetilde{\rho})=\mathrm{e}^{\rho_{\widetilde{\rho}}\widetilde{ \rho}};\] (12)
* hyperbolic tangent-based \[\rho=f_{3}(\widetilde{\rho})=\frac{r_{\rho}+1+\big{(}r_{\rho}+1\big{)}\tanh \big{(}\sigma_{\mathrm{s}}\widetilde{\rho}\big{)}}{2r_{\rho}}.\] (13)  
In Equation (13), \(r_{\rho}=\rho_{+}/\rho_{-}\) is the maximum density ratio occurring in the density field, with \(\rho_{+}\) and \(\rho_{-}\) the maximum and minimum density values for the flow problem at hand, respectively. In Equations (12) and (13), \(\sigma_{\mathrm{s}}\in\mathbb{R}\backslash\{0\}\) is a scaling factor defined as follows  
Figure 1: Structure of the Riemann problem solution. \[\sigma_{\rm s}=\max\left[\frac{1}{\Delta_{\widetilde{p}}}\ln\left(r_{p}\right),1 \right],\]  
where \(\Delta_{\widetilde{p}}=\widetilde{\rho}(\rho_{+})-\widetilde{\rho}(\rho_{-})\) is a user-defined range. In this work, we consider \(\Delta_{\widetilde{p}}=1\). The max function is introduced to gracefully deal with unit density ratios. Notice that in (11)-(13) the density has to be indented dimensionless with respect to the reference density \(\rho_{0}=\rho_{+}\). According to the linear transformation (11) the set of working variables falls back to the set of primitive variables, hence no particular properties are enforced at the discrete level. Instead, the exponential (12) and the hyperbolic tangent (13) based transformations restrict the range of admissible density values as follows  
\[f_{2}(\widetilde{\rho}):\mathbb{R}\to(0,\,+\infty)\qquad\text{and}\qquad f_{ 3}(\widetilde{p}):\mathbb{R}\to\left(0,\,1+\frac{1}{r_{p}}\right).\]  
thereby ensuring density positivity at the discrete level. Interestingly, by controlling the maximum density value, the hyperbolic tangent transformation in (13) also provides the ability to limit density over/under-shoots associated with high-order approximations of steep density gradients in multi-component flow problems. Nevertheless, it must be clear that Gibbs-type phenomena promoting such over/under-shoots are not eliminated. As detailed in Section 3.4, proper treatment is required to avoid simulation blow up. Moreover, when employing the hyperbolic tangent transformation, in order to avoid failure of the inverse transformation \(\widetilde{\rho}=f_{3}^{-1}(\rho)\), it must be ensured that \(\rho\in\left(0,\,1+1/r_{p}\right)\). Since the inverse transformation is applied to the density obtained from the Riemann solver in order to retrieve the working density to be plugged into the convective numerical flux, it is mandatory to rely on the SDRS. Indeed, as detailed in Section 3.2, only SDRS ensures that the star region density falls within the external states limits. Note that external density values computed from \(\widetilde{\rho}\) are within the limits by construction. Note that ERS is well suited to be used in combination with \(f_{1}(\widetilde{\rho})\) and \(f_{2}(\widetilde{\rho})\) transformations.