# Exciton Radiative Lifetime and Diffusion Length from First-Principles Atomic Trajectories: Method Assessment for P3HT

Cheick Oumar Diarra, Julien Taillieu, Mauro Boero, Thomas Heiser, Evelyne Martin

► To cite this version:

Cheick Oumar Diarra, Julien Taillieu, Mauro Boero, Thomas Heiser, Evelyne Martin. Exciton Radiative Lifetime and Diffusion Length from First-Principles Atomic Trajectories: Method Assessment for P3HT. Journal of Physical Chemistry B, 2025, 10.1021/acs.jpcb.5c04516 . hal-05227333

HAL Id: hal-05227333

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

Submitted on 28 Aug 2025

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Exciton radiative lifetime and diffusion length
from first-principles atomic trajectories:
Method assessment for P3HT

Cheick Oumar Diarra, Julien Taillieu, Mauro Boero, Thomas Heiser, and Evelyne

Martin*

University of Strasbourg, CNRS, ICube Laboratory UMR 7357, 300 Boulevard Sébastien Brant,
F-67412 Illkirch, France

E-mail: evelyne.martin@cnrs.fr

Abstract

Organic photovoltaic (OPV) technology offers a solution to the global demand for sus-
tainable energy production and harvesting due to its affordability, versatility, and moderate
environmental impact. Recent advancements have led to increased efficiency in bulk hetero-
junction (BHJ) solar cells, making practical applications feasible. However, these advance-
ments have also raised fundamental questions regarding exciton transport mechanisms. A
precise understanding of the relationship between the exciton diffusion length and the mate-
rial structure is crucial for improvements and applications. We formerly assessed a method to
extract the exciton diffusion coefficient [Diarra et al., Phys. Chem. Chem. Phys., **2023**, 25,
15539-15546] via first principles dynamical simulations in the excited state. Here, we propose
a novel approach to compute the radiative lifetime within the same method. This provided an
additional crucial parameter contributing to the exciton diffusion length. While the method

<footnote>*To whom correspondence should be addressed</footnote>
is limited to radiative contribution, we benchmark its validity on a specific semiconducting
polymer poly(3-hexylthiophene) (P3HT). The resulting exciton diffusion length is found to be
in close agreement with available experimental data, which might indicate that the radiative
process is dominant with respect to non-radiative contributions.

INTRODUCTION

The global demand for renewable energies with high efficiency, low production costs, and minimal environmental impact has propelled organic solar cells (OSCs) as a promising solution. These materials are inherently environmentally friendly due to their organic nature, and their relatively low manufacturing costs and exceptional flexibility make them attractive for various applications beyond photovoltaics. Recent breakthroughs¹⁻³ have overcome major obstacles in practical OSC use, with significant power conversion efficiencies (PCEs) now attaining efficiency values of 20.82% in bulk heterojunctions⁴ (BHJs). However, these advancements rely heavily on empirical verifications and repeated experimental efforts. To further enhance OSC efficiency, a deeper understanding of the underlying mechanisms at the microscopic level is crucial, along with reliable methods to extract key features.

Organic semiconductors generate bound electron-hole pairs (excitons) upon photon absorption. A critical parameter influencing power conversion efficiency is the exciton diffusion length (L), which determines the optimal dimensions of donor and acceptor domains in both BHJs and planar heterojunctions (PHJs). The diffusion length is given by the diffusion coefficient (D) and the intrinsic radiative lifetime (τ) of excitons, related by the equation $L = \sqrt{D\tau}$ for the case of one-dimensional diffusion as the one on which this article is focused and realized by our targeted system. If the domains where the exciton is generated exceed the critical length L, the exciton recombines and is unable to contribute to energy generation. The typical diffusion length L is equal to 10-20 nm, which means that BHJs should be preferred, despite instability issues over time due to the rearrangement of the two components in the mixture. The incorporation of bilayers in PHJs could potentially lift this problem. Nevertheless, the absorption of solar light by such thin
layers would be insufficient. Consequently, the competitiveness of the OPV sector hinges on the discovery of organic semiconductors with extended diffusion lengths. Our proposal entails the use of first-principles molecular dynamics (FPMD) simulations to unravel the behavior of materials in general, and of organic compounds in particular, at finite temperature at the atomic scale and with quantitative predictivity.

In the field of organic semiconductors, several theoretical approaches have been proposed over the years$^{5-8}$ to tackle the problem of excitons in organic semiconductors and to describe excitations from the ground state to the first excited state. These approaches, allowing for non-adiabatic transitions (e. g. surface hopping schemes), have been able to assess the validity of the Frenkel exciton theory and to probe possible contribution from lower-lying excited states beyond the first excitation. Often, to reduce the high computational cost of propagating dynamically the ground and excited states, electronic structure calculations are done on the fly during a first-principles dynamics but rather done on an accompanying classical molecular dynamics.$^{8}$

