A n-vector model for charge transport in molecular semiconductors

The rational control of charge or excitation transport underlies all organic semiconducting technologies.^1–3 While inorganic active layers rely upon the precise atomic positioning of a crystal lattice to transport charge over $μ$m length scales, organic semiconductors must establish percolative charge transport networks in structurally disordered morphologies controlled by van der Waals interactions.^4,5 Complicating matters further, solution-processability in organic semiconductors is often achieved via the attachment of non-conducting alkane chains to the π-electron system, disrupting intermolecular aggregation resulting from cohesive $π−π$ stacking interactions.⁶ Unfortunately, while interrupting π-stacking improves solubility, it is detrimental to optoelectronic functionality.⁴ Despite these issues, many organic semiconductors are robust to structural disorder and effectively transport charge over device-relevant length scales, as has been demonstrated by the recent proliferation of high-performance, non-crystalline materials.^7–10 

In many instances, the rationale behind why one organic semiconductor readily transports charge better than another is at best inconclusive. Strong local order of the molecules or even crystallinity has been recognized as a quality that would lead to efficient charge transport,^5,10,11 but there is scant evidence on the relationship between the manner of order (e.g., dimensionality, topology, periodicity) and device performance. In addition, the defect tolerance of the ordered material may play the most important role in the transport of charge.¹² 

Molecular planarity is a common feature among many organic semiconducting molecules due to the π-conjugation it affords. Conventional planar molecules, e.g., polyacenes and pyrenes showed promise as high electron mobility materials because of their capacity to strongly couple to neighboring molecules due to the π-orbital overlap. On the other hand, derivative molecules such as rubrene,^13,14 TIPS pentacene,¹⁵ and benzothiophene heterocyclic oligomers^16,17 rely less on just co-facial π-stacking along a single direction due to their herringbone molecular packing leading to 2D transport networks.^14,18 For more information on the 2D nature of many organic single crystal transport networks, we direct the reader to excellent discussion on the factors influencing charge mobility in the work of Coropceanu et al.¹ 

Due to molecular planarity, many organic semiconductors are capable of strong intermolecular coupling in only one spatial dimension, though much thought has gone into bypassing this limitation, i.e., 2D and 3D couplings.¹⁹ Conventional planar molecules in an ordered, condensed phase can be described as having “1D coupling topologies” (Figure 1) since there is no appreciable electronic coupling between valence orbitals in any direction other than co-facial stacking. An elementary result of bond percolation theory states that a one-dimensional network of randomly connected sites will only percolate at a percolation threshold (p_c) of p_c = 1; a one-dimensional network will only percolate if all sites are connected. Consequently, molecules with 1D coupling topologies, which form 1D columnar stacks,^20,21 can be compromised by a single break in the chain of connections, and their ability to transport charge is not robust to structural, specifically rotational, disorder. Higher dimensional coupling topologies mitigate the fragility of the 1D networks with the additional charge transport pathways they create, yet the design principles have not been articulated. For instance, grain size and surface energy of higher dimensional coupled molecules are a limiting factor for electron mobilities.²² Using molecular networks to describe π-electron coupling has been effective in describing charge transfer pathways and predicting percolation thresholds for different molecular topologies.^4,5,23

Given the weakness of van der Waals interactions (⁠$∼kBT$⁠), the prospect of ordering molecules over macroscopic length scales at room temperature using solution-processing appears daunting. Regardless, there has been a concerted effort to solve this problem by increasing π-system ordering via the incorporation of extended fused-ring π-systems to enhance $π−π$ stacking interactions.^24,25 Often these crystallinity-inducing modifications can have a detrimental effect on molecular solubility, as the forces conducive to crystallinity in the film simultaneously increase the cohesive energy density of the material.²⁶ To make these molecules more soluble, alkyl side-chains must be added to disrupt aggregation, which often negates the improved charge transport properties. While alkyl chains, by their nature insulating, are broadly considered detrimental to transport, it is important to note that it has recently been shown by Illig et al.¹⁷ that dynamic disorder can be fine-tuned and reduced via the proper selection of side-chains, enabling the development of champion materials.

The rational design goal of creating more ordered morphologies of organic semiconductors is contingent upon the assumption that crystalline materials should exhibit superior charge transport when compared to non-crystalline materials. While this wisdom holds true for many materials, numerous case studies have emerged where amorphous materials exhibit exceptional charge transport properties.⁹ In the field of non-fullerene acceptors for organic photovoltaic (OPV) applications, a transformation in molecular design principles appears to be occurring along these lines.^7,27,28

In this paper, we introduce a new methodology which maps molecular structure to a n-vector spin model for use in charge transport simulations of molecular semiconductors. Using the relative orientations between two 3D vectors as a proxy for the π-orbital overlap (and consequently the transfer integrals), we explicitly compute charge transport as a function of molecular topology and bulk structural disorder. We focus on the roles of long-range orientational order in charge transport and conclude with a summary of design goals for organic semiconductors.

