An explanation that we postulated for the origin of the mismatch between MD simulations and TvNS solutions was the non-commutability of the order of taking limits: the solution of the Euler equations after taking the scaling limit (as in TvNS solution) is not the same as applying the scaling limit to the solution of the Navier Stokes equations that include heat conduction and viscosity.
In particular, heat conduction terms impose the boundary condition of \(\nabla T=0\)  at \(r=0\) , \(T\)  being the temperature, which the TvNS solution can be seen to not obey [1]}. Recently, it was shown [2]}, [3]} that the results from direct numerical simulations (DNS) of the Navier Stokes equations, which include heat conduction and viscosity terms, compared very well with the results from MD simulations of shock propagation in a one dimensional gas consisting of two kinds of point particles. These studies also argue that the heat conduction term becomes important at distances very close to the shock center.
