Doping of semicrystalline conjugated polymers: dopants within alkyl chains do it better

Abstract

Depending on the sample preparation protocol, various structures for doped semi-crystalline polymers can be achieved, characterized by dopants either inserted in the alkyl side chains or packed with the conjugated backbone, which ultimately results in very different charge-transport and thermoelectric properties. This work targets such an intricate relationship between structure and properties with accurate hybrid quantum-classical calculations fully accounting for the effect of the environment. By considering representative structures for the crystalline domains of the F4TCNQ-doped PBTTT polymer, our calculations reveal that: (i) The electron affinity (EA) of the dopant is highly sensitive to the position occupied by the molecule in the polymer lamellae, with dopants inserted in the alkyl regions being much stronger electron acceptors than those stacked in the π-conjugated backbones (EA difference > 0.5 eV). (ii) The tiny orbital overlap between dopants in the alkyl regions and the polymer favors integer charge-transfer ground states, while dopants packed with conjugated chains are more inclined to fractional charge transfer. (iii) The Coulomb interaction between the charge carrier on the polymer and the ionized dopants is considerably (∼30%) smaller for dopants in the alkyl regions, pointing to less bound carriers. These findings rationalize the fact that record conductivities are generally associated with dopants inserted in the alkyl chains, raising awareness on the importance of controlling the dopant position in the polymer structure.

---

Gabriele D’Avino

Gabriele D’Avino received his PhD in Materials Science at the University of Parma in 2010. After a first postdoct at the University of Bologna and two Marie Curie fellowships at the University of Liege and at the University of Mons, he joined CNRS with a researcher position at the Institut Neél in Grenoble. His research work focuses on the theoretical modeling of organic functional materials for advanced applications in electronics and energy generation, synergistically combining quantum and classical techniques in a multiscale fashion. The goal of his research is the establishment of the relationships between structure and (electronic, optical, vibrational, dielectric) properties in complex molecular systems with special emphasis on environmental and cooperative phenomena in the solid state.

---

1 Introduction

Molecular doping is arguably the main technique to control charge carriers density and transport properties in organic semiconductors, enabling a large variety of technological applications from optoelectronics to thermoelectricity.^1–3 The microscopic mechanism for molecular doping has been depicted as a two-step process,⁴ consisting of (i) the transfer of an electron between the dopant and the host material, and (ii) the subsequent separation of the charge in the semiconductor from the ionized dopant, resulting in a free carrier. Both steps have attracted considerable interest, motivated by the aspiration for a deeper fundamental understanding of the doping process and of enhancing its efficiency in various applications. The second charge-release step is particularly puzzling if one considers the strong Coulomb forces keeping charges bound in low-dielectric constant organic materials. Explanations in terms of energetic disorder^4–6 and collective screening phenomena⁷ have been recently proposed.

Regarding the first charge-transfer (CT) step, the nature of the ground state of host-dopant CT complexes has been studied experimentally in a large variety of systems, pointing to the two qualitatively different scenarios of fractional versus integer CT.² In the former, only a fraction of an electron is transferred from or to the semiconductor. This condition, sometimes referred to as orbital hybridization scenario, seems far from being ideal for efficient doping as compared to integer CT, which leads to a fully-ionized dopant, with the charge being transferred to the host semiconductor. It has been argued that integer CT is predominant in polymer systems while fractional CT is usually encountered in small molecules,⁸ with some notable exceptions, e.g. pentacene.^9,10

It is presently unclear which among the many possible factors in place (e.g. charge delocalization, nanostructure, disorder, Coulomb interactions) determine the formation of fractional versus integer-CT host–dopant complexes. A common guideline consists in assessing the CT energetics by comparing the ionization potential (IP) and the electron affinity (EA) of the donor and the acceptor materials,¹¹i.e. host semiconductor and dopant in the case of p-type doping considered henceforth. This criterion should be, however, revised in view of the major role played by electrostatic interactions in the solid state,¹² which can affect the EA of common molecular dopants by up to 1 eV depending on the host semiconductor.¹³ Moreover, the crucial role of the excitonic interaction between electron and hole in the formation of integer CT complexes has been highlighted, along with polaronic effects.^9,14 All these factors do critically depend on the molecular organization in the doped semiconductor, calling for a deep understanding of the relationship between structure and electronic properties.

This is particularly challenging in conjugated polymers, due to the variety of possible nanostructures in pristine systems and even more so in doped samples. Many studies thus focused on highly-ordered (semicrystalline) conjugated polymers such as poly(2,5-bis(3-tetradecylthiophen-2-yl)thieno[3,2-b]thiophene) (PBTTT) or poly(3-hexylthiophene) (P3HT), typically doped with strong electron acceptors as 2,3,5,6-tetrafluoro-7,7,8,8-tetracyanoquinodimethane (F4TCNQ) or 1,3,4,5,7,8-hexafluorotetracyanonaphthoquinodimethane (F6TCNNQ). These polymers are characterized by a crystalline lamellar structure, consisting of layers formed by stacked chains of conjugated backbones, alternating with planar regions hosting the alkyl side chains (see Fig. 1). Experimental studies reported various nanostructures formed depending on the sample processing, with dopants inserting either within the π-conjugated backbones^8,15–17 or in the inert alkyl chains,^18–22 or in both regions.^23,24 An overview of available experimental data relating the position of molecular dopants in the polymer structure to the nature of the host–dopant complexes, along with the associated electrical conductivities, in provided in Fig. 2.

Fig. 1 (a) Chemical structures of the PBTTT polymer and F4TCNQ dopant investigated in this work. (b) Crystalline PBTTT lamella viewed from the top of stacked polymer-backbone chains, characterized by alternating conjugated and alkyl regions. The two F4TCNQ molecules illustrate the two positions of the dopant molecule considered in this work. (c) Rendering of the DFT/MM-optimized complexes with dopant in the conjugated region (π structure, top) and in the alkyl chains (α structure, bottom) that have been employed in the calculations. The subsystems described at the quantum level (QM region of our QM/MM calculations) are shown in bold, the classical embedding environment (MM region) is drawn with thin lines.

Fig. 2 Compilation of available experimental data for PBTTT and P3HT, relating the position of molecular dopants in the polymer lamellae (alkyl chains vs. conjugated backbone) to the nature of the host–dopant complexes (integer vs. fractional CT). The electrical conductivities of each sample is expressed through the size of the marker, see legend. The position of different points in each quadrant is arbitrary. We note that the large majority of doped polythiophenes feature integer CT, except for two systems corresponding to dopants inserted in the conjugated backbones. The systems with the largest conductivities are characterized by dopants inserted in the alkyl region and integer CT (upper-right quadrant).

