Generalized Kinetic Monte Carlo Framework for Organic Electronics

Abstract

In this paper, we present our generalized kinetic Monte Carlo (kMC) framework for the simulation of organic semiconductors and electronic devices such as solar cells (OSCs) and light-emitting diodes (OLEDs). Our model generalizes the geometrical representation of the multifaceted properties of the organic material by the use of a non-cubic, generalized Voronoi tessellation and a model that connects sites to polymer chains. Herewith, we obtain a realistic model for both amorphous and crystalline domains of small molecules and polymers. Furthermore, we generalize the excitonic processes and include triplet exciton dynamics, which allows an enhanced investigation of OSCs and OLEDs. We outline the developed methods of our generalized kMC framework and give two exemplary studies of electrical and optical properties inside an organic semiconductor.

1. Introduction

Organic semiconductors represent a promising candidate for cost-efficient, environmentally-friendly materials and renewable energy conversion devices. The advantage compared to currently widely-used inorganic semiconductor-based devices is captured in their mechanical flexibility, light weight and electrical and optical controllability at the molecular level. Organic semiconductors have exhibited a considerable potential for application in organic thin-film transistors (OTFT) [1], light-emitting diodes (OLEDs) [2] and solar cells (OSCs) [3,4]. Experimental research on organic-based devices is accompanied by computational methods providing a better understanding of the underlying processes in energy conversion devices made from nanostructured organic materials.

In the simulation of organic materials and devices, several methods on different length scales are of interest. While drift-diffusion (DD) simulations yield a good understanding at a macroscopic continuum level of several $\mu \mathrm{m}$, they can neither provide insight into the microscopic, internal properties, nor track the particle dynamics. Tools solving molecular dynamics give detailed insight into the processes occurring on an atomistic length scale ($\le 1\mathrm{nm}$), but the computational demand is not feasible for the simulation of bulk materials or thin films of a $100\mathrm{nm}$ thickness. Kinetic Monte Carlo (kMC) simulations have gained more and more attention in the field of materials science in recent years. Due to strong localization of charge carriers on molecules or molecular subunits, kMC is highly feasible for the modeling of organic disordered semiconductors. The aim of understanding the internal processes in organic devices and the rising available computational power has driven the development of the above microscopic numerical methods. Existing kMC models can reproduce the main macroscopic characteristics of organic optoelectronic devices and provide insight into fundamental physical properties considering charge transport and optical processes. Thus, kMC provides a bridge between different length scales from the atomistic level up to a macroscopic device level.

Ries and Bässler have presented one of the first applications of kMC in the study of charge transport in organic materials [5]. From that point on, kMC has been employed for many theoretical studies covering charge and exciton transport [6,7,8], recombination of charge carriers [9,10] and the impact of all processes on the photocurrent generation [11,12,13,14,15]. kMC simulations have been frequently utilized to study full devices such as bulk-heterojunction (BHJ) organic solar cells [10,11,12,13,15,16]. There are several groups working on kMC models describing all properties mentioned above to a certain extent. Many of them consist of a cubic lattice model, which cannot represent all geometrical features and the impact of the structure of organic materials on the physical processes. Furthermore, many do not consider the directionality of charge transport along, e.g., polymer chains. Another limitation of existing models is the simplified treatment of excitonic processes, which are of high importance regarding organic devices. Most kMC models either account for singlet [10,17,18,19] or triplet exciton dynamics [20,21,22,23] for the use in OSCs or OLEDs, respectively. However, the interplay between singlets and triplets can be used to modify exciton lifetimes and diffusion lengths, which is a promising way to tailor efficiencies in organic devices and therefore requires a better understanding [24,25,26].

In this work, we present our in-house-developed generalized kMC framework for the simulation of organic semiconductors and devices. Our model overcomes the mentioned limitations by the inclusion of a non-cubic lattice model, based on the Voronoi tessellation, and a model of polymer chains. Both yield a better geometrical description of organic materials and herewith provide new insight into the electrical properties. Charge transport is a strongly percolative process [27,28]. Thus, a generalized lattice model is of high importance in the kMC method. Considering the excitonic properties, we present a detailed model of exciton dynamics including both singlet and triplet excitons and their interaction. We present two exemplary studies on electrical and optical properties of organic materials using results obtained with our generalized kMC model.

In Section 2, we give an overview of the kMC method and show its applicability for the simulation of organic materials. A detailed description of the non-cubic Voronoi set, polymer chains, the model of the energetic landscape and excitonic processes is presented in Section 3. Section 4 discusses two exemplary studies. Firstly, the charge transport properties of organic materials with the inclusion of polymer chains is presented. We investigate the impact of the directionality of hopping processes on the mobility. Secondly, the effective lifetime of excitons inside an organic material is studied. We focus here on the impact of the dopant volume concentration on the effective exciton lifetime.

2. Kinetic Monte Carlo Method

Monte Carlo (MC) simulations refer to a broad class of numerical algorithms utilizing random numbers to solve real-world problems. The first Monte Carlo methods were developed during the 1940s within the Manhattan project to study neutron diffusion in fissionable materials. During that time, the Metropolis algorithm [29] originated, which is until today a convenient technique to study equilibrium properties of physical systems. However, the Metropolis algorithm is not capable of modeling the time-dependent evolution of a system. The kinetic Monte Carlo method is a specific MC method that allows propagating a physical system from state to state including the corresponding time information of the propagation steps. Therefore, kMC simulations are often termed simulated experiments because they essentially mimic the explicit evolution of a real-world system. The step by step evolution, where each time-step depends only on the previous one, fulfills the properties of a Markov chain and makes the kMC algorithm a sequential algorithm. kMC methods have gained large popularity in a variety of fields due to their ease of implementation and the increasing computational power of modern computers. The modern version of the kMC method is based on the work of Gillespie [30] and Bortz et al. [31]. For further information for the reader, an illustrative description of the method was provided by Voter [32].

To model a system, one has to identify the spatial dimensions and the time scales on which the processes determining the properties of the system occur. The kMC method commonly is referred to as a simulation technique at the mesoscale, in between molecular dynamics (MD) on the microscopic and hydrodynamic models on the macroscopic scale. In comparison to MD, which typically resolves systems on the timescale on which molecular vibrations and rotations occur (of the order of femtoseconds), the kMC method is based on the concept of coarse-graining the system dynamics. A system is characterized by a set of long-term states with discrete transitions between them (so-called infrequent events); all transitions on smaller timescales are neglected as long as they do not change a long-term state (illustrated in Figure 1a). Coarse-graining allows one to neglect a sufficient amount of (fast) processes and, thus, reach much longer overall simulation times up to seconds. Only transitions between long-term states determine the dynamic evolution of the system. It needs to be carefully considered which processes are implemented in the model and which ones are combined in long-term states. This way, the kMC can accomplish a tradeoff between accuracy and simulation time.