Recently, we introduced a methodology for computing the diffusion coefficient ($D$) by resorting to FPMD, showcasing its predictive capabilities with specific application to poly(3-hexylthiophene) (P3HT).$^{9}$ P3HT, a semiconducting polymer conventionally used as an electron donor in photo-voltaic cells, has been extensively studied for its structural and excitonic properties, making it an ideal material for evaluating the effectiveness of our approach. While the diffusion coefficient $D$ is more challenging to extract experimentally, the lifetime $\tau$ typically obtained by time-resolved photoluminescence (PL) may vary depending on the preparation conditions,$^{10,11}$ and several lifetimes may emerge from the PL decay due to the presence of molecular defects or impurities or the coexistence of crystalline and amorphous domains.$^{12}$ It is therefore challenging to establish a correlation between a controlled structural configuration and an exciton diffusion length. In this study, we compute the radiative exciton lifetime ($\tau_{rad}$) of P3HT using a FPMD-based approach on a well-established atomic model, the same used to calculate $D$. This eliminates uncertainties about the structure and provides a framework to connect a given and well-known structure of the material. We stress the fact that we are not aiming at providing a whole paradigm to include non-adiabatic
transitions. In fact, our approach makes use of adiabatic dynamics on the ground and excited state.
This will be explored in a forthcoming work. Our purpose is to check to which extent radiative
transitions contribute to the overall lifetime $\tau$, within the known relationship $\tau^{-1} = \tau_{rad}^{-1} + \tau_{nonrad}^{-1}$.
The case of P3HT, targeted in the present work and extensively studied in the literature$^{13}$ repre-
sents an exceptionally versatile reference system. In particular, the crystalline domains have been
accurately characterized$^{14-16}$ to unravel their specific role. The diffusion coefficient $D$ in these
crystalline domains could be extracted with great accuracy, and the associated lifetime $\tau$ was mea-
sured at the same high accuracy level, thus allowing to extract a reliable value of the diffusion
length $L$. This system is therefore an ideal template for benchmarking the methodology proposed
here.

METHODS

The P3HT model use in the present study is the one constructed and benchmarked in our previous
work.⁹ It consists of two chains containing six thiophene units each within a periodic box. The
box dimensions are adjusted to mimic an infinite-length polymer crystal, with full hexyl chains
included. This configuration results in a structure composed of 300 atoms, including 12 sulfur
(S), 120 carbon (C), and 168 hydrogen (H) atoms, in a simulation cell of dimensions $L_x = 24.0$
Å, $L_y = 16.0$ Å, and $L_z = 7.6$ Å. Here, the x direction aligns with the polymer axis, the y direction
accommodates the alkyl chains, regularly terminated by H atoms and not carrying any covalent
bonds across the edges of the periodically repeated simulation cell, and the z direction represents
the stacking direction. We have already shown elsewhere⁹ that this model is able to mimic the
regio-regular arrangement of the P3HT polymer chains and the periodicity along the z direction is
the correct one to simulate a periodically repeated system, as opposed to a simple dimer. This com-
putational set-up yields a structure closely ressembling the experimental one denoted as "Form I"
at room temperature ($T = 300$ K). The $\pi-\pi$ distance ($d_{\pi-\pi}$) measures 3.5 Å, the stacking distance
($d_{stack}$) is 3.7 Å, and the tilt angle ($\alpha$) is 23° (see Fig. 1), in good agreement with experimental
<figure><img src="image_2.png" /><figcaption>Figure 1: Model of the P3HT crystal at T = 300 K. **a.** view of the model along the *x* axis; the stacking and $\pi-\pi$ distances are indicated by the black and red lines, respectively, with numerical data provided in the text. **b.** z-axis view of the same polymers and their alkyl chains. In both view, the blue dashed line represents the unit cell.</figcaption></figure>

FPMD simulations were performed within the Car-Parrinello method¹⁹ frameworks implemented in the developers version 4.3 of the CPMD code.²⁰ The functional used to describe the exchange-correlation (XC) is the one proposed by Perdew, Burke and Ernzerhof (PBE).²¹ We are aware of both the inherent gap underestimation of the PBE functional and its tendency to give rise to more delocalized orbitals. Corrections are then made to cope with these shortcomings as specified in the ongoing discussion. Since long-range van der Waals (vdW) dispersion forces are not included in any XC functional, we complemented the PBE by adding the Grimme-D3 vdW semiempirical correction²²,²³ to cope with the interactions regulating the polymer chains stacking. Only valence electrons are treated explicitly, while core-valence interactions were described by norm-conserving Trouillier-Martins pseudopotentials.²⁴ Valence electron wavefunctions, in a spin unrestricted approach, were expanded on a plane wave basis set with a cutoff energy of 70 Ry (952.4 eV) and the sampling of the Brillouin zone was restricted to the $\Gamma$ point. To ensure good numerical control of the constants of motion, an integration time step of 0.01 fs (4.0 a.u.) was chosen for ion dynamics, and a fictitious mass of 400 a.u. was used to propagate the electronic degrees of freedom. The ionic temperature was controlled with a Nosé-Hoover²⁵⁻²⁷ thermostat chain²⁸ with a frequency of 2500 cm⁻¹. To preserve the adiabaticity in the Car-Parrinello approach, the fictitious kinetic energy of the electronic orbitals was controlled with a Blöchl-Parrinello thermo-
stat²⁹ with a frequency of 5000 cm⁻¹. The singlet excited state $S_1$ was accounted for by resorting to the restricted open-shell Kohn-Sham (ROKS) formalism.³⁰⁻³² We remark that this approach, being an extension of DFT-based FPMD, contains all interactions among electrons and ions, thus accounting for any harmonic and anharmonic effect by construction, hence not relying on any ad-hoc harmonic approximation hypothesis and implicitly accounting for the vibronic fine structure typical of the systems targeted.³³,³⁴ The accuracy of ROKS in comparison with time-dependent density functional theory (TDDFT) approaches, at least for small molecules, has been investigated in a series of former works on different benchmark systems³³⁻³⁶ showing its computational reliability.

