Kinetic Monte Carlo modelling of dipole blockade in Rydberg excitation experiment

We present a method to model the interaction and the dynamics of atoms excited to Rydberg states. We show a way to solve the optical Bloch equations for laser excitation of the frozen gas in good agreement with the experiment. A second method, the Kinetic Monte Carlo method gives an exact solution of rate equations. Using a simple N-body integrator (Verlet), we are able to describe dynamical processes in space and time. Unlike more sophisticated methods, the Kinetic Monte Carlo simulation offers the possibility of numerically following the evolution of tens of thousands of atoms within a reasonable computation time. The Kinetic Monte Carlo simulation gives good agreement with dipole-blockade type of experiment. The role of ions and the individual particle effects are investigated.


Introduction
From the seventies, physics of the Rydberg atoms has been an object of great interest. Most of the properties of Rydberg atoms are due to the dimension of the Rydberg orbit, typically in atomic units of the order of the square of the principal quantum number n. Possessing huge electric dipole moments, large lifetimes..., Rydberg atoms have offered the opportunity of studying atoms in extreme experimental conditions, for instance in presence of high electric, magnetic or electromagnetic fields or approximating the conditions which corresponds to (low n) atoms in the neighborhood of a star [1]. More recently the physics of Rydberg states of atoms in cold gases has stimulated interest since they are at the frontier of atomic, molecular, solid-state and plasma physics. In an ensemble of cold Rydberg atoms many-body phenomena have been observed in a Förster configuration where Rydberg atoms can exchange internal energy through long-range dipole-dipole interactions [2,3]. The possibility of controlling those strong interactions between atoms has been demonstrated by using an external controllable electric field [4].
A basic difference between experiments with cold Rydberg atoms and those with Rydberg atoms at room temperature is that cold Rydberg atoms on the (∼1 µs) timescale of the experiments can be approximately considered motionless. For example, in the case of cesium atoms they move 100 nm which is roughly the size of the atom, for n ∼ 30. Such a cold gas is expected not to exhibit any collisions and presents totally different characteristics than does a thermal gas at room temperature. No collisions does not mean no interactions, and a frozen Rydberg gas can present novel properties, close to those of an amorphous solid. The frozen Rydberg gas approximation leads to considering the ensemble of Rydberg atoms as interacting at large distances by the van der Waals interaction or the dipole-dipole one. The Förster configuration leads to a situation very similar to the migration of excitons [2,3]. Thus, fascinating perspectives are expected with cold Rydberg atoms. Controllable long-range interactions are particularly exciting for quantum information applications especially the so called dipole blockade mechanism of the excitation due to the strong interactions between Rydberg atoms [5,6]. The energy of a pair of interacting Rydberg atoms is shifted by dipole-dipole interactions and is not twice the energy of one Rydberg atom. A limitation of the excitation is expected when the dipole-dipole energy shift exceeds the resolution of the laser excitation. The use of a dipole blockade of the excitation constitutes an efficient way for the realization of a CNOT quantum gate [5,6,7,8,9,10]. The possibility of observing the dipole blockade of laser excitation has been demonstrated for the first time with the van der Waals case [7,8] and the Förster case [9]. Modelling the complex behavior of the dipole-dipole interaction in a frozen gas opens interesting ways of understanding the role of each particle by switching on and off the different interactions or effects. The advantage offered by our simulation is the possibility of selectively adding effects/interactions depending on their rates with up to thousands of particles under reproducible conditions within a computational time of a few minutes. After a review of our experimental conditions, we describe the different methods we have been using to model the dipole blockade effect observed in [9,10]. We briefly explain a first method based on the solution of the optical Bloch equations, then we discuss the use of a Kinetic Monte Carlo (KMC) simulation. We then present the results we obtained with the KMC model for different experimental situations. Due to the wide utility of algorithms used, we have presented in appendix a review of KMC method itself and the algorithm for the motion of the particles.

