Experimentally Calibrated Kinetic Monte Carlo Model Reproduces Organic Solar Cell Current-Voltage Curve

Kinetic Monte Carlo (KMC) simulations are a powerful tool to study the dynamics of charge carriers in organic photovoltaics. However, the key characteristic of any photovoltaic device, its current-voltage ($J$-$V$) curve under solar illumination, has proven challenging to simulate using KMC. The main challenges arise from the presence of injecting contacts and the importance of charge recombination when the internal electric field is low, i.e., close to open-circuit conditions. In this work, an experimentally calibrated KMC model is presented that can fully predict the $J$-$V$ curve of a disordered organic solar cell. It is shown that it is crucial to make experimentally justified assumptions on the injection barriers, the blend morphology, and the kinetics of the charge transfer state involved in geminate and nongeminate recombination. All of these properties are independently calibrated using charge extraction, electron microscopy, and transient absorption measurements, respectively. Clear evidence is provided that the conclusions drawn from microscopic and transient KMC modeling are indeed relevant for real operating organic solar cell devices.

organic solar cell. It is shown that it is crucial to make experimentally justified assumptions on the injection barriers, the blend morphology, and the kinetics of the charge transfer state involved in geminate and nongeminate recombination. All of these properties are independently calibrated using charge extraction, electron microscopy, and transient absorption measurements, respectively. Clear evidence is provided that the conclusions drawn from microscopic and transient KMC modeling are indeed relevant for real operating organic solar cell devices.

Introduction
Kinetic Monte Carlo (KMC) simulations have successfully been used to model the charge carrier dynamics in organic photovoltaics (OPVs) on the ps to µs time scale. For instance, it was shown that in thin-film OPV devices, thermalization in the disorder-broadened density of states (DOS) does not complete before charges are extracted. [1][2][3][4] The conclusions from these studies are drawn from the fitting of time-resolved experiments performed under certain bias conditions such as short circuit or open circuit. Other authors used KMC modeling to focus on the process of charge recombination and its dependence on the morphology in slabs of material, i.e., in absence of contacts. [5][6][7][8][9] However, it is still an open question to which extent nonequilibrium phenomena and other aspects that are not accounted for in macroscopic simulations such as quasi-equilibrium drift-diffusion (DD) models, govern the steady-state operation of complete OPV devices. To answer the question, it would be highly desirable to have a microscopic model that is also able to describe the current-voltage (J-V ) curve, particularly the open-circuit voltage (V OC ) and the fill factor.
Modeling J-V curves with KMC has so far proven nearly impossible. One of the main challenges is the presence of two injecting contacts. While it may be acceptable to consider the contacts as simple sinks for electrons and holes in transient extraction experiments (performed at V V OC ), this simplification does not work for situations closer to V OC . When the internal field is low, contacts inject many charge carriers into the active layer. This high carrier density is demanding from the computational point of view and challenging to correctly account for. Even though a few concepts exist how contacts can be implemented in KMC, literature studies have so far failed to fully describe J-V data of real devices or are based on assumptions that are not justified experimentally. [10][11][12] Besides computational challenges, the injected charge density also sets the boundary conditions for the recombination of photogenerated carriers. [13] Charge recombination generally becomes more important when going from short circuit to open circuit because transport will slow down. Indeed, the competition between charge extraction and recombination has been demonstrated to be the main determinant of the device fill factor. [14][15][16] For a device model to be reliable it must therefore capture the hopping transport characteristics and the recombination kinetics at the same time. Even though the mechanisms of charge recombination are highly disputed, it is commonly accepted that the morphology plays a key role. [6,17,18] For instance, it is well documented that aggregated donor or acceptor domains may lower the recombination rate. [19][20][21][22] However, although the morphology of many donor/acceptor blends is well characterized by electron microscopy and other techniques, the nanostructure is often neglected in KMC and an effective medium is assumed instead. [23] Here, we present a KMC model that successfully predicts device J-V curves while simultaneously accounting for nonequilibrium hopping transport and recombination dynamics.
We show that this is only possible when correct assumptions are made on the injection barriers, the morphology of the active layer, and the charge recombination rate. All these properties are calibrated by independent experimental techniques such as charge extraction, electron microscopy and transient absorption. We are thereby introducing a device model that works on a multitude of length and time scales. As such it will be useful for future investigations on the interplay between elementary processes and device characteristics of organic solar cells and other optoelectronic devices.

