Kinetic–Diffusion Monte Carlo (KDMC) Method

• KDMC is a Monte Carlo transport method that blends explicit kinetic particle tracking with a diffusion approximation for unresolved collisions, automatically adjusting to varying collisionality.
• It employs a hybrid particle trajectory that transitions from kinetic behavior until the first collision to a diffusive step within the same time frame, preserving kinetic moments and asymptotic accuracy.
• KDMC offers significant computational gains by reducing explicit collision counts and integrating with multilevel Monte Carlo techniques and fluid estimators for applications like fusion and radiotherapy.

Kinetic–Diffusion Monte Carlo (KDMC) denotes a class of Monte Carlo transport methods that combine explicitly kinetic particle motion with a diffusion-type representation of unresolved collisional dynamics. In its most specific and technically developed sense, KDMC refers to the asymptotic-preserving hybrid particle algorithm introduced for the Boltzmann–BGK equation in diffusive scaling, in which particles evolve kinetically until a first collision inside a fixed time step and diffusively thereafter, so that the method reverts to a standard velocity-jump process in low-collisional regimes and collapses to a random walk in high-collisional regimes (Mortier et al., 2020). Closely related literature extends this program to neutral transport in fusion reactors, multidimensional implementations, boundary-aware diffusive steps, source-term estimation, and particle-level fluid–kinetic decomposition (Steel et al., 4 Sep 2025, Lappi et al., 23 Sep 2025, Tang et al., 29 Dec 2025, Tang et al., 22 Jun 2026).

1. Terminology and scope

The most precise current usage identifies KDMC with the hybrid asymptotic-preserving Monte Carlo method for Boltzmann–BGK transport in the diffusion limit (Mortier et al., 2020). In that lineage, KDMC is a single particle algorithm valid across low- and high-collisional regimes, with no explicit domain decomposition and no explicit switch between kinetic and diffusive solvers; the method “automatically ‘adjusts’ via the collision frequency” (Mortier et al., 2020).

The designation also appears in neighboring literatures for event-driven Monte Carlo descriptions of transport processes whose macroscopic manifestation is diffusive. This includes proton transport governed by reorientation and hopping events in anhydrous proton conductors (Masłowski et al., 2014), self-learning kinetic Monte Carlo for surface self-diffusion of Ag islands on Ag(111) (Shah et al., 2015), first-passage kinetic Monte Carlo for reaction–drift–diffusion processes (Mauro et al., 2013), and efficient protective-domain Monte Carlo for diffusion with spatially varying annihilation rates (Schwarz et al., 2012). This suggests a broader family resemblance centered on stochastic transport, rare-event simulation, and diffusion as an emergent or asymptotic description.

Usage	Representative characterization	arXiv
Boltzmann–BGK AP hybrid	first collision kinetic, remainder diffusive	(Mortier et al., 2020)
Fusion neutral transport	hybrid fluid–kinetic Monte Carlo for neutrals	(Steel et al., 4 Sep 2025)
2D neutral transport	first 2D extension implemented in Eiron	(Lappi et al., 23 Sep 2025)
Source-term estimation	KDMC combined with associated fluid estimation	(Tang et al., 29 Dec 2025)
Particle-level decomposition	fluid–kinetic decomposition constructed at the particle level	(Tang et al., 22 Jun 2026)
Broader KMC-diffusion usage	proton, surface, and reaction–diffusion transport	(Masłowski et al., 2014)

A recurrent misconception is that KDMC is merely a spatial domain decomposition between kinetic and fluid regions. In the Mortier formulation, the hybridization occurs within each particle trajectory and within each fixed time step, rather than by assigning different spatial regions to different solvers (Mortier et al., 2020).

2. Boltzmann–BGK formulation and the diffusion limit

The foundational KDMC paper considers the one-dimensional Boltzmann–BGK equation in diffusive scaling,

$$\frac{\partial f(x, v,t)}{\partial t} + \frac{v}{\varepsilon}\frac{\partial f(x, v,t)}{\partial x} = \sigma(x)\frac{1}{\varepsilon^2} \left(M(v;x)\int f(x, v,t)\,dv - f(x, v,t)\right),$$