To understand the relationship between the structural characteristics of the organic material and the values of the exciton diffusion length $L$, it is essential to have a methodology allowing to extract simultaneously both $D$ and $\tau$, along with the requirement of being able to survey the underlying dynamical behavior at an atomic-scale level for a rationalization of the quantities extracted. Although we do not have, at this stage, a way to obtain the full $\tau$, we have a practical algorithm to compute the radiative contribution $\tau_{rad}$ and to attempt at evaluating the contribution that this would imply on the diffusion length. The method proposed here offers then a practical way to get $D$ and $\tau_{rad}$ from first principles molecular dynamics trajectories, avoiding computationally expensive workload of TDDFT or ad-hoc phenomenological approaches. As an initial step, we started by relaxing the polymer structure to its local minimum in the ground state $S_0$. Subsequently, the system is transitioned to the $S_1$ state, and a secondary relaxation was conducted within the ROKS framework. Then, the temperature was gradually raised, beginning with a 5 ps interval at $T = 100$ K before reaching 300 K. Since spontaneous emission occurs from the singlet excited state $S_1$ and that radiative transitions are vertical, we can safely assume that the emission spectrum is mainly dominated by the dynamics of the ions in $S_1$ (Fig. 2). Consequently, a 40 ps trajectory is generated in the $S_1$ state, from which a finite number $N$ of temporally equidistant (uncorrelated) configurations are sampled. We stress that within the ROKS-FPMD formalism, all interactions intra- and inter-chain, their coupling to the exciton and to the atomic motion responsible for the vibrational
spectrum are implicitly accounted for, along with the related energy fluctuations observable during the finite temperature run.

Auxiliary calculations on a tetramer system consisting of 600 atoms in a periodically repeated simulation cell of sizes 24.0 × 16.0 × 15.2 3 were done to check (i) for possible size effects in the z stacking and to inspect the role of the addition of exact exchange within a hybrid functional framework in terms of (ii) wavefunctions localization, (iii) energy separation between the ground and excited states and (iv) energy gap estimation. For this auxiliary simulation, a dynamical trajectory was generated at the PBE level for 5.86 ps and on this trajectory a configuration was sampled for refined calculations with a hybrid HSE06 37,38 functional in which the amount of exact exchange included was 25%. The specific configuration selected corresponds to the closest $E_{S_0}$ and $E_{S_1}$ energy separation at the PBE level. These results are summarized in the Supporting Information (SI, Figs. SI.1, SI.2, SI.3 and Table SI.1). In practice, a hybrid functional approach has the net effect of giving slightly less dispersed wavefunctions, as expected for electronic states treated at this level, 39 although preserving the overall PBE electron distribution of the main lobes. Moreover, they confirm that a diffusion of the excited electron occurs between adjacent polymer chains without involving the whole system and its periodicity induced by the computational set-up On an energetic standpoint, the $E_{S_0}$ and $E_{S_1}$ energy separation increases in the case of HSE06 estimations, thus making the two states $S_0$ and $S_1$ adiabatically well separated.

To characterize the optical transitions, we followed an approach assessed in former works 40,41 to compute the PL(E) spectrum versus transition energy E which includes the calculation of the transition moments $|\langle\psi_i|\mathbf{p}|\psi_j\rangle|^2$ between initial states i and final ones j:

$$ \mathrm{PL}(E) = \mathrm{constant} \cdot \frac{1}{E} \sum_{i,j} (f_i - f_j) |\langle\psi_i|\mathbf{p}|\psi_j\rangle|^2 \delta(E_i - E_j - E) \quad (1) $$

In Eq. 1, $\psi_i$ and $E_i$ are the wavefunctions and eigenvalues, $f_i = 0, 1$ are the occupation numbers and $\mathbf{p}$ the (quantum) momentum operator. The nature and localization of the excited state has already been reported in our former work. 9 Namely, both the electron and the hole are initially localized on
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration depicting the method employed to compute the emission spectrum and exci- ton lifetime. The red line represents the potential energy surface of the S$_{1}$ state, while the black line corresponds to the ground state S$_{0}$. The circular line with an arrowhead symbolizes the trajectory within the S$_{1}$ state. Black straight arrows mark the S$_{1}$ → S$_{0}$ transition at specific points along the trajectory, where configurations are sampled and transition matrices are calculated.</figcaption></figure>