Coming to electronic and transport properties, the few examples of polymers for which a fractional CT has been reported are all characterized by dopants packed within the conjugated backbone.^17,24 On the other hand, independent experimental reports have shown that high conductivities (σ > 100 S cm⁻¹) have been obtained only for samples with F4TCNQ or F6TCNQ inserting in the alkyl chains,^{18,19,21,22,25} while morphologies with dopants packed in between π-backbones attain values not exceeding σ = 2 S cm⁻¹.^15–17 We further note that other systems featuring exceptionally high conductivities, such as those obtained by ion-exchange doping (σ up to 1000 S cm⁻¹)^25,26 or upon incremental concentration doping with FeCl₃ (σ up to 22 000 S cm⁻¹),²⁷ are all characterized by dopants placed in the alkyl region. This set of experimental studies overall indicates a higher doping efficiency when dopants insert in the alkyl chains, as can be appreciated from Fig. 2.

In the present paper, we address the relationship between the dopant position in the polymer lamellar structure and the resulting electronic properties, considering F4TCNQ-doped PBTTT as a prototypical case study (see Fig. 1). Our analysis, based on accurate hybrid quantum/classical (QM/MM) calculations, considers two representative host–dopant morphologies, characterized by dopants inserted in the conjugated and in the alkyl region, showing that the latter configuration is the most favorable one for attaining efficient doping. In particular, we highlight the key role of electrostatic interactions in the crystalline polymer, determining a strong dependence of dopant electron affinity on the position in lamellae, de facto making dopants inserted in the alkyl chains stronger electron acceptors than those in the conjugated backbone.

The paper is organized as follows. Our multiscale QM/MM methodology, merging density functional theory (DFT) and many-body formalisms with classical microelectrostatics (ME), is described in Section 2. The presentation of the results in Section 3 follows a gradual approach, first discussing the energy levels of non-overlapping host and dopant units, and then introducing the intermolecular coupling. In doing so, we illustrate the limitations of conventional closed-shell calculations in describing integer CT ground states, which brought us to switch to spin-unrestricted calculations. The final section discusses our finding in the light of the experimental literature on doped conjugated polymers.

2 Materials and methods

The relationship between structure and electronic properties in doped PBTTT is here investigated with a set of complementary computational techniques, going from force field simulations, to QM/MM techniques combining DFT and Green's function many-body formalisms with a classical ME scheme of atomistic resolution.

2.1 QM/MM setup

The doped samples have been built by introducing a F4TCNQ impurity molecule in the crystal structure geometry of pristine PBTTT with 14 C atoms-long alkyl chains. The lamellar crystal structure of PBTTT and the force field have been reported previously.^28,29 Following an educated guess inspired by structural data,^8,17,18,21 dopants have been introduced either within the conjugated backbone (π structure), or in the region occupied by the alkyl side chains (α structure) in a large periodic system (three layers of nine π-stacked dodecamers) made by replicating the simulated PBTTT crystalline structure. In the π structure, a PBTTT monomer has been removed to allow the introduction of the dopant molecule, featuring a similar size. The resulting structure has then been optimized and subject to a 100 ps-long quenched dynamics (NPT ensemble; p = 1 atm, T = 300 K) from which structures have been optimized every 20 ps. The simulation of the structure for the α system has been more complex due to a larger number of possible intercalation sites. Here, tens of structures have been built by modifying the relative position of the dopant with respect to the alkyl side chains. All these structures have been optimized and subject to successive 500 ps-long quenched dynamics (NPT ensemble; p = 1 atm, T = 300 K; quench every 25 ps) until the energy of the lowest energy structure between two consecutive runs does not decrease anymore. The sample considered in the following corresponds to the lowest-energy structure among all the simulated conformers.

QM/MM approaches require to split the system into a QM region comprising one dopant and part of the polymer host, which we will describe by either DFT or many-body methods, and a MM region that comprises the surrounding organic environment, here described at ME level. The QM/MM partitioning of a system requires the QM and MM regions to be non-overlapping, a criterion which can plausibly be satisfied for systems composed of disjoint molecular units, but that is somewhat arbitrary for a virtually infinite conjugated polymer. Here we opt for a practical solution that consists in retaining within the QM region a finite segment of the polymer chain close to the dopant. In order to determine the actual length of this segment, we run a series of GW^30–33 and Bethe–Salpeter equation (BSE)^34–36 calculations for oligomer lengths ranging up to 5 repeating units, with and without the presence of the dopant. This analysis allowed us to conclude that a PBTTT segment of three repeating units (3-BTTT) represents a good model system, featuring reasonably-converged quasiparticle GW gap and optical excitations. The full detail of these preliminary calculations is reported in the ESI.†

We note that the alkyl chains of the 3-BTTT within the QM region have been replaced by methyl groups and the rest of those side chains has been included in the MM system, except for the bridging CH₂ unit that has been removed to avoid artifacts due to the too short distance between QM and MM atoms. Similarly, the two thiophene rings close to the 3-BTTT, one on each side of the trimer and originally belonging to the same polymer chain, have also been removed. After the QM/MM partitioning, the geometries of the QM region have been relaxed with dispersion-corrected DFT (PBEh*-D3/6-31G* level,^37,38 see below for the definition of the functional), accounting for the interactions with the MM region, described with Lennard-Jones potentials and atomic charges. The two structures employed in QM/MM calculations are shown in Fig. 1c.

2.2 Electronic structure calculations

The electronic structure of the host-dopant complexes has been investigated with our original implementation of QM/MM techniques, merging DFT and the Green's function many-body GW and BSE formalisms,^30,33 with classical polarizable models.^39,40 The main features of the employed methods are summarized below, the details can be found in the original papers.^9,41,42

The ground state of the system is obtained from Kohn–Sham DFT calculations on the QM subsystem, performed in the self-consistent electrostatic field of the polarizable MM environment, described at the atomistic ME level (see below). Self-consistent DFT/ME calculations follow an iterative scheme of cross-coupled quantum and classical calculations, in which each subsystem is computed in the electrostatic potential of the other until convergence.⁴² As we will show in the following, the starting guess of this iterative process may have an impact on the self-consistent solution, signaling a multistable behavior.

The discrimination between stable and metastable states requires comparing the total energy of the entire QM/MM system for the different self-consistent solutions. The total energy has been computed as

