Organic solar cells

As climate change and energy sufficiency are progressively becoming more pressing issues, environmentally friendly and cheap energy sources need to be developed. Solar energy production, with inorganic solar cells, has progressed greatly in recent years and is already able to compete with traditional energy sources. Photovoltaic solar cells convert photons into electric current. Although traditional photovoltaic technology based on inorganic materials has become commercially successful, it faces some limitations that can be overcome by organic solar cells. Organic solar cells are low-cost and easy to process, furthermore they possess innovative properties being potentially lightweight, flexible, and transparent. However, organic solar cells still have lower efficiencies and shorter lifetimes than traditional inorganic solar cells. While the best organic solar cells have reached around 11% efficiency, the best single junction crystalline silicon solar cells and thin film CdTe cells have efficiencies of around 25% and 22%, respectively.^1,2  Furthermore, the lifetime of organic solar cells is still short in comparison to the lifetimes of inorganic solar cells, so stability challenges must be addressed for organic solar cells to compete with conventional photovoltaics on the market.^3–5  The efficiency of an organic solar cell is determined by the efficiency of the different steps from photon absorption to charge collection, as shown in Fig. 1. Photons are absorbed in either the acceptor or the donor phase with efficiency η_abs, generating excitons (1). These excitons then either diffuse toward the interface and form a Charge Transfer (CT) state (2) or decay (6). The CT state is defined as an electron–hole pair Coulombically bound, composed usually by an electron on an acceptor molecule and a hole on a donor molecule. In an organic solar cell CT states (3) are formed at the interface from exciton dissociation. The efficiency of the excitons forming CT states is η_ex, and it is high when the distance to the interface is small. The CT states then split (4) into free electrons and holes with efficiency η_diss. After the charge carriers are freed, they may still move back to the interface and recombine in a CT state (8), or even to the ground state (7) and (9). η_trans is the efficiency of the free charge carriers being collected at the electrodes (5). The product of these different efficiencies gives the external quantum efficiency (EQE), the ratio of charge carriers collected to the number of incoming photons.^6–12 

The efficiency of organic solar cells is greatly impacted by the processes that occur at the interface between the donor and acceptor. The electron and hole have a strong Coulombic attraction, which must be overcome for the charges to separate. Organic solar cells have much lower dielectric constants and more localized electronic states than inorganic solar cells, so excited states are localized and the Coulombic barrier to overcome for the charges to dissociate is large.

The amorphous nature of these materials make analytical studies difficult and theoretical investigations of mobility and charge transport are to a large extent based on simulations: drift diffusion and Kinetic Monte-Carlo simulations (KMC). Drift–diffusion models^13–16  employ a classical picture, they are very suitable to simulate the whole solar cell and give some macroscopic information. Kinetic Monte-Carlo simulations^17–28  have potentially the capability to model the quantum phenomena happening at the nanoscale, but the more details are included in the KMC scheme the more computational effort is required. Detailed KMC schemes are thus mainly limited in studying some interesting portions of a solar cell. Several studies employing both approaches showed the significant contribution of morphology to the efficiency of solar cells.^14–16,18  Different structures for organic solar cells have been researched, the simplest consisting of a single layer of an organic semiconductor between two electrodes. However, solar cells with this structure have low efficiency, and the performance of the device is improved by a bilayer structure of two organic materials: a donor and an acceptor materials.²⁹  With this design, a trade-off must be made between light absorption and exciton dissociation. The thicker the layers are, the more light will be absorbed, increasing the efficiency of the solar cell. However, the excitons formed in the single material, will have in average a greater distance to the donor–acceptor interface, where they can dissociate in a CT state. If the domain size is too large, the exciton will decay before reaching the interface, thereby lowering the efficiency of the solar cell. A solution to this problem is to combine the donor and acceptor phases so that the distance to the interface is short even with thick films. This can be best achieved with interdigitated structures,^14,16  providing good exciton dissociation and at the same time a clear path to the electrodes for the free charge carriers. The interdigitated interface is very promising but also very difficult to obtain in practice. Another type of interface commonly used in the lab employs a blend of donor and acceptor materials forming a three-dimensional interpenetrated structure called bulk heterojunction (BHJ) solar cell.³⁰ 

