Automated calculation and convergence of defect transport tensors

Defect transport is a key process in materials science and catalysis, but as migration mechanisms are often too complex to enumerate a priori, calculation of transport tensors typically have no measure of convergence and require significant end user intervention. These two bottlenecks prevent high-throughput implementations essential to propagate model-form uncertainty from interatomic interactions to predictive simulations. In order to address these issues, we extend a massively parallel accelerated sampling scheme, autonomously controlled by Bayesian estimators of statewise sampling completeness, to build atomistic kinetic Monte Carlo models on a state space irreducible under exchange and space group symmetries. Focusing on isolated defects, we derive analytic expressions for defect transport tensors and provide a convergence metric by calculating the Kullback-Leiber divergence across the ensemble of diffusion processes consistent with the sampling uncertainty. The autonomy and efficacy of the method is demonstrated on surface trimers in tungsten and hexa-interstitials in magnesium oxide, both of which exhibit complex, correlated migration mechanisms.


Introduction
The migration and transformation of intrinsic and extrinsic complex crystal defects plays a central role in numerous materials science and chemistry phenomena such as post-irradiation annealing 1 , plasma surface interactions 2 , active site formation for heterogeneous catalysis 3,4 or mechanical properties of concentrated solid solution alloys 5 .
The atomistic mechanisms available to nanoscale defects are highly heterogeneous with defect size and impossible to divine a priori due to the routine presence of complex multiatom transformations [6][7][8][9][10] . Whilst these effects typically becomes less volatile with increasing defect size, phenomenological or higher-scale models can only capture the defect dynamics if suitable parameters can be calculated [11][12][13] . In addition, nanoscale clusters are typically the most mobile and thus often have a much greater influence on macroscopic phenomena than slower-moving larger defects. In the general case, atomistic mechanisms must be discovered through unbiased dynamic 9,14-18 or static 10,19,20 sampling approaches. When the true dynamics can be characterized as rare transitions between metastable basins on the energy landscape 21 , the basin-to-basin dynamics can be mapped to a continuous time Markov chain 22,23 , which forms the theoretical foundation of atomistic kinetic Monte Carlo (akMC) methods 24 , of which on-lattice kMC is a subclass. The resulting model can then be stochastically or in some cases (such as that presented here) analytically integrated to extract observables of interest.
A well recognized problem is that an akMC model will in general have an incomplete catalogue of available mechanisms due to finite amount of sampling, and this can produce catastrophically erroneous predictions if important mechanisms are omitted 17,18,25,26 . Sampling adequacy is often assesed qualitatively using the domain expertise of simulation practicioners. Whilst this approach has undoubtably yeilded significant successes, it requires significant end user analysis for each system under study. This has a punitive impact on the feasibility of automating complex materials simulations on massively parallel computational resources, where the required decision frequency rapidly exceeds practical human limits. Further, in absence of quantitative uncertainty quantification approaches, assessing the reliability of meso or macro-scale predictions is extremely challenging.
For example, the system sizes required even for small defect clusters render ab initio calculation unfeasible. As a result, interatomic potential models must be used which induces an additional model form uncertainty. The development of interatomic potentials has been revolutionized in recent years through the use of linear-in-descriptor or neural network regression techniques [27][28][29] . These approaches offer a natural encoding of model form uncertainty through isosurfaces of the cost function used for potential parametrization. Highthroughput calculations are essential in this context to enable the systematic propagation of this uncertainty on interatomic interactions to the observables of scientific interest, without prohibitive hours of end-user analysis.
In this contribution we present an autonomous, highly scalable sampling scheme to efficiently calculate defect diffusion tensors with quantified uncertainty on the sampling completeness. We demonstrate the ability to discover complex migration behaviours of defects in tungsten and magnesium oxide, and show how the quantified uncertainty can be used to rapidly yeild well defined convergence measures. Our approach enables high-throughput workflows to rapidly discover, converge and analyze complex kinetic properties of defect structures in surrogate energy landscapes, with minimal end user involvement.