Material System
The aim of this work is to develop and experimentally calibrate a KMC model that fits both transient experiments and device J-V curves. Our material system for experimental calibration is TQ1:PC 71 BM, [24] an archetypal polymer/fullerene blend. The reason for choosing TQ1:PC 71 BM is that for this specific system a clear picture of the carrier dynamics has emerged from time-resolved measurements and previous modeling, which is summarized in a recent review article. [23] Hence, many of the parameters for the KMC model are already known; in particular, it has been shown that the charge extraction in thin devices with an active-layer thickness ≈ 100 nm is strongly affected by nonequilibrium effects. Figure 1 shows

KMC Describes Device Current-Voltage Curve
The KMC model, which is extended and experimentally calibrated in this work to fully describe OPV devices, has been introduced previously. [ Table 1. Open symbols are simulations with the same parameter set, but assuming too low injection barriers of 0.1 eV (triangles) or only an effective medium without PC 71 BM aggregates (squares). hopping rate ν ij from site i to site j separated by a distance r ij is given by where ν 0 is the attempt-to-hop frequency, α the inverse localization length, ∆E ij = E j − E i the energy difference between the sites, and kT the thermal energy. Hopping is assumed to take place in a Gaussian DOS, where E is the single particle energy, E 0 the mean energy, and σ the width of the Gaussian DOS or the energetic disorder. We note that without of loss of generality, also other energy distributions could be assumed in the model, such as an exponential DOS. From previous studies, however, it is known that a Gaussian DOS gives the most appropriate description for the present TQ1:PC 71 BM system both when describing transient and steadystate experiments. [23,25] In this work, only hopping between nearest neighbors on a regular, six-fold coordinated lattice was considered. In this configuration, the localization length α is unimportant; the first exponential term of Equation (1) was implicitly included in ν 0 , that is, the rate of downward nearest-neighbor hops.
The core working principle of a KMC model is to simulate the time evolution of a system As mentioned above, the presence of injecting contacts causes computational challenges. Charge injection is mediated by the injection barriers, i.e., the energy offset between the contact Fermi level and the respective molecular orbital of the semiconductor. Especially for low barriers, carriers may oscillate multiple times across the contact interface before injection/extraction finally takes place. We mitigated this 'small barrier' problem by only allowing for a transfer if the number of charges next to the contact interface deviates from its equilibrium value. The transfer is modeled as hopping event with an attempt frequency ν 0,cont of the same order as for the transport of the faster carrier (here: electrons) in the semiconductor. This ensures that charge collection is not limited by the contacts. Both the cathode and anode were considered nonselective; hence, possible losses due to diffusion of carriers into the 'wrong' contact are implicitly accounted for.
An advantage of KMC simulations is that no explicit assumptions about the formalism of charge recombination need to be made. Recombination of free charges involves the formation of a CT pair as intermediate. Exciton formation is explicitly allowed, but requires overcoming the relevant energy level offset between the TQ1 and PC 71 BM; as such, it can be interpreted as the inverse of charge separation, i.e., the splitting of (CT) excitons into free electrons and holes. As discussed in more detail below, it is then the inverse lifetime of the CT state that determines the recombination rate and must be calibrated experimentally.
The filled circles in Figure 1a show that after the calibration discussed below, the KMC model fits the J-V curve of the TQ1:PC 71 BM solar cell well within experimental accuracy and matches both the device V OC and fill factor. Table 1 lists the key parameters used for the simulations. We note that these values are not the result of a fitting routine but come from independent characterizations. The hopping parameters were chosen in such a way that they represent earlier experiments, such as time-resolved electric-field-induced second harmonic generation (TREFISH) [4] and temperature-dependent space-charge-limited currents (SCLC), [4,26] but at the same time allow efficient calculations. This was done by assuming a single disorder for electrons and holes (σ e = σ h ≡ σ) and adjusting the attempt frequencies ν 0 such that the macroscopic transport characteristics of TQ1:PC 71 BM, e.g., the contrast between electron and hole mobility, are still captured (see Supporting Information for details). Figure 1b shows that also with the symmetrized hopping parameters, relaxation in the DOS is far from being complete when photogenerated carriers are extracted. This is true for both short-circuit and open-circuit conditions, which indicates that nonequilibrium effects may affect charge extraction along the entire J-V curve. A detailed discussion of how the nonequilibrium effects influence the individual performance parameters will be the topic of another publication.
The main result of this study is that a KMC model that can describe full J-V characteristics requires an appropriate description and calibration of the injection barriers and the morphology in the active layer. If wrong or too simple assumptions are made on these properties, our otherwise well validated KMC model can no longer describe the device (Figure 1a, open symbols). Because this mainly concerns V OC and the fill factor, these observations are closely related to the charge recombination. In the following sections we will therefore focus on the factors that determine the shape of the J-V curves in the fourth quadrant, that is the injection barrier height, the blend morphology, and the recombination rate.

Calibration of Injection Barriers
The injection barriers set the carrier density in the device around the built-in voltage. To get a realistic estimate of the barrier height, we compare the results of charge extraction experiments in the dark with device simulations. As KMC calculations are computationally too expensive for this approach, we used a DD model instead. [27,28] This is justified because the charges treated here were not photogenerated, but injected from the contacts, so that the complexities of exciton/charge separation are bypassed. Furthermore, charges are injected from thermalized reservoirs (contacts), so that it is reasonable to describe them by a quasi-equilibrium mobility. The mobility values were estimated by inserting the hopping parameters in Table 1 in the mobility functional by Pasveer et al. [29] Charge recombination is assumed to be strictly bimolecular with the steady-state recombination coefficient (6 × 10 −18 m 3 s −1 ) taken from experimental studies on TQ1:PC 71 BM. [30,31] Figure 2a illustrates the effect of the injection barrier height on the average carrier density. Here, we chose devices with an active-layer thickness of 150 nm; only at these larger thicknesses a 'bulk' region is established, which makes the comparison with charge-extraction experiments more reliable. [32,33] Note that especially at higher densities the carrier profiles are not perfectly symmetric, which is due to the imbalanced electron and hole transport. [28,34] The experiments to be simulated are charge extraction by linearly increasing voltage (CE-LIV) and bias-assisted charge extraction (BACE). In both techniques, the device is held at a certain pre-bias (V pre ) until a steady state is reached; the charges in the device are then extracted by applying a triangular (CELIV) or rectangular (BACE) voltage pulse. The dark carrier density is calculated from the transient current J(t) via where q is the elementary charge, d the active-layer thickness, J 0 the displacement current measured at V pre = 0, and t f the time at which charge extraction is completed. Note that in the form of Equation (3), the carrier density represents the average of electrons and holes, as pointed out by Hawks et al. [35] Figure 2b shows that CELIV and BACE give a consistent picture of the carrier density as a function of voltage. At V pre = 0.9 V, which approximately corresponds to open-circuit conditions under 1-sun illumination, n dark is about 1 × 10 22 m −3 . This is the same order of magnitude as for the photogenerated carrier density and indicates the importance of injected carriers for charge recombination. As can be seen, the best description of the dark carrier density and its voltage dependence is obtained for a barrier height of 0.2 eV; with this value, the KMC model reproduces the experimental J-V curve (see Figure 1). We note that the discrepancy between CELIV/BACE and DD simulation at voltages well below the built-in voltage is merely due to experimental limitations. In this regime, most carriers are situated in the thin space-charge regions close to the contacts, which makes them only partly visible to charge-extraction experiments. [14,36] If instead too small injection barriers are selected as input for the KMC model, it can no longer describe both V OC and the fill factor. The open triangles in Figure 1 illustrate this for a barrier height of 0.1 eV. Although this is not pursued further in this work, we would like to stress that this finding shows that defining a contact as 'Ohmic', in the sense that it does not limit injection and extraction in a particular experiment, is insufficient. Here, injection barriers of 0.1 and 0.2 eV both give rise to 'Ohmic' injection, implying bulk-limited transport under forward bias, but these barriers are not equivalent in terms of the resulting photovoltaic behavior.
Another interesting observation is that, as one would expect, lowering the injection barriers from 0.2 to 0.1 eV leads to an increase in V OC . But at the same time the fill factor becomes reduced, so that the overall power conversion efficiency stays roughly the same. Hence, we can deduce from our KMC simulations that reducing the injection barriers does not per se lead to a better performing solar cell device. Closer examination of this aspect, however, requires more extensive parameter studies, which are beyond the scope of the present paper and will be the subject of future work.

Morphology Governs Charge Recombination
In our previous KMC studies the photoactive blend was assumed as an effective hopping medium without any morphological features. [1,3,4] This zero-order approximation is reasonable when describing experiments on the ps-µs timescale where charge recombination is insignificant. However, we find that the effective-medium approach fails to fully describe the device J-V curve (Figure 1, open squares). In order to obtain a more realistic picture of the morphology, we performed transmission electron microscopy (TEM). Figure 3a  blends are commonly attributed to fullerene domains because of their higher density. However, this assignment is not unambiguous; the different intensities could also be caused by phase contrast due to local crystallinity differences. For comparison, we investigated the same sample in scanning transmission electron microscopy (STEM) mode using a high-angle annular dark field (HAADF) detector. [37,38] In the Supporting Information we show that HAADF-STEM reveals very similar structures as in Figure 3a, but of inverted contrast.
This clearly confirms that the clusters seen in TEM are PC 71 BM aggregates, in agreement with earlier work on similar blend systems. [39] The main effect of aggregation is to reduce the energy gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) compared to the amorphous material. This creates an energy cascade with a driving force for carriers to move from the (molecularly mixed) amorphous regions towards the (material-pure) aggregates and will affect the way how charges separate and recombine. [19][20][21][22] We implemented the aggregates in the KMC model as 7 × 7 inclusions in a 10 × 10 unit cell describing the mixed donor/acceptor phase (Figure 3a, inset). Inclusions were assumed to consist of pure PC 71 BM with a 0.2 eV lower-lying LUMO compared to the mixed phase; all other properties were left unchanged to keep the number of unknown parameters at a minimum. We did not consider pure TQ1 domains, as our TEM experiments do not provide any evidence for them. This is reasonable, since TQ1 is a relatively amorphous polymer that has no strong tendency to form aggregates, in particular in blends with excess fullerene. [40] Note that the aggregate size in the KMC model is smaller than what is suggested from the electron microscopy images. This was done to keep the simulation box computationally tractable while still getting reasonable statistics. The size of the inclusions and the unit cell were chosen such that the donor/acceptor ratio of the blend is maintained. A detailed examination of the structure size on the device performance is beyond the scope of this work; however, first tests indicate that the actual size of the aggregates is much less important than their presence. Likewise, a 0.1 eV lower-lying LUMO for the aggregate phase did not make any relevant difference as compared to the used 0.2 eV.
Only with the inclusions in the effective hopping medium we were able to match the fill factor of the experimental devices. Figure 3b shows that this is due to a reduction of the charge recombination. Importantly, the presence of aggregates simultaneously reduces the yields of geminate and nongeminate recombination. This confirms earlier suggestions that the generation and recombination of free charges are coupled via the ability of CT pairs to separate. [41,42] In other words, the possibility for carriers (here: electrons) to lower their energy by moving to the aggregates will not only increase the charge separation yield, but also reduce the nongeminate recombination. This is a clear hint that the different ability to form aggregates/phase-pure domains may explain why different OPV materials show so different recombination rates compared to the Langevin model. In the context of this work,  Table 1 but varied decay rate of CT states into the ground state. Inset: Illustration of the relevant states and transitions for charge recombination. (b) Measurements with various initial carrier densities by varying the pump fluency from 2.5 × 10 16 to 6 × 10 17 photons · m −2 and simulations for a fixed decay rate of k CT = 3 × 10 7 s −1 . Experimental data from Andersson et al. [43] however, it means that it is the kinetics of the CT states, i.e., how they dissociate and (re-)associate, that must be calibrated experimentally.