Guided by the relationships between molecular topology, orientational order, and bulk charge transport, this article accomplishes the following goals:

• Introducing concepts known as the molecular interaction topology (MIT) and the molecular coupling topology (MCT), we delineate when long-range orientational order is beneficial and deleterious to charge transport.

• We articulate and quantify the contributions of MCT to establishing percolative networks in structurally disordered environments. Additionally, we show that MCT with $d> 1$ (2D, 3D) coupling dimensionality leads to higher charge diffusion lengths.

• We identify a critical interaction strength (Ec) between molecules of $Ec≈ 4− 5kBT$⁠, where order-disorder transitions significantly influence charge transport.

• We explicitly outline the roles of local (orientational) and global (domain boundaries) disorder and their impact on charge transport. We identify two distinct contributions to charge transport, charge transport network topology, and charge transport connection strength (charge transfer rates), and deconvolve them accordingly.

• We develop an understanding of the relative impacts of the local molecular environment on site energies and transfer integrals of molecular materials. Specifically, we find that the average site-energy difference is not a monotonically decreasing function of long-range crystalline ordering.

We begin by proposing that an arbitrary molecule can be represented by two sets of spatial 3-vectors sharing both a fixed origin and relative orientation, but possessing complete rotational degrees of freedom along all three Euler angles (Figures 1 and 2). One set of spatial vectors is named the “molecular coupling topology” (MCT), denoted by ${a1𝐮1^,a2𝐮2^,a3𝐮3^}$⁠. The other is termed the “molecular interaction topology” (MIT), denoted as ${b1𝐯1^,b2𝐯2^,b3𝐯3^}$⁠. In these definitions, all elements $𝐮i^,𝐯i^$ are orthogonal basis unit vectors, whereas all coefficients a_i,b_i are defined to be positive real numbers. The MCT represents how many spatial dimensions a molecule’s relevant orbital(s) can efficiently couple to a neighbor’s orbital; practically speaking, this represents the direction of the molecule with the greatest v-electron density, or, for the majority of cases, the vector normal to the π-system of a planar molecule (see Figure 1). As an example, a₁ = 1 and a₂ = a₃ = 0 would represent a 1D coupling topology (pyrene). The MIT analogously represents the spatial directions in which energy lowering intermolecular interactions are strongest: again, in the case of a simple π-system that possesses no strongly interacting electrostatic groups, this is represented by the vector normal to the π-system of the molecule, indicating that $π−π$ dispersion interactions are the dominant ordering force for the molecule. In the case that some other molecular element influences the relative orientation of lowest energy, the MIT will point in the direction of that spatial preference (see Figure 1), with the MCT rotating along with the MIT. When the MCT and the MIT point in the same direction, we say they are coincident. A set of examples clarifying these definitions is shown in Figure 1, with more detailed examples provided in the supplementary material.

Having defined a representation of individual molecules using the MCT and MIT, expressions for the intermolecular energy and intermolecular electronic coupling between two arbitrarily oriented molecules must be established. To this end, we choose a simple generalized dot-product expression that incorporates the lattice translation vector t_L corresponding to the spatial separation between each molecule’s center of mass (COM) (the origin of the MIT/MCT vectors). These expressions have been used in previous work,²³ and qualitatively reproduce the behavior expected for two interacting π-electron systems,

The expression for the intermolecular electronic coupling (⁠$τij$⁠), Eq. (1), is rationalized by the fact that the interaction between two MCTs is the same regardless of a 180° rotation of either MCT; the form of Eq. (1) is bounded by 0 and $τmax$⁠, and each term in the sum of Eq. (1) behaves like $cos2[θ1]cos2[ϕ1]$ in accordance with the angle dependence of the overlap of two p-orbitals (Figure 2). This behavior has been recently confirmed by quantum-chemical studies on representative molecules at relevant intermolecular separations.²⁹ The expression for the interaction energy between two MITs (Eq. (2)) is validated by the similarity of its orientational dependence to the common expression for the energy of interaction between two dipoles.³⁰ The form of Eq. (2) models the effects of static disorder resulting from the misalignment of nearest-neighbor electrostatic dipoles, convoluted with short-range, directional dispersion interactions between π-electron systems.

A bulk material is assembled by copying the coarse-grained MCT/MIT representation of our molecule onto a lattice. We use the simplest possible bulk representation with a periodic cubic lattice of length L, each site having six nearest neighbors. Formally, the model presented is a version of a n = 3 O (n) spin model on a cubic lattice.³¹ It is informative to write down the Hamiltonian for the O(n) spin model, which can then be solved for the partition function. A critical temperature can then be derived from the partition function, which will serve for us as a prediction of the order-disorder transition of the model.

The Hamiltonian of a n-vector spin system is captured by