For the application of the kMC method, the system to be modeled needs to be characterized by a set of long-time system states $\left\{i\right\}$. These can be, for example, the entirety of positions and energies of all particles currently in the system. For a system being in a certain state i, a set of possible states j exists to which a system may evolve (illustrated in Figure 1b). Every transition is accompanied by a transition rate $a_{ij}$ (given in units of $s^{-1}$), describing the physical process leading to the transition $i\to j$. The magnitude of the rates is proportional to their probability of execution, i.e., larger rates (faster processes) are more likely to occur than smaller rates (slower processes).

The temporal evolution of discretized transitions between states is phenomenologically described by the master equation [33]:

$$\frac{\mathrm{d}{p}_i\left(t\right)}{\mathrm{d}t}=-\sum _j{a}_{ij}{p}_i\left(t\right)+\sum _j{a}_{ji}{p}_j\left(t\right),$$

(1)

where $p_{i/j}\left(t\right)$ gives the probability to find the system in state $i/j$ at time t. An analytic or numerical solution of the master equation is often not feasible. The kMC method provides a numerical approach to essentially execute the master equation to stochastically propagate the system between its long-term states. If the coarse-graining is chosen accurately such that the correlation time of the environment in the continuous-time dynamics is much shorter than the typical rate constants, the discretized system evolution $i_1\to i_2\to \cdots \to i_N$ can be regarded as Markovian. A single trajectory of state transitions does not yield the same system behavior as described by the master equation. Averaging over several trajectories provides an equivalent system behavior one would obtain by solving the master equation.

The dynamical evolution of the system is governed by the transitions rates. The rates are inputs to the kMC model and can either come from experimental measurements or are based on underlying theoretical models. Given rates, the kMC algorithm is straightforward and can be reduced to the following scheme. For a detailed description of the derivation, we refer to [30,34].

• Initialize the system, and define all parameters that are required for the transition rates of all processes in the kMC. Define termination conditions such as the maximum simulation time $t_{\max }$.
• Identify the activated transitions for the current system state i and calculate the transition rates $\left\{a_{ij}\right\}$ of all physical processes that can occur. The rates $a_{ij}$ mirror the probability of the occurrence of the transition to state j. The sum $a={\sum }_j{a}_{ij}$ of all rates is calculated.
• Two uniformly-distributed random numbers $r_1\in \left(0,a\right)$ and $r_2\in \left(0,1\right)$ are chosen. The first number $r_1$ chooses the transition $\mu$ from the set of events $\left\{a_{ij}\right\}$, which fulfills the condition:

  $$\sum _{j=1}^{\mu -1}{a}_{ij}<r_1\le \sum _{j=1}^{\mu }{a}_{ij}.$$

  (2)
  The second number $r_2$ gives the time passed before the transition is performed, given by:

  $$\tau =-\frac{ln\left(r_2\right)}{a}.$$

  (3)
  The derivation of Equations (2) and (3) is based on the direct method [30], which divides the probability density function (PDF) for the transition $\mu$ and the corresponding time step $\tau$ into two separate PDFs, each of which only depends on either the time $\tau$ or the transition $\mu$.
• The simulation time t is advanced by $\tau$, and the system configuration is updated following the transition $\mu$.
• If the simulation time exceeds $t_{\max }$, the simulation stops, and all particle trajectories and the counters tracking the events are evaluated. If $t<t_{\max }$ , the kMC loop (Step 2–4) is repeated.

For simple systems such as unipolar single layer devices, the master equation has been applied successfully [35,36,37]. However, critical properties such as the electrostatic Coulomb interaction between charge carriers can only be captured by an effective potential in the framework of the master equation. With the particular application to OSCs in mind, particles like excitons, electrons and holes reside in localized states of either monomers of a polymer or small organic molecules. Charge motion can occur by processes on shorter timescales, but as long as they do not perform a hopping process to another localized state, the particles are spatially assigned to that state. The overall dynamical behavior is determined by hopping events between localized states. This concept has been applied to model charge transport in disordered organic semiconductors, as well as organic devices such as OSCs or OLEDs. A significant contribution of describing charge transport has been conducted by Bässler [6]. Several models for OSCs have been set up investigating active layer morphology and exciton separation [11,17], recombination [10,14,38], Coulomb interaction [12] and charge densities [13].

3. kMC Framework for Organic Materials and Devices

In this section, we focus on the new features of our kMC framework extending the previously shown model [10,13]. A generalized representation of the organic material is presented, which covers the model of a non-cubic Voronoi set and polymer chains. Furthermore, dopants are included, which need to be taken into account for the simulation of all excitonic processes in organic devices. Correlations of the energy landscape are of importance for many materials. We present a simple weighted average approach in Section 3.2. The Voronoi tessellation and the correlation procedure have been discussed recently in more detail [39]. As they are part of our generalized kMC model, we provide a concise summary. We conclude this section with a detailed description of the particles and transitions containing the information of the physical properties.

3.1. Modeling of the Organic Material

3.1.1. Voronoi Tessellation

With the rising amount of computational power and powerful clusters, investigations on the detailed structure of organic materials and their influence on the physical properties using molecular dynamics (MD) have gained rising interest in the community [40,41,42]. This raises the demand of a generalized lattice representation for kMC studies [6,43,44,45].

For a representative model of the active organic layer, the use of the Voronoi tessellation generalizes the cubic lattice model [39,46]. Using the Voronoi cell library of Rycroft [47], next-neighbors of each lattice site, as well as the distance and all face-areas are computed. In Figure 2, the procedure is shown for the cubic lattice (left) and the non-cubic, Voronoi set (right) in a simple 2D representation. Starting from an arbitrary grid point, the perpendicular bisectors to all neighboring sites in proximity are computed. The algorithm obtains the volume closest to the considered starting point, named the Voronoi cell, as known from the Wigner–Seitz cell in periodic crystals. All neighboring sites sharing one face with the considered Voronoi cell are taken as next-neighbors. The knowledge of these is fundamental for the computation of all physical rates.

Comparing the cubic with the off-lattice approach, some obvious differences and their importance are to be highlighted. Due to the spatial disorder of the sites, the number of next-neighbors increases. Moreover, the major drawback of an equidistant spacing between lattice points in the cubic model is overcome with the use of the Voronoi set [48]. The local distribution of lattice sites can be imported from MD simulations using the center of masses of the molecules, providing a sophisticated multi-scale simulation method comparable to [42,43]. In the case of available quantum chemical simulations, the localized quantum state could directly be used as the input of the sites. Herewith, the information depth of the structural alignment of, e.g., polymers or small molecules down to the nanoscale can be represented on a mesoscopic level. Besides, amorphous domains in both polymers and small molecules can be explicitly investigated, which makes the study of the physical properties more realistic. In the case of the exciton diffusion (see Section 3.3.2), non-nearest neighbors are required. The Voronoi tessellation supports a fast iteration method to obtain all neighbors within a certain radius around the centroid.

