Spin-conserving Boltzmann theory for carriers and excitons in organic semiconductors

The rise of organic electronics calls for versatile modeling tools. In this context, we develop a semiclassical Boltzmann theory that describes transport and excitonic processes in crystalline organic semiconductors on equal footing. The generation of singlet and triplet excitons out of the ground state, their formation from free electrons and holes, the reverse processes, as well as the fusion and fission of excitons are included. The corresponding scattering integrals respect spin conservation, which requires matrix-valued distribution functions. They also include fermionic and bosonic many-particle effects such as Pauli blocking. We employ a multipole expansion of the distribution functions, where quadrupolar terms turn out to be essential for the triplet excitons. This work provides a basis for the modeling of organic solar cells, in which excitonic processes are crucial for the performance. Moreover, the theory is of general interest for transport and transitions of multiple (quasi-) particle species carrying spin in nonequilibrium systems.

However, the design of electronic properties requires a deeper theoretical understanding of the charge-carriertransport mechanisms. Several aspects make this a challenging goal for organic materials: On the one hand, the lattice is rather soft, i.e., the predominantly intermolecular phononic modes have low energies. They are also relatively strongly coupled to the charge carriers. This interaction is of short range with the largest contribution being the coupling of the electronic occupation of a certain molecular orbital to vibrations of the same molecule. This diagonal or Holstein coupling [13] results in the formation of so-called small polarons [14,15]. There is also a sizable coupling of intermolecular tunneling amplitudes to vibrations, called nondiagonal or Peierls coupling [16][17][18][19][20]. For rubrene, Ordejón et al. [21] have found that the Peierls coupling is indeed important. These couplings lead to strong polaronic effects [20,[22][23][24][25][26], which limit the carrier mobility.
On the other hand, creation, decay, and recombination of excitons as well as the fission of spin-singlet excitons into triplet excitons and the reverse fusion process are crucial for organic photovoltaic and light-emitting devices [27][28][29][30][31][32]: Lifetimes of triplet excitons can be six orders of magnitude longer than the ones of singlet excitons because their transition into the ground state, i.e., the recombination process, requires a spin flip [33]. Due to the * carsten.timm@tu-dresden.de spin exchange energy, the singlet-exciton energy can be higher than twice the triplet-exciton energy [34]. This enables singlet fission into two triplet excitons, which is spin allowed. This fission process can be the key mechanism to induce transport of excitation energy by triplet excitons over long distances, which is of major importance for photovoltaic applications [29,30]. Note that the photogeneration in such devices is dominated by internal interfaces while the exciton diffusion is a bulk effect. We return to this point in the conclusions.
We aim at a unified description of transport and the aforementioned local excitonic processes. A relevant property of organic materials is the typically weak spinorbit coupling, which ensures that spin is conserved to a good approximation. The description of transport thus has to respect spin rotation symmetry, in particular if singlet and triplet excitons are involved. We note that any violation of spin conservation, i.e., any relaxation of spin, can easily be incorporated by an additional scattering integral.
The excitons in typical devices are far from equilibrium. The quantum master equation constitutes a powerful tool for the description of nonequilibrium processes in quantum systems [35]. It is the equation of motion for the reduced density matrix or statistical operator ρ of an open system. This equation is linear in the density matrix since it is derived from the Schrödinger equation, which is linear in the state vector. Linearity is necessary to preserve the unit trace of the density matrix. However, one can relax this trace condition and interpret Tr ρ as a concentration of quasiparticles. The matrix structure of ρ then encodes the distribution of spin orientations, where ρ is a (2s + 1) × (2s + 1) matrix for spin-s quasiparticles. Furthermore, the equation of motion can then be nonlinear, for example describing generation and decay processes. A generalized quantum master equation in this sense was applied to local excitonic processes, such as fission and fusion, in organic crystals as early as 1969 by Johnson and Merrifield [36]. Essentially the same equation was proposed, under the name of stochastic Liouville equation, in the context of quantum optics [37] and for the description of chemical reactions involving spin [38][39][40]. More recently, Schellekens et al. [41] have applied it to excitonic processes in order to describe organic magnetoresistance.
In case of the stochastic Liouville equation, the nonlinearity is given by a ρ-independent source term [36,37,41]. The equation takes the form The first term on the right-hand side describes the unitary time evolution determined by the Hamiltonian H. The second term implements decay processes and is naturally linear in the density. Here W is a decay rate and P is a projection operator onto the spin state from which a decay is allowed by spin conservation. The specific anticommutator structure was shown by Haberkorn [39] to ensure the positivity of ρ-since the eigenvalues of ρ are interpreted as concentrations they must not become negative [42]. The third term describes generation processes and is of order zero in ρ since the generation rate is assumed to be independent of the density. The interpretation of ρ as a density also permits terms of higher order. Consider, for example, a pair-annihilation process of the type A+A → 0. For spinless particles, ρ just becomes the number density n and the process is naturally described by the equation This is of course well known in chemical kinetics [43,44]. Nonlinear terms also result from quantum statistics, which show that the rate of a scattering process also depends on the occupation of the final state, e.g., due to the Pauli principle.
The question arises what the corresponding equation looks like for particles with spin. The right-hand side then contains a product of the matrix ρ with itself and the question is how this product must be constructed so that spin is conserved. More complicated processes involving multiple species, e.g., A + B → C, can be described using multiple matrix-valued densities and we have to construct spin-conserving terms out of matrices of generally different sizes (2s + 1) × (2s + 1). This issue is central for the description of excitonic processes. For spin-1/2 quasiparticles, i.e., electrons and holes, the matrix-valued densities are 2 × 2 matrices acting on their spin Hilbert space. For triplet excitons with spin 1, they are 3 × 3 matrices, whereas for singlet excitons they reduce to scalars.
We combine the spin-rotation-invariant description of local excitonic processes with a spin-rotation-invariant description of carrier and exciton transport. We are interested in bandlike transport in clean organic materials at not very low temperatures. Characteristic quantum effects such as weak localization and universal conductance fluctuations are thus expected to be irrelevant and a semiclassical description is justified. Our goals thus require the derivation of SU(2) spin-rotation-invariant Boltzmann-type kinetic equations for the various quasiparticle species and of scattering integrals describing transitions between them.
As noted above, polaronic effects are typically strong in organic materials [20,[22][23][24][25][26]. The coupling to the phonons causes the charge carriers and also excitons to dress as small polarons [15]. Unlike for inorganic materials, there is no well-developed semiclassical transport theory for small polarons in organic semiconductors. This is a promising topic for future research, where one obstacle will likely be the inclusion of sizable Peierls coupling. Numerical results by Hutsch et al. [26] indicate that the phononic modes can be usefully separated into fast modes, which lead to polaron formation similar to inorganic materials, and slow modes, which can be treated as quasistatic disorder. This suggests that a semiclassical description should be possible, where the effect of the slow modes is incorporated by scattering integrals. Since polaronic effects are not at the focus of the present paper we assume that the polarons are adiabatically connected to the undressed quasiparticles, i.e., that polaronic effects only renormalize the parameters appearing in our description but do not change its structure.
Spin-rotation invariance requires the Boltzmann theory to contain the matrix-valued densities [45][46][47][48][49][50][51][52][53]. We emphasize that it is not sufficient to use kinetic equations for the densities in each spin channel, say spin-up and spin-down electrons, since this is equivalent to assuming the matrix-valued densities to be diagonal in the standard spin basis. This precludes the description of spin polarization in the x and y directions. It is also clear that the scattering integrals are significantly more complicated than for spinless particles since they depend on matrix-valued functions of various dimensions and are themselves matrices.
Our program is of general interest beyond the field of organic electronics as the principles developed here are useful for any system showing scattering of and reactions between spin-carrying particles, for example in spintronics and spin chemistry. This paper is organized as follows. In Sec. II, we set up the spin-conversing kinetic equations for charge carriers and excitons in uniform electric and magnetic fields. The treatment of scattering and transitions between quasiparticle species (excitonic processes) is then performed in Sec. III. The theoretical framework is illustrated for a simple model in Sec. IV. Conclusions are drawn in Sec. V.

II. SPIN-CONSERVING KINETIC EQUATIONS
In this section, we present the Boltzmann transport equations for electrons, holes, and singlet and triplet excitons including static, uniform electric and magnetic fields. We concentrate on the drift terms conventionally appearing on the left-hand side of the Boltzmann equation and leave the collision terms for Sec. III. Following Ref. [48], we write the general Boltzmann equation for the density f (r, k, t) as where q ∈ {−e, 0, e} is the charge of the quasiparticles under consideration and the hat • denotes matrices acting on spin space. The density f (r, k, t) is a (2s+1)×(2s+1) matrix, where s is the quasiparticle spin. Note that for s > 0, the matrix structure leads to an additional commutator term, which describes the dynamics in spin space [48]. For example, for spin-1/2 quasiparticles in a magnetic field B, the Hamiltonian H contains a Zeeman term proportional to B · σ, where σ is the vector of Pauli matrices σ 1 , σ 2 , σ 3 . We suppress the arguments r, k, and t from now on. I[ f ] contains all scattering integrals affecting the distribution. Since we will later consider nonlinear collision integrals, the normalization of f is important. We take the eigenvalues of f to be the occupation numbers per quantum state. For example, the density for all r and k describes a filled electron band. This convention will be useful for the description of Fermi and Bose statistics in the collision integrals. As noted above, Tr f is not restricted to unity and f is not a density matrix in the sense of a statistical operator [41]. The velocity in Eq. (3) is given by the derivative of the Hamiltonian acting on spin space with respect to mo- and thus, in principle, obtains a matrix structure. This structure is made nontrivial by those terms in the Hamiltonian that contain both the momentum and the spin, i.e., by spin-orbit coupling. Since we are primarily interested in organic materials, which show weak spin-orbit coupling, we neglect them. The velocity v then becomes proportional to the identity matrix and can be treated as a scalar [54]. The dispersion relation is the only relevant momentum-dependent term in the Hamiltonian so that we can express v as the group velocity v = 1 ∇ k (k) (6) in terms of the quasiparticle dispersion (k).

A. Electrons
The charge carriers in organic semiconductors are electrons and holes. For simplicity, we restrict ourselves to a model with one valence band and one conduction band; the generalization to multiple bands is straightforward. The spin 1/2 of the electrons and holes is essential to understand the related excitonic processes.
The identity matrix σ 0 and the three Pauli matrices σ x , σ y , and σ z form a basis of the space of Hermitian 2 × 2 matrices. Hence, for spin-1/2 quasiparticles, the 2 × 2 matrix f can be expanded as Then n = Tr f is the concentration of quasiparticles, summed over spin orientations, and c = Tr f σ is twice the spin density (setting = 1). Note that f must be Hermitian to ensure real densities. We will call relations of the form of Eq. (7) multipole expansions.
We start with the description of electrons in the conduction band with the Hamiltonian with the electronic Landé factor g el and the Bohr magneton µ B . Neglecting spin-orbit coupling, the bare electron Hamiltonian is H e0 = e (k) σ 0 , with the dispersion e (k). The Zeeman term remains the only part with nontrivial matrix structure. We decompose the electronic distribution function into the equilibrium distribution f e0 and a deviation g e according to We assume the distribution function to be close to equilibrium because the equilibrium distribution is perturbed only little by (a) a weak electric field E and (b) decomposition of a low concentration of excitons. The weak electric field also means that terms involving the product of E and g e are small of higher order and can be neglected. This is the assumption of linear response. It would be straightforward to drop this approximation. The multipole expansion of the deviation is similar to Eq. (7). In equilibrium, the distribution is spatially uniform since the external fields and the scattering term are assumed not to depend on position r, and can thus be written as This yields f e = 1 2 (n e0 +ñ e ) σ 0 + (c e0 B +c e ) · σ .
The equilibrium distribution can be derived starting from the specific case of B e 3 , where e 3 is the unit vector in the z direction. Here the Hamiltonian takes a diagonal form and so does the equilibrium function. It is thus given by a linear combination of σ 0 and B σ z . The rotation to an arbitrary direction of the magnetic field is then described by the transformation B σ z → B · σ. The induced magnetization is along B since we assume a magnetically isotropic medium [55]. The equilibrium term n e0 corresponds to the Fermi-Dirac distribution of electrons with energy e (k), chemical potential µ, and temperature T , In the linear-response approximation, the Boltzmann equation for the electrons reads as we obtain (30)

C. Singlet excitons
Singlet excitons are spin-0 quasiparticles. Consequently, the density f s and Hamiltonian H s are scalars so that the commutator term in Eq. (3) vanishes. This simplifies the Boltzmann equation to Again, we use the group velocity with regard to the dispersion relation of the singlet excitons, s (k). The description of excitons by a band structure s (k) and an associated velocity has a long history [56][57][58] and has also been applied to organic semiconductors [59][60][61]. The excitonic band structure in pentacene has been studied experimentally by Schuster et al. [62]. One should note that the band width and velocity of excitons can be larger than the ones of the charge carriers since the motion of excitons can take place through pure energy (Förster) transfer, without tunneling of electrons.
Similarly to the linear-response theory for electrons, we decompose the distribution function into an equilibrium distribution f s0 and a deviation g s , However, here the deviation is not related to the electric field and is not necessarily small since we also consider situations with optically excited excitons far from equilibrium. We assume f s0 to be uniform. In equilibrium, excitons in semiconductors follow the Bose-Einstein distribution where is the exciton energy and µ is their chemical potential [63]. Inserting Eq. (32) into Eq. (31) gives The scalar form of the Boltzmann equation for singlet excitons allows one to obtain the solutions analytically. This requires specific assumptions for the scattering term. In a simple case, we can neglect interactions with other particle species and use the relaxation-time describe the scattering. R s is a possibly momentum-dependent relaxation rate. This transforms the Boltzmann equation into a homogeneous partial differential equation for g s , It is of first order and the coefficients are constant for fixed k. The analytical solution can be derived using Fourier transformation with respect to r. The result is with the integration constant g 0 s . This equation expresses an exponential dampening of the deviation that is superimposed onto a ballistic divergence in position space. This dynamics agrees with the expected behavior. We can take this as a confirmation of the form of the Boltzmann equation derived above.

D. Triplet excitons
Triplet excitons are spin-1 quasiparticles. Hence, their distribution function f t is a 3 × 3 matrix acting on the spin-1 Hilbert space and their Boltzmann equation takes the form with the triplet-exciton dispersion t (k). As before, we decompose the distribution function into an equilibrium function and a deviation, f t = f t0 + g t , and assume f t0 to be uniform. Inserting this into Eq. (37) yields where the equilibrium commutator [ H t , f t0 ] again vanishes. We choose a basis of the space of Hermitian 3 × 3 matrices in a similar fashion as for the 2 × 2 case. This includes as the monopole term and the three-dimensional representation of SU(2), as dipole terms. The normalization condition for all basis matrices is Tr S 2 i = 2. The remaining five matrices can be interpreted as quadrupole terms and are arranged as a five-component vector Q. We use the Hermitian matrices [64] which are also normalized such that Tr Q 2 i = 2. The matrices Q 1 , Q 2 , and Q 3 transform as t 2g under the cubic point group O h , while Q 4 and Q 5 correspond to e g .
The multipole expansion of f t is with the coefficients n t ∈ R, c t ∈ R 3 , q t ∈ R 5 . The related physical quantities are the number density √ 6 n t , the spin density c t , and the quadrupole density q t . This is a good place to also introduce the Cartesian components of the quadrupole tensor operator, which form the symmetric and traceless matrix Importantly, this matrix acts on the space of Cartesian spin components and each component is a 3 × 3 matrix on the Hilbert space of a spin of length 1. Using the Cartesian form, the quadrupolar term in the distribution function f t can be rewritten as where the trace is over Cartesian components and Next, we motivate the form of the Boltzmann equation for these physical quantities. The Hamiltonian for the triplet excitons reads as where t (k) is the dispersion in the absence of a magnetic field (see the comments on excitonic band structures in Sec. II C) and g tr is the Landé factor of the triplet excitons, which is assumed to be scalar. A triplet exciton is of course charge neutral but since it generically consists of an electron and a hole with different Landé factors its total magnetic moment is nonzero. Additional terms linear in quadrupole matrices Q i have been proposed in the literature [65][66][67]. This magneto-crystalline anisotropy is due to spin-orbit coupling and is therefore weak in organic semiconductors. We omit these terms for simplicity but implementing them in the Hamiltonian is straightforward.
If the external magnetic field is aligned along the z axis, the Zeeman term simplifies to g tr µ B B S 3 with B = |B|. The Hamiltonian is diagonal in this case. We use this specific situation to derive the form of the equilibrium distribution function f t0 . The equilibrium distribution is the Bose-Einstein distribution and can immediately be evaluated for the diagonal components. The resulting f t0 has a diagonal matrix structure as well, A rotation yields the general description for an arbitrary direction of B. The matrix S 3 generalizes to S projected onto the direction of B. Thus, the general form of the equilibrium distribution function is given by The deviation g t is of the most general form of the multipole expansion, The velocity is given by a derivative of the Hamiltonian of triplet excitons in analogy to Eq. (5) and thus is a 3 × 3 matrix. Since we assume spin-orbit coupling to be negligible, we can approximate it by a scalar Kinetic equations for the coefficients can be obtained but the derivation is made complicated by quadrupolar terms. The details are relegated to Appendix A. Defining the differential operator and expanding the scattering integral into multipoles, we find the coupled equations This complicated system cannot be simplified without further knowledge about the scattering terms. A detailed investigation of these integrals is presented in the following section.

III. SPIN-CONSERVING SCATTERING INTEGRALS
In this section, we construct spin-conserving scattering integrals for the Boltzmann equations for electrons, holes, and excitons. These scattering integrals should not only describe collisions between quasiparticles and scattering off disorder but, importantly, also generation, decay, and transitions between excitons. An overview of the processes studied here is given in Table I. The selection and order of these processes as well as the presentation in this section are partly pedagogical. We progress from conceptually simpler to more complicated cases, where essentially every subsection introduces an additional aspect, usually having to do with conservation of spin or indistinguishability of quasiparticles. Additional processes can be treated in an analogous manner, e.g., the reverse of some of the processes in Table I and processes for holes instead of electrons. Which processes are important of course depends on specific materials. Disorder scattering is expected to be always present. Bipolar semiconductors in which the carrier energies are sufficiently high in comparison to the exciton binding energies will also show electron-hole binding and unbinding, e + h ↔ s or e + h ↔ t. For systems with photogenerated singlet excitons, e.g., photovoltaic devices, the fission s → 2t and the fusion 2t → s are crucial, as noted in Sec. I. Note that the derivations in this section do not require any quasiparticle species to be close to equilibrium. Scattering off disorder is described by collision integrals that are linear in the distribution function [68,69]. Possanner and Negulescu [52] have studied such terms for both spin-conserving and spin-flip scattering, from a mathematical point of view. El Hajj [53] has studied spin-diffusion models derived from the resulting linear matrix Boltzmann equation in detail.
Our main interest is in the description of transitions between quasiparticle species. The required collision integrals are necessarily nonlinear. Quantum statistics imply that the transition rates depend on the occupation of the final states, contrary to what is assumed in the stochastic Liouville equation. For fermions, complete occupation of the final state prevents a process due to the Pauli principle, whereas for bosons, a large occupation of the final state enhances the rate. In the context of excitonic processes, this has been studied by Bisquert [70], albeit not using an SU(2)-invariant formalism.
It is useful to first outline the principles behind the construction of the collision integrals. Generation, decay, and transitions between excitons involve quasiparticles on different energy levels. This is reminiscent of laser theory, in which transition rates are derived by considering the occupation of the initial and final states. However, unlike in laser theory, these transition rates depend on momentum. Thus, they take a form similar to the basic scattering term in Boltzmann theory, where f is a distribution function. The two integrals describe in-scattering and out-scattering, respectively. The processes considered here are more complicated in a number of ways. First, our Boltzmann equations contain matrix-valued density functions. The idea is to derive the matrix form from the special case of z -polarized particles by means of a rotation into a general direction, like discussed for the kinetic term in Sec. II. Second, the initial and final states can consist of more than one particle. We assume the initial state A to contain particles of species a 1 , a 2 , . . . (some of which may be the same) and the final state B to contain particles of species b 1 , b 2 , . . . (some of which may be the same). The integrand of the scattering integral is then proportional to where N a1 etc. are the occupation numbers of states and the sign in In the following, we will progress from simple to more complicated cases. To fix the notation, we first discuss the essentially trivial case of elastic scattering of electrons off nonmagnetic impurities. It will be beneficial to write the matrix-valued distribution function of the electrons in components, Since nonmagnetic scattering does not flip the spin, the matrix-valued collision integral appearing in the Boltzmann equation for f e is diagonal, The collision integral for spin-up electrons can be written as [68,69] I e↑ e→e (k) = and analogously for spin-down electrons. Detailed balance for elastic scattering implies that W e →e (k , k) = W e→e (k, k ) and thus (72) and analogously for spin down.
In the following, we employ a short-hand notation that suppresses (a) integrals over momenta k , k , . . . and (b) momentum arguments. Functions depending on k, k , and k are denoted by no, one, and two primes, respectively. Hence, we write the collision integrals for elastic impurity scattering as From Eq. (20), we then obtain the collision integrals for the electron number density, and for the density of the z -component of spin, The crucial next step is to reconstruct the SU(2)-invariant form of the collision integral for the spin by rotating all quantities into an arbitrary direction. In the present case, this is simple. We recognize that the lefthand side and both summands on the right-hand side are z -components of vectors. The SU(2)-invariant form is thus We can now use the multipole expansions in Eqs. (7) and (20) to express the collision rate in terms of the matrix densities as The case of holes is of course analogous. The result for excitons is also analogous: For example for singlet excitons, the collision integral is where we have used that excitons are bosons. The second-order terms again cancel and we obtain A. Generation and decay of singlet excitons Singlet excitons are described by the scalar distribution function f s , so SU(2) invariance is automatically maintained. The transition rate describing the generation of singlet excitons only depends on the occupation of the final state. Since singlet excitons are bosons the resulting transition rate reads as Here |0 → s is a short-hand way to denote that the singlet exciton is excited starting from the ground state with all single-electron states below (above) the Fermi energy occupied (empty). This is of course not a spontaneously occurring process but requires coupling to the radiation field, which is suppressed in this notation. We use the symbol W |0 →s with a tilde for the rate; this is an effective rate with properties of the radiation field such as occupation numbers of photonic states absorbed. For the decay of a singlet exciton into the Fermi sea, we analogously write These transition rates can be added to the right-hand side of the Boltzmann equation (31) for the singlet excitons. The corresponding processes for triplet excitons are of course described analogously.

B. Binding and unbinding of electrons and holes in a singlet state
The transition between unbound electron-hole pairs and singlet excitons involves three different particle species, which results in a description from three different perspectives. We discuss these in turn.

Point of view of the singlet exciton
In the case of z -polarized particles, there are only two different orientations of the electron and hole spins. To form a singlet exciton, a spin-up electron has to be paired with a spin-down hole or vice versa. Hence, the transition rate takes the form Here functions depending on the momenta k = k + k and k = k − k are described by two and three primes, respectively. Note that the momentum arguments of the electron and the hole distribution functions could be interchanged as all possible combinations are covered by integrating over k . The origin of the factor 1/2 is that the electron-hole states |↑↓ and |↓↑ can be written as |↑↓ = (|ψ s + |ψ t0 )/ √ 2 and |↓↑ = (− |ψ s + |ψ t0 )/ √ 2, respectively, where |ψ s = (|↑↓ − |↓↑ )/ √ 2 is the singlet state and |ψ t0 = (|↑↓ + |↓↑ )/ √ 2 is the m = 0 triplet state. Hence, the probability of, say, |↑↓ being in the singlet state is 1/2. (1 + f s ).
Now it is possible to generalize this transition rate to arbitrary spin directions as there is a unique SU(2)-invariant expression that reduces to Eq. (84) for z -polarized spins: This transition rate I s e+h→s ≡ I fs e+h→s appears in the Boltzmann equation (31) for the scalar distribution function f s of singlet excitons and is thus itself a scalar. It can also be expressed in a basis-independent form in terms of matrix-valued distribution functions as [71] I fs e+h→s = where S e and S h are the spin operators of the electron and the hole, respectively. It expresses the fact that an electron and a hole can only form a singlet exciton if they are in a relative spin-singlet state. The structure of the momentum-dependent transition rate is reminiscent of the Haberkorn approach [36,38,39]. As mentioned in Sec. I, the anticommutator structure ensures positivity of the density matrix on the electron-hole product space. The trace in Eq. (86) sums over all contributions allowed by spin conservation and leads to a scalar transition term.
It is useful to restore the momentum arguments and integrals at this point. For the process e + h → s from the point of view of the exciton, the outer momentum k is the one of the exciton, the momentum k of the electron needs to be integrated over and the momentum k = k − k of the hole is then fixed by momentum conservation. This leads to The breaking of singlet excitons into unbound electrons and holes can be described in the same way. For this reverse process only the sign and the occupation numbers have to be changed. Now the final state is fermionic, which yields in the multipole form and as the basis-independent expression; compare Eqs. (7) and (20).

Point of view of the electron
We now consider the perspective of the electron (the hole case is of course analogous). The electron's momentum is the momentum argument k of the solution of the Boltzmann equation. The transition rates for spin-up and spin-down electrons read as respectively. The transition rate for the electron number density is the sum of these two transition rates, while the transition rate of the density of the z-component of the spin is their difference. After generalizing to an arbitrary polarization direction, we obtain and The transition rate of the electron number density is the negative of the transition rate from the point of view of the singlet exciton given in Eq. (85). According to Eqs. (7) and (20), the matrix-valued transition term then reads as where Tr h is the trace over the hole sector. To restore the momenta and integrals, we note that the outer momentum k is the one of the electron, the hole momentum k is integrated over, and the exciton momentum k = k + k is then fixed. The rate reads as These transition rates can be compared to the transition rates from the point of view of the singlet exciton. This yields as the creation of one singlet exciton results in the annihilation of a free electron and a free hole. Furthermore, the spin transition rates satisfy which shows that spin is conserved, as the singlet exciton does not posses an overall spin. The reverse unbinding process can be described analogously. It can also be read off from Eq. (95), essentially by reversing the sign and interchanging occupied and empty states. The resulting rate reads as

C. Binding and unbinding of electrons and holes in a triplet state
Like for electrons and holes, it is beneficial to write the matrix-valued distribution function of the triplet excitons in components, The indices of the matrix elements denote the magnetic quantum number of the triplet exciton, where we usē 1 = −1. For a spin along the z axis, only the diagonal elements of the matrix are nonzero so that in the multipole representation it takes the form cf. Eq. (48). This allows us to express the elements of the density matrix in terms of the diagonal elements of the multipole expansion and vice versa,

Point of view of the triplet exciton
The transition term for the triplet-exciton distribution function f t also has to be a 3 × 3 matrix on the spin-1 Hilbert space. For the spin along the z axis, spin conservation implies that we only need to consider the transition rates of the diagonal components, which read as Recall that three primes denote a dependence on k = k − k . The next step is again to construct the SU(2)invariant transition rates by rotating the spin into a general direction. In analogy to the multipole expansion of the distribution function in Eq. (48), we expand the transition term as In these equations, the terms of second order in the coefficients represent the classical part since they do not depend on the occupation of the final state. The terms of third order are the quantum corrections due to the Bose-Einstein statistics of the excitons. Unlike for the previous examples, the SU(2)-invariant generalization of Eqs. (109)-(111) is not obvious since these transition rates contain the quadrupole coefficient q t5 . The transformation of the coefficient vector q t under rotations is not clear. On the other hand, the Cartesian quadrupole density q t defined in Eq. (52) transforms like a matrix under rotations. To make use of this, we have to identify the factors multiplying q t5 as components of Cartesian matrices and write Eqs. (109)-(111) as components of proper matrix products. This can be done by considering the two-particle state of an electron and a hole. As the unbound electron-hole pair can be of triplet character, it generally has a quadrupole moment. The projection operator projects onto the triplet subspace, in analogy to the singlet projection operator P s defined in Eq. (87). Since this subspace is three dimensional the projected density matrix P t f e ⊗ f h P t can be written as a 3 × 3 matrix with respect to a suitable basis of this subspace.
To be consistent, we choose the canonical spin-1 basis The notation [•] 3 refers to the corresponding matrix. The density matrix for an electron-hole pair in a triplet state is thus written as (113) The resulting matrix f eh is an operator on the Hilbert space of spin-1 particles and can thus be expressed in terms of the nine basis matrices of the multipole expansion. For a spin along the z axis, only three coefficients are nonzero, f eh = n eh S 0 + c z eh S 3 + q eh5 Q 5 , where As electron and hole together have the overall momentum k due to momentum conservation, f eh and the corresponding multipole coefficients do not carry any primes.
These relations allow us to construct the SU(2)-invariant generalizations without ambiguities. For terms containing quadrupole moments, Eq. (52) should be noted. The transition rate of the triplet exciton number density is where the classical part is proportional to n eh . For the last term, the identity Tr q eh q t = 2 q eh · q t has been used. The spin transition rate becomes where the classical part is proportional to c eh . The transition rate of the quadrupole part has to be traceless and symmetric in order to ensure that the full triplet-exciton transition rate is Hermitian. This requires the subtraction of the diagonal parts of matrices obtained by a naive restoration of SU(2) symmetry and yields with the classical part being proportional to q eh . It is also useful to express the triplet-exciton transition rate in a basis-independent form. In contrast to the formation of a singlet exciton, the final state is now also described by a matrix-valued distribution function. The resulting transition rate takes the form where f eh as defined in Eq. (113) appears. For the reverse process t → e + h, the matrix occurs instead, which results in a similar structure. In particular, the basis-independent transition rate from the point of view of the triplet excitons reads as from which the rates for the number, spin, and quadrupole densities can be inferred.

Point of view of electron and hole
In the following, we discuss the process e + h → t from the point of view of the electrons; the results for the holes are analogous. The transition rates of the different spin orientations for z -polarized electrons, can be combined to obtain the transition rates for the electron number and spin densities in analogy to Eqs.
(75) and (76). The resulting expressions are not symmetric in the coefficients of f e and f h for the simple reason that we are writing down transition rates for the electrons, not for the holes. The rates are not sufficient to read off the SU(2)-invariant form, as we will see shortly.
The basis-independent form of the transition rate is particularly useful for this reason. In analogy to the previous cases, the transition rate in terms of density matrices can be written as The notation (•) 4 means that a 3 × 3 matrix on the triplet subspace is extended to the full twoparticle space by adding a 1 × 1 null matrix on the singlet subspace. The matrix is then transformed into the product basis of electron and hole spins, The partial trace over the hole sector is defined by The resulting SU(2)-invariant transition rates given in the multipole expansion are thus for the electron number density and I ce e+h→t = − W e+h→t 12 9c e n h + 3n e c h + 3 √ 6 c e n h n t + √ 6 n e c h n t + 6n e n h c t + 6c e (c h · c t ) + 6n e q t c h (131) for the spin density. This rate contains a product of three spin densities. The corresponding term is proportional to c z e c z h c z t in the limit of z -polarized spins. It is thus clear that this limit does not allow to infer the full SU(2)-invariant expression (although one might guess correctly) so that the basis-independent expression (128) is required.
The results for the holes are analogous. We thus find that the rates satisfy the consistency relations Note that the factors √ 6 and 4 in these two equations stem from the different normalization conditions concerning electrons and triplet excitons. The reverse process can be treated similarly.

D. Fusion and fission of singlet excitons
In preparation for the important transitions between a singlet exciton and two triplet excitons, s ↔ 2t, we first consider the simpler process s ↔ 2s. In contrast to the processes investigated before, the fusion and fission of singlet excitons involves the same particle species in the initial and final states. Furthermore, the two singlet excitons on the same side of the reaction have to be treated as indistinguishable bosons.
Starting with the fusion, this process can either create or fuse a singlet exciton with the outer momentum k. This leads to in-scattering and out-scattering terms in the transition rate, which can be written as (134) Here the short-hand notations W ≡ W (k , k) and W ≡ W (k , k − k ) ≡ W (k , k ) are used. The factor 2 in the out-scattering term results from the two equal contributions of indistinguishable singlet excitons in the initial state. Conversely, the fission of one singlet exciton into two singlet excitons is described by Note that due to the spin-singlet character of all involved quasiparticles, SU(2) symmetry is trivially satisfied.

E. Absorption of a singlet exciton by a triplet exciton
In contrast to the previous process, the particles in the initial state are distinguishable. In the limit of z polarization, this yields, from the point of view of the singlet exciton, In terms of the multipole expansion this becomes The singlet exciton is described by a scalar distribution function so that the basis-independent transition rate is The transition term from the perspective of the triplet exciton features in-scattering and out-scattering terms as the triplet exciton is present on both sides of the reaction. The transition rates for the different spin orientations are Since the (positive) in-scattering and the (negative) outscattering contributions have the same structure it is sufficient to investigate one type. For the out-scattering contributions to the SU(2)-invariant transition rates of the multipole coefficients we obtain where the quadrupole transition rate is again symmetric and traceless. In a basis-independent form, these relations can be summarized as The treatment of singlet-exciton emission is analogous.

F. Fusion and fission involving two triplet excitons and one singlet exciton
The important fusion of two triplet excitons forming a singlet exciton [40] is similar to the fusion of electron and hole into a singlet exciton. In contrast to the latter process, the two particles forming the singlet exciton are indistinguishable bosons.

Point of view of the singlet exciton
Considering the Clebsch-Gordan coefficients, the transition rate in terms of diagonal elements of the distribution function is (146) Written in the multipole expansion and generalized to the SU(2)-invariant form this yields where the relation Tr q t q t = 2 q t · q t has been used. The only differences compared to the creation of singlet excitons from free charge carriers are the prefactor and the additional quadrupole term. The basis-independent form is which is analogous to the basis-independent form in the process e + h → s; see Eq. (86). Here P s = |0 2t 0 2t | projects the two-triplet-exciton state onto the singlet sector, where is the singlet state of the two triplet excitons. The rate for the reverse fission process can be obtained analogously or be read off from Eq. (148),

Point of view of the triplet exciton
For the diagonal elements of the distribution function, we obtain the three transition rates where the factor 2 is due to the indistinguishable triplet excitons in the initial state. The transition rate referring to the triplet-exciton number density is (154) Since the actual density is √ 6 n t , where n t is the expansion coefficient in Eq. (48), the transition rate for the actual density is twice the negative transition rate for the singlet-exciton density in Eq. (147). For the spin transition rate we obtain and the traceless and symmetric quadrupole transition rate is The basis-independent transition rate then assumes the form where Tr t is the partial trace over the sector of the triplet exciton with momentum k . The reverse fission process is then described by

G. Absorption of a triplet exciton by an electron
The process e + t ↔ e is more complicated than the previous cases because both particle species carry spin, in contrast to, e.g., t + s ↔ t. For the process e + t ↔ e to be possible, the total spin of the electron and triplet exciton on the left-hand side must be 1/2, i.e., they must form a doublet. The two states are These relations allow us to read off the numerical prefactors in the rates.

Point of view of the triplet excitons
In the limit of z polarization, the three transition rates describing triplet excitons are The SU(2)-invariant transition rates in terms of multipole coefficients are obtained in analogy to the previous cases. For the triplet-exciton number density we find for the spin density and for the symmetric traceless quadrupole density (c e · c e ) 1 3 n t Like for the process e + h → t, one might guess at this SU(2)-invariant form based on the limit of z -polarized spins but cannot infer it rigorously. The problem lies in the terms containing three c vectors and in the terms with two c vectors and the triplet-exciton quadrupole tensor q t . The above expressions have thus been checked by comparing them with the basis-independent expression describing this process, which takes the form where P d is the projection operator onto the doublet subspace. This is a 6 × 6 matrix. The notation (•) 6 means that a 2 × 2 matrix on the doublet subspace is extended to the full two-particle space by adding a 4 × 4 null matrix on the quartet subspace. The matrix is then transformed into the product basis for electron and tripletexciton spin to be consistent with f e ⊗ f t .

Point of view of the electron
The perspective of the electron leads to in-scattering and out-scattering terms. The resulting transition rates for the two spin orientations are where the Clebsch-Gordan coefficients have been used. As the in-scattering and out-scattering terms have the same structure it is sufficient to investigate out-scattering. For the multipole transition rates this yields 6 n e n t n e − 6n e c e · c t + √ 6 n t c e · c e + 6n e c e · c t + 6c e · q t c e , I ct, out e+t→e = − W e+t→e 18 12n e c t − 6 √ 6 c e n t − √ 6 n e n t c e − 6n e c t n e + 3 √ 6 c e n t n e + 6c e (c t · c e ) − 3 2 n e q t c e , which can be combined to obtain the full transition rate from the point of view of the electron. Together with the in-scattering terms the basis-independent form is then The absorption and emission of a triplet exciton by a hole is analogous.

H. Fusion of triplet excitons into a triplet exciton
The fusion of two triplet excitons into a triplet exciton [40] is the most complex process considered here. In the case of z polarization, the transition terms of the components of the density matrix read as In this process, the in-scattering and out-scattering terms differ in structure. This makes it necessary to investigate both terms separately.
1. In-scattering for 2t → t Combining Eqs. (173)-(175), we can obtain the transition rates expressed in the multipole expansion. Their generalization to the SU(2)-invariant form is ambiguous at first glance. However, the SU(2)-invariant form can be obtained in analogy to the formation of triplet excitons out of free charge carriers. We construct the density matrix f t ⊗ f t on the product space of the two tripletexciton spins and project it onto the triplet subspace. The resulting 3 × 3 matrix, which will be referred to as f tt , can be expressed in terms of coefficients of the nine basis matrices S i , i = 1, 2, 3, and Q j , j = 1, 2, 3, 4, 5. For a spin along the z axis, only three coefficients do not vanish: Combining these equations and using Eqs. (173)-(175), we can unambiguously generalize the rates to their SU(2)-invariant form: The basis-independent in-scattering term reads as . (182) The notation (•) 9 is analogous to (•) 4 in Sec. III C 2 and to (•) 6 in Sec. III G 1.

Out-scattering for 2t → t
For out-scattering, the strategy employed for inscattering does not work since the rates for the case of z polarization are not symmetric in the two incoming triplet excitons at the outer momentum k and the running momentum k . Hence, the projected density matrix on the product space, f tt , is not useful. The situation is analogous to the process e + h → t from the point of view of the electron or hole; see Sec. III C 2. Like for e+h → t, the matrix-valued transition rate in terms of projection operators and density matrices can be used to read off the correct SU(2)-invariant multipole transition rates. For the triplet exciton number density this yields Tr q t q t + 4n t n t n t + 2 n t c t · c t + n t c t · c t + n t c t · c t − n t Tr q t q t + n t Tr q t q t + n t Tr q t q t The spin transition rate becomes The quadrupole part reads as 2n t Tr q t q t − n t Tr q t q t − n t Tr q t q t 1 3 The derivation of the rates for the spin and for the quadrupole moment again require the basis-independent form in Eq. (183).
The reverse fission process t → 2t can be investigated in the same way. In the out-scattering term, the intro-duction of f 1−tt leads to the same structure as in the in-scattering term for 2t → t.

IV. EXAMPLES
In this section, we illustrate the theoretical framework by applying it to a simple model. In order to reduce numerical complications, we consider a one-dimensional toy model containing fission and fusion processes, disorder scattering of singlet and triplet excitons, as well as relaxation of singlet and triplet excitons. All bare rates W are assumed to be independent of momentum, except as dictated by scattering being elastic. The temperature is assumed to be sufficiently low so that the equilibrium thermal population of excitonic states is negligible. A magnetic field, if present, is taken to be uniform, constant in time, and oriented along the z direction.
The dispersion relations of singlet and triplet excitons are written as respectively, where m = 1, 0, −1 specifies the spin state of triplet excitons and is the energy difference between singlet and triplet excitons, which defines a momentum scale k 0 .
Choosing our units such that /m t = 1 and defining the parameters the Boltzmann equation for the singlet-exciton distribu-tion function f s (x, k, t) can be written as and the corresponding equation for the matrix-valued triplet-exciton distribution function f t (x, k, t) as Here I fs and I ft are the full scattering integrals, which contain the following contributions, where we only make the momentum dependence explicit: (a) Decay is discussed in Sec. III A. It is described by the scattering integrals (b) Disorder scattering is discussion in the introduction to Sec. III and is described by Note that the singlet-triplet energy difference ∆ and the magnetic energy of triplet excitons cancel. (c) Fusion is discussion in Sec. III F. From the point of view of the singlet excitons, fusion is described by with with κ given by Eq. (203), if κ is real and zero otherwise.
One can now expand f t into multipoles using Eq. (48). This is particularly useful if some of the multipoles vanish by symmetry. In any case, the equations to be solved are coupled first-order nonlinear partial differential equations for f s and f t as functions of time t and position x. The momentum k acts as a parameter labeling different coupled functions. None of the coefficients depend on x or t. Nevertheless, due to the nonlinearity of the equations, an analytical solution seems out of reach. We therefore use numerical forward propagation of the initial functions with discrete time steps and also discretize the x and k dependence.
As an example, we consider a localized distribution of singlet excitons, photogenerated close to the energy threshold. The distribution is modeled by a Gaussian function with mean zero and width σ x in real space multiplied by a Gaussian function with mean zero and width σ k in momentum space. To reflect some of the features of real systems, the singlet excitons are assumed to decay rapidly. We consider the singlet-excition distribution function in real space, and the distribution function of the occupation of triplet excitons in real space, For reference, we plot the time evolution of the singletexciton distribution function ρ s in Fig. 1 for the situation without fission. In this case, the singlets decay rapidly and triplet excitons are not generated. Note that the densities ρ s and ρ t in all plots are presented in arbitrary units but on the same color scale. If fission (but not fusion) is switched on the singlet excitons quickly transition into triplet excitatons, as shown in Fig. 2(a). The triplet excitons show predominant diffusive motion due to disorder scattering, see Fig. 2(b). The cone from ballistically moving, i.e., not scattered, excitons is visible as a weaker feature.
If we now also switch on fusion with a high bare rate we observe recovery of singlet excitons and loss of triplet excitons, as expected, see Fig. 3. However, the distribution of recreated singlet excitons is much narrower than the one of the triplet excitons. We attribute this to the fact that for our parameters two triplet excitons must collide head on in order to fuse, which is more likely in the center of the triplet cloud. It is an interesting question for the future to what extent such an effect is also present in real three-dimensional systems.
In the previous examples, the spin of the triplet excitons does not matter. The initial singlet distribution is of course nonmagnetic and there is no mechanism that generates any spin polarization. To exhibit a characteristic magnetic effect, we consider an initial state with a distribution of triplet excitons with momenta close to ±k 0 (a) (b) Figure 4. Spatial distribution function of (a) the occupation number of triplet excitons and (b) the x component of their spin polarization as functions of time in an applied magnetic field in the z direction. The initial distribution of triplet excitons has helical spin polarization with the spin parallel (antiparallel) to x for k > 0 (k < 0). The dimensionsless Larmor frequency is ωL = 10, the other parameters are the same as in Fig. 3. and helical spin polarization: the spin is polarized parallel (antiparallel) to x for k > 0 (k < 0). Except for the spin polarization, this distribution is similar to the one obtained by singlet fission in the previous examples. The occupation number ρ t of triplet excitons, shown in Fig.  4(a), is independent of the applied magnetic field. The x component of the spin polarization given by and plotted in Fig. 4(b) shows two effects: First, due to the helical spin polarization of the initial distribution, triplet excitons moving to the right (left) have positive (negative) spin. Second, their spin precesses about the magnetic field, i.e., the z direction, with a frequency given by the Larmor frequency ω L in Eq. (193).

V. SUMMARY AND CONCLUSIONS
The understanding of transport of charge carriers and excitons in organic semiconductors is crucial for applications. Semiclassical Boltzmann theory is a valuable tool for this since it is formulated in terms of physically intuitive quantities and processes, thereby simplifying the interpretation. In this paper, we have addressed the spinconserving charge transport and scattering in a semiclassical framework, including the creation, decay, recombination, fission, and fusion of excitons, which are crucial for organic solar cells. However, our results are of more general interest for any system showing scattering of and reactions between spin-carrying particles, for example in spintronics and spin chemistry.
To start with, we have set up kinetic equations for charge carriers as well as for singlet and triplet excitons in uniform electric and magnetic fields, taking care to preserve SU(2) spin-rotation invariance. The case of triplet excitons is the most interesting in that these spin-1 quasiparticles permit and require a quadrupolar contribution to their distribution function. We have obtained coupled equations for the concentration, spin-density, and quadrupole-density components. The generalization to nonuniform systems is conceptually straightforward.
In the second step, we have constructed SU(2)-invariant transition terms or scattering integrals describing disorder scattering, singlet-exciton generation, binding of electrons and holes into singlet and into triplet excitons, and exciton fusion, as well as the reverse processes. The occupation of final states enhances a process for bosons and suppresses it for fermions. These quantum effects are naturally included and lead to higher-order terms in the scattering integrals, which are relevant for large occupation numbers. It proved useful to employ two different representations of distribution functions, namely a multipole expansion and a description as matrix-valued functions acting on the spin Hilbert space. For the simpler processes, the multipole representation allows one to directly write down the SU(2)-invariant terms based on the transition terms under the (symmetry-breaking) assumption of all moments being parallel or antiparallel to the z axis. The results can be used to describe real materials when realistic models for the bare transition rates W are available. Moreover, our work serves as a guide on how to construct scattering integrals for complex transitions involving several quasiparticles with spin.
Future directions following from this work are quite clear. First, it will be necessary to derive models for the bare transition rates W for real materials and link the equations to observables, in order to obtain quantitative descriptions. For the application to photovoltaic devices, it is important to note that the generation of excitons takes place at internal interfaces in organic blends, while the diffusion of (triplet) excitons is a bulk effect. Our framework offers two approaches to this problem: On long length scales relative to the typical scale of the phases in the blend, coarse graining leads to an effective-medium description with an effective exciton-generation rate. On the other hand, on shorter length scales-but still long compared to the single-molecule scale-the spatial structure of the phases can be modeled. The semiclassical approach is in principle suitable for this but of course relies on a good structural model.
Second, additional processes are relevant in certain materials. For example, bound states of two electrons and one hole or vice versa, i.e., trions, can exist [72,73]. The formation and decay of trions can easily be described in our framework. Third, we have here assumed that polaronic effects, which are strong in organic materials, can be treated simply be renormalizing model parameters for charge carriers and excitons. This is certainly simplistic, in particular due to the broad distribution of timescales of phonons in organic materials. The connection of transport and excitonic processes with polaronic effects in a semiclassical framework is thus an important goal for the future.