where K_ij is the nearest-neighbor (positive) spin coupling constant and S_i is the spin state vector. The partition function can be constructed in the usual manner, but here we show the partition function that has been integrated over orientational phase space and averaged over all spin states

where $Ω$ is the phase space volume $Ω=∏i∫dΩi$⁠. From the moment theorem and Eq. (6), we can show the characteristic function of the n moments of the spin vector $⟨Sn⟩$ as $f(x)=sin(3)x)(3)x$⁠, which shows oscillatory behavior and non-zero 4th-order (and higher) moments.³² The fact that such higher-order moments exist in the characteristic function of Eq. (6) shows that the correlations in molecular orientation can intersect resulting in interactions of the fluctuations in the clusters’ molecular orientation. Finally, the n-vector model predicts the order-disorder critical temperature T_c. For n = 3 and up to 4th-order³³ we find

with $Tc0=2dEmax/kB$ and d = 3 being the dimension of the lattice. We translate Eq. (7) to a critical coupling energy E_c in units of k_BT,

This leads to a value of E_c = 4.58 k_BT, which (as we show in Section III) predicts well the onset of molecular ordering in our system.

A primary deficiency of our presented model is the inability of a rigid lattice to describe local density fluctuations. However, we believe that the insight generated via the orientational degrees of freedom obtained by fixing the COM of each molecule to a lattice site is highly valuable. In the future, density fluctuations could be incorporated by introducing distance-dependent interactions (i.e., a Lennard-Jones potential), and performing off-lattice simulations akin to conventional molecular dynamics simulations.

We utilize a cubic lattice length of 25 sites with periodic boundaries (∼15 000 sites). Larger lattices (10⁶ sites) were also studied, but produced results identical to those conducted at a lattice size of L = 25 sites. The L = 25 lattice is large enough to approach experimentally relevant film sizes (∼ 25 nm assuming a lattice spacing of ∼ 1 nm), while being tractable enough to have adequate sampling. With the descriptions of the lattice outlined in Eqs. (1)–(4), we have the ability to introduce an arbitrary degree of correlation in the disorder of both the energies and transfer integrals as a function of the molecular orientations of each “molecule.” Since our lattice system lacks natural dynamics due to the MIT/MCT abstraction, we utilize a Monte Carlo (MC) scheme to sample the canonical distribution of molecular orientations on the lattice.³⁴ The physical state of the system, defined by the set of {3L³} Euler angles on the lattice, is generated by MC moves to propagate a Markov chain of states. Random trial moves of the molecular orientation are proposed and accepted using the Metropolis criteria. This Markov chain is constructed by proposing local, single-site configuration changes between an old configuration “o” and a new configuration “n.” These transitions are accepted according to the probability of the following equation:

where $ΔEiint=Eiint(θxn,θyn,θzn)−Eiint(θxo,θyo,θzo)$ is the energy difference between the new and old configurations, and $β$ is the inverse temperature. Trial displacements are proposed by applying three consecutive random rotations to the MCT and MIT coordinate systems around fixed x, y, z axes,

Equation (10) shows the MC updating scheme of the three angles for each site; the MC algorithm performs trial moves for each of the angles. The reverse trial move according to this scheme is equally probable, guaranteeing that this algorithm fulfills detailed balance. This algorithm for trial orientations results in a uniform sampling of all spherical orientations, proof of which is provided in the supplementary material. Due to the relatively small lattice sizes examined in this work, our sampling method appears sufficient. In systems examining significantly larger lattices, Swendsen-Wang,³⁵ Wolff cluster,³⁶ or orientational biasing algorithms would be preferable,^37,38 and will be incorporated in future work. Lattice crystallinity is controlled by the dimensionless value of $E¯max$$=|Emax|β$ at which the lattice is simulated, which governs the acceptance of trial moves.

To account for the presence of both equilibrium and non-equilibrium morphologies (as are common in solution-processed molecular materials), we have applied two separate MC schemes with different initial conditions. Specifically, the presence of non-equilibrium morphologies is often indicated by the presence of domain boundaries (DB), which are deleterious to charge transport.^1,39–42 To simulate systems that are DB free, we begin our MC orientational sampling from the initial state of a perfectly aligned crystal. This allows us to avoid low-energy trap states at interfaces between large ordered regions of different relative orientations. To incorporate systems with DB, we start the MC orientational sampling from the initial state of a distribution of randomly oriented molecules on the lattice. This leads to rapid quenching of the lattice in the MC simulation, resulting in the development of a polycrystalline system with many large interfaces between ordered domains. While in principle enough MC steps should converge both systems to the same end-state, we limit our number of steps at 10⁷ MC steps to have well-equilibrated polycrystalline systems in deep local energy minima in the case of systems with DB.

