In simulations of non-equilibrium gases and plasmas, direct discretization of the Boltzmann equation for the particle distribution function, f(x,v,t), in its full six-dimensional phase space is often impractical. Instead, particle-based kinetic methods such as direct simulation Monte Carlo (DSMC) are commonly used, which sample f with N marker particles that are “pushed” under forces at each timestep and undergo stochastic collisions [1]. However, capturing non-equilibrium dynamics with particle-based methods face two major drawbacks: Monte Carlo estimates of moments converge only as O(1/N^{1/2}), corrupting higher-order statistics, and resolving short collisional scales requires prohibitively large particle counts. Two variance-reduction techniques have been proposed to mitigate statistical noise in particle-based discretizations. \delta f methods evolve only the deviation from a backgrounds equilibrium distribution, f^{eq}, thereby concentrating sampling effort on perturbations of interest [2]. The variance-reduced DSMC (VRDSMC) algorithm uses an importance-sampling formulation to evolve the ratio of f^{eq}/f at particle locations directly [3,4]. While these methods have shown success near equilibrium, they encounter difficulties in strongly non-equilibrium regimes. \delta f methods can produce non-physical particle weights when deviations become large, and VRDSMC's reliance on a fixed equilibrium state saturates variance-reduction once the true distribution departs significantly from that reference. In this talk, we propose to use a variance reduction framework that relies on a multifidelity representation of the full particle distribution function. For instance, such an approach might blend the high‐fidelity DSMC approximation of f with a computationally cheaper surrogate model, f^{NR}, coming from a 5-moment Euler fluid model. Moments can then be computed using a variance reduction formulation where f^{NR} acts akin to a control variate. By relying on multiple representations of the full particle distribution function, we hope to push beyond the difficulties of existing methods and enable accurate, low‐noise simulation of non‐equilibrium dynamics at particle counts suitable for engineering applications