the same polymer chain before evolving dynamically in a situation in which they separate and one of the two diffuses on the nearby polymer. Here the initial state is an atomic configuration of the S$_{1}$ trajectory, and the final state is the S$_{0}$ state of the same atomic configuration. In the framework of ROKS, the S$_{1}$ state is constructed as a linear combination of the antiparallel spin states m$_{1}$ and m$_{2}$ (Fig. 3), namely

$$
|S_1\rangle = \frac{1}{\sqrt{2}} (|m_1\rangle + |m_2\rangle). \qquad (2)
$$

Consequently, the matrix elements of the S$_{1}$ → S$_{0}$ transition in Eq. 1 can be expressed in terms

<figure><img src="image_4.png" /><figcaption>Figure 3: Spin occupation of the highest occupied (HOMO) and lowest unoccupied (LUMO) molecular orbitals in the m₁ (left) and m₂ (right) states.</figcaption></figure>

of the matrix elements of the $m_1 \rightarrow S_0$ and $m_2 \rightarrow S_0$ transitions. In the absence of a magnetic field, the $m_1 \rightarrow S_0$ and $m_2 \rightarrow S_0$ recombinations are symmetrical and $|\langle m_1 | \mathbf{r} | S_0 \rangle| = |\langle m_2 | \mathbf{r} | S_0 \rangle|$
and is noted $|\langle m|\mathbf{r}|S_0\rangle|$. By keeping into account all these factors, the rate $A_{S_1 \to S_0}^n$ characterizing the $S_1 \to S_0$ transition via spontaneous emission can be expressed as (writing the proportionality constant explicitly):

$$A_{S_1 \to S_0}^n = \frac{2e^2 E_n^3}{3\pi\epsilon_0\hbar^4 c^3} |\langle m|\mathbf{r}|S_0\rangle_n|^2 \quad (3)$$

for each sampled configuration $n$, where $\mathbf{r}$ is the (formal) position operator, actually rewritten in a gauge invariant form as in Eq. 1 by exploiting the canonical commutation relation $\mathbf{p}/m = (i/\hbar) [H, \mathbf{x}]$ being $H$ the Hamiltonian operator and $\hbar$ the reduced Planck's constant. The remaining parameters are $e$ is the electron charge, $c$ the speed of light in a vacuum and $\epsilon_0$ the vacuum permittivity. The radiative exciton lifetime $\tau_{rad}$, within our approach, is computed from the average of the transition rates over the dynamically generated configurations given by

$$\tau_{rad} = 1 / \left( \frac{1}{N} \sum_{n=1}^{N} A_{S_1 \to S_0}^{n} \right) \quad (4)$$

where the sum is on the $N$ sampled configurations, and $A_{S_1 \to S_0}^n$ is the transition rate in the $n^{th}$ configuration.

Transition moments and their associated energies are the quantities that determine the PL emission spectra. To be specific, when illustrating the emission spectra, the Dirac delta function in Eq. 1 is substituted by a Gaussian function. The contributions of these Gaussians, centered on the transition energy of the configuration, are then summed. The amplitude of each Gaussian, centered on a transition energy $E_n$, is weighted by the oscillator strength $f_n^{osc}$ given by:

$$f_n^{osc} = \frac{2 m_e}{3 \hbar^2} E_n |\langle m | \mathbf{r} | S_0 \rangle_n|^2. \quad (5)$$

$(m_e$ is the electron mass). This can be regarded as the probability of a direct transition between states $S_1$ and $S_0$. The Gaussian superposition becomes then

$$g(E) = \sum_{n=1}^{N} \frac{f_n^{osc}}{\sigma \sqrt{2\pi}} \exp \left( - \frac{(E - E_n)^2}{2\sigma^2} \right) \quad (6)$$
where the Gaussian standard deviation $\sigma$ is consistently set to $k_B T$.

Although most of the time the two surfaces $S_0$ and $S_1$ are energetically well separated, justifying the interpretation of $\tau_{rad}$ in Eq. 4 as essentially the radiative lifetime, we observe instantaneous events in which the two total energies $E_{S_0}$ and $E_{S_1}$ come close to each other by less than 0.5 eV, within one time step of the numerical integration of the equations of motion. It is difficult to interpret such a relatively small energetic difference as a trigger of (rare) non-radiative transitions between the two states.$^{42,43}$ This is further supported by the increased value that a hybrid functional would provide as shown in Table SI.1 of the SI. Although our approach does not account for non-adiabatic transitions, at least on the time scale of the simulations it seems that non-radiative events should play a marginal role for the targeted P3HT system. This is supported by two observations. The first one concerns the electronic structure of the system. We checked all along the trajectory that no sudden change of the electronic structure occurs, a fact that could jeopardize the Car-Parrinello-like propagation of the orbitals. The second observation is summarized in Fig. 4 where a non-zero separation between $E_{S_0}$ and $E_{S_1}$ is kept at each moment of the simulation.

<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of the two total energies $E_{S_0}$ and $E_{S_1}$ along the dynamics and their relative difference as a function of the simulation time.</figcaption></figure>
RESULTS AND DISCUSSION