In the present chapter we will focus on the KMC approach used to model charge transport in organic materials. The properties of organic materials are different from those of crystals; their intrinsic disorder tends to localize the charge carriers on one or few molecules. The conduction in these materials is therefore temperature activated, in contrast to the band conduction of crystalline materials. Two methods are predominantly used in literature to study charge transport using KMC simulations: the Miller-Abrahams and the Marcus hopping rates. The Miller-Abrahams³¹  formula is one of the simplest ways to couple temperature to the hopping rate. This is achieved by means of a Boltzmann factor helping to overcome the energy barrier for the transport. The Marcus hopping rate³²  can be derived from the Fermi golden rule and takes into account also the reorganization energy after each hop.

In this book chapter, we will show how to model the steps (3), (4) and (5) of Fig. 1, namely, how the CT state split and the charges are subsequently transported to the electrodes. We will thus not consider the transport of excitons¹²  or the coupling between the excited state, the CT state and the ground state^33–35  and we refer the reader to the appropriate literature. After a brief discussion on the method to obtain realistic interfaces with atomistic details, we will focus on the transport through Marcus equation and on the way to calculate the different parameters involved in it. We will also analyze the effect of the electric field on CT state splitting and conduction of free charge carriers. In particular, in order to optimize the device efficiency, we will illustrate how to choose the electric field to favor simultaneously both processes, if this is possible.

The morphology of the donor–acceptor interface in a BHJ organic solar cell has a large impact on the efficiency of the solar cell. Excitons can only diffuse 10–20 nm before decaying, so the donor and acceptor should be sufficiently mixed so that the excitons can reach the interface before decaying. However, once separated, the charge carriers need pathways to their respective electrodes. If, for example, an electron is in an acceptor domain that is completely surrounded by the donor, then the electron will have no path to the electrode and it will eventually decay through recombination. In a bilayer structure, there is always a pathway for free charges to travel to the electrode, but many excitons will decay before they can reach the interface and thus have a chance to split. In the bulk heterojunction structure, however, it is much easier for excitons to reach the interface before decaying, but some of the charges will get trapped in domains that have no path to the electrode or recombine with other charges when traveling by an interface. Two types of recombination can occur: geminate or nongeminate. Geminate recombination is when two charge carriers resulting from the absorption of the same photon recombine (7 in Fig. 1). Nongeminate recombination (9 in Fig. 1) occurs when two free charge carriers originating from different photons recombine with each other at the interface. When there is high interpenetration, free charge carriers are likely to travel close to an interface and recombine with other free charge carriers, increasing the nongeminate recombination rate. It is important to optimize the interpenetration of the donor and acceptor materials in order to extract as many charge carriers as possible from the photogenerated excitons.

Charges hop between molecules or polymer segments, and this hopping occurs with a certain probability, making probabilistic methods such as KMC appropriate for modeling charge transport. In early works^17,18  the molecules of an organic material were modeled only as a lattice of molecular sites, and the structure of the molecules was neglected. More recent works have improved upon these models by including an atomistic description of the molecules.^{19–28,36–39}  To study the complex interplay of geminate and non-geminate recombination and thus obtain a meaningful simulation of charge transport in the solar cell, a realistic morphology is needed. Several methods can be used to model the interface between donor and acceptor materials. As a first approximation, the crystal structure of the two materials can be brought into contact and the interdistance can be optimized without relaxing the molecule inside each phase.^38,39  In order to add more disorder, it is then possible to relax the interface by using molecular dynamics simulations like it has been done in several studies.^35,40  However, if this technique will allow to introduce disorder and create mixing at the interface, this latter will remain quite planar and it will not reach a real three-dimensional interpenetrated structure. In order to obtain a realistic BHJ interface we are currently working on another approach. The idea is to start from volumes obtained with the Ising model developed by Heiber et al. with different degrees of intermixing,⁴¹  and subsequently fill these volumes with donor and acceptor molecules as illustrated in Fig. 2. From those initial structures, molecular dynamics in the NPT ensemble can be performed to relax the different boxes that can be then used for Kinetic Monte Carlo simulations.

Kinetic Monte Carlo simulations are one of the prominent tools for the simulation of charge carriers (single or multiple) in organic materials. In the notation used in this chapter, a particular molecular orbital level in the system is identified by an upper case letter M=(i, m), i.e. the orbital level m of the molecule i. The molecular orbital M will always implicitly belong to the molecule i, unless explicitly stated. In the same way, the orbital N=(j, n) belongs to the molecule j (where i and j are different molecules). The flowchart of the KMC algorithm is presented in Fig. 3. Since organic materials exhibit usually a disordered amorphous structure, it is common practice in the literature to obtain a realistic structure as outcome of a MD simulation (see previous Section 2). Once the structure has been defined, every charge carrier a involved in the simulation is placed in its respective initial positions (step 1). At every time step the hopping rates for all the charge carriers are calculated with Marcus formula³²  (step 2). If a charge carrier is situated on the molecular orbital M, its probability to hop to the orbital N is