Experimental setup
In many experiments with hot or cold Rydberg atoms, the experimental procedure is the following. The atoms are excited by a short laser pulse (∼ 10 ns) to a Rydberg state, nl (l = s, p, d). Then after a duration of a few microseconds the Rydberg gas sample is selectively state analyzed by using a high voltage pulsed electric field with a risetime of the order of a microsecond. An important difference is observed between experiments realized at room temperature, using for instance a thermal atomic beam, and those realized with a cold atomic sample provided by a magneto-optical trap (MOT). In the case of cold atoms large fluctuations of the number of Rydberg atoms are generally observed between laser-shots. The reason is the very narrow linewidth (∼ 1 kHz if excited from the ground state, ∼ 5 MHz if excited from the 6p in cesium) of the Rydberg resonance, compared to the broad bandwidth multimode laser which cavity modes oscillates randomly (multiple cavity modes spread over a few GHz). It leads to uncontrollable frequency shifts of ∼ 500 MHz. In an atomic beam, the Doppler effect can be up to 1 GHz which limits the fluctuations of the Rydberg population. In the case of the atoms in a MOT, there is no Doppler effect, which explains the strong fluctuations. In broadband experiments the excitation of an ensemble of Rydberg atoms interacting altogether corresponds to the excitation of a band of energy levels, which can be excited by a short, thus broadband, laser pulse. The width of the band versus the Rydberg atomic density has been investigated by microwave spectroscopy [8] and laser spectroscopy [11]. Using monomode lasers for the excitation of cold atoms is a way to avoid the fluctuations of the Rydberg population from shot to shot.
Another important difference is expected between a broadband excitation and a high-resolution one. With narrow band, low power excitation, only a small part of the band of levels can be excited, leading to the limitation of the excitation and corresponding to a van der Waals [7,12] or dipole [10] blockade. The first excited Rydberg atoms shifts the resonance of the non-excited neighbors and prevent their excitation in a narrow-bandwidth laser excitation.
The details of our experimental setup have been described in several papers [2,9,4,11]. The Rydberg atoms are excited in a cloud of up to 5 × 10 7 Cs atoms (temperature 200 µK, characteristic radius ∼ 300 µm, peak density 1.2 × 10 11 cm −3 ) produced in a standard vapor-loaded MOT at residual gas pressure of 3 × 10 −10 mbar [2,4]. At the trap position, a static electric field and a pulsed high voltage field can be applied by means of a pair of electric field grids spaced by 15.7 mm. We consider an ensemble of cold cesium atoms excited in the Rydberg state, np 1/2 or np 3/2 . Three cw lasers provide a high resolution multistep scheme of excitation, as depicted in figure 1 A).   The first step of the excitation, 6s, F = 4 → 6p 3/2 , F = 5, is provided by the trapping lasers (wavelength: λ 1 = 852 nm) or a diode laser to avoid the excitation of hot atoms. The density of excited, 6p 3/2 , atoms can be modified for instance by switching off the repumping lasers before the excitation sequence. The second step, 6p 3/2 , F = 5 → 7s, F = 4, is provided by an infrared diode laser in an extended cavity device (wavelength: λ 2 = 1.47 µm, bandwidth: 100 kHz and available power: 20 mW). The average experimental intensity is ∼ 3 mW/cm 2 , twice the saturation one. The last step of the excitation, 7s, F = 4 → np 1/2,3/2 (with n = 25 − 300), is provided by a Titanium:Sapphire (Ti:Sa) laser. The wavelength λ 3 ranges from 770 to 800 nm, the bandwidth is 1 MHz, and the available power is 400 mW. The Ti:Sa laser is switched on for a time, τ = 0.3 µs, by means of an acousto-optic modulator at an 80Hz repetition rate. The beams of the infrared diode laser and of the Ti:Sa laser cross with an angle of 67.5 degrees and are focused into the atomic cloud with waists of 105 and 75 µm, respectively. Their polarizations are both linear and parallel to the direction of the applied electric field, leading to the excitation of the magnetic sublevel, np 1/2 or np 3/2 |m| = 1/2. The spectral resolution, ∆ν L , of the excitation is 5 − 6 MHz, limited by the lifetime, 56.5 ns, of the 7s state and by the duration, i.e by the spectral width of the Ti:Sa laser pulse. The magnetic quadrupole field of the MOT is not switched off during the Rydberg excitation phase, but it contributes less than 1 MHz to the observed linewidths. Just after the Ti:Sa laser pulse (between 0 and 1µs) the Rydberg atoms are selectively ionized by applying a pulsed high-voltage field with a rise time of 700 ns and detected on a Micro Channel Plate (MCP) detector. The experimental procedure is based on spectroscopy of Stark np states for different atomic densities, and for different Ti:Sa laser intensities.

