Kinetic-Diffusion Monte Carlo

• Kinetic-diffusion Monte Carlo is a hybrid simulation method that integrates kinetic free flight with diffusion updates to capture both low- and high-collisional regimes.
• It employs an event-driven approach where a kinetic free-flight step is followed by a diffusive remainder designed to match the mean, variance, and correlations of the true process.
• KDMC significantly reduces computational costs, provides second-order accuracy, and adapts well to multilevel Monte Carlo frameworks for complex transport simulations.

Searching arXiv for recent and foundational papers on kinetic-diffusion Monte Carlo to ground the article in the relevant literature. Kinetic-diffusion Monte Carlo denotes a family of stochastic simulation methods that couple explicit kinetic motion with diffusion-scale updates, or, more broadly, kinetic Monte Carlo formulations used to study diffusion-dominated phenomena. In the most specific contemporary usage, especially for Boltzmann-BGK neutral-particle transport, it refers to an asymptotic-preserving hybrid particle method in which each particle trajectory contains a kinetic free-flight segment and, when collisionality is high, a diffusive remainder within the same macroscopic time step. In this sense the method is designed to remain accurate in both low-collisional and high-collisional regimes, avoiding the exploding simulation cost of fully kinetic Monte Carlo in the diffusion limit (Mortier et al., 2020, Mortier et al., 2020). The expression is also used more loosely across lattice KMC, reaction-diffusion, and materials-diffusion literatures, where the common theme is event-driven stochastic evolution of diffusive or diffusion-reaction processes (Leetmaa et al., 2014, Schwarz et al., 2012).

1. Terminological scope and historical placement

The phrase “kinetic-diffusion Monte Carlo” is not confined to a single algorithmic lineage. In plasma-edge and kinetic-theory work it usually denotes a hybrid kinetic/diffusive particle tracer for transport equations with BGK-type relaxation. In reaction-diffusion and lattice KMC work, related ideas appear as event-driven diffusion-reaction simulation, first-passage acceleration, residence-time KMC, and diffusion-focused lattice hopping models (Mortier et al., 2020, Mauro et al., 2013).

Usage in the literature	Core model	Representative papers
Hybrid KDMC for kinetic equations	Boltzmann-BGK or linear Boltzmann transport with diffusion limit	(Mortier et al., 2020, Mortier et al., 2020, Tang et al., 29 Dec 2025)
Multilevel kinetic-diffusion Monte Carlo	AP particle schemes embedded in MLMC	(Løvbak et al., 2019, Mortier et al., 2020)
Event-driven reaction-diffusion KMC	Diffusion-annihilation or reaction-drift-diffusion in protective domains	(Schwarz et al., 2012, Mauro et al., 2013)
Lattice diffusion KMC in materials	Lattice hopping, reactions, coalescence, tracer diffusion	(Leetmaa et al., 2014, Rosenthal et al., 2011, Du et al., 2012)

Within the narrow KDMC sense, the key conceptual step is not merely replacing transport by diffusion, but hybridizing both within each particle update. The 2020 Boltzmann-BGK formulation describes “hybridized particles that exhibit both kinetic behaviour and diffusive behaviour depending on the local collisionality,” and proves that the method reverts to a standard velocity-jump process in the low-collisional regime while collapsing to a standard random walk in the high-collisional regime (Mortier et al., 2020). This distinguishes KDMC from purely kinetic Monte Carlo, purely fluid diffusion solvers, and domain-decomposition hybrids that impose a fixed interface between kinetic and fluid regions. A plausible implication is that KDMC occupies an intermediate methodological position: it is particle based like kinetic MC, but asymptotically diffusion consistent like a fluid approximation.

2. Mathematical basis: from kinetic transport to diffusion

The contemporary KDMC literature is grounded in kinetic equations under diffusive scaling. A canonical starting point is the Boltzmann-BGK model

$$\partial_t f + v\cdot \nabla_x f = \frac{1}{\varepsilon}(M[f]-f),$$

or, in diffusive scaling,

$$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$