Existing off-lattice models take into account all adjacent molecules within a certain radius [43,49]. Depending on the density of sites, this may cause a strong increase in simulation time. Using the Voronoi tessellation, only adjacent sites are considered, which keeps the computation time feasibly low. In the case of a strong disorder, a common radius might lead to the case that certain sites are fully disconnected from the surrounding material. This is avoided with the Voronoi tessellation. The tessellation furthermore allows the computation of local current and charge densities as a well-defined area of the faces, and the volume of the Voronoi cells is obtained [39]. Besides all the advantages, it still needs to be considered that the cell boundaries are purely based on the distance between centroids and do not contain any physical information, which, depending on the modeled system, has to be evaluated carefully for the computation of local quantities.

3.1.2. Polymer Chains

Organic materials used in organic devices can be distinguished into two separate categories, small molecules and polymers. Polymers are composed of repeating macromolecules called monomers forming long chains. One prominent example for polymers used in OSCs is poly(3-hexylthiophene) (P3HT) [3]. In this section, we briefly emphasize the importance of the correct model of polymer chains in our kMC model and present the algorithm for the chain morphology generation.

The spatial overlap of molecular orbitals leads to weaker bound and less localized electron states along polymer chains [50]. Therefore, charge transport along polymer chains, called intrachain transport, is faster than charge transport between two separate polymer chains, which we refer to as inter-chain transport [51,52,53]. Vanlaeke et al. showed that crystallization of P3HT increases the mobility [54]. Sirringhaus et al. have demonstrated that the mobility strongly depends on the orientation of the micro-crystalline domains in semi-crystalline P3HT [55,56]. In situ grazing incidence X-ray scattering techniques can be exploited to gain information on the polymer crystallization and ordering [57].

Many present kMC models treat every site as an independent localized state, and charge transport to all of its neighbors is described with the same mechanism [6,58]. Multi-scale simulation models combining molecular dynamics and kMC include the distance dependence and directionality captured in the coupling integrals calculated from MD simulations, but did not investigate the effect of polymer chains on mobility [42,43]. Mollinger et al. model single charge transport between crystalline regions based on kMC along a generated polymer chain and investigate the impact of the electric field on the single charge mobility [44]. However, the probability of leaving the chain is modeled with an effective, field-dependent exit rate and does not account for neighboring chains. To account for both inter-chain and intra-chain transport, it is necessary to generate a morphology in which all the sites belonging to a polymer chain are connected.

The algorithm for the chain morphology generation is based on a self-avoiding random walk [59]. It starts at a random free site, a node that is not yet part of a polymer chain. Then, all free neighbors of the site are determined. Free neighbors are next-neighbors, which are not already part of another polymer chain. One of these free neighbors is randomly chosen and connected to the chain. Every node stores a link to the previous and next chain element. From this newly-added node, again, the free neighbors are determined, and one is added to the chain. This procedure is repeated until the desired chain length is reached or there are no longer any free neighbors. Then, a new free site is searched, and the next polymer chain is generated until the whole simulation box is filled with chains. In our model, we assume that the chain length in organic materials follows a Gaussian distribution. Therefore, the average chain length $l_c$ and its standard deviation $\sigma$ are necessary input parameters for the algorithm. The desired chain length for every single chain will then be calculated according to the Gaussian distribution function. Chain length and distribution can be measured experimentally using, e.g., gel permeation chromatography [60] or static light scattering [61].

According to Carbone and Troisi [51], the orientation of the polymer chains is of great importance for the overall transport characteristics of organic materials. Therefore, we account for a preferred direction and the amount of orientation of the polymer chains, including three distinguished chain models. Gaussian chains are flexible and show no preferred orientation. They are obtained by the random walk algorithm described above. Rigid chains are entirely directed polymers pointing in one direction. The crossover between these two modes is described by the worm-like chain model, where a local stiffness is added to the polymer yielding to semi-flexible chains with a preferred direction. The parameter $\eta \in \left[0,1\right]$ controls the amount of stiffness, where a value of zero results in flexible Gaussian chains and a value of one leads to rigid chains where the polymers can only be generated along the preferred axis. For values in between, the probability for connecting sites to the chain along the preferred axis $p_{+}$ is increased, while in-plain probability $p_{\parallel }$ or backward growth $p_{-}$ is reduced. We use the equations:

$$\begin{array}{cc}\hfill p_{+}& =1+\eta ,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill p_{\parallel }& =1-\eta ,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill p_{-}& ={\left(\frac{1-\eta }{1+\eta }\right)}^2\hfill \end{array}$$

(6)

to account for these probabilities. Figure 3 shows the morphologies generated with our algorithm. The two single material morphologies show the Gaussian chain model ($\eta =0$) on the left and the worm-like chain model with high directionality of $\eta =0.9$ in the middle. The 2D case with regular grid is chosen for a better visibility, while the algorithm can be used to generate chains in 3D. For an irregular Voronoi set, the distance vector to the neighboring site is decomposed into in-plane and parallel components. The probability for adding a site to the polymer is obtained by averaging its components calculated with Equations (4)–(6).

It is also possible to generate polymer chains in BHJs. The desired chain length and orientation can be specified for each material individually. In the right part of Figure 3, we visualize a blend of one polymeric material (e.g., P3HT) and one small molecule material (e.g., PCBM). This tool for morphology generation leads to a more realistic representation of a polymeric material and herewith provides better insight into charge transport properties in organic materials and devices on a mesoscopic scale.

3.1.3. Phosphorescent Dopants

OSCs are one representative of the class of excitonic devices. In contrast to their inorganic counterparts, an absorbed photon does not directly excite a free electron-hole pair, but forms an exciton with binding energies exceeding $500\mathrm{meV}$. A morphology that introduces a high exciton splitting efficiency is the bulk-heterojunction (BHJ) [3,62]. The BHJ reduces the distance between the generation position of an exciton and the donor-acceptor interface. This distance should be as short as possible, as singlets have short diffusion lengths of 3–6 $\mathrm{nm}$ [63]. If the distance between origin and interface exceeds the diffusion length, singlet excitons probably decay before reaching an interface and cannot be converted into usable electrical energy.

Another possibility of enhancing the splitting efficiency is to enhance the exciton lifetime and diffusion length. Triplet excitons exhibit an effective lifetime of $10\mu \mathrm{s}$ for P3HT, which corresponds to a diffusion length of $100\mathrm{nm}$ [63]. The spin selection rules forbid the radiative decay of triplet excitons. However, due to the spin-orbit coupling triplets may radiatively decay, but this is very improbable for commonly-used organic materials. Triplet excitons cannot be generated by the photon absorption as it is forbidden by spin-selection rules [64]. In the presence of strong spin-orbit coupling, the spin-selection rules are lifted. Thus, singlet excitons can be converted into triplet excitons. This process is known as intersystem crossing (ISC) [64].