where ρ_QM is the charge density of the QM system, q_ni is the permanent charge of the MM atom at position r_ni, and E_n (V_n) is the DFT energy of the fragment n (either the QM system or a molecule in the MM region) subject to the potential V_n (r) due to all the other fragments:where the sum extends over both QM and MM molecules, excluding the fragment n itself. The charge density ρ_n of MM molecules is described in terms of discrete permanent atomic charges and induced dipoles, consistent with the ME model. The present approach can be classified as an additive QM/MM scheme,⁴³ in which the continuous charge density of the QM region interacts with permanent atomic charges and induced dipoles of the polarizable MM environment, without any discretization of the QM charge density and avoiding the double counting of interactions. Data (charge carrier environmental energies of molecular crystals) validating our approach to the calculation of extensive energies in a QM/MM framework are provided in the ESI,† Table S1.

In order to properly describe the relative alignment between the polymer highest occupied molecular orbital (HOMO) and the dopant lowest unoccupied molecular orbital (LUMO) in the solid state, as well as the resulting intermolecular orbital hybridization,^7,9 we adopt a tuned DFT functional whose gap matches that of reference embedded GW calculations. Specifically, we employ the single-parameter global hybrid PBEh* = PBEh(a* = 0.45) functional,³⁷ with 45% of Hartree–Fock exchange, in conjunction with the 6-311G* basis set. The details of the tuning protocol are reported in the ESI.† DFT calculations have been performed with the ORCA software package.⁴⁴

Quasiparticle excitations, corresponding to electron addition and removal energies, have been obtained within the framework of the embedded many-body GW formalism.^41,42 In order to describe the effect of the medium on quasiparticle energy levels, both the ground-state DFT and the subsequent GW calculation, should properly take into account the embedding environment. In this framework, the nth energy level of the embedded QM systems can be conveniently written as

E_n (GW_e|DFT_e) = E_n (GW_g|DFT_g) + Δ^E_n + Δ^P_n (3)

where the subscript g(e) labels a DFT or GW calculation performed for the QM molecule in the gas phase (embedded in the solid). E_n (GW_g|DFT_g) hence corresponds to the gas-phase energy which is evaluated with self-consistent scheme on the eigenvalues, ev GW. The latter represents an optimal compromise between computational cost and accuracy.⁴⁵ The other two terms in eqn (3) correspond to the two environmental contributions. Δ^E_n is the electrostatic term, sourced from the potential exerted by the MM subsystem on the QM region, which is captured at the ground-state DFT level. The second environmental contribution, Δ^P_n, is the polarization term, mirroring the dielectric screening of charged excitations by the MM environment.

Neutral singlet (optical) excitations have been obtained with BSE^34–36,46 calculations performed within the Tamm–Dancoff approximation (TDA). The BSE Hamiltonian has been built including occupied/unoccupied states located within 10 eV from the HOMO/LUMO. Earlier work has shown that these settings produce accurate excitation energies for CT excitations.^39,47,48

The dielectric response of the MM environment to charged and neutral excitations within the QM region has been described in the reaction field matrix framework,⁴⁰ explicitly considering an infinite dielectric medium described at the ME level.⁴¹ This approach assumes a frequency-independent dielectric response in the optical region, an approximation that has been validated in previous studies.^41,42 Many-body GW and BSE calculations have been performed with the Fiesta package,^{40–42,46,49,50} employing the Coulomb-fitting resolution-of-the-identity (RI-V) scheme,^51,52 and using the universal Weigend Coulomb fitting set of auxiliary basis functions.⁵³

The MM environment has been described with a ME model of atomistic resolution. The ME model has been parameterized by computing the partial atomic charges and polarizability tensors of its molecular constituents at the DFT (B3LYP/6-311+G**) level. The polymer chains have been partitioned into different units in order to achieve such a parameterization, treating thiophene, thienothiophene and alkyl chains as individual fragments. We note that H atoms have been made implicit within the MM region, leading to a drastic reduction of the number of atoms and interactions to compute, without significantly compromising the accuracy of the ME model. ME and DFT/ME calculations considered PBTTT films with edge-on alignment, consistent with grazing-incidence wide-angle X-ray scattering (GIWAXS) data.^18,54 These calculations employed periodic boundary conditions in the film plane, including interactions with replica up to a circular cutoff of 50 nm, ensuring the convergence of the electrostatic sums. ME calculations have been performed with the MESCal code.⁵⁵

3 Results

3.1 Electronic structure in the non-overlapping limit

Having defined the two systems under study (see Fig. 1) as well as the investigation methodology, we are now in the position to discuss the electronic properties of F4TCNQ–PBTTT complexes, starting from the photoemission energy levels in the zero overlap limit. In this preliminary step, we will hence determine the energy associated to adding an electron to the dopant and removing an electron from PBTTT, allowing us to focus on the effect of environmental interactions in the absence of intermolecular orbital hybridization, which will be introduced later. The non-overlapping limit is practically realized by performing two independent calculations including either the PBTTT segment or the dopant in the QM region of our hybrid QM/MM calculations.

The evolution of the energy levels of the F4TCNQ–PBTTT α and π complexes from the gas phase (GW_g|DFT_g) to the solid state (GW_e|DFT_e) is shown in Fig. 3. A first observation is that the energies in the gas phase are essentially equal for both systems, leading to a gap between the dopant LUMO and the polymer HOMO of about 2 eV. This is indeed expected since without the environment around, the molecules in the π and α systems reduce to the same chemical species, with small differences in the energy levels arising from the slightly different geometries, following from the independent QM/MM optimizations, and impacting the energy levels by less that 50 meV.

Fig. 3 Environmental contributions to the GW energy levels in 3-BTTT and F4TCNQ treated as non-overlapping units in the π (top) and in the α structure (bottom). The electrostatic contribution affects in a very different way the gap between the two systems, due to the very different electrostatic potential felt by the dopant in the two positions (see also Fig. 4). This reflects in pronounced position-dependence of the impurity electron affinity (LUMO level), which makes the dopant ionization more favorable for dopants placed in the alkyl chains (α structure) with respect to the conjugated region (π structure). The induction term similarly shrinks the gap by approximately 2 eV in both systems.

