# Mass transport computations via correlation splitting and a law of total
diffusion

Manuel Athènes  Université Paris-Saclay, CEA, Service de Recherches de
Métallurgie Physique, 91191, Gif-sur-Yvette, France  Gilles Adjanor  Groupe
Métallurgie, MMC, EDF, Les Renardières, 77818 Moret-sur-Loing, France  Jérôme
Creuze  Université Paris-Saclay, ICMMO, CNRS UMR 8182, F91405 Orsay Cedex,
France

###### Abstract

Directly computing mass transport coefficients in stochastic models requires
integrating over time the equilibrium correlations between atomic
displacements. Here, we show how to accelerate the computations via
correlation splitting and conditioning, which statistically amounts to
estimating the mass transport coefficients through a law of total diffusion.
We illustrate the approach with kinetic path sampling simulations of atomic
diffusion in a random alloy model in which percolating solute clusters trap
the mediating vacancy. There, Green functions serve to generate first-passage
paths escaping the traps and to propagate the long-time dynamics. When they
also serve to estimate mean-squared displacements via conditioning, colossal
reductions of statistical errors are achieved.

###### pacs:

05.10.Ln

###### pacs:

07.05.Tp

###### pacs:

05.20.Jj

##  I Introduction

Mass transport is the natural phenomenon that governs both the time evolution
of thermodynamic systems towards equilibrium and the diffusion of chemical
species at equilibrium. Its understanding is fundamentally important and
challenging in many engineering applications, ranging from biophysics to
materials science. Equilibrium correlations between particle displacements
play a crucial role because they allow characterizing transport coefficients
close to equilibrium. To give a typical example, chemical currents in linear
response theory Chandler (1987); Evans and Morriss (2007) are expressed as
products of diffusion matrices and negative gradients of chemical
concentrations considered as small thermodynamic forces. Kinetic Monte Carlo
(kMC) Bortz et al. (1975) provides a direct method to measure transport
coefficients in multi-dimensional systems governed by master equations.
However, the method is hindered by the low occurrence of important events. A
considerable amount of computations is often needed to collect sufficient
statistics in many systems of interest Trochet et al. (2017, 2019), including
basic models of defect migration through homogeneous media Swinburne and Perez
(2020). Besides, low-complexity models are amenable to non-stochastic master-
equation approaches Koiwa and Ishioka (1983); Nastar (2014); Bocquet (2014);
Vaks et al. (2016) that may yield accurate predictions at much lower cost.
These methods involve integrating fundamental matrices Koiwa and Ishioka
(1983); Bocquet (2014); Schuler et al. (2020) or equivalently evaluating
lattice Green functions Trinkle (2017); Agarwal and Trinkle (2017) that either
minimize a variational problem Arita et al. (2017, 2018); Trinkle (2018) or
satisfy a Poisson equation Duflo (1997); Veretennikov (2017); Delmas and
Jourdain (2009); Athènes (2018). The dimension of the linear systems to solve
and the viability of approximations made to handle exponentially increasing
numbers of configurations limit the applicability of non-stochastic
approaches.

Here, we develop a correlation splitting and conditioning scheme aiming to
facilitating the computations of mass transport coefficients. The approach
alleviates the aforementioned issues: scarcity of harvested important events
and combinatorial explosion in numerical algebra. It furthermore leads us to
formulate a law of total diffusion (LTD) that relates to the laws of total
expectation (LTE) and variance (LTV) considered so far to compute
thermodynamic expectations and their statistical variance via conditioning
Frenkel (2004, 2006); Athènes (2007); Delmas and Jourdain (2009); Athènes and
Marinica (2010); Adjanor et al. (2011); Cao et al. (2014); Athènes and Terrier
(2017); Frenkel et al. (2017); Athènes (2018). We illustrate the approach on a
random alloy model Trinkle (2018) exhibiting dynamical trapping and
percolation by estimating the diffusion matrix, denoted below by
${\mathbb{D}}\left[\mathbf{d}\right]$ and defined as half the asymptotic
variance of the vector $\mathbf{d}$ of chemical displacements.

##  II Diffusion matrix for reversible Markov chains

To monitor the displacements of each of the $c$ chemical species along each of
the $s$ space dimensions, any state $\chi$ is mapped onto a descriptor
$\mathbf{r}\in\mathbb{R}^{sc}$. Component $r_{i}$ represents here the
coordinate sum of all $\alpha$-type atoms in space direction $a$, the
descriptor index being encoded as $i=a+s(\alpha-1)$, where
$a\in\mathbb{N}^{\star}_{s}$ and $\alpha\in\mathbb{N}^{\star}_{c}$. Then,
displacement $\mathbf{d}(\chi,\chi^{\prime})$ from state $\chi$ to state
$\chi^{\prime}$ is computed using the minimum image representation of
$\mathbf{r}^{\prime}-\mathbf{r}$ 111$\mathbf{d}(\chi,\chi^{\prime})$ equals
$\mathbf{r}^{\prime}-\mathbf{r}+\arg\min\|\mathbf{r}^{\prime}-\mathbf{r}+\mathbf{m}\|$,
where the components $m_{i}$ of argument $\mathbf{m}$ run over the multiples
of the simulation supercell periods.. The variance matrix of vector
$\mathbf{r}$ over a Markov chain reads
${\mathbb{V}}\left[\mathbf{r}\right]=\mathbb{E}\left[\mathbf{r}^{2}\right]-\mathbb{E}\left[\mathbf{r}\right]^{2}$,
with $\mathbf{r}^{2}$ standing for tensorial square
$\mathbf{r}\otimes\mathbf{r}$ and $\mathbb{E}\left[\cdot\right]$ denoting the
expectation operator. The successive states of the chain are denoted by
$\chi_{h}$ with $h\in\mathbb{N}$ or $\mathbb{Z}$. The displacement vector
after $\ell$ transitions,
$\sum_{h=0}^{\ell-1}\mathbf{d}\left(\chi_{h},\chi_{h+1}\right)$, is simply
written $\mathbf{r}_{0\shortrightarrow\ell}$ in the following. The expected
square displacement divided by twice the elapsed time is written

${\mathbb{D}_{\ell}}\left[\mathbf{d}\right]=\frac{1}{2}\frac{\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow\ell}\otimes\mathbf{r}_{0\shortrightarrow\ell}\right]}{\mathbb{E}\left[t_{0\shortrightarrow\ell}\right]}=\frac{1}{2}\frac{{\mathbb{V}}\left[\mathbf{r}_{0\shortrightarrow\ell}\right]}{\mathbb{E}\left[t_{0\shortrightarrow\ell}\right]}.$
(1)

The asymptotic limit yields the diffusion matrix, i.e.,
${\mathbb{D}}\left[\mathbf{d}\right]=\lim_{\ell\shortrightarrow\infty}{\mathbb{D}_{\ell}}\left[\mathbf{d}\right]$.
The variance amounts to squaring here because the expected displacement
$\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow\ell}\right]$ is zero at
equilibrium. The Markov chain obeys detailed balance with respect to an
equilibrium stationary distribution $\rho^{\mathrm{eq}}_{\chi}$. Fulfillment
of this strong condition entails the invariance of expectations under
arbitrary translation of chain indexes and under chain reversal. Expected
autocorrelations are thus invariant after interchange of displacements:
$\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}\otimes\mathbf{r}_{h\shortrightarrow h+1}\right]$ is successively equal to
$\mathbb{E}\left[\mathbf{r}_{-1-h\shortrightarrow-h}\otimes\mathbf{r}_{-1\rightarrow
0}\right]$, $\mathbb{E}\left[\mathbf{r}_{h+1\shortrightarrow
h}\otimes\mathbf{r}_{1\shortrightarrow 0}\right]$ and eventually
$\mathbb{E}\left[\mathbf{r}_{h\shortrightarrow
h+1}\otimes\mathbf{r}_{0\shortrightarrow 1}\right]$, displacements being anti-
symmetrical under chain reversal: $\mathbf{r}_{h+1\shortrightarrow
h}=-\mathbf{r}_{h\shortrightarrow h+1}$, $\forall h$. Hence,
${\mathbb{D}}\left[\mathbf{d}\right]$ is symmetric nonnegative. Translational
invariance also entails that the expected elapsed time is $\ell$ multiplied by
the mean elapsed time before the next-event
$\bar{\tau}=\mathbb{E}\left[t_{0\shortrightarrow 1}\right]$. The
autocorrelation symmetry properties enable us to split the diffusion matrix
into two parts:

$\displaystyle{\mathbb{D}}\left[\mathbf{d}\right]={\mathbb{D}_{1}}\left[\mathbf{d}\right]+\frac{1}{\bar{\tau}}\sum\nolimits_{h=1}^{\infty}\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}\otimes\mathbf{r}_{h\shortrightarrow h+1}\right].$ (2)

For uncorrelated chains, the _uncorrelated_ part
${\mathbb{D}_{1}}\left[\mathbf{d}\right]=\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}\otimes\mathbf{r}_{0\shortrightarrow 1}\right]/\left[2\bar{\tau}\right]$
contributes to diffusion exclusively, since the _correlated_ part encompassing
the summation vanishes. For low-dimensional spaces, $\bar{\tau}$ and
${\mathbb{D}_{1}}\left[\mathbf{d}\right]$ can be readily evaluated from the
knowledge of $\rho_{\chi}^{\mathrm{eq}}$ via conditioning on the states. Let
$\mathbb{E}\left[f(\chi_{k},\chi_{k+1})|\chi_{0}\right]$ denote the
conditional expectation of function $f(\chi_{k},\chi_{k+1})$ given $\chi_{0}$
and integer $k\geq 0$. The mean first-passage time (MFPT) from $\chi_{0}$ then
writes $\tau_{0}=\mathbb{E}\left[t_{0\shortrightarrow 1}|\chi_{0}\right]$ and
its mean yields $\bar{\tau}=\mathbb{E}\left[\tau_{0}\right]$. Term
$\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}\otimes\mathbf{r}_{0\shortrightarrow 1}\right]$ can be similarly evaluated,
yielding the uncorrelated part. Evaluating the correlated part is however more
difficult because it requires formulating and solving a Poisson equation. For
this purpose, define the mean displacement from $\chi_{1}$ as

$\displaystyle\mathbf{e}(\chi_{1})=\mathbb{E}\left[\mathbf{r}_{1\shortrightarrow
2}|\chi_{1}\right]=-\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}|\chi_{1}\right],$ (3)

and the associated relaxation vector
$\bm{\epsilon}(\chi_{1})=\sum\nolimits_{h=1}^{+\infty}\mathbb{E}\left[\mathbf{e}(\chi_{h})|\chi_{1}\right]$
whose knowledge will give access to the correlated part of the diffusion
matrix. The relaxation vector is a particular solution of the following
discrete Poisson equation

$\displaystyle\bm{\epsilon}(\chi_{1})=\mathbb{E}\left[\bm{\epsilon}(\chi_{2})|\chi_{1}\right]+\mathbf{e}(\chi_{1})$
(4)

where $\chi_{1}$ runs over the state space. The solution $\bm{\epsilon}$ is
fully determined by additionally imposing
$\mathbb{E}\left[\bm{\epsilon}\right]=\mathbb{E}\left[\mathbf{e}\right]=0$. It
characterizes the expected correlations:

$\displaystyle\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}\otimes\bm{\epsilon}(\chi_{1})\right]$
$\displaystyle=\sum\nolimits_{h=1}^{\infty}\mathbb{E}\left[\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}\otimes\mathbf{r}_{h\shortrightarrow h+1}|\chi_{0},\chi_{1}\right]\right],$
$\displaystyle=\sum\nolimits_{h=1}^{\infty}\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}\otimes\mathbf{r}_{h\shortrightarrow h+1}\right].$ (5)

Combining (2) and (5) then yields the diffusion matrix:

$\displaystyle\begin{aligned} {\mathbb{D}}\left[\mathbf{r}_{0\shortrightarrow
1}\right]&={\mathbb{D}_{1}}\left[\mathbf{r}_{0\shortrightarrow
1}\right]+\frac{1}{\bar{\tau}}\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}\otimes\bm{\epsilon}(\chi_{1})\right]\\\
&=\frac{1}{2\bar{\tau}}\mathbb{E}\left[\left(\mathbf{r}_{0\shortrightarrow
1}+\bm{\epsilon}(\chi_{1})\right)^{2}-\bm{\epsilon}^{2}\right],\end{aligned}$
(6)

where first equality is similar to Green-Kubo formulae Arita et al. (2017,
2018) and symmetry is formally recovered in second equality. Expressions in
(6) are the cornerstone of variational approaches to mass transport Arita et
al. (2017, 2018); Trinkle (2018) and serve here the purpose of improving kMC
measurements via correlation splitting and conditioning.

##  III Correlation splitting and conditioning

We proceed by first writing the LTV for displacement
$\mathbf{r}_{1\shortrightarrow 2}$ with conditioning on $\chi_{1}$ and
rescaling with $2\bar{\tau}$,

$\displaystyle{\mathbb{D}_{1}}\left[\mathbf{r}_{1\shortrightarrow
2}\right]=\mathbb{E}\left[{\mathbb{D}_{1}}\left[\mathbf{r}_{1\shortrightarrow
2}|\chi_{1}\right]\right]+{\mathbb{D}_{1}}\left[\mathbb{E}\left[\mathbf{r}_{1\shortrightarrow
2}|\chi_{1}\right]\right],$ (7)

where ${\mathbb{D}_{1}}\left[\mathbf{r}_{1\shortrightarrow
2}|\chi_{1}\right]={\mathbb{V}}\left[\mathbf{r}_{1\shortrightarrow
2}|\chi_{1}\right]/\left[2\bar{\tau}\right]$ and the conditional variance is
${\mathbb{V}}\left[\mathbf{r}_{1\shortrightarrow
2}|\chi_{1}\right]=\mathbb{E}\left[(\mathbf{r}_{1\shortrightarrow
2})^{2}|\chi_{1}\right]-\mathbf{e}(\chi_{1})^{2}$. Further splitting the
diffusion matrix (6) is then obtained through: (i) inserting the LTV (7) in
(6) and conditioning the expectation of $\mathbf{r}_{0\shortrightarrow
1}\otimes\bm{\epsilon}(\chi_{1})$ on $\chi_{1}$; (ii) inserting the
reversibility property (3); (iii) plugging the Poisson equation (4); (iv)
simplifying the plugged expectation
$\mathbb{E}\left[\mathbb{E}\left[\mathbf{e}(\chi_{1})\otimes\bm{\epsilon}(\chi_{2})|\chi_{1}\right]\right]$
into
$\mathbb{E}\left[\mathbf{e}(\chi_{1})\otimes\bm{\epsilon}(\chi_{2})\right]$
using the LTE; (v) regrouping like quantities
${\mathbb{D}_{1}}\left[\mathbf{e}\right]=\mathbb{E}\left[\mathbf{e}^{2}\right]/\left[2\bar{\tau}\right]$;
(vi) identifying the diffusion matrix of the conditionally expected
displacements,

$\displaystyle\begin{aligned}
{\mathbb{D}}\left[\mathbf{e}(\chi_{1})\right]&={\mathbb{D}_{1}}\left[\mathbf{e}(\chi_{1})\right]+\frac{1}{\bar{\tau}}\mathbb{E}\left[\mathbf{e}(\chi_{1})\otimes\bm{\epsilon}(\chi_{2})\right]\\\
&=\frac{1}{2\bar{\tau}}\mathbb{E}\left[\left(\mathbf{e}(\chi_{1})+\bm{\epsilon}(\chi_{2})\right)^{2}-\bm{\epsilon}^{2}\right].\end{aligned}$
(8)

by analogy to expressions of Eq. (6). Both contributions (6) and (8) are
measurable by kMC simulations. The resulting splitting yields the LTD:

$\displaystyle{\mathbb{D}}\left[\mathbf{r}_{1\shortrightarrow 2}\right]$
$\displaystyle=\mathbb{E}\left[{\mathbb{D}_{1}}\left[\mathbf{r}_{1\shortrightarrow
2}|\chi_{1}\right]\right]-{\mathbb{D}}\left[\mathbb{E}\left[\mathbf{r}_{1\shortrightarrow
2}|\chi_{1}\right]\right],$ (9)

in which the diffusion matrix is expressed as the difference between two
symmetric nonnegative contributions, the _intra_ \- and _extra_ -correlated
diffusion matrices, respectively. Intra-correlations involve consecutive
displacements exclusively, since
$\mathbb{E}\left[{\mathbb{D}_{1}}\left[\mathbf{r}_{1\shortrightarrow
2}|\chi_{1}\right]\right]={\mathbb{D}_{2}}\left[\mathbf{r}_{0\shortrightarrow
1}\right]$. The remaining correlations contribute to the extra-correlated
part,
${\mathbb{D}}\left[\mathbf{e}(\chi_{0})\right]=-\mathbb{E}\left[\mathbf{r}_{0\shortrightarrow
1}\otimes\left(\mathbf{r}_{1\shortrightarrow
2}+2\sum\nolimits_{h=2}^{\infty}\mathbf{r}_{h\shortrightarrow
h+1}\right)\right]/[2\bar{\tau}]$. The former equality arises as the
particular case $\ell=2$ of a more general relationship,