To explicitly characterize the mobility of charges in this lattice model, we utilize a kinetic Monte Carlo (KMC) algorithm to solve the master equation for the lattice-directed graph, with edges defined by charge transfer rates derived from suitable rate theories. KMC charge hopping simulations are initialized by starting a charge on a random molecule in the lattice, and propagating it through the series of lattice sites using the KMC algorithm. Specifically, we use two forms for the rate expression, a Fermi’s Golden Rule (FGR) expression (Eq. (11)) which isolates only the topological (transfer integral) contributions to charge transport on the lattice, as well as a Marcus-Jortner expression⁴³ (Eq. (12)) which explicitly incorporates local energetic disorder, polaronic contributions, and specific molecular identity. In Eqs. (11) and (12), $τij$ is the transfer integral between sites i and j, $ℏ$ is Planck’s constant, $ρ$ is the density of states, $λext$ is the external reorganization energy, $ν$ is the number of vibrational quanta associated with the effective intramolecular mode governing internal reorganization $ωeff$⁠, S_eff is the effective Huang-Rhys factor associated with $ωeff$⁠, k_B is Boltzmann’s constant, T is temperature, and $ΔEijint$ is the energy difference between sites i and j (⁠$ΔEijint=α(Ejint−Eiint)$⁠). $α$ is a scale factor that sets the maximum value of $ΔEijint$ to be substituted into the rate theory,

To characterize the mobility of charges on the lattice, we examine the ensemble-averaged mean charge diffusion distance (Eq. (13)),

where r(x,y,z;t) is the lattice position of the charge carrier at time t. The first brackets represent the average value of $|Δr|$ over an ensemble of KMC trajectories, and the second brackets represent an ensemble average over independent morphologies conducted at the same value of $E¯max$⁠. For each data point, the simulations are averaged over 200 independent lattices, with each lattice annealed for 10⁷ MC steps, yielding acceptance rates between 10% and 50%. KMC trajectories are run for 10 ns and averaged over 20 000 independent trajectories per independent lattice configuration. A maximum value of the intersite transfer integral (⁠$τmax$⁠) is set to 0.025 eV, which satisfies the perturbative criteria of Eqs. (11) and (12). As the model described in this work is primarily phenomenological, i.e., we are using it to explore the general properties of charge transport in molecular materials, the authors by no means expect to quantitatively reproduce experimental charge mobility data with our simulations. Of course, if one were to release the lattice constraints and simulate a more realistic system, one could relate $⟨⟨Δr2⟩⟩$ to the zero-field mobility via an Einstein relation for the 3-dimensional diffusion coefficient. For FGR parameters, we assume a density of states value, $ρ$⁠, of 1 state/eV. For Marcus-Jortner parameters, we utilize physically reasonable values for small molecular semiconductors consistent with the literature (⁠$Seff= 1.4,λi= 0.2eV,λe= 0.05eV,ωeff=0.18eV$⁠).^44–46 In the calculation of the site-energy differences using Eq. (3), we scale the value of $Eiint$ by $α$ to ensure a range of site-energy values between 0 and 0.3 eV as input into the rate-theory expressions.

We begin with the examination of charge diffusion lengths in disordered systems parameterized with FGR rates (Eq. (11)). As FGR assumes energetic alignment, energetic site disorder (Eq. (3)) is not taken into account, and all sites are assumed to be isoenergetic. As a result, all intermolecular charge transfer rates are only dependent upon the value of the intermolecular transfer integral, which is related to the relative orientations of two lattice sites via Eq. (1). To state concretely, this approach does not take into account energetic disorder-limited transport, but will help deconvolve contributions to charge transport when energetic disorder is included in Sec. III B.

To examine how charge transport is influenced by the strength of intermolecular interactions, we plot $⟨⟨|Δr|⟩⟩$ for coincident 1D systems in Figure 3 (top) as a function of $E¯max$⁠. For clarity, we note that $E¯max$ is used as a measure of the propensity of a system to crystallize, including both temperature and the inherent strength of intermolecular interactions. Additionally, we emphasize that in the case of coincident MCT and MIT, the orientation of lowest energy is also that of strongest electronic coupling.

In Figure 3 (top) an increase in $⟨⟨|Δr|⟩⟩$ is observed as $E¯max$ increases from 0 to 4. As $E¯max$ approaches E_c, a stark change in the behavior of $⟨⟨|Δr|⟩⟩$ is observed in both 1D systems with and without DB. The system in which DB formation is forbidden exhibits a sharp increase in $⟨⟨|Δr|⟩⟩$ of nearly 100%, indicating a transition from disordered to well-aligned regions of molecules, whereas the system in which DB formation is allowed only slowly rises as $E¯max$ increases past E_c, indicating increased local alignment, followed by a drop-off at $≈Ec$ due to the development of strong DBs which limit long-range charge transport. This result is worth focusing on briefly as it indicates that the charge diffusion capability of a molecular material exhibits a sharp transition at a critical value of the crystallinity-inducing forces (E_c), and is not simply a linearly increasing function of intermolecular interaction strength. The ramifications of this fact will be mentioned in Sec. IV.