Following our previous work,⁴² we progressively introduce environmental effects in order to disentangle the electrostatic and induction contributions (see eqn (3)). We hence present first an intermediate result obtained upon considering the embedding in the DFT ground state, GW_g|DFT_e, in which energy levels account only for the electrostatic contribution Δ^E_n. The electrostatic term affects in a strikingly different way the energy levels of the two species in the π and α system. Δ_E can be seen as a measure of the electrostatic potential, weighted by the charge density of the molecular orbital involved in the ionization process. Since the polymer HOMO and HOMO−1 levels share the same spatial region, they similarly probe the same potential and are thus equally shifted in a rigid way in each structure (the intent for looking at the HOMO−1 will become clear below). However, the value of the electrostatic energy shift for the polymer levels is appreciably different for the π (Δ_E = −0.19 eV for the HOMO) and α (Δ_E = −0.35 eV for the HOMO) systems. This can attributed to the different electrostatic interactions with the nearby dopant in the MM region, the latter presenting different position and orientation with respect to the polymer chain in the two structures.⁵⁶

Most importantly, the electrostatic contribution to the dopant LUMO is quite different from that of the host levels, and actually very sensitive to position of the impurity in the polymer structure (see Fig. 3). Indeed, Δ_E for the dopant LUMO is found to be 0.73 eV larger in the π structure than in the α one, ultimately leading to very different gap for the two systems.

Such a disparity can be rationalized on the basis of the electrostatic potential in the crystalline polymer matrix, shown in Fig. 4. The potential, sourced from the quadrupole moments of the conjugated polymer backbone, presents a striped pattern characterized by an alternation of high and low potential regions in the alkyl and conjugated chains, respectively. This provides a visual explanation for the difference in the LUMO level for a dopant in the π and α structures. Note that a truly quantitative comparison between the electrostatic potential map in Fig. 4 and QM/MM calculations in Fig. 3 is hindered by the fact that the map is obtained as the convolution of the atomic potentials with a spherical kernel of radius 5 Å (microscopic probe), which is then averaged along the c crystal axis, perpendicular to the figure plane. As such, the map does not reflect the inhomogeneity below the kernel resolution and along the c axis, underestimating the electrostatic potential difference actually felt by dopants inserted in the two regions.

Fig. 4 Electrostatic potential map of the PBTTT crystal obtained from classical microelectrostatics calculations, illustrating the remarkable difference in potential experienced by dopants placed in the region of the conjugated backbone with respect to the alkyl side chains. The local potential has been obtained using a 5 Å-radius spherical probe and averaged over the c crystal axis, perpendicular to the figure plane. The calculation considers periodic boundary conditions in 2D in the ab crystal plane.

We now turn to the fully-embedded calculation results (GW_g|DFT_e in Fig. 3), also including the polarization term Δ_I (see eqn (3)), describing the dielectric screening of charged excitations by the environment. Similar to what observed in other systems,^41,42 this contribution results in approximately symmetric closure of the gap by about 1 eV per level, finally yielding a host-dopant gap of 0.83 eV in the π system and an almost vanishing gap of 70 meV in the α structure.

We emphasize that the very different gap of the two host-dopant structures is entirely due to the electrostatic energy landscape of the polymer lamellae, which heavily affects the EA of the dopant molecule. While recalling that these results have been obtained upon neglecting the intermolecular overlap between the polymer chains and the dopant, we remark that the electrostatic landscape of the crystalline polymer environment is a genuine feature that will also affect the subsequent calculations where the zero-overlap approximation will be lifted.

Before proceeding with the introduction of intermolecular overlap, we report on the calculation of the screened Coulomb interaction V_eh between an electron on the dopant and a hole on the polymer segment, which, together with the intramolecular relaxation energy,⁵⁷ contributes to the energy barrier to set free the doping-induced charge in the semiconductor. V_eh has been computed at the classical ME level, considering charges completely localized on the molecular fragments according to the zero-overlap assumption underpinning ME, and accounting for the frequency-independent dielectric screening provided by the polarizable environment. Our calculations yield V_eh = −0.59 and −0.88 eV for the α and π structure, respectively. The Coulomb binding is hence substantially larger for the π complex, as expected from the shorter distance between dopant and polymer chain, as compared to the α structure. On this basis, we expect that dopants inserted in the alkyl regions would be more effective in increasing the density of mobile carriers participating to conduction.

3.2 Intermolecular orbital hybridization

The energetics of non-overlapping donor and acceptor units provide important insights on the charge-transfer process that is at heart of the doping mechanism. It is, however, essential to introduce the quantum overlap between polymer and dopant, in order to account for intermolecular CT, which in principle can be either of integer or fractional nature.

Charged and neutral (optical) excitations in our many-body framework are calculated starting from the single-particle (supra)molecular orbital obtained at embedded DFT level, i.e. accounting for the embedding of the MM environment. We recall that in the partially self-consistent ev GW scheme here adopted, many-body corrections are applied only to the eigenvalues, not affecting the orbitals that are de facto obtained at DFT level. As customary for closed-shell systems, DFT calculations for the 3-BTTT–F4TCNQ complexes have been performed within the framework of spin-restricted Kohn–Sham theory (RKS).

The energy levels diagram obtained at RKS DFT/MM level, shown in Fig. 5, display the typical frontier orbital repulsion expected for a charge-transfer system upon introducing the intermolecular overlap. Such an opening of the HOMO–LUMO gap of the complex is the result of an intermolecular orbital hybridization, which leads to a moderate fractional charge transfer in both systems. The latter is here measured by the charge on the dopant molecule that amounts to Q_dop = −0.22 and −0.24 in the α and π structure, respectively.

Fig. 5 Frontier energy levels diagram for the π (top) and the α (bottom) 3-BTTT–F4TCNQ complexes, as obtained from spin-restricted (RKS) DFT/MM calculations. Intermolecular orbital hybridization leads to an opening of the gap with respect to non-overlapping fragments (left- and right-most levels), and to a fractional CT ground-state in both systems (Q_dop ≈ −0.2, see text).

Table 1 reports the charge transfer integrals controlling the intermolecular orbital hybridization, along with the overlap, calculated with the dimer projection method⁵⁸ at the PBEh*/6-311G* level of theory. The transfer integrals J between the frontier occupied levels of the polymer and the dopant LUMO are more than two times larger in the π structure than in the α one, an effect that can be attributed to the closest spatial proximity between the two molecular units in the former with stacked molecules. This trend is even more pronounced for the overlap S, which is two orders or magnitude smaller in the α complex.

Table 1 Charge transfer integrals (J) and overlap matrix elements (S) between the occupied molecular orbitals of the polymer, HOMO (H) and HOMO−1 (H−1), and dopant LUMO (L). The overlap between polymer and dopant orbitals is much smaller when the latter is located among alkyl chains (α structure) than in the conjugated region (π structure)