Dipole blockade model
Different approaches have been followed to study the problem of the excitation to a Rydberg state in presence of already excited atoms [13,14,15]. We discuss hereafter some hypotheses and simplifications we made to model the blockade effect. A first simplification is made by considering the excitation to np states only. The dipoledipole interactions are calculated for first and second orders between the np and all the neighboring states in the Stark diagram. When looking at a specific atom to be excited, the shift in energy is due to already excited neighboring atoms but not ground state atoms. As shown in the next section, the sum of each individual atom's contribution can be studied and the main effect is due to the nearest neighbor.
If a static electric field is present, either external or due to ions, Rydberg states are mixed, creating a permanent dipole moment for the Rydberg atoms. For instance in Cesium in the presence of an electric field F , the np state is mainly mixed with the (n − 1)d state. We denote by µ pd = np j , m j = 1/2 |q e z| (n − 1)d j+1 , m j+1 = 1/2 the transition dipole moment of an atom in state np (j, m j = 1/2) toward (n − 1)d (j + 1, m j+1 = 1/2). We introduce the scale parameter θ characterizing the dipole coupling for each level np defined by tan θ = µ pd .F h∆ pd /2 , where h∆ pd is the zero field energy difference between the (n−1)d and np levels. Energies and dipoles are obtained following [16], and are calculated only for |m j | = 1/2 states. The dipole of an atom in a np state aligned along the local electric field ( F ) is given by (here z is the coordinate along the vector defined by F ): Where the basis (|np( F ) >,|(n − 1)d( F ) >) are the eigenstates given by the diagonalization of the Hamiltonian matrix where E p and E d are the energies of states np and nd in absence of an electric field. The resulting shift in energy for np is h∆ p ( F ) = h∆ pd One can calculate the dipole-dipole interaction term V ij between two atoms labelled i and j separated by R ij = R ij n ij . The first order dipole-dipole interaction is : where h∆ pk is the energy difference between states np and n ′ k, and µ( is the classical permanent electric dipole of atom i which is aligned along the local electric field F i . In the absence of an electric field, θ = 0 and there is no permanent dipole moment and only the second order, so called van der Waals interaction is non zero. The potential energy of a np Rydberg atom is the sum of the energy of the state without electric field E p plus the shift of the state in the local electric field ∆ p ( F i ) plus the sum of the dipole-dipole interactions with all the atoms. For the atom i we then note An important part of the computation time is used to calculate the local electric fields and the potentials.

Reduced density matrix, mean field simulation
The details of this work can be found in [17]. We just briefly review the main results. The process describing the three excitation steps (see figure 1 A)) for one atom can be described using the optical Bloch equations.
Where the time evolution of the density matrix ρ is decomposed into three terms. The first term contains an Hamiltonian H being the sum of the individual potential energies, and the interaction of an atom with all the others V ij in presence of electric fields The second term accounts for the relaxation of the populations and coherences where Γ gives the lifetime of the considered state. The third term γ takes into account the radiative relaxation to state l with energy E l from states k with energy E k > E l due to spontaneous emission. Taking the trace over all the atoms except the one labeled i in the optical Bloch equations, gives the evolution of the density matrix for the particle i. The interaction term or shift in energy for the atom i due to the interaction with its neighboring atoms is j =i T r j [V ij , ρ i,j ], with V ij the dipole-dipole interaction and ρ i,j the two-body density matrix for atoms i and j. The coupling with all the other atoms becomes a mean field term proportional to the atom density. A similar treatment has been performed by [18].
As correlations appear during the excitation, the state of the system does not remain a product state. However the probability of excitation of a ground state atom into a Rydberg state being on the order of a few percent, and as long as the product of the individual density matrices is small, we can use the Hartree-Fock approximation. In this approximation the two-body density matrix ρ i,j can be developed as the product of single atom density matrices. We start with one atom i in ground state and atom j in the excited state, denoted (ge). After excitation, the system ends in a double excitation noted (ee). The Hartree-Fock approximation allows us to write, ρ (i,j)ge,ee = ρ ige × ρ jee .
Thus the interaction term can be written as which is simply a shift of the Rydberg level for the atom i. As an illustrative example, we look at a weak interaction with tan 2 (θ)=0.05, and with a 70p 3/2 state. The dipole blockade effect induces a shift of 6MHz (exactly our excitation linewidth) which prevents the excitation of two atoms at a distance of 5µm. At a density of 10 11 cm −3 a sphere of radius 5µm contains 50 atoms. This means that only one excitation could be present for the 50 atoms, and the probability to excite the atom j would be uniform within this sphere. The population in the excited state ρ jee is then replaced by a mean value ρ ee . Due to the inhomogeneity of the atomic density and laser intensity a local ρ ee ( r) is considered at different positions r over the whole atomic cloud. A naive (mean field) estimation for ρ jee could lead to wrong estimations. Indeed, a mean field interaction for an atom at the center of the cold atomic could naively be written as the integral of the interaction term V ij over all the possible directions Θ (the angle between the internuclear axis and the direction of the dipole i): π 0 V ij sinΘdΘ which is equal to zero, but is not the real value. In order to overcome this problem, a better way to evaluate the local interaction potential is to consider separately the nearest neighbor Rydberg atom from the other atoms. The shift in energy δ dd (i) is decomposed into a sum over all the atoms j treated as a continuous distribution out of a sphere containing only one excited atom (the nearest neighbor of i) in its center, plus the contribution of the nearest neighbor Rydberg atom. The result from this calculation is that the local field contribution associated to the nearest neighbor is dominant over the mean field contribution if an electric field is present. The nearest neighbor Rydberg atom contribution is considered at the most probable distance given by the mean of the Erlang distribution ‡ [19] from the atom i. The shift in energy relies on the local density (gaussian distributed) ρ 0 ( r) of the atoms in ground state. We finally find that the shift for the atom i is given by We then solve equation (5) for an atom i using the result from equation (7). The result given in figure 2 reproduces well the experiment. We take into account multiphotonic excitations as well as the finite coherence time of the lasers in the model through a temporal phase variation of the electric field of the lasers in the three step excitation represented in figure 1 A). Two results are given in figure 2, where the reduced density matrix approach is plotted versus the experimental data for different electric fields A) and different intensities B).