Following the methodology outlined in the former paragraph, we generated the emission spectrum
by sampling initial configurations from the excited state S<sub>1</sub> along a dynamic trajectory (Fig. 2). The
vertical transition to the S<sub>0</sub> state at each of these configuration is characterized by computing the
corresponding recombination rate using transition moments. The resulting spectra, depicted in Fig.
5, include three different configuration counts: N =45, 100, and 200. These specific choices are
motivated by the necessity for uncorrelated configurations and the evaluation we aim to conduct
regarding the sampling approach. Therefore, given our 40 ps dynamic simulation duration, we
confidently sampled the three configuration sets mentioned earlier, corresponding to time intervals
of 0.88, 0.40, and 0.20 ps along the trajectory. In each of the three samplings, the transition matrix
element is computed all along the generated dynamical trajectory at each sampled instant, where
the system is in the atomic configuration resulting from the ROKS-FPMD evolving at 300 K, thus
subject to its dynamical evolution. Fig. 5 shows that while the overall spectral width is nearly
constant (approximately 0.6 eV) and independent of the sampling, the features and peak positions
of the different spectra reflect a dependence on the sampling.

It's worth noting that any DFT approach requires an energetic correction due to the inher-
ent underestimation of energy gaps.⁴⁶ This challenge can be addressed by refining computational
methods with techniques better suited for gap evaluation, such as hybrid functionals, the GW ap-
proximation, or higher-level quantum chemical methods. As far as the hybrid functionals are con-
cerned, we wish to bring to the attention the fact that the amount of exact exchange included has to
be carefully tuned, as shown in Table SI.1 of the SI, where a 25% of Hartree-Fock exchange results
in an overestimation of the HOMO-LUMO gap of about 0.1 eV with respect to the experimental
value. Here we rather chose to utilize experimental energy gap information, as it is accessible,
to adjust the calculated PL spectrum by shifting the energy by a value equal to $\epsilon = E_{0-0}^{\text{exp}} - E_{0-0}^{\text{DFT}}$
to align the computed and experimental optical gaps. $E_{0-0}^{\text{exp}}$ represents the experimental value of
the optical gap reported in Ref. ⁴⁵ (1.88 eV), while $E_{0-0}^{\text{DFT}}$ is the corresponding quantity obtained
from our DFT-based calculation, i.e. the vertical excitation energy from S₀ to S₁ from the optimal
<figure><img src="image_6.png" /><figcaption>Figure 5: Emission spectra of P3HT. The blue line corresponds to a sampling of *N*=45 configurations from the dynamic trajectory, the green line to *N*=100 and the orange line to *N*=200. Experimental spectra of RR-P3HT thin films at *T* = 10 and 300 K, as reported in the literature, are depicted as black dotted (labeled "10 K")<sup>44</sup> and dashed (labeled "300 K")<sup>45</sup> lines, respectively, with a vertical shift of 0.7. The labels 0-0, 0-1, and 0-2 correspond to transitions assisted by 0, 1, or 2 phonons based on the rationale provided in Ref.<sup>45</sup></figcaption></figure>

geometry of *S*₀, equal to 0.72 eV. When we average the peak positions of the three shifted emission spectra (Table 1), we find a significant agreement with the experimental data obtained from reference<sup>45</sup> by photoluminescence (PL) on thin films of regio-regular P3HT (RR-P3HT).

An important observation is that our computed photoluminescence spectra qualitatively match with the experimental one in terms of overall shape, with well-reproduced peak positions. Notably, our calculation reveals that the intensity of the 0-0 vibronic peak is lower than that of the 0-1 peak, consistent with the experimental spectrum shown in Fig. 5. The relative intensities of the peaks in the PL spectrum of P3HT have been the subject of various studies, as discussed in reference.<sup>47</sup> Depending on the conditions of the film deposition (chloroform solution vs. toluene, film vs. nanoribbon), either the 0-0 peak (J-aggregate hypothesis) or the 0-1 peak (H-aggregate hypothesis) can be the most intense. J-aggregate refers to intrachain coupling among chromophores, while in
<table><caption>Table 1: Position in energy of the peaks in the calculated and measured emission spectra. The transition without phonon (0-0 peak) is adjusted to the experimental value (at T=0 K) to overcome the bandgap underestimation in DFT.</caption><thead><tr><th rowspan="2">Label</th><th colspan="2">Present calculation</th><th colspan="2">Measurement<sup>45</sup></th></tr><tr><th>Average(eV)</th><th>Error (eV)</th><th><i>E</i>(eV)@T=10K</th><th><i>E</i> (eV)@T=300K</th></tr></thead><tbody><tr><td>0-0</td><td>Aligned</td><td>—</td><td>1.88</td><td>1.83</td></tr><tr><td>0-1</td><td>1.66</td><td>0.03</td><td>1.70</td><td>1.68</td></tr><tr><td>0-2</td><td>1.55</td><td>0.02</td><td>1.53</td><td>1.54</td></tr></tbody></table>

H-aggregate, interchain coupling dominates. The fact that our model's PL spectrum resembles an H-aggregate type can be rationalized by considering the typical coherence length of P3HT, estimated to be about ~ 6.5 thiophene units.⁴⁸ In our simulated system, only a single chromophore of this length is present along the P3HT chain, limiting interactions to those between chromophores located on different chains along the stacking direction, i. e. an H-aggregate.