with $f(x,v,0)=Q(x,v)$ (Mortier et al., 2020). Here $\varepsilon$ encodes the diffusive scaling, $\sigma(x)$ is the collision frequency, and $M(v;x)$ is the local Maxwellian equilibrium. The low-collisional regime corresponds to moderate or large $\varepsilon$, whereas the high-collisional regime arises as $\varepsilon\to 0$, when the collision rate $\sigma/\varepsilon^2$ becomes large (Mortier et al., 2020).

A Hilbert-type expansion,

$f(x,v,t)=M(v;x)\rho(x,t)+\varepsilon g(x,v,t),$

yields the limiting macroscopic advection–diffusion equation

$$\frac{\partial \rho(x,t)}{\partial t} + \frac{\partial}{\partial x}\big(u(x)\rho(x,t)\big) - \frac{\partial}{\partial x}\left(\frac{T(x)}{\sigma(x)}\frac{\partial \rho(x,t)}{\partial x}\right)=0$$

in the $f(x,v,0)=Q(x,v)$0 limit (Mortier et al., 2020). In the spatially homogeneous case, this corresponds to the SDE

$f(x,v,0)=Q(x,v)$1

with diffusion coefficient $f(x,v,0)=Q(x,v)$2 (Mortier et al., 2020).

This kinetic-to-diffusive structure reappears in later neutral-transport formulations. For fusion applications, a linear kinetic equation with source, ionization, and charge exchange,

$f(x,v,0)=Q(x,v)$3

is treated by the same broad logic: kinetic transport in low-collisional regimes, and fluid or SDE approximations when charge exchange is large (Steel et al., 4 Sep 2025). A 2D BGK model with homogeneous $f(x,v,0)=Q(x,v)$4,

$f(x,v,0)=Q(x,v)$5

was used for the first two-dimensional KDMC implementation in Eiron (Lappi et al., 23 Sep 2025).

3. Hybrid particle construction

The defining KDMC construction is a hybrid particle that behaves kinetically until its first collision within a fixed time step $f(x,v,0)=Q(x,v)$6, and diffusively for the remaining part of that same step (Mortier et al., 2020). The kinetic segment is the standard velocity-jump process: positions and velocities evolve with free flight between collisions, collision times are sampled from an exponential law in homogeneous media or from

$f(x,v,0)=Q(x,v)$7

in the inhomogeneous case, and post-collisional velocities are sampled from the local Maxwellian (Mortier et al., 2020).

Within one global step of length $f(x,v,0)=Q(x,v)$8, the algorithm proceeds as follows. It starts from the current particle velocity, simulates the first kinetic collision time, advances ballistically until either that collision or the end of the time step, and, if the collision occurs before the end of the step, fills the remaining time with a single Gaussian increment whose mean and variance are chosen from exact kinetic moments conditioned on the final velocity (Mortier et al., 2020). The final velocity is then reused as the initial velocity for the next time step, which preserves inter-step correlation (Mortier et al., 2020).

A central design requirement is that, under spatially homogeneous conditions during the time step, KDMC maintains the correct mean, variance, and inter-step correlation of particle displacements for all $f(x,v,0)=Q(x,v)$9 (Mortier et al., 2020). The diffusive substep is written as

$\varepsilon$0

with $\varepsilon$1-dependent coefficients $\varepsilon$2 and $\varepsilon$3 chosen from exact kinetic moments, rather than from the asymptotic coefficients alone (Mortier et al., 2020).

The absence of an explicit regime switch is fundamental. When $\varepsilon$4, typically no collision occurs during $\varepsilon$5, so the algorithm is essentially purely kinetic. When $\varepsilon$6, the first collision occurs very early, the kinetic phase becomes negligible, and the method resembles a pure diffusion random walk (Mortier et al., 2020).

Later variants preserve this structure while adapting it to specific applications. The KDKMC variant used for fusion boundaries employs a kinetic segment, a diffusive segment, and then a small additional kinetic free-flight step after the diffusive step (Steel et al., 4 Sep 2025). The 2D implementation in Eiron generalizes the diffusive increment to a Gaussian with mean vector $\varepsilon$7 and covariance matrix $\varepsilon$8 in $\varepsilon$9 (Lappi et al., 23 Sep 2025). In electron beam therapy, the “Kinetic-Diffusion-Rotation” algorithm adds a rotation stage and a multiple scattering distribution; the mean kinetic motion is derived analytically, while the variance is supplied through a lookup table (Willems et al., 2024).