Results
In previous work 26 we introduced TAMMBER, a massively parallel accelerated sampling scheme whose formal objective is to exhaustively sample the set of all metastable minima M and all interminima transitions T for a given system. TAMMBER builds a matrix K ji of j ← i rates in the known state space and a vector k u i which estimates the 'unknown' (as yet undiscovered) escape rate from each state i, encoded as transitions to an absorbing sink 17,18,26,30,31 (methods). States are identified by constructing a connectivity graph from a minimized configuration, which is then hashed to produce a pseudounique integer label for each configuration 10,32 . The presence of nonzero k u causes trajectories to leave the known state space, giving a natual measure of model validity, the expected residence time 26 in the known state space where K tot ji = δ ji [k u i + l K li ], P 0 is the initial distribution and 1 is a row vector of ones. In practice, truly comprehensive sampling of M and T is often impractical as their size grows exponentially with the number of atoms in the system. However, M and T are both highly reducible under exchanges of indistinguishable atoms R, lattice vector translations S and space group symmetries G, especially in the case of quasi-zerodimensional defects. It is known that configurations which are degenerate (in the periodic minimum image sense) under these operations will have isomorphic connectivity graphs 10,32 , with isomorphisms that can be efficiently calculated 33 . We exploit this reducibility to make a partition of M, T into equally sized subsets where every member in a partition M p is identical to all other members within a reindexing in R and translation by t ∈ S. Similarly, the subset of transition T p0q ⊂ T pq that take a configuration p 0 ∈ M p to members in M q is identical to any other subset T piq under application of R and S. To exploit space group symmetries, we first note that every defect structure has a symmetry group G p ⊆ G, a subgroup of G containing at least the identity operation, whose order #G p (number of elements) factorizes the order #G of G. When acting on every member of M p with some space group symmetry G ∈ G to produce a new set of configurations GM p , there are two possible outcomes: In the latter case we have acheived a further reduction in the number of partitions we need to consider, as M q can be generated from M p through the action of G. However, this only defined up to postmultiplication by G p ∈ G p , as GG p M p = GM p . To resolve this ambiguity we generate all possible (left) cosets 34 with respect to G p , defined for some G ∈ G as {GG p : G p ∈ G p }, i.e. the set of operations formed by postmultplying G by every element in G p 35 . We can always form a setḠ p of exactly #G/#G p nonequivalent cosets of size #G p , where each member of G appears in exactly one coset 34 . Clearly, the elements of G p form one of these cosets, whilst acting on M p with any element from a given coset will give the same result. We can then index the #Ḡ p nonequivalent partitions {M q } which can be generated from M p by the coset which generates them; drawing some operation G qp from each coset we can thus make the further compression M = P G p=1 q G qp M p , where P G ≤ P , with analagous partitioning of T pq .
We therefore only need to sample the exit transitions from a maximum of P G configurations, as all others can be generated by the known symmetric relations. By focusing sampling resources on a much smaller number of possible states, the unknown rates decrease at a much faster rate, with a subsequent increase in the validity timescales for the constructed Markov chain. As in TAMMBER the unknown rates at low temperatures decrease as an inverse power of the MD time invested at high temperatures, state-space compression can yield extremely large benefits in practice. This maximally compressed representation is thus used for sampling, which is then partially decompressed into the partitioning (2) of P distinct states irreducible under translation for the calculation of transport tensors, as illustrated in figure (1). Whilst connectivity graph isomorphisms have been used to identitify already seen defect structures 10 , the production of Markov chains with complete or partial irreducibility for sampling or model construction is novel to the best of our knowledge. As an example of the efficiencies this approach affords, in previous work we investigated the breakup of two interstitial defects in bcc iron 26 , obtaining a Markov model on the uncompressed (reducible) state space with τ res ∼ 80s at 300K, an insufficent duration to make confident predictions on the breakup mechanism. Despite the relatively low symmetry of this defect system, we show in the supplementary material that the compressed sampling scheme described above yields a Markov model with τ res ∼ 5 × 10 6 s with only 75% of the computational effort, allowing convergence in the model predictions. The isolated point defects we consider in this work are typically of even higher symmetry, giving correspondingly greater efficiencies.