The radiative rate of the exciton in the $S_1$ state is determined by averaging the emission spectrum over the entire energy interval sampled during transitions from $S_1$ to $S_0$. The corresponding exciton lifetime is obtained as the inverse of this rate. Since the emission spectrum includes various vibronic peaks (0-0, 0-1, etc.), we compute the variance for each peak and then average them to estimate the total uncertainty on $\tau_{rad}$. The resulting values are presented in Table 2. With an

<table><caption>Table 2: Exciton radiative lifetime calculated for three sampling of the $S_1$ trajectory.</caption><thead><tr><th><i>N</i></th><th>τ<sub>rad</sub> (ps)</th><th>error (ps)</th></tr></thead><tbody><tr><td>45</td><td>682</td><td>172</td></tr><tr><td>100</td><td>496</td><td>129</td></tr><tr><td>200</td><td>464</td><td>118</td></tr></tbody></table>

increase in sampling size (represented by $N$), the error bar diminishes slightly. The calculation show a clear trend resulting to a value of $\tau_{rad}$ of 464± 118 ps for the densest sampling ($N = 200$). Despite visible variations of the overall line shape, the peaks positions can be clearly identified and it is evident that their relative amplitude depend on the sampling. In our case, 200 configurations is the best we can achieve being the upper limit of uncorrelated configuration sampling. Tables 1 and 2 complement the information of Fig. 5, confirming the quality of the estimation of peaks positions and radiative lifetime with associated error bar.
Beside the results reported in the SI mentioned above, an additional check for possible size effects was done by performing calculations on a larger system containing eight thiophenes per strand in the simulation cell at the PBE level. This result is shown Figure SI.4 reported in the supporting information. Namely, the peaks position and the intensity for the 0-0 transition lower than the 0-1 are well reproduced features independent of the system size. For the radiative recombination time, we get a value of 490 ± 103 ps, still compatible with the estimation of 464 ± 118 ps obtained for the smaller model, thus supporting the reliability of our approach and a numerically negligible dependence on the system size.

<figure><img src="image_7.png" /><figcaption>Figure 6: Transition matrix elements calculated according to Eq. 3 along the dynamics. The blue line indicates the instantaneous value $A_{S_1 \to S_0}^n(t)$ in the case of a sampling with $N=200$ configurations. The green dotted line is the corresponding average, resulting in a recombination time of $\tau_{rad} = 475$ ps (Table 2). The red line corresponds to the transition matrix element in the relaxed geometry of the $S_1$ state ($A_{S_1 \to S_0}^{\text{geo } S_1}$).</figcaption></figure>

Figure 6 illustrates the importance of considering the time evolution of the system in the $S_1$ state. The instantaneous value of the transition matrix element $A_{S_1 \to S_0}^n(t)$ to the $S_0$ state (in blue) for each of the $N = 200$ configurations changes considerably during the dynamics. The mean value $\langle A_{S_1 \to S_0}^n \rangle$, employed in the calculation of the value of $\tau_{rad}$ as summarized in Table 2, is represented by the dotted green line. In this figure, the red line corresponds to the transition matrix element of the system when it is in the relaxed geometry of state $S_1$ ($A_{S_1 \to S_0}^{\text{geo } S_1}$). The value is considerably smaller, with a corresponding radiative lifetime $\tau_{rad} = 12$ µs, and would result in a diffusion length
$L=3$ µm drastically larger than the experimental data. It is therefore imperative to avoid limiting the calculations to static descriptions of a single configuration, including the $S_1$ minimum energy one, even when DFT approaches are used to delineate the potential energy of the atomic system. Instead, one should properly account for the dynamics of the system at finite temperature to ensure quantitative agreement with the experimentally determined exciton diffusion length.

As mentioned in the introduction, the exciton lifetime (or photoluminescence lifetime) depends on a number of experimental parameters. In particular, its value can vary depending on whether the P3HT is isolated or blended with another material. When P3HT alone is considered, as in our study, there is the additional difficulty of the simultaneous presence of crystalline and amorphous domains in the experimental samples. In this context, the experimental efforts reported in Refs. 14–16 target crystalline domains specifically, and are therefore the ideal reference for validating our calculation. In the works of Ohkita et al., 15,16 the exciton diffusion coefficient $D$ and diffusion length $L$ are carefully reported, but the lifetime $\tau$ is not given explicitly, as well as the specific radiative contribution $\tau_{rad}$. In a former study, 9 we were able to get an excellent agreement of the exciton diffusion coefficient calculated for the atomic model used in the present work ($D = 8.0 \times 10^{-3} \text{ cm}^2 \text{ s}^{-1}$) with their measurement. In the absence of an analogous precise value of $\tau$ in Refs., 15,16 we use as a reference the exciton diffusion length $L$. Our calculation for this quantity yields $L = \sqrt{D\tau_{rad}} = 19.3 \pm 2.3 \text{ nm}$, assuming that only radiative processes occur and being fully aware that we are ignoring non-adiabatic processes. Despite our declared limits, we notice that our estimation in fairly good agreement with their value of $L = 20 \text{ nm}$. Furthermore, our result is also in line with other values available in the literature 14 ($L = 27 \pm 12 \text{ nm}$). This agreement, although satisfactory in terms of the validation of our calculation protocol, is also intriguing because non-radiative recombinations were not accounted for in our theoretical approach. This might indicate that these processes are negligible in the crystalline part. This hypothesis could be verified by determining the quantum yield specifically in the crystalline part, and not in the P3HT material as a whole. Nonetheless, such a fine discrimination seems difficult to implement and currently escape experimental probes. If it turns out that a scattering length of around 20 nm is indeed typical of
exciton scattering and recombination by essentially radiative processes in crystalline zones, this
would imply that the reduction of *L* to values of a few nm in blends¹⁰ could be due to its crossing
of amorphous zones where non-radiative phenomena would occur. This remains an open issue
calling for further experimental and computational developments.

