Diffusion of skyrmions: the role of topology and anisotropy

We study the diffusive motion of metastable localized spin structures with various topological charges by combining atomistic spin model simulations with an extension to the Thiele equation taking into account translational and rotational degrees of freedom. The Brownian motion of skyrmions is shown to obey markedly different diffusion laws depending on their topological properties. We demonstrate that skyrmions with topological charges Q = −1 and Q = 0 exhibit anisotropic translational diffusion due to their non-circular shape. The coupling between rotational and translational Brownian motion, however, causes memory of the initial orientation to be lost on a timescale governed by the effective rotational diffusion coefficient, ultimately leading to isotropic diffusion in the long time limit.


Introduction
Magnetic skyrmions are localized spin configurations where the directions of the magnetic moments span the whole unit sphere [1]. They are small in size and can be displaced by applied currents, which can be significantly lower than needed for moving domain walls [2][3][4], and are thus promising candidates for future magnetic logic and memory devices [5][6][7][8].
Although skyrmions exist as metastable excitations in the two-dimensional Heisenberg model it was shown by Belavin and Polyakov [9] that the presence of a finite magnetic anisotropy leads to a collapse of the skyrmion. Therefore, an additional interaction term in the Hamiltonian is needed in order to stabilize the radius of magnetic skyrmions. So far, several mechanisms for such an interaction have been identified: the Dzyaloshinsky-Moriya interaction (DMI) [1,10,11], four-spin interactions [12] and a frustrated Heisenberg-type exchange interaction [13]. In past investigations the focus was mostly on materials with DMI, where inversion symmetry is required to be broken. After the first detection of a skyrmion lattice in MnSi [14], skyrmions have been found experimentally in a variety of other bulk systems [15][16][17][18][19][20] and in thin films [21][22][23][24].
In recent theoretical works [25,26] on a (Pt 1-x Ir x )/Fe bilayer on Pd(111) it was demonstrated that the interplay between DMI and the frustration of the isotropic exchange interaction is responsible for the formation and stability of isolated skyrmionic spin structures of various types. These spin structures are characterized by their topological charge Q = 1/4π d 2 r S × (∂ x S × ∂ y S) where S is a unit vector describing the direction of the magnetic moment. So far, skyrmionic spin structures with topological charges ranging from −3 to 3 have been identified [27]. Usually only the spin structure with Q = 1 is called skyrmion but in the following all localized, metastable spin structures with integer topological charge will be referred to as skyrmions. Note, that what we call the Q = −1 skyrmion is frequently referred to as antiskyrmion in the literature. Furthermore it was demonstrated that skyrmions with Q = 1 are rendered non-circular by the DMI [26] leading to anisotropic current-driven dynamics [28].
Recently, stochastic skyrmion motion has attracted considerable attention because of its potential application for probabilistic computing [29]. Analytical calculations based on a rigid body approach predict that skyrmions exhibit Brownian motion [30]. So far, this has been demonstrated numerically [30][31][32] only for Q = 1 skyrmions. Skyrmion Brownian motion has also been observed experimentally [32][33][34][35] and it was demonstrated that its strength can be controlled by voltage application [36] and that a field-induced symmetry breaking leads to anisotropic skyrmion diffusion [37].
In this paper, we discuss the diffusive motion of skyrmions with various topological charges in the model system mentioned above, a (Pt 0.95 Ir 0.05 )/Fe bilayer system on a Pd(111) surface. We perform atomistic spin simulations based on the stochastic Landau-Lifshitz-Gilbert (LLG) equation [38] using a model Hamiltonian for the system parametrized by ab initio calculations in [25,32]. We demonstrate that the diffusion is strongly dependent on the respective topological charge and that skyrmions with broken rotational symmetry exhibit anisotropic diffusive dynamics. Moreover we study rotational Brownian motion and find a transition from short-time anisotropic to long-time isotropic diffusion.
The numerical results are compared with theoretical predictions based on a rigid body approach: we derive an equation of motion for the rotational degree of freedom from a Lagrangian approach and obtain an analytical prediction for the time-dependent diffusion coefficient tensor that is evaluated numerically making use of the Smoluchowski equation.