Marcus formula involves temperature T (through the Boltzmann constant k_B) and three main physical quantities depending on the particular hop: the transfer integral H_NM, the site energy difference ΔE_NM and the reorganization energy λ_NM. The evolution of the system is determined thanks to a weighted random selection algorithm (WRS), in which the probability of selecting a particular hop depends on its associated weight. First, the hopping charge carrier is selected using as weight for every charge carrier a the total escape rates R_Ma (step 3). This selects the charge that is jumping, thus identifying the origin of the hop (step 4). Let us call c the selected hopping charge and let us assume such charge is placed on the orbital M. Now, among all the possible hops that charge c can make, we will select the one determining the evolution of the system through another WRS using as weights the Marcus rates w_MN (step 5). This hop will be considered in the simulation as happening after a time Δt sampled from an exponential distribution ∼R_Mexp(−R_MΔt), with R_M being the total escape rate from the orbital M, i.e. (step 6). At step 7 the system is updated with the selected hop and the time of the simulation is increased by the estimated time for such hop calculated at step 6. The simulation can be stopped due to several stopping criteria. A stopping criterion can be the simulation time or the number of hops. Also positional criteria can be used, like for example checking if a charge carrier is passing through a particular plane in space. If no stopping criterion is respected (step 8), we go back at step 2 and continue the simulation. Periodic boundary conditions can be used in x, y and/or z. In the following of this section, we will explain how the three main parameters in Marcus equation (the transfer integral H_NM, the site energy difference ΔE_NM and the reorganization energy λ_NM) can be calculated with a combination of quantum chemistry and classical mechanics tools.

The site energy difference is the difference in energy of the charge carrier before and after the jump. The energy of a charge carrier can be calculated with several methods. For instance, it can be sampled from a Gaussian distribution with a certain width reflecting the energetical disorder of the system. This technique was initially proposed in the pioneering work by Bässler¹⁷  and it allows to phenomenologically recover the broadening of the density of states (DOS) of a charge carrier in the bulk. This approach allowed to gain great insight into the charge transport mechanisms of organic materials and also to achieve qualitatively the correct mobility temperature dependence. However, to achieve the correct mobility field dependence instead, the correlation of the energy landscape is mandatory as demonstrated by Novikov and Vannikov^42,43  and subsequently showed in a KMC scheme.⁴⁴ 

To model a realistic correlation of the energy landscape in the bulk, it is possible to consider in a more sophisticated way the energy of a charge placed on the molecular orbital M as composed by four terms

ϵ_M is the energy of the orbital M=(i, m) and the±sign depends on the type of charge carrier, if the energy of the orbital has to be added (electron) or subtracted (hole). This can be evaluated either at the semi-empirical or quantum mechanical level for each and every molecule in the box in order to consider the energetic shift due to specific conformation. The other terms represent the effect of the environment and they are depending on the positions and the quantum states of the molecules around. The array s̲=(s₁, s₂,…, s_i, s_j,…) represents the state of the system and it is composed by the quantum states of all the molecules. After every charge hopping, the state of the system will change. If for example an electron jumps from the LUMO (Lowest Unoccupied Molecular Orbital) of molecule i to the LUMO of molecule j, the state passes from s̲_1̲=(0, 0, …, LUMO, 0,…) before the jump, to s̲_2̲=(0, 0,…, 0, LUMO,…) after the jump. The term E$iperm-el$(s̲) considers the Coulombic interactions with permanent atomic charges within a cut-off distance, so that a charge carrier located on molecule i interacts only with the molecules inside its interaction set I(i), i.e. the set composed by the molecules nearer than the cut-off (Fig. 4). The interaction of the charge carrier with the permanent atomic charges of the surrounding molecules can be expressed as

where ϵ₀ is the vacuum dielectric constant, h is a molecule belonging to the interaction set I(i) and a, b are atoms. Every atom a has a position r⃑_a and a permanent atomic charge q_a(s_i) (which depends on the quantum state of molecule i).

