avel a rich phenomenology, including out of equilibrium phase transitions between homogeneous ($M_{QSS}=0$) and non-homogeneous ($M_{QSS} \ne 0$) QSS states. Second and first order transition lines are found that separate homogeneous and non homogeneous states and merge into a tricritical point approximately located in $(M_{0},e)=(0.2,0.61)$. When the transition is second order two extrema of the Lynden-Bell entropy are identified in the inhomogeneous phase: the solution $M=0$ corresponds to a saddle point, being therefore unstable; the global maximum is instead associated to $M \neq 0$, which represents the equilibrium predicted by the theory. This argument is important for what will be discussed in the following.

Let us now turn to direct simulations, with the aim of testing the above scenario, and focus first on the kinetic model (\[eq:VlasovHMF\])–(\[eq:pot\_magn\]). The algorithm solves the Vlasov equation in phase space and uses the so-called “splitting scheme", a widely adopted strategy in numerical fluid dynamics. Such a scheme, pioneered by Cheng and Knorr [@Cheng], was first applied to the study of the Vlasov-Poisson equations in the electrostatic limit and then employed for a wide spectrum of problems [@califano]. For different values of the pair $(M_{0},e)$, which sets the widths of the initial water-bag profile, we performed a direct integration of the Vlasov system (\[eq:VlasovHMF\])–(\[eq:pot\_magn\]). After a transient, magnetization is shown to eventually attain a constant value, which corresponds to the QSS value observed in the HMF, discrete, framework. The asymptotic magnetizations are hence recorded when varying the initial condition. Results (stars) are reported in figure \[fig1\](a) where $M_{QSS}$ is plotted as function of $e$. A comparison is drawn with the predictions of our theory (solid line) and with the outcome of N-body simulation (squares) based on the Hamiltonian (\[eq:ham\]), finding an excellent agreement. This observation enables us to conclude that (i) the Vlasov equati