Simulation and Theory of Classical Spin Hopping on a Lattice

: The behavior of spin for incoherently hopping carriers is critical to understand in a variety of systems such as organic semiconductors, amorphous semiconductors, and muon-implanted materials. This work speciﬁcally examined the spin relaxation of hopping spin/charge carriers through a cubic lattice in the presence of varying degrees of energy disorder when the carrier spin is treated classically and random spin rotations are suffered during the hopping process (to mimic spin–orbit coupling effects) instead of during the wait time period (which would be more appropriate for hyperﬁne coupling). The problem was studied under a variety of different assumptions regarding the hopping rates and the random local ﬁelds. In some cases, analytic solutions for the spin relaxation rate were obtained. In all the models, we found that exponentially distributed energy disorder led to a drastic reduction in spin polarization losses that fell nonexponentially.


Introduction
Spin relaxation of carriers has garnered much interest in the past few decades due to the prospect of spin electronic or spintronic devices that use spin functionality in addition to charge functionality. However, most work has been concentrated on band transport [1][2][3][4][5]. Less effort has been expended on the spin relaxation of localized electron spins where charge transport occurs via incoherent hopping motion. Noncrystalline, organic semiconductors have recently seen a surge in interest [6][7][8][9][10][11][12] as spintronic candidates due to their small spin-orbit coupling and expected long spin lifetimes of 10 −5 -10 −7 s [13]. Within the field of spin chemistry, spin relaxation plays a role in radical pair reactions [14], which has applications varying from organic electronic devices [15] to avian navigation [16]. Localized spin relaxation also plays a role in transport through disordered crystalline semiconductors [17,18] Spin relaxation of hopping spins has an older genesis: it was first studied in relation to muon spectroscopy [19][20][21][22][23][24]. The technique has seen much application in the study of magnetic materials. Muon spin spectroscopy is used to study kinetic isotope effects in chemical reactions [25]. The method involves positively charged polarized (spin 1/2) muons, µ + , being implanted into a material. The interaction of the muon spin with local magnetic moments allows the technique to be an excellent probe of magnetic properties. The muon's lifetime is only 2.1970 µs, after which it decays into two neutrinos and a positron. The utility of the mechanism is due to the positron being emitted preferentially along the direction of its parent's spin. By detecting emitted positrons, the spin evolution of the implanted muon is revealed. When that evolution is impacted by spin interactions with the host material, temporal muon decay measurements may yield information about the host material; time scales that can be probed by this technique range from picoseconds to microseconds. Due to the positron emission being sensitive to magnetic interactions, probing spontaneous magnetic fields in magnetic materials was an early application of muon spectroscopy [26]. In dielectrics, semiconductors, and organic compounds, a lowenergy muon can form a stable state with an electron: muonium. Muonium offers access to hydrogen-like states in these materials and can be especially beneficial in the study of chemical reactions where muonium serves as a radical [27].
Information from muon spectroscopy is obtained from the measured polarization of the implanted muons. For static fields, stationary muon spins experience a distribution of environmental magnetic fields whose ensemble average yields the measured response. When muons diffuse throughout the host, the fields fluctuate, which alters the spin decay function of the ensemble [24]. Calculations of the polarization in time yield information on the fields, as well as the hopping diffusion of the muons. Figure 1a shows how the classical spin rotation is coherent at each location, but then abruptly changes when an incoherent hop occurs; at the new location, the evolution is again coherent, and the process repeats until the spin polarization decays to zero. Since the spin is rotated while situated at a site, we call this an intrasite process of spin relaxation. The intersite process involves the spin rotating during the hopping process itself, while the spin is stationary when sitting on a site. Intrasite processes dominate the polarization functions of muon spectroscopy and therefore have been well studied for many decades [19,28,29]. Previous studies have shown that quantum and classical calculations of spin diffusion length agree well [8,30]. The time to hop, 1/k, varies widely from system to system. As discussed later, the time is strongly dependent on energetic disorder. Times can be as short as femtoseconds [31] Intersite spin evolution, depicted in Figure 1b, occurs when the spin experiences some rotation (or flip in a quantum picture) during the transition between two localizing centers. This type of mechanism is similar to the process of spin relaxation by spin-orbit interactions in organic semiconductors [12,[32][33][34][35][36]. Here, we investigated in depth a classical model of spin relaxation from a spin-orbit-type interaction.
rotate by γ 2 Figure 1. (Color online) (a) Intrasite spin evolution where the precession angle (ωt) depends on the wait time at the site (t i at the ith site) and the random local field. Arrow shading represents the chronology of spin evolution: the oldest times are represented by lighter grays with shade darkening to the current time, t, which is black. (b) Intersite spin evolution, which is independent of the wait time. At each hop, the spin undergoes a spin rotation, γ, due to a random local field.

Theory
In this section, an analytic theory of classical spin relaxation is reviewed for intrasite mechanisms [34].

Spin Evolution
The unit vector, S, is a classical electron spin, which in an arbitrary uniform, static magnetic field, B, has a precessional frequency −ω = −|γ e |B/2 ≈ µ B Bn/h, where γ e = gµ B /h, n is the unit vector of B, and the Lande factor for spin is taken to be two. The evolution of S is described by the following equation: where Ω is the skew-symmetric matrix: Intrasite spin relaxation (hopping independent) was not considered here. The solution to Equation (1) is: where S 0 is the initial spin vector and R is the following rotation matrix: which is equivalent to this vectorial expression [37]: If ω = 0, there is no spin evolution at a site, soR(t) =1 and S(t) = S 0 . Spin relaxation occurs when carriers hop between different sites with random local magnetic fields. This was explored in the context of hyperfine interactions and inhomogeneous g-factors in previous articles [12,34]. Intrasite spin evolution will be ignored from hereon out.
The focus here is on an alternate spin relaxation mechanism where spin rotations occur between hops-the so-called intersite process of spin relaxation. The spin changes abruptly at each hop instead of coherently evolving while fixed at a localizing center. The spin, after a single hop, is updated as: where:T γ =1 + sin γΩ + 2 sin 2 γ 2ΩΩ (7) or: The process can be repeated to model the evolution of the spin. For a single spin, the spin will randomly walk on the unit sphere as the spin hops to various sites with random γn fields. Later in this article, we will simulate this process for an ensemble of spins. However, analytic theories can be devised for the ensemble if appropriate approximations are made. This is the subject of the next section.

Strong Collision Approximation
We used an approximation known as the strong-collision approximation [38]. The approximation is characterized by the following assumptions: (1) local field changes are abrupt and not slow at each hop (non-Gaussian); (2) local fields during each hop are uncorrelated from any previous hop (Markovian).
With the strong-collision approximation, the theory of continuous-time random walks [39,40] can be used to calculate the polarization functions of hopping carrier spins. Two important functions of continuous-time random walk theory are the survival probability, Φ(t), and the wait time distribution, ψ(t). Φ(t) is the probability that a random walker starting on a site at t = 0 has not hopped by time t. ψ(t) is the probability distribution from which the hopping times are chosen. The two quantities are related by: The polarization of spins that have not hopped is then P 0 (t) = Φ(t)S 0 , since only those spins that have "survived" on their initial site contribute the polarization of stationary spins.
The polarization of spins that have hopped exactly one time by time t is P 1 (t) = R 1 (t)S 0 , where: whereR 1 (t) is the rotation matrix for the sub-ensemble of spins that have rotated one time by time t. This equation is interpreted term-by-term going from right to left: the probability of making a single hop in the interval dt is ψ(t )dt ; the rotation of the spin as it hops one time, due to the intersite process, isT γ,n , where n is the nth hop in the sequence of hops; once the hop is completed, the spin remains at the new site according to the survival probability. Since there is an ensemble of spins of which many will hop a single time within the interval t, we must integrate over the possible hopping time dt . Lastly, the angular brackets denote an average over the angular coordinates of theT γ,1 rotation matrix Equation (7), which yields: The magnitude of the rotation angle, γ, is held constant for the time being. The cumulative rotation matrix for spins that have hopped twice is: where the angular averaging of the two rotations,T γ,1 andT γ,2 , is independent since there is no correlation between local fields; the two average rotation matrices are also equal; we then denote their average as T γ,1 = T γ,2 ≡T γ . The polarization of these two-hoppers is P 2 (t) =R 2 (t)S 0 . The generalization to higher order hops is straightforward, which leads to a total polarization P(t) = ∑ ∞ n=0R n (t)S 0 . The net polarization function becomes tractable by Laplace transformations and the convolution theorem.
Analytic solutions are obtained by Laplace transforming the rotation matrices. For instance: feed intoP(s) = ∑ ∞ n=0R n (s) · S 0 , which is a geometric series summable to:

Non-Disordered Model
A simple example of the theory is demonstrated by using a exponential wait time distribution, ψ(t) = ke −kt , where k is the average hopping rate. The relevant quantities feeding into Equation (14) areψ(s) = k/(k + s) andΨ(s) = 1/(k + s), which yield first: which is reduced to:P The desired polarization function in time is ascertained by the inverse Laplace transform to be: When the spin rotation angle γ is small, the polarization function is: A comparison of the full polarization function, the approximate polarization function suitable for small rotation angles, and the multiple trapping simulations (to be described later) are displayed in Figure 2. The relaxation rate's kγ 2 dependence can be correlated with a quantum spin model for spin-orbit-induced spin relaxation, where γ is interpreted as a spin admixture parameter instead of a rotation angle. First principle calculations of γ for a common organic semiconductors yield 10 −2 -10 −4 [32].

Hopping Models
The spin relaxation of hopping spins is modeled in two ways. The first is Multiple Trapping (MT), where each hops is completely uncorrelated from previous hops. The advantage of this model is the ability to obtain analytic solutions. The disadvantage is that in a real system, hops are in fact correlated, which means the random fields giving rise to relaxation are correlated. The Multiple Hopping (MH) model accounts for these correlations. We used Monte Carlo methods to simulate both models.

Multiple Trapping Model
Realistic wait time distributions, appropriate for amorphous and organic semiconductors, are harder to handle than what was treated in Section 2.3. The two aforementioned disordered systems are typically characterized by an exponential density of states and a Gaussian density of states, respectively. The width of the Gaussian is typically taken to be on the order of tenths of an electron volt in organic semiconductors [41]. The Multiple Trapping model (MT) constructs wait time distributions ψ(t) = k(ε)e −k(ε)t from hopping rates that are of the form of trap release rates, k(ε) = k 0 e ε/k B T , where ε is the trap energy and is distributed according to the density of states. The energy levels are uncorrelated between hops. The strong-collision approximation is still assumed, so there is no correlation in rotational fields either. For an exponential density of states, g e (ε) = (k B T 0 ) −1 exp(ε/k B T 0 ), ψ(t) is described by 0 −∞ g e (ε)k(ε)e −k(ε)t dε or in Laplace space as: with x = ε/k B T 0 and α = T/T 0 . The result can be written as a hypergeometric function: The long time, or asymptotic, behaviors of ψ(t) and Φ(t) are of interest; they are determined to be [34]: The polarization at any time can be ascertained by the numerical Laplace transform inversion [42] of Equation (14) [43,44]. However, in the long-time case, asymptotic analysis determines [34]: which demonstrates an algebraic decay of spin polarization. n = 3 for γ distributed with a constant magnitude over a sphere and n = 1 for a randomly directed Gaussian γ. For a Gaussian density of states, g g (ε) = (2πσ 2 ) −1/2 exp(−ε 2 /2σ 2 ), the problem is less analytically and computationally tractable as neither pure exponential nor power law behavior are observed for finite σ/k B T [8,45].

Multiple Hopping Model
The MT model treats each hop as an independent trap release with release rate k(ε) = k 0 e ε/k B T . MT ignores correlations between sites, which would cause the hopping to not follow the carrier release from distribution sites with some density of states. For example, two sites lying near in energy and space would allow rapid back-and-forth hopping before the charge carrier escapes the dyad. Therefore, MT is inappropriate in general, but is useful due to its simplicity and the availability of analytic solutions to some problems.
To address more realistic hopping, we adopted a Monte Carlo approach where charge hopping and concomitant spin evolution are tracked on a cubic lattice where the rate of hopping from site i to site j is determined by the Miller-Abrahams rates: Each site on the lattice is assigned a fixed energy drawn from the chosen density of states (either the delta function, exponential, or Gaussian). The field vector responsible for rotating the spin is γn. Later on, different assignations for γ are discussed. Spin carriers were injected randomly onto the lattice of sites. The spin of each carrier was sampled at a chosen time interval and averaged over many different configurations of the disorder (typically 10,000-50,000). A single disorder configuration possesses a fixed landscape of site energies and local fields.
Various models for the local fields, γ, can be employed. Below, we describe two models for the correlation of the fields between sites (applicable only to MH) and two models of the local field vector (applicable to either MT or MH).

Local Field Correlations
The previous site is always "forgotten" after each hop in the MT model, so models of site correlations only apply to the MH model. Two models were employed: the vertex and edge local fields models.
Vertex local fields: Each site possesses a single local field vector γ, which exerts a torque on the spin hopping to that site, as shown in Figure 3. This model was employed in [12,34].
Local Field (γ) Edge local fields (commuting): The local field exists between two nearest neighbor sites. The local field for a forward hop is opposite that of a backward hop. Therefore, a spin that undergoes a forward-backward hop combination ends up with its initial spin state as shown in Figure 4. This model mimics spin-orbit interactions for sub-barrier carrier motion [35].

2
Local field gamma ( ) 1 → 2 2 → 1 Figure 4. (Color online) Edge model of local fields. A spin rotates a set amount γ about a random direction (green arrows) when hopping to a site (e.g., a hop from Site 1 to Site 2). The return hop (e.g., a hop from Site 2 to Site 1) rotates through the opposite direction (red arrows) an amount γ. Therefore, the net rotation for a hopping sequence 1→2→1 is zero.

Local Field Magnitude
Two different models of the local field magnitude were employed (demonstrated in Figure 5): (1) the magnitude was constant (the direction was random); (2) the magnitude was drawn from the Gaussian distribution (the direction was again random).
Constant local fields: Each local field has a fixed magnitude, but a random direction, such that the rotation angle is Γn, where Γ is constant. This leads to a spin relaxation function P(t) = e − 1 3 kΓ 2 t S 0 . Gaussian local fields: If each component of the rotational vector Γn is isotropically, normally distributed with width γ, then in the small-angle approximation for Γ, the relaxation rates change to γ 2 k such that P(t) = e −kγ 2 t S 0 .
Details on the implementation of these models are described in Section 5.

Results
Now, the results from the aforementioned simulations are presented. First, the results using zero disorder (σ = 0) are shown, and then, the results for finite disorder (using an exponential density of states) are displayed.

Zero Disorder
Simulations of MT and MH using each of the previously described models are demonstrated in Figures 6 and 7. Exponential fits for nondisordered simulations were carried out with the function: where P 0 and α serve as fitting parameters. Table 1 shows the values of α obtained for each model. The MT values were derived from the theory presented earlier in this article.
In view of all the MH simulations in the figures below, one can readily understand why the edges model led to slower spin relaxation: since to-and-fro hops commuted, such hops led to no spin losses where every hop in the vertex model had the potential of incurring a random spin rotation.   (24)). Additional parameters used in the simulations are γ = 0.025 (for the Gaussian approximation, γ is the width of the distribution) and k = 1.
An interesting finding demonstrated in Table 1 is that the MH-edges rate is found to be 2/3 the MT rate (for both field magnitude models). We can understand the factor of 2/3 in the following way. In MT, all six hopping directions lead to a spin rotation (and hence, contribute to spin relaxation). For the MH-edges model, only five hopping options contribute to spin relaxation; the sixth option is a return to the spin's previous site, which will actually suppress the spin relaxation due to the assumed commutation of back-andforth hops. By compiling the fraction of hops that lead to relaxation or reduce relaxation, on a cubic lattice, we have 5 6 (+1) + 1 6 (−1) = 4 6 = 2 3 . In general, with coordination number q, the relaxation suppression fraction is q−1 q (+1) + 1 q (−1) = q−2 q , which we confirmed for the square lattice.