Model
As an example for a thin, ferromagnetic film with frustrated exchange and DMI, we model a (Pt 0.95 Ir 0.05 )/Fe bilayer on a Pd(111) surface via a spin Hamiltonian H considering only the magnetic Fe moments, The tensorial exchange coefficients J ij , the on-site anisotropy tensor K and the magnetic spin moment μ s were taken from [25,32]. B is an external magnetic field which is applied perpendicularly to the surface and has a fixed value of B ⊥ = 0.5 T throughout this study. The systems ground state is a spin spiral state which transforms into a collinear state upon applying a field of B ⊥ 0.21 [26]. Skyrmions with different topological charges Q can then occur as metastable excitations. Skyrmions with topological charges ranging from −3 to 3 have been identified. For the sake of simplicity, however, the following discussion is restricted to skyrmions with Q = −2, −1, 0, 1 which are shown in figure 1. While in this system the DMI, which is modeled via the antisymmetric part of the interaction tensor J ij , favors the formation of skyrmions with Q = 1, those with other topological charges are stabilized by the frustration of the Heisenberg-exchange (the diagonal part of J ij ). However, they are deformed by the DMI [26] and only the skyrmion with Q = 1 keeps its circular shape, while the other skyrmions become non-circular.
Furthermore it was discovered that in this system the DMI induces an orientation-dependent potential, due to the different energies of domain walls along different crystallographic axes [26]. This potential leads to preferred orientations with respect to the underlying lattice for skyrmions with Q = 1. We find that the shape of the potential can be well approximated via U(θ) = ΔU/2(1 − cos Mθ) with M = 3 for Q = 0 skyrmions and M = 6 for skyrmions with Q = −2, −1, similar to what was found in [39].
The dynamics of the spins is calculated by means of the stochastic LLG equation of motion [38], Here, γ is the absolute value of the gyromagnetic ratio and α is the Gilbert damping parameter. The effective field H i = −∂H/∂S i is generated by the external field, the on-site anisotropy and the interaction with neighboring spins. ζ i denotes a fluctuating, stochastic field which is included to model the effects of thermal fluctuations by creating random torques on the magnetic moments. Following [40], ζ i is defined by where μ and ν are Cartesian components, k B is the Boltzmann constant and · · · th denotes an average taken over different realizations of the thermal fluctuations.

Equations of motion
It is well established [30,31] that translational Brownian motion of skyrmions can be described by a Langevin-type equation of motion, which describes the motion of a localized rigid spin structure under the effect of thermal fluctuations. It can be obtained by supplementing the equation derived by Thiele [41] with a fluctuating, stochastic force F, the autocorrelation of which is given by . G is the gyrocoupling vector which is perpendicular to the film and its absolute value is given by G = 4πQμ s /(a 2 γ) where a is the lattice constant. D is the dissipation tensor, whose entries are given by In contrast to the gyrocoupling vector, which only depends on the topological charge, the dissipation tensor depends on the details of the spin structure. From its definition it follows immediately that the dissipation tensor is symmetric, i.e. D xy = D yx , and consequently D μν can be diagonalized by rotating the coordinate system where the diagonal entries are given by its eigenvalues, denoted by D a and D b . As their cylindrical symmetry is broken, skyrmions with Q = 1 possess an additional degree of freedom, namely their orientation, in the following described by the orientation angle θ. Using an ansatz similar to [39] we derive in appendix A an equation of motion for θ which is given by where is the rotational dissipation coefficient, U is the orientation-dependent potential as discussed above and F θ is a stochastic torque which was added to the equation of motion accounting for thermal fluctuations. The torque autocorrelation can be obtained analogous to the translational Brownian motion as shown in [31] and follows as F θ (t)F θ (t ) th = 2αk B TD θ δ(t − t ). As equation (7) is derived assuming a rigid spin structure it lacks an inertia-term, i.e. a term that depends onθ. In the following section, we demonstrate that this term is indeed missing and that it has a strong impact on the dynamics in the low α regime. For Q = −1, equation (7) in the absence of thermal forces is equivalent to a previously derived equation describing the helicity dynamics [39]. We further want to emphasize that equation (7) does not hold for Q = 0 skyrmions: for those, the equation describing their orientation contains additional terms that couple the rotational and translational velocities. Consequently the derivation of the time-dependent diffusion tensor in the following subsection does not apply to Q = 0 skyrmions. However, this does not pose a problem because our simulations show that rotational motion of those spin structures can be neglected. This is because the energy barrier for their destruction is significantly lower than the orientation-dependent potential making it impossible for them to rotate. We still simulate their translational diffusive motion and compare the results with theoretical predictions based on equation (5).