Evaluation of transport tensors
To define the drift and diffusion tensors, we first require a defect position in the supercell for each configuration. To remove ambiguities from periodic boundary conditions we carve out a defective region of some configuration in each partition by thresholding some structural descriptor, here the centrosymmetry, then take a descriptor-weighted center of mass. As we allocate sampling to increase the residence time, the threshold value can be freely determined in post-processing; typically a single value is suitable for the same material system. It is also possible to determine the defect position through analysis of a sufficiently large set of mutally isomorphic configurations, but we find the descriptor-based approach to be simple, efficient and reliable.
A periodic minimum image displacement vector d lm for any l ∈ M p → m ∈ M q transition can then be found through application of the known rigid transformation from the sampled configurations. Importantly, transitions between states within the same space group partition are then represented as self transitions l ∈ M p → m ∈ M p (closed loops in figure 1).
We can then construct an irreducible group of states from the space group partitioning of equation (2), with one state indexed by p ∈ [1, P ] for each set M p . All the transitions from a state p to states q in the uncompressed model are mapped to their compressed states, forming a (possibly empty) set C qp = {k l , d l } of transitions and displacement vectors. The compressed rate matrix is then defined as K c qp = l∈Cqp k l , with a statewise total escape rate ma- where the superscript c indicates we consider rates in the compressed model. The unknown rates k u,c p are taken from the maximally compressed representation, with the same unknown rate for symmetrically equivalent partitions. One could then generate akMC trajectories 37 with a matrix and vector of branching probabil- δτ l before absorbtion. In the supplementary material we show that the drift and diffusion coefficients µ, D can then be extracted from the relations x = τ res µ and x ⊗ x = 2 τ res D + τ 2 res µ ⊗ µ. However, it is also possible to analytically evaluate averages over all possible pathways using a 'displacement generating function' Moments of the total displacement can then be written x = −∂ λ Z(λ, P 0 )| λ=0 and x⊗x = ∂ λ ⊗∂ λ Z(λ, P 0 )| λ=0 (supplementary material). Whilst (3) could be used for any choice of initial condition, we note that in the well sampled limit k u,c p r K c rp , which is necessary but not sufficient for global convergence, the matrix [K tot,c − K c ] will have a spectral gap, with one eigenvalue 0 < ν 0 ν 1 < ν 2 ... much smaller than all others 22 . The right eigenvector for ν 0 is the quasistationary distribution (QSD) π QSD in the known state space 22 , the limiting distribution conditional on not absorbing for an arbitrary long time, which as k u,c → 0 becomes the Boltzmann distribution, π QSD →π. As the QSD is the longest-lived mode and transport coefficients are defined as the limit of infinitely long trajectories, it is natural to set P 0 = π QSD = π QSD /(1π QSD ) to eliminate the influence of initial conditions. With this choice, it is simple to show that τ n res = n!/ν n 0 and the expected drift and diffusion coefficients emerge as (supplementary material) Equations (4,5) are a central result of this contribution, expressions for the drift and diffusion tensors of an arbitrarily complex diffusion process in a periodic system, autonomously constructed in a massively parallel sampling scheme, which crucially are dependent on 'unknown' rates k u,c that robustly quantify sampling incompleteness. In the supplementary material we show the limiting expressions µ ≡ lim k u,c →0 µ(k u,c ) and D ≡ lim k u,c →0 D(k u,c ) reduce to expressions obtained in previous derivations on multistate diffusion in periodic media 38,39 , with all uncorrelated and correlated contributions that are essential to capture complex diffusion pathways. In the remainder we focus on systems obeying detailed balance, where µ = 0 and the correlated contribution to the diffusivity is always nonpositive.

