Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics

This review summarizes the current understanding of electrostatic phenomena in ordered and disordered organic semiconductors, outlines numerical schemes developed for quantitative evaluation of electrostatic and induction contributions to ionization potentials and electron affinities of organic molecules in a solid state, and illustrates two applications of these techniques: interpretation of photoelectron spectroscopy of thin films and energetics of heterointerfaces in organic solar cells.


Introduction
The success of electronic and optoelectronic technologies relies on the possibility to control the energies of transport levels and excitations, and organic semiconductors based on π-conjugated molecules and polymers make no exception [1][2][3]. Owing to their nature of synthetic and soft mat erials, organic semiconductors present several advantages over inorganics, such as their low cost or the possibility to realize flexible devices [4], but also have specific features. The latter mostly arise from the weak magnitude of non-covalent intermolecular interactions, which result in structural disorder, localized charge and energy carriers, accentuating the molecular character of electronic and optical properties [5][6][7][8].
Established theoretical frameworks for the calculation of transport levels and optical excitations in periodic systems, such as density functional theory (DFT) and time-dependent DFT, or more accurate many-body methods constitute an important reference for perfect crystals [9,10], but are of moderate utility for disordered materials at room temperature. A more appropriate and effective approach consists in taking the molecular electronic structure as a starting point, and considering electrostatic interactions between localized charges and neutral excitations with their polarizable molecular environ ment [11][12][13].

Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics
Electrostatic and induction interactions are responsible for shifts of the energy levels of charge carriers in organic semiconductors on the order of eV [14], and hence have a primary impact on charge transport [15], photoionization measurements [16,17], charge injection at metal-organic interfaces [1,18], and charge separation at organic-organic heterojunctions for photovoltaic applications [19][20][21]. Over the last years, an increasing attention has been devoted to modelling these phenomena with often classical but accurate methods that provide the necessary bridge between molecular and materials properties [22][23][24].
This topical review has the twofold objective of summarizing the state of the art in the numerical simulation of long range electrostatic and induction interactions in molecular systems, and of critically discussing the insight brought by this type of approaches in the context of organic semiconductors. A special emphasis is given to the relationship between molecular properties, supramolecular organization, disorder and energy landscape. Where possible, a connection with experiments providing a direct benchmark for the theory will be established, in particular with photoelectron spectroscopy measurements on molecular films. This paper is organized as follows. Section 2 introduces the fundamental quantities and concepts that will be then used throughout the review. The different theoretical approaches describing intermolecular electrostatic interactions in polarizable molecular systems are presented and compared in section 3. Here, we first attempt a critical survey of different methods, from simple induced dipoles schemes to fully quantum approaches, and then focus on the subtle issues arising from the long range nature of electrostatic interactions in extended systems. Section 4 is dedicated to applications, covering systems of increasing complexity, from bulk crystals to films and interfaces, considering both ordered and disordered systems. General conclusions and perspectives for future development are given in section 5.

Fundamental concepts and nomenclature
In this review, we are interested in the determination of the energy of single-particle charged excitations, corresponding to the energy to add or remove one electron to the system, or of electrically neutral two-particle electron-hole excitations. While most of the quantities introduced in this section can be defined for any insulating systems, the background idea of the following development is that of localized charges in molecular organic solids.
The energy levels relevant to hole and electron transport are the ionization potential (IP, or ionization energy) and the electron affinity (EA): where U 0 , U + , U − are the energies of the neutral, positively and negatively charged system. The EA can be also viewed as the IP of a negatively charged system. We point out that IP and EA, which can be experimentally measured with photoelectron spectroscopy, are both intensive quantities defined as the difference of extensive energies and they should be calculated as such. The asymmetry in the definition of IP and EA comes from historical reasons, as both were defined as positive quantities, although the process of removing one electron is always endoergonic, while the reverse one is usually exoergonic in molecular solids.
Other related quantities that can be accessed by experiments and calculations are the transport gap that corresponds to the energy required for creating an electron and a hole at infinite distance, and the energy of an electron-hole (e-h) excitation at distance r: where we introduced the energy of the system with an e-h pair U r ( ) ± and the e-h exciton binding energy E r 0 b ( ) < ± accounting for the Coulomb interaction between electron and hole. We remark that the term charge-transfer (CT) state or excitation is generally used to identify coherent quantum states where electron and hole reside on neighbouring molecules, such as those formed upon direct photoexcitation and hence accessible by optical spectroscopy. Space-separated e-h pairs can be instead detected with electroabsorption measurements.
Intermolecular interactions in condensed phases largely affect the energetics of charge carriers, and the description of these phenomena and their relation to molecular properties and supramolecular organization is a central subject of this review. The contribution of intermolecular interactions to the energy of charged excitations can be quantified as the shift of IP and EA values in a given medium (crystal or amorphous solid, interface, liquid) with respect to gas phase (isolated molecule): The differences P +/− , schematically defined in figure 1(a), are historically called polarization energies because the formation of induced dipoles reacting to a localized extra charge was first recognized as the main source of its stabilization in the solid state, as pioneered by the work by Munn [5,11,12,25] and Silinsh [13,26]. The expression electronic polaron indeed describes a charge carrier dressed by the electronic polarization of the surrounding medium, in analogy with the concept of polaron that is instead related to structural relaxation upon charging. The term polarization energy will be retained throughout this review, although it may sound misleading according to the current understanding. In fact, other effects contribute to the energy of transport levels, the most important of which is the electrostatic interaction between the excess charge and the charge densities of the other molecules in the systems, as discussed in detail in section 3. This electrostatic contribution can have a comparable, and sometimes even higher, magnitude than that of induction. Structural relaxation and charge delocalization do also contribute to a lower extent to P +/− , but unless differently stated these effects will be neglected henceforth.
Turning to heterosystems with electron donating (D) and accepting (A) units, such as those employed for organic photovoltaics or molecular doping, the relevant gap is that involving hole and electron on D and A molecule, respectively. The energy levels scheme for a D/A heterojunction is sketched in figure 1(b), where we introduce the so-called photovoltaic gap: The energy of e-h states with a hole on the D and an electron on the A, Γ, and the corresponding exciton binding energy, are a straightforward generalization of equation (4).