$\displaystyle{\mathbb{D}_{\ell}}\left[\mathbf{d}\right]=\mathbb{E}\left[{\mathbb{D}_{1}}\left[\mathbf{d}|\chi\right]\right]+\tfrac{1}{\ell}{\mathbb{D}_{1}}\left[\mathbb{E}\left[\mathbf{d}|\chi\right]\right]-\tfrac{\ell-1}{\ell}{\mathbb{D}_{\ell-1}}\left[\mathbb{E}\left[\mathbf{d}|\chi\right]\right],$
(10)

derived from the reversibility property (3) and the law of total covariance
not , where $\chi\equiv\chi_{1}$ and
$\mathbf{d}\equiv\mathbf{r}_{0\shortrightarrow 1}$. Relationship (10) bridges
between the LTV at the lower extremity and the LTD obtained for
$\ell\rightarrow\infty$. The LTD entails that the diffusion matrix can be
estimated from mean local quantities, by plugging the expected displacements
and their conditional variances given the visited states into its extra- and
intra-correlated parts. Besides, the LTD and bridging law are meaningful for
reversible Markov chains and for any stochastic variable that is antisymmetric
under chain reversal.  Laws (7) and (9) also yield a Löwner partial ordering:

$\displaystyle{\mathbb{D}}\left[\mathbf{d}\right]\leq\mathbb{E}\left[{\mathbb{D}_{1}}\left[\mathbf{d}|\chi\right]\right]\leq{\mathbb{D}_{1}}\left[\mathbf{d}\right],$

entailing that ${\mathbb{D}}\left[\mathbf{d}\right]$ is better approximated by
its intra-correlated part ${\mathbb{D}_{2}}\left[\mathbf{d}\right]$ than its
uncorrelated part ${\mathbb{D}_{1}}\left[\mathbf{d}\right]$. Further
increasing $\ell$ provides a decaying sequence
$\left({\mathbb{D}_{\ell}}\left[\mathbf{d}\right]\right)_{\ell\geq 1}$ of
upper-bound approximates of ${\mathbb{D}}\left[\mathbf{d}\right]$ Trinkle
(2018). We then estimate ${\mathbb{D}}\left[\mathbf{d}\right]$ via
${\mathbb{D}_{\ell}}\left[\mathbf{d}\right]$ in (10) over a sample of $l$
trajectories of $L=\ell_{\mathrm{max}}$ steps each. For statistical errors to
be small, $l$ must be large enough, standard deviations decaying as
$1/\sqrt{l}$ Caflisch (1998). Letting $\chi_{i,h}$ denote state $h$ of
trajectory $i$ and $\mathbf{V}(\chi_{i,h})$ stand for
${\mathbb{V}}\left[\mathbf{r}_{i,h\shortrightarrow h+1}|\chi_{i,h}\right]$,
the LTD-based conditioned estimator of
${\mathbb{D}_{\ell}}\left[\mathbf{d}\right]$ writes

$\displaystyle\mathbf{D}_{l,\ell}^{\mathrm{cnd}}$
$\displaystyle=\frac{1}{2\widehat{\tau}}\cdot\frac{1}{lL}\sum\nolimits_{i=1}^{l}\sum\nolimits_{h=0}^{L-1}\left[\mathbf{V}(\chi_{i,h})+\tfrac{1}{\ell}\mathbf{e}(\chi_{i,h})^{2}\right]$
$\displaystyle-\frac{1}{2\widehat{\tau}}\cdot\frac{1}{l\ell}\sum\nolimits_{i=1}^{l}\left(\sum\nolimits_{h=1}^{\ell-1}\mathbf{e}(\chi_{i,h})\right)^{2},$
(11)

where $\widehat{\tau}$ is an estimate of $\bar{\tau}$ obtained via LTE-based
conditioning, i.e., by averaging the MFPTs Athènes et al. (1997); Puchala et
al. (2010):

$\displaystyle\widehat{\tau}=\frac{1}{lL}\sum\nolimits_{i=1}^{l}\sum\nolimits_{h=0}^{L-1}\mathbb{E}\left[t_{i,h\shortrightarrow
h+1}|\chi_{i,h}\right].$ (12)

Note that information about the $lL$ gathered states is included in the
estimation of static expectations in (III) and (12). Conditioned estimator in
(12) has a statistical variance that is lower than that of the plain estimator
$\sum_{i=1}^{l}t_{i,0\shortrightarrow\ell}$, a property guaranteed by the LTV
Frenkel (2006); Athènes and Terrier (2017); Athènes (2018). The conditioning
over time is traditionally done in the standard estimator of
${\mathbb{D}_{\ell}}\left[\mathbf{d}\right]$:

$\displaystyle\mathbf{D}^{\mathrm{std}}_{l,\ell}=\frac{1}{2\widehat{\tau}}\cdot\frac{1}{l\ell}\sum\nolimits_{i=1}^{l}\mathbf{r}_{i,0\shortrightarrow\ell}\otimes\mathbf{r}_{i,0\shortrightarrow\ell}.$
(13)

While estimators (III) and (13) are both valid, their statistical variances
usually differ. Hereafter, the optimal linear combination of the standard and
conditioned estimators will be employed as a third estimator, as in waste-
recycling Monte Carlo Delmas and Jourdain (2009); Adjanor et al. (2011),
because it provides maximum reduction of the statistical variance via a
control variate Caflisch (1998).

##  IV Auxiliary absorbing Markov chains

Correlation splitting and conditioning can be employed in conventional or
advanced kMC simulations. Here, we present simulations of severe dynamical
trapping using the kinetic path sampling (kPS) algorithm Athènes and Bulatov
(2014); Athènes et al. (2019); Sharpe and Wales (2020). kPS is an accelerated
kMC algorithm developed upon theories of absorbing Markov chains (AMC) and
graph transformation Wales (2009). It performs non-local displacements
avoiding the most trapping states. Those are formally pooled together into a
prespecified set, called _transient set_ using AMC terminology and denoted by
$\mathtt{T}$. Moreover, kPS algorithm is able to generate first-passage paths
and times efficiently based on matrix factorizations or inversions. Other AMC-
kMC algorithms rely either on matrix diagonalisations Novotny (1995);
Opplestrup et al. (2006); Oppelstrup et al. (2009); Donev et al. (2010) or a
mean rate ersatz Athènes et al. (1997); Puchala et al. (2010) based on Eq.
(12). The kMC and kPS stochastic matrices, respectively denoted by
$P^{\mathrm{o}}$ and $P$, satisfy the conditions $P^{\mathrm{o}}_{\chi\chi}=0$
and $P_{\chi\zeta}=0$, $\forall\chi$,
$\forall\zeta\in\mathtt{T}\cup\\{\chi\\}$. The probabilities of not
transitioning are always zero from any state and the respective MFPTs, denoted
by $\tau^{\mathrm{o}}_{\chi}$ and $\tau_{\chi}$, are state-dependent. Matrix
$P$ is constructed from $P^{\mathrm{o}}$ and an absorbing stochastic matrix,
$P^{\mathrm{a}}$ coinciding with $P^{\mathrm{o}}$ over $\mathtt{T}$ and with
identity matrix $I$ over $\bar{\mathtt{T}}$, the complement of $\mathtt{T}$.
This entails that transitions from $\bar{\mathtt{T}}$ are cancelled, making
$\mathtt{T}$-states transient and $\bar{\mathtt{T}}$-states recurrent:

$\displaystyle
P^{\mathrm{a}}_{\chi\chi^{\prime}}=\bm{1}_{\mathtt{T}}(\chi)P^{\mathrm{o}}_{\chi\chi^{\prime}}+\bm{1}_{\mathtt{\bar{T}}}(\chi)I_{\chi\chi^{\prime}}$
(14)

where indicator function $\bm{1}_{\mathtt{S}}$ is 1 if $\chi\in\mathtt{S}$ and
0 otherwise.
Resorting to AMC theory Novotny (1995); Athènes et al. (1997); Wales (2009);
Opplestrup et al. (2006); Oppelstrup et al. (2009); Donev et al. (2010);
Puchala et al. (2010); Athènes and Bulatov (2014); Athènes et al. (2019);
Sharpe and Wales (2020), the mean number of visits of transient state
$\zeta^{\prime}$ is defined as the conditional expectation of the indicator
function $\bm{1}_{\\{\zeta^{\prime}\\}}$ given initial state
$\chi_{1}\in\mathtt{T}$. It corresponds to the Green function $\chi_{1}\mapsto
G^{\mathrm{a}}_{\chi_{1}\zeta^{\prime}}=\sum_{k=0}^{\infty}\mathbb{E}^{\mathrm{a}}\left[\bm{1}_{\\{\zeta^{\prime}\\}}(\chi_{k+1})|\chi_{1}\right]$
solution of the Poisson equation over $\mathtt{T}$:

$\displaystyle
G^{\mathrm{a}}_{\chi_{1}\zeta^{\prime}}=\mathbb{E}^{\mathrm{a}}\left[G^{\mathrm{a}}_{\chi_{2}\zeta^{\prime}}|\chi_{1}\right]+\bm{1}_{\\{\zeta^{\prime}\\}}(\chi_{1}),$
(15)

entailing
$G^{\mathrm{a}}_{\zeta\zeta^{\prime}}=\sum_{\xi\in\mathtt{T}}P^{\mathrm{o}}_{\zeta\xi}G^{\mathrm{a}}_{\xi\zeta^{\prime}}+I_{\zeta\zeta^{\prime}}$,
while $G^{\mathrm{a}}_{\chi\chi^{\prime}}$ is set to zero whenever $\chi$ or
$\chi^{\prime}$ lies outside $\mathtt{T}$. Thus, $G^{\mathrm{a}}$ and
$I-P^{\mathrm{a}}$ are two singular matrices, pseudoinverse of each other,
whereas $G^{\mathrm{a}}_{\zeta\zeta^{\prime}}$ is the regular inverse of
$I_{\zeta\zeta^{\prime}}-P^{\mathrm{a}}_{\zeta\zeta^{\prime}}$ for
$\zeta,\zeta^{\prime}\in\mathtt{T}$. For an absorbing Markov chain starting
from $\zeta$, the knowledge of the number of visits of state $\zeta^{\prime}$
enables one to deduce the probability of being absorbed at
$\chi^{\prime}\in\bar{\mathtt{T}}$ and the mean time for being absorbed in
$\bar{\mathtt{T}}$:

$\displaystyle\begin{aligned}
\varPi^{\mathrm{a}}_{\zeta\chi^{\prime}}&=\sum\nolimits_{\zeta^{\prime}}G_{\zeta\zeta^{\prime}}^{\mathrm{a}}P_{\zeta^{\prime}\chi^{\prime}}^{\mathrm{a}}+I_{\zeta\chi^{\prime}},\\\
\tau_{\zeta}^{\mathrm{a}}&=\sum\nolimits_{\zeta^{\prime}}G_{\zeta\zeta^{\prime}}^{\mathrm{a}}\tau_{\zeta^{\prime}}^{\mathrm{o}}.\end{aligned}$
(16)

The overlying transition probability from $\chi$ to $\chi^{\prime}\neq\chi$
and the MFPT from $\chi$ eventually write

$\displaystyle\begin{aligned}
P_{\chi\chi^{\prime}}&=\sum\nolimits_{\zeta}P^{\mathrm{o}}_{\chi\zeta}\varPi^{\mathrm{a}}_{\zeta\chi^{\prime}}\big{/}\pi_{\chi}\\\
\tau_{\chi}&=\sum\nolimits_{\zeta}P^{\mathrm{o}}_{\chi\zeta}\left[\tau^{\mathrm{o}}_{\chi}+\tau^{\mathrm{a}}_{\zeta}\right]\big{/}\pi_{\chi},\end{aligned}$
(17)

where rescaling of both transition probabilities and elapsed times by
probability
$\pi_{\chi}=1-\sum_{\zeta}P^{\mathrm{o}}_{\chi\zeta}\varPi^{\mathrm{a}}_{\zeta\chi}$
serves to cancel the flicker probability, thus $P_{\chi\chi}=0$. Besides, if
$P^{\mathrm{o}}_{\chi\chi^{\prime}}$ obeys detailed balance w.r.t its
stationary probability distribution $\rho^{\mathrm{o}}_{\chi}$ over
$\bar{\mathtt{T}}\cup\mathtt{T}$, then so does $P_{\chi\chi^{\prime}}$ w.r.t
$\rho_{\chi}\propto\rho^{\mathrm{o}}_{\chi}\pi_{\chi}$ over
$\bar{\mathtt{T}}$. The scaled probability distribution
$\rho^{\mathrm{eq}}_{\chi}\propto\rho^{\mathrm{o}}_{\chi}\tau^{\mathrm{o}}_{\chi}\propto\rho_{\chi}\tau^{\mathrm{o}}_{\chi}/\pi_{\chi}$
then corresponds to the thermodynamic equilibrium. Fulfilment of the
reversibility property allows conditioning to be performed. The conditioned
estimator (III) can thus be used. Note that the kPS algorithm reverts to the
conventional kMC algorithm when none of the states are made transient.

Figure 1: Transformation of transition networks: (a) physical network wherein
blue and black barriers map out exchanges of a vacancy (empty squares) with
bulk $A$-atoms (not displayed) or solute $B$-atoms (red disks); (b) augmented
network after inserting saddle states along vacancy-$A$ exchanges (empty
diamonds represent vacancies at saddle positions); (c) reduced network after
eliminating stable states.

##  V Random alloy model

We illustrate the approach by computing the mass transport coefficients in a
random binary alloy on a square lattice with periodic conditions. Diffusion is
mediated by a single vacancy $V$ exchanging with $A$ or $B$ nearest-neighbor
atoms Allnatt and Lidiard (1993); Barbe and Nastar (2006); Trinkle (2018).
Transition rates are usually computed in the framework of classical transition
state theory Hänggi et al. (1990). They have the form
$\nu_{X}=\nu^{\star}\exp[-E_{X}/(k_{B}T)]$ where $E_{X}>0$ is the energy
barrier that neighbouring atom $X$ jumping into the vacancy must cross, $T$
denotes temperature, $k_{B}$ stands for Boltzmann’s constant and $\nu^{\star}$
is the attempt frequency for both species. Here, transition rates are
environment-independent and only depend on the type of the jumping atom.
Consequently, all states are equiprobable, the system energy has same constant
value for all stable states. No thermodynamic transition occurs in the random
solid solution Allnatt and Lidiard (1993). In particular, the site percolation
threshold is independent of the jumping frequencies $\nu_{A}$ and $\nu_{B}$.
Despite its simplicity, the random alloy model is a nontrivial system and is
thus customarily used to study trapping and percolation Barbe and Nastar
(2006); Trinkle (2018). Assuming $E_{A}>E_{B}>0$, the frequency ratio
$\nu_{B}/\nu_{A}$ increases with decreasing the temperature. Hence, low
temperatures result in dynamical trapping. Here, the jump rate of $B$-atoms,
$\nu_{B}$, will be much higher than the one of $A$-atoms, $\nu_{A}$. For the
sake of simplicity, we adopt $\nu_{A}$ units by setting this value to 1 in the
following. To perform kPS simulations, the transient graph is to be crafted as
the union of many disconnected subgraphs, so that relatively small blocks of
$I-P^{\mathrm{a}}$ need being numerically manipulated on the fly (for
computing $\varPi^{\mathrm{a}}$ and $\tau^{\mathrm{a}}$). We illustrate the
construction of the simulated networks on a square lattice in Fig. 1. Trapping
is suppressed by making transient all the states where the vacancy is located
next to a $B$-atom. The difficulty is that vacancy-solute cluster shapes are
connected to each others (see Fig. 1.a), rendering the size of the local
transient block so huge that its enumeration is numerically impracticable
except for small and isolated $B$-clusters. The problem is mitigated by first
augmenting the transition network through the insertion of saddle states along
$AV$ exchanges to divide the trapping network into many disconnected subgraphs
(see Fig. 1.b). Transition probabilities from any added saddle state to its
two adjacent stable states both equal half. Besides, doubling the transition
rates from a stable state to its adjacent saddle states and nullifying the
residence times $\tau^{\mathrm{o}}$ at the added saddles compensates the
effect of the occasional flickers between stable and saddle states. Finally,
stable states are all turned transient, as diagrammed in Fig. 1.c. That is,
set $\mathtt{T}$ is the space of stable states and the simulated network
$\bar{\mathtt{T}}$ exclusively consists of the inserted saddle states. The
augmentation/reduction of the network leaves the distributions of the
sequences of simulated events and times unaffected.

Figure 2: A supercell (in bold) is periodically replicated horizontally. The
blue and green paths represent two distinct sequences of atomic displacements
through the percolating solute cluster (the transient sub-graph) between the
same starting state in (a) and ending state in (b).

