# Disorder-Induced Transition from Transient Quantum Delocalization to Charge Carrier Hopping Conduction in a Nonfullerene Acceptor Material

Ljiljana Stojanović, Jack Coker, Samuele Giannini, Giacomo Londi, Jens

Wenzel Andreasen, Jun Yan, Gabriele d'Avino, David Beljonne, Jenny Nelson, Jochen Blumberger

► To cite this version:

Ljiljana Stojanović, Jack Coker, Samuele Giannini, Giacomo Londi, Jens Wenzel Andreasen, et al.. Disorder-Induced Transition from Transient Quantum Delocalization to Charge Carrier Hopping Conduction in a Nonfullerene Acceptor Material. Physical Review X, 2024, 14 (2), pp.021021. 10.1103/physrevx.14.021021. hal-04800970

HAL Id: hal-04800970

https://hal.science/hal-04800970v1

Submitted on 24 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Disorder-Induced Transition from Transient Quantum Delocalization to Charge Carrier Hopping Conduction in a Nonfullerene Acceptor Material

Ljiljana Stojanović,¹ Jack Coker,² Samuele Giannini,³ Giacomo Londi,³ Anders S. Gertsen,⁴

Jens Wenzel Andreasen,⁴ Jun Yan,⁵,² Gabriele D'Avino,⁶ David Beljonne,³

Jenny Nelson,²* and Jochen Blumberger¹,†

¹Department of Physics and Astronomy and Thomas Young Centre, University College London,
London WC1E 6BT, United Kingdom

²Department of Physics, Imperial College London, United Kingdom

³Laboratory for Chemistry of Novel Materials, University of Mons, Mons 7000, Belgium

⁴Department of Energy Conversion and Storage, Technical University of Denmark,
Fysikvej 310, DK-2800 Kgs. Lyngby, Denmark

⁵School of Science and Engineering, The Chinese University of Hong Kong,
Shenzhen, Guangdong Province 518172, People's Republic of China

⁶Grenoble Alpes University, CNRS, Grenoble INP, Institut Neel, 38042 Grenoble, France

<figure><img src="image_2.png" /></figure>

(Received 3 August 2023; revised 13 November 2023; accepted 12 March 2024; published 29 April 2024)

Nonfullerene acceptors have caused a step change in organic optoelectronics research but little is known about the mechanism and factors limiting charge transport in these molecular materials. Here a joint computational-experimental investigation is presented to understand the impact of various sources of disorder on the electron transport in the nonfullerene acceptor O-IDTBR. We find that in single crystals of this material, electron transport occurs in the transient quantum delocalization regime with the excess charge delocalized over about three molecules on average, according to quantum-classical nonadiabatic molecular-dynamics simulations. In this regime, carrier delocalization and charge mobility ($μ_a$ = 7 cm² V⁻¹ s⁻¹) are limited by dynamical disorder of off-diagonal and diagonal electron-phonon coupling. In molecular assemblies representing disordered thin films, the additional static disorder of off-diagonal electron-phonon coupling is sufficient to fully localize the excess electron on single molecules, concomitant with a transition of transport mechanism from transient quantum delocalization to small polaron hopping and a drop in electron mobility by about 1 order of magnitude. Yet, inclusion of static diagonal disorder resulting from electrostatic interactions arising from the acceptor-donor-acceptor (A-D-A) structure of O-IDTBR, are found to have the most dramatic impact on carrier mobility, resulting in a further drop of electron mobility by about 4–5 orders of magnitude to 10⁻⁵ cm² V⁻¹ s⁻¹, in good agreement with thin-film electron mobility estimated from space-charge-limited-current measurements. Limitations due to diagonal disorder caused by electrostatic interactions are likely to apply to most nonfullerene acceptors. They imply that while A-D-A or A-DAD-A motifs are beneficial for photo-absorption and exciton transport, the electrostatic disorder they create can limit carrier transport in thin-film optoelectronic applications. This work shows the value of computational methods, in particular, nonadiabatic molecular-dynamics propagation of charge carriers, to distinguish different regimes of transport for different types of molecular packing.

DOI: 10.1103/PhysRevX.14.021021

Subject Areas: Chemical Physics, Materials Science, Physical Chemistry

## I. INTRODUCTION

The development of nonfullerene acceptor (NFA) materials marked a milestone in organic photovoltaics (OPV) research [1–12]. They alleviated the issue of large voltage losses in traditional fullerene-based junctions resulting in record power conversion efficiencies of currently >19% [13]. NFAs are typically composed of alternating electron-accepting and electron-donating groups (A-D-A or A-DAD-A), which results in beneficial photophysical

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*jenny.nelson@imperial.ac.uk</footnote>

<footnote>†j.blumberger@ucl.ac.uk</footnote>
properties. Several materials have recently been reported that support exciton diffusion that is longer ranged (20–47 nm) [14] than in classical organic semiconductors (OSs, 5–20 nm). This makes NFAs attractive for light-emitting devices, emission-based sensors, and luminescent solar concentrators because long diffusion lengths are usually associated with high radiative efficiency.

While optical properties of several NFAs have been characterized in detail, their charge-transport properties are less well studied by comparison. In typical applications, NFAs form disordered thin films that often exhibit poor electron transport, which may limit the efficiency in (opto) electronic applications. NFAs with higher charge mobility could further enhance the dissociation efficiency of excitons and suppress recombination reactions. A better understanding of charge transport in these materials is thus indispensable for rational design of improved NFAs. This starts with a fundamental understanding of the nature of charge carriers and their transport mechanism in single-crystal and thin-film samples. Do they form small polarons or transiently delocalized states as recently proposed for high-mobility OSs? What is the role of the various sources of disorder in thin films? What is the impact of electrostatics arising from alternating acceptor-donor structure on carrier mobility?

Transient localization theory (TLT) and quantum-classical nonadiabatic molecular-dynamics simulations have now established a new picture of charge transport in high-mobility organic crystals that are typically composed of relatively apolar (π)-conjugated molecules [15–31]. Charge carriers may delocalize in these materials fairly extensively and move along the crystal by a transient delocalization mechanism in a scenario that is “between” the two limiting cases that is often assumed, small polaron hopping of a fully localized charge carrier and band transport of delocalized Bloch waves. Yet, TLT and nonadiabatic molecular-dynamics simulations have mostly been carried out on archetypical apolar and rigid organic molecular crystals, typically acenes, rubrene, and thiophene-based compounds [32–34]. It is unclear if the transient delocalization mechanism carries over to other structures such as NFAs where electrostatic, torsional, and conformational disorder is expected to be much stronger than in the above materials.

To address these questions, we investigate here one of the first successful NFAs, 5,5'-[(4,9-dihydro-4,4,9,9-tetraoctyl-s-indaceno[1,2-b:5,6-b']dithiophene-2,7-diyl)bis(2,1,3-benzothiadiazole-7,4-diylmethylidyne)]bis[3-ethyl-2-thioxo-4-Thiazolidinone (O-IDTBR)] [35]. This molecule is composed of a central electron-rich indacenodithiophene (IDT) group connected to electron-accepting benzothiadiazole (BT) and rhodanine (RH) groups forming an acceptor-donor-acceptor motif (A-D-A); see Fig. 1. The RH-BT and BT-IDT groups are coplanar, and the molecules arrange in an interdigitated columnar structure according to single-crystal x-ray crystallography [36] (a more detailed description of

<figure><img src="image_3.png" /><figcaption>FIG. 1. Chemical structure of O-IDTBR. The RH is presented in blue, the BT in green, and the IDT in orange; $R = \text{octyl}$. Dihedral angles connecting the IDT and BT groups are denoted $\gamma_1$ and $\gamma_2$.</figcaption></figure>

packing structure is given in Sec. II. A). The potential energy barrier along the dihedral angles connecting the IDT and BT groups, $\gamma_1$ and $\gamma_2$ in Fig. 1, is rather shallow (8.6 kcal/mol) [37,38] compared to some high-performance NFAs, e.g., Y6 and ITIC [39]. We thus expect larger torsional or conformational disorder than in the latter molecules. Moreover, the presence of a number of heteroatoms in the three constituent groups gives rise to local dipole and higher multipole moments as well as electrostatic disorder, in addition to the structural disorder present in thin films. O-IDTBR is therefore an excellent model system to investigate the effect of different types of disorder on the charge-transport mechanism.

In previous works, it was shown that O-IDTBR combined with a variety of donor polymers in bulk heterojunction blends enables power conversion efficiencies above 10% [40] and low recombination loss rates with high efficiencies with small molecule donors [41]. Apart from applications in OPV devices, O-IDTBR is also suitable as an *n*-type material in organic thin-film transistors [36]. Some basic charge-transport properties of crystalline and disordered O-IDTBR systems were studied previously. Bristow *et al.* [36] reported field-effect-transistor (FET) measurements of electron mobilities and calculated electronic couplings of molecular dimers isolated from the crystal structure. Gertsen *et al.* parametrized a force field for O-IDTBR, generated model structures of disordered thin films, and discussed the dominant intermolecular interactions [37,38].

In this study, we carry out quantum-classical nonadiabatic molecular-dynamics simulations on O-IDTBR, specifically, fragment-orbital-based fewest-switches surface-hopping FOB-SH [25,27,42], as well as kinetic Monte Carlo (KMC) simulations [43] to propagate an excess electron in single crystals and in molecular-dynamics-generated amorphous assemblies representing models of thin films. This is complemented by space-charge-limited-current (SCLC) measurements of electron mobility in thin films. Particular attention is paid to the different types of disorder in the system, structural and electrostatic, diagonal and off-diagonal, and how they affect the electron delocalization and propagation through the material. Our results suggest that in single crystals electronic coupling is sufficiently strong that excess electrons delocalize, albeit modestly, with electron
mobilities rivaling those of high-charge-mobility OSs. In thin-film models, the calculated electron mobility drops by 5–6 orders of magnitude compared to the single crystal. The static off-diagonal disorder is sufficient to induce complete localization of carriers, but this explains only a small fraction of the decrease in mobility. The main loss mechanism is the strong electrostatic disorder in the thin film.

This paper is organized as follows. In the following Sec. II A, we describe the O-IDTBR crystal structure and present the electronic coupling network. In Sec. II B, we present and discuss the charge-transport mechanism in single crystals obtained from nonadiabatic molecular-dynamics simulations. In Sec. II C, simulation results are presented for structural models of disordered thin films and compared to experimental SCLC measurements. Particular focus is placed on the analysis of the various sources of disorder in the thin film on the transport mechanism and electron mobility. After conclusion of this work in Sec. III, we explain the methods and give simulation and experimental details in Sec. IV. More technical details and further analysis results are included in the Supplemental Material [44].

## II. RESULTS AND DISCUSSION

### A. Crystal structure and electronic couplings

The 3D structure of single-crystal O-IDTBR is visualized in Fig. 2 by two projections, one onto the *a*-b and another onto the *a*-c planes. The molecules arrange in an interdigitated columnar crystal structure with four molecules in a monoclinic unit cell (*M*1–*M*4). The four dimers with the highest electronic couplings (for excess electron transfer) are denoted *D*1–*D*4 in Fig. 2. They interact via *π*–*π* stacking between their acceptor (BT and RH) units. In the *D*2 (*M*2*M*3) and *D*3 (*M*1*M*4'), dimers the longest axes of the monomers are parallel, whereas in the *D*1 (*M*1*M*2) and *D*4 (*M*2*M*4' or *M*1*M*3) dimers, the monomers are rotated in plane with respect to each other. In the notation used, primed labels refer to molecules in an adjacent periodic image cell, so *M*1*M*4' and *M*2*M*4' cross a unit-cell boundary. The *D*1 and *D*2 dimers have the highest-magnitude electronic couplings in the ideal crystal structure, *H*<sub>*kl*</sub>(*D*1) = −61 meV and *H*<sub>*kl*</sub>(*D*2) = −55.7 meV. *D*1 dimers form a continuous pathway along the *a* direction connecting adjacent unit cells, and they also contribute to the charge transport along the *b* direction but only within a single unit cell. Similarly, the parallel slip-stacked *D*2 dimer contributes mainly to the charge transfer in the *a* direction and *b* direction within a single unit cell. The *D*3 dimer, a slip-stacked dimer composed of two neighboring molecules from adjacent unit cells, features significantly smaller couplings [*H*<sub>*kl*</sub>(*D*3) = −12.3 meV] than *D*1 and *D*2 and is the weakest electronic connection along the *b* direction. Hence, this dimer is expected to limit charge mobility in this direction. The *D*4 dimer is a rotated slip-stacked dimer composed of molecules which are

<figure><img src="image_4.png" /><figcaption>FIG. 2. Crystal structure of O-IDTBR. The unit cell is indicated in thick lines and is comprised of four molecules denoted *M*1–*M*4. Molecules in neighboring unit cells are annotated by a prime. Two crystallographic projections are shown, one onto the *a*-b plane (left panels) and one onto the *a*-c plane (right panels). Dimers with the largest electronic couplings denoted *D*1–*D*4 are shown in color. The lowest unoccupied molecular orbital (LUMO) for a O-IDTBR molecule is shown in the top panel. The structure was taken from Ref. [36].</figcaption></figure>