with $\rho=\int f\,dv$ and $M$ a Maxwellian determined by local moments or prescribed background fields (Mortier et al., 2020, Mortier et al., 2020). In plasma-edge neutral transport, the collision term is often written with a charge-exchange rate $R(x)$ or $R_{\mathrm{cx}}$, and the diffusion limit yields drift-diffusion equations for the density (Tang et al., 29 Dec 2025, Lappi et al., 23 Sep 2025).

The diffusion limit follows from a Hilbert-type expansion. Writing $f=M\rho+\varepsilon g+\cdots$, one obtains relaxation to local equilibrium at leading order and a constitutive correction at next order. In one formulation this gives

$$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$

with diffusion tensor

$D(x)=\int (v-U)\otimes (v-U)M\,dv / R(x),$

while in the one-dimensional zero-mean isothermal case it reduces to

$$\partial_t \rho = \partial_x\bigl(R(x)^{-1}\partial_x \rho\bigr)$$

(Mortier et al., 2020). In the AP particle formulation for kinetic equations in the diffusion limit, the limiting macroscopic equation is likewise derived as

$$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$0

for a BGK collision operator $$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$1 (Løvbak et al., 2019).

This asymptotic structure is the mathematical reason KDMC is possible. When the mean free path is small, explicit simulation of every collision is unnecessary if the particle increment over a fixed step can be replaced by a random variable with the correct limiting drift and variance. The method is therefore “asymptotic-preserving” in the specific sense that the algorithm remains stable and consistent as $$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$2, without requiring a time-step restriction such as $$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$3 that burdens classical particle schemes (Løvbak et al., 2019, Mortier et al., 2020).

3. Core algorithmic structure

The basic KDMC update begins kinetically and becomes diffusive only after a collision occurs inside the current macroscopic step. In the homogeneous Boltzmann-BGK formulation, given a particle with position $$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$4 and velocity $$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$5, one first samples a free-flight time $$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$6. In homogeneous media,

$$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$7

while in nonhomogeneous settings the flight time $$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$8 solves

$$\partial_t f + \frac{1}{\varepsilon} v \partial_x f = \frac{\sigma(x)}{\varepsilon^2}\bigl(M(v;x)\rho(x,t)-f\bigr),$$9

If no collision occurs before $\rho=\int f\,dv$0, the particle free-streams for the full step. If $\rho=\int f\,dv$1, the particle first free-flies to the collision point, then draws a post-collision velocity from the local Maxwellian, and finally executes a single diffusive remainder over $\rho=\int f\,dv$2 (Mortier et al., 2020, Mortier et al., 2020).

The defining feature is the construction of the diffusive remainder. In the 2020 KDMC method, the Gaussian substep is chosen so that the positional increment over one full time step has the correct mean, variance, and inter-step correlation under spatial homogeneity during the time step (Mortier et al., 2020). In the more general plasma-edge formulation, the remainder takes the form

$\rho=\int f\,dv$3

where $\rho=\int f\,dv$4 and $\rho=\int f\,dv$5 are chosen so that the mean and variance match the true conditioned kinetic process over the unresolved remainder (Mortier et al., 2020). In the AP particle scheme for the diffusive limit, a related hybridization appears through the coefficients

$\rho=\int f\,dv$6

and the transport-diffusion update

$\rho=\int f\,dv$7

followed by a collision step with probability

$\rho=\int f\,dv$8

(Løvbak et al., 2019).

Two common misconceptions are addressed directly by this literature. First, KDMC is not a purely diffusive random walk with an occasional collision correction; the update always begins with a genuine kinetic segment, and in low-collisional regimes it reverts to standard kinetic simulation (Mortier et al., 2020, Mortier et al., 2020). Second, asymptotic-preserving does not mean bias free. The AP time-discrete scheme introduces an $\rho=\int f\,dv$9 bias, and later analyses identify an intermediate regime around $M$0 in which KDMC bias is largest (Løvbak et al., 2019, Tang et al., 29 Dec 2025).

4. Error analysis, multilevel acceleration, and computational complexity

A central motivation for KDMC is computational cost. In a standard kinetic BGK Monte Carlo, the expected number of resolved collisions per time step scales like $M$1. In the KD scheme, at most one collision is simulated per $M$2, which yields an asymptotic speed-up of

$M$3