Theoretical methods
We introduce here the different methods for the calculation of the energetics of charged excitations in molecular solids and at their interfaces. We recall that a correct description of long-range electrostatic effects requires performing calculations on large systems, comprising from one hundred to several thousands of molecules, posing a severe requirement that considerably restricts the number of suitable theoretical approaches. A handy reference to the main features of the different theor etical methods discussed in this section is provided in table 1.
The following discussion mostly focuses on methods based on the assumption of zero intermolecular overlap, a common and often implicit approximation that greatly simplifies the description of large systems of mutually interacting molecules. Owing to this approximation, charge carriers are fully localized on molecular units that interact only through classical electrostatic forces, intermolecular covalency and exchange being discarded. A net integer charge can be therefore assigned to individual molecules building up any charge configuration of interest (typically zero, one or a few molecular ions in a system of many neutral molecules) in the spirit of a valence bond (VB) approach. Before proceeding, we remark that overlap must conversely be taken into account for describing electronic bands and charge transport processes. However, the relative magnitude of electrostatic effects (∼1 eV) and intermolecular charge transfer couplings (∼1-100 meV) suggests the use of perturbative approaches where the second effect is introduced a posteriori, for instance via hopping or band models.
The energy of a system of non-overlapping molecules can be written as a functional of the set of charge densities r m { ( )} ρ relative to the different molecules: is the energy of molecule (or molecular ion) m in the electrostatic potential exerted by the other molecules: The second term in the right-hand side of equation (8) prevents from double counting interactions.
can possibly include an external applied potential, usually introduced in the calculation of the dielectric susceptibility [36][37][38]. For the moment we keep general expressions for E m and m ρ , anticipating that these quantities can be made explicit at various  [35] levels of theory. The minimization of the total energy U in equation (8) with respect to r m { ( )} ρ leads to a self-consistent variational problem. Self-consistency arises from mutual intermolecular interactions, i.e. the fact that the potential acting on a given molecule depends on the charge density of the others and determines the secondary polarization field created by the molecule itself.
Because of the weak intermolecular interactions in the solid state it is convenient to rewrite the molecular charge density as is the density of isolated (gas phase) molecule and r m ( ) δρ is the differential charge density induced by the potential , which is supposed to be small. This partitioning leads to a physically sound interpretation of the charge carrier polarization energy, introduced in equations (5) and (6), that corresponds to the contribution of intermolecular interactions to the ionization energy of a molecule in a condensed phase. The polarization energy can be in fact partitioned in two different contributions: The variational problem defined by equations (8) and (9) is often solved with the help of molecular models that provide discrete approximations of r m g ( ) ρ and r m ( ) δρ , that lead to accurate and practical classical or semiclassical schemes that we discuss in the following. A simple and effective approximation for gas phase charge densities r m g ( ) ρ is to describe them as a finite set of point charges, typically located at the nuclear positions: where q mi g is a permanent (i.e. non polarizable) charge at position r mi and r ( ) δ is the Dirac delta function. Partial atomic charges are not unequivocally defined in electronic structure calculations, but the ideal option for our purposes consists in a set of charges optimized in order to reproduce the molecular potential. For that purpose, different electrostatic potential fitting schemes (for instance ESP [39], RESP [40], CHELPG [41]), are available in most common quantum chemistry packages. Single-centered multipole expansions of r m g ( ) ρ , which were popular in early works [5,12] are only justified at large intermolecular separation and certainly not at van der Waals distances in molecular crystals.
Perhaps the simplest method to describe molecular induction, referred to as microelectrostatic (ME), consists in distributing the molecular linear polarizability α on a finite number of points and describing the response to electric fields in terms of induced dipoles [5,11,42,43]. These type of approaches are also referred to as distributed polarizabilities or induced dipole models. The number and the location of induced dipoles mi µ on a given molecule is arbitrary. The use of few polarizable points representative of molecular subunits has been popular in the past also because of limited computational power, while nowadays distributed polarizability tensors mi α are placed at atoms, leading to more accurate results and practical implementations.
For a ME scheme combining induced dipoles and permanent charges the energy of a given molecule reads: where mi φ and F mi mi φ = ∇ are the potential and the field at atoms due to other molecules and possibly external sources. The correction to the isolated molecule energy E m g includes several contributions: the second (third) term in the right-hand side corresponds to the electrostatic energy for placing permanent charges (induced dipoles) in the potential (field) of other molecules, while the last term is the energy required to polarize the molecule.
The Applequist-Thole model [27,42] is an induced dipole scheme that has been largely adopted by the computational community. While conceptually equivalent to ME, it presents some technical differences, such as the use of isotropic atomic polarizabilities and the lack of distinction between inter-and intra-molecular dipolar interactions that are all weighted by a distant-dependent screening function to avoid numerical instabilities resulting from polarizable points too close in space [27,43]. The characteristic inverse screening distance, also called Thole factor, is an additional parameter of the model that, together with atomic polarizabilities, should be tuned to reach a good compromise between magnitude and anisotropy of the molecular polarizability [44,45].
Induced dipole schemes are most accurate when molecular polarization comes from the distortion of atomic orbitals and the response of atoms is additive. This is hardly the case in π-conjugated molecules and polymers that find application in organic electronics, whose polarizabilities are largely due to intramolecular charge flows. Also in this case some form of discretization is convenient as proposed by Stone who introduced the concept of charge flow between different molecular subunits [43,46]. Indeed modern approaches such as the charge response kernel theory (CRK) by Morita and Kato [28,47] and the charge redistribution model (CR) by Tsiper and Soos [22,36] where non-locality clearly emerges as the charge on atom i is determined by the potential acting on all the atoms j. The symmetric ij m Π tensor can be obtained from first principles Hartree-Fock [47] or density functional [28,37] calculations, or at the simpler INDO/S semiempirical level [22]. Intramolecular charge flows account for the rearrangement through bonds but miss the effect of the distortion of atomic orbitals in the direction perpendicular to the bond [48]. This is particularly evident in planar molecules such as oligoacenes, where the out-of-plane α component should be reintroduced in some way to correct this limitation. In CR this is accomplished by introducing atomic polarizabilities mi α exactly as in ME, accounting for the difference between the reference α tensor, computed at the desired level of accuracy, and its charge-flow contribution [22,36], while CRK resorts to virtual sites above and below the molecular plane [37]. ME and CR(K) describe mutually interacting nonoverlapping molecules from the atomistic structure of the solids and molecular inputs (polarizability and atomic charges or molecular potentials) that can be accurately obtained from electronic structure calculations at suitable level. Both ME and CR rely on a perturbative treatment of intermolecular interactions, a discretization of the molecular charge density, and assume linear response to electric fields. These assumptions make these two methods practical and computationally cheap while granting a good accuracy, and hence ideally suited to the description of systems of several thousands of molecules.
We now turn our attention to methods that rely on a non perturbative quantum mechanical description of the molecular response, at least for a part of the system. We indeed start from hybrid quantum mechanics/molecular mechanics (QM/MM) approaches, in which the system is partitioned into an electronically important QM region, which in this context consists of the charged molecule(s), and a larger region providing a polarizable embedding, described at a coarser classical level. This broad family of methods, originally developed for complex biological systems, received some attention also in the context of organic electronics [29,30,49], yet the promise offered by an optimal compromise between accuracy and computational cost has been not sustained so far by adequate benchmarking.
Besides established QM/MM schemes, in which a selfconsistent solution is worked out for coupled sub-systems described either at the quantum or classical level, a new type of hybrid approach appeared recently where each molecule of the system is relaxed quantum mechanically in the field of the atomic charges approximating the relaxed charge densities of the other molecules [31,32,50,51]. Accordingly, quant um patch [31,50] or self-consistent charge field [32] schemes require an iterative process in which individual molecules, one at a time, are computed quantum mechanically updating their energy and atomic charges until self consistency is achieved.
Valence Bond Hartree-Fock (VBHF) is instead a fully quantum mechanical method originally developed to describe charge transport in organic superconductors and molecules of biological interest [33,34]. VBHF builds on semiempirical Hartree-Fock method of the neglect of diatomic differential overlap (NDDO) family, and uses orbitals strictly localized on molecular fragments to build the electronic states of a supramolecular system. Charge carriers are then straightforwardly localized onto a given molecule, while allowing its neighbors to electronically relax through a self-consistent procedure, thus accounting for both electrostatic and induction effects. This approach reduces the electronic energy of a molecular cluster to the sum of isolated fragment energies and of electrostatic interactions between electron densities of fragment pairs, as expressed by equation (8). The VBHF method suffers, however, from the inherent limitations of semiempirical schemes, mostly arising from the minimal basis sets used in these models and resulting in a poor description of the molecular polarizability tensor.
A recent achievement is represented by the application of a fully DFT approach to the calculation of localized charge carriers in disordered molecular assemblies of several thousands of atoms [35]. This has been made possible by a fragment orbital DFT formalism employing optimized localized functions and featuring a linear scaling of the computing time with respect to the size of the system, in contrast to the cubic scaling of standard implementations [52]. Charge localization is here enforced by supplementing the Kohn-Sham energy func- ρ with a Lagrange multiplier constraint [53,54]: where the weighting function w r m ( ) specifies the region of space over which the charge is constrained (i.e. the molecular fragments m), N m the required charge within the specified region, and V m is the Lagrange multiplier to be determined upon optimization. To describe diabatic CT states, a concurrent constraint on the molecular spins is usually added, to force the electron donating and accepting molecules to carry an excess spin of 1 2 ± . Constrained DFT (CDFT) offers a technically different pathway to the computation of the energetics of localized charge carriers with respect to zero overlap methods, with the advantage over VBHF of a superior description of the molecular electronic structure and polarizability. Moreover CDFT, by explicitly accounting for intermolecular exchange interactions, missed in the zero overlap approximation, allows not only to discriminate between different spin configurations, but can also be applied to compute exchange couplings in magnetic systems [55]. However, the computational cost of the CDFT approach of [35] is considerably higher than, for instance, induced dipole or QM/MM schemes, still posing some limitations to the treatment of long-range electrostatic phenomena.