Finite Disorder
Figures 8-10 demonstrate the role nonzero disorder (σ > 0) plays in the spin polarization functions. Figure 8 is the MT model with σ = 2. The thin black lines are the analytic results from inverse Laplace transforms of Equation (20). The linear dashed lines are the analytic approximation at longer times given by Equation (22). The key features, when σ > 1, are that the polarization function decays algebraically at long times no matter the model being used (however, an exponential density of states was assumed here).
The much slower algebraic relaxation (as opposed to exponential) was due to the highdisorder regime of σ > 1. Spin rotations only occurred during hops, and large disorder means that carriers would fall into deep energy states and take a long time to hop out. Since relaxation requires hops, the relaxation function follows the same trend as the survival probability function. σ = 2 c o n s t a n t γ g a u s s i a n γ  (14) and (20). The dotted lines are the expression Equation (22). Additional parameters used in simulations are γ = 0.025 (for the Gaussian approximation, γ is the width of the distribution) and k = 1000. The simulations are averaged over 60 k disorder configurations. Figure 9 shows a similar behavior for the MH-edges model, and we found the MHvertex model (not shown) to be qualitatively identical as well.

Discussion
The simulations presented above modeled the instantaneous rotation of a classical carrier spin as it hopped across a lattice. While our model was a classical approximation to a quantum mechanical, the motivation was that this type of model should qualitatively pick up some features of hopping-induced spin-orbit spin relaxation, which is believed to be responsible for the spin losses in some amorphous and organic semiconductors. In this section, we discuss the aspects and shortcomings of our classical model, as well as future directions this work will take.