Computing the conditional cumulants given the sampled states in estimator
(III) is strenuous when a solute cluster, i.e. a sub-graph of the transient
network, percolates through the periodically replicated supercell, as
illustrated in Fig. 2. Cumulants $\mathbf{e}(\chi_{0})$ and
$\mathbf{V}(\chi_{0})$ must account for all the displacements
$\mathbf{r}_{0\shortrightarrow 1}$ possibly sampled by kPS algorithm. The
sampled displacements starting in $\chi_{0}$ and ending in
$\chi_{1}\in\bar{\mathtt{T}}$ are all congruent to
$\mathbf{d}(\chi_{0},\chi_{1})$ modulo the cell periods along the percolating
directions (see Fig. 2.a), and many among them may differ from
$\mathbf{d}(\chi_{0},\chi_{1})$. This consequent sampling burden is avoided by
analytically integrating the underlying conditional expected displacements
$\mathbf{e}^{\mathrm{a}}(\chi_{h})\triangleq\mathbb{E}^{\mathrm{a}}\left[\mathbf{r}^{\mathrm{u}}_{h\shortrightarrow
h+1}|\chi_{h}\right]$, in a similar way to Athènes et al. (2019); Swinburne
and Perez (2020). Letting $\mathbb{E}^{\mathrm{o}}\left[\cdot\right]$ denote
expectation w.r.t. $P^{\mathrm{o}}$, both cumulants of overlying displacement
$\mathbf{r}_{0\shortrightarrow 1}$ are determined by the underlying sequences
$\mathbf{r}^{\mathrm{u}}_{0\shortrightarrow
1}+\mathbf{r}^{\mathrm{u}}_{1\shortrightarrow\infty}\equiv\mathbf{r}_{0\shortrightarrow
1}$ associated with the successive states of the AMC,
$\\{\chi_{h}^{\mathrm{u}}\\}_{h\geq 1}$, and where
$\chi_{0}^{\mathrm{u}}=\chi_{0}$ and $\chi_{\infty}^{\mathrm{u}}=\chi_{1}$.

Conditioning after the $P^{\mathrm{o}}$-transition within decomposition (17)
yields

$\displaystyle\mathbf{e}(\chi_{0})$
$\displaystyle=\mathbb{E}^{\mathrm{o}}\left[\mathbf{r}^{\mathrm{u}}_{0\shortrightarrow
1}+\mathbb{E}^{\mathrm{a}}\left[\mathbf{r}^{\mathrm{u}}_{1\shortrightarrow\infty}|\chi_{1}\right]|\chi_{0}\right]/\pi_{\chi_{0}},$
(18a) $\displaystyle\mathbf{V}(\chi_{0})$
$\displaystyle=\mathbb{E}^{\mathrm{o}}\left[\left[\mathbf{r}^{\mathrm{u}}_{0\shortrightarrow
1}+\mathbb{E}^{\mathrm{a}}\left[\mathbf{r}^{\mathrm{u}}_{1\shortrightarrow\infty}|\chi_{1}\right]-\mathbf{e}(\chi_{0})\right]^{2}\big{|}\chi_{0}\right]/\pi_{\chi_{0}}.$
(18b)

The mean displacement
$\mathbb{E}^{\mathrm{a}}\left[\mathbf{e}^{\mathrm{a}}\right]$ being null, the
relaxation vector
$\bm{\epsilon}^{\mathrm{a}}(\chi_{1})\triangleq\mathbb{E}^{\mathrm{a}}\left[\mathbf{r}_{1\shortrightarrow\infty}|\chi_{1}\right]=\sum\nolimits_{h=1}^{+\infty}\mathbb{E}^{\mathrm{a}}\left[\mathbf{e}^{\mathrm{a}}(\chi_{h})|\chi_{1}\right]$,
is null outside $\mathtt{T}$ and satisfies a Poisson equation inside
$\mathtt{T}$

$\displaystyle\bm{\epsilon}^{\mathrm{a}}(\chi_{1})=\mathbb{E}^{\mathrm{a}}\left[\bm{\epsilon}^{\mathrm{a}}(\chi_{2})|\chi_{1}\right]+\mathbf{e}^{\mathrm{a}}(\chi_{1}).$
(19)

The solution is unique, reads
$\bm{\epsilon}^{\mathrm{a}}(\zeta)=\sum_{\zeta^{\prime}}G^{\mathrm{a}}_{\zeta\zeta^{\prime}}\mathbf{e}^{\mathrm{a}}(\zeta^{\prime})$
for $\zeta\in\mathtt{T}$ and cancels outside ${\mathtt{T}}$. It serves to
evaluate the two conditional cumulants (18) for estimator (III). We resort to
their respective LTE and LTV expressions with conditioning over the underlying
absorbing chain initiated from $\zeta\equiv\chi_{1}$ to obtain a numerically
tractable expression

$\displaystyle\mathbf{e}(\chi)$
$\displaystyle=\sum\nolimits_{\zeta}\left[\mathbf{d}(\chi,\zeta)+\bm{\epsilon}^{\mathrm{a}}(\zeta)\right]{P^{\mathrm{o}}_{\chi\zeta}}/{\pi_{\chi}},$
(20a) $\displaystyle\mathbf{V}(\chi)$
$\displaystyle=\sum\nolimits_{\zeta}\left\\{\mathbf{V}^{\mathrm{a}}(\zeta)+\left[\mathbf{d}(\chi,\zeta)+\bm{\epsilon}^{\mathrm{a}}(\zeta)-\mathbf{e}(\chi)\right]^{2}\right\\}{P^{\mathrm{o}}_{\chi\zeta}}/\pi_{\chi},$
(20b)

where $\chi\equiv\chi_{0}$ in (18). The variance of the absorbing chain
conditioned on $\zeta\equiv\chi_{1}$ is also computed using the Green function
and relaxation vectors, resorting to relation (22) below:

$\displaystyle\mathbf{V}^{\mathrm{a}}(\chi_{1})$
$\displaystyle=\mathbb{E}^{\mathrm{a}}\left[\left.\left(\sum\nolimits_{h=1}^{\infty}\mathbf{r}^{\mathrm{u}}_{h\shortrightarrow
h+1}\right)^{2}\right|\chi_{1}\right]-\bm{\epsilon}^{\mathrm{a}}(\chi_{1})^{2}$
(21)
$\displaystyle=\sum\nolimits_{h=1}^{\infty}\mathbb{E}^{\mathrm{a}}\left[\left.(\mathbf{r}^{\mathrm{u}}_{h\shortrightarrow
h+1}+\bm{\epsilon}^{\mathrm{a}}(\chi_{h+1}))^{2}-\bm{\epsilon}^{\mathrm{a}}(\chi_{h})^{2}\right|\chi_{1}\right]$
$\displaystyle=\sum\nolimits_{\chi\zeta}G^{\mathrm{a}}_{\chi_{1}\chi}\left[P^{\mathrm{a}}_{\chi\zeta}\left[\mathbf{d}(\chi,\zeta)+\bm{\epsilon}^{\mathrm{a}}(\zeta)\right]^{2}-\bm{\epsilon}^{\mathrm{a}}(\chi)^{2}\right].$
(22)

Note that quantities $\bm{\epsilon}^{\mathrm{a}}$ and
$\mathbf{V}^{\mathrm{a}}$ are first and second derivatives of a displacement
cumulant generating function at its origin, respectively Swinburne and Perez
(2020).

##  VI Measurements of diffusion coefficients