second-nearest neighbors, with rather small and potentially negligible electronic couplings $H_{kl}(D4) = -0.5$ meV in the *a*-c plane. All other dimers have even smaller couplings. Thus, the electronic coupling network is very anisotropic in O-IDTBR. The highest mobility is expected along the *a* direction due to uninterrupted pathways formed by strongly coupled *D*1 dimers. Smaller mobility values are likely along the *b* direction due to the interruption of the strongly coupled *D*1 and *D*2 pairs by weakly coupled *D*3 dimers. There are no dimers with significant couplings in the *c* direction, implying much smaller mobility values along this direction than in the *a*-b plane. Thus, we perform FOB-SH charge-transport simulations only for the *a*-b plane (see Sec. II B below).

The transport dynamics are determined not only by electronic couplings but, crucially, also by the magnitude of their thermal fluctuations. We extract the coupling distributions for the *D*1–*D*3 dimers from the FOB-SH
<figure><img src="image_5.png" /><figcaption>FIG. 3. Distributions of electronic couplings ($H_{kl}$ in meV) for D1–D3 dimers in the O-IDTBR single crystal. Couplings obtained for the crystal structure [36] shown in Fig. 2 are indicated in solid vertical lines, and average values obtained from FOB-SH simulation at 300 K are shown in dashed vertical lines with root-mean-square fluctuations reported in the figure legend (in meV).</figcaption></figure>

dynamics at 300 K (Fig. 3) and compute averages $\langle\langle H_{kl} \rangle\rangle$, root-mean-square averages $\langle\langle H_{kl}^2 \rangle\rangle^{1/2}$, and the corresponding root-mean-square fluctuations. D1 couplings exhibit a very broad bimodal distribution ranging from –200 to 100 meV giving rise to very large coupling fluctuations that are almost as large as the average coupling $\langle\langle H_{kl} \rangle\rangle = –78.5 \pm 55.0$ meV, $\langle\langle H_{kl}^2 \rangle\rangle^{1/2} = 96 \pm 84$ meV). The broad distribution and the change in coupling sign from negative to positive for a number of D1 configurations is due to a change in the relative orientations of BT and IDT groups of the monomers, induced by thermal motions of the corresponding dihedrals. These motions give rise to large off-diagonal disorder that can be expected to be detrimental to charge carrier delocalization and mobility. Thermal fluctuations

also broaden the D3 coupling distributions markedly ($\langle\langle H_{kl} \rangle\rangle = 0.2 \pm 33.0$ meV, $\langle\langle H_{kl}^2 \rangle\rangle^{1/2} = 32.5 \pm 38.0$ meV). By comparison, the distribution of D2 couplings is rather narrow with the mean couplings close to the value in the crystal ($\langle\langle H_{kl} \rangle\rangle = –44 \pm 14$ meV and $\langle\langle H_{kl}^2 \rangle\rangle^{1/2} = 46 \pm 32$ meV).

## B. Electron transport in single crystals

Electron-transport simulations are carried out with FOB-SH [25,27,42]. In this method, the valence- or conduction-band states of an organic semiconductor are described by an electronic Hamiltonian represented in a molecular orbital or site basis; Eq. (1) in Sec. IV. The wave function of an excess electron in the crystal, Eq. (2), is propagated in time according to the time-dependent Schrödinger equation, Eq. (3), coupled to the thermal motion of the nuclei. We refer to Sec. IV for further details of the method. The carrier wave function is chosen to be initially localized on one of the O-IDTBR molecules in the center of the box (index i), $\Psi(t=0) = \phi_i(t=0)$, and propagated in time for 1 ps. Within the first 200 fs of dynamics, the mean-square displacement (MSD) of the wave function averaged over all trajectories increases sharply and in a nonlinear fashion [Figs. 4(a) and 4(b)]. This is because the initially localized wave function is a linear combination of low- and high-lying conduction-band states that relaxes to low-lying states close to the bottom of the conduction-band edge. This relaxation process is accompanied by wave-function delocalization over about three molecules, $\langle\langle \text{IPR} \rangle\rangle = 3.3$ [Fig. 4(c)]. Thereafter, the MSD increases linearly with time ($R^2 = 0.96$) indicative of Einstein diffusion, and the average inverse participation ratio (IPR), $\langle\langle \text{IPR} \rangle\rangle$, fluctuates stably around 2.7. Very similar carrier dynamics were previously reported for molecular crystals of acenes [45]. We note that other choices of wave-function initialization

<figure><img src="image_6.png" /><figcaption>FIG. 4. Dynamics of the excess electron wave function in the O-IDTBR crystal. The MSD of the wave function is shown along the crystallographic direction *a* as obtained from FOB-SH using a 10 × 3 simulation cell [panel (a), Eq. (4)], $\text{MSD}_{\alpha\beta} = \text{MSD}_{aa} = \text{MSD}_a$ and along the *b* direction as obtained from FOB-SH using a 5 × 5 simulation cell [panel (b), Eq. (4)], $\text{MSD}_{\alpha\beta} = \text{MSD}_{bb} = \text{MSD}_b$. Standard deviations of the MSD due to averaging over FOB-SH trajectories are indicated by error bars. Linear fits of the MSDs are shown in thick lines and used for the calculation of diffusion constants and electron mobilities according to Eqs. (5) and (6). The corresponding trajectory-averaged IPRs [integrand of Eq. (8)] are shown in panel (c). The mean IPR values obtained by averaging over times $t > 200$ fs according to Eq. (8) are indicated by dashed horizontal lines.</figcaption></figure>
will result in different short-time relaxation dynamics in the first few 100 fs but in the same Einstein diffusion dynamics at longer times.

The computed MSDs are used to obtain the diffusion and mobility tensors according to Eqs. (5) and (6) in Sec. IV. The eigenvalues are highly anisotropic; the mobility along the *a* direction (approximately parallel to the long molecular axis) is $\mu_a = 7.0 \pm 0.6 \text{ cm}^2 \text{ V}^{-1} \text{ s}^{-1}$, approximately 20 times larger than the mobility along the *b* direction, $\mu_b = 0.30 \pm 0.03 \text{ cm}^2 \text{ V}^{-1} \text{ s}^{-1}$. This finding is in line with the anisotropy of electronic coupling values along *a* and *b* directions, as discussed above. It also supports a previous finding that the anisotropy of mobility increases with the crystallinity of disordered O-IDTBR samples with a maximum ratio of 10 between the two $\mu$ tensor eigenvalues [46]. For single-crystal systems, this ratio can be even larger, about 20 according to our FOB-SH simulations.

According to FOB-SH simulations, wave-function delocalization in O-IDTBR crystals happens mostly along the *a* direction, and this is because the dominant electronic coupling in this direction [$\langle\langle H_{kl}^2 \rangle\rangle^{1/2} = 96 \text{ meV, D1 pair}$] exceeds half of the intramolecular reorganization energy as obtained from the four-point method ($\lambda_{\text{intra}} = 180 \text{ meV}$). Note that if the coupling was much smaller than the reorganization energy, FOB-SH would yield a stable small polaron as the electronic ground state that is stabilized relative to the undistorted molecule by half of the reorganization energy. However, this is not the case here; the excess electron delocalizes, implying that small polaron hopping is not the physically correct model to describe electron transport in the crystal. Yet, the extent of electron delocalization in O-IDTBR is fairly small compared to, e.g., rubrene ($\langle\langle \text{IPR} \rangle\rangle \approx 13$) even though the two materials exhibit similar electronic coupling strengths in their high-mobility (*a*) directions. The difference is that in O-IDTBR the electronic coupling fluctuations induced by dihedral motion of the BR and IDT groups are much larger than the thermal fluctuations in rubrene, causing larger off-diagonal disorder that prevents more extensive electron delocalization.

The electron transport in O-IDTBR crystals proceeds via transient delocalization events, as exemplified in Fig. S6 in the Supplemental Material [44]. In the snapshot at 910 fs, the excess electron is localized over about two molecules, which is close to the average polaron size. At this point in time, a thermal excitation to an energetically higher-lying and more delocalized conduction-band electronic state occurs via a surface hop, resulting in an expansion of the polaron along the *a* direction from two to about six molecules at 920 fs. This is followed by a surface hop to a lower-lying conduction-band state concomitant with a contraction of the polaron to four molecules at 940 fs. The expansion and contraction events result in a shift of the center of charge of the polaron and thus to charge transport. These events occur very frequently, approximately every 40 fs and typically last less than 20 fs. In the presented

trajectory, the polaron is observed to expand transiently over up to about 12 molecules, extending over seven unit cells (approximately 20 nm) in the *a* direction, and only over one unit cell in the *b* direction. As noted above, the delocalization in the *a* direction is supported by a continuous network of strongly coupled D1 and D2 dimers. By contrast, delocalization along the *b* direction is prevented by the weakly coupled D3 dimer, thus explaining the strong mobility anisotropy.

A caveat of our FOB-SH simulation is that the site energies, which are parametrized to reproduce density functional theory (DFT) intramolecular reorganization energy of the molecule in vacuum (see Sec. IV for details), do not contain intermolecular electrostatic energy contributions. This is because the computation of the full Ewald electrostatics for each charge state and for each nuclear time step in nonadiabatic FOB-SH simulations (0.025 fs) is not computationally feasible. Standard MD simulations of the O-IDTBR crystal including Ewald electrostatics gives a time-averaged intermolecular electrostatic energy contribution to the site energy of $84 \pm 69, 87 \pm 70, 80 \pm 70$, and $86 \pm 70 \text{ meV}$ for the four molecules in the unit cell of O-IDTBR, M1 to M4. Thus, the time-averaged electrostatic potential that the excess electron feels when localized on any one of the four molecules is the same within the statistical error bar. This means that the static electrostatic disorder in the O-IDTBR crystal is virtually zero. Hence, neglecting the electrostatic site energy contribution in FOB-SH dynamics does not introduce errors in static electrostatic disorder. However, the dynamic contribution to electrostatic disorder, which we estimate in Sec. II C below, corresponds to an estimated intermolecular reorganization energy $\lambda_{\text{inter}} = 80 \text{ meV}$ (Table II), and this contribution is neglected in FOB-SH simulations. Here we note that this contribution is due to all lattice vibrations, fast and slow. Yet, only the part of the intermolecular reorganization energy that is due to the fast lattice vibrations occurring on the timescale of transient delocalization events ($< 100 \text{ fs}$) will increase the effective diagonal disorder and couple to the carrier dynamics in FOB-SH. While we expect that the fast contribution to the intermolecular reorganization energy is small leading to only a slight reduction in carrier mobility, electrostatics-embedded FOB-SH would be required to fully quantify this effect.

### C. Electron transport in thin film

An amorphous assembly representing a model of a disordered thin film was created by nonequilibrium solvent evaporation MD simulations on an $\text{SiO}_2$ substrate in previous work [38]. This structural model includes the types of disorder likely to be present in a solution-processed nonequilibrium system. The solvent environment enabled large-amplitude rotations of dihedral angles $\gamma_1$ and $\gamma_2$ connecting BT and IDT units resulting in bimodal distributions of molecular conformation in the solid state
<figure><img src="image_7.png" /><figcaption>FIG. 5. Distribution of dihedral angles in the thin-film model structure. The bimodal distributions for the dihedrals $\gamma_1$ and $\gamma_2$ connecting IDT and BT groups [defined in Fig. 1 and again indicated in panel (a)] are shown in red and gray, respectively, for 300 K. Mean values are indicated by vertical lines. For comparison, the corresponding values for the experimental crystal structure [36] are shown in black vertical lines. A typical dimer structure in the thin-film model structure is depicted in panel (b).</figcaption></figure>

(Fig. 5). The more intense peaks of the distributions peak at 180°, in contrast to the equilibrium values in the crystal, close to 0°. The broad distributions indicate that most molecules are in nonplanar conformations, i.e., have marked backbone distortions. The most common type of electronically coupled dimers are interacting through BT-BT $\pi$-$\pi$ stacking (Fig. 5), not seen in the crystal structure. Dimers interacting via BT-RH units as in the crystal structure, in particular D2 and D4 with parallel, and D1 and D3 with rotated longitudinal axes, are very rare.

The bimodal dihedral distributions and the diverse orientations of neighboring molecules create strong static off-diagonal and diagonal disorder in the electronic Hamiltonian Eq. (1). We first focus on the former disorder. We find that the dimers are much weaker coupled than in the crystal. For most dimers, $H_{kl} < \lambda_{\text{intra}}/10$, some have larger couplings, $\lambda_{\text{intra}}/10 < H_{kl} < \lambda_{\text{intra}}/2$, and no pairs are sufficiently strongly coupled to induce carrier delocalization ($H_{kl} > \lambda_{\text{intra}}/2$). We characterize electronic couplings in the thin-film model structure in terms of an electronic coupling network as shown in Fig. 6. The centers of mass of neighboring molecules with couplings