Time-dependent diffusion tensor
Equations (5) and (7) are coupled Langevin equations describing the motion of a rigid elliptical particle in a periodic rotational potential. The coupling between translational and rotational motion is due to the orientation-dependence of the dissipation tensor, i.e. D μν = D μν (θ). Similar to the derivation of a time-independent diffusion tensor in earlier works [30,31,37], an expression for the time-dependent diffusion tensor D μν can formally be derived from its connection to the mean squared displacement For a given initial orientation θ 0 it is given by (details of this derivation can be found in appendices B and C) which describes the transition from a possibly anistropic regime (if D a = D b ) to a regime, where memory of the initial orientation is lost and hence diffusion becomes isotropic. This transition has been extensively studied in the case of ellipsoidal colloids [42]. Assuming θ 0 = 0, equation (9) reduces to previously obtained expressions in the isotropic limit (long timescales) [30,31] and anisotropic limit (short timescales) [37] The timescale τ * which separates the two regimes can be calculated using the effective rotational diffusion The effective rotational diffusion coefficient can be calculated using the free rotational diffusion coefficient D θ = k B T/αD θ and a correction coefficient η[βU] that depends on the potential and on temperature (β = 1/k B T). A detailed discussion of the derivation of equation (12) and the calculation of η[βU] can be found in appendix C. For a periodic potential it can be proven that 0 η 1 [43] which is directly related to the fact that the potential leads to an increased dwell time of a Brownian particle in a minimum of the periodic potential. We obtained the effective rotational diffusion coefficient D * θ by solving the Smoluchowski equation [44] corresponding to equation (7).

Simulation results
The numerical integration of the stochastic LLG equation has been performed by means of a GPU-based implementation of Heun's algorithm [38]. We simulate 512 × 512 spins with a fixed timestep Δt = 1.4 × 10 −15 s and periodic boundary conditions. The trajectory Δx i (t) of a skyrmion was determined by monitoring the change in the magnetization which has proven robust against thermal fluctuations. The diffusion tensor D was calculated using its connection to the MSD given in equation (8).
The simulation results were compared with the theoretical predictions evaluated using the numerically calculated parameters D, G and D θ for skyrmions at rest and T = 0 K. These values are summarized in table 1. We found that the skyrmions with Q = −1, 0 have an anisotropic dissipation tensor (D a = D b ) in contrast to the skyrmions with Q = 1, −2. This anisotropy is caused by the DMI and was discussed in detail in reference [28]. As a consequence, we expect anisotropic diffusion only for skyrmionic structures with Q = −1, 0. Moreover it was found, that the Q = 0 spin structure can only overcome the rotational potential at temperatures that exceed its range of thermal stability and hence, the discussion of the anisotropic-to-isotropic transition is restricted to the Q = −1 skyrmions.

Dependence on topological charge and temperature
First, we will discuss the influence of the topological charge and temperature on the diffusive motion of the considered spin structures. Hereby, we neglect for a moment a possible orientation-dependence and only Table 1. Eigenvalues of the dissipation tensor D a and D b and rotational dissipation coefficients D θ calculated numerically for different types of skyrmions at rest as shown in figure 1. consider the total diffusion constant D tot = Tr[D] = D x + D y . From equation (9) it can be directly concluded that the total diffusion constant is time independent and given by The above expression predicts that diffusion decreases with increasing topological charge which implies, as long as the dissipation tensors for different skyrmions do not differ drastically, that in the system considered here skyrmions with Q = −2 exhibit decreased while the Q = 0 skyrmions exhibits increased diffusive motion as compared to the skyrmions with Q = ±1. This was verified by our simulations as depicted in figure 2 where the total diffusion constant D tot at T = 1 K is shown. Overall we found good agreement between theoretical predictions (dashed lines) and simulation results (symbols). D tot for skyrmions with Q = ±1 are not exactly equal, due to the difference in the respective dissipation tensors. Moreover, the results depicted in figure 2 clearly illustrate the qualitatively different dependence on α: while the diffusion of topologically trivial structures, for which Q = 0, is inversely proportional to friction (αD describes friction in equation (5)), the gyrocoupling leads to linear dependence on friction for structures with finite topological charge for α → 0. Figure 3 depicts the dependence of D tot of Q = ±1 skyrmions on temperature at α = 0.1. The two insets correspond to snapshots of profiles of a Q = 1 skyrmion at T = 1 K and T = 30 K, illustrating the impact of thermal fluctuations on the profile. We found that at low temperatures (T 12 K for the Q = −1 skyrmion and T 15 K for the Q = 1 skyrmion) the simulation results agree very well with the theoretical predictions (dashed lines) and yield the expected linear dependence of the diffusive motion on temperature with an increasing difference between the obtained values for D tot for Q = ±1 skyrmions as a result of the different dissipation tensors. However, at temperatures above that threshold, deviations from that behavior are clearly visible and we found that theory tends to underestimate the diffusion constant for both types of skyrmions shown in figure 3. As discussed in [31] this can be attributed to the heavily deformed profile which is not accounted for in the theoretical predictions, as the Skymion profiles at zero temperature are used to calculate the dissipation tensor D. The deviations from the theory are higher for the Q = −1 skyrmion, suggesting a stronger impact of thermal fluctuations on its spin configuration, which results in a decrease in the difference between the obtained values for D tot for Q = ±1 skyrmions at high temperatures.
As a naive approach, one could use a finite temperature profile in equation (6) but this leads to non-vanishing diagonal elements of the dissipation tensor even in the absence of a skyrmion. This is due to the fact that the calculation of the diagonal elements of D involves an integration of the squared derivative of the profile which is strictly positive. But there is no clear connection of a thermal fluctuation in the vicinity of a skyrmion and the skyrmion velocity, as required in the derivation of equation (5) and hence,  there is no straightforward way of generalizing this equation by introducing temperature-dependent parameters.