(Mortier et al., 2020). In high-collisional regimes this is the basis for the observed orders-of-magnitude savings reported in later implementations (Mortier et al., 2020).

The analytical results distinguish low-collisional and diffusive limits. For fixed $M$4, the 2020 KDMC analysis proves a local $M$5 error of order $M$6 and a global $M$7 error of order $M$8, so the scheme is second-order accurate in $M$9 in the low-collisional regime (Mortier et al., 2020). In the high-collisional limit, the same work shows that as $R(x)$0, the hybrid step converges to the random-walk discretization

$R(x)$1

with a global Wasserstein error bounded by $R(x)$2, establishing the asymptotic-preserving property (Mortier et al., 2020).

A separate but closely related development places AP kinetic-diffusion particle methods inside a multilevel Monte Carlo hierarchy. For a quantity of interest $R(x)$3, the MLMC estimator uses levels $R(x)$4 with the telescoping decomposition

$R(x)$5

Under the reported variance and cost scalings,

$R(x)$6

one obtains overall cost $R(x)$7 in the case $R(x)$8, and in the representative Goldstein-Taylor experiment at RMSE $R(x)$9, the CPU cost drops from $R_{\mathrm{cx}}$0 for single-level AP Monte Carlo to $R_{\mathrm{cx}}$1 for MLMC, a $R_{\mathrm{cx}}$2 speed-up (Løvbak et al., 2019). In the BGK-specific ML-KDMC formulation, high-collision behavior is especially favorable: $R_{\mathrm{cx}}$3 and $R_{\mathrm{cx}}$4, leading to optimal $R_{\mathrm{cx}}$5 overall work in the diffusive regime and several-orders-of-magnitude speed-ups relative to single-level KDMC (Mortier et al., 2020).

Fusion-oriented analyses extend these results from trajectory simulation to moment estimation. For time-integrated moments $R_{\mathrm{cx}}$6, the 2025 analysis proves one-step and global error bounds, showing that the KDMC simulation bias obeys two asymptotic regimes,

$R_{\mathrm{cx}}$7

with $R_{\mathrm{cx}}$8 the statistical term. The same work reports that the hybrid estimator consistently achieves lower error than a purely fluid-based method, and even one order of magnitude lower in a fusion-relevant test case, while exhibiting speed-ups up to $R_{\mathrm{cx}}$9 compared to fully kinetic MC (Tang et al., 29 Dec 2025).

5. Boundary conditions, multidimensional extensions, and fluid-kinetic decomposition

Boundary treatment is a major technical issue because the diffusion approximation is least reliable near walls and interfaces. A traditional workaround in boundary-aware KDMC is to abandon the fluid step and revert to purely kinetic substepping whenever a diffusive move would cross a physical boundary, but this sacrifices the main computational benefit near the wall (Steel et al., 4 Sep 2025). The 2025 boundary-condition extension replaces that fallback by an exact analytical solution of the one-dimensional Fokker-Planck equation with a Robin boundary,

$f=M\rho+\varepsilon g+\cdots$0

and derives a boundary-aware sampling strategy based on the decomposition $f=M\rho+\varepsilon g+\cdots$1. In the reported reflecting-wall test, this “KDMC_Fluid” treatment yields speed-ups of up to $f=M\rho+\varepsilon g+\cdots$2 times over a KDMC method that switches to a purely kinetic method, while maintaining the characteristic fluid-approximation bias rather than the erratic switching-induced behavior of the older boundary treatment (Steel et al., 4 Sep 2025).

The move to higher dimensions preserves the same hybrid principle but complicates the diffusive increment. In the first two-dimensional implementation, the kinetic equation is

$f=M\rho+\varepsilon g+\cdots$3

and after one kinetic free-stream of length $f=M\rho+\varepsilon g+\cdots$4, the remaining time is covered by a two-dimensional Gaussian increment $f=M\rho+\varepsilon g+\cdots$5, where $f=M\rho+\varepsilon g+\cdots$6 contains both an isotropic part and an anisotropic rank-one part involving $f=M\rho+\varepsilon g+\cdots$7 (Lappi et al., 23 Sep 2025). In high-collisional benchmarks on a $f=M\rho+\varepsilon g+\cdots$8 slab, KMC cost scales with $f=M\rho+\varepsilon g+\cdots$9 while KDMC runtime is essentially constant, producing speed-ups up to $$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$0–$$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$1 (Lappi et al., 23 Sep 2025).

