Effects of Delocalized Charge Carriers in Organic Solar Cells: Predicting Nanoscale Device Performance from Morphology

Monte Carlo simulations of charge transport in organic solar cells are performed for ideal and isotropic bulk heterojunction morphologies while altering the delocalization length of charge carriers. Previous device simulations have either treated carriers as point charges or with a highly delocalized mean‐field treatment. This new model of charge delocalization leads to weakening of Coulomb interactions and more realistic predicted current and fill factors at moderate delocalization, relative to point charges. It is found that charge delocalization leads to significantly increased likelihood of escaping interface traps. In isotopic two‐phase morphologies, increasing the domain sizes leads to slight decreases in predicted device efficiencies. It was previously shown that tortuous pathways in systems with small domain sizes can decrease device performance in thin film systems. However, the diminishing effects of Coulomb interactions with delocalization and efficient separations of excitons by small domains make morphological effects less pronounced. The importance of delocalization, which has largely been ignored in past simulations, as a parameter to consider and optimize when choosing materials for organic solar cells is emphasized.


Introduction
A strong theoretical framework of charge transport in organic materials [ 1,2 ] is needed to improve the effi ciencies of organic fi eld effect transistors (OFETs) [3][4][5][6][7] organic light emitting diodes (OLEDs) [ 8,9 ] and organic photovoltaics (OPVs). [10][11][12][13][14][15] The promise is that these devices will be cheaper and more environmentally friendly than their inorganic counterparts. Costs are may be low because devices can be manufactured using existing solution processing techniques [16][17][18] and synthetically tailored [ 19 ] to specifi c applications. However, the understanding of charge transport in disordered organic semiconductors is far from complete. Mechanisms of charge transport in organic materials are very different from crystalline inorganic materials, where band transport occurs. While band transport can be observed in pure molecular crystals at low temperature, charge