For the polarization several approaches have been considered. For instance, it is possible to retrieve polarizabilities from polarizable force fields, in particular from works of van Duijnen and Swart (ref. 45) and the more recent work of Wang, Cieplak et al. (ref. 46) In our approach, we use another cut-off, leading to a different interaction set denoted by I_p(i), see Fig. 4. Atoms belonging to the molecules in this set will be polarized according to Thole model,⁴⁷  using linear atomic polarizabilities. These polarizabilities can for example be calculated²⁸  at the QM level using the LoProp localization method.⁴⁸  The charge carrier interaction with the nearby induced atomic dipoles is computed

where p_b(s̲) represents the induced dipole at the atom b, given that the molecules of the system are in the state s̲. The dipoles inside I_p(i) are calculated depending on the permanent charges of all the molecules inside the interaction set I(i). To avoid polarization border effects, the cut-off chosen for polarization has to be significatively smaller than the cut-off for Coulombic interactions.

The last term of eqn (3) considers the field contribution to the site energies

where E⃑_ext is the external field considered.

Transfer integrals can be evaluated in several different ways. We will present in this subsection a non-exhaustive list of methods to calculate them.

The simplest method to calculate the transfer integral is the exponential decay or distance dependent (DD) transfer integral

where the transfer integral is considered to be exponentially decaying with the distance r_NM between the two molecular orbitals M and N. H₀ and α are two parameters that need to be chosen carefully for the molecules at hand. This method was widely used in the early studies where molecules were represented only by a site. An obvious limitation of this method is of course the fact that the relative orientation of the two molecules does not influence the value of the transfer integral.

It is possible to improve the Distance Dependent model through a weighted Mulliken formula⁴⁹  for carbon atoms²⁷  to take into account the relative orientation of the molecules.

where S_ab is the overlap between two 2pπ orbitals and C$mai$ is the 2pπ expansion coefficient related to the atom a for the orbital level m of molecule i. S_ab can be calculated with a Mulliken formula modified to consider the relative orientation of the two 2pπ atomic orbitals on the atoms a and b.^27,50  This approach is based on the approximation that the transfer integral is proportional to the molecular wavefunctions overlap where only the p electrons contributes to the transfer integral.

This method is valid for a symmetric dimer. In this case the extra charge can be considered roughly equally distributed over the two molecules (having equal conformation) and the transfer integral H_MN is equal to half the energy difference between the energetic levels of molecules i and j. Therefore, once the geometry of both the molecules is known, a single point QM calculation can be performed and the transfer integral can be calculated as

respectively for electron (LUMO) and hole (HOMO) transport. This method is called the energy-splitting-dimer (ESD) method.⁵¹  It has been used widely because of its simplicity and low computational cost. However, the fact that this method breaks for non-symmetric dimers makes it difficult to use in real 3D bulk case studies.

It is of course also possible to calculate the transfer integral directly from its definition, as non-diagonal element of the QM Hamiltonian. For every pair of nearby molecules (nearer than a chosen cut-off), three QM calculations are performed: one for the pair and two for the single molecules. The Fock matrix of the pair is then projected on the single molecular orbitals basis (the basis of the molecular orbitals of the two single molecules) and the transfer integrals |H_NM| can be extracted as the non-diagonal elements of this projected Fock matrix in the single molecular orbitals basis. This is known in literature as the projection method. It was originally formulated for ZINDO⁵²  and it has been subsequently extended⁵³  for more sophisticated QM calculations that do not assume orthogonality of atomic orbitals. In its quite general formulation it considers three QM calculations: for the molecule i, for the molecule j and for the dimer ij (also indicated as d). For each of these calculations we obtain a set of wavefunctions

expanded in the atomic basis set {φ_v} and identified by the quantum numbers of the particular orbital level m. The transfer integral from M=(i, m) to N=(j, n) can be written as

where we have used the completeness of the eigenfunctions of the dimer {|ψ$rd$〉}. If we now assume the dimer system d in the void, we can write

where ϵ$rd$ is the eigenvalue associated to the eigenfunction |ψ$rd$〉 and δ_rs is the Kronecker delta. Substituting eqn (10) and (12) in eqn (11) we arrive at

where F^d is the Fock matrix of the dimer, composed by the energies {ϵ$rd$}. Let us assume there are n_i atomic basis functions {φ$vi$} for molecule i and n_j atomic basis functions {φ$μj$} for molecule j, then there will be n_d=n_i+n_j atomic basis functions {φ$χd$} for the dimer and these basis functions can be rearranged such that the first n_i functions of {φ$χd$} are equal to {φ$vi$}, while the remaining n_j atomic basis functions for the dimer are equal to {φ$μj$}. The matrices П^id and П^jd have dimensions n_i × n_d and n_j × n_d respectively, and they take care of the projection of the eigenfunction of the dimer on the eigenfunctions of the single molecules and vice versa. They can be defined element-wise as