Incorporating heavy metal atoms into the organic layer induces a strong spin-orbit coupling. Doping of a BHJ with nanoparticles can increase the effective exciton lifetime and thus improve the power conversion efficiency [24,25,65,66]. Furthermore, in phosphorescent OLEDs, organometallic complexes are used to improve the efficiency. The electrically-generated triplet excitons decay radiatively and help thus to enhance the efficiency of the OLED [67,68].

Dopant sites are implemented to investigate the impact of doping with strong spin-orbit materials in OSCs or OLEDs. Dopant sites are assigned randomly to the lattice sites in the simulation box. The ratio of dopant to organic sites gives the volume concentration of the dopant. The spin-orbit coupling of an atom affects excitons only in its proximity. Thus, a radius is introduced, in which the spin-orbit coupling influences the excitons. As explained in Section 3.3, the intersystem crossing is activated within this radius.

3.2. Energetic Landscape

The energetic landscape of an organic semiconductor plays a significant role in the charge transport and exciton dissociation [10,27,28]. In contrast to inorganic semiconductors, the molecular orbital energies inside an organic material are strongly disordered due to the non-crystalline arrangement of the molecules and polymers [28]. Energetic disorder values of $\sigma =0.1\mathrm{eV}$ are present in many materials [6]. It is widely accepted that the shape of the density of states (DOS) can be modeled using a Gaussian distribution [15,27,36,58,69,70]. Some experimental observations require a correlation in the site energies, as described in the correlated disorder model (CDM) [15,36,71,72], while many studies are based on the uncorrelated, Gaussian disorder proposed by the Gaussian disorder model (GDM) [6,58,72].

We implement a correlated disorder model based on a simple averaging method of site energies. Figure 4 visualizes the CDM for an organic system consisting of PCBM (black dots) and P3HT (red dots). First, for every site, we choose an energy randomly from a Gaussian distribution:

$$g\left(ϵ\right)=\frac{1}{\sigma \sqrt{2\pi }}exp\left(-\frac{{(ϵ-E^0)}^2}{2{\sigma }^2}\right),$$

(7)

with the molecular orbital energy $E^0$ and the global energetic disorder $\sigma$. Herewith, the energy distribution in the GDM is obtained. For the CDM, the following procedure is performed.

• Starting from the randomly Gaussian-distributed energetic landscape, we compute the weighted average of the energies $E_k$ of all sites k with $r_{ik}=\left|{\mathbf{r}}_i-{\mathbf{r}}_k\right|\le r_{\mathrm{corr}}$ for every site i using:

  $$E_i^{*}=\frac{1}{{\sum }_k\frac{1}{\left|r_{ik}\right|}}\sum _k{E}_k\frac{1}{\left|r_{ik}\right|}.$$

  (8)
  Note that the average is only performed over sites k representing the same material as site i. The average energies $E_i^{*}$ are stored in an intermediate vector.
• Replace the site energies by the average energies $E_i\to E_i^{*}$.
• Scaling of the site energies: The averaging procedure reduces the energetic disorder. To keep the initially chosen $\sigma$, we compute the energetic disorder ${\sigma }^{*}$ after the averaging process. Then, the site energies are scaled by:

  $$E_i\to E_i\sqrt{\frac{\sigma }{{\sigma }^{*}}}-\left(\sqrt{\frac{\sigma }{{\sigma }^{*}}}-1\right)E_i^0.$$

  (9)
• The level of correlation can be varied by either changing $r_{\mathrm{corr}}$ or repeating Steps 1–3.

A measure of the level of correlation is defined by the correlation function $c\left(r\right)$. The correlation of sites i and j can be obtained with:

$$c_{i,j}=\frac{\left(E_i-E_i^0\right)\left(E_j-E_j^0\right)}{{\sigma }^2}.$$

(10)

If we take the average of the correlation of the ensemble $〈·〉$ of pairs ($i,j$) with the same distance $r_{ij}$, we obtain a measure of the correlation of the whole organic material. Using this correlation mechanism, we have shown that a high correlation in the site energies gives rise to filamentary structures of the local currents [39]. For distances below $2\mathrm{nm}$, the correlation function obtained by atomistic simulations [73] has been reproduced with our mechanism, while for long-range hopping, other procedures, e.g., a dipole-based correlation [36], give better results. While a dipole-based correlation is obtained by summation over all dipoles within the simulation box [36,71], the presented method is computationally favorable as only sites within a certain radius are considered.

3.3. Particles and Transitions

3.3.1. Charge Carriers

In disordered organic materials, charge carriers are spatially localized due to the lack of long-range order. Charge transport is modeled as a hopping process from one localized state to a neighboring localized state. This state can be seen as local energy valleys where charges locate in a long-term, which e.g., represent a small molecule or a few monomers along a polymer chain.

Miller and Abrahams describe hopping as a combination of thermal activation and tunneling [74]. The hopping rate between states i and j reads:

$$a_{ij}=a_{0}exp(-2\gamma r_{ij})\left\{\begin{array}{cc}exp\left(-\frac{\Delta E_{ij}}{k_{B}T}\right)\hfill & \mathrm{for}\Delta E_{ij}\ge 0,\hfill \\ 1\hfill & \mathrm{for}\Delta E_{ij}\le 0,\hfill \end{array}\right.$$

(11)

where $a_0$ is the attempt-to-escape frequency, $\gamma$ is the inverse localization length, $k_B$ is the Boltzmann constant, T gives the temperature, $r_{ij}$ the distance and $\Delta E_{ij}=E_i-E_j$ the energy difference between the initial configuration i and the final configuration j, before and after a hopping process, respectively. When hopping upwards in energy, the Boltzmann term accounts for thermal activation [74]. This rate model is usually used when the coupling between electrons and phonons is weaker [75,76].

Marcus has introduced a model that accounts for stronger coupling between lattice and charges [77]. Because of the low permittivity of organic materials, charges lead to local polarizations of the lattice, which are described as quasi-particles called polarons. The formula for polaron transitions reads:

$$a_{ij}=\frac{{\left|I_{ij}\right|}^2}{\hslash }\sqrt{\frac{\pi }{\lambda k_{B}T}}exp\left(-\frac{{(\Delta E_{ij}-\lambda )}^2}{4\lambda k_{B}T}\right),$$

(12)

with the reduced Planck’s constant ℏ, transfer integral $I_{ij}$ and the reorganization energy $\lambda$ describing the reconfiguration of the molecule due to the deformation by the charge. The effect of a spatial disorder in an irregular Voronoi set is included via the site to site distance $r_{ij}$ in Miller–Abrahams rates (Equation (11)) and via the distance-dependent transfer integral $I_{ij}$ in Marcus rates (Equation (12)). The hopping rates are calculated for every particle and all free neighbors. If a particle already occupies a neighboring site, the hopping rate is zero. When a hopping transition occurs, the occupation information of origin and destination sites and the particle position have to be updated.

To calculate hopping rates, it is necessary to compute the electric energy $E_i$ for every site i. The energy at site i is given by:

$$E_i=E_i^0+E_i^{\sigma }+E_i^F+E_i^C,$$

(13)

with $E_i^0$ being the energy level of the molecular orbital and the energetic disorder $E_i^{\sigma }$. The energy arising from the external applied field is referred to as $E_i^F$ and $E_i^C$ represents the Coulomb energy. The Coulomb interaction between all particles is calculated using the Ewald sum [12,78]. The Ewald sum makes use of the principle of image charges. As the potential at the contact layers has to be constant, all the field vectors have to end perpendicular to contacts. This can be achieved by removing the contact and introducing an image charge of opposite sign and the same distance behind the contact. Since an OSC device has two contacts, every image charge behind one contact also induces an image charge behind the other contact, leading to a periodic scheme. The electrostatic potential at node i can be calculated by:

$$U_i=\frac{1}{4\pi {ϵ}_0{ϵ}_r}\sum _{\overline{n}}\sum _{j=1}^{2N}{}^{{}^{\prime }}\frac{q_j}{\left|{\overline{r}}_{ij}+\overline{n}\right|},$$

(14)

where the outer sum corresponds to the periodic replica in the x-, y- and z-direction and the inner one sums over all charged particles $q_j$ inside the box $\overline{n}$. The prime in the inner sum notes that self-interactions are excluded. Here, ${ϵ}_0$ is the vacuum permittivity and ${ϵ}_r$ is the permittivity of the material.

The Miller–Abrahams rate provides a generic model for the hopping process modeling the contribution of tunneling and thermal activation by the parameters $a_0$, $\gamma$ and the Boltzmann term, respectively. It is used if hopping is mainly caused by single phonon-assisted tunneling between sites of energy $E_i$ and $E_j$. Using Marcus theory, one implicitly accounts for molecular details including $I_{ij}$ and $\lambda$. It is frequently used in multiscale studies, as for given organic molecules, the molecular positions and transfer integrals need to be calculated by molecular dynamics and quantum chemical calculations, respectively [42,43,79,80]. Marcus theory accounts for multiphonon hopping processes.

Recombination of charge carriers is possible when an electron and a hole are located at neighboring sites. Charge recombination is one of the primary loss processes in OSCs besides exciton decay. When recombination occurs, the two particles are removed, and the sites are free again.

3.3.2. Excitons

An exciton is an intramolecular excited state between a Coulomb-bound electron-hole pair [64]. Compared to the weakly-bound Wannier–Mott excitons present in inorganic semiconductors, strongly-bound Frenkel excitons are formed in organic materials due to the low dielectric constant [3,81]. The excited state is a singlet (triplet) state, if the spins of the electron and the hole are anti-parallel (parallel). We will now describe all physical rates modeled for both singlets and triplets, as well as the conversion rates between these excitons. For a detailed discussion on the photophysical processes, please refer to [64,82]. A schematic overview of the rates and the corresponding energetic levels are depicted in Figure 5. The singlet (blue) and triplet excitons (red), residing at the energetic levels $S_1$ and $T_1$, respectively, as well as the separated electron-hole pair forming a charge-transfer (CT) state and the ground state $S_0$ are visualized.

Singlet excitons can be generated upon light absorption. Absorption of a photon excites one electron, resulting in a singlet exciton. The optical generation profile of singlets is modeled using a transfer matrix model (TMM) [83,84,85]. The TMM is used to calculate the generation profile in the active layer under an AM1.5 illumination spectrum. For every layer in penetration direction z, a constant generation rate $a_{\mathrm{gen}}^s$ is obtained. As the TMM does not provide the generation profile perpendicular to the penetration direction, the x and y position of the generated singlet is chosen randomly [10,13].

The generation of both singlets and triplets can occur by electrical excitation. The rate $a_{\mathrm{exg}}^{s/t}$ for electrical generation of a singlet/triplet is enabled if two charges of different polarity reside at neighboring sites and if the respective spins, a constant property of the charge carriers, show a singlet/triplet configuration. Depending on the origin of charges, we assign their spin states following a distinguished procedure. If a charge is injected, we allocate the spin following the expected singlet to triplet ratio of 1:3. If an exciton dissociates, the spin states of the charges reflect the spins of the singlet or triplet excitation. We use the convention that holes always have a spin with eigenvalue $-\hslash /2$ (spin-down), while the spin of the electrons is assigned such that the above criteria are fulfilled. Following Mesta et al. [23], the electrical generation rate can be calculated as a hopping process of the electron on a hole-occupied neighboring site using the Miller–Abrahams equation (Equation (11)) [74]. In comparison to Mesta et al., our model only considers electron hops for the electrical exciton generation. The energy difference for the singlet/triplet generation is given by:

$$\Delta E^{s/t}=E_{\mathrm{B}}^{s/t}+\Delta E_{\mathrm{LUMO}}-E_{\mathrm{CT},\mathrm{B}},$$

(15)

where $E_{\mathrm{B}}^{s/t}$ denotes the exciton binding energy. $E_{\mathrm{CT},\mathrm{B}}$ gives the charge transfer binding energy, corresponding to the Coulomb binding energy of an electron-hole pair residing at neighboring sites, and $E_{\mathrm{LUMO}}$ defines the LUMO energy level difference of the neighboring sites. Coulomb interactions with the surrounding charge carriers are omitted. If an electrical generation event is chosen, the electron and hole are removed, and a singlet/triplet exciton is generated at the site the former hole was located.

All excited states have a certain lifetime before they decay to the ground state. The lifetime of triplets exceeds that of singlets due to the spin selection rules by far. The rate for the exciton decay is chosen as the inverse of average lifetime $a_{\mathrm{exd}}^{s/t}=1{\tau }^{s/t}$. Due to the spin-orbit coupling of a phosphorescent dopant, triplets may decay due to the lifted spin-selection rules. Thus, the triplet lifetime is reduced within the radius $r_{\mathrm{ISC}}$ around dopant sites.

Excitons are quasi-particles with a net charge of zero. Thus, they diffuse through the organic semiconductor by energy transfer between sites. Existing kMC models account for exciton diffusion using either Förster resonance energy transfer (FRET), Dexter type transfer or are based on a random walk [19,86,87,88]. The diffusion of singlets and triplets can be modeled using FRET [89] or by a random walk based on the Einstein–Smoluchowksi model for the Brownian motion with a diffusion coefficient D [90,91]. Förster described a non-radiative energy transfer between an exciton donor and acceptor as a dipole-dipole interaction [89]. As the Coulomb interaction causes this energy transfer, FRET is a long-distance mechanism. For the random walk (RW), the hopping rate is computed as the inverse of the mean time interval for a hopping process ${\tau }_{\mathrm{hop}}$ between neighboring sites, which gives the rate:

$$a_{\mathrm{exh},\mathrm{RW}}^{s/t}=\frac{1}{{\tau }_{\mathrm{hop}}^{s/t}}=\frac{6D^{s/t}}{l^2},$$

(16)

where l is the next-neighbor distance. Diffusion rates by FRET are modeled with:

$$a_{\mathrm{exh},\mathrm{F}}^{s/t}=\Gamma {\left(\frac{r_F}{r_{ij}}\right)}^6\left\{\begin{array}{cc}exp\left(-\frac{\Delta E_{ij}}{k_{B}T}\right)\hfill & \mathrm{for}\Delta E_{ij}\ge 0,\hfill \\ 1\hfill & \mathrm{for}\Delta E_{ij}\le 0,\hfill \end{array}\right.$$

(17)

where $r_F$ is the Förster radius, $\Gamma$ is the total decay rate, $\Delta E_{ij}=E_i-E_j$ is the exciton energy difference between the initial and final configuration and $r_{ij}$ denotes the distance between two sites. The computational effort of random walk is typically 2–3 orders of magnitude smaller compared to FRET as the number of sites considered for the FRET diffusion increases with $r_F$, while RW diffusion only considers next-neighbors [19]. Especially for a low energetic disorder, RW and FRET give similar results [19], with rising energetic disorder FRET needing to be considered.

In case of the triplet diffusion, FRET is a spin-forbidden process and can only occur between phosphorescent dopants within the Förster radius $r_F$. A second mechanism for the diffusion comes into play. In contrast to the Förster rate, the Dexter type energy transfer rises due to the exchange interaction and may occur between next neighbors [92]. Following [22], the Dexter type hopping is modeled as a Miller–Abrahams hopping defined by:

$$a_{\mathrm{exh},\mathrm{D}}^t=a_{\mathrm{D}}exp\left(-2\gamma r_{ij}\right)\left\{\begin{array}{cc}exp\left(-\frac{\Delta E_{ij}}{k_{B}T}\right)\hfill & \mathrm{for}\Delta E_{ij}\ge 0,\hfill \\ 1\hfill & \mathrm{for}\Delta E_{ij}\le 0,\hfill \end{array}\right.$$

(18)

with the exciton energy difference $\Delta E_{ij}$ of initial state i and final state j, the Dexter hop rate $a_{\mathrm{D}}$ in the zero distance limit and the inverse localization length $\gamma$ of donor and acceptor sites [64]. The triplet binding energy $E_{\mathrm{B}}^t$ is usually assumed to be equal for every site, but can be modified by a disorder following a Gaussian distribution. The Boltzmann factor vanishes in the case of no energetic disorder in the binding energy. Both transfer mechanisms are shown schematically in Figure 6.

To extract electrical energy from the absorbed photons in organic solar cells, singlets need to be dissociated. This process is only possible if singlets reach an acceptor-donor interface. The dissociation process is reported to happen on extremely fast timescales ${\tau }_{\mathrm{exs}}^s$ in BHJs [93]. A separation rate $a_{\mathrm{exs}}^s=1/{\tau }_{\mathrm{exs}}^s$ is enabled if the singlet is located at an interface. Dissociation of a triplet exciton is interpreted as an electron hopping from its current position towards an unoccupied neighboring site [21,22,23]. The separation rate $a_{\mathrm{exs}}^t$ is calculated with the Miller–Abrahams rate using an energy difference of $\Delta E=\Delta E_{\mathrm{LUMO}}+E_{\mathrm{CT},\mathrm{B}}-E_{\mathrm{B}}^t$. The dissociation is probable if the state following a dissociation process is energetically lower than the triplet binding energy.

Regarding high spin-orbit coupling, intersystem crossing from a singlet to a triplet state is possible. The reverse process is not considered here. In our model, a transition from singlets to triplets is only possible within a certain radius $r_{\mathrm{ISC}}$ around a dopant molecule. The radius of intersystem crossing is a material dependent property. If an exciton resides at a site within $r_{\mathrm{ISC}}$, the intersystem crossing rate $a_{\mathrm{ISC}}$ is enabled. The intersystem crossing rate is given by:

$$a_{\mathrm{ISC}}=a_{\mathrm{ISC}}^{0}exp\left(-\frac{2r_{i\to d}}{r_{\mathrm{ISC}}}\right),$$

(19)

with an intersystem crossing frequency $a_{\mathrm{ISC}}^0$. The intersystem crossing decreases exponentially with the distance $r_{i\to d}$ to a dopant site d. If the intersystem crossing occurs, the singlet is converted into a triplet exciton, residing at the corresponding site.

With rising concentration of triplets, another undesired process occurs, the triplet-triplet annihilation (TTA). If two triplet excitons meet, they either create a singlet or a triplet exciton [20,21,94,95]. The spin state of the resulting excitons depends on the spin angular momentum of the two involved triplet excitons. If a triplet exciton resides at a neighboring site or two involved triplet excitons reside at two dopant sites within the Förster radius, TTA is enabled. As explained before, TTA is interpreted as a triplet exciton hopping on a site occupied by a second triplet. If this transition is performed, at least one triplet exciton is removed. If the TTA process is Dexter-mediated, with a ratio of 1:3, the particle at the final site is either a singlet or a triplet, respectively. In the case of the Förster mediated TTA, one triplet is excited to a higher triplet level $T_{\mathrm{n}}$, followed by a relaxation to the lowest triplet exciton level $T_1$. We assume the relaxation to happen immediately. The second triplet decays to the singlet ground state.

4. Applications

An inherent feature of the kMC method in comparison to, e.g., continuum models is the availability of the underlying time-dependent evolution of the modeled system. The kMC method allows the simulation of entire devices including effects at the nanometer scale. Single-particle dynamics, from which both the steady-state quantities and the transient behavior of the system can be extracted, are computed, and particle ensemble characteristics are evaluated.

In this section, we give two exemplary studies on organic materials using kMC. Firstly, the electrical properties of organic materials are investigated. We simulate the charge carrier mobility and current distribution in organic materials showing the effect of a correlated energetic landscape and the inclusion of polymer chains. Secondly, a numerical study of the impact of the dopant volume concentration on the exciton dynamics is presented. An analysis of the conversion efficiency between singlets and triplets is shown in single organic material, neglecting any charge transport properties.

4.1. Charge Carrier Transport through Polymer Chains

The charge carrier mobility is a crucial figure of merit for the applicability of organic semiconductors. While for some devices such as OLEDs, a low mobility is critical to optimize the performance, OSCs and OTFTs demand a high mobility. In semiconducting polymers, charge carrier mobility strongly depends on the morphology [54,55], but its interplay is not well understood [52,96]. Investigating the impact of different transport mechanisms, such as intrachain hopping, on carrier mobility can help the development of organic materials and improve device performance.

A single material cubic lattice of $50\mathrm{nm}\times 50\mathrm{nm}\times 50\mathrm{nm}$ with periodic boundary conditions in all three spatial dimensions is used for the following investigation of the charge carrier mobility under the influence of intrachain hopping. For the energetic disorder, a Gaussian density of states with $\sigma =0.1\mathrm{eV}$ is used. The Miller–Abrahams rate given in Equation (11) is used for the hopping process with an attempt-to-escape frequency $a_0=3\mathrm{e}12{\mathrm{s}}^{-1}$ for inter-chain transitions and an inverse localization length of $\gamma =2{\mathrm{nm}}^{-1}$. Charges are allowed to hop to the six nearest neighbors on the cubic lattice. We run the simulation with four electrons placed in the simulation box resulting in a charge density of $n=3.2\mathrm{e}16{\mathrm{cm}}^{-3}$. This charge density is constant throughout the simulation because of the periodic boundary conditions. To account for the faster intrachain hopping, an attempt-to-escape frequency of $a_{\mathrm{chain}}=m\times a_0$, $m\in {\mathbb{R}}^{+}$, is used to calculate the rates. Further simulation parameters are summarized in

Table 1. Simulation parameters for mobility considering polymer chains.

Description	Parameter	Value
Gaussian energetic disorder	$\sigma$	$0.1$ eV
Attempt-to-escape frequency	$a_0$	3 × 1012 s−1
Inverse localization length	$\gamma$	2 nm−1
Charge density	n	3.2 × 1016 cm−3
Temperature	T	298 K
Simulation time	$t_{\mathrm{sim}}$	6 μs
Applied field	F	5 × 105 − 5 × 106 V cm−1
Dielectric constant	${ϵ}_r$	3.5
Intersite distance	a	1 nm

. For each electric field F, three simulation runs for a simulation time of $6\mu \mathrm{s}$ are performed. The charge carrier mobility in the field direction is calculated by:

$${\mu }_z=\frac{\langle \overrightarrow{v}·\overrightarrow{F}\rangle }{F^2}=\frac{d_z}{t_{\mathrm{sim}}F},$$

(20)

with the electric field $\overrightarrow{F}=F\overrightarrow{e_z}$ in the z-direction $\overrightarrow{e_z}$ and the average carrier velocity $\overrightarrow{v}$ given by the distance $d_z$ traveled during the simulation time $t_{\mathrm{sim}}$ in the field direction. The results represent the mean value and the standard deviation of the mobility.

For this study, different morphologies with chain length $l_c=20\mathrm{nm}$ are generated. We study both Gaussian distributed chains with $\eta =0$ and worm-like chains with $\eta =0.7,0.9$ oriented in the z-direction. The multiplication factor m is varied between one and 100. A factor of $m=1$ means that hopping along polymer chains is not preferred and the chain morphology does not have any impact. Thus, the blue curve in Figure 7a is a reference for comparing the impact of intrachain hopping. Because of the relatively high energetic disorder of $\sigma =0.1\mathrm{eV}$, the carriers are slowed down due to the high energetic barriers between sites. At low fields, this slowing effect is dominant, and transport is mainly determined by diffusion. For higher fields, the impact of the energetic barriers becomes weaker, and transport is dominated by the drift field. For high fields, the carrier velocity saturates, and the mobility decreases with $1/F$.

Figure 7a shows the carrier mobility for completely undirected Gaussian chains. For higher intrachain hopping rates, the mobility increases slightly, but when increasing intrachain hopping rates by a factor of 100, mobility only increases by a factor of 2.2. When including very fast hopping along undirected polymer chains, only a very small amount contributes to the carrier mobility along the field direction. The limiting factor for carrier mobility is the slower inter-chain hopping. Polymer chain orientation the in field direction, on the other hand, shows a larger influence on carrier mobility. In Figure 7b, mobility curves with constant attempt-to-escape frequency $a_{\mathrm{chain}}=100·a_0$ and different orientation are shown. Mobility increases from Gaussian chains ($\eta =0$) to highly oriented chains ($\eta =0.9$) by a factor of 3.7. In total, the impact of polymer chains on carrier mobility adds up to a factor of 8.1 for highly oriented chains with $\eta =0.9$ and increased hopping rates of $a_{\mathrm{chain}}=100\times a_0$.

The combination of polymer chains and the non-cubic Voronoi set allows a realistic representation of the organic materials. In future work, the chains shall be generated from molecular dynamics simulations. With this, kMC can be used to gain insight into the device characteristics based on a realistic lattice model down to the sub-nanoscale.

4.2. Impact of Dopant Volume Concentration on Exciton Dynamics

Using the presented kMC framework, we study the impact of the dopant concentration on the exciton dynamics in an organic bulk material. We focus on the diffusion length and lifetime of excitons and analyze the number of singlets converted into triplets. With this simulation, which is rather theoretical, we aim to mimic the exciton dynamics and interplay of excitons within a single organic material without any interfaces, defects and contacts.

As for the charge transport simulation, a single material cubic lattice of $50\mathrm{nm}\times 50\mathrm{nm}\times 50\mathrm{nm}$ with a lattice constant of $1\mathrm{nm}$ and periodic boundaries in all spatial dimensions is generated. Rates for Dexter energy transfer and random walk are computed for the six nearest-neighbors, while FRET is considered for all neighbors within the Förster radius. Simulations are executed for different dopant volume concentrations varying from $0\mathrm{vol}$% up to $10\mathrm{vol}$%. As no contacts are included, charges cannot be injected or collected. Due to the lack of interfaces, exciton splitting cannot occur. Singlet excitons are generated with a constant optical generation rate $a_{\mathrm{gen}}^s=7{\mathrm{s}}^{-1}{\mathrm{nm}}^{-3}$. Triplet excitons can be induced upon intersystem crossing close to the phosphorescent dopants. In

Table 2. Simulation parameters for the exciton dynamics.

Description	Parameter	Value	Ref.
Common Parameters
Phosphorescent dopant concentration	c	$0\mathrm{vol}\%-10\mathrm{vol}$%	-
Energetic disorder excitons	${\sigma }^{s/t}$	$0\mathrm{eV}$	-
Inverse localization length	$\gamma$	$2{\mathrm{nm}}^{-1}$	[11,12]
Simulation time	$t_{\mathrm{sim}}$	$1\mathrm{ms}$	-
Exciton parameters
Optical generation rate	$a_{\mathrm{gen}}^s$	$7{\mathrm{s}}^{-1}{\mathrm{nm}}^{-3}$	-
Hopping rate (RW)	$a_{\mathrm{exh},\mathrm{RW}}^s$	$2.0\times {10}^{10}{\mathrm{s}}^{-1}$	[13]
Singlet decay rate	$a_{\mathrm{exd}}^s$	$2.0\times {10}^9{\mathrm{s}}^{-1}$	[12,17,97]
Singlet binding energy	$E_{\mathrm{B}}^s$	$0.5\mathrm{eV}$	-
ISC radius phosphorescent dopant	$r_{\mathrm{ISC}}$	$2\mathrm{nm}$	-
ISC rate	$a_{\mathrm{ISC}}^0$	$1.0\times {10}^{12}{\mathrm{s}}^{-1}$	-
Förster radius dopant	$r_{\mathrm{F},\mathrm{dopant}}$	$1.5\mathrm{nm}$	[98]
Dexter rate	$a_{\mathrm{D}}$	$2.9\times {10}^8{\mathrm{s}}^{-1}$	[22]
Decay rate dopant	$a_{\mathrm{exd},\mathrm{dopant}}^t$	$1.0\times {10}^6{\mathrm{s}}^{-1}$	-
Decay rate host material	$a_{\mathrm{exd},\mathrm{host}}^t$	$1.0\times {10}^6{\mathrm{s}}^{-1}$	[24,63]
Triplet binding energy	$E_{\mathrm{B}}^t$	$1.0\mathrm{eV}$	-
TTA rate Dexter type	$a_{\mathrm{D},\mathrm{TTA}}$	$2.9\times {10}^8{\mathrm{s}}^{-1}$	[22]
TTA Förster radius	$r_{\mathrm{F},\mathrm{TTA}}$	$1.5\mathrm{nm}$	[98]

, the simulation parameters used are summarized.

In Figure 8, the dependence of the average lifetime $\tau$ and effective diffusion length $L_{\mathrm{diff},\mathrm{eff}}$ on the dopant volume concentration c is illustrated for both singlet (a) and triplet excitons (b). The effective exciton diffusion length is determined using:

$$L_{\mathrm{diff},\mathrm{eff}}^{s/t}=\sqrt{\frac{\sum {\left(\mathrm{d}{L}_i^{s/t}\right)}^2}{N^{s/t}}},$$

(21)

where $\mathrm{d}{L}_i$ is the Euclidean distance between the generation and final site of exciton i and $N^{s/t}$ denotes the number of either singlet or triplet excitons [82]. The presented results show the average effective diffusion lengths of eight configurations simulated for each dopant concentration, while the error bars indicate the standard deviation. For a concentration of $c=0\mathrm{vol}$%, an average singlet exciton lifetime of ${\tau }^s=501\pm 24\mathrm{ps}$ is obtained, which corresponds to an effective diffusion length of $L_{\mathrm{diff},\mathrm{eff}}^s=7.6\pm 0.1\mathrm{nm}$. As can be seen in Figure 8, both ${\tau }^s$ and $L_{\mathrm{diff},\mathrm{eff}}^s$ decrease with rising dopant concentrations. This is caused by the increased intersystem crossing of singlets into triplets within simulation times below the decay time. Effective triplet diffusion lengths of $L_{\mathrm{diff},\mathrm{eff}}^t=101\pm 3\mathrm{nm}$ are obtained for a low dopant concentration, whereas for increasing dopant concentrations, the diffusion length decreases down to $L_{\mathrm{diff},\mathrm{eff}}^t\approx 65\mathrm{nm}$ for $c=10\mathrm{vol}$%. The average lifetime of triplets decreases from $5.5\pm 0.1\mathrm{s}$ for $c=0.2\mathrm{vol}$% to $2.3\pm 0.1\mathrm{s}$ at $c=10\mathrm{vol}$%. This strong decrease in both effective diffusion length and lifetime is caused by the higher decay rates of triplets within the radius $r_{\mathrm{ISC}}$ around a dopant. With rising volume concentration, the influence of dopants becomes increasingly visible.

Figure 9 depicts the dependence of the percentage of singlets undergoing intersystem crossing and triplets annihilated by TTA on the dopant concentrations. Multiple counting of singlets back-converted from triplets upon a TTA event is omitted. For a dopant concentration of $c=0.2\mathrm{vol}$%, the proportion of converted triplets to optical generated singlets is $30\%$. For higher concentrations, more than $80\%$ of singlets are converted into triplets before decaying. Figure 9 displays the percentage of triplets annihilated upon TTA. Starting with a ratio of $8.1\%$ at $c=0.2\mathrm{vol}$%, the proportion of triplet-triplet annihilation decreases to a value of $2.3\%$ for $c=10\mathrm{vol}$%. Triplet-triplet annihilation becomes unlikely for higher concentrations. The decrease in TTA can be related to the decreasing lifetime of triplets and the reduced diffusion length. Thus, the probability of triplets to meet each other decreases tremendously. This is explicitly caused by the lowering of the spin-selection rules within the strong spin-orbit coupling domains of the dopants. Higher triplet densities would significantly increase the impact of TTA. In Figure 9a, the effective exciton lifetime defined by:

$${\tau }_{\mathrm{eff}}=\mathrm{ISC}\times {\tau }^t+(1-\mathrm{ISC})\times {\tau }^s$$

(22)

is visualized. For a low dopant concentration, the ${\tau }_{\mathrm{eff}}$ increases strongly as most of the excitons are not converted into triplets. With rising c, the effective lifetime rises and reaches its maximum at $c=1\mathrm{vol}$%. Beyond the maximum, the effective lifetime decreases and saturates at a value of ${\tau }_{\mathrm{eff}}\approx (2.09\pm 0.06)\mu \mathrm{s}$. A similar trend is observed by Gonzales et al. for doped P3HT:PCBM solar cells [24].

The results for exciton diffusion lengths and lifetime conform with the experimental values within ideal material setups, documented in [63,64,86,99]. Similar to the trend observed in [24,26], the lifetime of excitons decreases with dopant concentration. Firstly, for a small dopant concentration, the lifetime increases tremendously in comparison to singlet lifetimes caused by the conversion of singlets to the long-living triplets. However, for a higher dopant concentration, a quenching mechanism is observable due to the TTA and a higher decay rate. Simulations with higher triplet exciton decay rates within the ISC radius due to stronger spin-orbit coupling matches observed results for exciton lifetime trends [26]. Yang et al. measured the decrease of the triplet exciton lifetime by a factor of 2.8 upon doping with Ir(mppy)${}_3$, which is in the same order of magnitude of our simulation results [26].

5. Conclusions

To summarize, we have presented a generalized kinetic Monte Carlo framework for the study of organic semiconductors and optoelectronic devices. New developments on a realistic representation of the organic material by a generalized non-cubic Voronoi set and polymer chains have been emphasized. We have included dopant sites and improved our description of the exciton dynamics by a model of the triplet dynamics for the detailed study of OSCs and OLEDs. A simplified spatial correlated density of states accounting for a correct treatment of the energetic landscape has been presented. Our kMC framework overcomes limitations of existing kMC models due to both a generalized spatial and energetic representation of the underlying multifaceted organic materials, as well as improved treatment of the particles involved in the photocurrent generation of OSCs and the light emission of OLEDs. We presented two applications concerning electrical and optical processes inside an organic semiconductor. The impact of the intra- to inter-chain hopping rate and the directionality of the polymer chains on the charge carrier mobility has been studied. We showed that an enhanced intrachain hopping rate itself does not improve the mobility of charges, but directionality is necessary. Furthermore, the impact of the dopant volume concentration on the effective lifetime and diffusion length of excitons has been studied. The investigated exciton dynamics coincides with reported values of both singlet and triplet lifetimes and diffusion lengths.