<figure><img src="image_8.png" /><figcaption>FIG. 6. Electronic coupling map for the thin-film model structure. (a) Dimers are connected by a red line if the electronic coupling is large, in the range $\lambda_{\text{intra}}/10 \le H_{kl} \le \lambda_{\text{intra}}/2$. (b) Distribution of continuous coupling paths of length $l$ where $l$ is the number of molecules with electronic couplings to their adjacent neighbors in the range $\lambda_{\text{intra}}/10 \le H_{kl} \le \lambda_{\text{intra}}/2$. The lengths of the shortest, mean, and longest coupling paths in each MD-generated structure ($l_S, l_M$, and $l_L$) are averaged over all 100 MD structures that are used as initial configurations for FOB-SH simulation (averages indicated by vertical lines).</figcaption></figure>

exceeding $\lambda_{\text{intra}}/10$ are connected by a red line. These molecular pairs are expected to contribute most to charge transport. The coupling network spans in three directions and is composed of a set of graphs of variable length. Here, a graph of length $l$ is defined as a continuous path of $l$ molecules where the electronic coupling between next-nearest neighbors exceeds $\lambda_{\text{intra}}/10$. To include thermal effects, we run a 1-ns MD simulation at 300 K, 0.5 ns in the canonical (NVT) ensemble followed by 0.5 ns in the microcanonical (NVE) ensemble, starting from the structure of Ref. [38]. We then extract the graphs for 100 geometries sampled equidistantly from the last 100 ps of the trajectory. The distribution of the graph lengths, i.e., number of molecules forming a graph, is shown in Fig. 6(b). The minimum and maximum graph lengths averaged over all 100 structures are $\langle l_S \rangle = 2$ and $\langle l_L \rangle = 92$, respectively, and the mean is $\langle l_M \rangle = 8$ (note, for the crystal $\langle l_S \rangle = \langle l_L \rangle = \langle l_M \rangle \to \infty$). Our analysis suggests that electron transport in the thin-film model will occur via small polaron hopping since $H_{kl} < \lambda_{\text{intra}}/2$, and that transport is limited by slow hops between weakly coupled molecules $H_{kl} < \lambda_{\text{intra}}/10$ connecting two different graphs.

FOB-SH simulations of charge transport in the thin-film model structure may be sensitive to where the excess electron is initially placed, in a region that is electronically well connected (i.e., long graph) or poorly connected (i.e., short graph). Thus, we carry out FOB-SH dynamics simulations for three sets of trajectories, where the initial wave function $\Psi(t=0) = \phi_i(t=0)$ is located on molecules $i$ that are part of short ($l=2$), medium-sized ($l=8-10$), and long graphs ($l=50-150$). For each case, 250-350 FOB-SH trajectories are run for 1 ps. We find that in all three cases the static off-diagonal disorder in the thin-film model is sufficient to keep the excess electron localized on a single molecule, except in hopping events where the electron transits between molecules, resulting in
a mean IPR of 1.1–1.2. Thus, the static off-diagonal disorder is sufficient to change the transport mechanism from transient delocalization to small polaron hopping.

As the thin-film model structure is grown on a substrate, it is of interest to quantify the MSD of the carrier wave function in the directions parallel to the substrate, denoted in the following $\text{MSD}_{xy}$, in the z direction orthogonal to the x-y plane, $\text{MSD}_z$, and in full 3D space, $\text{MSD}_{xyz}$; see Eq. (7) and Sec. IV for details. Note, we use $(x,y,z)$ to indicate the directions for the thin-film model structure because the crystallographic directions $(a,b,c)$ of the single crystal are no longer meaningful for the disordered thin film. All MSDs are shown in Fig. 7 for initializations of the wave function in regions of different electronic coupling strengths. Initiated from molecules in the middle of long graphs (denoted as $l_L$), $\text{MSD}_{xy}$ and $\text{MSD}_{xyz}$ increase rapidly in the first 200 fs, reaching approximately two-fold-larger values compared to the MSD from dynamics initiated from sites in medium and short graphs ($l_M$ and $l_S$). However, at longer simulation times when diffusive dynamics sets in, the MSD slopes are rather insensitive to the initialization. The in-plane, out-of-plane, and 3D mobility of the film are obtained from the corresponding MSDs using the Einstein relation for 1D, 2D, and 3D

transport, respectively, Eq. (7), and they are summarized in Table I for all three initial conditions. The mobility values are all very similar, in the range 0.3–0.4 $\text{cm}^2 \text{ V}^{-1} \text{ s}^{-1}$, indicating that the charge transport in the generated thin-film model is approximately isotropic. These values are 20- to 30-fold smaller than along the high-mobility direction in the crystal, and they are very similar to the mobility in the b direction of the crystal. This agreement is likely to be accidental since the crystal and thin-film model structures are very different.

In the following, we examine the transport mechanism in the thin-film model in more detail by considering a representative FOB-SH trajectory; Fig. S7 in the Supplemental Material [44]. According to the IPR($t$), the excess electron is mostly localized on a single molecule and occasionally expands to two to three molecules for short durations of time (<10 fs), often in response to a surface hop to a higher-lying band state. Only a small fraction of these events are successful, meaning the wave function becomes displaced to the next molecules after relocalization, contributing to mobility. More extensive and prolonged wave-function delocalizations, as frequently seen in the crystal, are very rare, and this is why the mobility in the thin film is much smaller than in the crystal. One such rare event is shown in

<figure><img src="image_9.png" /><figcaption>FIG. 7. Dynamics of the excess electron wave function in the structural model of an O-IDTBR thin film. MSDs of the wave function (upper panels) and inverse participation ratios (lower panels) are obtained from FOB-SH simulations where the excess electron is initially localized on a molecule in the thin-film model that is part of a short [$l_S$, panels (a) and (d)], mean [$l_M$, panels (b) and (e)], and long electronic coupling path [$l_L$, panels (c) and (f)]. The MSD is shown for the x-y plane parallel to the substrate ($\text{MSD}_{xy}$, data in blue), along the z direction orthogonal to the substrate ($\text{MSD}_z$, data in green), and in 3D ($\text{MSD}_{xyz}$, data in black); see Sec. IV for details on the calculations. Standard deviations of the MSD due to averaging over FOB-SH trajectories are indicated by error bars. Linear fits of the MSDs are shown in thick lines and used for the calculation of diffusion constants and electron mobilities according to Eqs. (7) and (6); see Table I for numerical values. The corresponding trajectory-averaged IPRs [integrand of Eq. (8)] are shown in black lines, and the mean IPR values obtained by averaging over times $t > 200$ fs according to Eq. (8) are indicated by horizontal red lines.</figcaption></figure>
<table><caption>TABLE I. Charge mobilities and average IPR values for the structural model of an O-IDTBR thin film from FOB-SH simulations. The mobilities ($\mu_z$, $\mu_{xy}$, and $\mu_{xyz}$ in cm$^2$ V$^{-1}$ s$^{-1}$) are obtained from linear fits of the MSDs shown in Figs. 7(a)–7(c) according to Eqs. (7) and (6). The R$^2$ values of the linear fits are given in parentheses. The average IPR values ($\langle \text{IPR} \rangle$) are obtained by averaging over times $t > 200$ fs according to Eq. (8). They are indicated by horizontal red lines in Figs. 7(d)–7(f). Values are given for three different initial conditions where the excess electron is initially localized on a molecule that is part of a short ($l_S$), mean ($l_M$), and long electronic coupling path ($l_L$).</caption><thead><tr><th></th><th><i>l</i><sub>S</sub></th><th><i>l</i><sub>M</sub></th><th><i>l</i><sub>L</sub></th></tr></thead><tbody><tr><td>μ<sub>z</sub></td><td>0.42 (0.78)</td><td>0.34 (0.94)</td><td>0.42 (0.96)</td></tr><tr><td>μ<sub>xy</sub></td><td>0.42 (0.89)</td><td>0.28 (0.96)</td><td>0.39 (0.98)</td></tr><tr><td>μ<sub>xyz</sub></td><td>0.30 (0.96)</td><td>0.24 (0.98)</td><td>0.28 (0.98)</td></tr><tr><td>⟨IPR⟩</td><td>1.15</td><td>1.15</td><td>1.21</td></tr></tbody></table>

Fig. S7(b) in the Supplemental Material [44]. At 750 fs, the wave function is localized on a single site, then delocalizes, first over two (at 770 fs) and then over four to five sites at 780 fs. In the subsequent time steps, the wave function remains delocalized, but the center of the charge is shifted, followed by relocalization of the charge on two sites (at 800 fs) and a single site (at 810 fs). Since these events are very rare and transport is dominated by smaller expansions described above, we refer to the transport scenario in the thin-film model as mostly small polaron hopping.

The FOB-SH thin-film simulations above include static and dynamic off-diagonal disorder, diagonal disorder from the intramolecular force-field terms, but no static or dynamic electrostatic diagonal disorder. For disordered thin films it can be expected that the electrostatic diagonal disorder is large and can no longer be neglected, in particular, for molecules with significant intramolecular variations of charge distribution and high conformational flexibility as is the case for O-IDTBR. Here we take advantage of the fact that the carrier dynamics in the O-IDTBR thin-film model are by small polaron hopping according to FOB-SH quantum dynamics, and the inclusion of electrostatic disorder will not change this condition (it will even further increase carrier localization). Thus, it is appropriate to model electron transport in the thin-film structure as a sequence of electron hopping events sampled with, e.g., KMC methods in combination with hopping rates from semiclassical electron-transfer theories such as Marcus theory. The electrostatic disorder can then be included in the calculation of the electron-transfer parameters determining the rate, in particular, the free-energy difference and the reorganization free energy.

To validate the KMC approach, we first calculate the mobility for the thin-film model without the electrostatic diagonal disorder and compare to the results of FOB-SH. In KMC, we use Marcus theory to define electron-transfer rates using thermally averaged electronic couplings and a

reorganization free energy that is set equal to the intramolecular contribution, as in FOB-SH, $\lambda = \lambda_{\text{intra}}$. The free-energy difference or driving force for electron transfer in this model is negligibly small and is set to zero for each pair because the fast intramolecular terms thermally average out to values very close to zero. Further details on the KMC simulations are given in Sec. IV. The resulting mobilities are found to be $\mu_{xy} = 0.31$ cm$^2$ V$^{-1}$ s$^{-1}$, in good agreement with $\mu_{xy}$ obtained from FOB-SH after averaging over initial conditions, $\mu_{xy} = 0.36$ cm$^2$ V$^{-1}$ s$^{-1}$. This suggests that the charge carrier wave-function dynamics obtained from FOB-SH can indeed be well approximated by a small polaron hopping model with standard Marcus rates, the regime assumed by our KMC approach. It also provides validation for the use of FOB-SH as a tool for ascertaining the regime of transport in a given system.

As a further calibration of the KMC mobility calculation, we use the same approach to calculate the mobility in the O-IDTBR single-crystal structure using thermally averaged electronic coupling values. Again, the reorganization energy is set equal to the intramolecular contribution, and no electrostatic disorder contributions are included in the site energies, similar to the FOB-SH simulations on the single crystal. The KMC-derived mobility values are not expected to be representative of single crystal because that system is not in the small polaron hopping regime (see Sec. II B above), and indeed the KMC-derived mobilities overestimate the FOB-SH single-crystal mobilities by a factor of 3 (see Fig. S9 in the Supplemental Material [44]). However, the orientation dependence obtained in the FOB-SH simulations of the crystal is preserved in the KMC values.