Anisotropic diffusion and anisotropic-to-isotropic transition
The timescale τ * , at which the anisotropic-to-isotropic transition occurs, is inversely proportional to the effective rotational diffusion coefficient given in equation (12). For free rotational diffusion, the correction factor η = 1 and thus τ * is inversely proportional to temperature. In the presence of a potential however, η depends on temperature and U(θ) in a non-trivial way. In [39] it was shown that the rotational potential can be approximated by a cosine function. Using the Smoluchowski equation under this assumption we show that the correction factor can be well approximated via η(βΔU) = (I 0 (βΔU/2)) −2 , where I n (x) are the modified Bessel functions of the first kind (details in appendix C).
Together with equation (12) one can then obtain ΔU by fitting to the mean switching time τ sw ∼ τ * obtained from simulation data, which is the mean time to get from one minimum of U(θ) to another (figure 4), yielding ΔU/k B ≈ 8.5 K. In order to ensure consistency we also calculated the potential directly by using the relation U(θ) = −k B T ln p(θ) [45], with p being the probability of finding a skyrmion with orientation θ, which was calculated from simulation data. The resulting potential for T = 5 K is shown in the inset of figure 4. It indeed has a cosine shape and the potential strength follows as ΔU/k B ≈ 9 K, which compares well with what was obtained above. Note, that this value is specific to the Q = −1 skyrmions and is not necessarily equal for the other skyrmions.
From figure 4 we concluded that at T = 1 K the skyrmions are unable to overcome the rotational pinning in the timescale of our simulation ( 2 ns) and hence we deemed this value of T suitable for the study of the anisotropic diffusive regime. The corresponding simulation results for the diffusion constants D x and D y , the diffusion constants along x-and y-axis, are depicted in figure 5 for skyrmions with Q = 0, −1 initially oriented as shown in figure 1. Our results clearly verify that both types of spin structures   (14). The red dotted line is a fit of equation (14) to the data.
Upon increasing temperature, τ sw becomes comparable to the simulation time and thus, one can certainly expect deviations from the anisotropic limit of diffusive motion, ultimately leading to isotropic diffusion in the long-term limit. This anisotropic-to-isotropic transition for α = 0.1 and T = 3 K is depicted in figure 6 for the Q = −1 skyrmion. Figure 6(a) shows the MSD along x-and y-axis and the dotted lines correspond to linear fits for t 0.1 ns (illustrated by the gray area) which serve as a guide to the eye. It is clearly visible that the two lines bend toward each other and away from their initial slopes. In order to analyze this transition in a more quantitative way, one possibility would be to calculate the ratio of diffusion coefficient. However, this poses the problem that this ratio only converges to 1 in the limit t → ∞, even though the slopes of the MSD curves are parallel, which indicates that diffusive motion along both axes is equal. Hence we concluded that it is preferable to calculate the ratio of the slopes of the MSD curves (for a general, time dependent diffusion coefficient 2D(t) = ∂ t MSD). In an effort to smooth out numerically obtained derivatives of the MSD, this was done by performing linear fits within intervals with a fixed length of 0.1 ns starting every 0.02 ns. The ratio of the resulting slopes, i.e. ∂ t MSD y /∂ t MSD x is shown in figure 6(b) (blue circles) exhibiting a transition from the anisotropic limit (D a /D b ) to isotropic diffusion. The blue dashed line corresponds to the theoretical prediction which follows from equation (9), if assuming that θ 0 = 0, as The above expression was evaluated by using equation (12) and the parameters listed in table 1. The correction factor was calculated using the Smoluchowski equation with a potential strength of ΔU/k B = 8.5 K, as estimated above, and yields η ≈ 0.28 (details in appendix C). Clearly, the prediction based on our theory agrees poorly with the simulation. The reason for this is not that equation (14) is incorrect, but rather that the value D * θ is clearly overestimated: the red dotted line corresponds to a fit of equation (14) to the simulation data yielding a roughly 4 times lower value.
By calculating D * θ from the simulation data for several values of α ranging from 0.01 to 1.0 (shown in figure 7) following the above procedure we found that accuracy of the theoretical prediction depends on α. For high α, we found quantitative agreement while for values lower than 0.4, there are large deviations from the theoretically expected inverse proportionality to α as instead the effective rotational diffusion goes to zero for vanishing damping. This discrepancy between simulation data and theoretical predictions can be attributed to the fact that an inertial term is missing in equation (7). Similar to the superparamagnetic switching process (see for example [46]) or the Brownian motion in an external periodic potential [47] an inertial term has a strong impact in the low damping regime. In the high damping regime on the other hand, the friction term dominates over an inertial term and the latter can be neglected. Thus, equation (7) can be understood as the high-damping approximation to the correct equation of motion for the orientation. In order to treat the full dynamics correctly, one would have to extend the derivation in appendix A and then solve the corresponding Fokker-Planck equation, which is not the Smoluchowski equation but the Klein-Kramers equation [48].