CONCLUSIONS

In summary, the method introduced in this study for determining the radiative lifetime of the S<sub>1</sub> exciton state proves its reliability through comparison with experimental results. Additionally, it demonstrates versatility by not depending on specific model characteristics, thereby making it widely applicable to various semiconducting polymers. When combined with the estimation of the diffusion coefficient (as discussed in<sup>9</sup>), our approach also allows for the extraction of the diffusion length, which demonstrates favorable agreement with existing experimental data. Hence, the approach is ideally suited for uncovering the relationship between structural features, such as molecular stacking degree, and exciton diffusion prior to its recombination via light emission. Furthermore, it can be applied to various organic materials, including both polymers and assemblies of small molecules. The overall picture that emerges from the work presented here, theoretically well grounded, is a viable alternative to computationally more expensive methods such as TDDFT. Not only the method to extract the lifetime is generalizable to other systems, but it provides also a direct inspection of the atomistic and electronic mechanism not accessible to empirical or phenomenological methods as, for instance, the Förster resonance energy transfer. Having a dynamical atomistic approach, as opposed to static ones, to get the complete picture in terms of energetics, dynamics of the underlying mechanisms and both lifetime and diffusion constants seems to be crucial. In fact, a proper account of the underlying dynamics, accessible by the proposed approach, is essential for a process that is by its own nature dynamical and we have shown that a simple static picture can lead to results off of orders of magnitudes in terms of diffusion length and related lifetime.
# ACKNOWLEDGEMENTS