We perform a series of kPS simulations of extreme dynamical trapping by
setting $\nu_{B}=10^{5}$ and $\nu_{A}=1$ in a periodically replicated square
lattice of size $16\times 16$. Compositions are gradually increased from 1 to
$254$ $B$-atoms. The atomic concentration of $X$ is denoted by $C_{X}$,
entailing that $C_{A+B}=C_{A}+C_{B}$ is one. The diffusion matrix
$\mathbb{D}(\mathbf{d})$ is estimated using $\mathbf{D}_{l,L}^{\mathrm{std}}$
and $\mathbf{D}_{l,L}^{\mathrm{cnd}}$ estimators with
$L=\ell_{\mathrm{max}}=200$, $l=5\cdot 10^{5}$ for $C_{B}<0.5$ and $l=10^{4}$
for $C_{B}\geq 0.5$. Estimates $D_{X}^{\mathrm{std}}$ and
$D_{X}^{\mathrm{cnd}}$ of the diffusion coefficients $D_{X}$ for both species
and the vacancy ($X=A$, $B$ and $A+B$) are then obtained via
$\mathbf{D}_{l,L}^{\mathrm{std}}$ and $\mathbf{D}_{l,L}^{\mathrm{cnd}}$ by
averaging the corresponding elements over the two space directions. The
optimal estimate is then obtained as
$D_{X}^{\mathrm{opt}}=(1-c^{\star})D_{X}^{\mathrm{std}}+c^{\star}D_{X}^{\mathrm{cnd}}$
where the optimal control variate is
$c^{\star}=-\mathrm{cov}\left(D_{X}^{\mathrm{std}},D_{X}^{\mathrm{cnd}}-D_{X}^{\mathrm{std}}\right)/\mathrm{var}\left(D_{X}^{\mathrm{cnd}}-D_{X}^{\mathrm{std}}\right)$.
We furthermore record the transient regimes
$D^{\mathrm{cnd}}_{X}(\ell)=\varSigma_{\ell}+\frac{1}{\ell}\varGamma_{1}+\frac{1-\ell}{\ell}\varGamma_{\ell-1}$
from (10) where $\varSigma_{\ell}$ and $\varGamma_{\ell-1}$ are estimates of
the corresponding intra- and extra-correlated parts. Last, we define and
evaluate the intra-to-extra correlation ratio
$\phi=\varSigma_{L}/\varGamma_{L}$ and the reduction factor of the statistical
variance
$\eta=\mathrm{var}\left(D^{\mathrm{std}}_{X}\right)/\mathrm{var}\left(D^{\mathrm{opt}}_{X}\right)$.

Figure 3: Atomic diffusivities $D^{\star}_{X}$ ($X=A$, $B$ and $A+B$) for
extreme dynamical trapping ($\nu_{A}=1$ and $\nu_{B}=10^{5}$), as a function
of $C_{B}$, with corresponding variance reduction factor $\eta$, intra-to-
extra correlation ratio $\phi$ and optimal variate $c^{\star}$. The supercell
size is $16\times 16$.

The atomic diffusivities $D^{\star}_{X}=D_{X}/C_{X}$ evaluated using the
standard and optimal estimators for the entire composition range are displayed
in Fig. 3, together with $\eta$, $\phi$ and $c^{\star}$. We observe that the
speedup measured in term of $\eta$ is considerable whenever the intra-
correlated contribution to diffusivity dominates ($\phi>>1$). In this
situation, the optimal estimator perfectly matches the conditioned estimator
($c^{\star}=1$). The optimal estimator is relatively less efficient for the
fast-diffusing $B$-atoms when the extra-correlated contribution to
$B$-diffusivity is significant, at $B$-concentration lower than 40%. The
lowest reductions of variance are about a factor of 2. To investigate the
origin of the relatively lower performance of the conditioned estimator, we
display in Fig. 4 the transient regime estimated from
$D^{\mathrm{std}}_{X}(\ell)$, $D^{\mathrm{cnd}}_{X}(\ell)$ and
$\varGamma_{\ell-1}$ with estimates of their standard errors for
$0\leq\ell\leq\ell_{\mathrm{max}}$. The errors associated with
$\varGamma_{\ell-1}$ terms are negligible compared to those with
$\varSigma_{\ell}$ terms and are attributed to the intermittent binding of the
vacancy to the $B$-clusters. Removing this intermittency in the kPS algorithm
would require making additional states transient and handling $A$-connected
subgraphs similarly to $B$-clusters.

Interestingly, the negative contributions
$\tfrac{1-\ell}{\ell}\varGamma_{\ell}$ in Fig. 4 undergo fast algebraic decays
towards their plateau values $\delta$. The three curves are very well fitted
by a law of the form $g(\ell+\ell_{0})^{\alpha}+\delta$. Implementing the
nonlinear least-squares Marquardt-Levenberg algorithm with the four fitting
parameters $g$, $\ell_{0}$, $\alpha$ and $\delta$ yields the power-law
exponents $\alpha$ associated with $A$, $B$ and $A+B$ diffusivities. The
obtained $\alpha$ values lie in the following $68\%$-confidence intervals:
$-1.302\pm 0.091$, $-1.015\pm 0.008$ and $-1.015\pm 0.008$, respectively. The
convergence features of the mean-squared mean-displacements entails that a
standard diffusive regime exhibiting short-range dependence is quickly reached
as $\ell$ increases. This property is related to the fast algebraic decay of
the autocorrelation function, as detailed in not . Note that the larger
exponent associated with $A$-diffusivity decreases as the frequency ratio is
decreased not . This physical trend results from a nontrivial interplay
between composition, transition rates, site interactions and geometry. We
refer to the recent work of Settino _et al._ Settino et al. (2020) who studied
anomalous diffusion in the Aubry-André model via the autocorrelation functions
and investigated the dependence of their power-law decay on the model physical
parameters.

Figure 4: Time dependence of the estimated $D_{X}(\ell)$ with
$\ell=t/\bar{\tau}$, using standard (13) and conditioned (III) estimators, for
$X=A$, $B$ and $A+B$. Filled areas around the curves represent $95\%$
confidence intervals (CI). CIs around $\tfrac{1}{\ell}\varGamma_{1}$ and
$\tfrac{1-\ell}{\ell}\varGamma_{\ell}$ curves are too small to be visible.

Here, the extremely severe dynamical trapping of the vacancy with $B$-atoms
strongly impacts the monotonous increase of $B$-diffusivity with $C_{B}$. The
sudden sharp increase of $D_{B}^{\star}$ observed at $C_{B}\approx 0.35$
occurs at a concentration much lower than $p_{s}=0.592745$, the site
percolation threshold Ziff (1992). This dynamical transition is due to the
cooperative motion of small $B$-clusters whose effectiveness increases
considerably with $C_{B}$ in small supercells. With increasing supercell
sizes, the transition further sharpens and occurs closer to the percolation
threshold, as shown in Fig. 5.

Figure 5: Diffusion coefficient $D_{A+B}$ as a function of composition for
various lattice sizes S using the kPS and kMC algorithms. The optimised
estimator is used in all measurements. The kPS simulations are performed with
same set-up as in Fig. 3, but four additional supercell sizes, $S=8$, 12, 24
and 32. The conventional kMC simulations employ $10^{4}$ trajectories of
$10^{6}$ steps. Their CPU costs are almost independent of the system sizes.
Note that the speed-ups $\eta$ obtained through conditioning are considerable
in conventional kMC simulations of dilute $A_{x}B_{1-x}$ systems with $x<0.1$.

The kPS simulations for the two largest supercell sizes in Fig. 5 became
inefficient and were stopped before terminating. The inefficiency results from
the larger sizes of the percolating solute clusters and to the cost of
computing the Green functions for the associated transient sub-graphs, a task
involving matrix inversion or factorisation. Here, a direct non-optimized
dense solver was used whose complexity scales cubically with the sizes of the
transient sub-graphs. As a result, the approach starts consuming CPU time
enormously when supercell sizes and $B$-concentrations are too large.
Nevertheless, dynamical trapping disappearing beyond the percolation
threshold, conventional kMC simulations can be used instead. One then obtains
reliable estimates of the diffusion coefficients at much larger system sizes,
as shown in Fig. 5. The variance reduction factor $\eta$ allowed by the
optimal estimator in kMC simulations increases from 1 below the percolation
threshold to $10^{5}$ in pure $B$. For concentrated alloy system
$A_{0.1}B_{0.9}$, the factor of variance reduction $\eta$ is $10^{2}$,
yielding an appreciable speed-up. However, the conditioned estimator
underperforms the standard one at the low-$B$ concentrations, when trapping
becomes important in conventional kMC simulations. Results obtained using the
conventional kMC algorithm in systems exhibiting moderate dynamical trapping
are also reported in not , where it is shown that significant variance
reductions via conditioning are systematically observed. The conventional kMC
algorithm outperforms kPS algorithm with excessively large transient sub-
graphs because the trapping strength measured in number of steps the vacancy
remains attached to a solute cluster scales linearly with the cluster size,
while the factorization cost scales cubically here. A smaller pre-factor for
the factorization costs makes kPS advantageous for transient sub-graphs
smaller than a certain threshold value. Hence, there are two complementary
ways to improve kPS simulations: (i) avoid turning a large vacancy-solute
cluster into a transient sub-graph when its size is above the threshold, (ii)
resort to a sparse linear solver to reduce the factorization/inversion costs.
The latter costs grow quadratically with the sub-graph sizes for the direct
multi-frontal linear solver used in Ref. Athènes et al. (2019).

