Weighted-ensemble Brownian dynamics simulation: Sampling of rare events in non-equilibrium systems

We provide an algorithm based on weighted-ensemble (WE) methods, to accurately sample systems at steady state. Applying our method to different one- and two-dimensional models, we succeed to calculate steady state probabilities of order $10^{-300}$ and reproduce Arrhenius law for rates of order $10^{-280}$. Special attention is payed to the simulation of non-potential systems where no detailed balance assumption exists. For this large class of stochastic systems, the stationary probability distribution density is often unknown and cannot be used as preknowledge during the simulation. We compare the algorithms efficiency with standard Brownian dynamics simulations and other WE methods.

Rare events are ubiquitous in many biological, chemical and physical processes [1,2]. Whereas the density of states is known in systems at thermal equilibrium, interesting phenomena often occur in non-equilibrium systems [3]. Unfortunately, many such problems are inaccessible to analytic methods. Therefore computer simulations are a widely used tool to estimate the density of states or transition rates between them [4,5]. Since standard Brownian dynamic simulation [6,7] provides computational costs that are inversely proportional to the state's probability, specialized methods [8][9][10] have to be used to adequately sample rare events, i.e. states with low probability or low transition rates.
In the last decades, flat histogram algorithms [11] have been developed, allowing one to evenly sample states with highly different probabilities. These algorithms are implementations of the umbrella sampling [12], where each state is sampled according to a given probability distribution, the so-called umbrella distribution. Within nonequilibrium umbrella sampling (NEUS) [13] the space of interest is divided into different but almost evenly sampled subregions. The interaction between different regions occurs solely due to probability currents between then. Whereby the probability distribution within a region is then calculated by performing Monte Carlo simulations.
In order to calculate low rates between a starting and a final state, forward flux sampling methods can be used (for a review see) [14]. These methods introduce a sequence of surfaces between these states and introduce walkers (copies of the system) to perform weighted trajectories according to the underlying dynamics. If walkers cross one of the surfaces, getting closer to the final state, new walkers with smaller weights are introduced. Finally, many walkers with particular small weights reach the final state. The consideration of the particular weights allows one to calculate very low rates in a finite simu- * justuskr@physik.hu-berlin.de lation time. Recently, extensions to these methods have been developed to calculate both, transition rates using umbrella sampling [15] and probability distributions using forward flux sampling [16] algorithms.
In this work, we present an algorithm, based on the previously developed weighted-ensemble (WE) Brownian dynamics simulations [17][18][19][20], that allows one to calculate the stationary probability density function (SPDF) as well as transition rates between particular states. Like in WE simulations the space of interest is divided into several subregions and the probability for finding the system in them is calculated by generating equally weighted walkers in each region. By moving to the underlying dynamics, the walkers transport probability between the subregions. Thus, WE methods are usually applied to systems of Brownian particles moving in a potential landscape [18,21].
We are interested in an algorithm which allows simulations of stochastic dynamical systems which, apriori do not obey detailed balance for probability fluxes or suppose some special topology of the flow [22][23][24][25][26]. Such are given for Brownian particles in conservative force fields under the influence of additive noise. Even canonic dissipative systems possess a vanishing probability flow if transformed to the energy as dynamic variable [24,27,28]. In consequence, both systems allow exact analytic solutions of the SPDF. In contrary, we aim to develop an algorithm which does not assume that neither the deterministic nor the stochastic items (see Eq. (1) below) underly such conditions. Thus, no information on the SPDF can be used for the simulations.
In general the algorithm can be applied to arbitrary dynamical systems of the form: x n = f n (x) + g n (x)ξ n (t), n = 1, ..., d, where d is the number of stochastic time-dependent degrees of freedom x n (t), n = 1, ..., d; the functions f n (x) describe the deterministic velocities for the n-th direction; ξ n (t) represents zero-mean Gaussian noise with delta-like correlation function ξ n (t)ξ m (t ′ ) = δ nm δ(t − t ′ ). The noise intensity along the n-th direction is scaled by the functions g n (x), which in general depend on the vector x = (x 1 , x 2 , ..., x d ) T . We are interested in high precision sampling of the stationary probability current J st (x) and the SPDF p st (x) of finding the system in the d-dimensional cube [ with a finite resolution. We will specify the resolution by the number M n,res of evenly spaced supporting points along the n-th direction, for which we will determine p st (x). This article is organized as follows: In section I we introduce an algorithm, based on WE methods [17], that allows one to calculate low probabilities and rates. Afterwards, we study one-and two dimensional model systems and analyze the algorithms efficiency compared to Brownian dynamics simulation (BDS) and WE techniques.

I. THE ALGORITHM
First, for sake of clarity, we restrict ourselves here to particle motion in one dimension, d = 1, under additive noise. The deterministic part of the dynamics Eq. (1) can be always represented as a conservative force f (x) = −U ′ (x). The noise strength is scaled by the parameter D, we put g(x) = √ 2D. For such systems, the SPDF is known to be where Z st is a normalization constant. We are interested in finding numerically the system's SPDF, p st (X j ), at a set of M res evenly spaced supporting points X j in a finite part of the physical space, given by Supporting points are given explicitly by X j = L − + (j − 1 2 )∆X res , j = 1, 2, . . . , M res , with ∆X res = ⌈ M Mres ⌉∆x, see Fig. 1. Here ⌈z⌉ denotes the largest integer smaller than or equal to z.
Let us introduce the probability P i (t) for finding a particle in the i-th subregion at time t, and the corresponding set P(t) = (P 0 (t), P 1 (t), ..., P M−1 (t)). Initially, no information on the system is available, thus, each subregion is given an arbitrary amount of probability P i (0), simply fulfilling the normalization condition M−1 i=0 P i (0) = 1. Naturally, equilibration can be accelerated if one already has information on the SPDF of the system (Eq. (1)).
In that case, one can choose P(0) close to the set P st , which optimally approximates the SPDF: However, in general no such information is required.

A. Time evolution
After setting the initial set P(t = 0), the time evolution of the P(t) → P(t + h) is performed using three different steps.
We start with a redistribution step, in which N walkers (copies of the system) are uniformly distributed in each subregion. Besides their individual positions x k i (t), where i = 0, ..., M − 1 denotes the particular subregion and k = 1, ..., N the individual walkers, each walker possesses a given amount of weight q k i (t). This is nothing but the present probability in the i-th subregion distributed to the N walkers, which yields Note that one does not need to introduce walkers in subregions with P i (t) = 0. After the redistribution step has been performed, Eq. (1) is integrated for all walkers, using a Brownian dynamic simulation step h and an arbitrary integration scheme. This integration step realizes the time evolution Here walkers transport probability between the subregions. As walkers are independent of each other, it is of importance to note that the particular time evolution of each one of the N × M walkers is due to different sample paths in the stochastic parts of the Langevin equation.
Lastly, an updating step is performed, in which the new probabilities P i (t) → P i (t+h) are calculated by summing up the weights of all walkers that are currently located in the particular subregion, In what follows, we will name the sequence of redistribution, integration and updating step as running step. After an equilibration time T therm , the set P(t) reaches a stationary regime, where the P i (t)'s fluctuate around their mean values P i .

B. Calculating the stationary probability
The individual P i are estimated by averaging over a total amount of N T sets P(t ℓ ), ℓ = 1, 2, ..., N T , taken, after the system has reached the stationary regime, every n av running steps: t ℓ = T therm + (ℓ − 1)n av h. It turns out that the mean probabilities P i coincide with the stationary probability P st,i (see Eq. (3)), for compatibly chosen time step h and the size of a subregion ∆x (see Sec. I E).
Finally the SPDF on the supporting points p st (X j ) is calculated by adding the adjacent P i and dividing by the size ∆X res (in order to have a properly normalized PDF): C. Calculation of the probability current The stationary probability current J st (x) at position x can be easily calculated by adding up (with the right sign) the weights of all walkers, passing x to the right and to the left per unit time. In practice x should be the boundary of a subregion. If x = x i the current J(x i )(t) is given by: and R i indicate these walkers which cross the boundary moving rightwards, i.e.
Averaging over N T such estimates, taken in the stationary regime, leads to the average current J(x i ) , which converges towards J st (x i ) for N T → ∞.

D. Implementation of boundary conditions
The implementation of boundary conditions for the probability current or the SPDF is straightforward. Right now, absorbing boundaries are already implemented at L − and L + , since walkers that pass these boundaries are not located in any subregion. Therefore, their weights will get lost in the next updating step. Reflecting boundary conditions at L + can be implemented by setting Hence, the probability current at L + will be zero. Reflecting boundaries at L − can be implemented analogously.

One-dimensional systems
In order to ensure, that P i and J(x i ) converge towards the stationary probability distribution and the probability current of Eq. (1) for N T → ∞, the time step h and the size of a subregion ∆x have to fulfill specific criteria. This is due to the redistribution step, where walkers are uniformly distributed in each subregion. This implicates statistical errors, since they can reach positions in a subregion that, are inaccessible or at least more improbable. Therefore, walkers can more easily escape from potential minimums or reach regions of low probability. This effectively flats the probability distribution, leading to more probability in regions of low probability, for instance, around local maximums of U (x), and less probability in the potentials minimums. In order to overcome this problem, earlier works [17,18] have stored the positions and weights of all walkers. In the next redistribution step, walkers were only spaced on the stored positions according to the weights belonging to them. This requires a lot of computer memory, especially for large N and M .
However, we found that one does not need to store these information, if the subregions are small enough to ensure that walkers have a non-negligible probability to leave them during one integration step. As a measure of how far a walker can step, due to the fluctuations, in one time step, we use the diffusion length L dif = 2 √ Dh. Thus, the size of a subregion ∆x should be small com- The distance a walker can pass during an integration step is not only determined by L dif , but also by the deterministic dynamics, leading to a step length L det = f (x)h in first order. Usually f (x) changes very fast at the boundaries of the simulation area, which produces high deterministic velocities and regions of low probability. Walkers can only reach these regions, if the fluctuation are strong enough to balance the deterministic force, i.e. if for all x ∈ [L − , L + [. This leads to a condition for the time step h: Hence, lower time steps allow one to sample regions, far from the potential extrema and for instance, the tails of the SPDF. If a larger time step is chosen, the P i will run to zero in subregions with larger deterministic force.
Since it is often difficult to fulfill Eq. (10) in the entire simulation area, one should choose a time step, which allows one to fulfill Eq. (8).

Multidimensional systems
In general, our method can be applied to stochastic dynamical systems in arbitrary dimension d ≥ 1 (Eq. (1)). For the foundation of the algorithm, we refer to the Appendix Sec. B. However, in order to ensure convergence in a finite region A, the criteria for the time step h and the size of a subregion ∆x i along the i-th direction should be fulfilled. For noise dominated directions the criteria hold, therefore, the condition for the time step (Eq. (10)) becomes: and the criteria for the size of a subregion along the n-th directions (Eq. (8)) reads: However, there might be directions without any noise (D n = 0). In that case the length of a subregion should ensure, that walkers can leave it due to the deterministic term f n (x), leading to otherwise information on the deterministic dynamics gets lost during the redistribution step, since walker that stay into a subregion do not produce any change in the P i and are again randomly placed in their subregion during the next redistribution step. For equally sized subregions ∆x n should be the minimum value of Eq. (13) with respect to all x in the simulation area for each of the d directions. These criteria can lead to a huge number of subregions, especially in high dimensional spaces.
To appropriately reduce the number of subregions, the ∆x n should be chosen in order to locally fulfill the criteria. This was implemented by a grouping algorithm, that groups original subregions into "larger" ones as long as walkers can leave these due to the deterministic term (Eq. 13). Depending on the system, this procedure highly reduces the total number of "larger" subregions M group . If the grouping algorithm was used, N walkers are randomly placed in each of these "larger" subregions and a probability P i of finding a walker in the corresponding area was introduced. We find that such grouping highly reduces the computational costs, since less subregions and therefore less walkers are required. Since walkers jump out of these regions until the next redistribution step starts, the algorithm still approximates the correct SPDF.

F. Simulation techniques
Simulations were performed on a Intel R Xeon R CPU E31245 @ 3.30GHz processor with 16 Gb DDR-3 RAM. The algorithm described above was implemented in a C ++ program for one-and two-dimensional systems. Runs of the algorithm are specified by the time step h, the size of a subregion ∆x, (∆y, in two-dimensional problems), the simulation area, given by L − and L + (L ± x , L ± y ), the number of walkers per subregion N , the thermalization time T therm , and the number of running steps between two sets of P denoted as n av . The numerical integration of the Langevin Eq. (1) was done using a Heun scheme. A resolution of M res = 200 was used in any direction.
To compare the results with other methods, we also perform Brownian dynamics simulation (BDS) using N Brown initially uniformly placed particles in the simulation area. Integration was done using the Heun scheme [7,29] with integration time step h Brown . After a thermalization time T therm,Brown the particles positions were recorded after time intervals ∆t Brown . The BDS was given a running time T run (real CPU time) which usually equals the time our algorithm needs to produce its results. After T run the BDS was stopped and the SPDF was calculated using the recorded particle positions. We set T therm,Brown = T therm , h Brown = h, ∆t Brown = n av h to make results comparable.

A. One dimensional system
In order to demonstrate the implementation of the algorithm, we study overdamped Brownian motion in a bistable potential (1) which results in a bistable system which is often used to study bistable systems or stochastic resonance therein [25,30]. The two stable states come up to the potentials minimums, located at x = −1 and x = 1, respectively. The corresponding SPDF is given by Eq. (2). For low noise strength, the SPDF attains sharp peaks at the potentials minimums and decreases down to low values at the borders and the local maximum, for instance for D = 0.01 p st (0) ≈ 10 −11 .

Equilibration
At first, we study the equilibration process, performing simulations with different numbers of walkers per subregion N . Results are shown in Fig. 2. Analyzing the time dependence of the probability P i (t), we find that longest thermalization time occurs at the local maximum of U (x). Note that runs with larger N thermalize at lower t, but one needs more integration steps.
After thermalization has been achieved, we evaluate the coefficient of variation of the probability, given by where averages · · · are performed for a fixed number N T of sets P and different N and x i . We find it to scale according to 1 √ N (data not shown).

Stationary probability density function
We start to calculate the SPDF p st in the region [L − , L + [ for a small noise strength (D = 0.01). The Using the results shown in Fig. 2, we set the thermalization time T therm = 50 for a run with N = 2. Time averages after thermalization were calculated over an ensemble of N T = 10 4 sets P. Results for the SPDF are shown in Fig. 3. The algorithm calculates the tails of the distribution down to 10 −300 correctly, after a running time T run ≈ 27 hours. We also calculate the SPDF using Brownian dynamics simulation using N Brown = 10 4 , which stops estimating at a level of 10 −6 after the same running time. Further runs were performed (see Tab. I (run 2) and (run 3)), approximating the tails down to 10 −10 (M group = 1136) and 10 −48 (M group = 6789) after a running time of ≈ 30 sec and ≈ 10 min , respectively (data not shown).
Simulations for different values of ∆x L dif indicate that insignificant deviations from the analytic SPDF occur for ∆x > 1 20 .

Probability current
Next, we present that our algorithm can be used to calculate the escape rate to pass the energy barrier at x max = 0. Such problems are typical for chemical reac- tions [31] and in the field of neuroscience [32]. Initially, only N particles are assigned at the subregion including x min = 1, so approximating an initial deltalike probability distribution for t = 0. Furthermore, an absorbing boundary right behind the local maximum (x abs = −0.01) is included. To fulfill normalization of the SPDF, walkers that reach x abs are reinjected immediately at x min . The escape rate to pass the barrier is given by the probability current J(x max ). For small noise intensities, the probability current on top of the potential barrier can be described using Arrhenius law, namely: where ∆U = U (x max ) − U (x min ) = 0.25. Since strong fluctuations are rare, but possible, we will use J(x abs ) to approximate J(x max ). Probability currents J(x abs ) were recorded for each time step and averaged over a sequence of n av = ⌈ 0.1 h ⌉ running steps, resulting in J(x abs ) . Figure 4 shows the time dependence of ln J(x abs ) . After a relaxation regime, where the current decays exponentially, the current reaches its stationary value. The values of ln J(x abs ) , averaged over the stationary regime, are shown in Fig. 5 for different noise intensities. Fulfilling the criteria described above, the algorithm reproduces well Arrhenius law down to ln J(x max ) ≈ −650 corresponding to a current J(x max ) ≈ 10 −286 .

Efficiency compared to weighted-ensemble Brownian dynamics simulation
In order to compare the efficiency of two algorithms, important quantities are the transient time required to first reach the steady state. Once the algorithm reaches the steady state, we quantify the size of the fluctuations by the coefficient of variation Eq. (14) at the potential's    local maximum x = 0. The computation time mainly depends on the number of integrations N int needed, to reach the stationary regime. In order to compare the efficiency, the size of fluctuations (Eq. (14)) in the local minimum relative to p st (0) ≈ 3.88717 × 10 −11 during the stationary regime is plotted over N int in Fig. 6. The most effective algorithm would be located close to the origin. Comparing the efficiency of standard WE simulations and our algorithm, we find that WE simulations with low N equilibrate faster. However, the precision highly depends on the fluctuations during the stationary regime. To produce results of same precision (same c(0)) both algorithms approximately need the same N int . Usually WE simulations were performed using thousands of walkers per subregions, resulting in a low value of c(0).
Here runs of our algorithm producing the same c(0) need much less N and have memory requirements independent of N .
B. Two-dimensional systems

Poincaré Oscillator
As an example of a two-dimensional system with known SPDF, we consider the Poincaré oscillator [27,33,34], represented by the dynamical system: where ξ(t) represent delta correlated white Gaussian noise with zero mean. Using the energy function H(x, y) = 1 2 (x 2 +y 2 ), which only depends on the distance to the origin, one can calculate the associated SPDF: (see appendix Sec. A). Since noise only applies to the ydirection, the lengths of a subregion ∆x and ∆y in xand y-direction are calculated by Eq. (13) and (8), respectively. The minimum of |f x (x, y)| = |y| is equal to 0, therefore, we choose ∆x = ∆y h, which corresponds to a first order approximation of |f x (x, y)| in the next subregion. Results for the SPDF are shown in Fig. 7. Interestingly, the algorithm approximates better the SPDF along the direction where no noise was applied. Here the SPDF is sampled down to 10 −30 . We found that the algorithm slightly oversamples the analytic SPDF in the tails. This is due to the statistical errors, made during the redistribution step. By reducing the size of a subregion, this error can be reduced further. Along the y-direction, noise is applied. Here the behavior is similar as in the one dimensional example (see above). For runs with larger ∆y (results not shown) the algorithm slightly oversamples the SPDF in the origin.

Bistable system with colored noise
As a further example, we calculate the SPDF of the two dimensional system: where τ denotes the time scale separation between x and the colored noise y. The white Gaussian noise ξ(t) has been already described above. This system has been studied previously in [35,36]. Like in the bistable system we have studied above, the SPDF has maximums at (x, y) = (1, 0) and (x, y) = (−1, 0). However, for some combinations of τ and D, the SPDF possesses a local minimum at (x, y) = (0, 0). Plots of the SPDF are depicted in Fig. 8. Our algorithm samples the SPDF down to 10 −12 , whereas BDS breaks down at a level of 10 −6 . The minimum, occurring for τ = 2.5 was clearly found by the algorithm.

FitzHugh-Nagumo-system
As a last example, we consider the widely used FitzHugh-Nagumo-system [37], which is often used in the field of Neuroscience [38] or to study synchronization [39,40] and coherence phenomena [41], represented by:ẋ Here ǫ denotes the timescale separation between the activator variable x and the inhibitor variable y. ξ x (t), ξ y (t), represent independent zero-mean delta-correlated Gaussian white noises. We want to study the stationary probability density in the case of D x = D y = D for a time scale separation ǫ = 0.1. We set the parameters according to Ref. [42] to b = 1.4, ǫ = 0.1 and γ = 2. Thus, the system is in the excitable regime. Since the deterministic part of the equation for the activator variable increases very fast if x is increased, we have to choose a time step h = 0.01, which is small enough, that the walkers' steps are small compared to 1, but allows us to fulfill the criteria for the size of the subregions. Contour plots of the SPDF are shown in Fig. 9. Especially regions of low probability are much better sampled, using the algorithm. In the case of low diffusion (D = 0.01, Fig. 9 (top)) the algorithm runs down to 10 −16 , whereas BDS stops at a level of 10 −6 . Especially the minimum is much better sampled by our algorithm. Note that the local maximum located in the surroundings of (x, y) = (0.7, 0.5) was not found by BDS. For higher diffusion values (D = 0.1, Fig. 9 (top)), significant differences can be found only in the tails.

A. One dimensional system
In the one dimensional system (see section II A) we succeeded to approximate the SPDF down to 10 −300 . If equally sized subregions are used, a huge number of subregions has to be implemented, in order to fulfill the convergence criteria (see section I E). However, by evaluating the SPDF according to Eq. (6) the additional computational costs also reduce the fluctuations of the estimated SPDF. If the subregions are to large, the theoretical SPDF is overestimated by the algorithm in potential minimums and underestimated in the potentials maximums.
We also managed to calculate escape rates of size 10 −286 . Although the finite size of a subregion leads to small errors, in the calculated SPDF, Arrhenius law is well reproduced for such small probability currents.
At the stationary regime the P i fluctuate around their mean value, estimating the SPDF. By either increasing the number of averages N T or the number of walkers per subregion N these fluctuations can be reduced, leading to higher precision. The estimation error scales with √ N T and √ N , respectively. However, increasing the number of walkers per subregion also effects the computational costs during the thermalization. Simulations for different N show that runs with higher N become stationary faster, but this does not compensate for the additional computational costs. We also find, that increasing N slightly improves the sampling of the SPDF's tails. Once the system is in the stationary regime the increase of the computational costs scale linearly with N and N T . An advantage of simulations with small N is, that one does not need to perform running steps for subregions with P i = 0 and the number of such regions naturally increases for small N . We usually use the smallest possible N = 2 and scale the estimation error by increasing N T .

B. Comparison with weighted-ensemble Brownian dynamics
Since the general idea of our method was adapted from prior simulation techniques known as weighted-ensemble (WE) sampling [17,18,20], we want to discuss advantages and disadvantages of our algorithm in comparison with these techniques. The main difference in WE techniques is the redistribution step. Here, using WE techniques, positions and weights of all walkers are stored, and new walkers are introduced on the stored positions considering their particular weights.
In our algorithm, there is no need to store any position or weight, since walkers are randomly placed in each subregion. The resulting statistical errors can be neglected, if the size of the subregions fulfills conditions which, unfortunately, lead to much larger numbers of subregions. By averaging the probabilities, these extra computational costs contribute to an reduction of fluctuations in the stationary regime.
The comparison of the computational costs until equilibration and the achieved precision shows, that both algorithms posses the same efficiency for high precision runs.

C. Two dimensional systems
In the case of two dimensional systems, we found that the algorithm outperforms Brownian dynamics simulations. However, the number of subregions needed according to the criteria can be really high, especially if there are some directions without any noise. Here it is necessary to size the subregions to fulfill the criteria locally. A first step in that direction has already been done by the implementation of the grouping algorithm. We are quite confident that it is possible to reduce the computational costs by optimizing the Mesh. The analyzed examples show, that the running times highly depend on the investigated system. For the bistable system with colored noise and the FitzHugh-Nagumo-system, we obtained good results within ≈ 20 minutes, whereas the algorithm needed about 5 hours for the Van-der-Pol oscillator.

IV. CONCLUSION
We provided and tested an algorithm that allows the calculation of low probabilities and low rates. The algorithm is based on WE Brownian dynamic simulations, but uses a uniform distribution of walkers within each subregion. To our findings, the resulting statistical errors can be neglected if one uses subregions small compared to the diffusion length. In contrast to WE methods, the required memory does not depend on the number of walkers, which leads to less memory requirements for runs with large numbers of walkers.
Special attention was payed to non-equilibrium dynamical systems. Applying the method to one-and twodimensional model systems, we analyze its efficiency compared to standard Brownian dynamics simulation. Our method outperforms Brownian dynamics simulation by several orders of magnitude and its efficiency is comparable to weighted-ensemble Brownian dynamic simulations in all studied systems and lead to impressive results in regions of low probability and small rates.

ACKNOWLEDGMENTS
This paper was developed within the scope of the IRTG 1740 / TRP 2011/50151-0, funded by the DFG / FAPESP. LSG acknowledges IFISC of the University of Balearic Island for cordial hospitality and support. He also acknowledges Dr. Volkhard Bucholtz from Logos-Verlag Berlin for earlier participation and work on the project. RT acknowledges financial support from MINECO (Spain), Comunitat Autònoma de les Illes Balears, FEDER, and the European Commission under project FIS2007-60327, as well as the warm hospitality at Humboldt University.
In this Appendix we show that the presented algorithm is described by a corresponding Master equation for the probability distribution density P i (t) for the case of a Markovian hopping process between boxes. We use a single index i to label the boxes, but the argument applies to any spatial dimension d. We identify the dynamics of the stochastic system which shall be simulated with the discrete stochastic dynamics of the walkers. The latter is defined via the matrices of probabilities per unit time w(i → i ′ ) which describe the hopping in the given discretized space. We assume that it shall converge for sufficiently small time scales and box lengths to the outgoing dynamics.
Let assume that we have the probability distribution density given at time t in every box with index i. It holds i P i (t) = 1. (B1) We redistribute the probability in every box to N walkers. In result any walker k = 1, . . . , N of the same box gets an identical weight until the next new redistribution. At later time t + h the walkers may be still inside the box or may have jumped to other boxes. Let U i denote all possible box-indices which can be reached during a single step from the i-box. Then in accordance with the assumption above, the probability is w(i → i ′ )h that during h a single walker leaves i-box and jumps to the box with index i ′ ∈ U i . By this hopping the walker k carries the weight q k i (t) to the new box, which step is, obviously, connected with a lost in the outgoing box. Therefore, the lost per particle for this specific hopping from i → i ′ can be expressed as The whole lost of weight will be realized on all possible hopping channels. It is identical for all N particles being located in the present box. Hence, the full lost becomes Alternatively, one can also introduce U ′ i as the boxes where from walkers can reach the i-box. Then, gain of weight transferred by every walker arriving at the i-box is expressed, if i ′ ∈ U ′ i , by Again, summing over the different hopping steps and considering that all the N walkers inside the box with i ′ ∈ U ′ i reach the i box, yields the gain Therefore, the balance of transported weight results in the following shift of the full weight in the i-box at time t + h compared to the former one plus corrections of order O(h 2 ) corresponding to cases in which two or more walkers jump from one box to another during the time interval h. After dividing by the time step h and taking the limit h → 0 we obtain the wanted Master equation for the discretized dynamics in the boxed phase space