Effects of Delocalized Charge Carriers in Organic Solar Cells: Predicting Nanoscale Device Performance from Morphology
Adam G. Gagorik , Jacob W. Mohin , Tomasz Kowalewski , and Geoffrey R. Hutchison * Monte Carlo simulations of charge transport in organic solar cells are performed for ideal and isotropic bulk heterojunction morphologies while altering the delocalization length of charge carriers. Previous device simulations have either treated carriers as point charges or with a highly delocalized mean-fi eld treatment. This new model of charge delocalization leads to weakening of Coulomb interactions and more realistic predicted current and fi ll factors at moderate delocalization, relative to point charges. It is found that charge delocalization leads to signifi cantly increased likelihood of escaping interface traps. In isotopic two-phase morphologies, increasing the domain sizes leads to slight decreases in predicted device effi ciencies. It was previously shown that tortuous pathways in systems with small domain sizes can decrease device performance in thin fi lm systems. However, the diminishing effects of Coulomb interactions with delocalization and effi cient separations of excitons by small domains make morphological effects less pronounced. The importance of delocalization, which has largely been ignored in past simulations, as a parameter to consider and optimize when choosing materials for organic solar cells is emphasized. transport is typically described by localized carrier hopping at room temperature. [ 20 ] The charge carriers are localized to molecular sites through disorder of the morphology and electron-phonon coupling. In contrast, highly delocalized carriers move through the energy bands of an inorganic material, yielding a high charge carrier mobility.
More effi cient charge delocalization is therefore a major goal of molecular design in conjugated materials. Although the carriers are localized, some degree of delocalization can still play a critical role, especially in OPVs, where Coulomb bound interface traps easily form between "separated" carriers. For example, charge delocalization between π-stacked benzenes has been observed . [ 21 ] Delocalization has also been suggested as a mechanism for effi cient long-range charge separation in organic semiconductors. [ 22,23 ] For example, well-ordered, semicrystalline domains likely lead to charge delocalization and improved device performance in OPVs. Single-junction devices with well ordered morphologies have been reported with very large fi ll factors (76-80%) and power conversion effi ciencies of up to 8.7%. [ 24 ] A complete picture of the role of crystallinity is not yet clear. While charge transport is enhanced by crystalline domains, exciton transport may be inhibited in well ordered systems. [ 25 ] This additional level of complexity is due to the dominance of exciton diffusion occurring by an intermolecular mechanism. In fact, exciton delocalization in regioregular P3HT is estimated to be as low as 1 to 2 nm, and therefore may not play a major role in charge separation. [ 26 ] Studying the interplay of delocalization with electrostatic interactions between carriers may help elucidate these effects.
The combination of electrostatics and delocalization can play a major role in the charge separation process. For example, along with hot charge transfer (CT) states, calculations show that charge delocalization can lead to ultra fast (<100 fs) charge separation via the lowering of Coulomb barriers. [ 27 ] While delocalization can be treated very accurately in quantum calculations [26][27][28] the effect is generally ignored in meso-scale simulations of charge transport. In such simulations, highly accurate treament of electrostatic effects is achieved by largely ignoring the effects of delocalization and treating carriers as point charges [29][30][31][32] Since previous simulations have treated carriers exclusively as point charges or not at all, our primary goal was to add delocalization to Monte Carlo charge transport methods to more accurately represent carrier electrostatic interactions. To address these questions, we have incorporated delocalization into our existing Monte Carlo model [29][30][31][32] of charge transport in OPVs. We have used the Monte Carlo model to examine a series of four isotropic two-phase morphologies with increasing domain size, three ideal morphologies based on mathematical minimal surfaces, and morphologies based on bands. We have not yet considered how delocalization effects transport for every possible mixed phase morphology. Such effects are likely worth further investigating.
In this work, the effects of delocalization have been examined by varying the delocalization length of the carriers, from 0.75 to 1.50 nm, which are compared to simulated point charges and non-interacting charges. These length scales are comparable to those found from quantum calculations. [26][27][28] As a result, transport behavior is affected on the length scale of the exciton dissociation radius (≈20-30 nm). We have not yet addressed various questions such as how delocalization may be different in the donor and acceptor phases, anisotropy of delocalization, or how delocalization may vary during transport. For example, upon initial absorption of a photon, an exciton may be highly delocalized, and then localize after a very short time period (a few ps) due to electron-phonon coupling. Such topics are clearly relevant for future simulation efforts.
We have used our Monte Carlo simulations to predict current-voltage curves (IV Curves) of a solar cell device. From the IV curves we compute the effi ciency of each system by calculating fi ll factors, a common fi gure of merit for solar cells. The fi ll factor (FF) is the ratio of the maximum power that can be extracted from the cell to the theoretical power, as shown in Equation 1 : Here, i mp is the current at maximum power, v mp is the voltage at maximum power, i sc is the short circuit current, and v oc is the open circuit voltage. The fi ll factor ranges between 0 and 100%, and can be viewed as the ratio of the areas of two rectangles-a smaller one formed by i mp and v mp , and a large one formed by i sc and v oc . The closer the fi ll factor is to 100%, the more square (ideal) the IV curve, and the more effi cient the performance of the solar cell. In a real device, the solar cell is connected to an external load, which has a resistance, R . Only when R is the ratio v mp / i mp can the maximum power be extracted from the device. The short circuit current i sc corresponds to R = 0, while the open circuit voltage corresponds to R = ∞. While the fi ll factor gives information on IV curve shape, solar cells with similar fi ll factors can have different p mp . For example, a device that poorly passes current can have a small i sc and i mp , yet appear to have a good fi ll factor. Therefore, it is important to consider multiple measures of performance, as discussed below.

Lattice
The model consists of a lattice of cubic sites, with each site representing neutral molecules, 1.0 nm 3 in size, approximately the same size as an isolated fullerene molecule. This is therefore the smallest possible lattice without treatment of atomicscale details.
In this work, we have used a lattice of 192 nm × 192 nm × 32 nm. Charge carriers, such as holes and electrons, can occupy the sites. Holes are to be viewed as an electron missing from the highest occupied molecular orbital (HOMO) of a molecule. Electrons are to be viewed as an excess electron in the lowest unoccupied molecular orbital (LUMO). When an electron and hole occupy the same site, they are an exciton. Excitons can recombine to form unoccupied sites.