A more ambitious extension is particle-level fluid-kinetic decomposition. Instead of imposing a spatial interface between a fluid solver and a kinetic solver, this approach lets each KDMC trajectory split naturally into kinetic and diffusive segments inside each time step. The diffusive segment is then represented by a fluid reconstruction. In the 2026 formulation, the kinetic model

$$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$2

is coupled to a Navier-Stokes-type fluid system derived via Hilbert-Chapman-Enskog expansions and tailored for KDMC (Tang et al., 22 Jun 2026). The reported one-dimensional tests show at least $$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$3 times speedup over kinetic MC while maintaining relative $$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$4 errors around $$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$5 in a charge-exchange-dominant case. The same work emphasizes a limitation that is also visible in earlier boundary studies: in non-CX-dominant regimes, accuracy 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. Relation to broader diffusion-oriented KMC traditions

Outside the BGK/AP literature, kinetic Monte Carlo has long been a standard tool for diffusion and diffusion-reaction systems, but these methods usually resolve a different mathematical object: a Markov jump process on a lattice, graph, or protective-domain decomposition rather than a transport equation with an explicit diffusion limit. The connection is methodological rather than terminological identity.

In lattice KMC, the classical variable-step-size method or $$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$6-fold way selects events using cumulative rates

$$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$7

followed by the stochastic time increment

$$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$8

KMCLib implements this framework for diffusion and reaction of millions of particles in one, two, or three dimensions, supports site-specific Arrhenius rates through a Python rate-calculator plugin, and maintains per-particle coordinate lists for on-the-fly mean-square-displacement analysis (Leetmaa et al., 2014). In materials science this basic architecture underlies KMC studies of coupled surface and bulk diffusion of metal clusters in polymers (Rosenthal et al., 2011), hydrogen diffusion in idealized bcc-Fe grains with DFT-derived site energies and migration barriers (Du et al., 2012), Ag diffusion in high-energy grain boundaries of SiC using Gaussian-distributed site energies and KRA barriers (Ko et al., 2016), and vacancy diffusion in non-dilute Ni-X alloys using a diffusion-window model for local barriers (Grabowski et al., 2018).

Event-driven reaction-diffusion methods form another adjacent tradition. For diffusion-annihilation with spatially varying annihilation rates, an exact protective-domain algorithm samples a trial annihilation time from a maximal rate $$\partial_t \rho + \nabla_x\cdot(\rho U) = \nabla_x\cdot(D\nabla_x \rho),$$9, a first-exit time from a free-diffusion Dirichlet propagator, and accepts annihilation with probability $D(x)=\int (v-U)\otimes (v-U)M\,dv / R(x),$0; summing over rejected trials reproduces the Dyson expansion of $D(x)=\int (v-U)\otimes (v-U)M\,dv / R(x),$1 (Schwarz et al., 2012). First-Passage KMC and Dynamic-Lattice FPKMC similarly propagate particles between analytically or locally discretized protective-domain events, extending exact diffusion sampling to drift-diffusion under background potentials (Mauro et al., 2013). These methods share with KDMC the aim of replacing many small diffusive hops by larger statistically correct events, but their analytical basis is first-passage decomposition rather than asymptotic preservation for kinetic transport.

This broader view clarifies a final point. “Kinetic-diffusion Monte Carlo” is best understood not as one universally standardized algorithm, but as a cluster of methods organized around a common problem: how to simulate systems that are microscopically kinetic yet macroscopically diffusive without paying the full cost of resolving every microscopic event. In the most developed present usage, KDMC for Boltzmann-BGK transport provides a particularly explicit solution to that problem by embedding kinetic free flight, Maxwellian velocity refreshment, and moment-matched diffusion in a single particle update, then extending that construction to multilevel estimators, exact fluid boundary conditions, multidimensional transport, and particle-level fluid-kinetic decomposition (Mortier et al., 2020, Mortier et al., 2020, Steel et al., 4 Sep 2025, Tang et al., 22 Jun 2026).