The authors acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsac calcul/Big Data. Calculations were also performed by using resources from GENCI (Grand Equipement National de Calcul Intensif) (grant no. 0910296 and 0905071). This work of the Interdisciplinary Institute HiFunMat, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm, was supported by IdExUnistra (ANR-10-IDEX-0002) and SFRI (STRAT'US project, ANR-20-SFRI-0012) under the framework of the French Investments for the Future Program.

# References

(1) Zhu, L. et al. *Nature Materials* **2022**, *21*, 656–663.

(2) Chong, K.; Xu, X.; Meng, H.; Xue, J.; Yu, L.; Ma, W.; Peng, Q. *Advanced Materials* **2022**, *34*, 2109516.

(3) Chen, H.; Jeong, S. Y.; Tian, J.; Zhang, Y.; Naphade, D. R.; Alsufyani, M.; Zhang, W.; Griggs, S.; Hu, H.; Barlow, S.; Woo, H. Y.; Marder, S. R.; Anthopoulos, T. D.; McCulloch, I.; Lin, Y. *Energy & Environmental Science* **2023**, *16*, 1062–1070.

(4) Chen, H. et al. *Nature Materials* **2025**, *25*, https://doi.org/10.1038/s41563-024-02062-0.

(5) Hernández, F. J.; Crespo-Otero, R. *Annual Review of Physical Chemistry* **2023**, *74*, 547–571.

(6) Binder, R.; Burghardt, I. *Faraday Discuss.* **2020**, *221*, 406–427.

(7) Peng, W.-T.; Brey, D.; Giannini, S.; Dell'Angelo, D.; Burghardt, I.; Blumberger, J. *J. Phys. Chem. Lett.* **2022**, *13*, 7105–7112.
(8) Giannini, S.; Peng, W.-T.; Cupellini, L.; Padula, D.; Carof, A.; Blumberger, J. *Nature Comm.* **2022**, *13*, 2755.

(9) Diarra, C. O.; Boero, M.; Steveler, E.; Heiser, T.; Martin, E. *Phys. Chem. Chem. Phys.* **2023**, *25*, 15539–15546.

(10) Shaw, P. E.; Ruseckas, A.; Samuel, I. D. W. *Advanced Materials* **2008**, *20*, 3516–3520.

(11) Dimitriev, O. P.; Blank, D. A.; Ganser, C.; Teichert, C. *The Journal of Physical Chemistry C* **2018**, *122*, 17096–17109.

(12) Classen, A.; Chochos, C. L.; Lüer, L.; Gregoriou, V. G.; Wortmann, J.; Osvet, A.; Forberich, K.; McCulloch, I.; Heumüller, T.; Brabec, C. J. *Nature Energy* **2020**, *5*, 711–719.

(13) Dimitriev, O. P. *Chemical Reviews* **2022**, *122*, 8487–8593.

(14) Cook, S.; Liyuan, H.; Furube, A.; Katoh, R. *The Journal of Physical Chemistry C* **2010**, *114*, 10962–10968.

(15) Tamai, Y.; Matsuura, Y.; Ohkita, H.; Benten, H.; Ito, S. *The Journal of Physical Chemistry Letters* **2014**, *5*, 399–403.

(16) Ohkita, H.; Tamai, Y.; Benten, H.; Ito, S. *IEEE Journal of Selected Topics in Quantum Electronics* **2016**, *22*, 100–111.

(17) Kayunkid, N.; Uttiya, S.; Brinkmann, M. *Macromolecules* **2010**, *43*, 4961–4967.

(18) Brinkmann, M. *Journal of Polymer Science Part B: Polymer Physics* **2011**, *49*, 1218–1233.

(19) Car, R.; Parrinello, M. *Physical Review Letters* **1985**, *55*, 2471–2474.

(20) CPMD, Copyright IBM Corp. 1990–2022, copyright MPI für Festkörperforschung Stuttgart 1997–2001. Available on GitHub under MIT License.

(21) Perdew, J. P.; Burke, K.; Ernzerhof, M. *Phys. Rev. Lett.* **1996**, *77*, 3865–3868.
(22) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. *The Journal of Chemical Physics* **2010**, *132*,
154104.

(23) Grimme, S. *Journal of Computational Chemistry* **2006**, *27*, 1787–1799.

(24) Troullier, N.; Martins, J. L. *Physical Review B* **1991**, *43*, 1993–2006.

(25) Nosé, S. *Molecular Physics* **1984**, *52*, 255–268.

(26) Nosé, S. *The Journal of Chemical Physics* **1984**, *81*, 511–519.

(27) Hoover, W. G. *Physical Review A* **1985**, *31*, 1695–1697.

(28) Martyna, G. J.; Klein, M. L.; Tuckerman, M. *The Journal of Chemical Physics* **1992**, *97*,
2635–2643.

(29) Blöchl, P. E.; Parrinello, M. *Physical Review B* **1992**, *45*, 9413–9416.

(30) Frank, I.; Hutter, J.; Marx, D.; Parrinello, M. *The Journal of Chemical Physics* **1998**, *108*,
4060–4069.

(31) Grimm, S.; Nonnenberg, C.; Frank, I. *The Journal of Chemical Physics* **2003**, *119*, 11574–11584.

(32) Nonnenberg, C.; Grimm, S.; Frank, I. *The Journal of Chemical Physics* **2003**, *119*, 11585–11590.

(33) Röhrig, U. F.; Frank, I.; Hutter, J.; Laio, A.; VandeVondele, J.; Rothlisberger, U.
*ChemPhysChem* **2003**, *4*, 1177–1182.

(34) Odelius, M.; Laikov, D.; Hutter, J. *Journal of Molecular Structure: THEOCHEM* **2003**, *630*,
163–175.

(35) Kowalczyk, T.; Tsuchimochi, T.; Chen, P.-T.; Top, L.; Van Voorhis, T. *The Journal of Chemical Physics* **2013**, *138*, 164101.
(36) Büchel, R.; Álvarez, L.; Grage, J.; Maniscalco, D.; Frank, I. *Molecules* **2024**, 29.

(37) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. *Journal of Chemical Physics* **2003**, 118, 8207–8215.

(38) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. *Journal of Chemical Physics* **2006**, 124, 219906.

(39) Falletta, S.; Pasquarello, A. *Phys. Rev. B* **2022**, 106, 125119.

(40) Boero, M.; Parrinello, M.; Terakura, K.; Ikeshoji, T.; Liew, C. C. *Physical Review Letters* **2003**, 90, 226403.

(41) Boero, M. *The Journal of Physical Chemistry A* **2007**, 111, 12248–12256.

(42) Jones, G. A.; Acocella, A.; Zerbetto, F. *The Journal of Physical Chemistry A* **2008**, 112, 9650–9656, PMID: 18767783.

(43) Ibele, L. M.; Curchod, B. F. E. *Phys. Chem. Chem. Phys.* **2020**, 22, 15183–15196.

(44) Paquin, F.; Yamagata, H.; Hestand, N. J.; Sakowicz, M.; Bérubé, N.; Côté, M.; Reynolds, L. X.; Haque, S. A.; Stingelin, N.; Spano, F. C.; Silva, C. *Phys. Rev. B* **2013**, 88, 155202.

(45) Clark, J.; Silva, C.; Friend, R. H.; Spano, F. C. *Phys. Rev. Lett.* **2007**, 98, 206406.

(46) Perdew, J. P. *Int. J. Quantum Chem.* **1985**, 28, 497–523.

(47) Spano, F. C.; Silva, C. *Annual Review of Physical Chemistry* **2014**, 65, 477–500.

(48) Raithel, D.; Simine, L.; Pickel, S.; Schötz, K.; Panzer, F.; Baderschneider, S.; Schiefer, D.; Lohwasser, R.; Köhler, J.; Thelakkat, M.; Sommer, M.; Köhler, A.; Rossky, P. J.; Hildner, R. *Proceedings of the National Academy of Sciences* **2018**, 115, 2699–2704.
Graphical TOC Entry

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