Nonetheless, an interesting trend emerges from the results of the kPS
simulations with severe trapping. We indeed observe that the extra-to-intra
correlation ratio $\phi^{-1}$ vanishes for all diffusion measurements reported
in Fig. 3 when the solute concentrations are larger than the dynamical
percolation threshold values. This property means that the correlations of the
underlying process are well captured by the auxiliary absorbing Markov chains.
It entails that the extra-correlated part of the diffusion matrix is
negligible and justifies the use of the kinetic cluster expansion (KineCluE)
method Schuler et al. (2020) to solve this particularly difficult diffusion
problem. This method neglects the extra-correlated part and similarly resorts
to auxiliary absorbing Markov chains, but these ones are constructed from a
sub-graph expansion. To adapt the KineCluE method to concentrated alloys, the
absorbing Markov chains should be solved for a sample of alloy configurations
generated at thermodynamic equilibrium. Whenever the extra-correlated part
containing the cluster correlations is negligible, the transport coefficients
may be correctly estimated from their intra-correlated contributions
exclusively.

Concerning the transport properties of the slowly migrating solvent species,
we observe that $A$-diffusivity is almost not impacted by dynamical trapping
and remains constant. Reproducing a sluggish diffusion behaviour where
$A$-diffusivity exhibits a minimum at intermediate compositions, as observed
experimentally in several concentrated solid solution alloys Athènes et al.
(1997) and highly debated in the literature on high entropy alloys Tsai et al.
(2013); Paul (2017); Miracle and Senkov (2017); Osetsky et al. (2018); Kottke
et al. (2020); Daw and Chandross (2021), will necessitate to parameterize the
exchange rates $\nu_{A}$ and $\nu_{B}$, for instance by accounting for the
binding energies and entropies between alloying elements and vacancies Trochet
et al. (2017, 2019); Swinburne and Perez (2018); Lapointe et al. (2020). This
would affect short-range ordering or clustering, introduce a dependency of
vacancy concentration on alloy composition and impact atomic diffusion.

##  VII Conclusion

Splitting the diffusion matrix into its intra- and extra-correlated parts and
conditioning the correlations enables one to formulate a law of total
diffusion. A conditioned estimator then allows measuring atomic diffusivities
with considerably improved accuracy in absence of dynamical trapping. The use
of an optimized control variate guarantees systematic variance reductions
compared with standard computations performed with conventional or advanced
kinetic Monte Carlo methods. Conditioning is particularly useful when the
latter methods employ auxiliary absorbing Markov chains able to suppress
trapping at the underlying scale, like kinetic path sampling, as it leverages
the Green functions that are computed to generate the escaping first-passage
paths. However, the cost of computing these Green functions must remain
reasonable so that kinetic path sampling simulations always outperforms
conventional kinetic Monte Carlo simulations. If this condition is not met, it
is preferable to limit the sizes of the transient sub-graphs or to revert to
conventional kMC simulations. The conditioning approach is in principle
amenable to the numerous kMC simulations dealing with realistic potential
energy surface accounting for lattice distortions Trochet et al. (2019) or
describing amorphous structures Li and Ando (2020). In these simulations, the
computation of the transition rates is often so costly that the number of
simulated events is quite limited. The variance reduction technique proposed
here will certainly be useful in this particular context. As the conditioning
approach represents a small overhead, it can be readily applied to measure
transport properties of any particles migrating through thermally activated
events in crystalline or non-crystalline solids. Its range of applicability
encompasses surface and bulk diffusion in presence of extended defects
Grabowski et al. (2018); Athènes et al. (2019); Landeiro Dos Reis et al.
(2020); Rahman et al. (2021) like cavities, clusters, dislocations or grain
boundaries, within complex alloys or any material compounds.

###### Acknowledgements.

Stimulating discussions with Vasily Bulatov, Benjamin Jourdain, Maylise
Nastar, Kirone Mallick and Dallas Trinkle are gratefully acknowledged.