Kinetic Monte-Carlo (KMC) simulations
The previous approach based on the reduced density matrix has some limitations. Despite correctly describing the excitation it was not possible to look for the dynamics ‡ The probability to find a k th nearest neighbor at a distance r of a an atom is given by of the system, the orientation of the dipoles in a local electric field or the individual interactions between atoms instead of a mean field term or ionization. For these reasons we developed a KMC simulation. All the above mentioned limitations can then be overcome. However in KMC simulations the excitation has to be based on the solution of rate equations. A more detailed description of the KMC algorithm, and more generally of possible numerical solution of any kind of master or rate equations, is given in Appendix A. Briefly if a system is driven by a master equation describing the time evolution of the probability P k of a system to occupy each one of a discrete set of states numbered by k. Each process occurs at a certain average rate Γ lk , which may either be constant in time, or dependent on how the system has evolved up to that time.
The KMC algorithm is then the iteration of the following steps.
• Initializing the system to its given state called k at the actual time t.
• Creating the new rate list Γ lk for the system, l = 1, . . . , N.
• Choosing a unit-interval uniform random number generator § r: 0 < r ≤ 1 and calculating the first reaction rate time t ′ by solving • Choosing a unit-interval uniform random number generator r ′ : 0 < r ′ ≤ 1 and searching for the integer l for which and R 0 = 0. This can be done efficiently using a binary search algorithm.
• Setting the system to state l and modifying the time to t ′ . Then go back to the fist step.
One fundamental result is that the KMC method makes exact numerical calculations and cannot be distinguished from an exact molecular dynamics simulation, but is orders of magnitude faster. It is therefore indistinguishable from the behavior of the real system (if evolving through a master equation), reproducing for instance all possible data in an experiment including its statistical noise.

KMC simulation of dipole-dipole interaction
We consider a cloud with a gaussian spatial density. Initially the thousands of atoms are in their ground state with a maxwellian distribution for the velocity σ v = k b T /m at , where T, k b , m at are respectively the temperature of the frozen gas in the MOT, the Boltzmann constant and the mass of the atomic species under consideration. In such a system, considering coherent excitations would lead to solving the Schrödinger equation with ≈ 2 100 states, which is obviously beyond our capability. Nevertheless, it is possible to obtain good agreement for the shift in energy and the dynamics of the system with a § In our case we use the free implementations by GSL (GNU Scientific Library) of the Mersenne twister unit-interval uniform random number generator of Matsumoto and Nishimura.
reasonable computation time using a simplified numerical treatment. We consider a two level system for each atom so they can be either in their ground state (6s 1/2 ) or laser excited to a given (np 3/2 ) Rydberg state. Initially, the electric dipoles µ of the atoms are aligned along the direction of polarization z of the exciting laser and the applied DC electric field, which is the quantization axis. During the evolution (ionization especially) the electric dipoles will be aligned along the local electric field. The shift δ dd (i) of the Rydberg state of atom i is then j =i V ij where V ij is given by equation (3). We start from the sets of Bloch equations for our two level system [20]: where ρ gg and ρ ee are the populations of ground state and excited state atoms, Γ spont is the spontaneous decay rate, Γ laser is the FWMH of the exciting laser(Γ spont << Γ laser ), Ω is the local Rabi frequency of the transition and δ laser is the detuning from the isolated and field free atom resonance. In order to end up with rate equations for each atom, we neglect the coherencesρ ge =ρ eg = 0. This is obviously the biggest assumption made here. Coherences are in fact tractable with a Monte Carlo method as described in [21]. We finally get pure rate equations for ρ ee and ρ gg : where the rate Γ exc of excitation of atom i depends on the projection of |np( F i ) > onto |np >, which is cos 2 (θ i /2), times the stimulated laser rate. The deexcitation rate is similar to the excitation rate plus the spontaneous decay of the np state. The great advantage offered by this 'kinetic' method, is the simplicity of adding phenomena or particles with evolution based on rates. Contrary to the previous model the potential energy of an atom is the sum of the energies due to the Stark shift and the interactions with all the other atoms. Numerical methods such as the ordinary Runge-Kutta methods are not ideal for integrating Hamiltonian systems because they do not conserve energy. On the contrary symplectic integrators such as the Verlet integrator does conserve energy. Mechanical effects are then treated via classical movements of the atoms, and dipolar forces between atoms are taken into account using the (leapfrog) Verlet-Störmer-Delambre algorithm (see Appendix B). The computation is realized as follows. Ion formation is possible and creates an electric field felt by other atoms, thus the direction and strength of the dipoles can vary depending on local electric fields. Two main mechanisms for ionization exist. The first one is due to laser ionization or blackbody ionization and has a rate of ionization proportional to the atomic density. The second one happens if two Rydberg atoms move toward each other and reach a smaller internuclear distance than 4n 2 a 0 [22]. In the latter case one Rydberg atom is ionized, the second atom falls to a lower state. Due to energy conservation, its binding energy is at least twice as large after the ionization, and as the final atomic states are often different from s or d, it does not interact with np atoms. Consequently we assume for simplicity in the simulation that the atom state is changed to a non interacting state. Due to its time and spatial resolution the KMC simulation takes into account all the dipolar interactions developed during the excitation and gives access to individual atoms. A dynamical evolution of the system is made except if a collision between two atoms in np state is detected or if a reaction is detected. The time evolution of the simulation is incremented either with the KMC timestep or with a small fraction of the collisional timestep. After each change in the position or change of any particle state, the fields and potentials are recalculated over all the atoms. Then operations are repeated until the end of the excitation.