Electrostatic interactions in extended systems
It is well established that the evaluation of electrostatic interaction energies requires some care. These interactions are long-ranged, that is one cannot truncate the interaction sum by simply neglecting pairs of molecules with separations exceeding a certain cutoff distance. Moreover, these sums in some cases converge only conditionally, i.e. the result of the summation depends on the dimensionality and on the macroscopic shape of the system. This problem has already been acknowledged by Madelung when evaluating electrostatic sums of ionic crystals, where convergence in 3D can only be achieved under specific conditions. The situation is even more complex in the case of charged excitations, where instead of dealing with electrostatic sums over charged atoms or molecules in a neutral unit (super)cell, one has to compute the interaction between an excess charge and surrounding neutral molecules, see equation 12. In ordered organic semiconductors, the first non-zero molecular multipole of neutral molecules is typically the quadrupole (see figure 2). Of course dipolar molecules also exist but in a crystal phase they often arrange in an apolar cell that can be still treated as a (super) quadrupole. The leading term in equation (12) is then the interaction between the excess charge and quadrupole moments of the surrounding neutral molecules. The charge-quadrupole interaction energy decays as 1/r 3 , leading to a conditional convergence in 3D and extremely slow convergence in 2D films.
Let us illustrate the conditional convergence of electrostatic sums with the help of a simple lattice of permanent multipole moments Q lm of order l (l = 0, 1, for monopole, dipole; m l l , , = − … ). The electrostatic energy of a charge q placed for convenience at the origin lattice site can be written as where the information on the lattice structure and on the multipoles orientation have been adsorbed in the lattice sum where R a b c i j k ijk = + + points to lattice sites defined by vectors a, b and c, and Y l m are spherical harmonic functions. The primed sum excludes summation over the site bearing the charge.
In the case where only molecular quadrupoles (l = 2) are concerned, the sum of terms R lm ijk ijk 3 ξ ∼ − converges in infinite 3D systems thanks to its angular dependence (spherical harmonics), in the absence of which a logarithmic divergence occurs. In turn, the value of lm ξ depends on the order of summation and hence on the macroscopic shape of the sample, making the electrostatic energy of a charge E ∆ at the center of a sphere of infinite radius (bulk limit) different from the value for an infinite slab ( film limit), irrespectively on its thickness. It is worth noticing that, conversely, the sum m shape independent for a lattice of quadrupoles and acquires a weak shape dependence in realistic molecular systems.
In infinite 2D systems lm ξ is instead absolutely convergent, although such a convergence may be very slow for thick films. The lack of translational symmetry along the film normal z leads to position-dependent energies, although energy profiles are rather flat with sizeable differences only for sites at the film surfaces. Moreover, due to the angular dependence of lm ijk ξ , the value of E ∆ depends on the crystal facet along which the film is cut, or the orientation of a molecular quadrupolar tensor with respect to the film normal. These general considerations apply also to off-lattice calculations of realistic morphologies with extended molecular charge distributions, for which compact expressions like equations (17) and (18) are not available.
We will now briefly review numerical methods for the calcul ation of long-range electrostatic interactions in molecular systems, whose description requires specific summation techniques, especially when periodic boundary conditions are used to mimic macroscopic systems. Notice that this is not only the case of units cells of perfect crystals, but also (super) cells containing a large number of molecules as typically obtained from finite temperature simulations of disordered systems. The direct evaluation of Coulomb sums represents a simple and generally viable option that requires calculations on systems of increasing size until convergence within a specified tolerance, or extrapolations to the infinite system limit.
In neutral periodic systems one can take advantage of translational symmetry and resort to reciprocal space-techniques, such as the Ewald [56] or the Parry [57] algorithms for periodic boundary conditions in 3D and 2D, respectively.
If one or more charged species are present, it is no longer possible to use standard summation techniques, since charges are not repeated in periodic replicas. Again, one can resort to a direct evaluation and extrapolation of sums, paying however attention to the shape of the system that is crucial for a correct extrapolation, as illustrated in figure 3 for the pentacene crystal and films. The bulk limit is recovered for spheres centered around the charged molecule, while reaching the film limit requires cylinders with axis normal to the film and radius much larger than film thickness. In both cases chargequadrupole interactions scale as R −1 , where R is the radius of the sphere or cylinder, making direct sums prohibitive in the case of thick films. Alternatively, reciprocal space techniques can be extended to charged systems, as done by Tsiper and Soos [22] and Poelking and Andrienko [24,58,59] for periodic boundary conditions in three and two dimensions, respectively. These approaches, by partitioning the charged system into a periodic background plus one or more aperiodic excitations, combine the accurate evaluation of long range electrostatic interactions with the computational efficiency of reciprocal space sums.
Induction interactions can be also treated exactly in the case of periodic neutral systems, while when aperiodic excitations break translational symmetry the dipoles induced by the charge carrier(s) have to be computed within a given cutoff, yet accounting for their interactions with the polarized background of the infinite neutral system [22,24]. The extrapolation approach to induction term I ∆ is instead shown in figure 3, where appropriate power laws are used to obtain the values of infinite bulk and films.
We anticipate that the aforementioned effects are important for interpreting photoelectron spectra of ordered organic films. Indeed, here both the ionization energy and electron affinity have a sizable electrostatic contribution which originates from long-range interactions between the excess charge and quadrupolar moments of neutral molecules. This is further discussed in sections 4.2 and 4.3.

Bulk organic semiconductors: crystal packing and disorder
We discuss here the application of methods introduced in section 3.1 to the calculation of charge transport levels in bulk organic semiconductors as obtained with electrostatic sums in 3D according to the discussion in section 3.2. We will not attempt a comparison with experimental values of ionization energies since these are surface quantities as extensively discussed in the next section. The following discussion focuses instead on methodological aspects, the relationship between molecular properties, crystal packing and polarization energies in bulk systems, and disorder. ∆ shows an oscillating convergence in the bulk and R −1 behavior in films. I ∆ follows R −1 and R −2 trends in bulk and films, respectively.