Exciton Injection
The simulation begins by injecting excitons randomly in the system, as a model for photo-generated electronic excitations. We have assumed that the system is ideal in that all sites have the same probability to form an exciton, regardless of thickness. In the future, we will relax these assumptions, to allow for non-uniform thickness or self-adsorption in deeper layers. For now, we have used an injection probability estimated from the AM1.5 solar spectrum, and is proportional to the surface area of the device. For a 256 nm × 256 nm area, we have estimated an injection probability of − 10 3 ps −1 . This injection probability is scaled by a factor of 192 /256 0.5625 for an area of 192 nm × 192 nm. While optical effects may change the intensity distribution (e.g., due to absorbance attenuating light in the fi lm, and differences in refractive indices between the p-type and n-type phase) these are currently ignored.

Recombination
Recombination is also allowed to occur at a fi xed probability for excitons in the system. Unlike injections, a fi xed recombination probability does not directly correspond to a fi xed rate because excitons are only considered when holes and electrons occupy the same site. A constant probability is a simplifi cation; recombination in organic solar cells can happen by many mechanisms with different timescales. [ 33 ] For example, in a bimolecular recombination mechanism, the Langevin expression predicts recombination rate to depend upon the hole/ electron mobilities, the electron/hole concentrations, and the intrinsic carrier concentration. [ 33 ] We have used a fi xed recombination probability of − 10 4 ps , inspired by kinetic Monte Carlo simulations. [34][35][36] Recombination is performed at the end of each step. All electrons and holes occupying the same site are allowed to recombine at the fi xed rate given above. The units on − 10 4 ps −1 mean that, if, on average, only one exciton is present at the end of every Monte Carlo step, then it would take 10 ns for an exciton to recombine. There are frequently more than one exciton at the end of each step. In practice, this fi xed rate is not realized because the average number of excitons present in the simulation is highly dependent upon charge concentration. The concentration will change with voltage. This means charges are less likely to combine at high voltages (carriers are being pulled apart) and more likely to recombine at v oc . At v oc , the external bias balances the internal bias. With only Coulomb interactions present, pulling separated carriers apart becomes diffi cult; it is likely for excitons to reform and recombine.

Energy Evaluations
Charges hop between adjacent sites using the Metropolis criterion. The probability to hop between sites is , independent of hopping direction. The value used determines the overall magnitude of the current calculated, and approximates experimental current in OFETs. [ 32 ] Directional independence is an approximation that may not be representative of many systems, where charge transfer rates may be anisotropic (i.e., along π-π stacks). For example, charge transfer effi ciency will be preferred along stacked domains, although in an amorphous device, such directional effects would be expected to average out. Additionally, coupling would be different among differing materials. H AB is made exponentially smaller, , for a two-site hop. This is made to simulate an exponentially less likely tunneling event in variable-range hopping. Carriers are considered individually for trial moves. Sites are assigned energies representing the HOMO and LUMO energies of the molecules. There are three main contributions: the donor/acceptor offset, the overall device potential between electrodes (including differences between metal work functions), and the Coulomb interaction between charges.
The sites are classifi ed as donor or acceptor, according to the system morphology chosen, as discussed below. The donor/ acceptor LUMO and HOMO are offset by a constant energy ΔE LL = ΔE HH = 0.5 eV. This value is slightly larger than the Coulomb interaction between two opposite charges spaced 1 nm apart 0.411 eV. It provides a driving force to separate carriers at the interface. and the open-circuit voltage v oc is idealized at an optimum of ≈1.5 V, to consider the possible effi ciencies of ideal materials. Realistic experimental v oc values are usually 0.5-0.8 V. Next, a linear electrostatic potential is applied between the electrodes, using the potential difference Δ = − V V V f i and width of the device, L . This potential includes the external and intrinsic bias. To construct an IV curve, this potential is varied from −2.0 to 2.0 V to pass the v oc .