System	H(3-BTTT):L(F4TCNQ)		H−1(3-BTTT):L(F4TCNQ)
	J (meV)	S	J (meV)	S
π	435	1.7 × 10^−2	265	1.4 × 10^−2
α	191	1.8 × 10^−4	29	2.2 × 10^−4

---

We remark that within the RKS framework adopted so far, intermolecular charge transfer follows directly from the formation of doubly-occupied supramolecular orbitals from the quantum mixing (hybridization) between the occupied orbitals of the donor and the unoccupied ones of the acceptor.

Having characterized the ground state of the 3-BTTT/F4TCNQ complex, our analysis proceeds with quasiparticle excitations. Table 2 compares the GW/MM gap obtained for non-overlapping fragments and for the complex, the latter being 0.1–0.2 eV larger, in line with what discussed above for DFT calculations. We underline the very small gap (<0.2 eV) obtained for the α structure. This value is 0.9 eV smaller than in the π system, a result that is essentially due to the electrostatic landscape of polymer lamellae, as discussed in the previous Section.

Table 2 Photoemission and optical gap (eV units) of the two 3-BTTT/F4TCNQ complexes obtained with embedded GW and BSE calculations. The GW/MM gap increases upon considering intermolecular orbital hybridization in the complex, consistent with DFT/MM results in Fig. 5. A very small photoemission gap is found in the α structure, which then results in a negative-energy singlet optical excitation (S1), according to BSE/MM calculations. This nonphysical result calls into question the reliability of the ground state obtained from spin-restricted (RKS) calculations, see text

System	Non-overlapping	Complex
	GW/MM gap	GW/MM gap	BSE/MM ES1
π	0.83	1.08	0.51
α	0.07	0.18	−0.16

---

We have finally calculated the optical excitation energies with embedded BSE calculations, also reported in Table 2. The lowest-energy singlet excitation, S₁, has a neat intermolecular charge-transfer character in both structures, as can be evinced from the analysis of the BSE wavefunction, shown in the ESI,† Fig. S3. The very small gap for the α structure leads to negative optical excitation energies, an unsettling result that can be rationalized, at least formally, upon considering that the gap is considerably smaller in magnitude that the screened electron–hole interaction.

Such a nonphysical negative-energy exciton is here interpreted as a signal of the instability of the ground state obtained from RKS theory, which has been used as a starting point for the many-body perturbative treatment of excitations. The fact that charge-transfer states are found at energies below the moderately-hybridized ground state suggests that the constraint of assigning the same spatial wavefunction to two electrons with opposite spins is too severe, strongly affecting the degree of charge transfer attainable in the ground state of the complex. We further observe that in the present RKS picture, the transfer of one electron from the donor to the acceptor can be realized only on average, i.e. by equally sharing two electrons in the fully-hybridized complex HOMO, i.e. a supramolecular orbital featuring equal contributions from the frontier orbitals of the two molecular fragments. Transferring two half-electrons between two molecules in a complex is something conceptually and practically very different from transferring a full electron, as it would be described in a rigorous multi-reference picture of intermolecular CT.

3.3 Ground-state charge transfer from spin-unrestricted DFT

The unexpected outcome of many-body perturbation theory applied to the α structure, resulting in a negative-energy excitation, has raised a fundamental issue on the reliability of the ground state obtained with spin-restricted (RKS) DFT/MM calculations. In the following, we attempt to overcome the double-occupancy spin constraint by switching to spin-unrestricted Kohn–Sham (UKS) theory, always accounting for the effect of the MM environment on the ground-state properties. It has been shown that UKS theory can improve the description of the ground state of charge-transfer complexes with respect to a RKS scheme,⁵⁹ while remaining within a cost-effective DFT framework. The use of UKS theory, however, comes at the price of introducing some degree of spin contamination in the ground state.

Coming back to the dopant-polymer systems under investigation, we report in Table 3 the net charge on the dopant molecule Q_dop, calculated from the Löwdin atomic charges, as obtained from different RKS and UKS calculations performed on the α and π structures. Mulliken atomic charges yield Q_dop values differing by less than that 0.03 from the Löwdin estimates for all RKS and UKS calculations in Table 3. For each QM/MM calculation, we further report the relative ground-state energy of the system ΔU = U–U (RKS), in order to assess the relative stability of a given solution with respect to the RKS one. We recall that the energy U, computed according to eqn (1), represents the extensive energy of the entire QM/MM system, including the contribution from MM environment and QM–MM interactions.

Table 3 Ground-state properties of 3-BTTT/F4TCNQ complexes obtained in DFT calculations, which differ by the presence/absence of MM embedding, the treatment of the spin degrees of freedom, and the density guess used to initiate the self-consistent field process. ΔU is the total energy of the QM/MM system, referenced to the energy of the spin-restricted (RKS) calculations of the corresponding, π or α, system. The charge on the dopant Q_dop measures the degree of CT in each system. The most stable ground state in the α complex, obtained with spin-unrestricted (UKS) calculations, features integer CT (α₄). Two different fractional CT ground states (π₁ or π₂, and π₄) of similar energy have been instead found for the π structure. ρ_at and ρ(α₃/π₃) identifies the density guess obtained as a superposition of atomic densities and from the α₃/π₃ calculation, respectively. ΔU ∼ 0 corresponds to energy differences smaller than 1 meV in magnitude, well below the accuracy of the method

Label	Method	Spin	Density guess	ΔU (eV)	Q dop (e)
α1	RKS/MM	Singlet	ρ at	0	−0.22
α2	UKS/MM	Singlet	ρ at	∼0	−0.22
α3	UKS	Triplet	ρ at	—	−0.95
α4	UKS/MM	Singlet	ρ (α3)	−0.58	−0.95
α5	UKS	Singlet	ρ (α3)	—	0.03
π1	RKS/MM	Singlet	ρ at	0	−0.24
π2	UKS/MM	Singlet	ρ at	∼0	−0.24
π3	UKS	Triplet	ρ at	—	−0.90
π4	UKS/MM	Singlet	ρ (π3)	−0.07	−0.68
π5	UKS	Singlet	ρ (π3)	—	−0.22

---

We start by discussing the α system. Table 3 compares the embedded UKS calculation, labelled α₂, with the reference calculation, α₁, which corresponds to the embedded RKS that has been discussed in the previous Section. The reference α₁ is characterized by a moderate charge transfer in the ground state (Q_dop = −0.22) and serves as a reference for the total energy (ΔU = 0). The same calculation performed within the UKS scheme, α₂, results in the same ground state within numerical accuracy.