Calibration of the Recombination Rate
The inset in Figure 4 illustrates the kinetic model of charge recombination that has emerged from literature. [2,19,30,41,42] As we discuss to some detail in the Supporting Information, recombination in TQ1:PC 71 BM is not limited by the rate k enc at which free carriers meet to form an interfacial CT complex. This implies that the probability for the CT pair to dissociate is much higher than to decay to the ground state (k d k CT ). It has been shown that in such a situation an equilibrium between CT states and free charge carriers is established. [2,41] The position of the equilibrium is determined by the rate k CT , which is the relevant parameter in the KMC model to calibrate the recombination.
In order to do so, we use the results of transient absorption (TA) experiments. TA is a pump-probe technique that optically tracks a carrier population created by a short light pulse over time. As the experiment is carried out under flat-band conditions, the measured decay solely reflects the recombination kinetics. Figure 4a shows the TA decay of a TQ1:PC 71 BM device for a pump fluency of 4 × 10 16 photons · m −2 taken from literature. [43] The traces are attempts to describe the experiment with our KMC model. One can clearly see that the (inverse) CT state lifetime is the crucial parameter for the decay dynamics. The best fit on short time scales is obtained for k CT = 3 × 10 7 s −1 . Figure 4b demonstrates that with the calibrated value for k CT , we are able to reasonably describe transient absorption data for a range of initial carrier densities.
On longer time scales, however, the fit between TA experiment and KMC model is not as good. The reason for this are the symmetrized transport parameters we use for computational effectiveness. As discussed in the Supporting Information, the disorder σ and attempt-to-hop frequency ν 0 are largely interchangeable, i.e., increasing the one parameter can be compensated by decreasing the other and vice versa. This interchangeability allows us to use the values given in Table 1, which keep the KMC calculations manageable, while still reproducing the measured quasi-steady state mobilities. Nevertheless, using symmetrized transport parameters remains a simplification, so that some of the details necessary to describe the full TA traces are lost. In Figure S2 we show that a better fit can be obtained when the 'real', non-symmetrized values for σ and ν 0 are used in the simulation. However, significant differences between the parameter sets are only noticeable at very high initial carrier densities (∼ 10 24 m −3 ) and on the time scale of µs and beyond. At those times, most of the carriers have already been extracted, as can be seen from the histograms in Figure 3b and from previous experiments. [1,4,23] Hence, the use of the simplified transport parameters is well justified when describing a solar cell under standard operating conditions.

