---
title: Kinetic–Diffusion Monte Carlo (KDMC) Method
url: https://www.emergentmind.com/topics/kinetic-diffusion-monte-carlo-kdmc
type: topic
---

# Kinetic–Diffusion Monte Carlo (KDMC) Method

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 [2012.08985]. 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 [2509.03942], [2509.19140], [2512.23580], [2606.23368].

## 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 [2012.08985]. 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” [2012.08985].

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 [1408.6165], self-learning kinetic Monte Carlo for surface self-diffusion of Ag islands on Ag(111) [1508.00260], first-passage kinetic Monte Carlo for reaction–drift–diffusion processes [1302.0793], and efficient protective-domain Monte Carlo for diffusion with spatially varying annihilation rates [1206.2203]. 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 | [2012.08985] |
| Fusion neutral transport | hybrid fluid–kinetic Monte Carlo for neutrals | [2509.03942] |
| 2D neutral transport | first 2D extension implemented in Eiron | [2509.19140] |
| Source-term estimation | KDMC combined with associated fluid estimation | [2512.23580] |
| Particle-level decomposition | fluid–kinetic decomposition constructed at the particle level | [2606.23368] |
| Broader KMC-diffusion usage | proton, surface, and reaction–diffusion transport | [1408.6165] |

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 [2012.08985].

## 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)\) [2012.08985]. 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 [2012.08985].

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 \(\varepsilon\to 0\) limit [2012.08985]. In the spatially homogeneous case, this corresponds to the SDE
\[
dX = u\,dt + \sqrt{\frac{2T}{\sigma}}\,dW,
\]
with diffusion coefficient \(D=T/\sigma\) [2012.08985].