However, we shall note that the DFT calculations proceed through an iterative self-consistent field (SCF) procedure in which the starting Hamiltonian is formed from an initial guess of the density, usually a superposition of atomic densities. Such a choice might impact the result of the SCF process, driving the system towards an local energy minimum, i.e. a metastable state. This observation drove us to attempt initiating the SCF process from a qualitatively different starting point.

We have hence performed triplet gas-phase calculations (α₃ in Table 3), in order to build a quasi full-CT electron density (Q_dop = −0.95) to be employed as a starting guess. When such an ion-pair density is used to initiate a singlet UKS DFT/MM calculation (α₄), we indeed observe its convergence towards an almost integer CT ground state, characterized by Q_dop = −0.95. Interestingly, we observe that the energy of this state is 0.58 eV lower than the one of the reference RKS DFT/MM calculation α₁, which means that this fully-ionized ground state represents the most stable electronic configuration for the α structure.

Additional insight on this integer-CT ground state, α₄, is obtained from the inspection of the spin-polarized Kohn–Sham orbitals, shown in the lower panel of Fig. 6. The figure reveals a negligible hybridization between the molecular orbitals of the polymer and dopant, resulting from the very small overlap between the wave functions of the two units (see Table 1). The LUMO of the isolated dopant is found to be split between the singly-unoccupied molecular orbital (SUMO) and the singly-occupied occupied orbital SOMO−3, which elucidates the electronic configurations resulting in a full-CT ground state.

Fig. 6 Frontier energy levels diagram for the π (top) and the α (bottom) 3-BTTT–F4TCNQ complexes, as obtained from spin-unrestricted (UKS) DFT/MM calculations. A large factional CT resulting from a strong orbital hybridization is obtained for the π structure (π₄ in Table 3), characterized by dopants in the conjugated backbone. A negligible orbital hybridization and an integer CT is instead obtained for the α complex (α₄ in Table 3) with dopants in the alkyl-chains regions. Note the qualitative difference with respect to the RKS results in Fig. 5, the latter considered unreliable for the α system, see text.

It is also interesting to note that while the α₄ UKS DFT/MM calculation starting from the density of the triplet state, α₃, converges towards an ionized ground state, the same guess leads to a quasi neutral state in a gas-phase UKS calculation (α₅). This testifies once more the crucial role played by the electrostatic landscape of the molecular environment for the ionization of dopant impurities in organic semiconductors.

Pursuing our analysis to the π structure, we see that also in this case UKS or RKS DFT/MM calculations starting from the same atomic-density guess lead to the same ground state with moderate charge transfer, Q_dop = −0.24 (see π₁ and π₂ in Table 3). Again, upon initiating the UKS DFT/MM calculation with a triplet quasi full-CT density guess leads to a qualitatively different ground state with respect to the RKS one. In this case (π₄), the ground state presents a large, but not integer, charge transfer with a net charge on the dopant Q_dop = −0.68. We remark that this solution is only slightly more stable that the RKS one, the 70 meV energy difference being comparable to the accuracy of the method. The sizable molecular orbital hybridization in the π₄ solution can be appreciated in the upper panel of Fig. 6, especially in the SOMO and SOMO−2 orbitals, both featuring contributions from both the dopant and polymer unit.

Before concluding this section, we comment on the eigenvalues squared spin operator 〈S²〉 = 1.16 and 0.78 in the lowest-energy UKS solutions α₄ and π₄, respectively. These values close to unity signal a similar proportion of singlet and triplet electronic configurations, which is not really surprising, nor worrying. Indeed, because of the very small intermolecular overlap, singlet and triplet configurations are generally nearly degenerate in pure intermolecular CT states, determining an unavoidable spurious mixing in the UKS ground state. Such a spin contamination is however not expected to have an appreciable impact on the charge density, nor or the ground state energy, being the two spin configurations nearly isoenergetic and both characterized by a hole on the donor and an electron on the acceptor.

To summarize, both the UKS calculation and the instability at the GW/BSE level starting from a RKS ground state point to an integer CT in the ground-state for the α structure. The environmental embedding is found to be crucial to describe such a CT character. The π complex is instead more incline to partial CT, yet our calculations did not allow us to conclude whether the fractional charge transferred to the dopant is small (Q_dop = −0.24, as obtained with RKS calculations) or large (Q_dop = −0.68, as from UKS).

4 Discussion and conclusions

In this paper, we have investigated the relationships between nanostructure and electronic properties in the prototypical crystalline polymer PBTTT doped by the strong electron acceptor F4TCNQ. By applying theoretical tools such as density functional theory and many-body perturbation theory combined with classical microelectrostatics, we have elucidated the differences in the electronic structure of host–dopant complexes characterized by impurity molecules placed either in the conjugated backbone (π complex) or in between the alkyl chains (α complex) of PBTTT crystalline lamellae.

On the methodological side, we remark two crucial, though often-overlooked, factors for the proper description of host–dopant complexes. First, the importance of accounting for the electrostatic fields and the dielectric screening provided by the molecular environment, strongly affecting the energetics of the charge-transfer phenomena at the core of doping. Second, we highlight the difficulty of routine spin-restricted calculations in describing integer charge transfer states, with most of the calculations on host–dopant complexes found in the literature reporting fractional ionization.^8,17,60,61 Spin-unrestricted calculations represent a convenient framework for describing the crossover between neutral and ionized ground states, without resorting to constrained formalisms.¹⁴

Our analysis shows that dopants inserting into the alkyl chains can undergo an integer charge transfer, while impurities placed in the conjugated region are instead more prone to fractional ionization. The different propensity for an integer or fractional ground-state charge transfer in the two systems can be ascribed to two different factors. First, the electrostatic landscape within the polymer lamellae, which results in a strong position dependence of the LUMO level (electron affinity), making F4TCNQ a much stronger acceptor when it is placed in the alkyl chains than in the conjugated region. We emphasize that this is a major effect, with a difference in electron affinity exceeding 0.5 eV, comparable to the variations for a given impurity in different organic hosts calculated by some of us in a previous study.¹³ Second, the larger overlap between the host and polymer wave functions favors a stronger intermolecular hybridization in the π complex with respect to the α one.