Next, we investigate the electrostatic disorder in the thin-film model structure by way of a classical microelectrostatic (ME) model, which some of us have previously developed [47,48]. We use the same terms and notations as in the original ME studies; in particular, we refer to the electrostatic contribution to the site energy as “polarization” contribution $\Delta^{-}$. The latter energy is composed of three terms: the fixed-point charge ($\Delta_{\bar{E}}^{-}$), induced-point dipole ($\Delta_{\bar{\text{In}}}^{-}$) and induction contributions ($\Delta_{\bar{\text{Ie}}}^{-}$); see Eqs. (15)–(18) in Sec. IV for their definitions. We find that the fixed-point charge and dipole contribution are relatively small, while the induction term has the largest contribution to the average polarization energy $\langle \Delta^{-} \rangle$; see Table II. Yet, the relevant property affecting electron transfer rates is the electrostatic disorder of the site energies $\sigma_{\text{tot}}$ defined as the root-mean-square fluctuations of the polarization energies averaged over time and all sites, Eq. (19). We note that the electrostatic disorder of the site energies is a collective effect of all molecules in the crystal due to the long-range nature of this interaction and cannot be related to some local structural parameters. We find that the largest contribution to $\sigma_{\text{tot}}$ is due to the fixed-point charge contribution ($\Delta_{\bar{E}}^{-}$), while inclusion of induced dipoles and induction leads to a reduction in the
<table><caption>TABLE II. Summary of classical ME calculations for an excess electron in a thin-film model structure of O-IDTBR. The polarization energy contribution to the site energy [Δ⁻, Eq. (15)] is broken down in fixed-point charge [Δ⁻*, Eq. (16)], induced-point dipole [Δ⁻*, Eq. (17)], and induction contributions [Δ⁻*, Eq. (18)], in accord with Refs. [47,48]. The average site energy or site energy contributions are denoted ⟨X⟩, X = Δ⁻, Δ⁻*, Δ⁻*, Δ⁻*, and the corresponding root-mean-square fluctuations [σtot, Eq. (19)] are broken down in static [σstat, Eq. (20)] and dynamic [σdyn, Eq. (21)] contributions. σcross, Eq. (22), is negligibly small. The polarization energy contribution to the site energies at the energy minimum of a negatively charged and neutral molecule in the crystalline solid are denoted X⁻ and X₀, respectively, and the corresponding intermolecular or external reorganization energy obtained from these values is denoted λinter. All energies are in eV.</caption><thead><tr><th></th><th>⟨X⟩</th><th>σ<sub>tot</sub></th><th>σ<sub>stat</sub></th><th>σ<sub>dyn</sub></th><th>X<sub>-</sub></th><th>X<sub>0</sub></th><th>λ<sub>inter</sub></th></tr></thead><tbody><tr><td>Δ<sub>E</sub><sup>-</sup></td><td>-0.231</td><td>0.209</td><td>0.198</td><td>0.067</td><td>-0.345</td><td>-0.301</td><td>0.088</td></tr><tr><td>Δ<sub>E</sub><sup>-</sup> + Δ<sub>In</sub><sup>-</sup></td><td>-0.161</td><td>0.130</td><td>0.114</td><td>0.063</td><td>-0.475</td><td>-0.439</td><td>0.071</td></tr><tr><td>Δ<sub>Te</sub><sup>-</sup></td><td>-0.801</td><td>0.096</td><td>0.095</td><td>0.016</td><td>-0.696</td><td>-0.691</td><td>0.010</td></tr><tr><td>Δ<sup>-</sup></td><td>-0.962</td><td>0.172</td><td>0.159</td><td>0.065</td><td>-1.171</td><td>-1.130</td><td>0.080</td></tr></tbody></table>

electrostatic disorder. This is in accord with many previous studies that have shown that inclusion of electronic polarization decreases the energy fluctuations for charge insertion or removal as well as reorganization energies [49–52], as we show further below.

Here, we further analyze σtot in terms of the static contribution σstat, dynamic contribution σdyn, and cross contribution σcross, where σ²tot = σ²stat + σ²dyn + σ²cross. The static contribution is due to the variation of the time-averaged site energies across the different sites with respect to the site energy averaged over all times and sites, Eq. (20). It originates from the packing disorder in the thin film. The dynamic contribution is due to the time variation of the site energy of a single site, averaged over all sites, Eq. (21). The latter energy is due to the intermolecular and conformational thermal motion of surrounding molecules. The results of this analysis are summarized in Table II. We find that the electrostatic disorder in the thin film is dominated by the static contribution, whereas this contribution is equal to zero in the single crystal (see also discussion at the end of Sec. II B). The dynamic contribution is relatively small by comparison indicating that thermal motion is fairly restricted by the surrounding molecules.

The intermolecular or external reorganization energy for electron transfer between two O-IDTBR molecules (λinter, Table II) is obtained as twice the polarization energy difference between the state where a molecule in the solid is negatively charged in the optimized geometry of the negatively charged state (X⁻ in Table II) and the state where the same molecule is negatively charged but in the geometry optimized for the neutral state (X₀; see Sec. IV for details of this calculation). The polarization energies are again decomposed into the three contributions as above. The total value λinter = 80 meV is due to a small rearrangement of the environmental charge density in response to a change of the intramolecular geometry from neutral to charged. It originates mostly from fixed-point charge and induced-dipole contribution (71 meV), whereas

the induction contribution is smaller by comparison (9 meV). Thus, the intermolecular reorganization energy (80 meV) is significantly smaller than the intramolecular contribution (180 meV), as is usually assumed for organic semiconductors, but it is not negligible resulting in a total reorganization energy of 260 meV.

The electrostatic disorder described above is accounted for in KMC simulations by including the polarization energy in the calculation of the site energies from which the driving force and reorganization energy are obtained that enter the Marcus rates; see Sec. IV for details. This inclusion has a major effect on the electron mobilities. As shown in Fig. 8, the diagonal disorder introduced by the polarization energy causes the mobility predictions to decrease by several orders of magnitude. The largest drop results from the ΔE contribution, the inclusion of which results in mobilities μx = 1.9 × 10⁻⁷ cm² V⁻¹ s⁻¹ and μy = 5.5 × 10⁻⁷ cm² V⁻¹ s⁻¹. Adding the ΔIₙ term acts to partially smooth out the

<figure><img src="image_10.png" /><figcaption>FIG. 8. Summary of computed and experimental electron mobilities for O-IDTBR. Notice the large drop in electron mobility upon inclusion of electrostatic disorder in the thin film.</figcaption></figure>
energetic landscape leading to an increase in mobility, while the $\Delta I_e$ term reduces it again. When all polarization contributions are included, mobilities of $\mu_x = 7.6 \times 10^{-6} \text{ cm}^2 \text{ V}^{-1} \text{ s}^{-1}$ and $\mu_y = 1.2 \times 10^{-5} \text{ cm}^2 \text{ V}^{-1} \text{ s}^{-1}$ averaging to $\mu_{xy} = 9.6 \times 10^{-6} \text{ cm}^2 \text{ V}^{-1} \text{ s}^{-1}$ are predicted. As we see below, this value is in very good agreement with the experimental results obtained from SCLC measurements, thus providing validation for our modeling approach.

The above calculations suggest a strong effect of static energetic disorder on the electron mobility in thin films of O-IDTBR. Such effects can be observed in experiments where charge carriers occupy states lying in the tail of the energetic distribution, typically referred to as trap states. Most previous measurements of electron mobility of thin-film O-IDTBR were obtained using the FET method, which is not sensitive to trap distributions as it tends to operate in a high charge carrier density regime. Consequently, O-IDTBR electron mobilities extracted from FET measurements tend to be relatively high (approximately $0.1 \text{ cm}^2 \text{ V}^{-1} \text{ s}^{-1}$) [36,53]. SCLC measurements, on the other hand, are usually carried out under low injection conditions and are very sensitive to the degree of disorder in the site energies, and therefore are more relevant to the low-charge-density simulations done with FOB-SH and KMC above. In a nonideal device, an SCLC experiment does not yield an unambiguous single value for mobility, but by doing measurements at different temperatures one can pin down some of the nonideal factors and in particular the density of states [54], i.e., degree of static energetic disorder. Having an estimate of the density of states then allows us to simulate the effective mobility under different levels of charge injection.

Such a method was applied to a disordered thin film of O-IDTBR. Samples were spin coated without adding

$1,8$-diiodooctane and without thermal postannealing (details in Sec. IV) to minimize crystallization and to better represent the disordered assembly used in simulations [37]. It was observed that the temperature-dependent SCLC data were best fit by a system with a Gaussian-type density-of-states distribution centered on the LUMO energy and with a Gaussian width of 0.150 eV, along with a trap-free mobility of $1.83 \times 10^{-3} \text{ cm}^2 \text{ V}^{-1} \text{ s}^{-1}$ [see Fig. 9(a) and Supplemental Material Tables S2 and S3 [44] for fit parameters]. Note that the Gaussian density-of-states distribution is of comparable width to the energetic disorder of 0.159 eV arising from the electrostatic calculations when all contributions are included (Table II). With a typical free-charge-carrier density for OPVs at room temperature at 1 sun, which is between $10^{22}$ and $10^{23} \text{ m}^{-3}$, we obtained an effective mobility lying between $1 \times 10^{-5}$ and $2 \times 10^{-5} \text{ cm}^2 \text{ V}^{-1} \text{ s}^{-1}$ at 0 V applied bias at room temperature, as shown in Fig. 9(b), which agrees well with results from KMC calculations when electrostatic disorder is included.

When comparing the SCLC mobilities with the KMC results, it is important to consider the different transport conditions in experiment and the KMC model. In particular, the higher charge density present in the device during SCLC experiments can have the effect of increasing mobilities due to trap-filling effects. For this reason, we provide the charge density at which the mobility obtained from SCLC is estimated. Finally, we note that in more ordered thin films as, e.g., obtained after thermal postannealing, the presence of crystalline domains will lead to an increase in charge mobility (see Ref. [34] and references therein for the well-studied system pentacene). The higher SCLC mobility values reported previously for O-IDTBR thin films could have arisen from a more crystalline thin film [55].

<figure><img src="image_11.png" /><figcaption>FIG. 9. Charge-transport mobility measured using temperature-dependent SCLC experiments. (a) Current-voltage characteristics at forward bias at different temperatures (220 to 300 K). Colored circles are experimental results, and black lines are fitting results using a free-charge-carrier drift-diffusion model via OghmaNano [56,57]. (b) Calculated effective mobility as a function of the applied bias assuming different free-charge-carrier densities (no injection, $n_f = 10^{22} \text{ m}^{-3}$, or $n_f = 10^{23} \text{ m}^{-3}$), where $n_f$ is the density for free-charge carriers. Illumination at 1 sun typically results in $n_f$ between $10^{22}$ and $10^{23} \text{ m}^{-3}$ in organic photovoltaics.</figcaption></figure>
### III. CONCLUSIONS

In this work, we carried out state-of-the-art molecular simulations and experiments to characterize electron transport in the nonfullerene acceptor material O-IDTBR. We found that both the electron-transport mechanism and the electron mobilities change drastically when going from single crystals to thin-film morphologies. In the crystal, the excess electron delocalizes modestly, over about three molecules on average, due to substantial electronic couplings in the a-b plane. More extended delocalization is hampered by off-diagonal and diagonal electron-phonon coupling. The FOB-SH nonadiabatic molecular-dynamics simulations predicted that the carrier propagates through the material via a transient delocalization mechanism resulting in relatively high electron mobilities that rival those of high-mobility acenes. In the thin-film model structures investigated, the excess electron localizes on a single molecule and moves across the material via site-to-site hopping. Experimentally validated mobility values for the thin film were predicted to be 5-6 orders of magnitude lower than those in the single crystal.

Experimental measurements of the electron mobility in a single crystal of O-IDTBR would help support this prediction. However, the relatively large size of the molecule and the impact of the solubilizing side chains make these molecules, like most solution-processible OSs, notoriously difficult to crystallize. Direct experimental evidence for the transient delocalization mechanism has not yet been reported in the literature. However, indirect experimental evidence for this mechanism exists for charge [18,28] and exciton transport in OSs [58,59]. Given that the FOB-SH method successfully predicted a number of important transport properties in OSs, e.g., the average delocalization length of electron holes in pentacene crystals [26], the different temperature dependences of mobility in two structurally similar OS crystals [28], as well as the transition from bandlike to small polaron hopping transport [60], we are confident that the predictions for crystalline O-IDTBR are robust, in particular, that carrier transport is in the regime of transient delocalization.

We systematically investigated the impact of the different types of disorder in the thin-film model that resulted in this dramatic loss in mobility. The packing and conformational disorder in the thin film and the resultant smaller electronic coupling values were sufficient to cause localization of the excess electron on a single molecule resulting in a decrease in electron mobility by about 1 order of magnitude. Yet, electrostatic disorder in the thin-film model had the most dramatic impact decreasing mobilities by an additional 4-5 orders of magnitude to values that are typically reported for structurally disordered organic semiconductors. The electrostatic disorder is mainly caused by the static contribution, that is, the different mean electrostatic potential at the different molecular sites. They in turn are generated by the different orientations of the molecules

with respect to one another and the different conformations of the molecules that form as a result of the soft dihedral potential between IDT and BT groups as well as the different partial charges on those groups that enhance intermolecular electrostatic interactions.

The results obtained here for O-IDTBR are likely to apply to several other NFAs featuring an A-D-A or A-DAD-A motif with soft dihedral potentials between A and D moieties. Our study suggests that while these motifs are beneficial for photoabsorption and exciton transport [42], the electrostatic disorder they create can strongly limit carrier transport in thin-film optoelectronic applications. Greater conformational rigidity between A and D groups and/or higher crystallinity are expected to reduce electrostatic disorder and increase charge mobility in thin films. This could be a reason why other NFAs, e.g., Y6, exhibit higher and often ambipolar mobilities. It will be of interest to investigate these aspects in future work alongside the question of how the structural and electrostatic disorder affect charge generation and recombination processes in optoelectronic donor-NFA interfaces.

### IV. METHODS

#### A. FOB-SH simulation of electron transport and mobility