Convergence of the diffusivity
The central novelty of (5) is that D(k u,c ) estimates the diffusion tensor over all possible trajectories in the known state space before exit due to the unknown rates. Possible changes to the diffusivity under the discovery of additional transition rates between known states can then be bounded; the rate matrix K c can be modified by an additional rates matrix δK c which must satisfy detailed balance and not increase the total exit rate from each known state p by more than the unknown rate k u,c p , meaning [1δK c ] p ≤ k u,c p . We have desiged a Monte Carlo procedure to sample the space of permissible δK c (see supplementary material) which typically requires less than a core-minute for the systems studied here and is trivally parallelizable. We discuss the sensitivity to the discovery of additional states at the end of this section.
To analyze the ensemble of diffusion tensors produced in the Monte Carlo procedure, we diagonalize the 3 × 3 diffusion matrix, producing eigenvalues {D l }, l ∈ [1,3]. The Monte Carlo procedure yeilds upper and lower bounds D l± for each of the eigenvalues from which a convergence metric can be obtained; in practice, we find these bounds to be highly asymmetric, with D l− → D l when we constrain the Monte Carlo procedure to consider only the addition of new transition rates without any additional state. In addition, the eigenbasis undergoes negligible changes close to convergence, meaning that to a high degree of accuracy the matrices D ± (with eigenvalues D l± ) can be simultaneously diagonalized. To produce a dimensionless convergence measure, consider the fundamental solutions ρ(x, t), ρ ± (x, t) to the three dimensional diffusion equation with a diffusion tensor D, D ± . A natural measure is the Kullback-Leibler divergence 40 R ± ≡ dxρ ln ρ/ρ ± , which for diffusion is time independent, reading R ± = Tr D −1 D ± /2 − 3/2 + ln D −1 D ± . Our convergence measure is the spread δR ≡ R + − R − ≥ 0, which when D ± and D can be simultaneously diagonalized reads where the limit applies close to convergence, being the leading order expansion in D l± − D l . As δR = 0 for D + = D − , we have a well defined dimensionless convergence metric ideal for autonomous implementation, which has an informative limit, namely the relative spread in eigenvalues consistent with the sampling uncertainty. It is clear that the diffusivity cannot be globally bounded against the discovery of some new set of states which are free to possess arbitrary transport properties. However, as transition rates to such a set of states are bounded by the k u , the Monte Carlo prodedure outined above does characterizes transport behaviour in the known state space over timescales of order τ res . Convergence to the large τ res limit can be accelerated by seeding the TAMMBER procedure with as many states as possible, which are free to be completely disconnected 41 . The integration of automated structural search algorithms 21,42 into the present workflow will be the subject of a future work. Of course, as TAMMBER generates thousands of high temperature molecular dynamics trajectories across the known state space it is an automated global minimum search method, with τ res acting as quantitative measure of sampling quality when the system is ergodic. As a result, whilst in common with all theoretical studies on high dimensional landscapes we cannot provide bounds on global minima, we can provide a key uncertainty quantification on the validity of our findings, namely a rigorous prediction timescale from an arbitrary distribution on the known state space with a corresponding bound on transport coefficients. We emphasize that the inability to assign bounds on global minima searches applies equally to human guided sampling or the present approach. The convergence metrics we provide are thus a valuable analytical tool which removes ambiguities inherent to traditional methods in addition to enabling a fully automated workflow suitable for high-throughput computation. Application to trimer diffusion on W(110) The diffusion of adatom clusters is a fundamental process in surface science, and has recently been conjectured to play a crucial role in the formation of complex 'fuzzy' surface morphologies during plasma exposure in nuclear fusion reactors 43 . We focus here on the trimer defect a demonstrative case study; comprehensive high-throughput investigations for which the present approach is designed will be presented elsewhere. In the present case TAMMBER was initialized with state A 1 in figure  1, then run for 8 hours on 144 cores using the EAM4 embedded atom method potential by Marinica et al. 44 , covering a temperature range of 400K-1400K. The corresponding values of τ res as a function of temperature are also shown in Fig.  1.
The resultant diffusion behaviour is highly correlated, with many A 1 ↔ A 2 and D 1 ↔ D 2 transitions in particular. However, the overall system does not exhibit a clear 'superbasinto-superbasin' diffusion mechanism, meaning access to the full highly correlated trajectory ensemble is essential to extract accurate transport coefficients. As shown in Fig. 1 To look for dominant pathways, we form a weighted graph from the connectivity of four primitive unit cells, where the graph edges are weighted by the corrsponding saddle point energy 21,26 . Dijstra's shortest path algorithm 45 was then used to identify the dominant pathways. In agreement with the found Arrhenius slopes, migration along v 2 is dominated by A 1 → D 1 → C 1 → D 2 → A 5 paths, with a well defined activation energy of ∆E 1 = 1.32eV at the lower temperatures. Migation along v 1 is similarly dominated by A 1 → A 4 paths at lower temperatures, but A 1 → A 3 paths have a growing contribution with temperature, giving a weak nonlinearity to the Arrhenius gradient ∆E(β) = −∂ β ln |D|. An in depth study of this procedure, and its role in a fully automated workflow, will be the subject of future work. We note that a recent study 43 of surface island diffusion on W (110) at 1000K using a different interatomic potential 46 reported trimer migration via A 1 → A 3 , which is accounted for in the present study but is not found to be the dominant mechanism. The ability to efficiently resolve such ambiguities, without prohibitive person hours, is a key advantage of methodology presented here.