These results, together with the considerably smaller (∼30%) Coulomb binding between electron and hole in the α complex, point towards the structures characterized by dopants inserted in the alkyl side chains of polymer lamellae being the most favorable for high doping efficiency and hence large electrical conductivities. Our findings rationalize a series of experimental observations, such as the fact that, to the best of our knowledge, fractional charge transfer has never been observed in semi-crystalline polymers where it has been ascertained that dopants are inserted in the alkyl side chains.^17,24 Most importantly, conductivities exceeding 100 S cm⁻¹ have been reported only for systems with dopants intercalated in the alkyl chains,^{18,21,22,26,27,62} an observation that is in line with the general picture emerging from our electronic structure calculations. The present study, by illustrating the key factors controlling the doping efficiency as a function of the dopant position in a polymer matrix, brings forward our understanding of the often elusive relationship between structure and electronic properties in doped conjugated polymers.

Author contributions

M. C. performed electronic structure calculations on doped polymers, using morphologies obtained by V. L. with force-field simulations; A. G. developed and validated our approach to the calculation of total energies in a DFT/MM framework; G. D. supervised the work and drafted a first version of the manuscript. All authors discussed the results and contributed to the redaction of the manuscript in its final form.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank Ivan Duchemin and Jing Li for providing codes implementing many-body Green's function techniques and tools for QM/MM calculations. M. C. acknowledges PhD scholarship from Grenoble Quantum Engineering (GreQuE) program, funded by Fondation Nanosciences and the European Unions Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 754303. The work in Grenoble has been supported by the French “Agence Nationale de la Recherche” through the project RAPTORS (ANR-21-CE24-0004-01). High-performance computing resources from GENCI-TGCC (Grant no. 2020-A0090910016) are acknowledged. The work in Mons has received funding from the European Unions Horizon 2020 research and innovation programme under grant agreement No 964677 (MITICS). The molecular modeling activities in Mons are supported by FNRS (Consortium des Équipements de Calcul Intensif - CÉCI, under Grant 2.5020.11) and by the Walloon Region (ZENOBE Tier-1 supercomputer, under grant 1117545). DB is a FNRS research director.

References