##### 1. Theory

The FOB-SH approach is described in detail in a series of publications [25,42,60-62]. In this method, the valence-or conduction-band states of organic semiconductors are described by the following Hamiltonian:

$$ \hat{H} = \sum_{k} \epsilon_{k} |\phi_{k}\rangle\langle\phi_{k}| + \sum_{k \neq l} H_{kl} |\phi_{k}\rangle\langle\phi_{l}|, \quad (1) $$

where $\phi_k = \phi_k(\mathbf{r}, \mathbf{R}(t))$ are the quasidiabatic or charge localized wave functions approximated here for electron transport by the LUMOs of the isolated molecules, $\mathbf{r}$ describes the position of the excess electron, $\mathbf{R}(t)$ is the vector containing all nuclear positions, $\epsilon_k = \epsilon_k(\mathbf{R}(t)) = H_{kk}(\mathbf{R}(t))$ is the site energy of site $k$, i.e., the total potential energy of the system with site $k$ occupied by an excess electron and all other sites in their neutral state, and $H_{kl} = H_{kl}(\mathbf{R}(t))$ is the electronic coupling between $\phi_k$ and $\phi_l$. Note that the time dependence of diagonal and off-diagonal Hamiltonian matrix elements $\epsilon_k$ and $H_{kl}$, respectively, is due to the time dependence of the nuclear coordinates. The total excess electron wave function $[\Psi(t)]$ is expanded in the basis spanned by localized wave functions $[\phi_k(\mathbf{r}, \mathbf{R}(t))]$ used to represent the Hamiltonian

$$ \Psi(t) = \sum_{k} u_{k}(t) \phi_{k}(\mathbf{r}, \mathbf{R}(t)), \quad (2) $$

where $u_k(t)$ are the expansion coefficients. When this wave-function form is inserted into the time-dependent
Schrödinger equation, a system of time-dependent equations for the wavefunction coefficients is obtained

$$i\hbar\dot{u}_k(t) = \sum_{l=1}^{M} u_l(t)[H_{kl}(\mathbf{R}(t)) - i\hbar d_{kl}(\mathbf{R}(t))], \quad (3)$$

where $d_{kl}(\mathbf{R}(t)) = \langle \phi_k | \dot{\phi}_l \rangle$ are the time-dependent non-adiabatic coupling elements between localized states $k$ and $l$. The nuclei are propagated on a single adiabatic electronic state (termed “active state”), which is obtained by diagonalization of the Hamiltonian Eq. (1). The electronic and nuclear motions are coupled by the time dependence of all Hamiltonian matrix elements. Nonadiabatic transitions of the nuclear dynamics from the potential energy surface of the active state and the one of other adiabatic electronic eigenstates are determined stochastically according to Tully’s surface-hopping probability [63]. The FOB-SH method includes improvements over the original surface-hopping algorithm, such as decoherence correction, trivial crossing detection by a state-tracking algorithm based on tracking the overlaps between states in subsequent time steps, and eliminating spurious long-range charge-transfer states induced by the decoherence correction [62]. Recently, a similar surface-hopping methodology has been applied to exciton transport [64] in organic crystals and to exciton dissociation in a simple donor-acceptor system [65] and was found to agree well with multiconfigurational time-dependent Hartree simulations [65] considered to be the gold standard for nonadiabatic quantum dynamics [66].

The mean-square displacement MSD$_{αβ}$ of the charge carrier wave function $\Psi(t)$, Eq. (2), obtained by propagation via Eq. (3), is calculated as

$$\text{MSD}_{\alpha\beta} = \frac{1}{N_{\text{trj}}} \sum_{n=1}^{N_{\text{trj}}} \langle \Psi_n(t) | (\alpha - \alpha_0)(\beta - \beta_0) | \Psi_n(t) \rangle, \quad (4)$$

where $\alpha, \beta$ are the Cartesian coordinates of the excess electron in the high-mobility plane of the O-IDTBR crystal, $\alpha_0, \beta_0$ the initial positions of the center of charge of the excess electron wave function, and $N_{trj}$ the number of FOB-SH trajectories. Notice that the MSD Eq. (4) accounts for both spreading and translational or center of charge motion of the wave function [26]. After the initial (quantum) relaxation, $\Psi(t)$ enters a diffusive regime, where MSD components increase approximately linearly with time. In this regime, linear fits of the MSD components give the diffusion tensor via the Einstein relation,

$$D_{\alpha\beta} = \frac{1}{2} \lim_{t \to \infty} \frac{d\text{MSD}_{\alpha\beta}(t)}{dt}, \quad (5)$$

and the mobility tensor

$$\mu_{\alpha\beta} = \frac{eD_{\alpha\beta}}{k_B T}. \quad (6)$$

The diffusion constants and corresponding mobilities for charge carrier motion in the thin film are obtained from the MSD as

$$D_{\gamma} = \frac{1}{m} \lim_{t \to \infty} \frac{d\text{MSD}_{\gamma}(t)}{dt}, \quad (7)$$