This kinetic-to-diffusive structure reappears in later neutral-transport formulations. For fusion applications, a linear kinetic equation with source, ionization, and charge exchange,
\[
\partial_t f + v \cdot \nabla_x f = S - R_i f + R_{cx}\Big(M(v|x,t)\int f\,dv' - f\Big),
\]
is treated by the same broad logic: kinetic transport in low-collisional regimes, and fluid or SDE approximations when charge exchange is large [2509.03942]. A 2D BGK model with homogeneous \(R_{\mathrm{cx}}\),
\[
\partial_t f(x,v,t) + v \cdot \nabla_x f(x,v,t) = R_\mathrm{cx}\left(\rho(x,t)\mathcal{M}(v)-f(x,v,t)\right),
\]
was used for the first two-dimensional KDMC implementation in Eiron [2509.19140].

## 3. Hybrid particle construction

The defining KDMC construction is a hybrid particle that behaves kinetically until its first collision within a fixed time step \(\Delta t\), and diffusively for the remaining part of that same step [2012.08985]. 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
\[
\int_0^{\Delta T_k} \frac{\sigma\big(x(\mathcal T_k) + V_k s/\varepsilon\big)}{\varepsilon^2}\,ds = E_k,\qquad E_k\sim \mathrm{Exp}(1),
\]
in the inhomogeneous case, and post-collisional velocities are sampled from the local Maxwellian [2012.08985].

Within one global step of length \(\Delta t\), 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 [2012.08985]. The final velocity is then reused as the initial velocity for the next time step, which preserves inter-step correlation [2012.08985].

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 \(\varepsilon\) [2012.08985]. The diffusive substep is written as
\[
X(t+\theta)=X(t)+A_\varepsilon(\theta)\,\theta+\sqrt{2D_\varepsilon(\theta)\,\theta}\,\xi,\qquad \xi\sim\mathcal N(0,1),
\]
with \(\varepsilon\)-dependent coefficients \(A_\varepsilon\) and \(D_\varepsilon\) chosen from exact kinetic moments, rather than from the asymptotic coefficients alone [2012.08985].

The absence of an explicit regime switch is fundamental. When \(\sigma\Delta t/\varepsilon^2\ll 1\), typically no collision occurs during \(\Delta t\), so the algorithm is essentially purely kinetic. When \(\sigma\Delta t/\varepsilon^2\gg 1\), the first collision occurs very early, the kinetic phase becomes negligible, and the method resembles a pure diffusion random walk [2012.08985].

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 [2509.03942]. The 2D implementation in Eiron generalizes the diffusive increment to a Gaussian with mean vector \(\mu\) and covariance matrix \(\Sigma\) in \(\mathbb R^2\) [2509.19140]. 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 [2412.05063].

## 4. Accuracy, asymptotic preservation, and computational scaling

In the foundational BGK setting, the scheme is asymptotic-preserving if two limits hold simultaneously: for fixed \(\varepsilon>0\) and \(\Delta t\to 0\), the method converges to the exact kinetic dynamics; and for fixed \(\Delta t\) and \(\varepsilon\to 0\), it converges without time-step refinement to a consistent Monte Carlo discretization of the limiting diffusion equation [2012.08985]. The practical significance is that \(\Delta t\) is not forced to resolve the microscopic mean free time \(\varepsilon^2/\sigma\) [2012.08985].

The low-collisional error analysis is carried out in 1-Wasserstein distance. With \(S_{\Delta t}\) the exact kinetic evolution operator and \(\hat S_{\Delta t}\) the KDMC operator, the local one-step error satisfies
\[
W_1(S_{\Delta t}(f)-\hat S_{\Delta t}(f)) \le 0.24959\, T \sigma^3 \Delta t^5 / \varepsilon^8 + O(\Delta t^7),
\]
under a spatially homogeneous background, which yields convergence to standard kinetic Monte Carlo as \(\Delta t\to 0\) [2012.08985].

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
\[
W_1(S'_{\Delta t}(f)-\tilde S_{\Delta t}(f)) \le 0.58 \sqrt{T}\,\varepsilon^3 \sqrt{\sigma^3 \Delta t},
\]
which vanishes for fixed \(\Delta t\) as \(\varepsilon\to 0\) [2012.08985]. Summed over \(N\) steps, this gives
\[
W_1(f(t_N),f^{(N)}) \le 0.58 \sqrt{T}\,\varepsilon^3 \sqrt{\sigma^3 \Delta t^3},
\]
establishing the asymptotic-preserving property [2012.08985].

The computational advantage follows directly from the collision count. In standard kinetic MC, the expected number of collisions per unit time is \(\sigma/\varepsilon^2\), 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
\[
\mathbb P(\text{≥1 collision in } \Delta t)=1-e^{-\sigma\Delta t/\varepsilon^2},
\]
and the ratio of expected collision counts between standard MC and KDMC scales like
\[
\frac{\sigma \Delta t/\varepsilon^2}{1-e^{-\sigma\Delta t/\varepsilon^2}},
\]
which grows like \(\sigma \Delta t/\varepsilon^2\) at high collisionality [2012.08985].

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” [2004.04071]. 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 [1902.04347]. In the high-collisional case, ML-KDMC “outperforms the single-level KDMC method by several orders of magnitude” [2004.04071].

## 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 [2509.03942]. 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 \((-\infty,L]\) with Robin or absorbing boundary conditions and sampling the diffusive endpoint directly from that boundary-aware density [2509.03942]. 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” [2509.03942]. 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 [2512.23580]. In the analyzed 1D setting, the relevant moments are
\[
m_0(x)=\int_0^{\bar t}\int f(x,v,t)\,dv\,dt,\qquad
m_1(x)=\int_0^{\bar t}\int v f(x,v,t)\,dv\,dt,\qquad
m_2(x)=\int_0^{\bar t}\int \frac{v^2}{2} f(x,v,t)\,dv\,dt,
\]
and the algorithm combines a kinetic estimator with a fluid estimator derived from the diffusion limit [2512.23580]. 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” [2512.23580]. 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” [2509.11883].

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 [2606.23368]. 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 [2606.23368]. 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” [2606.23368]. 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” [2606.23368].

## 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” [2509.03942], and the 2D implementation in Eiron identifies neutral particle simulations as “one of the major computational bottlenecks in tokamak scrape-off layer simulations” [2509.19140]. In the highest-collisional tests of that 2D implementation, KDMC is “between 2 and 3 orders of magnitude faster than fully kinetic simulation” [2509.19140].

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 [2412.05063]. 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” [2412.05063].

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 \(\gamma_R\) and hopping events with rate \(\gamma_T\), with the ratio \(\gamma_T/\gamma_R\) controlling conductivity regimes [1408.6165]. The self-learning KMC study of Ag islands on Ag(111) treats diffusion as thermally activated hopping with rates
\[
k_i = \nu_i \exp\!\left(-\frac{E_i}{k_B T}\right),
\]
and uses center-of-mass mean-square displacement to extract diffusion coefficients and effective barriers [1508.00260]. 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 [1302.0793]. Protective-domain algorithms for spatially varying annihilation rates use an upper-bound rate \(k_m(t)\) and rejection sampling to generate exact annihilation or boundary-event statistics without small diffusive hops [1206.2203].

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 [2012.08985].

Source: https://www.emergentmind.com/topics/kinetic-diffusion-monte-carlo-kdmc