Localized charges in perfect crystals.
Pioneering work on the calculation of charge carriers polarization energies in molecular crystals dates back to the late sixties/early seventies, with seminal contributions by Lyons and co-workers [60,61], Munn and co-workers [11,12,25] and Silinsh and co-workers [13,26]. The merits of these very first attempts have been several, from the identification of charge-induced dipole and dipole-dipole interactions as the leading contributions to the polarization energy, up to the exploration of different partitioning schemes of the excess charge over the molecular ion [60,61]. Self-consistent polarization field approaches following the classical work by Mott and Littleton for atomic lattices [62] were next developed and applied to acene crystals [13], including the Fourier transform method by Bounds and Munn [11], in which the shape and orientation of π-conjugated molecules started to be taken into account through distributed polarizabilities on molecular regions (submolecules). These calculations were however blind with respect to the nature of the charge (positive or negative) [12]. An important advancement in this respect was the observation that the so-called charge-quadrupole interaction breaks the symmetry of the electronic polarization for opposite charge carriers assumed in previous studies [12,26]. Before proceeding, it is important to provide an order of magnitude assessment of the two dominant, electrostatic and induction, contributions. The induction term ( I ∆ in equation (11)), referred to as charge-induced dipole in early works [5,11], can be roughly estimated by the Born equation for the polarization energy of a charge in a spherical cavity of radius R inside a dielectric medium: where r is the dielectric constant of the medium. A generalization of equation (19) to anisotropic dielectrics is due to Bounds and Munn [11]. Typical values for molecular crystals (R = 5 Å, 3.5 r ε = ) yield 1 I ∆ ≈ − eV for both holes and electrons, a value that is consistent with submolecules or atomistic calculations that instead explicitly account for the details of the structure and molecular anisotropy. The Born equation provides a general understanding of the decrease in magnitude of I ∆ with the molecular size, a trend that has been confirmed by the UPS measurements on 44 organic semiconducting crystals by Sato et al [26], and by state of the art ME and CR calculations (see figure 4(a)) [23]. The polarization of the medium is the main responsible for the large (about 2 eV) reduction of the transport gap E t = IP − EA in organic materials with localized carriers with respect to its gas-phase counterpart IP EA g g − . The electrostatic contribution, E ∆ in equation (11), accounts for the interaction of the charge carrier with the charge distributions of neighboring molecules. Although this term has been originally introduced to describe quadrupolar molecules such as oligoacenes [12], the definition in equation (12) applies to arbitrary charge densities r g ( ) ρ . The electrostatic term has opposite sign and approximately the same magnitude for positive and negative charges, hence representing the main source of difference between hole and electron polarization energies (see below). The magnitude of E ∆ essentially depends on and on the supramolecular packing, much more than I ∆ does. For instance, along the well-studied series of linear oligoacenes, the absolute value of E ∆ increases from 0.2 to 0.35 eV, correlating with the strength of the molecular quadrupole (figure 4(a), middle panels). These values, obtained with ESP atomic charges [23], are consistent with earlier estimates based on different partitioning of the excess charge and point quadrupoles over molecular units [5,12] and agree well with values obtained from DFT electron densities [63]. We stress that E |∆ | can be in some cases comparable or even larger than the induction term, as in the cases of perylenetetracarboxylic acid dianhydride (PTCDA) [22,63] or a perylene bisimide derivative with fluorinated side chains [23]. Ryno et al focused instead on the role of molecular packing by considering oligoacenes and their TIPS [6, 13-bis(2-(tri-isopropylsilyl)ethynyl)] derivatives [64]. This is an interesting case study as the two series of molecules, acenes and TIPS-acenes, feature comparable quadrupole and polarizability but they form very different crystals characterized by herringbone and brickwork packing, respectively. Different structures lead to different intermolecular interactions, resulting in strikingly different polarization energies (see figure 4(b)), and smaller electron-hole asymmetry in TIPS compounds. For instance, P P P 0.23 ∆ = − = + − eV in pentacene against 0.10 eV for the TIPS analogue [64].
Various modern computational approaches have been applied to study electrostatic and induction interactions in molecular crystals [22,30,33,64,65]. We conclude this part by analysing in depth how the P + versus P − asymmetry builds up as different approximations are progressively lifted in the modelling, specifically referring to the work presented in [23]. This is illustrated by figure 4(c), showing hole and electron polarization energies of anthracene computed with CR and molecular models of increasing sophistication numbered from I to V. At level I, we neglect the effect of permanent charges of neutral molecules and assume a uniform charge distribution for cation and anion, which bear the same polarizability of the neutral molecule. In this case, P P 1.2 eV, in agreement with the original ME work [11]. At level II we add the electrostatic interaction between the charge carrier and the permanent charge distribution of the surrounding neutral molecules, which results in a large asymmetry, eV, as predicted in [12]. At level III, a fully self-consistent calculation of induced multipoles in the field of permanent charges is performed. Because induced dipoles screen the field generated by the charge carriers, their interaction with the surrounding neutral molecules is reduced and so is the polarization energy asymmetry. Further refinements consist in using the atomic ESP charges (IV) and polarizability (V) specific to molecular cation and anion, which contribute in this case to slightly lower P P − − + .