## References

  * Chandler (1987) D. Chandler, _Introduction to modern statistical mechanics_ (Oxford University Press, 1987).
  * Evans and Morriss (2007) D. J. Evans and G. P. Morriss, _Statistical Mechanics of Nonequilibrium Liquids_ (ANU Press, 2007), URL http://www.jstor.org/stable/j.ctt24h99q.
  * Bortz et al. (1975) A. Bortz, M. Kalos, and J. Lebowitz, _A new algorithm for monte carlo simulation of ising spin systems_ , Journal of Computational Physics 17, 10 (1975).
  * Trochet et al. (2017) M. Trochet, A. Sauvé-Lacoursière, and N. Mousseau, _Algorithmic developments of the kinetic activation-relaxation technique: Accessing long-time kinetics of larger and more complex systems_ , The Journal of Chemical Physics 147, 152712 (2017).
  * Trochet et al. (2019) M. Trochet, N. Mousseau, L. K. Béland, and G. Henkelman, _Off-Lattice Kinetic Monte Carlo Methods_ (Springer International Publishing, Cham, 2019), pp. 1–29, URL https://doi.org/10.1007/978-3-319-42913-7_29-1.
  * Swinburne and Perez (2020) T. D. Swinburne and D. Perez, _Automated calculation and convergence of defect transport tensors_ , npj Computational Materials 6, 190 (2020).
  * Koiwa and Ishioka (1983) M. Koiwa and S. Ishioka, _Integral methods in the calculation of correlation factors for impurity diffusion_ , Philosophical Magazine A 47, 927 (1983).
  * Nastar (2014) M. Nastar, _Atomic diffusion theory challenging the cahn-hilliard method_ , Phys. Rev. B 90, 144101 (2014).
  * Bocquet (2014) J.-L. Bocquet, _Correlation factor for diffusion in cubic crystals with solute–vacancy interactions of arbitrary range_ , Philosophical Magazine 94, 3603 (2014).
  * Vaks et al. (2016) V. G. Vaks, K. Y. Khromov, I. R. Pankratov, and V. V. Popov, _Statistical theory of diffusion in concentrated bcc and fcc alloys and concentration dependencies of diffusion coefficients in bcc alloys FeCu, FeMn, FeNi, and FeCr_ , Journal of Experimental and Theoretical Physics 123, 59 (2016).
  * Schuler et al. (2020) T. Schuler, L. Messina, and M. Nastar, _Kineclue: A kinetic cluster expansion code to compute transport coefficients beyond the dilute limit_ , Computational Materials Science 172, 109191 (2020).
  * Trinkle (2017) D. R. Trinkle, _Automatic numerical evaluation of vacancy-mediated transport for arbitrary crystals: Onsager coefficients in the dilute limit using a Green function approach_ , Philosophical Magazine 97, 2514 (2017).
  * Agarwal and Trinkle (2017) R. Agarwal and D. R. Trinkle, _Exact model of vacancy-mediated solute transport in magnesium_ , Phys. Rev. Lett. 118, 105901 (2017).
  * Arita et al. (2017) C. Arita, P. L. Krapivsky, and K. Mallick, _Variational calculation of transport coefficients in diffusive lattice gases_ , Phys. Rev. E 95, 032121 (2017).
  * Arita et al. (2018) C. Arita, P. L. Krapivsky, and K. Mallick, _Bulk diffusion in a kinetically constrained lattice gas_ , Journal of Physics A: Mathematical and Theoretical 51, 125002 (2018).
  * Trinkle (2018) D. R. Trinkle, _Variational principle for mass transport_ , Phys. Rev. Lett. 121, 235901 (2018).
  * Duflo (1997) M. Duflo, _Random Iterative Models_ , vol. 34 of _Applications of Mathematics_ (Springer-Verlag, Berlin Heidelberg, 1997).
  * Veretennikov (2017) A. Veretennikov, in _Modern Problems of Stochastic Analysis and Statistics_ , edited by V. Panov (Springer International Publishing, Cham, 2017), pp. 457–511.
  * Delmas and Jourdain (2009) J.-F. Delmas and B. Jourdain, _Does waste recycling really improve the multi-proposal Metropolis–Hastings algorithm? An analysis based on control variates_ , J. Appl. Probab. 46, 938 (2009).
  * Athènes (2018) M. Athènes, _Conditioning and enhanced sampling schemes for simulating thermodynamic and kinetic properties of condensed matter_ , Habilitation à diriger des recherches, Université Paris Saclay ; Université Paris Sud (2018), URL https://hal-cea.archives-ouvertes.fr/tel-01851686.
  * Frenkel (2004) D. Frenkel, _Speed-up of monte carlo simulations by sampling of rejected states_ , Proceedings of the National Academy of Sciences of the United States of America 101, 17571 (2004).
  * Frenkel (2006) D. Frenkel, _Waste-Recycling Monte Carlo_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006), pp. 127–137.
  * Athènes (2007) M. Athènes, _Web ensemble averages for retrieving relevant information from rejected monte carlo moves_ , The European Physical Journal B 58, 83 (2007).
  * Athènes and Marinica (2010) M. Athènes and M.-C. Marinica, _Free energy reconstruction from steered dynamics without post-processing_ , Journal of Computational Physics 229, 7129 (2010).
  * Adjanor et al. (2011) G. Adjanor, M. Athènes, and J. Rodgers, _Waste-recycling monte carlo with optimal estimates: Application to free energy calculations in alloys_ , The Journal of Chemical Physics 135, 044127 (2011).
  * Cao et al. (2014) L. Cao, G. Stoltz, T. Lelièvre, M.-C. Marinica, and M. Athènes, _Free energy calculations from adaptive molecular dynamics simulations with adiabatic reweighting_ , The Journal of Chemical Physics 140, 104108 (2014).
  * Athènes and Terrier (2017) M. Athènes and P. Terrier, _Estimating thermodynamic expectations and free energies in expanded ensemble simulations: Systematic variance reduction through conditioning_ , The Journal of Chemical Physics 146, 194101 (2017).
  * Frenkel et al. (2017) D. Frenkel, K. J. Schrenk, and S. Martiniani, _Monte carlo sampling for stochastic weight functions_ , Proceedings of the National Academy of Sciences 114, 6924 (2017).
  * Note (1) Note1, $\mathbf{d}(\chi,\chi^{\prime})$ equals $\mathbf{r}^{\prime}-\mathbf{r}+\mathop{arg}\nolimits\mathop{min}\displaylimits\|\mathbf{r}^{\prime}-\mathbf{r}+\mathbf{m}\|$, where the components $m_{i}$ of argument $\mathbf{m}$ run over the multiples of the simulation supercell periods.
  * (30) See Supplemental Material at [URL will be inserted by publisher].
  * Caflisch (1998) R. E. Caflisch, _Monte Carlo and quasi-Monte Carlo methods_ , Acta Numerica 7, 1–49 (1998).
  * Athènes et al. (1997) M. Athènes, P. Bellon, and G. Martin, _Identification of novel diffusion cycles in B2 ordered phases by Monte Carlo simulation_ , Philosophical Magazine A 76, 565 (1997).
  * Puchala et al. (2010) B. Puchala, M. L. Falk, and K. Garikipati, _An energy basin finding algorithm for kinetic monte carlo acceleration_ , The Journal of Chemical Physics 132, 134104 (2010).
  * Athènes and Bulatov (2014) M. Athènes and V. V. Bulatov, _Path factorization approach to stochastic simulations_ , Physical Review Letters 113, 230601 (2014).
  * Athènes et al. (2019) M. Athènes, S. Kaur, G. Adjanor, T. Vanacker, and T. Jourdan, _Elastodiffusion and cluster mobilities using kinetic monte carlo simulations: Fast first-passage algorithms for reversible diffusion processes_ , Phys. Rev. Materials 3, 103802 (2019).
  * Sharpe and Wales (2020) D. J. Sharpe and D. J. Wales, _Efficient and exact sampling of transition path ensembles on markovian networks_ , The Journal of Chemical Physics 153, 024121 (2020).
  * Wales (2009) D. J. Wales, _Calculating rate constants and committor probabilities for transition networks by graph transformation_ , The Journal of Chemical Physics 130, 204111 (2009).
  * Novotny (1995) M. A. Novotny, _Monte carlo algorithms with absorbing markov chains: Fast local algorithms for slow dynamics_ , Phys. Rev. Lett. 74, 1 (1995).
  * Opplestrup et al. (2006) T. Opplestrup, V. V. Bulatov, G. H. Gilmer, M. H. Kalos, and B. Sadigh, _First-passage monte carlo algorithm: Diffusion without all the hops_ , Phys. Rev. Lett. 97, 230602 (2006).
  * Oppelstrup et al. (2009) T. Oppelstrup, V. V. Bulatov, A. Donev, M. H. Kalos, G. H. Gilmer, and B. Sadigh, _First-passage kinetic monte carlo method_ , Phys. Rev. E 80, 066701 (2009).
  * Donev et al. (2010) A. Donev, V. V. Bulatov, T. Oppelstrup, G. H. Gilmer, B. Sadigh, and M. H. Kalos, _A first-passage kinetic monte carlo algorithm for complex diffusion–reaction systems_ , Journal of Computational Physics 229, 3214 (2010).
  * Allnatt and Lidiard (1993) A. R. Allnatt and A. B. Lidiard, _Atomic Transport in Solids_ (Cambridge University Press, 1993).
  * Barbe and Nastar (2006) V. Barbe and M. Nastar, _A self-consistent mean field calculation of the phenomenological coefficients in a multicomponent alloy with high jump frequency ratios_ , Philosophical Magazine 86, 1513 (2006).
  * Hänggi et al. (1990) P. Hänggi, P. Talkner, and M. Borkovec, _Reaction-rate theory: fifty years after kramers_ , Rev. Mod. Phys. 62, 251 (1990).
  * Settino et al. (2020) J. Settino, N. W. Talarico, F. Cosco, F. Plastina, S. Maniscalco, and N. Lo Gullo, _Emergence of anomalous dynamics from the underlying singular continuous spectrum in interacting many-body systems_ , Phys. Rev. B 101, 144303 (2020).
  * Ziff (1992) R. M. Ziff, _Spanning probability in 2d percolation_ , Phys. Rev. Lett. 69, 2670 (1992).
  * Tsai et al. (2013) K.-Y. Tsai, M.-H. Tsai, and J.-W. Yeh, _Sluggish diffusion in Co–Cr–Fe–Mn–Ni high-entropy alloys_ , Acta Materialia 61, 4887 (2013).
  * Paul (2017) A. Paul, _Comments on “Sluggish diffusion in Co–Cr–Fe–Mn–Ni high-entropy alloys” by K.Y. Tsai, M.H. Tsai and J.W. Yeh, Acta Materialia 61 (2013) 4887–4897_ , Scripta Materialia 135, 153 (2017).
  * Miracle and Senkov (2017) D. Miracle and O. Senkov, _A critical review of high entropy alloys and related concepts_ , Acta Materialia 122, 448 (2017).
  * Osetsky et al. (2018) Y. N. Osetsky, L. K. Béland, A. V. Barashev, and Y. Zhang, _On the existence and origin of sluggish diffusion in chemically disordered concentrated alloys_ , Current Opinion in Solid State and Materials Science 22, 65 (2018).
  * Kottke et al. (2020) J. Kottke, D. Utt, M. Laurent-Brocq, A. Fareed, D. Gaertner, L. Perrière, Łukasz Rogal, A. Stukowski, K. Albe, S. V. Divinski, et al., _Experimental and theoretical study of tracer diffusion in a series of (CoCrFeMn) 100-xNix alloys_, Acta Materialia 194, 236 (2020).
  * Daw and Chandross (2021) M. S. Daw and M. Chandross, _Sluggish diffusion in random equimolar fcc alloys_ , Phys. Rev. Materials 5, 043603 (2021).
  * Swinburne and Perez (2018) T. D. Swinburne and D. Perez, _Self-optimized construction of transition rate matrices from accelerated atomistic simulations with bayesian uncertainty quantification_ , Phys. Rev. Materials 2, 053802 (2018).
  * Lapointe et al. (2020) C. Lapointe, T. D. Swinburne, L. Thiry, S. Mallat, L. Proville, C. S. Becquart, and M.-C. Marinica, _Machine learning surrogate models for prediction of point defect vibrational entropy_ , Phys. Rev. Materials 4, 063802 (2020).
  * Li and Ando (2020) W. Li and Y. Ando, _Effect of local structural disorder on lithium diffusion behavior in amorphous silicon_ , Phys. Rev. Materials 4, 045602 (2020).
  * Grabowski et al. (2018) M. Grabowski, J. Rogal, and R. Drautz, _Kinetic monte carlo simulations of vacancy diffusion in nondilute ni- $x$ ($x$ = re, w, ta) alloys_, Phys. Rev. Materials 2, 123403 (2018).
  * Landeiro Dos Reis et al. (2020) M. Landeiro Dos Reis, L. Proville, M.-C. Marinica, and M. Sauzay, _Atomic scale simulations for the diffusion-assisted crossing of dislocation anchored by vacancy clusters_ , Phys. Rev. Materials 4, 103603 (2020).
  * Rahman et al. (2021) M. M. Rahman, F. El-Mellouhi, O. Bouhali, C. S. Becquart, and N. Mousseau, _Pressure effect on diffusion of carbon at the $85.{91}^{\circ}\phantom{\rule{1.60004pt}{0.0pt}}\langle 100\rangle$ symmetric tilt grain boundary of $\alpha$-iron_, Phys. Rev. Materials 5, 043605 (2021).