For $E¯max$ $>Ec$⁠, both systems with DBs and without DBs steadily decrease as a function of increasing $E¯max$⁠, implying that as the system approaches complete crystallinity, $⟨⟨|Δr|⟩⟩$ decreases. While the explanation for this result is straightforward for the system with DB (higher crystallinity-inducing forces leads to sharper, stronger DB), in the DB-free system this outcome is counterintuitive. To understand this result, it is important to consider the dimensionality of $⟨⟨|Δr|⟩⟩$ relative to that of the percolation pathways in the system. $⟨⟨|Δr|⟩⟩$ is a measure of the 3-dimensional distance a charge travels over a specified time-interval. In our 1D MCT/MIT, the strong charge transport pathways that will form are inherently one-dimensional⁴ due to the topology of the molecule (planar π-electron systems form strong 1D percolation pathways when they are stacked in 1D). Consequently, as $E¯max$ increases, the pathways that are formed become strongly one-dimensional. The results of Figure 3 (top) indicate that the charge diffusion capabilities of a system are increased if there is some disorder in the ordered 1D structural domains, allowing charges to hop between columns in the other two spatial dimensions, increasing the effective spatial dimensionality of the percolation pathways, and increasing $⟨⟨|Δr|⟩⟩$⁠. Molecules with higher spatial dimensionality charge transport networks will exhibit increased $⟨⟨|Δr|⟩⟩$⁠. We next test this hypothesis by increasing the number of spatial coupling dimensions in the MCT to 2.

In Figure 3 (bottom) we examine the dependence of $⟨⟨|Δr|⟩⟩$ on $E¯max$ in systems with 2D MCT and MIT. The most immediately observable feature of transport in these systems is that the baseline value of $⟨⟨|Δr|⟩⟩$ is an order of magnitude larger than in the case of molecular materials with 1D MCT and MIT (Figure 3 (top)). This suggests that the ability of molecules to electronically couple in multiple spatial directions is beneficial to transport, as it increases the spatial dimensionality of the percolating charge transport network, in agreement with our previous work.⁴ Moreover, the baseline value of $⟨⟨|Δr|⟩⟩$ is significantly higher than that of even the fully aligned 1D MCT/MIT systems, suggesting that the ability of molecules to couple in multiple spatial directions, even if their ability to crystallize is poor, is essential to achieving high charge carrier diffusion distances/charge mobilities and may in some circumstances counteract small transfer integrals between molecules.

In 2D systems without DBs, a large increase in $⟨⟨|Δr|⟩⟩$ is observed as $E¯max$ increases past a value of E_c, while in 2D systems with DBs $⟨⟨|Δr|⟩⟩$ slowly decreases after a $E¯max$ value of $≈Ec$ is reached. The relatively sharp increase of $⟨⟨|Δr|⟩⟩$ in DB-free systems at $E¯max$ values of 4.5 corroborates the same critical phenomenon observed previously for the 1D systems. The decreasing value of $⟨⟨|Δr|⟩⟩$ in the 2D system with DBs falls in line with the idea that DBs continue to limit transport in strongly ordered systems. However, the sudden drop-off in $⟨⟨|Δr|⟩⟩$ at $E¯max>Ec$ is not as severe because 2D MCT resists DB limited transport due to their ability to couple in multiple spatial directions; whereas 1D MCT/MIT molecules form columns as their characteristic percolation pathways, 2D MCT/MIT molecules form sheets, and all that is needed to connect two unconnected 2D networks is a single disordered segment between the two. These results support the idea that molecules with higher-dimensional MCT should not only exhibit increased charge diffusion, but also should more effectively resist DB-induced mobility decreases than 1D MCT molecules.

The case of non-coincident 1D MCT/MIT is especially interesting in the context of crystallinity (Figure 3 (top)). When the direction in which a molecule prefers to align is not coincident with the direction in which it electronically couples most strongly, increasing crystallinity (⁠$E¯max$⁠) does not lead to increasing $⟨⟨|Δr|⟩⟩$⁠, as observed in Figure 3 (top). The reason for this is simple: the lowest energy morphology of the non-coincident MCT/MIT molecular representation is not the morphology with the highest intermolecular couplings. Alternatively put, the MIT induces the formation of columns of molecules; however the transfer integrals within these columns of oriented molecules are all essentially zero. This results as a nearly flatlined value of $⟨⟨|Δr|⟩⟩$ in Figure 3 as a function of $E¯max$ that is relatively insensitive to the details of the morphology; since short-range charge transport is already compromised, DB’s do not have a noticeable effect on the system’s charge transport properties over the time scale of our simulations (10 ns).