As already pointed out, in all the methods introduced in this section the transfer integrals are evaluated excluding the surrounding molecules. However, several studies showed that the polarization can have a non negligible effect on the transfer integral.^54–59  For instance, Castet et al. showed by using a valence bond approach that the inclusion of polarization can decrease the value of transfer integrals by 10–15% for a hole in anthracene.⁵⁹ 

The reorganization energy λ is the energy required for all the structural adjustments of the system, after the charge hopping has occurred. It can be decomposed in two parts, an internal and an external reorganization contribution

The internal reorganization energy (λ_int) is related to the internal geometry rearrangement of the molecule, after a charge is added or subtracted to the molecule. Let us consider a charge hopping from molecule i to molecule j (Fig. 5). The internal reorganization energy is the sum of the reorganization energy of the molecule i that is passing from charged to neutral state and the reorganization energy of the molecule j that is passing from neutral to charged state

where

E$inC$, for example, is the energy of molecule i in the geometrical conformation of the charged molecule (C), but in the neutral electronic configuration (n). All the other terms of eqn (17) can be interpreted analogously. These energies can be calculated by performing quantum mechanical geometry optimizations for the neutral and charged molecules followed by two single point calculations.^25,51 

The external reorganization energy (λ_out) is related to the rearrangement of the environment after the charge carrier's jump. It can be calculated for example from the electric displacement fields created by the charge transfer complex²⁶ 

ϵ₀ is the permittivity of free space, D⃑_I(r⃑) and D⃑_F(r⃑) are the electric displacement fields created by the charge transfer complex in the initial (charge on molecule i) and final (charge transferred to molecule j) states, respectively. V^out is the volume outside the complex, and

is the Pekar factor, determined by the low (ϵ_s) and high (ϵ_opt) frequency dielectric permittivities.

Another possibility to calculate the external reorganization energy is by using atomic induced dipoles and computing the energy needed for a re-polarization of the environment.²⁸  In these calculations we neglect the rearrangement of the permanent nearby multipoles, i.e. we do not perform molecular dynamics during the KMC simulation. The MD structural rearrangement is much slower than both the wavefunctions polarization and the internal reorganization of the molecule (fast reorganization of covalent bonds), thus it is happening on a different timescale respect to charge conduction. The polarization contribution is taken into account by means of induced dipoles on the nearby molecules. Let us consider the energy contribution

where as above a and b are atoms, orbital M belongs to molecule i, and h is a molecule in the polarization interaction set of i. This is the interaction energy between the charges on i as given by the quantum state s̲_1̲ and the nearby dipoles induced by the quantum state s̲_2̲. When a jump from M=(i, m) to N=(j, n) occurs, the quantum states are modified from s̲_1̲ to s̲_2̲, i.e. the state of i changes from charged to uncharged and the state of j changes from uncharged to charged. An illustration of the states s̲_1̲ and s̲_2̲ with their respective induced dipoles are presented in Fig. 6. The external reorganization energy can be expressed as the difference between the energy of the new charge configuration in the old dipoles arrangement (the jump happened but the structure did not relax yet) and the energy of the new charge configuration in the new dipoles arrangement (system has re-polarized through dipoles reorganization)

The energetic contribution between the electric field and the induced dipoles is not considered. For more details on this method of calculating the external reorganization energy we refer the reader to ref. 28.

The CT state splitting into free charge carriers is one of the main processes determining the efficiency of an organic solar cell. Kinetic Monte Carlo simulations can be a very useful tool for the study of such a system. Nevertheless, even once all the material parameters have been calculated, the electric field has to be given as input of the simulation. An iterative process based on trial and error is needed to establish which are the fields at which the CT state splits. This process of determination of the meaningful set of electric fields for the specific case study is slow to perform since many simulations are needed at each field value to obtain a statistically significant ensemble to analyze. Here we describe a method allowing to analyze the qualitative behaviour of a bilayer interface, based on the parameters calculated for Marcus formula, but without having to run KMC simulations. It is therefore possible to determine a priori the interesting sets of fields for the system under consideration and if the two chosen donor and acceptor materials are suitable to be used together in an OPV device.