Electrostatic Interactions
Disorder is added to the system using the Coulomb interaction, as shown in Equation 2 . In Equation 2 , N is the number of charges (carriers or charged defects), q is the charge, r ij is the distance between charges i and j , and = 3.5 ε is the dielectric constant. For every Monte Carlo move, the Coulomb interaction must be calculated for the initial and fi nal state -around N N 2 steps × × times. Even with a cutoff of 50 nm and lookup table for inter-cell distances, the brute force evaluation of Equation 2 is the most time consuming step of the simulation, because it scales as N 2 . Therefore, a parallel GPU implementation of the Coulomb sum was used to make the calculation feasible. With the GPU code, carrier numbers of 10 000 or more can be handled.
The error function only modifi es the Coulomb interaction at short range. For example, when σ = 1 nm, at ≥ r 2 nm. However, at short distances, the Coulomb potential is signifi cantly diminished. For example, when σ = 1 nm, at ≤ r 1 nm. This effect can be seen in Figure 1 A, where the potential at each value of σ converges to the Coulomb potential at large r (i.e., ≥ 3 nm). The charge density as a function of r is shown in Figure 1 B. The charge density was chosen such that each Gaussian integrates to a unit charge. Note that this solution to Poisson's equation can only be obtained in 3D, so could not be applied to our previous simulations [ 30 ] on 2D monolayer OPVs. The use of spherically symmetric charge distributions is an oversimplifi cation. In reality, molecular orbitals with angular momentum are not spheres. Anisotropy will be added in the future, for example, by considering Gaussians ellipsoids.
The choice of the σ delocalization lengths of 0.75 nm to 1.5 nm, as discussed above, corresponds to a range found through quantum mechanical calcualtions. [26][27][28] Moreover, as seen in Figure 1 A, with the 1.0 nm grid spacing between molecules, σ ranging from 0.75 nm, 1.0 nm, and 1.25 nm, roughly span the energetic range between point charges and σ = 1.5 nm. Smaller σ values experience almost the same electrostatic interactions at this grid spacing as point charges, and vanishingly small probabilities of acceptable Monte Carlo moves, as indicated in Table S1 and illustrated in Figure S5 (Supporting Information). Since the spatial resolution of our Monte Carlo grid roughly corresponds to a molecular size, while decreased delocalization lengths may be interesting, such considerations would require treatment of atomic-level (i.e., quantum) details.