2D non-coincident MCT/MIT exhibits a similar crystallinity dependence as the 1D system for the DB-free case: increasing crystallinity leads to the formation of ordered sheets of molecules, where increasing $E¯max$ beyond E_c does not lead to increased charge transport. Curiously, in non-coincident 2D systems with DBs, an increase in $⟨⟨|Δr|⟩⟩$ is observed as a function of $E¯max$ $>Ec$ (Figure 3). This effect is a side-effect of our definitions of the coupling topologies. In the case of 2D non-coincident systems without DB, the sheets of molecules align via the MIT, but since the MCT is partially coincident with the MIT, the formation of DB actually leads to coupling within and between sheets, leading to a 3D transport network. This behavior highlights the non-trivial behavior of charge transport as a function of molecular topology and interaction strengths. This feature corroborates the sometimes non-intuitive effects of crystallinity in molecular systems. If the MCT and MIT are aligned in a certain fashion, DB development can theoretically be advantageous to long-range charge transport. It is for this reason that we encourage the reader to think carefully about the impact of crystallinity in their systems.

While Sec. III A was concerned with topological contributions to charge transport (both the strength of interactions and the spatial dimensionality of the charge transport network), now we utilize a Marcus-Jortner rate theory⁴³ incorporating all energetic and polaronic contributions to charge transport. In Figure 4 we plot $⟨⟨|Δr|⟩⟩$ using rates computed with Eq. (12).

In Figure 4 (top) it is immediately obvious that the inclusion of energetic disorder decreases the values of $⟨⟨|Δr|⟩⟩$ in all systems relative to those in Sec. III A (Figure 3), confirming the dominant role that energetic and polaronic contributions play in non-crystalline systems. Notably, we observe similar increases in the value of $⟨⟨|Δr|⟩⟩$ when $E¯max$ $>Ec$⁠. This suggests that in systems where the formation of DBs can be prevented, increasing the ability of a 1D MCT/MIT molecule to crystallize can be very beneficial to charge mobilities. In the case of the 1D system with allowed DB formation, $⟨⟨|Δr|⟩⟩$ exhibits similar behavior as in the FGR cases, with a slight increase in $⟨⟨|Δr|⟩⟩$ at $E¯max$ $<Ec$ due to an increase in local ordering, following by a fall off due to the formation of rigid, transport-limiting DB.

Though mitigated in the case of including energetic disorder, the slight fall off in the values of $⟨⟨|Δr|⟩⟩$ as $E¯max$ increases significantly above 5 is still observed. The mitigation of $⟨⟨|Δr|⟩⟩$’s decrease beyond $Ec≈5$ relative to the FGR system is likely due to the fact that energetic and polaronic contributions will slow down transport along the 1D percolation paths (columns), making hops between columns more likely to occur than in the purely FGR case, increasing the effective spatial dimensionality of the charge transport network. In other words, if k_intra (hopping within columns) $>>$ k_inter (hopping between columns), charges only move within 1D columns, and the transport network dimensionality is 1. However, if $kintra≈kinter$⁠, the column is not the preferred means of transport, and the transport network dimensionality becomes $>$ 1. As a consequence, while the individual site connections in the network might be weaker than in the FGR case, the effective spatial dimensionality of the network has increased, and these two effects move $⟨⟨|Δr|⟩⟩$ in opposite directions.

To deconvolve the energetic and coupling contributions to $⟨⟨|Δr|⟩⟩$⁠, we plot the ensemble-averaged standard deviation of the lattice site energies as a function of $E¯max$ in Figure 5. Here, we are using

as a metric for the average site-energy difference, where E^avg is the ensemble-averaged value of the site energy, which is consistent with the band tail description of charge transport.³ In the case of the 1D coincident DB-free system, a steady increase in $⟨σE⟩$ is observed for $0<E¯max<Ec$⁠, implying that as the system becomes more ordered, the average site energy difference increases. In the case of 1D coincident systems with DBs, a less-smooth plot of the energetic deviation is observed, resulting from the fact that the system is trapped in many different high-energy, non-equilibrium morphologies. This is further showcased by the fact that $⟨σE⟩$ in the system with DBs appears to converge to a higher value relative to the system without DBs, confirming that the average site energy differences in systems with DBs limit transport.

The non-monotonic behavior of $⟨σE⟩$ in Figure 5 (top) warrants comment. The result that $⟨σE⟩$ exhibits a maximum as a function of $E¯max$ can be rationalized: at $E¯max$ = 0, the “effective dipoles” on every site in the system are orientationally randomly distributed. Consequently, every site exhibits an interaction energy averaged over 6 randomly oriented dipoles. As a result of this, every site energy is similarly random, and $⟨σE⟩$ is relatively small. However, as $E¯max$ increases, $⟨σE⟩$ increases because orientational correlations of dipoles begin to develop, and the local dipolar environment around each site becomes non-uniform; ordered regions form along with amorphous ones, leading to a broader distribution of dipolar environments, leading to a larger value of $⟨σE⟩$⁠. However, after a critical interaction strength, $E¯max$ = E_c, is reached, the system becomes able to fully orient itself into a low-energy, ordered state, where every site is in its lowest energy state, and thus $⟨σE⟩$ trends towards zero.