Field induced dipole blockade
As described in the experiment reported in [10] the dipole moment of np states increases with the strength of the coupling field and reaches a maximum when tan θ is equal to 1. Indeed, as the strength of the field increases the blockade radius increases and the number of excited atoms in a given volume gets smaller. It is worth noting that the dipole blockade condition in a 2-level approach is j =i V ij > hδ laser . Figure 3 A) represents the result of a Monte Carlo simulation for different np states, where the number of Rydberg atoms present at the end of the excitation is given as a function of the detuning. The parameters are close to those described in [10].

Effect of the ions
At a distance of 10µm the electric field due to an ion is 150mV/cm which represents a shift of 150MHz for an atom in state 70p. The contribution to the blockade of the excitation of such an ion in the sample during the excitation is so important that it completely hides the observation of the dipole blockade effect. In the experiment the ions present before and during the excitation can be discriminated from the ionized Rydberg atom in the time of flight signal. In the simulation we simply monitor the number of ions at the end of the the excitation time. This ionic effect is shown in figure  4. In figure 4 A) we increase the laser intensity to ten times the saturation intensity and we apply the excitation laser for 20µs in absence of external electric field. One can see that the number of excited atoms is important but also that the resonance is broadened to the low frequency side of the atomic resonance. This result is similar to the one obtained in the experiment [12]. Ions are formed due to the high laser power but also by collisions during the interaction time due to the long range attractive dipole force between atoms. The number of ions produced during the excitation is so important that it allows for the excitation of Rydberg atoms on a 100MHz range, red detuned from the center of the line. This is due to the fact that np states can only be shifted to the red of the resonance when no external electric field is present. In figure 4 B) is shown an experimental curve obtained in a combined pulsed and cw-excitation (see figure 1 B)). In this case the main source of broadening is the inhomogeneous electric fields due to

Role of the nearest neighbor
As the dipole-dipole interaction term V ij strongly depends on the distance between two atoms, the energy shift due to the nearest neighbor Rydberg atom has to be distinguished from the mean field shift due to all the other atoms in a Rydberg state. In some cases the contribution of the nearest neighbor can be dominant.
In figure 5 is represented for two different densities, but with no ions to avoid extra effects, the nearest neighbor Rydberg atom shift versus the sum of the shifts of all the Rydberg atoms in the sample. Above the solid line the contribution of the nearest neighbor Rydberg atom to the blockade effect is dominant. The nearest neighbor Rydberg atom shift is twice more important than the shift from all the others Rydberg atoms for 66% of the ground state atoms at a density of D=5 × 10 9 cm −3 ( figure 5 A)), whereas at D=5 × 10 10 cm −3 (figure 5 B)) it is 73%. This confirm the fact that the energy shift is dominated by the effect of the nearest neighbor. This result has been used to derive the local mean field approximation in the single atom density matrix model previously described.