Static versus dynamic disorder.
We next focus on the effect of disorder that strongly affects the energetics of localized charge carriers in soft organic semiconductors, ultimately hampering their charge transport properties. In fact, the different environments experienced by charge carriers at each molecular site result in a broad distribution of transport levels, a phenomenon that is akin to the inhomogeneous broadening brought by solute-solvent configurations in spectroscopy. The energy distributions of charge carriers often feature an approximately Gaussian shape, whose standard deviation σ is usually taken as a measure of the energetic disorder.
Most of the approaches presented in section 3 can be also applied to systems with no translational symmetry, yet with the additional hurdle that usually a large configurational space has to be sampled. This necessity makes the requirement of computational efficiency even more stringent, especially for amorphous systems. In the following we distinguish the two cases of static (or positional) disorder, related to energy differences among molecular sites, and dynamic disorder, when individual site energies vary in time as a result of thermal molecular motion.
Static disorder is typical of amorphous systems, but is also present near structural defects. The simplest point defect one can imagine is the vacancy of a molecule in a perfect crystal lattice. This case has been analyzed by Eisenstein and Munn in oligoacene crystals, reporting variations of site energies in the range of 0.1-0.3 eV of both positive (scattering centers) and negative (traps) sign [66]. Grain boundaries are unavoidable defects in molecular films. Verlaak and Heremans applied ME calculations to study how the energy levels of single holes are affected in proximity of four idealized grain boundaries [67]. It was found that electronic states near grain boundaries give rise in most cases to energy barriers though some structural arrangements also result in trapping centers. Both effects are sourced primarily by charge-quadrupole interactions and are expected to hinder severely charge transport.
The recent studies of energetic disorder in molecular solids mostly rely on the combination of the methodologies discussed so far with classical molecular dynamics (MD). The latter technique allows obtaining reliable atomistic morphologies of disordered large molecular systems, as well as sampling the effect of thermal molecular motion [68]. Tris-(8-hydroxyquinoline)aluminum (Alq 3 ) is a common n-type organic semiconductor studied in several theoretical works [29,31,50,69] that well illustrates the role of intermolecular interactions on energetic disorder. In figure 5(a) we show the distribution of site energy differences between neighboring Alq 3 molecules computed using different flavors of the quant um patch method proposed by Wenzel and co-workers [31,50]. These rather broad distributions ( 0.2 σ ≈ eV) due to the dipolar character of Alq 3 and to the positional disorder in the amorphous phase, compare well with those reported by Ruhle et al, shown in the inset of figure 5(b) and based on the simpler Thole model [69].
The energy landscape in fullerene derivatives of different molecular structure and morphology has been recently studied with ME calculations by D'Avino et al [71]. This work showed that while amorphous phases of the polar PC 61 BM and PC 71 BM derivatives sustain an energetic disorder comparable to Alq 3 ( 0.16 σ = eV for both), the disorder is strongly reduced in crystalline PC 61 BM ( 0.09 σ = eV) and almost negligible ( 5 σ = meV) in the highly symmetric and apolar C 60 fullerene, entirely due to the impact of the thermal motion on the induction term (see inset of figure 5(c)).
These and other studies allow drawing some general conclusion on the energy landscape in bulk organic semiconductors: (i) The energetic disorder is largely electrostatic in nature and originates from the potential probed by each molecule in its specific environment, as generated by the charge densities of the surrounding disordered neutral molecules. The dis order hence depends on the type and magnitude of molecular multipole moments as well as on the (positional, orientational, conformational) order in a given material. (ii) The induction of microscopic dipoles, mainly reacting to the localized carrier, stabilizes the transport levels but, most importantly, it strongly reduces the energetic disorder by 20-50% [65,69,71,72]. This smoothening of the energy lansdcape results from the dielectric screening provided by molecular polarizabilities and manifests in a pronounced anti-correlation between E ∆ and I ∆ in equation (11) [71]. (iii) The site energies on different molecules are not independent but spatially correlated and correlations appreciably extend up to 2-3 nm [69,71,73]. The decay with distance approximately follows the r −1 behavior expected for random dipolar disorder [74].
The distributions in figure 5 actually hinder the information whether the energetic disorder is static or dynamic in time, as the total variance is just the sum of the two 2 dyn 2 sta 2 σ σ σ = + . For fullerene derivatives it has been shown that actually the dipolar energetic disorder is static with respect to charge carriers dynamics, varying with a characteristic time scale of a few hundreds of ns, which corresponds to the slow rotational motion of this nearly spherical molecules [51,71]. Figure 6(a) takes the example of an amorphous PC 71 BM systems to illustrate how a broad polarization energy distribution builds up from large positional disorder and relatively small dynamic disorder. Disentangling static and dynamic disorder is difficult in reality but rather simple in the simulations, although the results are dependent on the length of the observation window. In fact, the dynamic disorder is just the average standard deviation of the site energy of single molecules with respect to its individual mean, while the static one is the standard deviation of the individual means with respect to the global mean.
Dynamic disorder associated with thermally-activated molecular motion represents the only source of dis order in ideal defect-free molecular crystals. Martinelli et al used both MD simulations and normal mode sampling coupled to VBHF  to evaluate the impact of the lattice and molecular vibrations on site energies in anthracene single crystal [75]. This study allowed the estimation of the purely dynamic contribution to the spread of transport levels ( 0.1 σ ≈ eV at 300 K) and, most interestingly, to assess the characteristic timescales of the modes modulating IP and EA. The Fourier transforms of the IP sampled at 5 fs intervals along the MD trajectory in figures 6(d)−(e) show that while there is a contribution from low-frequency intermolecular modes, energy fluctuations are mostly due to intramolecular modes associated with the carbon atoms backbone for which, however, classical treatment of nuclear dynamics should be inappropriate [75].