Conclusions
We were taken at an acceleration voltage of 300 kV in a FEI Titan 80-300.

Discussion of Symmetrized Hopping Parameters A well-known problem of kinetic
Monte Carlo (KMC) simulations of hopping transport in energetically disordered media is that calculation times explode when the normalized disorderσ = σ/kT becomes larger than 3-4. An additional problem arises in bipolar systems if the attempt to hop frequencies ν 0 of the two charge species (electrons and holes) differ by more than roughly an order of magnitude. In this case, the hopping events of the species with the highest ν 0 outnumber those of the slower species. This leads to poor statistics for the slower species that can only be cured by increasing the total number of hops considered in the simulation, i.e. by increasing the total calculation time. Unfortunately, the previously determined hopping parameters for the TQ1:PCBM system studied here cause both problems to arise.
For holes in TQ1 we found in Refs. S1-S5 σ h ≈ 50-100 meV and ν 0,h ≈ 0.1-1 × 10 10 s −1 and for electrons in PCBM σ e ≈ 120 meV and ν 0,e ≈ 1×10 13 s −1 . In order to keep calculation times to acceptable levels (days), we had to use partially symmetrized hopping and disorder parameters. In this, we could make use of the previously observed interchangeability of the disorder and the attempt-to-hop frequency, [S2] where increases in one of the two parameters can be largely compensated by a simultaneous increase in the other parameter. Bearing the above in mind, we used σ e = σ h = 75 meV for both electrons and holes, and ν 0,e = 1×10 11 s −1 and ν 0,h = 1 × 10 10 s −1 , keeping the corresponding steady state electron mobility significantly larger than that of the holes. The other parameter that has (a minor) influence on the hopping process in the used model is the nearest neighbor distance a NN = 1.8 nm that we fixed to the value that we obtained from temperature-dependent charge transport studies. [S5] A consequence of the symmetrized hopping parameters is that the other rates in the simulation, specifically the recombination rate of the CT exciton, become relative to these values. This explains why the used CT rate k CT = 3 × 10 7 s −1 is two to three times lower than the typical values used before. [S1,S3]