The hopping probabilities used in the KMC simulation are calculated with Marcus formula (eqn (2)). As detailed in Section 3 both the transfer integrals and the internal reorganization energies are calculated without any external field or extra charges. Theoretical studies support the independence of transfer integrals and internal reorganization energy from the external electric field.^36,60  Also the external reorganization energy (eqn (21)) is independent of the external field and additional charges, due to the linear polarization method used. The only parameter influenced by the field is thus the energy difference ΔE, at least in a first approximation.

Since we are interested in the role of the electric field E⃑ on the hopping rates, we study the Marcus rate w_MN (eqn (2)) while everything is kept constant except ΔE_NM, the only parameter affected by the electric field. Marcus rate w_MN, reported here for the reader's convenience,

or even

The symbol ≏ assumes thus a precise meaning in this context, it is used to specify that the first member of the equation belongs to the interval identified by the second member ±3σ_NM. In this way we will obtain shorter formulas, but the reader should keep in mind that the standard deviation σ_NM is implicit in the notation used and it depends on the particular hopping considered. The fields maintaining the energy difference inside the interval specified by eqn (23) favor conduction in Marcus theory.

Let us now consider an electron and a hole forming a CT state at a flat interface in the xy plane, with the electric field applied in the −z direction (Fig. 7). We focus first on the movement of the electron, keeping the hole still. Two interesting situations can be analyzed: when the electron is in its initial position at the interface and when the electron is far from the interface. When the electron is at the interface, the difference in energy for the electron hop from the first to the second layer of molecules in the opposite direction of the field is

where can be regarded as the contribution of the electron–interface interactions in passing from the first to the second layer (Fig. 7), and the 0 at the subscript specifies that we are considering the ground state. is composed by the shift of the energy levels of the molecules due to their geometrical deformation at the interface (Δϵ), and by the electrostatic interaction of the interface with the charge carrier due to permanent atomic charges (ΔE^perm) and induced dipoles (ΔE^ind), see terms in eqn (3). is the change in energy due to the Coulombic interaction of the hole and the last term in eqn (24) is the contribution of the field, with Δz$FSA$ being the distance between first and second layer of the acceptor in the z direction (that of course does not depend on which molecular energy levels we are considering). Using eqn (24) we can write the Marcus rate for an electron passing from first to second layer as a function of the external electric field

This is still a Gaussian but with standard deviation

Eqn (25) represents the efficiency of the process considered (the electron splitting the CT state by passing from first to second layer) as a function of the external electric field. With the symbol defined in eqn (23), this can be written concisely as

When the electron instead is far from the interface, the energy difference of a hop between two layers in the z direction can be simplified as

representing the difference in energy between the ground states of two consecutive layers in the bulk B and B′. The average distance of two layers in the bulk in the z direction Δz$bulkA$ is in general different from Δz$FSA$ due to the deformations at the interface. The interface will indeed usually introduce an extra stress on the molecule leading to an increased compression of the interfacial layers (for example see ref. 61). In this case, the hole is too far to interact with the electron, and the energy shifts due to molecular distortions (Δϵ), permanent charges (ΔE^perm) and induced dipoles (ΔE^ind) will be neglected since they average to 0 in the bulk. The Marcus rate for a free electron moving in the bulk is

whose standard deviation is

For the electron movement then we can finally obtain the pair of equations

these represent the system in the two situations when the electron is near to or far from the interface. Eqn (31) and (32) specifies the conditions for the field to obtain maximum conduction in the z direction, when the electron is near to and far from the interface, respectively. Eqn (31) is related to the efficiency of the CT state splitting initiated by the electron, while eqn (32) is related to the transport of the free electron to the electrode. If these two equations can be satisfied simultaneously, we will obtain both CT state dissociation (initiated by the electron) and free electron transport in the acceptor phase. Figure 8 illustrates these facts in what we will call the efficiency diagram for the electron. The 3σ tolerances are illustrated both on the energy and on the field axis.

Analogous considerations can be made for the hole, leading to other two equations regarding the efficient CT state splitting initiated by the hole and the transport of the free hole to the electrode