Evaluating the Models
In terms of capturing the relevant physics, the MH model is superior due to the accounting of correlations during hops. The edges version of MH is most applicable to underbarrier motion undergoing a spin-orbit interaction [35]. The vertex version of MH is most suitable for overbarrier hopping transport, where one would not expect commuting spin rotations between adjacent sites.
Determining which of all the models presented here (MT, MH, vertex, edges, Gaussian, constant) are governing a particular system is difficult due to the uncertainty in the parameters (γ, k) for a given system. This can be seen for the following reasons: (1) if σ = 0, all models predict spin relaxation rates, αγ 2 k with α differing by order unity; (2) if σ > 0, all models predict algebraic spin decay ∼t −1/σ . The predicted ∼γ 2 k dependence for spin relaxation was used in [46] to fit experimental electron spin resonance measurements in undoped amorphous silicon [47,48] and hydrogenated amorphous silicon [49,50]. While these amorphous systems are disordered (which entails algebraic decay), the ∼ γ 2 k dependence for the line shape width may be due to resonance experiments picking up short-time exponential dynamics. A better comparison with the theory presented here, to probe the influence of strong disorder, would require experiments performed using time-resolved electron spin resonance.
Time-of-flight experiments clearly demonstrate the role of disorder in charge transport (so-called dispersive transport) [51,52]; however, extending this method to spin transport is nontrivial, and the authors are not aware of any experiments of this type.