Conclusion
In this paper we have presented different methods used to model the dipole blockade effect in a manner as close as possible to an experimental situation. We have first shown that a nearest neighbor mean field analytical approach, based on the solution of the Bloch equations for the partial density matrix equations of one atom of interest gives results in good agreement with the experiment if no ions are present. Second the Kinetic Monte Carlo simulation based on rate equations is able to introduce all the electric dipole interactions. Furthermore, we have included N-body spatial dynamics using the Verlet integrator, but the overall code remains very simple. This allows us to look at the important role of the ions formed during the excitation. If ions are present during the excitation they may mimic the dipole blockade effect and lead in extreme cases to a broadening of the lines. We have shown that the two models reproduce well our experiments. It is observed that the number of excited atoms for a given ground state density decreases with the principal quantum number and with the intensity of the electric field which creates the permanent dipole. The effects associated with the experiment, such as the shift of the resonance, the density effect, the broadening of the line, the artifact due to the ion blockade are reproduced with the simulations we describe. For varying initial atomic densities the amplitude of the energy shift due to the nearest neighbor is monitored which confirms its dominant role in the dipole shift. We could also use the Kinetic Monte Carlo model to analyse other experiments than ours as in figure 4 A) [12]. As another simple example, the number of atoms excited to a Rydberg state in our simulation is analyzed using the the Mandel Q parameter and the statistics follow a sub-Poissonian distribution as described in [23,24,13,18]. The density matrix model can be a good tool to model recent experiments of coherent excitation of Rydberg atoms [25,26,27,28]. The next step will be the modelling of the Förster case [9] where the coherences between pairs of atoms are expected to play a dominant role.

Acknowledgments
This work is in the frame of "Institut francilien de recherche sur les atomes froids" (IFRAF). This theoretical development could aslo be applied to the CORYMOL experiment supported by an ANR grant (No. NT05-2 41884). The authors thank Thomas F. Gallagher and Jianming Zhao.

Appendix A. KMC method for solving rate equations
Appendix A.1. The master, kinetic or (reaction-)rate equation There are innumerable instances, in physics and in other sciences, where a system evolves in time through many competing internal stochastic processes. For instance, many classical and quantum physical problems can be reduced to the form of a master equation: This master equation, sometimes called kinetic or (reaction-)rate equation is a phenomenological set of coupled first-order differential equations describing the time evolution of the probability P k of a system to occupy each one of a discrete set of states numbered by k. In probability theory, this identifies the evolution as a continuous-time Markov process. Each process occurs at a certain average rate Γ lk , which may either be constant in time, or dependent on how the system has evolved up to that time. The goal of this appendix is to describe why Kinetic Monte Carlo methods are a standard means of modelling such problems, especially when one wishes to model the evolution of the system over periods of time much longer than those accessible by direct simulation.

Appendix A.2. Solving the master equation
Following Gillespie [29] we can distinguish, among several competing methods commonly used to solve the master equation, two major approaches: The deterministic approach which regards the time evolution as a continuous, wholly predictable process governed by a set of coupled, ordinary differential equations (the "reaction-rate equations") and the stochastic approach which regards the time evolution as a kind of random-walk process which is governed by a single differential equation (the "master equation") governing the time-dependent behavior rather than a fixed probability distribution.
Appendix A.2.1. Deterministic approach The deterministic approach is based on the fact that equation (A.1) can be written as a matrix ordinary differential equations dP dt = ΓP The formal solution can be obtained for instance by using Direct diagonalization algorithms [30]. A second deterministic approach to the time-dependent population distribution comes from the integrand form of the master equation: Explicit numerical integration can in principle be achieved by different numerical integration schemes. All these methods require the explicit formation of the matrix and the computational effort is dominated by N 3 terms [30]. The deterministic approach is simply the exact time-evolution for the function P. However, the stochastic probabilistic formulation has often a stronger physical basis, especially in the quantum world or for non equilibrium systems, than the deterministic formulation. Instead of the deterministic approach which deals only with one possible "reality" of how the process might evolve under time, in a stochastic or random process there is some indeterminacy in its future evolution described by probability distributions. This means that even if the initial condition (or starting point) is known, there are many possibilities where the process might go to, but some paths are more probable and others less. This approach correctly accounts for the inherent fluctuations and correlations that are necessarily ignored in the deterministic formulation. In addition, as we shall see, this point of view, opens the way to stochastic simulation algorithms, such as the Kinetic Monte Carlo one, making exact numerical calculations which are much faster O(N) than the deterministic reaction-rate algorithms.