Transport levels in molecular films and photoelectron spectroscopy
Photoelectron spectroscopies are the reference techniques for the determination of transport levels in organic materials [2,18,76]. Ultraviolet photoelectron spectroscopy (UPS) measures the first ionization potential (IP) and other valence occupied states, and similarly does x-ray photoelectron spectroscopy (XPS) with core levels. Inverse photoemission spectroscopy (IPES) probes instead unoccupied levels, most often the electron affinity (EA). The energy resolution of UPS and XPS are on the order of 0.1 eV, and a comparable resolution has been attained only recently in IPES through the use of near-ultraviolet photons (low-energy inverse photoemission spectroscopy, LEIPS) [77,78]. Photoelectron spectroscopies are surface sensitive techniques probing the topmost 1-3 nm of a sample, corresponding to a limited number of molecular layers. Hence UPS or IPES measure the transport levels at surfaces, either of thin films or of crystals taken as semi-infinite solids, and not the bulk transport levels discussed in the previous section.
Charge carriers energetics in molecular films can be also probed with other techniques, such as scanning tunneling spectrosopy (STS) [79,80] or two photon photoemission spectroscopy (2PPE) [81,82]. A simple and qualitative insight on transport levels can also be obtained from cyclic voltammetry (CV) on solid state films [83,84]. Nevertheless, each of these techniques has its own specificities that generally lead to different values for the same measured quantities. For instance, STS determines the density of electronic states by measuring the current flowing between the semiconductor surface and the tip of a scanning tunneling microscope, whose interaction with the molecules cannot be neglected. On the other hand, 2PPE is a pump-probe technique tracking nonequilibrium electronic dynamics, while the phenomena at the interface to the electrolite solution do impact the redox potentials measured in CV. The theoretical techniques discussed in this review can in principle account for the specialties of each experiment, although our primary interest is in the determination of IPs and EAs of organic films, as directly probed by UPS and IPES.
The following discussion mostly focuses on UPS/IPES spectra of molecular films on metallic or dielectric substrates, i.e. the transport levels relevant to organic electronic applications. Nevertheless, we start the discussion with XPS spectra of noble gases, an almost ideal example of van der Waals crystal that most clearly and simply illustrates surface effects without any additional effect complicating the interpretation and modelling of molecular films (see below). Figure 7(a) shows XPS spectra of Xenon clusters of different size as obtained by supersonic adiabatic expansion [85]. These spectra present well resolved peaks corresponding to core holes for isolated atoms, which serve as a reference, and to atoms at the surface and in the bulk of the clusters. The relative intensity of the surface to bulk peak decreases with increasing cluster size, according to the ratio between bulk and surface atoms. Electronic polarization reduces the binding energy by 1.2 eV when a localized hole is created in the inner cluster region and by 0.9 eV for holes at the surface. These effects can be quantitatively modelled with ME calculations relying on a cluster structure and atomic polarizabilities α, as discussed in detail in [86]. Our results in figure 7(b) for a spherical cluster of N = 887 Xe atoms (face centered cubic with lattice constant a = 6.2 Å, ) further show that the polarization energy for an atom at the center of the cluster is still 0.1 eV higher than that of an infinite 3D bulk (dashed line). Surface effects can be also detected in a ten-layer film with a difference of 0.2 eV between outer and inner layers, as shown in figure 7(c).
The modelling of transport levels of organic molecular films poses additional challenges, the most important being: (i) molecules, in contrast to atoms, possess an internal structure that tends to relax upon ionization. Nuclear intramolecular reorganization energies are on the order of 0.1 eV per charge and usually decrease with molecular size [15,87]. Nuclear reorganization is a fast process and its contribution is usually considered as a correction to the vertical electronic energy [5].
(ii) The anisotropy in shape and charge density of the large conjugated molecules used in organic electronics produces charge-multipole interactions that can be as large as induction effects and that in addition strongly depend on molecular orientation. Among those multipoles, the quadrupole moment is the leading contribution for molecules, or crystal cells, with an inversion center. (iii) The substrate plays a key role in thin films of a few molecular layers (ML) and requires an explicit modelling. In the following we address the specific features of molecular films, starting from the role of the substrate.
Molecular films on insulating substrates are usually modelled as free-standing films, leading to results in good agreement with experiments [24,88,89]. To illustrate the accuracy with which ionization energies of thin-film can be predicted in silico, figure 8 compares calculated (Thole model) and measured values for a set of molecules covering a wide spectrum of gas-phase IP, ranging between 5.2 and 6.8 eV [24]. To address the surface sensitivity of UPS, simulated hole energies are shown as a function of the penetration depth. The simulation results (blue bars) not only quantitatively reproduce measured trends for the orientation dependence, but also the absolute experimental IP values extracted from UPS (black bars).
Actually the electrostatic potential of a dielectric substrate, its polarizability and its roughness should all contribute to the transport levels of the organic semiconductor. For instance, Martinelli et al studied the impact of different polymer di electric substrates on the energy landscape of a pentacene film [90]. They found a larger energetic disorder in the first layers on poly(methyl methacrylate) than on polystyrene and ascribed this effect to the rougher electrostatic landscape of the former, featuring polar carbonyl groups. Similar effects can be produced by chemically different self-assembled mono layers (SAMs) coating a given dielectric, or even by coverage inhomogeneities of the SAM [91].
Metallic substrates require additional care with respect to dielectrics, due to possible shift of the vacuum level at the interface and because of the infinite polarizability of the metal. The two effects can be treated separately to a degree of accuracy that allows understanding most of the exper imental evidences. The typical energy-level diagram of a generic metal-organic interface is exemplified in figure 9(a): notably, a shift of the vacuum energy level, ∆, is a common feature of metal-organic interfaces that typically saturates at the coverage of 1 monolayer (ML) of molecules. Ishii et al identified (a) XPS spectra of 4d 9 core holes in (Xe) N clusters. The two peaks replica correspond to ionization processes from spin orbit coupling-resolved levels ( j = 5/2 and 3/2). Reprinted from [85] with the permission of AIP Publishing. Right-column panels show the position-dependent polarization energies (P) in a Xe 887 cluster (b) and in an infinite film of ten atomic layers (c) as obtained from ME calculations. Dashed lines correspond to the P extrapolated for an infinite bulk. In noble gas systems experiments and theory quantitatively agree on the presence of well-resolved peaks for ionization from the surface and from inner regions. several possible contributions to ∆ [1]. Among those, metalorganic charge transfer, and the consequent formation of an interface dipole layer, is considered the most important term. In this context, the terms metallization or metal-organic hybridization are sometimes used as synonymous of charge transfer. The theoretical description of this phenomenon requires accounting for the overlap between the two materials and goes beyond the scope of this review. We only mention that cooperative electrostatic interactions between metalorganic charge-transfer dipoles account for the nonlinear variation of ∆ as a function of sub-ML coverage [92].
For molecular layers not in direct contact with the metal the interaction with the substrate can be instead safely described in terms of image charges and image induced dipoles. This approach has been extensively applied by Soos and coworkers in computing transport levels at metal-organic interfaces with the charge redistribution model [80,88,93,94]. Figure 9(a) sketches their typical calculation setup where a few ML-thick organic film is placed on top of an inert metal surface. A hole (or an electron) is placed at a molecular site of a given layer and the polarization energy P + (P − ) can be computed selfconsistently as a function of the ion position. The approach was successful in reproducing the 0.3-0.4 eV increase of the transport gap E IP EA t = − measured for perylenetetracarboxylic dianhydride (PTCDA) on silver and gold substrates (see figure 9(b)) when going from one ML (0.4 nm) to a thick film (>5 nm) [80].
On general grounds, the asymmetry induced by the metal and vacuum interfaces, along with long-range electrostatic and induction interactions, leads to transport levels that vary from ML to ML, as shown for 10 ML thick films of different molecules in figure 9(c) [93]. The C 60 (1 1 1) film features P P = + − since high molecular symmetry leads to negligible charge-multipole interactions ( 0 ). P + is calculated to change by ∼0.1 eV between molecules at metal and vacuum interface. This result is consistent with the XPS measurements by Maxwell et al [95] that reported a 0.15 eV vacuum level shift upon deposition of the first ML, and then small changes in the C(1s) binding energy with increasing thickness. In contrast to what is observed in Xe N clusters, it was not possible to resolve specific XPS features from surface or subsurface C 60 layers. In the case of films of pronouncedly anisotropic molecules as anthracene and PTCDA, charge-quadrupole interactions adds leading to an electronhole asymmetry P P ≠ + − qualitatively similar to that discussed in section 4.1 for bulk systems. The predicted layer dependence of transport levels remains weak, exception made for P − in PTCDA that increases substantially from the surface to the metal [93]. The calculated change of ∼0.2 eV between the surface and subsurface layer may now be within the resolution of IPES spectra.
We conclude this section by discussing the influence of molecular anisotropy and orientation on photoelectron spectra of molecular films. As recognized by Koch and coworkers [96], pentacene (PEN) and perfluoropentacene (PFP) films represent an ideal case study for this purpose as they allow to sharply distinguish between electrostatic and polarization contributions to transport layers. The two molecules have in fact the same π-conjugated backbone (see figure 10(a)), which results in nearly equal molecular polarizabilities, while As shown in figures 10(c) and (d), the IP of films of standing PEN molecules is 0.55 eV lower than for lying PEN on metal, while PFP films feature an IP increase of 0.85 eV when going from standing to lying films. The experimental picture (see figure 10(d)) has been very recently completed by the LEIPS measurements by Yoshida et al [89], which also reported opposite trends of comparable magnitude for EAs of PEN and PFP.
The opposite trends in transport levels, summarized in figure 10(d), have an electrostatic origin that can be qualitatively interpreted as the difference of electrostatic potential between the position of the ionized molecule (e.g. at its centroid) and the vacuum level above the film [16,17,88]. Since molecules in films with standing or lying orientations expose different components of their quadrupole moment to the vacuum, they experience different electrostatic potentials in the two cases, resulting in orientation dependent IP and EA [16,89,96]. This effect is similar to the dependence of the work function of a metal on the crystal face, although in the latter case it is due to the electrostatic potential step at polar surface [1,97]. Opposite quadrupole components for PEN and PFP rationalize the different trends for the two molecules. On the computational side, Topham et al calculated the IP of PEN and PFP films with the charge redistribution model, considering bilayers of standing molecules and monolayers of lying molecules on a metallic substrate [88]. Theoretical calculations based on a very similar methodology have been successfully extended to EA by Yoshida et al [89]. The agreement of calculations with UPS and LEIPS data is within the experimental resolution, as shown in figure 10(d).
We conclude this section by remarking that the overall agreement between experiments and theory, as for instance evinced from figures (8)-(10), or from the very recent report of continuous band energy tuning in mixed films of hydrogenated and halogenated phthalocyanine derivatives [98], demonstrates the maturity of the current understanding of ionization processes in molecular films and the robustness of the modern theoretical approaches. The discrepancy between experiment and theory, comparable to experimental uncertainty, can be attributed to different phenomena, including intramolecular relaxation, charge delocalization and disorder. Theoretical calculations have therefore a predictive value that can establish a link between the electronic properties of the molecules and those of the materials.