where O and T represents the first and second layer for the hole (Fig. 7). Due to the fact that h is present in the superscript, now the subscript 0 in O₀ and T₀ refers to the ground state of the hole instead. The conditions imposed by eqn (33–34) can be illustrated in the efficiency diagram for the hole, Fig. 9. To have both an efficient CT state splitting and free charge carriers collection at the electrodes the field should respect eqn (31–34). This results in a combination of Figs. 8 and 9 in what will be called a CT state splitting diagram, Fig. 10. The CT state splitting diagram plots together the electric field strength, the energy of the field gained in jumping to the next layer and the efficiencies of the 4 processes of interest in eqn (31–34). From the CT state splitting diagram, it is possible to tell, before performing the KMC simulations, if two materials are suitable to be used together in an organic solar cell and which electric field strengths are required for the simulations to work. Also we can spot immediately from the hypothetical CT state splitting digram of Fig. 10, that the hole is more mobile than the electron. Ideally we would like to have the 4 Gaussians of a CT state splitting diagram to overlap as much as possible. This means that there is an electric field at which the 4 processes are efficient. Suboptimal situations are also possible, leading sometimes to the same efficient CT state splitting. When reading a CT state splitting diagram indeed it is mandatory to keep in mind that the two Gaussians corresponding to the splitting of the CT state (first to second layer) somehow represent the state of the system in the initial situation. We can think that the state of the system at every instant of time is represented by two Gaussians (since only two charges are present). As soon as the system starts to evolve, if the CT state splits and one of the two charge carriers starts moving, the two Gaussians representing the state of the system will be in an intermediate situation, between charge at the interface and charge in the bulk. In the example of Fig. 10 then, the field that is expected to work best for both CT state splitting and free charge carrier conduction is around 0.5×10⁷ V cm⁻¹. At this field indeed, the CT state splitting will always be initiated by the hole, but as soon as the hole moves away from the interface, the Gaussian of the electron will approach the Gaussian of the free conduction and also the electron will start moving efficiently.

The main question that can arise now is which are the factors that keep the Gaussians separated. Let us first focus on the two Gaussians of the same charge carrier, for this task we have to look at eqn (31) and (32). Usually the difference between the Δz near to and far from the interface is not predominant. Also, if considering donor and acceptor materials with similar polarizabilities, we can argue that the external reorganization energies do not change so much at the interface or in the bulk. Under these conditions the main factor shifting the Gaussian of the free conduction with respect to the Gaussian of the CT state splitting is the energy difference in passing from first to second layer without field contribution (noF), i.e.  and for electron and hole, respectively

In the initial CT state, the charge carriers will feel two contrasting forces, for example the electron will feel the attraction of the hole and the contribution of the interface . To have the two Gaussians as much overlapped as possible the interface should have a destabilization contribution balancing the hole attraction (Fig. 11)

This agrees somehow with the intuition that the interface should destabilize “sufficiently” the CT and thus helps the free charge carriers formation. Eqn (37) clarifies how much is “sufficiently”. Polarization helps for two reasons. It introduces a screening of the Coulombic attraction between electron and hole, and increases the destabilization of the interface, thus reducing . Secondly it increases the reorganization energy, adding the external reorganization contribution and thus increasing the Gaussians standard deviations and σ$bulke$. The reasoning for the hole is analogous, leading to the condition

Another question is if the two materials selected are compatible. Assuming that either one of eqn (37) or (38) are respected, we have established that it is possible to split the CT state and conduct at least one of the free charge carriers. To assure that both the charge carriers can be collected we have to analyze the factors that favor the overlap between the free electron conduction and the free hole conduction Gaussians. For this task we have to compare eqn (32) and (34). It turns out that we would like the ratio between reorganization energy and layer spacing in the direction of the field to be quite similar, meaning

If this condition is respected, it exists a field at which both electrons and holes can move. In combination with the fact that eqn (37) or (38) is respected we have just found two materials for which it is possible to find a field at which the CT state splits and both the free charge carriers can be collected at the electrodes.

The analysis done so far focuses on the ground CT state, in which the electron is always on the LUMO of an acceptor molecule and the hole is always on the HOMO of a donor molecule. There is a quite large consensus in the literature on the importance of the role of hot charge carriers for an efficient functioning of organic solar cells.^62–65  This can also be rationalized with the CT state splitting diagram just presented. Introducing excited levels for the electron and the holes introduces new Gaussians in the diagram and thus creates more possibilities for overlaps and complex path mixing of decay and transport. We can argue that in the bulk the main transport is given by the ground levels, thus we can assume that the free conduction Gaussians will not change (at least not significantly). The main effect of the excited molecular levels is expected at the interface in the hot CT states. The transfer integrals H_FmSn have an oscillatory behaviour but are expected to increase with increasing m and n due to the increasing delocalization of the molecular orbitals. This will result in Gaussians with higher amplitude. Neglecting the change of the reorganization contribution (λ_FmSn), we focus again on the energy difference (), that is the main factor determining the overlap. We had seen in eqn (37) and (38) which are the conditions for an overlap between the Gaussians of the charge carrier in the ground CT state and the free conduction Gaussian. Analogously, the condition for the overlap between the Gaussians of the charge carriers in the hot CT state and the free conduction is