Appendix A.2.2. Monte Carlo algorithms
Monte Carlo refers to a broad class of algorithms that solve problems through the use of random numbers [31,32]. They first emerged in the late 1940's and 1950's as electronic computers came into use. The most famous of the Monte Carlo methods is the Metropolis algorithm (sometimes called Monte Carlo Markov chain methods), offering an elegant and powerful way to generate equilibrium properties of physical systems [33]. For many years, researchers thought Monte Carlo methods could not be applied to molecular dynamics simulations because it seems necessary to follow individual motion and/or interactions [34,35,36]. However some Monte Carlo algorithms of this type exist; they are sometimes called Coarse-Grained methods [37]. The simplest algorithm of this kind, sometimes called fixed time step algorithm [38], is based on the first-order formula P(t + dt) = P(t) + Γ(t)P(t)dt i.e.
A similar scheme is, for instance, used in stochastic quantum simulation in the Quantum Monte Carlo wave-function approach [39] because the Lindblad equation in quantum mechanics is a generalization of the master equation describing the time evolution of a density matrix.
To illustrate the fixed time step method [38], let's assume that at time t the system is in state k: P k = 1 and P l = 0 for l = k. The algorithm consists of choosing a small dt, and for each possible reaction k → l generating a random number r between 0 and 1. If r < Γ lk (t)dt the system changes configuration and evolves to state l at time t + dt (a quantum jump occurs in the stochastic quantum simulation terminology). The main disadvantage is that dt has to be small enough to maintain accuracy and such that at most one reaction occurs during each time step: meaning Γ lk (t)dt ≪ 1. Several steps are then needed before effectively doing the evolution. As we shall see, especially for time independent Γ lk rates, this Monte Carlo algorithm is very inefficient compared to the kinetic Monte Carlo algorithm which ensures at each time step an evolution of the system.

. Derivation
The Kinetic Monte Carlo algorithm is also known as the residence-time algorithm, the n-fold way, the Bortz-Kalos-Liebowitz (or BKL) algorithm [40], the dynamical Monte Carlo method [41], the Gillespie algorithm [42,29], the Variable Step Size Methods (VSSM) [38] (by comparison with the fixed time step method) ..., depending on the physical or chemical context. For other reviews of the KMC algorithm see [38,32,33,43]. Some minor subtle changes between these algorithms exist, but we will describe here only 3 types of KMC algorithms, namely the Kinetic Monte Carlo (KMC) algorithm, the First Reaction Method (FRM) and the Random Selection Method (RSM).
The Kinetic Monte Carlo method is a Monte Carlo method intended to simulate the time evolution of independent (non correlated) Poisson processes [38,44]. This means that the KMC method solves the master equation and is therefore of great interest, for instance for relaxational processes and transport processes on mesoscopic to macroscopic time scales. Indeed, the KMC algorithms are able to model the evolution of the system over periods of time very much longer than those accessible by direct simulation such as molecular dynamics. Surprisingly enough, up to now it has been more or less limited to the study of chemical reactions, surface or cluster physics (diffusion, mobility, vacancy motion, transport process, epitaxial growth, dislocation, coarsening, ...).
Due to the Markovian behavior the system loses its memory of how it entered state k at time t. Therefore, in order to simulate the stochastic time evolution of such a reacting system, i.e. in order to move the system initially at time t in state k forward in time, we just need to know when the next reaction will occur, and what kind of reaction it will be [29]. We then need to determine the probability distribution function p(t) for the time of the first system change. From equation (A.2) the probability that the system has not yet escaped from state k at time t ′ is given by p survival (t ′ ) = exp t ′ t Γ k (τ )dτ , where Γ k = l Γ lk is the total rate from state k. Thus 1 − p survival (t ′ ) gives the probability that the system has been modify at time t ′ , which is exactly the integral t ′ t p(τ )dτ . The first-passage-time distribution in then found by differentiation: which is characteristic of the Poissonian nature of the process and is the starting point of the KMC algorithm. We then know when the next reaction will occur. We just have to find what kind of reaction it will be. At time t ′ a reaction takes place, so just before t ′ the system is still in state k and according to the master equation (A.1) the probability that the system will be in configuration l at time t ′ + dt ′ is Γ lk (t ′ )dt ′ , where dt ′ is a small time interval. We therefore have to generate a new configuration l by picking it out of all possible new configurations with a probability proportional to Γ lk (t ′ ). The algorithm is based on our ability to generate a time t ′ from equation (A.4) when the first reaction actually occurs and on our ability to choose the correct reaction l.
The t ′ choice can be done by the inverse transform method [45] which is based on the fact that if t ′ is chosen with probability density p(t ′ ), then the integrated probability up to point t ′ , is sometimes possible to improve the speed of the algorithm, for instance if several rates are similar, using a binning or weighted methods to reduce for instance to O(log N) the complexity [46,47,48,38].
It is of interest to mention at least the First Reaction Method (FRM) and the Random Selection Method (RSM) [38]. Depending of the type of reaction in the system these algorithms can be useful. The FRM, also called Discrete Event Simulation in computer science, consists of choosing the first occurring reaction, meaning choosing the smallest time t ′ l , and the corresponding reaction number l, from the formula t ′ l t Γ lk (τ )dτ = − ln r l where the r l ∈]0, 1] are N independent random numbers. As the KMC one this algorithm generates an exact evolution of the system but is usually less efficient because it necessitates a random number for each possible reaction, whereas the KMC advances the system to the next state with just two random numbers.
The RSM can be used only when the rates Γ lk (t) = Γ lk are time independent, as for Poisson processes (such as radiative lifetime, radiative decay rate, ...). The RSM consists of evolving the system up to the time t ′ = t − (ln r)/Γ max where r ∈]0, 1] is a random number and Γ max = max l Γ lk is the maximal possible rate. Then choosing randomly a possible reaction l ∈ [1, N] and accepting the reaction with probability Γ lk /Γ max . If the reaction is accepted, the configuration is changed. Contrary to the FRM or the KMC algorithms it does not necessarily imply a system evolution at each time step. But, here again this algorithm generates an exact evolution of the system. The RSM is optimized for system having just one (or a small number of) type of reaction because it is then of O(1) complexity ! Following [38], KMC is generally the best method to use unless the number of reaction types is very large. In that case use FRM. If you have a type of reaction that occurs almost everywhere, RSM should be considered. Simply doing the simulation with different methods and comparing is of course the best.