Shortcomings of the Classical Spin Model
Consider a simple model atom in the 2p state, where p z has an energy ∆ below that of p x and p y . For the spin-orbit Hamiltonian, H so = ξL · S, perturbation theory yields doubly degenerate eigenstates: and a spin-orbit energy shift of −ξ 2 /∆. |+ and |− are pseudospin states: since ξ/∆ is small, the |+ pseudospin state is mostly | ↑ with a small component of | ↓ . Spin-orbit coupling mixes up-and down-spin spin states, which entails that spin is not conserved. The hopping integral between two sites i and j is spin independent, but since the pseudospin states |+ and |− are admixtures of up and down states, there is a finite probability for transitions between |+ and |− states during hopping events.
If one takes the local fields to signify spin-orbit interactions, we can contrast the classical model with quantum models of carriers hopping between sites experiencing spinorbit interactions [32,33,36]. Due to the spin-orbit Hamiltonian not commuting with the spin operators, spin is not conserved during hops. Typically, in an organic semiconductor (e.g., Alq3), the spin-orbit interaction is weak, such that many spin-conserving hops may take place before a spin flip. The classical spin model presented herein describes the spin relaxation as a more continuous process, where small spin rotations occur at each hop and the accrual of such rotations leads to relaxation. The quantum model as described by [32,33,36] is inconsistent with the underbarrier motion modeled by MH-edges [35]. The deviation of the classical MH model with the quantum model (for simplicity, we assumed no disorder) is starkest in lower dimensional hopping systems, where commuting hops are assumed by the edges model. The commutation of oppositely directed hops suppresses spin relaxation, while we do not expect such a dramatic suppression in the quantum spin model (since spin flips occur probabilistically, there would be no spin flip reversals when returning to a site). A detailed study of the quantum model awaits our further focus.

Materials and Methods
Calculations were carried out in Python, C++, and CUDA. These three languages were chosen in order to confirm each other and also to parallelize the computations. MH codes initialized a classical spin vector and then allowed the spin to randomly walk over a cubic (for 3D simulations) lattice (typically 80 3 , which proved to be far beyond where finite-size effects would be observed) according to the Miller-Abrahams hopping rates. The spin vector would be sampled every ∆t. Many thousands of spins were evolved in this manner, and their vector at every ∆t would then be averaged to find the average spin decay function. When disorder was present, each run would randomize the site energies according to the density of states. The MT code operated in a similar fashion, except that the activation energy at each hop was randomly drawn according to the density of states. In this manner, each hop is completely uncorrelated from the previous hops.

Conclusions
For the first time, we provided a systematic treatment of intersite classical spin relaxation processes for hopping carrier spins. The variety of models examined, while classical, carried enough features to give insight into what may occur for a hopping quantum mechanical spin experiencing spin-orbit interactions during the hopping processes. Given a specific amount of disorder, the classical models do not deviate enough from each other to adequately distinguish between them since important spin-orbit and hopping parameters are not known to high precision. Future efforts will be to study spin relaxation from hopping using a fully quantum mechanical spin-orbit interaction.