1. K. Walzer, B. Maennig, M. Pfeiffer and K. Leo, Chem. Rev., 2007, 107, 1233–1271.
2. I. Salzmann, G. Heimel, M. Oehzelt, S. Winkler and N. Koch, Acc. Chem. Res., 2016, 49, 370–378.
3. I. E. Jacobs and A. J. Moulé, Adv. Mater., 2017, 29, 1703063.
4. M. L. Tietze, J. Benduhn, P. Pahner, B. Nell, M. Schwarze, H. Kleemann, M. Krammer, K. Zojer, K. Vandewal and K. Leo, Nat. Commun., 2018, 9, 1182.
5. A. Mityashin, Y. Olivier, T. Van Regemorter, C. Rolin, S. Verlaak, N. G. Martinelli, D. Beljonne, J. Cornil, J. Genoe and P. Heremans, Adv. Mater., 2012, 24, 1535–1539.
6. A. Fediai, F. Symalla, P. Friederich and W. Wenzel, Nat. Commun., 2019, 10, 4547.
7. M. Comin, S. Fratini, X. Blase and G. D’Avino, Adv. Mater., 2022, 34, 2105376.
8. H. Mendez, G. Heimel, S. Winkler, J. Frisch, A. Opitz, K. Sauer, B. Wegner, M. Oehzelt, C. Rothel, S. Duhm, D. Tobbens, N. Koch and I. Salzmann, Nat. Commun., 2015, 6, 8560.
9. J. Li, G. D’Avino, A. Pershin, D. Jacquemin, I. Duchemin, D. Beljonne and X. Blase, Phys. Rev. Mater., 2017, 1, 025602.
10. C. P. Theurer, M. Richter, D. Rana, G. Duva, D. Lepple, A. Hinderhofer, F. Schreiber, P. Tegeder and K. Broch, J. Phys. Chem. C, 2021, 125, 23952–23959.
11. B. Maennig, M. Pfeiffer, A. Nollau, X. Zhou, K. Leo and P. Simon, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 195208.
12. G. D’Avino, L. Muccioli, F. Castet, C. Poelking, D. Andrienko, Z. G. Soos, J. Cornil and D. Beljonne, J. Phys.: Condens. Matter, 2016, 28, 433002.
13. J. Li, I. Duchemin, O. Roscioni, P. Friederich, M. Anderson, E. Da Como, G. Kociok-Köhn, W. Wenzel, C. Zannoni, D. Beljonne, X. Blase and G. D’Avino, Mater. Horiz., 2019, 6, 107.
14. C. Gaul, S. Hutsch, M. Schwarze, K. S. Schellhammer, F. Bussolotti, S. Kera, G. Cuniberti, K. Leo and F. Ortmann, Nat. Mater., 2018, 17, 439.
15. D. T. Duong, C. Wang, E. Antono, M. F. Toney and A. Salleo, Org. Electron., 2013, 14, 1330–1336.
16. J. E. Cochran, M. J.-N. Junk, A. M. Glaudell, P. L. Miller, J. S. Cowart, M. F. Toney, C. J. Hawker, B. F. Chmelka and M. L. Chabinyc, Macromolecules, 2019, 47, 6836–6846.
17. I. E. Jacobs, C. Cendra, T. F. Harrelson, Z. I. Bedolla Valdez, R. Faller, A. Salleo and A. J. Moulé, Mater. Horiz., 2018, 5, 655.
18. K. Kang, S. Watanabe, K. Broch, A. Sepe, A. Brown, I. Nasrallah, M. Nikolka, Z. Fei, M. Heeney and D. Matsumoto, Nat. Mater., 2016, 15, 896–902.
19. S. N. Patel, A. M. Glaudell, K. A. Peterson, E. M. Thomas, K. A. OHara, E. Lim and M. L. Chabinyc, Sci. Adv., 2017, 3, e1700434.
20. D. T. Scholes, P. Y. Yee, J. R. Lindemuth, H. Kang, J. Onorato, R. Ghosh, C. K. Luscombe, F. C. Spano, S. H. Tolbert and B. J. Schwartz, Adv. Funct. Mater., 2017, 27, 1702654.
21. V. Vijayakumar, P. Durand, H. Zeng, V. Untilova, L. Herrmann, P. Algayer, N. Leclerc and M. Brinkmann, J. Mater. Chem. C, 2020, 8, 16470.
22. V. Untilova, H. Zeng, P. Durand, L. Herrmann, N. Leclerc and M. Brinkmann, Macromolecules, 2021, 54, 6073–6084.
23. H. Hase, K. ONeill, J. Frisch, A. Opitz, N. Koch and I. Salzmann, J. Phys. Chem. C, 2018, 122, 25893–25899.
24. D. A. Stanfield, Y. Wu, S. H. Tolbert and B. J. Schwartz, Chem. Mater., 2021, 33, 2343–2356.
25. Y. Yamashita, J. Tsurumi, M. Ohno, R. Fujimoto, S. Kumagai, T. Kurosawa, T. Okamoto, J. Takeya and S. Watanabe, Nature, 2019, 572, 634–638.
26. I. E. Jacobs, Y. Lin, Y. Huang, X. Ren, D. Simatos, C. Chen, D. Tjhe, M. Statz, L. Lai, P. A. Finn, W. G. Neal, G. D’Avino, V. Lemaur, S. Fratini, D. Beljonne, J. Strzalka, C. B. Nielsen, S. Barlow, S. R. Marder, I. McCulloch and H. Sirringhaus, Adv. Mater., 2021, 2102988.
27. V. Vijayakumar, Y. Zhong, V. Untilova, M. Bahri, L. Herrmann, L. Biniek, N. Leclerc and M. Brinkmann, Adv. Energy Mater., 2019, 9, 1900266.
28. D. Venkateshvaran, M. Nikolka, A. Sadhanala, V. Lemaur, M. Zelazny, M. Kepa, M. Hurhangee, A. J. Kronemeijer, V. Pecunia, I. Nasrallah, I. Romanov, K. Broch, I. McCulloch, D. Emin, Y. Olivier, J. Cornil, D. Beljonne and H. Sirringhaus, Nature, 2014, 515, 384–388.
29. V. Lemaur, J. Cornil, R. Lazzaroni, H. Sirringhaus, D. Beljonne and Y. Olivier, Chem. Mater., 2019, 31, 6889–6899.
30. L. Hedin, Phys. Rev., 1965, 139, A796–A823.
31. G. Strinati, H. J. Mattausch and W. Hanke, Phys. Rev. Lett., 1980, 45, 290–294.
32. M. S. Hybertsen and S. G. Louie, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 5390–5413.
33. G. Onida, L. Reining and A. Rubio, Rev. Mod. Phys., 2002, 74, 601–659.
34. E. E. Salpeter and H. A. Bethe, Phys. Rev., 1951, 84, 1232–1242.
35. L. J. Sham and T. M. Rice, Phys. Rev., 1966, 144, 708–714.
36. W. Hanke and L. J. Sham, Phys. Rev. Lett., 1979, 43, 387–390.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
38. S. Grimme, WIREs Comput. Mol. Sci., 2011, 1, 211–228.
39. B. Baumeier, M. Rohlfing and D. Andrienko, J. Chem. Theory Comput., 2014, 10, 3104–3110.
40. I. Duchemin, D. Jacquemin and X. Blase, J. Chem. Phys., 2016, 144, 164106.
41. J. Li, G. DAvino, I. Duchemin, D. Beljonne and X. Blase, J. Phys. Chem. Lett., 2016, 7, 2814–2820.
42. J. Li, G. D’Avino, I. Duchemin, D. Beljonne and X. Blase, Phys. Rev. B, 2018, 97, 035108.
43. H. M. Senn and W. Thiel, Angew. Chem., Int. Ed., 2009, 48, 1198–1229.
44. F. Neese, WIREs Comput. Mol. Sci., 2012, 2, 73–78.
45. X. Blase, C. Attaccalite and V. Olevano, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 115103.
46. X. Blase, I. Duchemin and D. Jacquemin, Chem. Soc. Rev., 2018, 47, 1022–1043.
47. X. Blase and C. Attaccalite, Appl. Phys. Lett., 2011, 99, 171909.
48. I. Duchemin, T. Deutsch and X. Blase, Phys. Rev. Lett., 2012, 109, 167801.
49. C. Faber, P. Boulanger, C. Attaccalite, I. Duchemin and X. Blase, Philos. Trans. R. Soc., A, 2014, 372, 20130271.
50. D. Jacquemin, I. Duchemin and X. Blase, J. Chem. Theory Comput., 2015, 11, 3290–3304.
51. J. L. Whitten, J. Chem. Phys., 1973, 58, 4496–4501.
52. I. Duchemin, J. Li and X. Blase, J. Chem. Theory Comput., 2017, 13, 1199–1208.
53. F. Weigend, Phys. Chem. Chem. Phys., 2006, 8, 1057–1065.
54. H. Tanaka, A. Wakamatsu, M. Kondo, S. Kawamura, S.-I. Kuroda, Y. Shimoi, W.-T. Park, Y.-Y. Noh and T. Takenobu, Commun. Phys., 2019, 2, 1–10.
55. G. D’Avino, L. Muccioli, C. Zannoni, D. Beljonne and Z. Soos, J. Chem. Theory Comput., 2014, 10, 4959–4971.
56. A. Privitera, G. Londi, M. Riede, G. D’Avino and D. Beljonne, Adv. Funct. Mater., 2020, 30, 2004600.
57. M. Schwarze, C. Gaul, R. Scholz, F. Bussolotti, A. Hofacker, K. S. Schellhammer, B. Nell, B. D. Naab, Z. Bao, D. Spoltore, K. Vandewal, J. Widmer, S. Kera, N. Ueno, F. Ortmann and K. Leo, Nat. Mater., 2019, 18, 242.
58. E. F. Valeev, V. Coropceanu, D. A. da Silva Filho, S. Salman and J.-L. Brédas, J. Am. Chem. Soc., 2006, 128, 9882–9886.
59. V. Geskin, R. Stadler and J. Cornil, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 085411.
60. I. Salzmann, G. Heimel, S. Duhm, M. Oehzelt, P. Pingel, B. M. George, A. Schnegg, K. Lips, R.-P. Blum, A. Vollmer and N. Koch, Phys. Rev. Lett., 2012, 108, 035502.
61. A. M. Valencia and C. Cocchi, J. Phys. Chem. C, 2019, 123, 9617–9623.
62. D. Kiefer, R. Kroon, A. I. Hofmann, H. Sun, X. Liu, A. Giovannitti, D. Stegerer, A. Cano, J. Hynynen, L. Yu, Y. Zhang, D. Nai, T. F. Harrelson, M. Sommer, A. J. Moulé, M. Kemerink, S. R. Marder, I. McCulloch, M. Fahlman, S. Fabiano and C. Müller, Nat. Mater., 2019, 18, 149–155.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2tc01115h

---