Conclusion
We examine the diffusive motion of metastable skyrmions with a variety of topological charges in a (Pt 0.95 Ir 0.05 )/Fe bilayer on a Pd(111) surface by the means of spin model calculations based on the stochastic LLG equation.
In agreement with earlier analytical predictions [30] we find a linear dependence on temperature and a strong topology-dependence of the diffusive motion: while structures with finite topological charge exhibit decreasing diffusion in the limit α → 0, topologically trivial structures show an inverse proportionality to damping. Furthermore our simulation results reveal that skyrmions with topological charges 0 and −1 exhibit anisotropic diffusion, i.e. increased diffusive motion along one axis of their profile and reduced along the perpendicular axis. As a consequence, we observe a transition from an anisotropic regime to a regime, where memory of the initial orientation is lost and hence diffusion becomes isotropic.
In order to compare our numerical findings with analytical predictions based on a rigid body approach we extend the existing theory on skyrmion diffusion by adding the rotational degree of freedom for structures with broken rotational symmetry. In general, we find good agreement between theoretical predictions and numerical findings. We observe quantitative agreement for translational diffusion in the anisotropic and the isotropic limits as well as for the rotational diffusion in the high damping limit. However, if damping is small, we find large deviations from the theoretically expected inverse proportionality to α. These deviations can be attributed to the absence of an inertia term in the equation of motion for the orientation of skyrmions.
Our findings indicate that the strength of Brownian motion of skyrmions depends strongly on their topological charge. This could be used experimentally to determine their topological properties, similar to measurements of the Q dependent skyrmion Hall angle [28]. Especially, this dependence could serve as method to distinguish between structures with finite topological charge and topologically trivial structures since the latter are expected to be subject to stronger diffusive motion. Moreover, if the topological charge can be changed experimentally, one could use this to switch between high and low diffusive motion, which could for example be used in skyrmion reshuffler devices which are essential for the realization of probabilistic computing [32]. Furthermore, our results demonstrate that deformations of the spin configurations can give rise to anisotropic Brownian motion. Here, these deformations are intrinsic to the system via its interaction parameters. In reference [37], however, it was shown that deformations can also be induced by an external stimulus, namely via an external in-plane field, which gives rise to anisotropic diffusion with higher diffusive motion along the field direction than perpendicular to it. This opens another pathway for the efficient manipulation of skyrmion Brownian motion.
Following equation (6), the dissipation tensor is symmetric and real and hence can be diagonalized using rotation matrices. Inversely, the dissipation tensor for a skyrmion with any orientation angle θ can be expressed in terms of its eigenvalues D a and D b and rotation matrices as After a simple matrix inversion, equation (5), in components form, reads as Using the force autocorrelation, which is given by (derived in reference [31]) the velocity autocorrelation function for fixed θ follows as Lastly, performing the orientational average · · · θ and then time integration, where Δx The first, time-independent term describes isotropic diffusion and was derived earlier [30,31].
The calculation of c 2 (t) requires the solution of the Smoluchowski equation [44] corresponding to equation (7) which is shown in the next section.