Organic heterointerfaces for photovoltaics
In this last section of our critical survey, we discuss the latest literature results on the energy landscape at organic-organic interfaces. These studies are fostered by the importance of these heterojunctions in organic photovoltaic devices, where the active layer is constituted by a phase separated blend of two different organic semiconductors, a p-type electron donor (D), and a n-type electron acceptor (A). The active layer can be somehow regarded as a p-n junction where a large level offset of chemical origin, i.e. the photovoltaic gap (see equation (7) and figure 1) is applied to prompt the separation of tightly bound and localized photo-generated electron-hole pairs.
As a consequence of the large variety of materials, ranging from small molecules (see figure 2) to polymers, and of the strong dependence of structural organization on the device preparation technique, the morphology of organic heterojunctions is a rather ill-defined and system-dependent concept, difficult to generalize and even to probe experimentally. An optimal interface morphology should fulfil the following requirements: (i) be thick enough to absorb most of the incident light (≈100 nm); (ii) contain D and A domains of typical size comparable to the exciton diffusion length (≈10 nm), so that the exciton can reach the interface before decaying; (iii) have a large D/A interfacial area to maximize the probability exciton dissociation; (iv) D and A domains should all be in contact with anode and cathode, respectively, in order to collect the totality of generated charges.
Since all these requisites are difficult to achieve simultaneously, real morphologies can be quite different from the ideal picture, and range between two extreme situations, a completely planar bilayer heterojunction and a fully interpenetrated bulk heterojunction, as schematically depicted in figures 11(a) and (c). As we shall see below, also the simplicity of planar heterojunctions geometry is largely apparent, because different relative orientations of D and A molecules are possible at the interface (figures 11(d) and (e)), where the miscibility of the two components can also lead to some intermixing (figures 11(b) and (f)). Of course intermixing and interfacial disorder can occur both in the planar and in the bulk heterojunction case.
The theoretical methods for the calculation of the energetics of charged excitations presented in section 3.1 proved to be very useful in bridging the gap between molecular properties and device performances, giving insights on the working mechanism of organic solar cells. Charge carrier energy levels at organic heterojunctions provide information on the electrical work that can be extracted from an organic solar cell, i.e. its open circuit voltage V oc , whose low values represent a major bottleneck for efficient applications. In addition, the analysis of electrostatic forces at D/A interfaces rationalizes the high charge generation quantum yields and photocurrents of best performing devices, despite the strongly bound excitons characterizing these low dielectric media.

Energetics of multilayers and open circuit voltage.
In this section we will see how long-range electrostatic interactions in a polarizable environment can be applied to the study of charge carriers energy profiles in layered organic solar cells, resorting to the approach for infinite 2D systems by Poelking and Andrienko mentioned in section 3.2 [24,58]. Before turning to real materials, we propose the illustrative example of a bilayer composed of a bcc lattice of polarizable points. Lattice positions in the z < 0 half-space bear no permanent multipole, sites for z > 0 are characterized by a quadrupole moment. In the context of organic solar cells, the former corresponds to a fullerene (C 60 ) A unit, the latter to a D molecule with a finite quadrupole moment. Next to these molecular sites, we define ficticious sites on either side of the bilayer, which can be singly charged in order to probe the vacuum level above the thin film.
We parametrize the lattice model by choosing sites polarizabilities such that they result in a dielectric constant 3 r ε ≈ for a lattice constant c = 0.77 nm. As opposed to the non-polar A sites, D ones carry an axial quadrupole with largest principal component Q 2 whose axis is tilted by an angle θ with respect to the interface plane (see inset of figure 12(a)). Q 2 is chosen negative (−10 a.u.) as an abstraction of the electrostatic layout of typical D units. A negative quadrupole component is for instance found along the molecular axis of dicyano-substituted thiophenes, due to the strongly electron accepting terminal groups (see D5M in figure 2), but also normal to the molecular plane of common π-conjugated molecules as unsubstituted thiophenes or acenes. It is worth noticing that tuning the angle θ in such a model can be both interpreted either as a change of the chemical structure of the D molecule, or as a variation of the orientation of a given molecule. Figure 12 configuration, respectively. To interpret the profiles, we first note that the dielectric stabilization at this heterojunction would amount to almost 1 eV in the absence of permanent quadrupoles. Any deviation from this average value is hence attributable to the action of the quadrupolar sites in z > 0. Notably, the polarity of the interface changes from the 0 θ = to the 90 θ = case. In the former scenario, holes are more stabilized in the donor part than the acceptor part, and vice versa for electrons. In the latter, 90 θ = , scenario, the situation is reversed, with a twofold increase in D-A offset resulting from the alignment of the largest quadrupole component Q 2 parallel to the interface normal.
Beyond the dielectric stabilization of charge carriers, the polarization response of the lattice sites to the abrupt step in quadrupolar fields generates interfacial dipole layers [24,59]. As can be seen in figure 12(b), the dipole profile is sharply peaked across the interface, with a slight kickback in polarization density over the adjacent layers as a consequence of the bcc packing. These dipole layers are responsible for the vacuum-level shifts observed on both sides of the thin film ( figure 12(a)).
What happens if we now sandwich a 0 layer in between the non-polar layer and a 90   in the top-layer (right-hand panel). The lattice sites are occupied by apolar (z < 0) and quadrupolar (z > 0, Q 0 2 < ) polarizable particles with orientations θ as specified. (b) Interfacial dipoles generated by quadrupolar fields for the bilayer (left-hand panel) and trilayer (right-hand panel) configurations. Additivity of the peaks indicates that vacuum-level alignment holds for the trilayer system, as can also be seen from direct comparison of the vacuum levels of holes and electrons and the relative positioning of the interlayer, see the arrows in the right-hand panel of (a). In the right column panels, the data for bilayers are repeated in grey as a guide for the eye.