Systems
Morphologies were constructed with a domain of 192 nm × 192 nm × 32 nm, or 1,179,648 sites. Periodic boundary conditions were not used, to properly include fi nite size effects in real materials, such as surface trapping. All morphologies consisted of a donor and an acceptor phase whose frontier molecular orbital energies were offset arbitrarily by 0.5 eV, as described above. Device performance can be affected by many parameters, such as mobility and connectivity, and is not limited to just charge delocalization. Therefore, we have considered ideally connected morphologies and morphologies more representative of bulk heterojunctions. Realistic, isotropic morphologies were simulated, and compared with "ideal" morphologies constructed from equations of triply periodic minimal surfaces, [ 37 ] shown in In Equation 5 -7 , k is an angular wavenumber, defi ning the periodicity of the surface. We choose a k that corresponds to a wavelength λ = 8 nm, as shown in Equation 8 , where L is a dimension of the grid, and n is the number of periods along the dimension. Only systems where ( = = = k k k k x y z ) were studied, and the minimum possible distance between surface boundaries is λ 2 = 4 nm.
In practice, sites were assigned by evaluating f x y z ( , , ) on the domain and thresholding. Donor sites were placed where ≥ f x y z ( , , ) 0, while acceptor sites were placed where < f x y z ( , , ) 0. This assignment makes the zero-level contour-isosurface of f x y z ( , , ) the boundary between donor and acceptor sites. The last ideal morphology studied was the "band" morphology. The band morphology is a checkerboard pattern of alternating donor and acceptor sites in the yz-plane, extended along x -direction. Two band morphologies were used. The band4 morphology has a checker size of 4 × 4 nm 2 . The band8 morphology has a checker size of 8 nm × 8 nm. Figure 2 D shows a (32) 3 nm 3 slice of the band4 morphology.
Traditional isotropic morphologies, intended to simulate a random organic bulk heterojunction, are also shown in Figure 2 . Such isotropic morphologies were constructed using an algorithm discussed previously, [ 30 ] extended into three dimensions. In the algorithm, a grid of random noise is constructed in the domain and then convoluted with a threedimensional Gaussian of desired width, σ. These real valued sets are then converted into binary maps by a thresholding rule, with zeros corresponding to D (donor) and ones corresponding to A (acceptor) phases. The average domain size in such maps is set by the value of σ. To produce morphologies with equal volume fractions of both phases ( φ φ = =50 1 2 %), the thresholding rule was applied at the 50% level of the cumulated probability of the normalized distribution of intensities. In these morphologies, the domain size 〈 〉 L p is approximately twice the size of σ , such that the σ = 3, 4, 6, 9 nm leads to domain sizes of 6, 8, 12, 18 nm. Four isotropic systems were studied, named iso3 , iso4 , iso6 , and iso9 . An important metric of the morphologies is the inter-facial area between phases, which is correspondingly 196.04, 123.31, 104.39, 71.52 m 2 cm −3 . The iso3 morphology domain size is comparable to the minimum curvature morphologies' minimum distance between phases of 4 nm. This excludes the band8 morphology, which has a domain size more comparable to the iso4 system.
Although the isotropic two-phase morphologies represent a mixed phase system already, we probed the effect of adding inter-dispersed acceptor sites in the donor material of the iso9 system. Such "mixed phases" occur when individual fullerene molecules disperse into the polymer donors. In this case, called iso9p , 1179 sites (0.2%) of the donor phase was changed to acceptor material at randomly chosen locations. In the end, however, one must keep in mind that a real macroscopic device Adv. Funct. Mater. 2015, 25, 1996-2003 www.afm-journal.de www.MaterialsViews.com will be a combination of multiple regions with somewhat different morphologies, and thus should view our simulation as a subset of such a device.
While our model could consider more phases, that is not our current focus. In the future, a few avenues will be pursued to better differentiate phases. For example, multiple energetic offsets will be used to create a 3 or more domain system. Anisotropic transfer couplings (H AB , discussed above) will be used to consider multi-scale size effects, such as crystalline regions with larger transfer rates. While we use a rectangular lattice, others will be implemented, to better probe the effects of crystallization. Though this would modify transfer probability through changing the energetic calculation, the more direct approach would be changing H AB . Finally, we will consider a multi-region cell, where different regions contain differing order parameters, such as domain size.