4. Accuracy, asymptotic preservation, and computational scaling

In the foundational BGK setting, the scheme is asymptotic-preserving if two limits hold simultaneously: for fixed $\sigma(x)$0 and $\sigma(x)$1, the method converges to the exact kinetic dynamics; and for fixed $\sigma(x)$2 and $\sigma(x)$3, it converges without time-step refinement to a consistent Monte Carlo discretization of the limiting diffusion equation (Mortier et al., 2020). The practical significance is that $\sigma(x)$4 is not forced to resolve the microscopic mean free time $\sigma(x)$5 (Mortier et al., 2020).

The low-collisional error analysis is carried out in 1-Wasserstein distance. With $\sigma(x)$6 the exact kinetic evolution operator and $\sigma(x)$7 the KDMC operator, the local one-step error satisfies

$\sigma(x)$8

under a spatially homogeneous background, which yields convergence to standard kinetic Monte Carlo as $\sigma(x)$9 (Mortier et al., 2020).

In the high-collisional regime, the local difference between the kinetic displacement distribution and the Gaussian KDMC step is estimated through an Edgeworth expansion. The paper proves

$M(v;x)$0

which vanishes for fixed $M(v;x)$1 as $M(v;x)$2 (Mortier et al., 2020). Summed over $M(v;x)$3 steps, this gives

$M(v;x)$4

establishing the asymptotic-preserving property (Mortier et al., 2020).

The computational advantage follows directly from the collision count. In standard kinetic MC, the expected number of collisions per unit time is $M(v;x)$5, whereas in KDMC at most one collision is explicitly simulated per time step per particle. The probability of at least one collision in a time step is

$M(v;x)$6

and the ratio of expected collision counts between standard MC and KDMC scales like

$M(v;x)$7

which grows like $M(v;x)$8 at high collisionality (Mortier et al., 2020).

The multilevel extension ML-KDMC embeds KDMC in a Multilevel Monte Carlo hierarchy and introduces “a new and improved recipe for correlating particle trajectories with different time step sizes,” together with “a new and more general level selection strategy” (Mortier et al., 2020). Earlier multilevel asymptotic-preserving particle work on the Goldstein–Taylor model supplied the general principle that coarse diffusive simulations can act as control variates for fine kinetic simulations in the diffusion limit (Løvbak et al., 2019). In the high-collisional case, ML-KDMC “outperforms the single-level KDMC method by several orders of magnitude” (Mortier et al., 2020).

5. Boundaries, estimators, and hybrid fluid–kinetic closures

Boundary treatment is a major technical issue because the diffusive step in its original form is derived in an unbounded or periodic setting. In the fusion-neutral literature, the original practical response was to revert to a fully kinetic treatment near boundaries whenever the proposed diffusive trajectory crossed the wall or came sufficiently close to it (Steel et al., 4 Sep 2025). This preserves accuracy but can erase much of the computational gain.

A later extension addresses this by deriving the exact Fokker–Planck transition density on $M(v;x)$9 with Robin or absorbing boundary conditions and sampling the diffusive endpoint directly from that boundary-aware density (Steel et al., 4 Sep 2025). In one-dimensional reflecting-boundary experiments, this “can lead to a speedup of up to 500 times compared to a KDMC method that switches to a purely kinetic method, while not sacrificing too much accuracy” (Steel et al., 4 Sep 2025). A plausible implication is that boundary-aware diffusive sampling is indispensable if KDMC is to remain efficient in edge-plasma geometries where walls are dynamically important.

A second major theme is estimator design. In neutral transport, the quantities of interest are often time-integrated source terms or macroscopic moments rather than terminal particle positions. The associated “fluid estimation” procedure combines KDMC trajectories with a fluid post-processing step for the diffusive segments (Tang et al., 29 Dec 2025). In the analyzed 1D setting, the relevant moments are

$\varepsilon$0