Dissociation of charge transfer states.
Finally, we would like to discuss the current understanding of the mechanism that leads to an efficient dissociation of charge transfer (CT) states into free carriers at organic D/A interfaces. The assumption of localized charge carriers offers a relatively simple framework for grasping the different factors playing a role in the e-h separation in a cold scenario, i.e. when charges split without the need of an excess excitation energy. Recent experiments suggest that while hot processes involving excited states do occur leading to ultra-fast charge separation [101,102], high quantum yields do not require excess energy [103,104]. Following [105] the energy of an e-h pair reads: where we introduced the e-h polarization energy P ± measuring the contribution of intermolecular interactions to E ± , in analogy to its single charge analogues in equations (5) and (6) This equation highlights different contributions to P ± : single-charge electrostatic energies, E / ∆ + − , the unscreened e-h Coulomb attraction V ± that should be overcome to obtain free carriers, and finally the induction contribution I ∆ ± , which requires a self consistent calculation of the induced dipoles at the interface for each specific position of the two charges. Moreover, structural disorder leads also to IP g and EA g that depend on the specific geometry of each molecular site [105,106].
We next quantify the different contributions to P ± , starting from the Coulomb e-h attraction. Here, a simple unit point charge approximation provides a crude upper limit of V 2.9 ≈ − ± eV for an intermolecular distance of 5 Å. Better estimates can be obtained with extended charge distributions, leading to results between −2.0 and −1.1 eV for nearest neighbours with different packings and orientations [105][106][107][108], as shown in figure 14(b).
The polarization of the environment by the charge pairs gives rise to an important reduction of the energy barrier for charge separation that we again illustrate by means of a lattice model. In this case we consider a simple cubic lattice of polarizable points with no permanent multipoles except an electron and a hole that we take apart while computing the relevant interaction energies as a function of their distance r eh . This elementary model, entirely defined by the lattice spacing and dielectric constant, is not specific to interfaces and can be considered as a minimal picture of e-h splitting also in pristine materials [22] or in the presence of dopant molecules [88,109]. Figure 15 (see equations (4) and (7)). The limits of this macroscopic arguments are made explicit in figure 15(b), where we compare the microscopic and numerically exact E b ± and its macroscopic approx imation r . The latter appears justified when charges are far apart, while significant differences occur when e-h distances are comparable to molecular length even in the ideal case of an isotropic cubic lattice. This highlights the importance of a microscopic description of dielectric screening, especially in realistic offlattice morphologies of molecular systems.
Atomistic calculations of the induction energy of interfacial and separated e-h pairs for different model interfaces [107,108], shown in figure 14(c), confirm the qualitative trend in figure 15. These results obtained for spherical clusters of finite radii provide, however, only a lower limit to the effect of the dielectric medium. The large difference between the values of I ∆ ± obtained by Verlaak et al [107] and Gorczak et al [108] for the same pentacene (0 0 1)/C 60 (0 0 1) interface might be attributed to different polarizability inputs or to the size of molecular clusters employed in the calculation.
Even considering the contribution of the dielectric medium, a residual binding energy of 0.5-1 eV is still too large for an efficient dissociation of CT states at D/A interfaces. We hence focus on interface electrostatics, namely the dependence of E ∆ + and E ∆ − on the molecular position, which may also affect the energetics of charge separation at D/A interfaces. Actually, electrostatic energy profiles at interfaces suffer from the same ambiguity of transport levels of pure materials and depend on the dimensionality and macroscopic shape of the systems. For an ideal bilayer with neat D/A interface as the one sketched in figure 11(a) the superposition of quadrupolar fields results in a homogeneous electrostatic potential in both materials, thereby not providing any additional force for charge separation (see figure 12). Interfacial roughness is, however, unavoidable in real-life bilayers and a mild D-A intermixing consisting in nanometric protrusions of one phase into the other as sketched in figures 11(b) and (f ) can radically alter the local electrostatic landscape.
Poelking and Andrienko recently showed that charges in these protrusions actually experience electrostatic forces that amounts to tenths of an eV over the distance of about one nanometer. Depending on the orientation of molecular quadrupoles and on the curvature of these protrusions, these forces can both keep charges stuck at the interface or push them apart [58]. Where the second condition is met, CT states in the protruding region will have energies close to the photovoltaic gap, allowing the dissociation of cold excitons into free charges through a barrierless process. The same forces can also keep charges away from the interface, preventing non-geminate charge recombination [58]. For this to happen, specific orientation of molecular quadrupole moments at the interface, as well as reasonable intermixing of two phases are needed. An interesting observation is that large charge push-out forces prompting  an efficient CT-state dissociation do also result in a reduction of the open circuit voltage, since they reduce the photovoltaic gap Γ [58]. Thus, a good balance between energy level alignment and electrostatic forces assisting CT state dissociation and charge detrapping has to be achieved, making material and processing conditions such a formidable challenge.
An alternative approach to interfacial electrostatics consists in the use of finite size clusters, approximating the junction between two semi-infinite solids [105,107,108]. In this case single charge electrostatic energy profiles z E ( ) / ∆ + − are not flat as for the infinite bilayer in figure 12, but rather feature band bending effects [110], resulting in a net difference between the energy of a charge at the interface and in the bulk (see figures 14(d) and (e)). Also in this case molecular quadrupoles and orientations play a crucial role discriminating between interfaces where e-h pairs remain strongly bound, and others where they can dissociate spontaneously or almost so. For instance, Verlaak et al showed that a binding energy E 0.4 b ∼ ± eV should be overcome at the edge-on pentacene(1 0 0)/C 60 (0 0 1) interface, while the charge separation is nearly barrierless for the face-on pentacene(0 1 1)/C 60 (0 0 1), as shown in figure 14(f). Similar trends in the energy profiles have been also reported by Van Voorhis and coworkers in their studies on quadrupolar and dipolar fields at D/A interfaces [49,111].
The polymer/fullerene blend between poly(3-hexylthiophene) (P3HT) and PC 61 BM is probably the most investigated organic photovoltaic system, typically as bulk heterojunctions (see figures 11(c) and (g)) obtained from solution processing. D'Avino et al applied ME calcul ations, assuming ten monomer-long polymer segments as D units, to map the dense manifold on e-h states at a dis ordered interface between a P3HT crystallite and amorphous PC 61 BM obtained with molecular dynamics simulations [105]. The resulting energy landscape for e-h states at this interface is shown in figure 16, partitioned into the different contributions according to equation (21). This study has the merit of estimating critical materials parameters, such as the e-h capture radius r 3 c ∼ nm, corresponding to the distance at which charges can be considered as free, and the average energy barrier for charge separation of E 0.3 b ∼ ± eV, the latter being substantially reduced by the quadrupolar fields of polymer chains with face-on orientation [105]. Both values are in very good agreement with experimental estimates for different polymer-fullerene heterojunctions [102]. Moreover, as a consequence of the large energetic disorder, mostly sourced from the dipole moments of randomly oriented PC 61 BM molecules, approximately half of the interfacial e-h states are able to separate barrierless [105].
We finally remark that an explicit calculation of bulk heterojunctions with these approaches remains an open challenge that is still out of the reach of atomistic techniques. Nevertheless, the qualitatively similar results obtained for rough bilayers, or with clusters comparable in size to D and A domains, support the idea that induction and electrostatic interactions, acting in synergy with disorder, represent the driving force for thermally accessible charge separation pathways.

Conclusions and perspectives
In this topical review we provided a theoretician's perspective on the most recent advances in the understanding of longrange electrostatic phenomena in organic semiconductors. Modern computational approaches based on atomistic models parametrized from first principles are able to capture, in many cases quantitatively, subtle effects arising from the anisotropy of π-conjugated molecules, such as the strong dependence of ionization energies of molecular films on supramolecular organization. The same electrostatic phenomena have profound implications in the context of organic photovoltaics, where charge carriers energetics is key. Major advances in this context are represented by the accurate calculation of the open circuit voltage of solar cells and the rationalization of how microscopic fields at donor-acceptor interfaces can prompt the separation of tightly bound carriers.
Open challenges and interesting routes to pursue still exist, especially concerning the realistic description of dis ordered films on rough substrates, relaxed grain boundaries, molecularly intermixed D/A systems, and bulk heterojuntions. In these cases and in all other where atomistic structural information is generally not available from experiments, high quality molecular dynamics simulations become of highest importance for predicting reasonable morphologies and, consequently, electronic properties.
On the methodological side a promising perspective is represented by hybrid quantum/classical approaches merging state of the art many-body electronic structure methods based on the GW formalism and the Bethe-Salpeter equation, which are receiving growing attention by the molecular science community because of their accuracy and good computational performances [112], with classical discrete polarizable models [113,114].
We conclude by recalling that electrostatic and induction interactions profoundly impact the energetics of charge carriers in a way that is highly sensitive to the molecular anisotropy and to the supramolecular organization at the mesoscale. These type of effects have been largely overlooked in the past, Figure 16. Electron-hole polarization energy and its different contributions according to equation (21) computed with ME calculations for an interface between crystalline P3HT and an amorphous PC 61 BM obtained from molecular dynamics simulations [105]. The average P r eh ( ) ± profile and energetic disorder fluctuations are shown as white bullets and error bars, respectively. Adapted with permission from [105]. Copyright 2013 American Chemical Society. and we believe this to be one of the reasons of the often elusive structure-property relationships in organic materials. We hence hope that this topical review will contribute to raise awareness on electrostatic phenomena, which should be taken into account in the rational design of materials and devices.