Figure 3
A shows the computed fi ll factors for ideal morphologies. Aside from the band8 morphology, each ideal system has the same domain size of 4 nm. As the delocalization is increased from point charges ( σ = 0.00 nm) to delocalized charges ( σ = 0.75 nm), the average FF increases from 35% to 58%. The FF continues to increase, less dramatically, as the delocalization is increased to σ = 1.50 nm. Non-interacting charges, which can be viewed as infi nitely delocalized carriers in this model, yield devices with the best FF (70%). While there is no signifi cant trend in morphology with delocalized charges, the trend with point charges is gyroid < d-surface < p-surface < band4. Small differences exist, but are not statistically signifi cant with delocalized charges. Figure 3 B shows the maximum power ( p mp ) for ideal morphologies. The p mp indicates the maximum work per second that can be performed by the solar cell. The maximum power and power conversion effi ciency are both proportional to the fi ll factor, and thus increase with delocalization. The largest increases are seen when moving from point charges to delocalized charges, and from delocalized charges to non-interacting charges. With point charges, transport is hindered by energetic disorder, but as the energetic disorder is smoothed out by delocalization, the maximum power increases because more carriers can navigate to the electrodes. Interestingly, the predicted i sc is consistent across all ideal morphologies and delocalization lengths, but the i mp drives the increase in maximum power. This suggests that delocalization effects change the shape of the IV curve, not simply increasing the overall current due to an increased number of carriers escaping interfacial pinning. Figure 4 A shows computed fi ll factors for isotropic two phase morphologies as a function of sigma. The iso4 and iso6 Adv. Funct. Mater. 2015, 25, 1996-2003 www.afm-journal.de www.MaterialsViews.com  systems have very similar average domain size (8.1, 9.57 nm). Within error, their FF are very similar. As the delocalization of carriers increases, the device performance increases, much like the ideal morphologies. The most dramatic increase occurs between point charges and σ = 0.75 nm, where the fi ll factors increase from < 40% to > 50%. Since charge carriers in organic solar cells are not point charges, but delocalized, the former is unrealistic for such idealized morphologies, and the latter is more likely considering champion solar cells already achieve FF approaching 70%. However, the differences of FF for σ = 0.75-1.50 nm is much less pronounced. Only when non-interacting charges (σ = ∞) are used, does the FF increase signifi cantly. The fi ll factor decreases in iso9p relative to iso9 for σ = 1.00 nm and non-interacting charges. However, for point charges, the fi ll factor is not smaller. Figure 4 B shows p mp for isotropic morphologies, which increases with delocalization. The p mp is consistently smaller in iso9p relative to iso9 . In the iso9p system, acceptor sites, termed "pepper," were randomly added in the donor phase, representing PCBM inter-dispersed in the P3HT phase. This is due a smaller i mp when pepper defects are present ( v mp does not vary between iso9 and iso9p ). The pepper defects acts as electron traps, reducing the number of carriers that can escape to the electrodes. When we compare the idealized mathematical minimal surfaces, such as gyroid, p-surface, and d-surface to the two-phase isotropic morphologies, such as iso3 , iso4 , and so on, we note that the i sc and i mp are slightly lower. That is, the typical curved two-phase isotropic systems show greater predicted maximum power and thus power conversion effi ciency than idealized systems. We suspect the increased interfacial surface area in these bulk-heterojunction morphologies, combined with efficient exciton separation and escape from the Coulomb-pinned interfaces with delocalized charges, lead to such high maximum powers. These results suggest that in materials with high delocalization and effi cient exciton separation, traditional bulk-heterojunction morphologies may be at least as effi cient as band or "comb" idealized morphologies due to high surface area.
In a previous study, on simulations of 2D monolayer devices, we found that decreasing the domain size leads to poorer charge extraction in thin fi lms. [ 30 ] We argued that while smaller domain sizes lead to effi cient exciton separation, small tortuous pathways can hinder charge extraction. The close proximity of charges and tortuous morphology to the electrodes lead to a disordered energy landscape with a limited number of viable pathways to the electrodes. The trend appears to reverse in larger bulk 3D systems. That is, smaller domains show improved device performance at each value of the delocalization parameter σ. However, we note that the effects are relatively small and within statistical errors. More replicas need to be simulated and averaged to elucidate the effect. Given the interplay between tortuous pathways and effi cient exciton separation, the separation appears to win. That is, we propose that more pathways have opened up in the 3D environment than were present in our previous 2D thin fi lm simulations, due to the increased dimensionality. Moreover, the decreased electrostatic interactions due to the delocalized electrostatic model likely increases the infl uence of effi cient exciton separation. Figure 5 shows snapshots of the Monte Carlo simulation for the iso3 morphology at v oc . At v oc , the potential between the electrodes balances the internal intrinsic bias. Therefore, there is very little incentive for excitons to separate and free carriers to travel to the electrodes. Recombination is the most likely fate of the carriers in this regime. However, the effects of charge delocalization are easily visualized at v oc . When σ = 0.75 nm (Figure 5 A), holes and electrons remain bound in interface traps. However, as the delocalization is increased to σ = 1.50 nm (Figure 5 B), the number of bound interface traps decreases signifi cantly. The energetics of escaping the interface trap are more favorable when the carriers are delocalized.
Consider, for example, an electron-hole pair trapped at an interface. The exciton has separated, putting the electron in Adv. Funct. Mater. 2015, 25, 1996-2003 www.afm-journal.de www.MaterialsViews.com  the LUMO of the acceptor, and the hole in the HOMO of the donor. The Coulomb attraction between the hole and electron in dielectric continuum ( ≈ 3.5 ε ) separated by 1 nm is 0.411 eV, nearly 15.9 times larger than kT at 300 K (0.0259 eV). A potential difference of 1.0 V between electrodes separated by 256 nm only releases 0.00390 eV of energy during the 1.0 nm hop out of this state. Given that the Coulomb attraction is only halved at 2.0 nm, this means the hole and electron will "separate" only a minuscule fraction of the time ( β − Δ e E = 0.0407%). Even if the carriers reach this state, it is extremely likely the charges fall back into the charge trap state, and possibly recombine. There must be other mechanisms to explain the effi cient charge separation in organic materials.
The delocalization of charges can weaken the Coulomb interaction between trapped carriers. While a fully quantum treatment would be ideal, the interaction of spherically symmetric Gaussian charge distributions can be modeled using the standard Coulomb law and the error function (Equation 4 ). Using this model, the probability to separate carriers, of only modest delocalization lengths ( σ = 0.5-1.5 nm) can increase dramatically. For example, a delocalization length of 1.0 nm dampens the Coulomb interaction by erf (1.0 / 2) = 68% at an electron hole separation of 1.0 nm. The probability to escape the charge trap state increases by over a factor of 100 (ΔE = 0.0806 eV, P = 4.42%). These calculations are compiled in Table S1 (Supporting Information).