For different values of , we will have hot CT Gaussians centered at different positions on the field axis, determining the possibility of the overlap (analogous to the ground CT case of Fig. 11). It will thus be possible to find an energy difference for which the CT splitting Gaussian from F_m to S_n is overlapping with the free conduction Gaussian in the bulk, or with another CT splitting Gaussian that in turn overlap with the free conduction Gaussian in the bulk. In this way it will be possible to establish, at a particular field, which is the most probable chain of events leading to the free charge carriers conduction. If it is possible to find one.

Which chain of events (the system passing from one hot CT state to another) is more favourable depends on the specific case. The definition of is analogous to the ground state case (eqn (35))

i.e. the average energy difference in passing from the m orbital level in the first layer F to the n orbital level in the second layer S. Higher level wavefunctions are more delocalized and can produce a different momentum in the charge distribution on the molecule (dipoles, quadrupoles, etc.). These represent higher order corrections to the electrostatic energies and we can argue that in first approximation the electron–hole interactions do not vary significantly from the ground case, i.e.. Similar reasoning can be done for the interaction of the charge carrier with the surrounding in term of permanent charges (ΔE^perm) and induced dipoles (ΔE^ind). Consequently, the contribution of the interface will mainly differ caused by the orbital energy difference Δϵ. In first approximation the differences in energy without the field can be written as

An analogous argument can justify somehow the neglecting of the change in reorganization energy, leading to and to the simplified set of equations for the hot CT states

Selecting appropriately the initial and final levels of the hop m and n, it is possible to obtain and/or a respecting the eqn (45–46) for a given electric field. The analysis done so far is focused on the hopping rates of the electron or the hole composing the hot CT. If trying to determine the most probable path for the evolution of the system, also the pure decays of a single charge carrier from higher to lower level on the same layer have to be considered. For example for electron on the first layer, decaying rates F_m→F_p with p<m will compete with hopping rates F_m→S_n. Eqn (45–46) have been obtained neglecting higher order terms that could lead to slightly different results, but the message of the equations is still valid: higher levels help the CT state splitting and the overlap of the CT splitting Gaussians with the free conduction Gaussians. Also other paths can be thought, complicating further the CT state splitting diagram, for example considering intermediate steps in the CT state dissociation, i.e. F_m→S_n followed by S_m→R_p, where R represents here a consecutive layer.

We have presented in this book chapter the steps required to study charge transport in organic solar cells. Starting from a realistic morphology using Ising model and MD simulations, Kinetic Monte Carlo simulations based on Marcus formula can thus be performed on the obtained geometry. The parameters entering in the Marcus formula can be calculated with a combination of semi-empirical, QM and classical methods. A detailed analysis of those parameters allows the construction of the CT state splitting diagram, to determine the qualitative behavior of the materials and the electric fields which allow simultaneously the CT state splitting and the conduction of the free charge carriers in the bulk. This allows to select carefully the electric fields at which the KMC simulations should be performed. With the progress of supercomputing, it is becoming more and more feasible to screen a large variety of donor and acceptor systems and interpenetrating networks. However, the analysis of the KMC simulations is not always trivial, especially for non planar interfaces (e.g. BHJ). For these reasons we are currently collaborating with the Visualization center in Norrköping to develop a visualization tool. Figure 12 shows a set of 100 trajectories, where trajectories for the holes in the donor part and for the electrons in the acceptor part are colored in light gray [red] and dark gray [blue], respectively. [colors online] The spheres represent the end points of the trajectories and their colors indicate how long it took for the charge carrier, initially in the CT state, to arrive at this point (light for short time and dark for long time). It is also possible to highlight a specific trajectory in the different visualization panel in order to identify preferential paths in the BHJ morphology paving the way to establishing structure–property relationship.

M. L. thanks SERC (Swedish e-Science Research Center) for funding and SNIC (Swedish National Infrastructure for Computing) for providing computer resources. The authors would like to thank Sathish Kottravel and Ingrid Hotz from the Visualization center in Norrköping for a fruitful collaboration within the SeRC.