The case of 2D systems exhibits features similar to that in the discussion of Sec. III A 2 utilizing FGR rates: $⟨⟨|Δr|⟩⟩$ is overall larger than in the 1D case, and $⟨σE⟩$ is overall lower. These two effects can be attributed to the increase in the spatial dimension of the charge transport network mentioned previously as a result of increasing MCT dimensionality, and a more uniform distribution of randomly oriented dipoles (every site now has two orthogonal, randomly oriented dipoles on it, leading to a total of 12 randomly oriented dipole interactions per site), respectively. Consequently, the distribution of site-energies is more uniform, and $⟨σE⟩$ is lower. Regardless, a noticeable transition in the $E¯max$ $≈Ec$ region is still observed.

Finally, in Figure 6, we examine the coupling contributions to $⟨⟨|Δr|⟩⟩$ as a function of $E¯max$⁠. We note that this plot applies both to the Marcus-Jortner results of this section, as well as the FGR results of Sec. III A. Figure 6 exhibits expected behavior, with an increase in the average value of the transfer integral as $E¯max$ increases, indicating that in 1D coincident molecules, increasing molecular orientation leads to larger values of the transfer integral. Note that in non-coincident systems, such a trend will not be observed except in exceptional circumstances. Again, we observe a relatively sharp increase in the value of $τij$ as $E¯max$ approaches E_c, consistent with the previous suggestions of a critical value of $E¯max$ where the system begins to form strong orientational correlations.

It is a common maxim that crystallinity is beneficial for charge transport in molecular and polymeric systems. In this work, we coarse-grain molecular detail down to what we deem to be the simplest possible representation necessary for understanding the relationship between charge transport and crystallinity: the MCT and MIT. While the interactions in this work are modeled as simple spin vector couplings and interactions, the simulations provided suggest that the influence of crystallinity on charge transport in molecular materials is complex, and that to fully understand the impact of crystallinity requires an inclusion of basic molecular topological features and explicit structural disorder.

First, there are two important contributions to charge transport: the strength of connections within a charge transport network and the spatial dimensionality and extent of that charge transport network. When examining the simple metric of $⟨⟨|Δr|⟩⟩$ in this work, it becomes clear that large transfer rates between sites and increasing the dimensionality of a molecule’s coupling topology are both paths to increased $⟨⟨|Δr|⟩⟩$⁠. Moreover, the increased spatial dimensionality of 2D MCT/MIT allows for a more robust resistance to DB-induced mobility decreases (Figure 3(b)). In the case of coincident, DB-free, 1D MCT/MIT systems (Figure 3(a)), the effects of charge transport network spatial dimensionality can be clearly observed as a decrease in the value of $⟨⟨|Δr|⟩⟩$ as a function of $E¯max$ $>Ec$⁠. This somewhat counterintuitive result suggests that careful attention should be paid to the MCT of a molecule, as well as the topology of the formed charge transport network. Specifically, in Sec. III B we argued that slower rates within 1D columns can lead to an effective increase in the spatial dimensionality of the charge transport network, two effects that have opposite impacts on $⟨⟨|Δr|⟩⟩$⁠.

Second, we note that the value of the standard deviation of the mean site-energy (⁠$⟨σE⟩$⁠) as a metric for the energetic offset between sites does not uniformly decrease as a function of $E¯max$⁠. Instead, a rise is observed until a critical value of E_c = 4.58 k_BT, followed by a decrease to a lower value (see Figure 5). This is due to the fact that in an effective dipolar model for the interaction energies, the isotropically distributed random dipoles should have a lower average site-energy difference than a lattice that is a mixture of ordered and disordered regions. As the structure dominates the system and overcomes the random component, $⟨σE⟩$ decreases to its lowest value, in agreement with common intuition about the effects of crystallinity on energetic disorder.

Lastly, we note that the complex influences of molecular topology can lead to highly unusual behavior: take the case of the non-coincident 2D system with DBs. While one would expect DBs to always decrease charge mobilities, if a molecule’s MIT and MCT are aligned in a certain fashion, the presence of energetic DBs can actually lead to a more extended and connected charge transport network as shown in Figure 3.