Conclusion
In summary, our simulations show decreased effect of morphological details and that charge delocalization leads to significantly increased charge extraction in OPVs. We propose that the weakening of energetic barriers, particularly in Coulomb interface traps leads to the effi cient separation of charges. Multiple length scales are involved in assessing device performance, such as micron-scale morphological effects and the nanoscale hole-electron radius of interaction. Delocalization increases the interaction radius, allowing charges to escape the electrostatic pull at the interface. This occurs well under the length scale of the morphology. The escaped charge carriers can more readily exit the system, increasing device parameters such as current, maximum power, and fi ll factor.
The predicted fi ll factors in 3D are overall larger (average FF of 58 ± 12%) and closer to experimental champion devices (FF ≈ 60-70+%) than previous 2D thin fi lm simulations, [ 30 ] where the predicted FF were unrealistically low for idealized devices. The predicted fi ll factors also depend less on the morphology, suggesting new materials with highly delocalized carriers (i.e., along the polymer chain, between chains, and between acceptor molecules) may be less sensitive to morphological differences between devices.
In thin fi lms, the lattice easily becomes frustrated and the Coulomb interactions dominate the energy landscape. In 3D systems, more pathways open up, alleviating some of the frustration, and increasing the device effi ciency. While the differences in morphology (domain size for isotropic systems and geometrical structure for ideal systems) did not play an important role in device effi ciency, we note that a larger sample of morphologies, particularly those with mixed phases, must be examined in the future with more factors to differentiate them.
Considering charges in π-conjugated organic materials are known to be delocalized, this new delocalized electrostatic model is a much more realistic and computationally effi cient method to treat carrier-carrier interactions in organic solar cells and similar devices. Point charge models dramatically overestimate the interfacial charge pinning, and mean-fi eld methods signifi cantly underestimate such effects.
Seeking materials that can effectively allow charges to delocalize over multiple molecules, especially near an interface, is a promising route to create effi cient solar cells. While a material with an overall higher dielectric constant will allow for a more effi cient device (i.e., inorganic materials), it may be that the weakening of Coulomb interactions near an interface is the dominant player in organic materials and should be the focus of future materials design.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.