Application to hexa-interstitial in bulk MgO
To conclude this contribution we investigate the diffusion of a stochiometric hexainterstitial in MgO, whose lattice has space group F m3m with point group O h . In a cubic supercell with axes aligned with 100 directions we thus retain the full symmetry of the host lattice before the introduction of any defective structures. Connectivity graphs are constructed with vertices colored to indicate specie. Due to the high degree of symmetry, the maximally reduced state space contained only four states, whilst the state space irreducible under translation contained 56 states. TAMMBER was initialized in a relatively high energy state (red circle) that was found by Uberuaga et al. 7 to possess very low migration barriers; a much lower energy state was rapidly found, upon which sampling was subsequently concentrated, demonstrating the ability of TAMMBER to act as a massively parallel global minimum search routine. Using the same modified Buckingham potential as in that work, TAMMBER was run for 8 hours on 144 cores targeting a temperature range of 500-1000K. Whilst we find the very low migration barriers (¡0.1eV) in agreement with 7 , the full diffusion tensor was found to be isotropic D 1 = D 2 = D 3 , with a characteristic activation energy that converged to 0.887eV at the lower temperatures, which slowly increases at higher temperatures, as found for the above trimer example. Analyzing the dominant pathways with temperature reveals that this activation barrier corresponds to self-migration of the lowest energy state, with little effect of correlation on the diffusion paths; analyzing the state-by-state contribution to the uncorrelated diffusivity, shown in Fig. 2, we see that at higher temperatures the self-migration of other states have an increasingly large contribution. τ res as a function of temperature are also shown in Fig. 2. The Monte Carlo procedure indicated a high degree of convergence, with δR < 0.001 over all temperatures considered.