and the algorithm combines a kinetic estimator with a fluid estimator derived from the diffusion limit (Tang et al., 29 Dec 2025). In a fusion-relevant test case, the method “consistently achieves lower error than the fluid-based method, and even one order of magnitude lower” while also exhibiting “a significant speedup compared to the reference kinetic MC method” (Tang et al., 29 Dec 2025). A complementary numerical study reports that KDMC with fluid estimation achieves “at least one order of magnitude lower errors than the fluid method for both high- and low-collisional regimes” (Tang et al., 15 Sep 2025).

The 2026 particle-level decomposition makes this estimator logic structural. It constructs a “distribution-decomposition hybrid model constructed at the particle level based on the kinetic-diffusion Monte Carlo (KDMC) method,” with a fluid–kinetic split induced by the kinetic and diffusive stages of each particle trajectory (Tang et al., 22 Jun 2026). To improve fluid-part estimation, the authors derive a “Navier-Stokes-type fluid system” via “Hilbert-Chapman-Enskog expansions, tailored for KDMC,” and introduce a “tunable reflective boundary condition” that balances accuracy and efficiency (Tang et al., 22 Jun 2026). In the reported 1D tests, the method exhibits “at least 500 times speedup over the kinetic MC, while maintaining relative L2 errors around 10% in a charge exchange (CX)-dominant test case” (Tang et al., 22 Jun 2026). The same paper states that accuracy in non-CX-dominant regimes becomes increasingly sensitive to boundary treatment, because of “the inherent limitations of the fluid approximation near the boundary” (Tang et al., 22 Jun 2026).

6. Applications and adjacent lineages

Neutral transport in fusion reactors is presently the most developed application domain. KDMC has been described as “a powerful tool for simulating neutral particles in fusion reactors” (Steel et al., 4 Sep 2025), and the 2D implementation in Eiron identifies neutral particle simulations as “one of the major computational bottlenecks in tokamak scrape-off layer simulations” (Lappi et al., 23 Sep 2025). In the highest-collisional tests of that 2D implementation, KDMC is “between 2 and 3 orders of magnitude faster than fully kinetic simulation” (Lappi et al., 23 Sep 2025).

Radiotherapy provides a distinct but structurally related adaptation. The “Kinetic-Diffusion-Rotation Algorithm for Dose Estimation in Radiation Therapy” states that the underlying kinetic–diffusion particle tracing scheme was “first proposed in the context of neutral transport in fusion energy,” but is here modified by analytic mean matching, a lookup-table variance, and a multiple-scattering distribution (Willems et al., 2024). On a 2D CT scan of a lung patient, the Python implementation is “nearly 33 times faster than an equivalent kinetic simulation at the cost of a small modeling error” (Willems et al., 2024).

Outside the AP Boltzmann–BGK lineage, related event-driven Monte Carlo methods illuminate the wider conceptual neighborhood. The proton-conductivity study models the Grotthuss mechanism through reorientation events with rate $\varepsilon$1 and hopping events with rate $\varepsilon$2, with the ratio $\varepsilon$3 controlling conductivity regimes (Masłowski et al., 2014). The self-learning KMC study of Ag islands on Ag(111) treats diffusion as thermally activated hopping with rates

$\varepsilon$4

and uses center-of-mass mean-square displacement to extract diffusion coefficients and effective barriers (Shah et al., 2015). First-passage approaches such as DL-FPKMC for reaction–drift–diffusion decompose motion into protective domains and sample first-passage events or continuous-time random walks on dynamically generated lattices (Mauro et al., 2013). Protective-domain algorithms for spatially varying annihilation rates use an upper-bound rate $\varepsilon$5 and rejection sampling to generate exact annihilation or boundary-event statistics without small diffusive hops (Schwarz et al., 2012).

These neighboring methods are not identical to the Boltzmann–BGK KDMC scheme. A plausible interpretation is that they delineate a broader methodological zone in which kinetic Monte Carlo is used to realize transport processes whose macroscopic description is diffusive, while the Mortier-type KDMC supplies the most explicit asymptotic-preserving synthesis of kinetic and diffusion dynamics within a single particle algorithm (Mortier et al., 2020).