Appendix C. Effective rotational diffusion coefficient
The orientation dynamics are described by which corresponds to an overdamped Langevin equation as inertia terms are absent. The corresponding Fokker-Planck equation in the overdamped limit is the Smoluchowski equation with the free diffusion coefficient D θ = k B T/αD θ and the potential U(θ) which was found to be well approximated by a cosine as U(θ) = ΔU/2(1 − cos Mθ) with M = 3 for skyrmions with Q = 0 and M = 6 for Q = −2, −1. As skyrmions with Q = 1 are rotationally symmetric, this potential is absent for them. Note, that also the strength of the potential ΔU depends on the type of skyrmions.
where we introduced a dimensionless potential strengthŨ = MβΔU/4. In what follows, we assume that p(θ, t = 0) = δ(θ), where δ(θ) is the Dirac delta function which is symmetric at θ = 0. Therefore we can expand the probability density in a cosine basis as p(θ, t) = ∞ m=0 c m (t) cos mθ with c 0 (t = 0) = 1/2π and c m (t = 0) = 1/π in order to fulfill the initial condition. Inserting this expression into equation (C.3) leads to Multiplying with cos nθ and integration over θ yields, by using the fact that 2π 0 cos(nθ) cos(mθ)dθ = πδ nm for n, m > 0, a coupled set of linear differential equations for the c n (t), which can be summarized asċ with n = 0, 1, 2, . . . , or in matrix notation for M = 3 and similarly for M = 6. An exact calculation of c n (t) would require an exact diagonalization of the above matrix (without the first row and column) which is impossible, as it is infinite. However, the diagonal elements scale with n 2 whereas the off-diagonal elements increase linearly with n. Thus the relative coupling strength between different modes decreases with 1/n and becomes negligible above a certain n max . In other words, for n > n max the corresponding eigenvalue λ n can be well approximated by n 2 D θ . Thus the problem is reduced to the diagonalization of a n max × n max matrix. Moreover, we are only interested in the time evolution of c 2 (t), as equation (B.9) only depends on that single mode. Formally, where d (n) 2 is the corresponding element of the nth eigenvector and λ n is the nth eigenvalue of the matrix on the rhs of equation (C.6). From the fact that lim t→∞ p(θ, t) ∼ exp[−βΔU/2(1 − cos Mθ)] it follows that c n (t → ∞) = 0, ∀ n/M / ∈ N 0 and hence d (0) 2 = 0. Physically, this means that only a potential with M ∈ 1, 2 can lead to anisotropic diffusion at all times. The timescale at which c 2 (t) goes to zero is governed by the smallest eigenvalue λ min in equation (C.7), which can be expressed in terms of an effective diffusion coefficient D * θ as λ min = 4D * θ . Thereby, D * θ = D θ η[βU] ( C . 8 ) where η[βU] is a correction factor to the free rotational diffusion coefficient which depends on the potential. The factor 4 was chosen as such that λ min (Ũ → 0) = 4D θ , which immediately follows as the damping constant for c 2 (t) in the absence of a potential, and consequently 0 η 1. Lastly, we approximate c 2 (t) ≈ c 2 (0) exp(−4D * θ t) and together with c 2 (0) = 1/π we can rewrite equation (B.9), finally leading to the expression for the time-dependent diffusion tensor given by equation (9). The correction factor η was obtained by numerical diagonalization of the n max × n max matrix as discussed above and is shown in figure 8 for M = 6. It monotonically decreases to zero with increasingŨ following roughly the Lifson-Jackson (LJ) formula [50] η LJ = 1 e −βU(θ) e βU(θ) (C.9) where · · · = 1/2π 2π 0 (· · · )dθ. The LJ formula describes effective diffusion as long as ΔU is not too high. For the cosine potential considered here, it yields η LJ = (I 0 (βΔU/2)) −2 , irrespective of M, where I n (x) are the modified Bessel functions of the first kind. The inset of figure 8 shows the dependence of η on n max depicting quick convergence. Hence, the assumption that the off-diagonal elements for n > n max in the matrix in equation (C.6) can be neglected in the calculation of c 2 (t) is justified.