Appendix A.4. Conclusion
The KMC algorithm is a stochastic algorithm generating quasi classical trajectories, "i.e." creating a Markov chain representing the exact evolution of the system in the sense that it will be statistically indistinguishable from an exact dynamics simulation. Indeed, each system configuration l is reached with its real physical probability. Unlike most procedures such as the often used fixed time step method for numerically solving the deterministic reaction-rate equations, this algorithm never approximates infinitesimal time increments dt by finite time steps ∆t = t ′ − t.
The fact that the mechanisms and so the rates have to be known in advance is the main limitation of the use of the KMC method. If the rate have to be modified in an unpredictable fashion during the free time evolution of the system, i.e. between two reactions, the fixed time step method has to be preferred. However for several physical system this is not the case and the KMC or FRM algorithms can be used. As the interaction between particles of a system depends often of the distance between particles, the KMC methods can be advantageously used when the motion of the particles is slow. Ultimately when the rates are time independent (which does not mean that they are constant, because they often have to be recalculated after each system evolution) KMC and RSM approaches are very powerful.
Appendix B. Simple method to solve the N-body problem Appendix B.1. Introduction A large number of physical systems can be studied by simulating the interactions between the particles constituting the system. In a typical system each particle influences every other particle. The interaction is often based on an inverse square law such as Newton's law of gravitation or Coulomb's law of electrostatic interaction but, as in our case, a more complex anisotropic interaction with an inverse higher power law dependence might exist. Examples of such physical systems can be found in astrophysics, plasma physics, molecular dynamics and fluid dynamics. Since the simulation involves following the trajectories of motion of a collection of N particles, the problem is termed the N-body problem.
Since it is not possible to solve the equations of motion for a collection of many particles in closed form, iterative methods are used to solve the N-body problem. At each discrete time interval, the force on each particle is computed and this information is used to update the position and velocity of each particle. A straightforward computation of the forces requires O(N 2 ) work per iteration. The rapid growth with N limits the number of particles that can be simulated by this method. Several approaches, especially that by Aarseth [49], have been used to reduce the complexity per iteration and to speed up the calculation, for instance each particle is followed with its own integration step. Non full N-body codes also exist transforming the problem imposing for instance a grid on the system of particles and computing cell-cell interactions. These are known as hierarchical methods, or tree methods such as the Barnes-Hut one where a Ahmad-Cohen neighbor type of scheme is used which updates less frequently the non neighbor force than the neighboring one. Finally, multipole expansion methods have also been developed as well as Monte Carlo algorithms for very large number of particles using a set of representative "macro" particles (not point-masses) like in Fokker-Planck or gaseous methods with smooth potentials after the pioneering work of Hénon [50]. The required CPU time scales with the number N of particles as N ln(N) (for a given number of relaxation times), while the scaling is N k with k of order 2−3 for direct N-body codes. It is beyond the scope of our article to discuss all the methods .. Very good references, discussing also their stability or their complexity can be found on the web site http://www.manybody.org/ or in the of O((∆t) 4 ) for local truncation error for position and O((∆t) 2 ) for velocity. We are interested in having accuracy in position and velocity so we use the so called Velocity Verlet method which is of O((∆t) 3 ) accuracy for both position and velocity for a ∆t timestep. This leapfrog integrator often turns out to be more accurate than expected from a simple second-order integrator [57]. The 'unreasonable' accuracy stems from its symmetry properties under time invariance due to its simplectic structure. The scheme of the algorithm is the following.
This has also the big advantage that accuracy can be improved by using higher order symplectic integrators such as the one by Yoshida [58].