Computing statistical solutions of a Mach 2000 astrophysical jet

The simulation of extreme Mach astrophysical flows is traditionally viewed through the lens of deterministic positivity-preserving schemes. However, due to Kelvin--Helmholtz instabilities and shock anomalies, the multi-dimensional Euler equations admit a variety of non-unique entropy solutions in turbulent regimes. Here, we computationally explore the limits of weak-strong uniqueness of a Mach 2000 jet by defining the statistical solution as the pushforward of a probability measure through a vectorial lattice Boltzmann method operator. Utilizing optimized CUDA kernels, we compute an ensemble of 1000 Monte Carlo samples across a sequence of highly refined spatial grids of up to 3.2 million cells and subsequently post-process the empirical measures via memory-mapped CPU streaming. We contrast the strong sample-wise $L^1$ error divergence with the convergence of the probability measure in the 1-point Wasserstein distance via empirical Cauchy rates. Our results demonstrate that while individual flow realizations physically diverge due to chaotic shear-layer instabilities, the statistical solution converges to an admissible limit measure at a rate of 0.5. Consequently, we provide numerical evidence that the statistical solution to the considered problem is non-Dirac and remains stable in the extreme compressible regime.