Microscopic theory of OMAR based on kinetic equations for quantum spin correlations
aa r X i v : . [ c ond - m a t . d i s - nn ] M a r Microscopic theory of OMAR based on kinetic equations for quantum spincorrelations.
A.V. Shumilin Ioffe Institute, 194021 St.-Petersburg, Russia
The correlation kinetic equation approach is developed that allows describing spin correlations in amaterial with hopping transport. The quantum nature of spin is taken into account. The approachis applied to the problem of the bipolaron mechanism of organic magnetoresistance (OMAR) inthe limit of large Hubbard energy and small applied electric field. The spin relaxation that isimportant to magnetoresistance is considered to be due to hyperfine interaction with atomic nuclei.It is shown that the lineshape of magnetoresistance depends on short-range transport properties.Different model systems with identical hyperfine interaction but different statistics of electron hopslead to different lineshapes of magnetoresistance including the two empirical laws H / ( H + H )and H / ( | H | + H ) that are commonly used to fit experimental results. PACS numbers:
I. INTRODUCTION
Hopping conductivity is one of the fundamental typesof electron transport in solid-state materials. It existswhen the electron wavefunctions are localized. The con-ductivity is achieved due to the acts of hopping whenthe electron hops between localized functions (sites) withdifferent energies due to the emission or absorption ofphonons. The conventional theory of hopping conduc-tivity is closely related to doped semiconductors withcompensation. It is based on the mean-field approxima-tion and Miller-Abrahams resistor network, which fol-lows from this approximation in weak applied electricfield [1, 2]. The drawback of the mean-field approxima-tion is that it neglects correlations between occupationnumbers of different sites.One type of materials with hopping transport that isactively developed right now is the organic semiconduc-tors. They are already widely applied in OLEDs [3] andhave other possible applications, for ex. in organic solarcells [4]. These materials very often display an intrigu-ing property that is called “organic magnetoresistance”or “OMAR” [5–11]. It is quite a strong magnetoresis-tance observed in magnetic fields 10 − gs both at lowand room temperatures. Although the qualitative expla-nations [12, 13] and semi-qualitative theories [14–20] ofthis phenomenon started to appear ten years ago, the de-tailed microscopic theory of OMAR is not yet developed.One of the reasons for it is the close relation of OMARto non-equilibrium spin correlations. The magnetoresis-tance is equal to zero in the mean-field approximation[21] but re-appears when the correlations are includedin the theory even in an oversimplified model [22]. Thephysical reason for OMAR is the dependence of the relax-ation of spin correlations on the applied magnetic field.This relaxation is often associated with hyperfine inter-action with atomic nuclei. With some simplification, itcan be described as spin rotation around the so-called“hyperfine fields”. These fields are different at differentsites therefore random hopping of electron with rotation around these fields leads to spin relaxation. When theexternal magnetic field is large compared to hyperfinefields, the spin rotates around approximately the samedirection on all the sites and its relaxation is suppressed.The microscopic theory of OMAR requires a theoreti-cal approach that takes into account the non-equilibriumcorrelations including the spin correlations. Up to veryrecent times, practically the only theoretical tool to dothis was the Monte-Carlo numerical simulation. It wasused in one of the pioneering studies of OMAR to showthe possibility of its bipolaron mechanism (the mecha-nism related to double-occupation of a single site withtwo electrons in the spin-singlet state) [13]. However,this method has its drawbacks. It is a numerical methodnot suited for analytical theory. Also, it is based on thesemi-classical nature of hopping transport were all thequantum mechanics is included in the electron hoppingrates. It has some problems with spin correlations thatactually have quantum nature. In [13] the spin was de-scribed semi-classically as the two possibilities for an elec-tron: to have spin up or spin down. This approximationcan readily be used in Monte-Carlo simulation, however,it cannot describe actual spin rotation around hyperfinefields. To describe it the spin should be allowed to bedirected along any axis, not only “up” and “down”. Itwill be shown that it requires a more rigorous descriptionof spin correlations based on quantum mechanics.There is an approach to consider pair correlations inclose pairs of sites as a modification of Miller-Abrahamsresistors [22–24]. However, up to now, electron spin wasconsidered with this approach only in the semi-classical“up” and “down” model.Very recently the approach that allows to include cor-relations of arbitrary order into the analytical theoryof hopping transport was developed [25, 26]. The ap-proach operates with correlation kinetic equations (CKE)that relate occupation numbers and their correlations. Itis based on Bogolubov chain of equations [27, 28]. In[25, 26] this approach is developed only for charge cor-relations, i.e. it does not consider electron spin in anymodel.The goal of this study is to develop CKE theory thatincludes spin correlations and explicitly takes into ac-count their quantum nature and to apply this theory tothe problem of organic magnetoresistance. The studyis restricted to the bipolaron mechanism of OMAR inlow electric fields and large Hubbard energy. The spinrelaxation is considered to be provided by hyperfine in-teraction with atomic nuclei.Although the qualitative theories can give a generalunderstanding of OMAR it is desirable to have an ap-proach that can be used to calculate OMAR explicitly.The progress in general understanding of organic semi-conductors and in simulation technics leads to the pos-sibility to calculate the microscopic properties of organicmaterials: the energy of molecular orbitals and theiroverlap integrals [29]. Some additional study of electron-phonon interaction in organic materials may lead to thepossibility of direct calculation of hopping rates. If allthese properties will be known the CKE approach can beused to quantitatively calculate the magnetoresistance.In the present paper, the magnetoresistance is calculatedin model systems. It is shown that the so-called lineshapeof magnetoresistance (the shape of the dependence of re-sistivity on the applied magnetic field) depends on theproperties of short-range electron transport. Differentmodel systems with identical hyperfine interaction butdifferent statistics of hopping rates show different line-shapes of organic magnetoresistance. The obtained line-shapes include (but are not limited to) the two empiricallaws that are most often used to fit experimental data[6, 7]: H / ( H + H ) and H / ( | H | + H ) . Here H isthe applied magnetic field and H is a fitting parame-ter. It gives hope that the calculations of OMAR can beused to relate microscopic models of hopping transportin organics with experimental results.The paper is organized as follows. In Sec. II I dis-cuss the model that is used to describe organic semi-conductors. In Sec. III the general mathematical defi-nitions of quantum spin and charge correlations are in-troduced. In Sec. IV the kinetic equations that relatethese correlations to currents are derived. In Sec. V theobtained system of CKE is used to describe the possi-ble lineshapes of the bipolaron mechanism of OMAR. InSec. V A it is done analytically in the model of modifiedMiller-Abrahams resistors. In Sec. V B it is done in themore general case with the numerical solution of CKE.In Sec. VI the general discussion of the obtained resultis provided. In Sec. VII the conclusion is given. Somepart of quantum mechanical calculations that are madeto derive spin CKE is discussed in the appendix A. II. MODEL
The following model is considered in the present work.The material contains a number of hopping sites wereelectrons (or polarons) are localized. In principle a hop-ping site can contain two polarons, however, the energy of its double occupation is larger than the energy of itssingle occupation. If the energy of the single occupa-tion of site i is ε i , the energy of its double occupation is ε i + U h , were U h is the Hubbard energy. The spins of elec-trons on a double-occupied site should form spin-singlet.The possibility of double occupation with electrons in atriplet state is neglected. It is considered that the sys-tem has some concentration of electrons and the Fermilevel µ . The current in small applied electric fields in thelinear response regime is discussed.I consider the Hubbard energy to be much larger thantemperature and energy differences in the hopping pro-cess. In this situation all the sites participating in hop-ping transport can be divided into the two groups. A -type sites have the energy of single-occupation near Fermienergy ε i ∼ µ . They have 0 or 1 electrons but are neverdouble-occupied because ε i + U h is too large for an A -typesite i . B -type sites have the energy of double occupationnear chemical potential ε j + U h ∼ µ . They always haveat least one electron because ε j ≪ µ for a B -type site j . Therefore “unoccupied” B -type site is a B -type sitewith one electron and has the spin degree of freedom. Anoccupied B -type site has two electrons in the spin-singletstate. This model is most convenient for the descriptionof a situation when the distribution of site energies ε i is broad not only compared to temperature but also tothe Hubbard energy U h (Fig. 2 (A)). In this case, onlya small part of hopping sites effectively participate intransport. The consideration of sites without an energylevel near µ significantly complicates the numeric simu-lation and has little impact on the result. I adopt thesimplified model where there are two independent densi-ties of states for A -type sites and for B -type sites (Fig.2 (B)). (cid:1) g (cid:2)(cid:1)(cid:3)(cid:4)(cid:4) -U h (cid:1) g (cid:2)(cid:1)(cid:3)(cid:4)(cid:4) -U h ype (B)(A) -type-type FIG. 1: The density of states in the system with a broad dis-tribution of energies (A) and in the adopted simplified model(B). The A -type and B -type sites marked with red and bluecolor correspondingly. In the following part of the text, the energy of doubleoccupation of B -type site i is denoted as ε i because theenergy of its single occupation does not appear in thetheory.The hopping between sites is controlled by hoppingrates W ij . It is assumed that the electron spin is alwaysconserved in the hopping process. When sites i and j have the same type, the site j is occupied and site i isnot, the electron hops from j to i with rate W ij . W ij = W | t ij | exp (cid:18) − max[( ε i − ε j ) , T (cid:19) . (1)Here t ij is the overlap integral between localized stateson sites i and j . W describes the strength of electron-phonon interaction. In principle, it can have a power-law dependence of site energies ε i , ε j . However, thisdependence is not universal (is material-depend) and isnot considered in the present study. Therefore W istreated as a constant. When site j is of type B and i is A -type site ( B → A hop) the hopping occurs with therate 2 W ij because both of the two electrons on site j canhop to i . After the hop, the spins of electrons are in thesinglet state. In the situation of A → B hop, when spinsare in thermal equilibrium, the hopping occurs with therate W ij /
2. However, this rate increases to 2 W ij whenthe spins are in singlet state. It ensures that the detailedbalance holds in the thermal equilibrium in a pair of sitesof different kinds. The A → B hop is impossible whenthe spins are in a triplet state because the electron spinis conserved during the hop.The equilibrium occupation number of an A -type site i is equal to n (0) i = 1 / (1 + e ( ε i − µ ) /T / A -type site [2]. The A -type site i has twooccupied states with different spin. The joint probabilityof the occupied states is 2 e − ( ε i − µ ) /T /Z i where Z i is thestatistical sum of site i : Z i = 1 + 2 e − ( ε i − µ ) /T . It leads tothe mentioned expression for the occupation probability n (0) i . For the B -type site j the free state is degenerate.With similar arguments it leads to slightly different ex-pression for the equilibrium occupation number of site j : n (0) j = 1 / (1 + 2 e ( ε j − µ ) /T ).In any pair of sites i − j , the equality holds withoutrespect for site typesΓ ij = W ij p ( ij ) sp (1 − n (0) i ) n (0) j = W ji p ( ji ) sp (1 − n (0) j ) n (0) i = Γ ji . (2)Γ ij is the number of electrons that hops from site j tosite i in unit of time in thermal equilibrium. p ( ij ) sp is thespin term in the hopping probability. It is equal to 1 / j has type A and site i has type B . When site j is B -type site and i is A -type site, p ( ij ) sp = 2. Whenthe types of sites i and j are the same, p ( ij ) sp = 1. Theequation (2) shows that detailed balance holds in thermalequilibrium.The organic magnetoresistance is closely related to thedynamics of spin correlations. The mathematical de-scription of these correlations is introduced in Sec. III.Here I describe the physics that is considered. The mainreason for spin dynamics is the hyperfine interaction withatomic nuclei. It is described as effective on-site magneticfields H ( i )hf that are different on different sites. Thesefields should be added to the external magnetic field H . FIG. 2: Rotation of spins on sites i and j with Larmor fre-quencies Ω i and Ω j . The electron spin rotates around the total field with Lar-mor frequency Ω i = µ b g ( H + H ( i )hf ). This descriptionis valid when the slow dynamics of nuclear spins can beneglected.The rotation in hyperfine fields can modify spin cor-relations. Consider that at some point of time the spinsof electrons on sites i and j were parallel. After sometime due to rotation with different vector frequencies Ω i and Ω j there will be an angle between spin directions.In combination with electron hops, this rotation leads tospin relaxation [30–32].I also introduce a phenomenological time of spin re-laxation τ s . It can describe the time of on-site spin re-laxation related to (for example) the spin-phonon inter-action. However, the main reason to introduce τ s is thepossibility to compare the results with previous theorieswhere spin relaxation was treated in this way [13, 22]. III. SPIN CORRELATIONS
The site occupation numbers are insufficient for the de-scription of OMAR that is controlled by spin degrees offreedom. OMAR appears due to spin correlations. Start-ing from the first qualitative studies of OMAR [12, 13]the magnetoresistance is attributed to different proba-bilities of parallel and anti-parallel configurations of twospins. In terms of statistics the description of these prob-abilities is possible only beyond the mean-field approxi-mation. By definition the mean-field approximation dealswith occupation numbers on a single site. It can describeonly the situation when one spin direction on some siteis more probable than another without respect to othersites. It corresponds to the appearance of averaged on-site spin polarization that is impossible at room temper-atures and magnetic fields ∼ gs . OMAR is controlledby the probabilities of the relative directions of two spins.For example let us consider the situation when paral-lel configuration is more probable than anti-parallel. Inthis case the probability for both spins to have up direc-tion ↑↑ is equal to the probability of ↓↓ direction but islarger than the probability of ↑↓ and ↓↑ directions. Interms of statistics this situation is described with cor-relations. When the correlations are neglected in themean-field approximation, it is impossible to distinguishbetween parallel configuration of two spins (without aver-aged polarization) and the equilibrium statistics. There-fore a microscopic theory of OMAR should be developedbeyond the mean-field approximation and take at leastsome correlations into account. The approach adoptedin the present paper is to write the equations for arbi-trary correlations and neglect the insignificant ones atthe last step of the theory when the equations are solvednumerically.In this section, the general notations for spin correla-tions between electrons on different sites are introduced.The special attention is paid to the quantum nature ofspin. It allows to describe simultaneously the hoppingtransport and the spin rotation around local hyperfinefields. These fields are responsible for spin relaxation indifferent materials [30–32] including many organic semi-conductors.At first, I describe the electron state on a single A-typesite i . The electron has the two quantum-mechanicalstates: with spin up | ↑i or down | ↓i . The generaldescription of its state can be given by 2 × b ρ i . Its diagonal terms are real. Their sum is unitybecause the site i is considered to be occupied. The non-diagonal terms are complex and are conjugate. It leadsto three independent parameters describing the densitymatrix. These parameters can be selected to have clearphysical meaning [33]: the averaged spin polarizations ofsite in the directions x , y and z . The density matrix b ρ i can be re-constructed with this averaged polarizations. b ρ i = (1 / (cid:16)b s xi b σ x + s yi b σ y + s zi b σ z (cid:17) . Here s x,y,zi are theaveraged values of operators ˆ s x,y,zi of spin polarizationalong the axes x , y and z . b σ x,y,z are the Pauli matrices. b × b σ Now consider an A -type site i with finite occupationprobability n i ≡ s i . When the conductivity is due toelectron hops the terms of density matrix with uncertainoccupation numbers can be neglected between hops (al-though they should be treated with perturbation theorywhen the hopping rates are calculated). The density ma-trix b ρ i can be given in block-diagonal form. The free siteis described by the block b ρ (0) i that actually is a singlenumber equal to 1 − n i . The single-occupied site is de-scribed by the block b ρ (1) i that can be expressed as follows b ρ (1) i = P p s pi σ p (cid:0) s i b σ + s xi b σ x + s yi b σ y + s zi b σ z (cid:1) . (3)Here index p can have one of the four values (0 , x, y, z ).The averaged quantities s i , s xi , s yi and s zi can describeany state of the site i .It is possible to describe spin correlations in a similarmanner. The density matrix of two sites i and j can beexpressed as the four blocks with well-defined occupation numbers. Let i be an A -type site and j be a B -type site.In this case the site i can have 0 or 1 electrons and j canhave 1 or 2 ones. The joint density matrix of sites i and j can be expressed in terms of four blocks b ρ ij = b ρ (12) ij b ρ (11) ij b ρ (02) ij
00 0 0 b ρ (01) ij . (4)Here the upper indexes stand for the occupation numbersof the sites. For example, b ρ (12) ij describes the part of den-sity matrix related to single-occupied site i and double-occupied site j . b ρ (11) ij is 4 × b ρ (11) ij = 14 X pq s pi s qj (cid:16)b σ ( i ) p ⊗ b σ ( j ) q (cid:17) (5)where s pi s qj is the quantum mechanical average of opera-tor ˆ s pi ˆ s qj . Here ˆ s x,y,zi and ˆ s x,y,zj are the operators of spinpolarizations of sites i and j correspondingly. s i and s j are the operators of single occupations of sites i and j . They always have well-defined values between hops.Note that the relation between occupation number and s is different for different types of sites. For A-type site i , s i = n i . For a B-type site j , s j = 1 − n j because B -type sites have one electron in unoccupied state. Ikeep these double notations because the notation n i isuseful to track the charge conservation law while the no-tation ˆ s allows to give the expression for density matrixin terms of correlations in unified form for both A-typeand B-type sites.The sign ⊗ in Eq. (5) denotes the Cartesian product ofmatrices. The upper indexes ( i ) and ( j ) in b σ ( i ) p and b σ ( j ) q have no mathematical meaning but help to track whatPauli matrix is related to what site.The block b ρ (02) ij is a single number equal to (1 − n i ) n j The blocks b ρ (12) ij and b ρ (01) ij are 2 × b ρ (12) ij = X p (1 − s j ) s pi b σ ( i ) p , b ρ (01) ij = X p (1 − s i ) s pj b σ ( j ) p . (6)The whole matrix ρ ij can be parameterized with 24 av-eraged values: s pi , s pj and s pi s qj . These values have clearphysical meaning: they describe the occupation proba-bilities, mean spin polarizations and their correlations.The kinetic equations for these averaged values would al-low to describe all the dynamics of the density matrixfor a system with hopping transport. This result can begeneralized for arbitrary number of sites. Let I be someset of sites. The density matrix b ρ I of this set can be de-scribed by averaged products of ˆ s pi in all the subsets I n of the set I . s PI n = Y i ∈ I n ˆ s p i i , P = { p , ...p n ( I n ) } , I n ⊂ I. (7)Here P is the set of upper indexes p i equal to 0, x , y or z related to the sites i in the set I n .All the kinetics of a hopping system can be describedwith s PI . However, these values cannot be considered ascorrelations. When two A-type sites i and j are not corre-lated s i s j = n i n j = 0 in thermal equilibrium. Therefore s PI cannot be neglected even when correlations inside theset I are not important.The idea of correlation kinetic approach [25] is to writeequations for correlations themselves and neglect the cor-relations at large distance and of high order. To followthis idea the correlations c PI are introduced c PI = Y i ∈ I c p i i , c i = n i − n (0) i , c αi = ˆ s αi . (8)Here the Greek letter α stands for the spin projectionon x , y or z . n (0) i is the equilibrium filling number. Inthis notations c αi = 0 when we neglect averaged spin po-larization and c i = n i − n (0) i is the perturbation of theoccupation number due to applied electric field. Whenthe hopping system is close to equilibrium, eq. (8) de-scribes the correlations that can be neglected if occupa-tion numbers and spins in the set I are not considered tobe correlated. IV. KINETIC EQUATIONS FOR SPIN ANDCHARGE CORRELATIONS
In this section the kinetic equations are derived for thecorrelations defined in Sec. III. Consider the correlation c PI in some set I of sites. It can be changed due to one ofthe following processes: the hopping of electrons betweensites of the set and outer sites, the hopping of electronsinside set I and due to the internal spin dynamics. ddt c PI = X i ∈ I,k / ∈ I (cid:18) ddt (cid:19) ik c PI + X i,j ∈ I (cid:18) ddt (cid:19) ij c PI + X i ∈ I (cid:18) ddt (cid:19) i c PI (9)Here in r.h.s of Eq. (9) the symbolical expressions fordifferent terms are used. ( d/dt ) ij c PI means the changingrate of c PI due to the hops between sites i and j of the set.( d/dt ) ik c PI stands for the transitions between site i fromthe set and outer site k . ( d/dt ) i c PI describes the changingrate of c PI due to the internal spin dynamics (spin rotationand phenomenological spin relaxation) related to site i .The term ( d/dt ) ik describes the transition of correla-tions between sets of sites. It can be expressed as follows. (cid:18) ddt (cid:19) ik c pi c P ′ I ′ = T ( p ) ik c pk c P ′ I ′ − T ( p ) ki c pi c P ′ I ′ . (10) Here T ( p ) ik is the rate of charge or spin transition fromsite k to site i . c P ′ I ′ describes the part of correlation thatis not related to site i and is conserved during i ↔ k hops. I ′ = I \{ i } where the notation \ stands for the setdifference. P ′ is the set of indexes from P other than theindex p related to the site i . T (0) ik = W ik p ( ik ) sp (cid:16) − n (0) i (cid:17) + W ki p ( ki ) sp n (0) i (11) T (0) ik describes the rate of transition of small perturbationof charge density from site k to site i [25]. FIG. 3: The spin transfer from site i to site j in differentpairs of sites. On the left hand side of the figure the site i isspin-polarized and site j is not. On the right-hand side of thefigure (after the hop) the site j acquires spin polarization inthe direction of the initial polarization of site i . T ( α ) ik describes the process of spin polarization trans-fer between sites. This process is different for differenttypes of sites, as shown in Fig. 3. In a pair of A-typesites the spin transfer is achieved due to hops of spin-polarized electrons. In a pair of B-type sites the hops ofspin-polarized holes are responsible for the spin transfer.In a mixed AB pair the spin transfer occurs because elec-tron with one spin projection can hop from the A -typesite to the B -type site while the electron with other spinprojection cannot. It leads to the following expressionsfor spin transfer rates in different pairs of sites (see [21]for details). T ( α ) ik = W ik (1 − n (0) i ) , AAW ki n (0) i , BBW ki n (0) i / , ABW ik (1 − n (0) i ) / , BA (12)The term ( d/dt ) ij for charge correlations is derived in[25] using the fact that joint occupation of sites i and j cannot be changed due to i ↔ j hops, ( d/dt ) ij n i n j = 0.This argument should be generalized to include spin anddifferent types of sites. When both the sites i and j have the same type the hops between them are impos-sible when both of them are single occupied. It leadsto the expression ( d/dt ) ij s pi s qj s P ′ I ′ = 0 for AA and BB pairs of sites i , j . When the sites i and j have differ-ent types the hop between them is possible when bothare single-occupied. Actually it is the process that re-lates spin correlations and charge transport in bipolaronmechanism of OMAR . Rate equations for this processwhen i is A -type site and j is B -type site are derived inappendix A with quantum mechanic approach: (cid:18) ddt (cid:19) ij s i s j s P ′ I ′ = 2 W ij n j ↑ n j ↓ (1 − n i ↑ )(1 − n i ↓ ) s P ′ I ′ − W ji s i s j − X α s αi s αj ) s P ′ I ′ , (13) (cid:18) ddt (cid:19) ij s i s αj s P ′ I ′ = W ji s αi s j − s i s αj ) s P ′ I ′ , (14) (cid:18) ddt (cid:19) ij s αi s j s P ′ I ′ = W ji s i s αj − s αi s j ) s P ′ I ′ , (15) (cid:18) ddt (cid:19) ij s αi s βj s P ′ I ′ = − W ij δ αβ n j ↑ n j ↓ (1 − n i ↑ )(1 − n i ↓ ) s P ′ I ′ + W ji δ αβ ( s i s j − X α s γi s γj ) s P ′ I ′ − W ji s αi s βj − s βi s αj ) s P ′ I ′ (16)Here Greek indexes α , β and γ stand for spin polariza-tions along Cartesian axes.The physical meaning of Eqs. (13-16) is as follows.Eq. (13) and the first two terms in r.h.s. of Eq. (16) showthat the hop from site i to site j is possible only whenboth sites are single-occupied and electron spins on thesesites are in singlet state. It decreases the probability ofsingle occupation of both sites. The backward hop from j to i leads to the single occupation of sites with electronsin singlet state. The third term in Eq. (16) leads to re-laxation of the antisymmetric combinations s αi s βj − s βi s αj .It can be shown that these combinations are related toa coherent combination of singlet and triplet states ofthe two electrons. When the electrons are either in thesinglet or in a triplet state s αi s βj − s βi s αj = 0. However,when their state is a coherent combination of singlet andtriplet s αi s βj − s βi s αj = 0. Even if due to some reason theprobabilities of singlet and triplet states are conserved,their coherent combinations relax because the hopping i ↔ j can occur in the singlet state and cannot occur ina triplet state. Similar term appears in spin dynamic ofa double quantum dot [34]. Eqs. (14,15) show that spintransfer process in AB pairs of sites is correlated withoccupation numbers.The known relations between c PI and s PI allow to ob-tain the expressions for ( d/dt ) ij c PI similar to Eqs. (13-16). These expressions are not provided here becausethey are quite cumbersome but can be given in a muchmore compact form when the correlation potentials areintroduced.The term ( d/dt ) i c PI is responsible for the internal spindynamics. It is present only when the index p corre-sponding to site i is a spin index. (cid:18) ddt (cid:19) i c α,P ′ i,I ′ = ǫ αβγ Ω i,β c γ,P ′ i,I ′ − τ s c α,P ′ i,I ′ (17)Here Ω i,β is the projection of spin rotation frequency vec-tor on site i to the axis β . τ s is a phenomenological on-sitespin relaxation time. ǫ αβγ is Levi-Civita symbol.In the linear response regime it is useful to introduceeffective correlation potentials ϕ PI . ϕ i = c i / (1 − n (0) i ) n (0) i , ϕ αi = c αi / ( s i ) eq . (18)Here ( s i ) eq is the equilibrium probability of single occu-pation of site i . It is equal to n (0) i if site i has type A andto 1 − n (0) i if site i has type B . There is no averaging inthe definition (18). ϕ pi should be ensemble averaged insome combination to have the meaning of potential. Forexample ϕ xyijk = ϕ i ϕ xj ϕ yk can be considered as a potentialof correlation c xyijk .In these notations ( d/dt ) ik c p,P ′ i,I ′ corresponds to a “cor-relation flow” J p ; P ′ ik ; I ′ between correlations c p,P ′ i,I ′ and c p,P ′ k,I ′ : (cid:18) ddt (cid:19) ik c p,P ′ i,I ′ = − (cid:18) ddt (cid:19) ik c p,P ′ k,I ′ = J p ; P ′ ik ; I ′ . (19)When I ′ is the empty set J ik is the particle flow from site k to site i . The flow J p ; P ′ ik ; I ′ is expressed as follows J p ; P ′ ik ; I ′ = Γ ik Θ P ′ I ′ (cid:16) ϕ p,P ′ k,I ′ − ϕ p,P ′ i,I ′ + S pik,I ′ +Φ qq ′ p ( ik ) ϕ qq ′ ,P ′ ik,I ′ (cid:17) . (20)Here Γ ik is the average number of electrons that hopsfrom site k to site i in unit time in equilibrium (2). Θ P ′ I ′ isthe coefficient related to sites of I that do not participatein transition i ↔ k .Θ p j ,p l ,p m ,...j,l,m,... = θ p j j · θ p l l · θ p m m · ... (21) θ j = n (0) j (cid:16) − n (0) j (cid:17) , θ αj = (cid:0) s j (cid:1) eq . (22) S pik,I ′ is the source term related to the external electricalfield E . It is not equal to zero only when I ′ is empty setand p = 0. In this case S ik, ∅ = e Er ik where e is electroncharge, r ik is vector of distance between i and k .The term Φ qq ′ p ( ik ) ϕ qq ′ ,P ′ ik,I ′ describes the effect of higher-order correlations to the flow of lower-order correlations.For different indexes p , q and q ′ , the coefficients Φ qq ′ p ( ik )are Φ ( ik ) = n (0) k − n (0) i , Φ αα ( ik ) = τ k − τ i , Φ α α ( ik ) = (1 − τ i ) n (0) k − τ i (1 − n (0) k ) , Φ αα ( ik ) = τ k (1 − n (0) i ) − (1 − τ k ) n (0) i . (23)Here τ i = 0 for A -type site i and τ i = 1 if site i has type B . All the coefficients Φ qq ′ p not listed in (23) are equal tozero. For example when all the three indexes p , q and q ′ are spin indexes Φ qq ′ p = Φ βγα = 0.The term ( d/dt ) ij c pq,P ′ ij,I ′ is closely related to the corre-lation flow between correlations c q ′ ,P ′ i,I ′ and c q ′ ,P ′ k,I ′ : (cid:18) ddt (cid:19) ij c pq,P ′ ij,I ′ = − Φ pqq ′ ( ij ) J q ′ ; P ′ ij ; I ′ − Υ pqij Θ pq,P ′ ij,I ′ (cid:16) ϕ pq,P ′ ij,I ′ − ϕ qp,P ′ ij,I ′ (cid:17) (24)Note that the same coefficients Φ qq ′ p enter the equations(20) and (24). The second term in r.h.s. of Eq. (24) isrelated to the relaxation of coherent singlet-triplet com-binations. It contains the coefficient Υ pqij that is equal tounity when sites i and j have different types and p and q are spin indexes. Otherwise, Υ pqij = 0.The term ( d/dt ) i c α,P ′ i,I ′ related to the spin rotation andrelaxation should also be expressed in terms of potentials (cid:18) ddt (cid:19) i c α,P ′ i,I ′ = Θ α,P ′ i,I ′ ǫ αβγ Ω i,β ϕ γ,P ′ i,I ′ − ϕ α,P ′ i,I ′ τ s ! . (25)In a stationary system the derivatives dc PI /dt are equalto zero. Therefore the equations (9), (19),(20), (24) and(25) compose a closed system of linear equations for thecorrelations potentials. It incudes all the charge and spincorrelations. The total number of these correlations is ex-tremely large, 4 N where N is the number of sites. How-ever, one can hope that correlations between sites at verylarge distances and the correlations of very high orderare not relevant for the electron transport and can beneglected. Actually, to treat reasonably large systemssome of the correlations should be neglected to make thesystem of equations solvable. The idea of CKE approachis to write the equations in general form relevant for arbi-trary correlations and make the cutoff at the “final step”taking into account the structure of considered system(that defines the correlations that are really relevant)and the possibility to numerically solve the system ofequations of the desired size. When some correlations are neglected in this way, the potentials ϕ PI of these cor-relations are considered to be equal to zero in all theequations. V. MAGNETORESISTANCE DUE TO THERELAXATION OF SPIN CORRELATIONS
In this section, the discussed approach to the theory ofhopping transport with spin correlations is applied to thebipolaron mechanism of OMAR. As it was discussed inthe previous section, it is necessary to cut the system ofkinetic equations at some point. The most simple cutoffis the model when only the correlations in close pairs ofsites are considered. In this case, it is possible to reducethe problem to a network of modified Miller-Abrahamsresistors. This approach was used in [22, 24] with thesemi-classical model of spins up and down. However, inthe present study, the quantum nature of spin correla-tions is taken into account and the expressions for resis-tors are different from [22]. This model is discussed inSec. V A. It allows the analytical solution in the limitingcases of fast and slow hopping.The long-range and high-order correlations can betaken into account with the numerical solution of kineticequations. Such solutions are provided in Sec. V B andare compared with analytical results.
A. Modified resistor model
When the long-range spin correlations are neglectedthe rate equation for spin correlations c αβij includes onlythe spin generation due to the electron flow J ij and theinternal spin dynamics. The effect of other sites is re-duced to spin relaxation. When the correlation i − j isconsidered in this model the spins on other sites are as-sumed to be in thermal equilibrium. The transition ofcorrelation i − j to other sites can be formally includedinto internal spin dynamics as additional source of relax-ation. ddt c αβij = − R ( ij ) αβ ; α ′ β ′ c α ′ β ′ ij + δ αβ ( τ i − τ j ) J ij . (26)Here R ( ij ) αβ ; α ′ β ′ is a matrix that describes the dynamicsand relaxation of correlations. R ( ij ) αβ ; α ′ β ′ = γ ij δ αα ′ δ ββ ′ − ǫ αγα ′ δ ββ ′ Ω i,γ − ǫ βγβ ′ δ αα ′ Ω j,γ + Γ ij ( s i ) eq ( s j ) eq ( δ αα ′ δ ββ ′ − δ αβ ′ δ βα ′ ) . (27) γ ij is the effective rate of relaxation of spin correlationseither due to phenomenological on-site spin relaxationmechanism or due to electron transition to or from othersites. γ ij = 1 τ s + X k = i,j (cid:16) T ( α ) ki + T ( α ) kj (cid:17) . (28)Note that spin transition rates T αki do not depend onvalue of spin index α = x, y, z . The index is kept onlyto show that it is a spin index, not the charge index0. It is assumed in the present section that i − j is theresistor that controls the resistivity of some mesoscopicpart of the sample. It happens when the rate of hoppinginside the pair i − j is slow compared to the hoppingbetween this pair and other sites. It allows to assumethat Γ ij ≪ γ ij . Therefore Γ ij and the last term in r.h.s.part of Eq. (27) is neglected in the present section.It is possible to give a closed expression for J ij with ac-count to short-range pair correlations in terms of inversematrix ( R ( ij ) ) − . J ij = ϕ j − ϕ i Γ − ij + F s ( ij ) + F c ( ij ) , (29) F s ( ij ) = δ αβ δ α ′ β ′ (cid:0) R ( ij ) (cid:1) − αβ ; α ′ β ′ ( τ i − τ j ) ( s i ) eq ( s j ) eq , (30) F c ( ij ) = (cid:16) n (0) i − n (0) j (cid:17) (cid:16)P k = i,j T (0) ki + T (0) kj (cid:17) − n (0) i n (0) j (1 − n (0) i )(1 − n (0) j ) . (31)Here Γ − ij corresponds to ordinary Miller-Abrahams resis-tance. F s and F c describe the additional resistance thatappear due to spin and charge correlations correspond-ingly.The effect of the external magnetic field on the con-ductivity is incorporated in F s . As shown in eqs. (27,30)it depends on the vectors of on-site rotation frequen-cies Ω i and Ω j that are proportional to the sum of hy-perfine field and applied external magnetic field Ω i = µ b g ( H + H ( i )hf ) / ~ . However, even if the system is com-posed from only one resistor the averaging over Ω i and Ω j is required. The hyperfine fields are slowly changeddue to the nuclear spin dynamics. Although this dynam-ics is considered to be slow compared to electron hopsand electron spin rotation, it is usually fast compared tothe current measurement procedures. Therefore the fi-nal expression for dc resistivity should be averaged overhyperfine fields.I assume that hyperfine fields have normal distribution p ( H ( i )hf ,α ) = 1 √ πH hf exp − (cid:16) H ( i )hf ,α (cid:17) H . (32)Here H hf is the typical value of hyperfine fields. Thedifferent components of the hyperfine field on a givensite and the fields on different sites are considered not tobe correlated.The spin correlation part of resistance is related to thereverse relaxation function R ( γ ij , H, H hf ): F s = ( τ i − τ j ) R / ( s i ) eq ( s j ) eq . R = D δ αβ δ α ′ β ′ ( R ( ij ) ) − αβ ; α ′ β ′ E hf . (33) Here h ... i hf means the averaging over the hyperfine fields. R can be thought of as the time of relaxation of proba-bilities for the two spins to be in singlet or triplet state.In the general case R ( γ ij , H, H hf ) can be found numeri-cally. However, it is possible to find it analytically in thelimiting cases of slow and fast hops.In the limit of fast hops the rate of relaxation γ ij due to electron transition to other sites is fast com-pared to the typical rate of rotation in hyperfine fields h hf = µ b gH hf / ~ . The rate of rotation in the externalmagnetic field h ext = µ b gH/ ~ is arbitrary. In this caserelaxation matrix R ( ij ) can be divided into R ( ij )1 relatedto hopping and rotation in the external magnetic field R ( ij )1 = − γ ij δ αα ′ δ ββ ′ + µ b gH ( ǫ αzα ′ δ ββ ′ + ǫ βzβ ′ δ αα ′ ) / ~ and R ( ij )2 related to the rotation in hyperfine fields. The firstof this matrices R ( ij )1 can be inverted analytically. Thesecond one can be considered as a small perturbation.The total inverse relaxation matrix averaged over hyper-fine fields can be approximately expressed as (cid:28)(cid:16) R ( ij ) (cid:17) − (cid:29) hf = (cid:16) R ( ij )1 (cid:17) − + (cid:28)(cid:16) R ( ij )1 (cid:17) − R ( ij )2 (cid:16) R ( ij )1 (cid:17) − R ( ij )2 (cid:16) R ( ij )1 (cid:17) − (cid:29) hf . (34)In principle the expression for ( R ( ij ) ) − also includes thefirst order term ( R ( ij )1 ) − R ( ij )2 ( R ( ij )1 ) − , however, it be-comes equal to zero after the averaging over hyperfinefields.With straightforward calculations, Eq. (34) leads toexplicit expression for the function RR = 1 γ ij " − h γ ij h /γ ij ! . (35)The dependence of resistance on the magnetic field cor-responding to eq. (35) is described by Lorentz function.Note that its width is controlled not by the relation of ex-ternal magnetic and hyperfine fields h ext /h hf but by therelation of hopping rate ant the rate of rotation in the ex-ternal field h ext /γ ij . The magnetoresistance is relativelyweak due to the small prefactor h /γ ij .The opposite limit is the situation of slow hopping, γ ij ≪ h hf . In this case it is possible to use the ran-dom phase approximation. It is assumed that componentof electron spin on site i normal to the on-site effectivemagnetic field H + H ( i )hf relaxes very fast due to the ro-tation around this field. However, the component along H + H ( i )hf is conserved until the electron is transferred tosome other site. In this case R is proportional to the av-eraged squared cosinus of the angle between on-site fields R = D cos (cid:16) H ( i ) , H ( j ) (cid:17)E hf /γ ij . (36)Here H ( i ) = H + H ( i )hf . The averaging in (36) can be done FIG. 4: Reverse relaxation function R ( h ext /h hf ) corresponding to different relations γ ij /h hf . Solid blue curve in all the figuresis numeric calculation of R . Yellow dash-dot curve in Fig. (A) and (B) is approximation with Eq. (37). Red dashed curve inFig. (B) is the non-Lorentz fit described in text. Red dashed curve on Fig.(C)-(D) is the Lorentzian fit. Yellow dash-dot curvein Fig. (D)-(F) is approximation with Eq. (35). Inset in Fig. (C) shows the curves at small magnetic fields. analytically in terms of special functions R = 13 γ ij + 12 γ ij H × " √ H ( H − H ) + H D + (cid:18) H √ H hf (cid:19) . (37)Here D + ( x ) is the Dawson function D + ( x ) = e − x R x e t dt . Eq. (37) shows that in slow hopping limitthe dependance of resistance on magnetic field has non-Lorentz shape. It saturates when applied magnetic fieldis much larger than hyperfine fields, while γ ij controls itsoverall strength.In Fig. 4 I compare the approximate expressions (35)and (37) with the function R calculated numerically.It can be seen that Eq. (35) can quantitatively de-scribe R only for quite large values of hopping rate γ ij & h hf . However, the dependence of R on the ap-plied magnetic field can be described by the Lorentzian R = A + B/ ( h + e γ ) for significantly smaller γ ij & h hf .Here A , B and e γ are fitting parameters. At smaller hop-ping rates γ ij . . h hf the function R becomes non-Loretzian. It is most clearly seen when comparing thenumeric results for R with its Lorentzian fit at smallmagnetic fields as shown on the inset in Fig. 4(C).For small hopping rates γ ij . . h hf the non-Lorentzian fit related to Eq. (37) becomes relevant. Itis the fit R = A + B e R ( H, H hf , γ ij ) where e R is describedwith Eq. (37). For γ ij . . h hf Eq. (37) can describethe reverse relaxation function quantitatively.It is interesting to compare the results of the presentsection with results of [22] where the correlations inclose pairs of sites were considered with semi-classicalmodel of spin. In [22] OMAR was related to spin re- laxation time τ s ( H ) that had a phenomenological de-pendence on the applied magnetic field. The approachof the present section can describe quantum spin corre-lations with the same relaxation time. When the spinrelaxation is reduced to a single time τ s ( H ) the reverserelaxation function can be found analytically and the spinpart of resistance is equal to F s = 3( τ i − τ j ) τ s ( H ) / [(1 + γ ij τ s ( H ))( s i ) eq ( s j ) eq ]. It is exactly 3 times larger thanthe correction to Miller-Abrahams resistor due to spincorrelations obtained in [22]. Note that F c quantitativelyagree with the correction to Miller-Abrahams resistor dueto charge correlations obtained in [22]. The difference be-tween the results is related to the quantum nature of spincorrelations taken into account in the present paper. Itappears that semi-classical model of spin cannot quan-titatively describe OMAR even if the spin relaxation isreduced to a single time τ s ( H ). B. Numerical results
This section includes the results of the numerical so-lution of correlation kinetic equations for several disor-dered systems. The calculation is made with some ofthe long-range and high-order correlations taken into ac-count. The results are compared with the model de-scribed in Sec. V A to show when the theory that includesonly the correlations in close pairs of sites is applicableand when it is not.I start from a single numerical sample consisting of25 sites. The sample is shown in Fig. 5(A). The sites areplaced on a square lattice with the same overlap integrals t ij between neighbor sites. The site energies are selectedindependently with normal distribution with the stan-dard deviation ∆ E . The temperature is T = ∆ E . Some0 FIG. 5: The dependence of conductivity on the applied mag-netic field in mesoscopic numeric sample. (A) the structureof the sample. The color bar stands for site energies. (B)conductivity calculated with different approximations as de-scribed in the text. sites are randomly selected to have type B. They aremarked with crosses in Fig. 5(A). Other sites have typeA. Each site is ascribed with a random hyperfine field.The typical rate of rotation in random fields is equal tothe pre-exponential term in the hopping rate betweenneighbors h hf = w . Here I define the pre-exponentialterm in the hopping rates as w = W | t ij | . The peri-odic boundary conditions are applied.In Fig. 5(B) the calculated conductivity of this sampleis presented. The conductivity is calculated in four differ-ent approximations. The blue dashed curve correspondsto (1 ,
2) approximation where pair correlations betweensites i and j are included in the theory when the distancebetween sites i and j along the lattice bonds is 1. It is theapproximation used in Sec. V A. In general, the notation( p, q ) stands for the approximation when correlations upto the order q are considered provided that the distancebetween sites in the correlations along lattice bonds isno longer than p . The yellow dash-dot curve correspondsto (3 ,
2) approximation when most pair correlations aretaken into account. Green solid curve corresponds tojoint (3 ,
2) and (2 ,
4) approximation when most of paircorrelations and correlations up to the 4-th order in closecomplexes of sites are considered. The red dots stand forjoint (3 ,
3) and (2 ,
5) approximation where most of thecorrelations of the third order and the correlations up to5-th order in close complexes are included in the theory.The introduction of new correlations into the the- ory decreases the calculated conductivity and increasesits calculated dependence on the applied magnetic field.However, the difference between (3 , ,
4) and (3 , ,
5) approximations is rather small and one can hopethat (3 ,
2) + (2 ,
4) approximation adequately describesthe system. In the following analysis of larger numericalsamples, the correlations of order q > p ≤
2. Thepair correlations will be considered at slightly larger dis-tances. Unfortunately, the number of spin correlationsgrows extremely fast with the correlation order and itwas technically impossible to go beyond the (3 ,
3) + (2 , × σ ( H )for real samples is symmetric even if they are mesoscopic.In Fig. 6 I show results for larger 10 ×
10 samples. Theresults are averaged over 30 disorder configurations. Theproperties of the numerical samples (except their size) arethe same as in Fig. 5. Fig. 6(A) shows the structure ofone of these samples. Fig. 6(B) and (C) show the magne-toresistance ( R ( H ) − R (0)) /R (0) in different scales. Here R ( H ) is the sample resistance in the external magneticfield H . The blue points stand for the magnetoresistancecalculated with the numeric solution of CKE. The paircorrelations at the distance no longer than 4 and 4-thorder correlations with the distance 2 were taken intoaccount ((4 ,
2) + (2 ,
4) approximation). Green dashedcurve is the fit with Lorentzian. Yellow solid curve is thefit with function R described in Sec. V A. It means thatthe expression R ( H ) − R (0) R (0) = A × R ( e γ, H, H hf ) − R ( e γ, , H hf ) R ( e γ, , H hf ) (38)was used for fitting. Here A and e γ are the fitting param-eters. e γ can be considered as the effective rate for a cor-relation to leave the pair of sites where it appeared. Theresults shown in Fig. 6(B) correspond to e γ ≈ . h hf . A isthe general amplitude of magnetoresistance, it describesthe relative part of sample resistance that is related tospin correlations.The shape of resistance dependence on the magneticfield significantly deviates from Lorentzian. It can beeasily seen in Fig. 6(C) where the results for small mag-netic fields are shown on a close scale. The fitting (38)has better agreement with numerical simulation.In Fig. 7 the similar analysis is provided for similarnumeric samples with different rates of hopping betweenneighbors w = 0 . h hf and w = 5 h hf . Other charac-teristics of the samples are the same as in the previous1 (A) (B) (C) FIG. 6: Results of CKE solution for random 10 ×
10 numerical samples. (A) the structure of a single sample. (B) and (C)calculated magnetoresistance in different scales averaged over 10 random samples.FIG. 7: Magnetoresistance in numerical samples with slow w = 0 . h hf and fast w = 5 h hf hopping rates. Blue dotson all subplots correspond to results of numerical solutionof CKE. Yellow solid curves correspond to Eq. (38). Greendashed curves in Figs (A) and (B) is the fit with Eq. (38)where function R is described by Eq. (37) and fitting hy-perfine fields are artificially increased h hf → . h hf . Greendashed curve in Figs (C) and (D) is the fit with Lorentzian. numeric experiment including the averaging over 30 dis-order configurations. Fig. 7 (A) and (B) shows the resultsfor the samples were the hopping is slow, w = 0 . h hf .Blue points correspond to the numeric solution of CKE.Yellow solid curve corresponds to the fit with Eq. (38),the fitting parameter e γ was e γ = 1 . h hf . Green dashedcurve corresponds to fit with Eq. (38), where function R is described by Eq. (37) (that is valid in the limit e γ ≪ h hf ) but the strength of hyperfine fields was artifi-cially increased h hf → . h hf to achieve better fit withnumerical results.Note that the model from Sec. V A does not take intoaccount long-range and high-order correlations. There-fore, the possibility of the description of numerical resultswith this model could not be taken for granted. How-ever, for the considered numerical samples, this descrip- tion is possible and the effect of the simplifications madein Sec. V A is reduced to small modification of h hf andto some changes in the amplitude of magnetoresistance(described with fitting parameter A ).In Fig. 7 (C) and (D) the results for numeric sampleswith fast hopping w = 5 h hf are shown. The results ofthe simulation (blue dots) agree with fitting with Eq. (38)(yellow solid curve), where the value of fitting parameter e γ is 8 . h hf . At this e γ the function R can be described byLorentzian, as shown with green dashed curve in Fig. 7(C) and (D).The provided numeric results show that in some casesthe theory from Sec. V A can be used as a toy model forthe understanding more complex situations when long-range and high-order correlations are required to quan-titatively calculate the magnetoresistance. However, notall the lineshapes that appear in numeric experiments canbe described with this toy model. Let us consider the nu-meric sample shown in Fig. 8. It is constructed from 3 × w = 3 h hf . The blocksare connected with links with slow hopping w = 0 . h hf .The energies of sites are random with normal distribu-tion with the standard deviation ∆ E = T . The averageenergy of sites in A -type block is 1 . T and in B -typeblock it is − . T . The idea behind this sample is as fol-lows. The conductivity of the sample is controlled by theprocess of generation and recombination of electron-holepairs in AB pairs of sites. The structure of the sampleensures that the generated electron and hole will staynear the pair of sites were they are generated for quite along time. However, they are trapped not on a single sitebut on a cluster of nine sites, therefore their spins notonly rotate around local hyperfine fields but also relaxdue to hops between sites with different hyperfine fields.The positive average energy of A -type sites and the neg-ative one of B -type sites ensures that A -type clusterscontain a small number of electrons and B -type clusterscontain a small number of holes. Physically this situa-tion can correspond to small polymer molecules wherethe hopping between monomers of a single polymer is2 FIG. 8: Results for the numeric samples constructed with 3 × ,
2) approximation. (C) magnetoresistance in (7,2) approximation. fast while the hops between different molecules are slow.Note that for the relaxation of spin correlation in the dis-cussed electron-hole pair all the 18 hyperfine fields in twoneighbor blocks are relevant. It cannot be captured withthe toy model from Sec. V A.The structure of the described sample is shown inFig. 8 (A). Fig. 8 (B) shows the magnetoresistance ofsuch samples calculated in the (1 ,
2) approximation thatcorresponds to the model of modified Miller-Abrahamsresistor. In Fig. 8(C) the magnetoresistance calculatedin (7 ,
2) approximation is shown. In both of the cases,the magnetoresistance is averaged over 30 random disor-der configurations (i.e. the random hyperfine fields andsite energies).The magnetoresistance calculated in (1 ,
2) approxima-tion is quite weak ∼
2% and its lineshape is Lorentzian.However, (1 ,
2) approximation is not adequate for the de-scription of the magnetoresistance in these samples be-cause the correlations can easily leave the initial sitesbut are trapped in the clusters. To take this trappinginto account it is required to consider the correlations atthe inter-site distance equal to 7. When these correla-tions are taken into account the estimated magnetore-sistance increases ∼ R ( H ) − R (0) R (0) ∝ H ( | H | + H ) , (39)where H is a fitting parameter. The expression (39) wasused to describe the lineshape of OMAR in a number ofexperimental works [6–10]. I want to stress that the statistics of hyperfine fields isexactly the same in all the considered numerical samples.However, the obtained OMAR lineshapes are differentincluding the two shapes most commonly obtained in ex-periments: H / ( H + H ) and H / ( | H | + H ) . What isdifferent in the numerical samples is the statistics of elec-tron hops. One can conclude that the shape of OMARcontains information about short-range transport in or-ganic materials. VI. DISCUSSION
Up to the recent time, the most known method forthe calculation of hopping transport with the account tocorrelations of filling numbers was the numeric Monte-Carlo simulation. When only the charge correlations areimportant it leads to the correct description of transportprovided that the time of simulation is enough to achievethe averaging. In [25] the Monte-Carlo simulation wasused to prove that CKE approach also leads to correctresults when a sufficient number of correlations are takeninto account. Therefore the Monte-Carlo simulation wasconsidered to be “an arbiter” for CKE approximations.However, the situation is different when the spin cor-relations are relevant. In [13] the Monte-Carlo simula-tion with the semi-classical description of electron spinsin terms of “spin up” and “down” was used to showthe possibility of bipolaron mechanism of OMAR. How-ever, this simplified description cannot include spin rota-tion around hyperfine on-site fields. Naturally, “up” and“down” spins cannot rotate. Therefore the spin relax-3ation in [13] was described by a single relaxation time τ s with phenomenological dependence on the applied mag-netic field.Even when the spin relaxation is reduced to a singletime τ s the semi-classical description of spin does not leadto the correct qualitative estimate of the spin-correlationpart of resistance. In [22] the semi-classical spin correla-tions were considered in the approximation when onlypair correlations in close pairs of sites are taken intoaccount. It leads to correlation corrections of Miller-Abrahams resistors. These corrections are comparedwith similar corrections due to quantum spin correlationsin Sec. V A. The corrections due to charge correlationsobtained in [22] and in the present study agree. How-ever, the corrections due to spin correlations obtained inSec. V A are exactly three times larger than the spin cor-rections to resistances in [22]. I believe that the reasonfor this difference is related to the quantum nature of spinthat was not taken into account in [22]. Note that in [22]only the correlations in close pairs are taken into account.This approximation is insufficient for the description ofthe OMAR shape H / ( | H | + H ) that appears in thesample on Fig. 8 only when long-range pair correlationsare included into the theory.The author does not know about any way to makeMonte-Carlo simulations where electron spin is allowedto have arbitrary direction (instead of only “up” and“down”) and take into account quantum spin correla-tions. Any wavefunction of a single spin 1 / s α = c x ˆ s x + c y ˆ s y + c z ˆ s z , where c x + c y + c z = 1. Therefore itis tempting to describe the electron spins 1 / h s i s j i = s xi s xj + s yi s yj + s zi s zj . In classical statisticsof unit vectors − ≤ h s i s j i ≤
1. The value 1 describesthe vectors that always have the same direction. Thevalue − − ≤ h s i s j i ≤
1. Thevalue − s zi was selected to haveeigenvalues 1 and −
1, however, it corresponds to theactual spin 1 /
2. In quantum mechanics the eigenval-ues of the operator ˆ l of squared angular momentumare equal to l ( l + 1). It the case of the angular mo-mentum 1 / /
4. There-fore the operator (ˆ s i ) = ˆ s zi ˆ s zi + ˆ s xi ˆ s xi + ˆ s yi ˆ s yi is alwaysequal to 3. When two spins are in the singlet statetheir total spin is equal to zero (ˆ s i + ˆ s j ) = 0. Theproduct ˆ s i ˆ s j in this case is well defined and is equal to s i s j = [(ˆ s i + ˆ s j ) − (ˆ s i ) − (ˆ s i ) ] / −
3. When the spinsare in a triplet state, the total momentum of the systemis equal to l = 1 and (ˆ s i + ˆ s j ) = 4 l ( l + 1) = 8. It leadsto s i s j = 1. Therefore the quantum spin statistics is dif-ferent from the statistics of classical vectors, although itcan be described with the classical values: averaged spinpolarizations and their correlations. The Monte-Carlo calculations at least with a naive description of electronspins cannot be used to quantitatively calculate the mag-netoresistance related to the spin correlations and act as“an arbiter” for CKE approach.Although the properties of transport in organic semi-conductors are studied for some time they are notcompletely understood. The low-field mobility in or-ganic semiconductors is often extremely small ∼ − − − cm /V s [35–37]. However, these small values canbe related to the long-range correlations of electrostaticpotential produced by the unscreened molecular dipoles[38–40]. There is an evidence that at small length-scalesthe electron mobility can be much higher [41]. The pro-vided results show that some information about short-range mobility can be obtained from the measurementsof OMAR. The lineshape of organic magnetoresistancedepends on hopping rates. However, it is not related tothe time for an electron to cross the macroscopic sample.What is important is the time τ sep that is required fortwo spins to become separated with sufficient distancethat prevents their meeting before their spin correlationrelaxes.A theoretical estimate of τ sep can be possible in theframework of complex numerical modeling of organicsemiconductors similar to the one made in [29]. How-ever, even the analysis of the provided simplified modelscan give some hints on τ sep and nature of the short-rangetransport. When τ sep is small compared to the period ofprecession in hyperfine field τ sep h hf ≪ τ sep h hf & ∼
10% as shownon Fig. 7 and its lineshape is described by Lorentzianor by Eq. (37). When some of the hops are fast but τ sep h hf & ∼ H / ( | H | + H ) .The present study deals only with the most simplemodel with large Hubbard energy and small applied elec-tric field. In principle, it is possible to generalize CKEtheory to include other cases. In [25] the far from equi-librium CKE that can be applied for high electric fieldsare derived for charge correlations only. In [22] the sit-uation with arbitrary Hubbard energy is considered forclose-range pair correlations with the semi-classical spinmodel. It is shown in [22, 25] that the discussed gen-eralizations make the theory much more complex. Thepresent work shows that different lineshapes of OMARappear even in the simplest model due to the differentproperties of short-range transport. VII. CONCLUSIONS
The system of correlation kinetic equations (CKE) isderived for the spin correlations in materials with hop-ping transport with large Hubbard energy for a smallapplied electric field. The spins are assumed to be con-served in the hopping process and can rotate around on-4site hyperfine fields. The spin degrees of freedom aredescribed with quantum mechanics as averaged productsof spin operators. The derived CKE approach allows todescribe the bipolaron mechanism of OMAR. It is shownthat the shape of the magnetic field dependence of re-sistivity contains information about short-range electrontransport. Different statistics of hopping rates lead todifferent OMAR lineshapes including the empirical laws H / ( H + H ) and H / ( | H | + H ) that are often usedto describe experimental data.The author is grateful for many fruitful discussionsto Y.M. Beltukov, A.V. Nenashev, D.S. Smirnov, V.I.Kozub and V.V. Kabanov. The work was supportedby the Foundation for the Advancement of Theoreti-cal Physics and Mathematics “Basis”. A.V.S. acknowl-edges support from Russian Foundation for Basic Re-search (Grant No. 19-02-00184). Appendix A: Derivation of ( d/dt ) ij term in kineticequations In this appendix I derive the term ( d/dt ) ij s PI in thekinetic equation. It is supposed that sites i and j are in-cluded into set I . For definiteness the site i is consideredto be an A -type site and j to have type B . s PI can be expressed as quantum mechanical average of the operator s PI = (cid:10) ˆ s PI (cid:11) = D ˆ s pi ˆ s qj ˆ s P ′ I ′ E . (A1)Here when the index p corresponds to x , y and z , ˆ s pi isthe spin polarization operator ˆ s αi = a + i,n ( σ α ) nm a i,m . σ α is a Pauli matrix. a + i,n and a i,m are the creation anddestruction operators for electron on site i . Indexes n and m correspond to spin states ↑ or ↓ . The operator ˆ s i is the operator of single occupation and can be expressedas follows ˆ s i = a + i ↑ a i ↑ + a + i ↓ a i ↓ − a + i ↑ a i ↑ a + i ↓ a i ↓ . (A2)The expression (A2) is valid for any type of site, thereforeˆ s j can be expressed in a similar way.The transitions between sites i and j can be describedby the operator b H ij (it is the term in Hamiltoinan relatedto these transitions). b H ij = t ij b τ ij Φ ij + t ji b τ ji Φ ji , τ ij = a + i ↑ a j ↑ + a + i ↓ a j ↓ . (A3)The operator Φ ij describes the interaction with phononsthat appears in the transition term of Hamiltonian b H ij after the polaron transformation [42].With the approximations corresponding to hopping transport the time derivative of operator ˆ s PI can be expressedas follows (cid:18) ddt (cid:19) ij ˆ s PI = − ~ Z t −∞ Dh b H ij ( t ) , h b H ij ( t ′ )ˆ s PI ( t ) iiE ph dt ′ . (A4)Here square brackets denote commutator [ b A, b B ] = b A b B − b B b A . h ... i ph means the averaging over phonon variables withequilibrium distribution of phonons. Note that even with the simplifications made eq. (A4) contains not only hoppingterms but also terms corresponding to quantum mechanical perturbation of electron wavefunctions on sites i and j and to exchange interaction between electrons that is neglected in our study. In further calculations I keep only theterms related to hopping process that are proportional to the hopping rates W ij = 1 ~ | t ij | (cid:28)Z −∞ Φ ji ( t )Φ ij (0) e ( i/ ~ )( ε j − ε i ) t + Φ ji (0)Φ ij ( t ) e ( i/ ~ )( ε i − ε j ) t dt (cid:29) ph , (A5) W ji = 1 ~ | t ij | (cid:28)Z −∞ Φ ij ( t )Φ ji (0) e ( i/ ~ )( ε i − ε j ) t + Φ ij (0)Φ ji ( t ) e ( i/ ~ )( ε j − ε i ) t dt (cid:29) ph . (A6)The term in ( d/dt ) ij b s PI proportional to W ji and related to hops i → j is equal to − W ji ( b τ ij b τ ji b s PI + s PI b τ ij b τ ji ) /
2. Theterm related to j → i hops is W ij b τ ji b s PI b τ ij . Here I took into account that site i cannot be double-occupied and site j cannot have zero electrons.The following computation is quite cumbersome but straightforward operator algebra. There are two useful relationsthat make it simpler ˆ τ ij ˆ τ ji = ˆ s i ˆ s j − P α ˆ s αi ˆ s αj , ˆ s αi ˆ s βi = δ αβ ˆ s i + iǫ αβγ ˆ s γi , (A7)where ǫ αβγ is Levi-Civita symbol.5The operator algebra yields (cid:18) ddt (cid:19) ij ˆ s i ˆ s j ˆ s P ′ I ′ = 2 W ij ˆ n j ↑ ˆ n j ↓ (1 − ˆ n i ↑ )(1 − ˆ n i ↓ )ˆ s P ′ I ′ − W ji ˆ s i ˆ s j − X α ˆ s αi ˆ s αj ! ˆ s P ′ I ′ , (A8) (cid:18) ddt (cid:19) ij ˆ s i ˆ s αj ˆ s P ′ I ′ = W ji (cid:0) ˆ s αi ˆ s j − ˆ s i ˆ s αj (cid:1) ˆ s P ′ I ′ , (A9) (cid:18) ddt (cid:19) ij ˆ s αi ˆ s j ˆ s P ′ I ′ = W ji (cid:0) ˆ s i ˆ s αj − ˆ s αi s j (cid:1) ˆ s P ′ I ′ , (A10) (cid:18) ddt (cid:19) ij ˆ s αi ˆ s βj ˆ s P ′ I ′ = − W ij δ αβ ˆ n j ↑ ˆ n j ↓ (1 − ˆ n i ↑ )(1 − ˆ n i ↓ )ˆ s P ′ I ′ + W ji δ αβ ˆ s i ˆ s j − X α ˆ s γi ˆ s γj ! ˆ s P ′ I ′ − W ji (cid:16) ˆ s αi ˆ s βj − ˆ s βi ˆ s αj (cid:17) ˆ s P ′ I ′ (A11)The quantum mechanical averaging of these equations leads to eqs. (13-16) from the main text. [1] A. Miller and E. Abrahams, Impurity conduction at lowconcentrations, Phys. Rev. , 745 (1960).[2] B. I. Shklovskii and A.L. Efros, ”Electronic Properties ofDoped Semiconductors” (Springer, Berlin, 1984).[3] S. R. Forrest, Nature (London) , 911 (2004)[4] C.J. Brabec, N. S. Sariciftci, J.C. Hummelen, Adv.Funct. Mater. , 15 (2001)[5] J. Kalinowski, M. Cocchi, D. Virgili , P. Di Marco, V.Fattori, Chem. Phys. Lett.
710 (2003)[6] O. Mermer, G. Veeraraghavan, T. L. Francis, Y. Sheng,D. T. Nguyen, M. Wohlgenannt, A. Kohler, M. K. Al-Suti, M. S. Khan, Phys. Rev. B , 205202 (2005)[7] Y. Sheng, T.D. Nguyen, G. Veeraraghavan, O. Mermer,M. Wohlgenannt, S. Qiu, U. Scherf, Phys. Rev. B ,045213 (2006)[8] T.D. Nguen, Y.G. Sheng, J.Rybicki, G.Veeraraghavan,M.Wohlgenannt, Journ. of mat. chem. , 103715 (2008)[10] L. Yan, M. Wang, N.P. Raju, A. Epstein, Loon-Seng Tan,A. Urbas, Long Y. Chiang, B, Hu, J. Am. Chem. Soc. , 7 (2012)[11] Hongbo Gu, Xi Zhang, Huige Wei, Yudong Huang, Suy-ing Wei, Zhanhu Guo, Chem. Soc. Rev. , 5907 (2013)[12] V.N. Prigodin,J.D. Bergeson, D.M. Lincoln, A.J. Ep-stein, Synthetic Metals ,216801 (2007).[14] N. J. Harmon and M. E. Flatte, Phys. Rev. Lett. ,186602 (2012)[15] N. J. Harmon and M. E. Flatte Phys. Rev B , 075204(2012)[16] N. J. Harmon and M. E. Flatte Phys. Rev B , 245213 (2012)[17] F.J. Yang, W. Qin, S.J. Xie, The Journal of ChemicalPhysics 140, 144110 (2014)[18] N. Gao, L. Li, N. Lu, C. Xie, M. Liu, H. Bassler, Phys.Rev. B , 075201 (2016)[19] N.Lu, N. Gao, L. Li, M. Liu, Phys. Rev. B , 165205(2017)[20] A. Larabi, and D. Bourbie, J. Appl. Phys. , 085502(2017)[21] A.V. Shumilin, V.V. Kabanov, Phys. Rev. B , 014206(2015)[22] A.V. Shumilin, V.V. Kabanov, V.I. Dediu, Phys. Rev. B , 094201 (2018)[23] O. Agam, I.L. Aleiner, Phys. Rev. B , 224204 (2014)[24] O. Agam, I.L. Aleiner, B. Spivak, Phys. Rev. B ,100201(R) (2014)[25] A.V. Shumilin, Y.M. Beltukov, Phys. Rev. B , 014202(2019)[26] A.V. Shumilin, Y.M. Beltukov, Phys. Solid State, , 115204 (2017)[30] Z.G. Yu, F. Ding, H. Wang, Phys. Rev. B , 205446(2013)[31] N.J. Harmon, M.E. Flatte, Phys. Rev. Lett. , 176602(2013)[32] V.V. Mkhitaryan, V.V. Dobrovitski, Phys. Rev. B ,054204 (2015)[33] T. Damker, H. Bottger, V. V. Bryksin, Phys. Rev. B , ,075409 (2019).[35] L.B. Schein, A. Rosenberg, S.L. Rice, Journ. of Appl.Phys. , 4287 (1986)[36] P.M. Borsenberge, E.H. Magin, M. Van Der Auweraer,F.C. De Schryve, Phys. Stat. Sol. (a) , 9 (1993)[37] P.M. Borsenberger, L. Pautmeier, H. Bassler, J. Chem.Phys. , 1258 (1991)[38] S.V. Novikov, A.V. Vannikov, J. Phys. Chem., , 14573(1995) [39] D.H. Dunlap, P.E. Parris, V.M. Kenkre, Phys. Rev. Lett.,
542 (1996)[40] S.V. Novikov, Russian Journal of Electrochemistry,
388 (2012)[41] A.J. Drew, F.L. Pratt, J. Hoppler, L. Schulz, V. Malik-Kumar, N.A. Morley, P. Desai, P. Shakya, T. Kreouzis,W.P. Gillin, K.W. Kim, A. Dubroka, R. Scheuermann,Phys. Rev. Lett.100