Discussion
We have presented a fully automatable and efficient method to evaluate the transport tensors resulting from the arbitrarily complex diffusion processes of crystal defects, with a well defined convergence criteria based on quantitative measures of sampling uncertainty combined with a Monte Carlo procedure to sample the admissable diffusion tensors consistent with sampling uncertainty and detailed balance. The method was demonstrated on a surface trimer in in tungsten and a hex-interstitial in magnesium oxide. By effectively eliminating user input beyond the seeding of some initial state(s), the presented approach demonstrates sufficient computational and critically end user efficiencies to extend the phenomenological reach of high throughput computations to point defect kinetics. Future work will exploit these efficiencies to analyze defect transport over a wide range of material systems, and the influence of breaking detailed balance through external driving forces, giving nonzero limits for the drift vector µ. Acknowledgements TDS gratefully recognizes support from the Agence Nationale de Recherche, via the MEMOPAS project ANR-19-CE46-0006-1. This work was granted access to the HPC resources of IDRIS under the allocation AP010910718 attributed by GENCI. Work at Los Alamos National Laboratory was supported by the U. S. Department of Energy, Office of Nuclear Energy and Office of Science, Office of Advanced Scientific Computing Research through the Scientific Discovery through Advanced Computing (SciDAC) project on Fission Gas Behavior. Los Alamos National Laboratory is operated by Triad National Security LLC, for the National Nuclear Security administration of the U.S. DOE under Contract No. 89233218CNA0000001.

Competing Interests
The Authors declare no Competing Financial or Non-Financial Interests. Contributions TDS and DP designed the research program. TDS derived the theoretical results, designed the Monte Carlo procedure, and ran the simulations. TDS produced an initial manuscript, which was then discussed and refined with DP.

Data Availability
The datasets generated during and analysed during the current study are available from the corresponding author on reasonable request. Methods Sampling procedure Through the use of NEB calculations 47 and transition state theory 48 TAMMBER constructs a transition matrix K of rank equal to the number of discovered states, giving a continuous time Markov chainṖ(t) = K − K tot P(t). Exploiting the known Poissionian distribution of exit times from a suitably thermalized basin 22 , a Bayesian likelihood (and thus posterior distribution) for the {k u i } was derived using parallel trajectory data obtained through a modified temperature accelerated dynamics method 16 . Continuous self-optimization was acheived by calculating the derivative of each k u i with respect to additional computational work, allowing the degree of temperature acceleration to be statewise optimized and the massively parallel sampling effort distributed across states to differentially maximize a key measure of model validity, the expected residence time before absorbtion 26 , equation (1). Typically, a simulation starts with one state then the rank of P, K and K tot increases as states are discovered. If P0 is fixed to be unity for the initial state and zero otherwise, i.e.
[P0]j = δij for an initial state i, τres is monotonically increasing with sampling effort, a key consequence of the estimation procedure for {k u i }. An implementation of this method is available as an open source code 49 , whose near-ideal parallel efficiency has been demonstrated on massively parallel resources employing 1000 to more than 80,000 cores.

Isomorphically compressed representation
Isomorphic configurations are identified in two ways-the connectivity graph for every state is duplicated and reindexed into McKay's pseudounique canonical order 33 , which is identical for isomorphic states. Alternatively, the VF2 graph matching routine 50 is applied to a single state to find self-isomorphisms. To determine a given reindexing Rij between two isomorphic states i, j we first find the mappings Ric and Rjc to the McKay canonical order then obtain Rij = RicR −1 cj . As isomorphisms will be with respect to the simulation supercell, not the host crystal structure, the relevant point group is G = W ∩ G ∈ O h , the intersection of the point group of G and the supercell point group W, where W is a subgroup of the cubic group O h . Unlike an arbitrary element of G or W, any element ofḠ is guaranteed to leave both the perfect lattice and supercell unchanged up to a translation and reindexing. We then iterate through all 48 members of O h , applying the point transform and applying a constant displacement such that the first indexed atom from each configuration are minimum image coincident. We then check for minimium image coincidence atom-by-atom, rejecting each candidate member of O h at the first failure. Consider the set Pn of paths ξ = {ξi} n 0 which execute n transitions between n + 1 states before absorbtion. The probability weight P ξ of each path and all n-paths is given by Below: Probability P B<∆ that the system executes a A → B transition before absorbtion, residence times τres and various approximate projected breakup times 26 at a range of temperatures using a Markov model built using the previously presented reducible (left) or new irreducible (right) representation. We see the irreducible representation achieves a significant compression of phase space and yields much greater certainty, as encoded by P B<∆ , that our model has sufficient predictive power to capture the dominant breakup mechanism.
The sum over of path weights over the set of all possible paths P = ∪nPn can be confirmed to be unity as expected- The total displacement x ξ = l d ξ l ξ l−1 across each path can be found by simply constructing the pathwise probability as before, but replacing B(0) with B(λ) to yield The expected value (or first moment) of the total displacement over the set of all possible paths is therefore The second moment matrix of the total displacement follows naturally- where ⊗ is the outer (or dyadic) product.