where for $\gamma = z$ (out of plane) MSD$_z = (1/N_{trj}) \times \sum_{n=1}^{N_{trj}} \langle \Psi_n(t) | (z - z_0)^2 | \Psi_n(t) \rangle$ and $m = 2$, for $\gamma = xy$ (in plane) MSD$_{xy} = (1/N_{trj}) \sum_{n=1}^{N_{trj}} \langle \Psi_n(t) | (x - x_0)^2 + (y - y_0)^2 | \Psi_n(t) \rangle$ and $m = 4$, and for $\gamma = xyz$ (3D) MSD$_{xyz} = (1/N_{trj}) \sum_{n=1}^{N_{trj}} \langle \Psi_n(t) | (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 | \Psi_n(t) \rangle$ and $m = 6$.

The IPR of the charge carrier wave function $\langle \Psi(t) \rangle$ is a measure for the number of molecules over which the wave function is delocalized. It is computed as [26]

$$\langle \text{IPR} \rangle = \frac{1}{T_{\text{trj}}} \int_{0}^{T_{\text{trj}}} dt \frac{1}{N_{\text{trj}}} \sum_{n=1}^{N_{\text{trj}}} \frac{1}{\sum_{k=1}^{M} |u_{k,n}(t)|^4}, \quad (8)$$

where $u_{k,n}$ is the expansion coefficient $k$ of the wave function $\Psi(t)$ in trajectory $n$, $M$ is the number of electronically active molecules, $N_{trj}$ is the number of trajectories, and $T_{trj}$ is the trajectory length in time.

## 2. Parametrization of electronic Hamiltonian

For the calculation of the site energy $e_k$, the intramolecular and Lennard-Jones interaction terms for the neutral molecules are taken from the force field that was previously parametrized using B3LYP/6-311G(d,p) calculations [37,38]. For parametrization of the molecule in the anionic charge state, the equilibrium bond lengths are adjusted so that the intramolecular reorganization energy $\lambda_{intra}$ of O-IDTBR in vacuum obtained from the force field matches the corresponding value obtained from the B3LYP/6-311G(d,p) calculations,

$$\lambda_{\text{intra}} = [E_C(\mathbf{R}_N) + E_N(\mathbf{R}_C)] - [E_C(\mathbf{R}_C) + E_N(\mathbf{R}_N)], \quad (9)$$

where $E_C(\mathbf{R}_N)$ is the energy of the charged molecule $E_C$ at the minimum of neutral state $(\mathbf{R}_N)$, $\lambda_{intra} = 180$ meV; see Figs. S1 and S2 in the Supplemental Material [44]. We choose the B3LYP functional to remain consistent with the previous force-field parametrization and note that a range-separated functional CAM-B3LYP gives a similar reorganization energy of 213 meV. The difference (33 meV) does not have any significant effect on the ratio between thermally averaged electronic coupling along the high-mobility direction and reorganization energy that determines the extent of carrier delocalization, transport
mechanism, and mobility $\langle H_{ab}^2 \rangle^{1/2}/\lambda_{\text{intra}} = 0.53$ at B3LYP and 0.45 at the CAM-B3LYP level. In both cases, small polarons are no longer stable and charge carrier delocalization occurs. Force-field parameters for angular and dihedral terms are the same as for the neutral state, since differences are found to be negligibly small, and the same Lennard-Jones parameters are used in both charge states. As discussed in detail in Sec. II B, the site energies in the present FOB-SH simulations do not contain an electrostatic energy contribution. The off-diagonal electronic Hamiltonian matrix elements (electronic couplings $H_{kl}$) are calculated using the analytic overlap method (AOM) [67,68]. The orbital overlap $S_{kl}$ between the LUMO orbitals $\phi_k$ and $\phi_l$ are evaluated on the fly during dynamics and the couplings approximated as $H_{kl} = CS_{kl}$. The C constant is obtained from the $H_{kl}$ vs $S_{kl}$ correlation plot computed on a set of O-IDTBR dimer geometries taken from MD simulations of the crystal; Figs. S3 and S4 in the Supplemental Material [44]. The electronic couplings $H_{kl}$ are computed applying the projector-operator-based diabatization (POD) method [69,70] in combination with the PBE [71,72] functional and the double-zeta valence polarized basis set and the Goedecker-Teter-Hutter pseudopotentials for the core electrons. The POD couplings are scaled by a factor of 1.325 as recommended in Ref. [73] where POD was benchmarked against CASSCF/NEVPT2 electronic couplings for organic dimer anions [73,74]. For the crystal, it is necessary to introduce two C constants with the value depending on the dimer pose $C_1 = -0.4854$ Ha for the D1 and D4 dimers and $C_2 = -0.2548$ Ha for the D2 and D3 dimers; Fig. S3 in the Supplemental Material [44]. For the thin-film model, it is found that a single C constant equal to the value of $C_2$ is sufficient; Fig. S4 in the Supplemental Material [44]. Further details on the parametrization can be found in the Supplemental Material [44].

### 3. Simulation details

For the FOB-SH simulations of crystalline O-IDTBR, two supercells of different sizes were built starting from the experimental crystal structure 10 × 10 × 1 (400) and 14 × 7 × 1 (392), with the total number of molecules given in parentheses. Each supercell was first relaxed to their local minimum energy configuration and subsequently equilibrated to 300 K by running standard MD simulations for 0.5 ns in the NVT ensemble and applying periodic boundary conditions. The dynamics was run with a single molecule in the center of the cell in the negatively charged state and all other molecules in the neutral state. The same force field was used as described above for electronic Hamiltonian parametrization. The final configurations and velocities from these runs were used to initialize MD trajectories in the NVE ensemble, run for 0.5 ns. The structure of the O-IDTBR crystal was stable during these runs, root-mean-square-deviation = 0.7 Å with respect to the crystal structure. The initial conditions for the FOB-SH dynamics were

chosen from the equilibrated NVE structures (100 structures and velocities taken from the last 100 ps of NVE dynamics with 1 ps spacing). The FOB-SH simulations were performed for rectangular electronically active regions in the a-b crystallographic plane 5 × 5 × 1 (100) and 10 × 3 × 1 (120) with the number of electronically active molecules given in parentheses. Only the electronically active molecules were used for the construction of the electronic Hamiltonian, Eq. (1), and their LUMO orbitals were used for the expansion of the charge carrier wave function, Eq. (2). The remaining molecules (composed of one or two molecular layers embedding the active region) were electronically inactive but interacted with the electronically active region via their nonbonded force-field terms to ensure their structural integrity. The excess electron dynamics in the electronically active region was not periodically repeated; i.e., there were no interactions of the excess electrons with period images. In accordance with the prepared initial configuration, the charge carrier wave function [$\Psi(t)$] was chosen to be an eigenstate of the electronic Hamiltonian, localized on a central molecule in the active regions, i.e., a small polaron. The wave function and the nuclei were propagated in time by applying the FOB-SH algorithm in the NVE ensemble, applying the decoherence correction, removal of the spurious long-range charge-transfer transitions, trivial crossing detection, adjustment of the velocities in the direction of the nonadiabatic coupling vector, and applying the multiple-time-step algorithm. The nuclear time step was set to 0.025 fs. The convergence of the electron mobility with respect to the nuclear time step was investigated and summarized in Table S1 in the Supplemental Material [44]. The electronic time step was equal to 1/5 of the nuclear time step. The number of simulated trajectories varied between 250 and 500, each of length 1 ps. The number and length of trajectories were determined to be large enough to converge the MSD tensor. The FOB-SH simulations were found to recover important properties, specifically, the internal consistency between the electronic wave function and surface populations and a detailed balance of the electronic populations in the long-time limit; see Fig. S5 in the Supplemental Material [44].

An amorphous structural assembly of O-IDTBR was previously generated by nonequilibrium solvent evaporation simulation on an approximately 2 nm layer of amorphous silicon-dioxide substrate [37,38]. We adopted this structure as a model representative of a disordered thin film and used it as an initial structure for MD simulation in the NVT ensemble. The simulations were carried out for 0.5 ns applying periodic boundary conditions in two directions parallel to the substrate (x and y). In order to maintain the structure along the z direction where a periodic boundary condition (PBC) was not applied, we froze the coordinates of the molecules containing atoms within 3 Å thick layers on the surface of the sample along the z direction; i.e., all molecules with the atoms within this
region (214 molecules) were kept frozen in MD and treated electronically inactive in the FOB-SH simulations. The rest of the molecules (234) could move freely in MD simulations and were treated electronically active in the FOB-SH simulations. The equilibrated configuration was used as an initial structure for a 0.5 ns NVE dynamics. The last 100 configurations (1 ps spacing) were used to build electronic coupling maps and to determine the distribution of the length of the electronic coupling paths. For each of the 100 configurations, we simulated FOB-SH dynamics starting from three different initial wave functions of the excess electron in each of the 100 configurations, ones where it is localized on a molecule that belongs to a short, medium, and long electronic coupling path. The FOB-SH dynamics was performed with the same setup as in the case of crystalline systems. All MD and FOB-SH simulations were done with an in-house implementation of FOB-SH in the CP2K program package [75].

## B. KMC calculation of electron transport and mobility

### 1. Theory

In the KMC simulations, the excess electron is described by a small polaron that hops from one molecule to the next under the influence of an applied external electric field F. Hopping rates $\Gamma_{kl}$ between different molecules are defined by high-temperature nonadiabatic electron transfer (or Marcus) theory

$$ \Gamma_{kl} = \frac{\langle H_{kl}^2 \rangle}{\hbar} \sqrt{\frac{\pi}{\lambda k_B T}} \exp \left( -\frac{(\Delta E_{kl} + \lambda)^2}{4\lambda k_B T} \right), \quad (10) $$

where

$$ \Delta E_{kl} = \langle \epsilon_l \rangle - \langle \epsilon_k \rangle - e \mathbf{F} \cdot \mathbf{r}_{kl}, \quad (11) $$

$\langle \epsilon_k \rangle$ is the thermal average of the site energy $\epsilon_k$, $\lambda$ is the reorganization energy, $\langle H_{kl}^2 \rangle$ and $\mathbf{r}_{kl}$ are the thermally averaged electronic coupling and separation distance, respectively, between molecular sites k and l, and e is the elementary charge. Thus, different forward and backward hopping rates are present for sites with different positions in the field direction. During the KMC simulations, the time evolution of a single excess electron is propagated until a set simulation time is reached, at which point the simulation halts, and the average velocity of the electron in the field direction $\langle v \rangle$ is computed. The electron mobility along the field direction is then calculated from

$$ \mu = \frac{\langle v \rangle}{F}, \quad (12) $$

where $F = |\mathbf{F}|$ is the field strength.

### 2. Simulation details

The KMC simulations were carried out for a single excess electron on the same thin-film structures that were used as initial coordinates in the FOB-SH simulations. Each molecule was represented by a site with a single coordinate. Site positions were determined from the center of mass coordinates for each molecule, averaged over the 100 configurations previously described as the starting point for the FOB-SH simulations. Only the 234 active molecules were included. The site energies $\epsilon_k$ were evaluated at four levels of theory: without electrostatic interactions, with the electrostatic contribution arising from the external potential created by atomic charges of neutral surrounding molecules ($\Delta_E^{-}$), with the total electrostatic contribution ($\Delta_E^{-} + \Delta_{In}^{-}$) including also electrostatic potentials generated by dipoles and multipoles induced by the charge density of the neutral molecules ($\Delta_{In}^{-}$), and with the total electrostatic and induction contributions ($\Delta_E^{-} + \Delta_{In}^{-} + \Delta_{le}^{-}$); see Sec. IV C below for details on the calculation of these contributions and the thermal averaging. For each set of site energies considered, a corresponding reorganization energy $\lambda$ was determined, which was used for all hops in that simulation. The total $\lambda$ was taken as the sum of the previously described $\lambda_{intra}$, Eq. (9), and an additional intermolecular contribution $\lambda_{inter}$ determined based on the energetic landscape considered for that simulation run. For the case where we considered only $\epsilon_k$ without electrostatic interactions, thermal averaging led to site energies of zero, and only intramolecular reorganization energy was taken into account in charge hopping rate computations. The electronic couplings $H_{kl}$ between all dimer pairs in the thin film used in the KMC computations were averaged over several FOB-SH trajectories. In addition to those calculated in FOB-SH, $H_{kl}$ was also computed between sites separated across periodic boundaries in the x and y directions parallel to the substrate. This allowed the excess electron to cross between boundaries and loop over the system multiple times during a simulation run. To this end, the thin-film structures during 20 FOB-SH trajectories (100 structures from each of them) were replicated by a single periodic image in the x and y directions. In each of the sampled geometries, for each of the molecules in the central cell, the molecules from the eight surrounding cells with the centroid distance smaller than 40 Å were isolated. The AOM couplings of the created dimers were computed and averaged over all structures, yielding $\langle H_{kl}^2 \rangle$. The thermal averaging over the electron-transfer parameters requires some further comments. The timescales of intramolecular site energy fluctuations are on the typical time-scale of intramolecular vibrations ($10^{-14}-10^{-13}$ s), and they are sufficiently fast compared to the computed electron hopping rates. The time constant distributions of the latter peak between $10^{-12}$ and $10^{-11}$ s; see Fig. S8 in the Supplemental Material [44]. For this reason, thermally averaged values for intramolecular site energies were used.
The timescale of intermolecular (including electrostatic) site energy fluctuations and electronic coupling fluctuations is slower, typically on the timescale of intermolecular vibrations (10⁻¹³–10⁻¹² s). This timescale partially overlaps with the electron hopping rates. In keeping with standard practice, we opted to use also in this case thermally averaged values for the calculation of the Marcus rate, although this is less well justified than for the intramolecular site energies. The KMC runs were repeated applying an electric field along the x or y directions of the thin film until the corresponding electron mobilities in these directions, μₓ and μᵧ, converged. Large enough simulation times were considered so that the initial location of the mobile charge did not influence the mobility found for each run. The average KMC mobility in the x-y plane of the thin film was obtained by averaging μₓᵧ = (μₓ + μᵧ)/2. All simulations were carried out using the in-house KMC code ToFeT [43].

### C. Classical microelectrostatics model

#### 1. Theory

The interaction between the charge carrier and the molecular environment is responsible for the variations of site polarization energies Δₖ of molecule k in the sample. The disorder arising from different polarization energies experienced by different molecules is particularly important for O-IDTBR thin-film structures, where each molecule feels a different surrounding environment. The polarization (or environmental) energy is evaluated here with a classical ME model developed by D'Avino et al. [47,48] based on permanent atomic charges and distributed molecular polarizabilities that give rise to induced atomic dipoles. Self-consistent ME calculations are performed in order to account for the mutual interactions between induced dipoles. In a nutshell, ME allows for the calculation of the total energies of the neutral (indicated with 0) and, in this work, negatively charged (indicated with −) systems as

$$U^0 = \frac{1}{2} \sum_{ki} q_{ki}^0 (V_{ki}^0 + \tilde{V}_{ki}^0), \quad (13)$$

$$U^{-} = \frac{1}{2} \sum_{ki} (q_{ki}^{0} + q_{ki}^{\delta})(V_{ki}^{0} + V_{ki}^{\delta} + \tilde{V}_{ki}^{0} + \tilde{V}_{ki}^{\delta}), \quad (14)$$

where q₀ₖᵢ(q₀ₖᵢ + qᵢ) is the permanent charge on atom i of molecule k in the neutral (charged) state, V₀ₖᵢ(Vᵢ) is the electrostatic potential felt on atom i of molecule k from all other molecules in the neutral (charged) system. The differential charges (qᵢ) and potentials (Vᵢ) between the charged and neutral systems are denoted by a δ superscript. Exact partitioning of the site energy of the negatively charged molecule k gives [76]

$$\Delta^{-} = \Delta_{E}^{-} + \Delta_{\text{In}}^{-} + \Delta_{\text{le}}^{-}, \quad (15)$$

where the first two terms can be regarded as electrostatic contributions, while the third is the induction contribution to the site energy. The electrostatic term is a sum of a “bare” electrostatic contribution due to the electrostatic potential generated by permanent atomic charges of the neutral molecules (V₀ₖᵢ),

$$\Delta_{E}^{-} = \frac{1}{2} \sum_{k \in i} q_{ki}^{\delta} (V_{ki}^{0}), \quad (16)$$

and a contribution due to electrostatic potentials (Ṽ₀ₖᵢ) generated by dipoles and multipoles induced by the charge density of the neutral molecules:

$$\Delta_{\text{In}}^{-} = \frac{1}{2} \sum_{k \in i} q_{ki}^{\delta} (\tilde{V}_{ki}^{0}). \quad (17)$$

Therefore, this term accounts for the interaction between the molecular charge difference of the molecule of interest and the screened field of the surrounding neutral molecules. On the other hand, the induction term Δ⁻ᵢₑ is given by electric potentials (Ṽᵢ) generated by dipoles induced by the charge density of the negatively charged molecule itself:

$$\Delta_{\text{le}}^{-} = \frac{1}{2} \sum_{k \in i} q_{ki}^{\delta} \tilde{V}_{ki}^{\delta}. \quad (18)$$

#### 2. Calculation of electrostatic disorder and intermolecular reorganization energy

The ME model was parametrized by computing gas-phase electrostatic potential (ESP) atomic charges at the DFT ωB97X-D/6-311++G(d, p) level for the neutral and the charged (anionic) state. ESP charges were computed for all the individual active molecules of the thin-film sample (234 in total) extracted from the last MD frame of a 5-ns-long NVT trajectory. The polarizability tensor of the neutral species computed at the same level of theory as above was applied to describe both the neutral and the charged state. In order to evaluate the electrostatic interactions in Eq. (15), self-consistent calculations were performed on the thin-film sample, where (PBCs were applied in the x-y plane with a large cutoff of 800 Å to ensure convergence of the electrostatic energies. The induction term in Eq. (18) was evaluated separately on spherical clusters where each molecule in the thin film was placed at the center of the cluster. As the induction term scales as 1/R, where R is the cluster radius, extrapolation in the infinite-crystal limit was carried out. The polarization energy and its contributions Xₖᵢ, X = Δ⁻, Δ⁻ᵢ, Δ⁻ᵢₙ, Δ⁻ᵢₑ were calculated for all active sites k (total number of active sites K) and for L = 10 different MD configurations l, taken from the last 1 ns of the NVT trajectory for the thin film in equidistant steps of 100 ps.
The corresponding thermal averages were denoted $\langle X_{kl} \rangle$. The total electrostatic disorder and the corresponding static, dynamic, and cross contributions discussed in Sec. II C are defined by

$$ \sigma_{\text{tot}} = \left( \frac{1}{K} \sum_{k=1}^{K} \frac{1}{L} \sum_{l=1}^{L} (X_{kl} - \langle X_{kl} \rangle)^2 \right)^{1/2}, \quad (19) $$

$$ \sigma_{\text{stat}} = \left( \frac{1}{K} \sum_{k=1}^{K} (\langle X_{kl} \rangle_l - \langle X_{kl} \rangle)^2 \right)^{1/2}, \quad (20) $$

$$ \sigma_{\text{dyn}} = \left( \frac{1}{K} \sum_{k=1}^{K} \frac{1}{L} \sum_{l=1}^{L} (X_{kl} - \langle X_{kl} \rangle_l)^2 \right)^{1/2}, \quad (21) $$

$$ \sigma_{\text{cross}}^2 = \sigma_{\text{tot}}^2 - \sigma_{\text{stat}}^2 - \sigma_{\text{dyn}}^2, \quad (22) $$

where $\langle X_{kl} \rangle$ denotes averaging over all MD configurations $l$ and molecules $k$ and $\langle X_{kl} \rangle_l$ denote averaging over MD configurations $l$ of molecule $k$ only.

The intermolecular reorganization energy $\lambda_{\text{inter}}$ was calculated as explained in Sec. II C. The potential energy of the negatively charged state in the optimum geometry of the negatively charged state denoted $X_-$ in Table II was calculated by performing ME calculations on a crystalline sample, where a molecule in the center was optimized at the DFT level [ωB97X-D/6-31G(d, p)] in the negatively charged state by freezing angles and dihedrals of the alkyl chains. Atomic ESP charges were calculated for the neutral state of both the optimized molecule mentioned above and an unoptimized molecule extracted directly from the crystal. ME calculations to evaluate the electrostatic terms Eqs. (16) and (17) were performed with 3D PBC using a cutoff of 500 Å, while the induction term in Eq. (18) was evaluated for spherical clusters of increasing radius, with the molecule optimized in the anionic structure placed at its center. We assigned the charge difference between the negatively charged and the neutral state computed at the optimized geometry to this molecule. The energy of the system of the negatively charged state at the neutral geometry $X_0$ in Table II was calculated considering all the molecules in the crystal at their neutral geometry and repeating the steps described above for $X_-$. Here, the molecule in the center of the spherical cluster was assigned the charge difference between the negatively charged and the neutral state computed at the neutral geometry.

# D. Experimental mobility measurements

## 1. Device fabrication

The device structure for the O-IDTBR electron-transport device was ITO/ZnO(40 nm)/O-IDTBR(120 nm)/ BCP(10 nm)/Ag(80 nm). ITO substrates were precleaned using acetone, soap water, deionized water, and isopropanol in sequence (10 min each step). A 40 nm ZnO thin film was

then deposited on top of ITO from a zinc acetate anhydrous solution (110 mg/ml in 2-methoxyethanol with 30 µL ethanolamine per 1 mL) by spin coating at 4000 rpm, followed by 30 mins annealing at 200 deg. After that, pristine O-IDTBR film was spin coated from chlorobenzene solution (40 mg ml⁻¹) with a speed of 2000 rpm, yielding a thickness of approximately 120 nm. Bathocuproine was then evaporated using an organic source in a vacuum chamber with a thickness of 10 nm, followed by the evaporation of silver contact metal (80 nm).

## 2. Current-voltage measurement

The devices were measured with a voltage range of -5 to +5 V using Keithley 236 under dark. The temperature of the device was controlled by the Linkam stage (version HFS600E-PB4 PROBE STAGE) using liquid nitrogen. Each measurement was carried out after the temperature was stabilized for > 2 min.

## 3. Determination of electron mobilities

The electron mobility was estimated experimentally using the SCLC method at different temperatures assisted by a full free-charge drift-diffusion model fitting using OghmaNano [56,57]. The temperature-dependent SCLC experiment allowed us to better resolve the DOS of trap states (below gap states, or called tail states), and a Gaussian-type DOS for the conduction-band tail [$g(E)$] was assumed in the model fitting following

$$ g(E) = U^{\text{Gau}} \exp \left[ -\frac{1}{2} \left( \frac{(E - E_c)}{\sigma} \right)^2 \right], \quad (23) $$

where $U^{\text{Gau}}$ is the effective trap density, $E_c$ is the band edge of the conduction band, and $\sigma$ is the Gaussian width, i.e., static disorder of the conduction-band tail. The total trap density ($N_t$) is then determined via $N_t = U^{\text{Gau}}\sigma\sqrt{2\pi}$. In the fitting process, we aimed to fit the current-voltage (J-V) curves at forward bias at all temperatures in dark at once using one set of physical parameters. The inputs for the drift-diffusion model fitting are shown in Table S3 of the Supplemental Material [44]. We considered five main fitting parameters, within which the trap-free mobility ($\mu_0$), Gaussian width ($\sigma$), and effective trap density ($U^{\text{Gau}}$) were the most important. The trap-free mobility and DOS of the trap states extracted based on the Gaussian model allowed us to estimate the effective electron mobility at room temperature at a given charge density (injection level) for the free electrons, following the multiple trapping and detrapping framework via [57,77]

$$ \mu_{\text{eff}} = \mu_0 \frac{n_f}{n_f + n_t}, \quad (24) $$

where $\mu_0$ is the trap-free mobility, $n_f$ is the free-electron density above the conduction band, and $n_t$ is the trapped charge density in the Gaussian tail, i.e., trap states.
ACKNOWLEDGMENTS

L. S. was supported by the European Research Council (ERC) under the European Union, Horizon 2020 research and innovation program (Grant Agreement No. 682539/ SOFTCHARGE). J. N., J. Y., and J. F. C. thank the ERC for financial support under the European Union, Horizon 2020 research and innovation program (Grant Agreement No. 742708/CAPaCITy). J. N. thanks the Royal Society for the research professorship. J. Y. thanks Guangdong Province for support under Guangdong Provincial Natural Science Foundation General Project (Grant No. 2024A1515012318) and Youth Fund (Grant No. 2023A1515111140), the Shenzhen government for support under the special appointed position “Pengcheng Peacock Plan—C,” and The Chinese University of Hong Kong (Shenzhen) for support under The University Development Fund (Grant No. UDF01003117) and special funds for high-level universities-talent projects, Presidential Young Fellow (Grant No. UF02003117). S. G., G. L., and D. B. acknowledge funding from the Belgian National Fund for Scientific Research within the Consortium des Equipements de Calcul Intensif, under Grant No. 2.5020.11, and by the Walloon Region (ZENOBE Tier-1 supercomputer) under Grant No. 1117545. Via our membership of the UK’s HEC Materials Chemistry Consortium, which is funded by EPSRC (Grants No. EP/L000202 and No. EP/R029431), this work used the ARCHER UK National Supercomputing Service. We are grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC (Grant No. EP/P020194/1).

[1] J. Hou, O. Inganas, R. H. Friend, and F. Gao, Organic solar cells based on non-fullerene acceptors, Nat. Mater. **17**, 119 (2018).

[2] A. Karki, A. J. Gillett, R. H. Friend, and T.-Q. Nguyen, The path to 20% power conversion efficiencies in nonfullerene acceptor organic solar cells, Adv. Energy Mater. **11**, 2003441 (2020).

[3] K. Vandewal, S. Mertens, J. Benduhn, and Q. Liu, The cost of converting excitons into free charge carriers in organic solar cells, J. Phys. Chem. Lett. **11**, 129 (2020).

[4] R. Wang, C. Zhang, Q. Li, Z. Zhang, X. Wang, and M. Xiao, Charge separation from an intra-moiety intermediate state in the high-performance PM6:Y6 organic photovoltaic blend, J. Am. Chem. Soc. **142**, 12751 (2020).

[5] R. Jasiunas, H. Zhang, J. Yuan, X. Zhou, D. Qian, Y. Zou, A. Devizis, J. Sulskus, F. Gao, and V. Gulbinas, From generation to extraction: A time-resolved investigation of photophysical processes in non-fullerene organic solar cells, J. Phys. Chem. C **124**, 21283 (2020).

[6] C. Li, J. Zhou, J. Song, J. Xu, H. Zhang, X. Zhang, J. Guo, L. Zhu, D. Wei, G. Han, J. Min, Y. Zhang, Z. Xie, Y. Yi, H.

Yan, F. Gao, F. Liu, and Y. Sun, Non-fullerene acceptors with branched side chains and improved molecular packing to exceed 18% efficiency in organic solar cells, Nat. Energy **6**, 605 (2021).

[7] Y. Dong, H. Cha, H. L. Bristow, J. Lee, A. Kumar, P. S. Tuladhar, I. McCulloch, A. A. Bakulin, and J. R. Durrant, Correlating charge-transfer state lifetimes with material energetics in polymer:non-fullerene acceptor organic solar cells, J. Am. Chem. Soc. **143**, 7599 (2021).

[8] S. Karuthedath et al., Intrinsic efficiency limits in low-bandgap non-fullerene acceptor organic solar cells, Nat. Mater. **20**, 378 (2021).

[9] W. Deng, W. Liu, R. Qian, and H. Wu, Toward high-efficiency organic photovoltaics: Perspectives on the origin and role of energetic disorder, J. Phys. Chem. Lett. **13**, 544 (2022).

[10] Z. Chen and H. Zhu, Photoinduced charge transfer and recombination dynamics in star nonfullerene organic solar cells, J. Phys. Chem. Lett. **13**, 1123 (2022).

[11] C. Li, X. Gu, Z. Chen, X. Han, N. Yu, Y. Wei, J. Gao, H. Chen, M. Zhang, A. Wang, J. Zhang, Z. Wei, Q. Peng, Z. Tang, X. Hao, X. Zhang, and H. Huang, Achieving record-efficiency organic solar cells upon tuning the conformation of solid additives, J. Am. Chem. Soc. **144**, 14731 (2022).

[12] E. Saglamkaya, A. Musiienko, M. S. Shadabroo, B. Sun, S. Chandrabose, O. Shargaieva, G. Lo Gerfo, N. F. van Hulst, and S. Shoaee, What is special about Y6; the working mechanism of neat Y6 organic solar cells, Mater. Horiz. **10**, 1825 (2023).

[13] L. Zhu et al., Single-junction organic solar cells with over 19% efficiency enabled by a refined double-fibril network morphology, Nat. Mater. **21**, 656 (2022).

[14] Y. Firdaus, V. M. Le Corre, S. Karuthedath, W. Liu, A. Markina, W. Huang, S. Chattopadhyay, M. M. Nahid, M. I. Nugraha, Y. Lin, A. Seitkhan, A. Basu, W. Zhang, I. McCulloch, H. Ade, J. Labram, F. Laquai, D. Andrienko, L. Jan Anton Koster, and T. D. Anthopoulos, Long-range exciton diffusion in molecular non-fullerene acceptors, Nat. Commun. **11**, 5220 (2020).

[15] A. Troisi and G. Orlandi, Charge-transport regime of crystalline organic semiconductors: Diffusion limited by thermal off-diagonal electronic disorder, Phys. Rev. Lett. **96**, 086601 (2006).

[16] S. Fratini, D. Mayou, and S. Ciuchi, The transient localization scenario for charge transport in crystalline organic materials, Adv. Funct. Mater. **26**, 2292 (2016).

[17] S. Fratini, S. Ciuchi, D. Mayou, G. Trambly de Laissardiere, and A. Troisi, A map of high-mobility molecular semiconductors, Nat. Mater. **16**, 998 (2017).

[18] S. Fratini, M. Nikolka, A. Salleo, G. Schweicher, and H. Sirringhaus, Charge transport in high-mobility conjugated polymers and molecular semiconductors, Nat. Mater. **19**, 491 (2020).

[19] T. Kubař and M. Elstner, A hybrid approach to simulation of electron transfer in complex molecular systems, J. R. Soc. Interface **10**, 20130415 (2013).

[20] A. Heck, J. J. Kranz, T. Kubař, and M. Elstner, Multi-scale approach to non-adiabatic charge transport in high-mobility
organic semiconductors, J. Chem. Theory Comput. **11**, 5068 (2015).

[21] S. Roosta, F. Ghalami, M. Elstner, and W. Xie, *Efficient surface hopping approach for modeling charge transport in organic semiconductors*, J. Chem. Theory Comput. **18**, 1264 (2022).

[22] L. Wang and D. Beljonne, *Flexible surface hopping approach to model the crossover from hopping to band-like transport in organic crystals*, J. Phys. Chem. Lett. **4**, 1888 (2013).

[23] L. Wang, O. V. Prezhdo, and D. Beljonne, *Mixed quantum-classical dynamics for charge transport in organics*, Phys. Chem. Chem. Phys. **17**, 12395 (2015).

[24] Z. Wang, J. Dong, and L. Wang, *Large-scale surface hopping simulation of charge transport in hexagonal molecular crystals: role of electronic coupling signs*, J. Phys. Condens. Matter **35**, 345401 (2023).

[25] J. Spencer, F. Gajdos, and J. Blumberger, *FOB-SH: Fragment orbital-based surface hopping for charge carrier transport in organic and biological molecules and materials*, J. Chem. Phys. **145**, 064102 (2016).

[26] S. Giannini, A. Carof, M. Ellis, H. Yang, O. G. Ziogos, S. Ghosh, and J. Blumberger, *Quantum localization and delocalization of charge carriers in organic semiconducting crystals*, Nat. Commun. **10**, 3843 (2019).

[27] S. Giannini and J. Blumberger, *Charge transport in organic semiconductors: The perspective from nonadiabatic molecular dynamics*, Acc. Chem. Res. **55**, 819 (2022).

[28] S. Giannini, L. Di Virgilio, M. Bardini, J. Hausch, J. Geuchies, W. Zheng, M. Volpi, J. Elsner, K. Broch, Y. H. Geerts, F. Schreiber, G. Schweicher, H. Wang, J. Blumberger, M. Bonn, and D. Beljonne, *Transiently delocalized states enhance hole mobility in organic molecular semiconductors*, Nat. Mater. **22**, 1361 (2023).

[29] N. B. Taylor and I. Kassal, *Generalised Marcus theory for multi-molecular delocalised charge transfer*, Chem. Sci. **9**, 2942 (2018).

[30] D. Balzer, T. J. A. M. Smolders, D. Blyth, S. N. Hood, and I. Kassal, *Delocalised kinetic Monte Carlo for simulating delocalisation-enhanced charge and exciton transport in disordered materials*, Chem. Sci. **12**, 2276 (2021).

[31] J. T. Willson, W. Liu, D. Balzer, and I. Kassal, *Jumping kinetic Monte Carlo: Fast and accurate simulations of partially delocalized charge transport in organic semiconductors*, J. Phys. Chem. Lett. **14**, 3757 (2023).

[32] T. F. Harrelson, V. Dantanarayana, X. Xie, C. Koshnick, D. Nai, R. Fair, S. A. Nuñez, A. K. Thomas, T. L. Murrey, M. A. Hickner, J. K. Grey, J. E. Anthony, E. D. Gomez, A. Troisi, R. Faller, and A. J. Moulé, *Direct probe of the nuclear modes limiting charge mobility in molecular semiconductors*, Mater. Horiz. **6**, 182 (2019).

[33] O. G. Ziogos, S. Giannini, M. Ellis, and J. Blumberger, *Identifying high-mobility tetracene derivatives using a non-adiabatic molecular dynamics approach*, J. Mater. Chem. C **8**, 1054 (2020).

[34] M. Ellis, H. Yang, S. Giannini, O. G. Ziogos, and J. Blumberger, *Impact of nanoscale morphology on charge carrier delocalization and mobility in an organic semiconductor*, Adv. Mater. **33**, 2104852 (2021).

[35] S. Holliday, R. S. Ashraf, A. Wadsworth, D. Baran, S. A. Yousaf, C. B. Nielsen, C.-H. Tan, S. D. Dimitrov, Z. Shang, N. Gasparini, M. Alamoudi, F. Laquai, C. J. Brabec, A. Salleo, J. R. Durrant, and I. McCulloch, *High-efficiency and air-stable P3HT-based polymer solar cells with a new non-fullerene acceptor*, Nat. Commun. **7**, 11585 (2016).

[36] H. Bristow, K. J. Thorley, A. J. P. White, A. Wadsworth, M. Babics, Z. Hamid, W. Zhang, A. F. Paterson, J. Kosco, J. Panidi, T. D. Anthopoulos, and I. McCulloch, *Impact of nonfullerene acceptor side chain variation on transistor mobility*, Adv. Electron. Mater. **5**, 1900344 (2019).

[37] A. S. Gertsen, M. K. Sørensen, and J. W. Andreasen, *Nano-structure of organic semiconductor thin films: Molecular dynamics modeling with solvent evaporation*, Phys. Rev. Mater. **4**, 075405 (2020).

[38] A. S. Gertsen, *Molecular dynamics force fields for organic solar cell materials*, Ph.D. thesis, Technical University of Denmark, 2021.

[39] J. Wu, J. Lee, Y. C. Chin, H. Yao, H. Cha, J. Luke, J. Hou, J.-S. Kim, and J. R. Durrant, *Exceptionally low charge trapping enables highly efficient organic bulk heterojunction solar cells*, Energy Environ. Sci. **13**, 2422 (2020).

[40] A. Wadsworth, M. Moser, A. Marks, M. S. Little, N. Gasparini, C. J. Brabec, D. Baran, and I. McCulloch, *Critical review of the molecular design progress in non-fullerene electron acceptors towards commercially viable organic solar cells*, Chem. Soc. Rev. **48**, 1596 (2019).

[41] R.-Z. Liang, M. Babics, V. Savikhin, W. Zhang, V. M. Le Corre, S. Lopatin, Z. Kan, Y. Firdaus, S. Liu, I. McCulloch, M. F. Toney, and P. M. Beaujuge, *Carrier transport and recombination in efficient “all-small-molecule” solar cells with the nonfullerene acceptor IDTBR*, Adv. Energy Mater. **8**, 1800264 (2018).

[42] S. Giannini, A. Carof, M. Ellis, O. G. Ziogos, and J. Blumberger, *From atomic orbitals to nano-scale charge transport with mixed quantum/classical non-adiabatic dynamics: Method, implementation and application*, in *Multi-scale Dynamics Simulations: Nano and Nano-Bio Systems in Complex Environments* (The Royal Society of Chemistry, Cambridge, 2022), Chap. 6, pp. 172–202.

[43] J. J. Kwiatkowski, J. M. Frost, and J. Nelson, *The effect of morphology on electron field-effect mobility in disordered C60 thin films*, Nano Lett. **9**, 1085 (2009).

[44] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.14.021021 for further details on parametrization of electronic Hamiltonian, convergence of FOB-SH electron mobilities and IPR, detailed balance internal consistency and energy conservation, illustration of electron transport mechanism, KMC hopping rates, and parameters for drift diffusion model fitting.

[45] S. Giannini, O. G. Ziogos, A. Carof, M. Ellis, and J. Blumberger, *Flickering polarons extending over ten nanometres mediate charge transport in high-mobility organic crystals*, Adv. Theory Simul. **3**, 2000093 (2020).

[46] A. F. Paterson, R. Li, A. Markina, L. Tsetseris, S. MacPhee, H. Faber, A.-H. Emwas, J. Panidi, H. Bristow, A. Wadsworth, D. Baran, D. Andrienko, M. Heeney, I. McCulloch, and T. D. Anthopoulos, *N-doping improves charge transport and morphology in the organic semiconductors*, J. Phys. Chem. C **114**, 1054 (2009).
non-fullerene acceptor O-IDTBR, J. Mater. Chem. C **9**, 4486 (2021).

[47] G. D'Avino, L. Muccioli, C. Zannoni, D. Beljonne, and Z. G. Soos, *Electronic polarization in organic crystals: A comparative study of induced dipoles and intramolecular charge redistribution schemes*, J. Chem. Theory Comput. **10**, 4959 (2014).

[48] G. D'Avino, L. Muccioli, F. Castet, C. Poelking, D. Andrienko, Z. G. Soos, J. Cornil, and D. Beljonne, *Electronic phenomena in organic semiconductors: Fundamentals and implications for photovoltaics*, J. Phys. Condens. Matter **28**, 433002 (2016).

[49] G. King and A. Warshel, *Investigation of the free energy functions for electron transfer reactions*, J. Chem. Phys. **93**, 8682 (1990).

[50] S. Gupta and D. V. Matyushov, *Effects of solvent and solute polarizability on the reorganization energy of electron transfer*, J. Phys. Chem. A **108**, 2087 (2004).

[51] J. Blumberger and G. Lamoureux, *Reorganization free energies and quantum corrections for a model electron self-exchange reaction: Comparison of polarizable and non-polarizable solvent models*, Mol. Phys. **106**, 1597 (2008).

[52] X. Jiang, Z. Futera, and J. Blumberger, *Ergodicity-breaking in thermal biological electron transfer? Cytochrome c revisited*, J. Phys. Chem. B **123**, 7588 (2019).

[53] M. I. Nugraha, M. Gedda, Y. Firdaus, A. D. Scaccabarozzi, W. Zhang, S. Alshammari, F. Anies, B. Adilbekova, A.-H. Emwas, I. McCulloch, M. Heeney, L. Tsetseris, and T. D. Anthopoulos, *Addition of diquat enhances the electron mobility in various non-fullerene acceptor molecules*, Adv. Funct. Mater. **32**, 2202954 (2022).

[54] J. A. Röhr, X. Shi, S. A. Haque, T. Kirchartz, and J. Nelson, *Charge transport in spiro-OMeTAD investigated through space-charge-limited current measurements*, Phys. Rev. Appl. **9**, 044017 (2018).

[56] R. C. MacKenzie, https://www.oghma-nano.com/ (2022).

[57] B. Xiao, P. Calado, R. C. I. MacKenzie, T. Kirchartz, J. Yan, and J. Nelson, *Relationship between fill factor and light intensity in solar cells based on organic disordered semiconductors: The role of tail states*, Phys. Rev. Appl. **14**, 024034 (2020).

[55] D. Baran, R. S. Ashraf, D. A. Hanifi, M. Abdelsamie, N. Gasparini, J. A. Röhr, S. Holliday, A. Wadsworth, S. Lockett, M. Neophytou, C. J. M. Emmott, J. Nelson, C. J. Brabec, A. Amassian, A. Salleo, T. Kirchartz, J. R. Durrant, and I. McCulloch, *Reducing the efficiency-stability-cost gap of organic photovoltaics with highly efficient and stable small molecule acceptor ternary solar cells*, Nat. Mater. **16**, 363 (2017).

[58] A. J. Sneyd, T. Fukui, D. Paleček, S. Prodhan, I. Wagner, Y. Zhang, J. Sung, S. M. Collins, T. J. A. Slater, Z. Andaj-Garmaroudi, L. R. MacFarlane, J. Diego Garcia-Hernandez, L. Wang, G. R. Whittell, J. M. Hodgkiss, K. Chen, D. Beljonne, I. Manners, R. H. Friend, and A. Rao, *Efficient energy transport in an organic semiconductor mediated by transient exciton delocalization*, Sci. Adv. **7**, eabh4232 (2021).

[59] A. J. Sneyd, D. Beljonne, and A. A. Rao, *A new frontier in exciton transport: Transient delocalization*, J. Phys. Chem. Lett. **13**, 6820 (2022).

[60] S. Giannini, A. Carof, and J. Blumberger, *Crossover from hopping to band-like charge transport in an organic semiconductor model: Atomistic nonadiabatic molecular dynamics simulation*, J. Phys. Chem. Lett. **9**, 3116 (2018).

[61] A. Carof, S. Giannini, and J. Blumberger, *Detailed balance, internal consistency, and energy conservation in fragment orbital-based surface hopping*, J. Chem. Phys. **147**, 214113 (2017).

[62] A. Carof, S. Giannini, and J. Blumberger, *How to calculate charge mobility in molecular materials from surface hopping non-adiabatic molecular dynamics—beyond the hopping/band paradigm*, Phys. Chem. Chem. Phys. **21**, 26368 (2019).

[63] J. C. Tully, *Molecular dynamics with electronic transitions*, J. Chem. Phys. **93**, 1061 (1990).

[64] S. Giannini, W.-T. Peng, L. Cupellini, D. Padula, A. Carof, and J. Blumberger, *Exciton transport in molecular organic semiconductors boosted by transient quantum delocalization*, Nat. Commun. **13**, 2755 (2022).

[65] W. T. Peng, D. Brey, D. Dell'Angelo, S. Giannini, I. Burghardt, and J. Blumberger, *Exciton dissociation in a model organic interface: Excitonic state-based surface hopping versus multi-configurational time-dependent Hartree*, J. Phys. Chem. Lett. **13**, 7105 (2022).

[66] G. A. Worth, H.-D. Meyer, H. Koppel, L. S. Cederbaum, and I. Burghardt, *Using the MCTDH wavepacket propagation method to describe multimode non-adiabatic dynamics*, Int. Rev. Phys. Chem. **27**, 569 (2008).

[67] F. Gajdos, S. Valner, F. Hoffmann, J. Spencer, M. Breuer, A. Kubas, M. Dupuis, and J. Blumberger, *Ultrafast estimation of electronic couplings for electron transfer between conjugated organic molecules*, J. Chem. Theory Comput. **10**, 4653 (2014).

[68] O. G. Ziogos and J. Blumberger, *Ultrafast estimation of electronic couplings for electron transfer between pi-conjugated organic molecules*, II, J. Chem. Phys. **155**, 244110 (2021).

[69] Z. Futera and J. Blumberger, *Electronic couplings for charge transfer across molecule/metal and molecule/semiconductor interfaces: Performance of the projector operator-based diabatization approach*, J. Phys. Chem. C **121**, 19677 (2017).

[70] I. Kondov, M. Cizek, C. Benesch, H. Wang, and M. Thoss, *Quantum dynamics of photoinduced electron-transfer reactions in dyesemiconductor systems: First-principles description and application to coumarin 343 – TiO₂*, J. Phys. Chem. C **111**, 11970 (2007).

[71] J. P. Perdew, K. Burke, and M. Ernzerhof, *Generalized gradient approximation made simple*, Phys. Rev. Lett. **77**, 3865 (1996).

[72] J. P. Perdew, M. Ernzerhof, and K. Burke, *Rationale for mixing exact exchange with density functional approximations*, J. Chem. Phys. **105**, 9982 (1996).

[73] A. Kubas, F. Gajdos, A. Heck, H. Oberhofer, M. Elstner, and J. Blumberger, *Electronic couplings for molecular charge transfer: Benchmarking CDFT, FODFT and FODFTB against high-level ab initio calculations*, II, Phys. Chem. Chem. Phys. **17**, 14342 (2015).
[74] A. Kubas, F. Hoffmann, A. Heck, H. Oberhofer, M. Elstner, and J. Blumberger, *Electronic couplings for molecular charge transfer: Benchmarking CDFT, FODFT, and FODFTB against high-level ab initio calculations*, J. Chem. Phys. **140**, 104105 (2014).

[75] J. Hutter, M. Iannuzzi, F. Schiffmann, and J. VandeVondele, *CP2K: Atomistic simulations of condensed matter systems*, WIREs Comput. Mol. Sci. **4**, 15 (2014).

[76] G. Londi, S.-U.-Z. Khan, L. Muccioli, G. D'Avino, B. P. Rand, and D. Beljonne, *Fate of low-lying charge-transfer excited states in a donor:acceptor blend with a large energy offset*, J. Phys. Chem. Lett. **11**, 10219 (2020).

[77] R. C. I. MacKenzie, T. Kirchartz, G. F. A. Dibb, and J. Nelson, *Modeling nongeminate recombination in P3HT: PCBM solar cells*, J. Phys. Chem. C **115**, 9806 (2011).