One of the most striking results of these simulations is the sharp transition in the behavior at critical energy E_c. This suggests that in order to maximize transport in systems capable of reaching equilibrium, a critical interaction strength (forn = 3) of at least $Ec>4.58kBT$ is needed to establish the long-range charge transport networks. Restricting the interactions to a 2D surface leads to a more modest critical intermolecular energy $Ec2D=2.31kBT$⁠. The 3-dimensional value of the interaction strength is at the limit of dispersion interactions; to obtain this magnitude of dispersion interactions, large fused π-electron systems are needed, which implies rigid 1D MCT/MIT, and spatially 1D charge transport networks. To form higher-dimensional charge transport networks, one will need to utilize biaxial interactions which are stronger than the threshold of E_c = 4.58 k_BT, and which do not lead to strictly 1D transport networks. A possible solution is the use of hydrogen bonds; it is perhaps not a coincidence that biological molecules such as protein complexes utilize hydrogen bonds to establish robust charge transport networks capable of self-assembly.^47,48 Hydrogen bonds have been modestly explored to enhance the efficiency of organic semiconductor devices,^49,50 but of course present their own difficulties, such as increased electrostatic disorder or the “locking-in” of poor morphologies.⁵¹ Potential applications may involve weak, non-traditional hydrogen bonds studied by the authors.⁵² This point is less the advocation of the use of hydrogen bonds, and more the understanding of the fundamental limitations of dispersion interactions in achieving the self-assembly of desired morphologies.

In addition to considering the role of the magnitude of E_c in the charge transport of the molecular networks, it is also important to discuss the topological defects that form near $E¯max$ $≈Ec$⁠. In this study, we considered a line of dislocations that formed at low-angle DB. The role that DBs have on transport $⟨⟨|Δr|⟩⟩$ is based on the strength of the dislocations and shifts the order-disorder transitions to lower energies $E¯max$ $<4kBT$⁠. This energetic shift is based on the domain size, which is primarily due to the energetic cost of unbinding the dislocation pairs (higher T) and the energy coupling to the density changes in the DB. From the elastic theory of topological defects, the energetic cost of line dislocations in DBs scales as $∼log(R/b)$⁠, where R is the domain size and b is the magnitude of the Burgers’ vector.⁵³ 

The electronic transport behavior of coupled molecules was investigated with a coarse-grained n-vector lattice model. The impact of coupling and interaction topologies on the order/disorder of the molecules was related to charge transport via electron hopping models using a kinetic Monte Carlo algorithm. Using rational design constructs we concluded that five strategies emerged from our study. (i) Molecular topologies have two decomposable topologies: (inter)molecular electronic coupling topology (MCT) and (inter)molecular interaction topology (MIT). Molecular design strategies should look at both topologies when trying to optimize charge transport. (ii) MCT contributions in disordered molecular networks show higher dimensional coupling, e.g., 2D establish percolative networks more efficiently especially when not coincident to MIT. (iii) We derive a critical interaction strength E_c = 4.58 k_BT between neighboring molecules and we show that charge transport is not guaranteed to be enhanced at interactions above E_c. Our model shows interaction strengths greater than 3 kcal/mol can lead to worse transport especially in 1D coupled networks. (iv) We distinguish between local and long-range orientational (dis)order through the formation of domain boundaries. From the inclusion of domain boundaries, we identify the two multi-scale contributions to transport: network percolation and coupling strength (rates). (v) Local molecular environmental effects are quantified in the model via site-energy dispersion. Using site-energy differences, $⟨σE⟩$⁠, we find that energetic dispersion does not monotonically decay to zero as the system becomes more crystalline. The dispersion results show that optimal molecular design would be having interaction energies just above E_c.

We hope these five strategies can be utilized in the design of new materials for non-crystalline semi-conducting devices. Since defects are the rule and not the exception in these materials, we find that the type of defects and the tolerance need to be at the forefront during the design and in some cases should be embraced, not avoided. An important outcome of that is that the dimensionalities of the percolation pathways are vitally important and many of the current semi-conducting molecules only utilize 1D transport networks leading to easily broken pathways. One strategy to mitigate the fragility of 1D transport is micron-length conjugated nanowires,⁵⁴ a topic that future studies would include.

Charge transport calculations in disordered environments are distinct from crystalline materials in that a fixed atomic snapshot is not sufficient to operate a transport calculation on it. In this paper, we utilize a coarse-grained lattice model to calculate transport coefficients for a wide-range of molecular topologies. Although lattice models have constrained translational periodicity, orientational order should be sufficiently decoupled from translational order, especially in the condensed phases where organic electronics operate, though rigorous verification of this hypothesis will be conducted in future work.

See supplementary material for more examples of different molecular structure and the coarse grained mapping to MCT and MIT vectors. Higher dimensional couplings are also described. Additionally, MC orientational sampling data are presented.

The authors thank Brett Savoie for his thoughtful conversations. This research was supported through the ANSER Center, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award No. DE-SC0001059, and the lab equipment was supported through the Division of Chemical Sciences, Office of Basic Energy Sciences, the U.S. Department of Energy, under Contract No. DE-AC02-06CH11357. K.L.K. would also like to thank the AFOSR MURI Grant No. FA9550-11-1-0275 for their support.