Appendix C: Evaluating derivatives of the generating function
Consider the following application of the chain rule, using the above definition G(λ) = [I − B(λ)] −1 , but writing G λ instead of G(λ) for readability- . The second derivative of G(λ) follows by induction as Using the above identity We now evaluate the vector and matrix-valued first and second derivatives of the components of B λ , writing Appendix D: Transport coefficients from trajectories of variable duration As discussed in the main text and below, in the well sampled limit the matrix K tot − K has a spectral gap, with one eigenvalue ν0 > 0 much smaller than all the others. The corresponding right eigenvector is the quasistationary distribution 22 π QSD , which as ν0 → 0 becomes the Boltzmann distributionπ in equilibrium. The QSD is a natural choice of initial distribution, as it is the limiting distribution in the known state space conditional on not absorbing. As a result, setting P(0) =π QSD effectively eliminates initial transients in the extraction of transport coefficients. With this choice of initial condition the residence time distributionṖ(τ ) = 1Ṗ(τ ) is a single exponential with decay constant ν0, i.e.Ṗ(τ ) = ν0 exp(−ν0τ ). Consider an average . . . τ over the set of all trajectories Pτ with τres ∈ [τ, τ + dτ ]. It is clear that . . . = ∞ 0Ṗ (τ ) . . . τ dτ . The drift and diffusion estimators µτ and Dτ over such trajectories are defined through x ⊗ x τ = ξ∈Pτ x ξ ⊗ x ξ P ξ = 2τ Dτ + τ 2 µτ ⊗ µτ .
With increasing trajectory length, τ → ∞, the drift and diffusion coefficients are only well defined if the estimators tend to a constant value, i.e. limτ→∞ µτ = µ and limτ→∞ Dτ = D.
In the well sampled limit where K tot − K has a spectral gap, averages over the set of all possible paths P = ∪τ Pτ will be dominated by the large τ regime, where µτ and Dτ are constant. Furthermore, due the preparation of the initial conditions in the QSD (P(0) =π QSD ) the exit time distribution is the single exponential decayṖ(τ ) = ν0 exp(−ν0τ ), meaning τ = ν −1 0 and τ 2 = 2ν −2 0 = 2 τ 2 .
Exploiting the chain rule and the identity Z(0) = 1, the drift and diffusion constants can thus be succinctly written as quoted in the main text. We thus see that when trajectories are of variable length, the appropriate expression to extract the diffusivity in the possible presence of drift is not simply the variance in the total displacement, but must account for the properties of the distribution of trajectory durations.
As every jump vector d l , l ∈ δCpq will have a one-to-one correspondence with a vector d l = −d l for some l ∈ δCqp, the detailed balance symmetries imply in particular that These relations imply that the drift vector is always zero in equilibrium-using π = Π1 and the fact that 1M1 = 1M 1 for any matrix M we have µ = 1K ⊗ dΠ1 = 1 K ⊗ dΠ 1 = −1K ⊗ dΠ1 = −µ, ⇒ µ = 0.
We now split D into uncorrelated and correlated terms as As all rates are positive, it is clear that Tr Duc > 0. For the correlated contribution we first use the symmetry of W + to write where we introduce W to write the symmetric matrix Π −1 W + as W W. We then use detailed balance symmetries to rewrite the correlated contribution as which shows that the correlated contribution always reduces the diffusivity for a system obeying detailed balance.