S2
Charge Recombination in TQ1:PC 71 BM The three relevant transitions for the recombination of photogenerated charges are the separation rate of CT states (k d ), the encounter rate of free electrons and holes (k enc ), and the decay rate of CT states into the ground state (k CT ). Importantly, the encounter complex formed by two independent carriers has been identified as the same CT state as involved in charge separation. [S6,S7] Only if the rate at which CT states recombine were much faster than the rate at which they dissociate (k CT k d ), the recombination would be encounter-limited and Langevin theory applicable. In practice, however, virtually all OPV blends show recombination rates that are substantially reduced compared to the Langevin model. The apparent steady-state recombination rate constant of 6 × 10 −18 m 3 /s we assume Section 2.3 in the main text for TQ1:PC 71 BM implies a reduction of 2 orders of magnitude.
Such large reduction factors cannot be explained by geometrical confinement, i.e., the fact that electrons and holes are attributed to different material phases and only meet at the heterointerface. [S8] We must therefore assume that recombination in TQ1:PC 71 BM is not encounter-limited. This is reasonable considering the following: although the details of charge separation are not yet well understood, the rate k d will depend on the ability of the two involved carriers to move away from each other. Still being a hopping process, this will happen with a rate on the order of 10 10 to 10 13 s −1 . Because this is much faster than the rate at which the CT state decays into the ground state, i.e., its inverse lifetime (k CT = 10 7 to 10 8 s −1 ), there is enough time to establish an equilibrium between CT states and free electrons and holes. [S1, S9] In other words, the carriers forming the CT state have multiple attempts to escape from their mutual Coulomb attraction before they ultimately recombine. One can easily see that it is then the rate k CT which determines the position of the equilibrium; decreasing k CT will shift the equilibrium more towards dissociation, effectively decreasing the charge recombination rate detectable experimentally.
S3 Figure S1: Transmission electron micrograph in HAADF-STEM mode of the same TQ1:PC 71 BM blend film as discussed in the main text. The HAADF signal in STEM originates from electrons that have inelastically scattered to high angles when transmitted through the sample. The signal intensity is higher if the electrons scatter against molecules of higher average atomic number (Z-number). This is likely not relevant here, since the TQ1 monomer and the PC 71 BM molecule have a similar average Z-number. The signal intensity is also increased for regions of higher density or thickness. Since the film thickness is rather uniform, the bright regions in the HAADF-STEM image show areas of higher density. This correlates well with what is known about fullerenes and their aggregates. [S10-S12] For the electrons σ e = 120 meV and ν 0,e = 1 × 10 13 s −1 was assumed and for the holes σ h = 83 meV and ν 0,h = 5 × 10 9 s −1 , according to previous experimental work. [S3] As mentioned in the text above, using these parameters required us to assume a slightly higher CT rate of k CT = 8 × 10 7 s −1 .