Tracking the coherent generation of polaron pairs in conjugated polymers

The optical excitation of organic semiconductors not only generates charge-neutral electron-hole pairs (excitons), but also charge-separated polaron pairs with high yield. The microscopic mechanisms underlying this charge separation have been debated for many years. Here we use ultrafast two-dimensional electronic spectroscopy to study the dynamics of polaron pair formation in a prototypical polymer thin film on a sub-20-fs time scale. We observe multi-period peak oscillations persisting for up to about 1 ps as distinct signatures of vibronic quantum coherence at room temperature. The measured two-dimensional spectra show pronounced peak splittings revealing that the elementary optical excitations of this polymer are hybridized exciton-polaron-pairs, strongly coupled to a dominant underdamped vibrational mode. Coherent vibronic coupling induces ultrafast polaron pair formation, accelerates the charge separation dynamics and makes it insensitive to disorder. These findings open up new perspectives for tailoring light-to-current conversion in organic materials.


Supplementary Tables
Supplementary Table 1  3. Dynamics of two-dimensional electronic spectra of rr-P3HT thin films for waiting times up to 1 ps. To study the persistence of vibronic coherence on longer timescales, we measured the dynamics of absorptive two-dimensional optical spectra for waiting times of up to 900 fs with a stepsize of 5 fs. In the measurements, the maximum energy of the pump pulse pair was set to 5 nJ at τ = 0 fs, whereas the probe energy was 1. a-d insets, dashed lines). The peak at 1380 cm -1 corresponds to the symmetric C-C stretching mode and is also seen in the differential transmission measurements in Supplementary Figure   3. The peak at 1520 cm -1 may be assigned to the anti-symmetric C=C stretching mode of thiophene 7 . Due to the larger stepsize along the waiting time axis, the signal-to-noise ratio in these measurements is reduced in comparison to those shown in Fig. 1 of the main text.

Details of the experimental spectral splittings up to 60 fs
To better show that each of the exciton peaks predicted by the incoherent model 6  Additionally, to show that the subpeak structure is indeed visible up to about 60 fs, we also report 2DES maps at later waiting times (Supplementary Figure 6) with their respective cross sections at 2.25 eV. From these maps and cross sections one can clearly recognize that the peak splitting along the detection energy persists for up to T = 58 fs.

Theoretical model
The Hamiltonian. To analyze the experimentally observed linear absorption and 2DES spectra and to understand the basic physical mechanisms underpinning these observations, we consider a phenomenological model for P3HT consisting of four electronic states (ground state, exciton, polaron pair and excited polaron pair) coupled to a vibrational degree of freedom provided by the C=C stretch mode. This reduced model is described by the Hamiltonian of the form + P P (!Ω P + !ωa † a + !ω S P (a † + a) + !ωS P ) + ( X P + P X )(!J + !ω S XP (a † + a)) + P * P * (!Ω P* + !ωa † a + !ω S P* (a † + a) + !ωS P* ).
Here G , X , P , P * denote ground state, exciton, polaron pair and excited polaron pair state, respectively, with the purely electronic energies represented by Ω X , Ω P , Ω P* . The site energy of each electronic state is given by the sum of the purely electronic energy Ω k and the reorganization energy ωS k due to the vibronic coupling. The electronic coupling between exciton and polaron pair denoted by J is responsible for the coherent population transfer between exciton and polaron pair. The C=C stretch mode is modeled by an effective, single vibrational mode (annihilation and creation operators denoted by a and a † , respectively) with vibrational frequency of ω = 1450 cm -1 and a vibrational relaxation rate of γ vib = (650 fs) -1 , both of which are determined based on experimentally observed oscillatory pump-probe and 2D signals. The interaction between electronic states and the C=C stretch mode is modeled by vibronic couplings quantified by Huang-Rhys (HR) factors that we denote as S X , S P , S P* , S XP .
Here S X , S P , S P* describe the displacement of the equilibrium position of the vibrational mode in the electronic excited states with respect to the equilibrium position of the mode in the electronic ground state. These three vibronic couplings describe the modulation of the excited state energy ! " Ω k = !Ω k + !ωa † a + !ω S k (a † + a) + !ωS k (cf. k = X,P,P *) due to the vibrational motion. It is notable that when S X ≠ S P , the energy-level difference ! Ω X -! Ω P between exciton and polaron pair states depends on ω( S X -S P )(a † + a) , implying that nonequilibrium vibrational motion can dynamically modulate the resonance between exciton and polaron pair states significantly when the difference in the corresponding HR factors is sufficiently large 8 . On the other hand, S XP describes the modulation of the electronic coupling strength between exciton and polaron pair due to the vibrational motion. In the case of natural photosynthetic systems, it has been shown that typically the vibronic modulation of the energy levels leads to larger HR factors than that modulating electronic couplings 9 . In simulations, we found that larger S X , S P , S P* than S XP are needed to reproduce experimental observations within our model.

Decoherence.
The unavoidable coupling of the electronic states to their environments destroys electronic coherences and induces homogeneous broadening and incoherent features in the population transfer dynamics between excited states. In the absence of a detailed microscopic model for the dephasing of vibronic excitations in the polymer, the influence of the noisy environment on electronic states can be modeled phenomenologically by dissipative (i.e. non-unitary) contributions to the time evolution of the vibronic system that we assume can be expressed in the so-called Lindblad form: Here the first Lindblad super-operator describes pure dephasing noise with a rate of γ dep , which contributes to homogeneous broadening but does not induce population transfer between exciton and polaron pair states by itself. The second Lindblad super-operator describes the incoherent population relaxation from the exciton to the polaron pair state with a rate of γ rxn . The vibrational relaxation of the C=C stretch mode is modeled by a formally analogous term where the relaxation rate of the vibration is taken to be γ vib = (650 fs) -1 based on experimental results. To take into account the fact that the equilibrium position of the vibrational mode depends on the electronic states, we introduce a displaced mode operator  a = U † aU with a k=X,P,P* ∑ in terms of the so-called displacement operator D(x) = exp(x(a † -a)) . The Lindblad super-operator in terms of  a therefore describes the relaxation of the vibrational motion to different equilibrium positions depending on the electronic state of the system. When the HR factors are sufficiently small,  a is well approximated by a , leading to a widely used Lindblad operator for mode damping.
However, due to the large HR factors of the present system, such an approximation is not valid here.
For the Hamiltonian and noise model described above, the time evolution of the global, vibronic (electronic + vibrational) state ρ(t) is governed by the quantum master equation We note that the above phenomenological model is employed to demonstrate how coherent electronic and vibronic couplings affect 2D optical responses in the presence of noise and disorder. We found that, remarkably, several relevant experimental features, such as multifrequency transients and splitting of 2D peaks, can be reproduced even though the employed model is rather simple. A detailed microscopic model based on electronic wavefunctions and realistic phonon spectral densities may be able to reproduce more detailed features of experimental observations, such as the clean splitting of 2D peaks and the spectral diffusionlike effects discussed in the manuscript, and also provide a more accurate description of the transient dynamics of 2D spectra.
Absorptive 2D spectra. To calculate the absorptive 2D spectra of the present system, we consider the transition dipole between the ground and exciton state and that between the polaron pair and excited polaron pair states, denoted by ! µ X and ! µ P , respectively, with the assumption that the transition dipole strength between ground state and polaron pair is sufficiently weak. In this case, the optical response from the polaron pair state can be observed mediated by the optical transition between polaron pair and excited polaron pair states, leading to excited state absorption (ESA) signals in pump probe and 2D spectra.
The third-order polarization responsible for 2D spectra 10 is described by where the electric field E(t) is modeled by a sum of three Gaussian fields to fit the laser spectrum used in experiments, as shown in Supplementary Figure 7 and Supplementary Table   1, with a carrier frequency of υ and a slowly varying, time-dependent function ε(t) where c.c denotes the complex conjugate. Here τ, T , t represent the time delays between the first and second excitation pulses, between the second and third excitation pulses, and between the third excitation pulse and the emerging signal from the system, respectively. The Fourier transformation of the third-order polarization with respect to (τ , t) leads to 2D spectra as a function of excitation and detection energies (E X , E D ) for each waiting time T .
In this work, we calculated absorptive 2D spectra within the rotating wave approximation 7 by taking into account only the signal contributions emitted into two phase-matched directions We note that for our model, the influence of the transition dipoles and their rotational averaging (ensemble) on the third-order optical response S (3) (t 3 ,t 2 ,t 1 ) is characterized by two effective transition dipole strengths. More specifically, here we consider the ground state bleaching (GSB), stimulated emission (SE) and excited state absorption (ESA) components of the rephasing pathway, whose contributions to the third-order optical response are given by Here ρ eq denotes the equilibrium state of the vibronic system in the electronic ground state, u(t) is a formal expression of the time evolution operator describing both the Hamiltonian and noise, and µ X,P ± are the transition dipole operators defined by (13) with ê denoting the unit vector parallel to the polarization of the electric field. In this case, the third-order optical responses are reduced to Therefore the rotational average of the transition dipoles (ensemble) with respect to the polarization direction of the electric field is characterized by two parameters 2 , where µ GSB,SE and µ ESA determine the relative intensity between the GSB/SE contributions and ESA contribution, respectively. In this work, we take µ ESA ≈ 0.2 µ GSB,SE to take into account the fact that the experimentally observed (negative-valued) ESA signals at detection energies around 1.89 eV are weaker than the (positive-valued) GSB/SE signals at detection energies between 2.0 eV and 2.4 eV (cf. eV). We found that when the FWHM of the inhomogeneous broadening of the electronic energy difference Ω X -Ω P becomes larger than ~50 meV, there are notable changes in simulated absorption and 2D spectra, but we found decreasing agreement with the experimental results. We also considered an alternative inhomogeneous broadening model where each value of the electronic energies Ω X , Ω P , Ω P* is modeled by an independent Gaussian distribution, but could not find a better match between theory and experiments. We We found that these results do not depend on the choice of the rates of pure dephasing noise and relaxation and also the vibronic coupling strengths quantified by S X , S P , S P* . In addition, the incoherent model does not predict low-frequency components ( < 1000 cm -1 ) in transients, as shown in Supplementary Figure 10d and e, contrary to the results of the coherent model (cf. Fig. 3 in the manuscript and Supplementary Figure 9). These results support our conclusions that the splitting of 2D peaks and multi-frequency transients observed in experiments are the evidence of coherent couplings in the present system. We note that the model parameters estimated in this work have a significant uncertainty as there are many different sets of model parameters that reproduce the absorption spectrum quite well and also induce the splitting of homogeneously broadened 2D peaks in the absence of inhomogeneous broadening. For all cases, we found that coherent electronic and vibronic couplings are required with different values of S X and S P to obtain the above mentioned features within our phenomenological model: a sufficiently large HR factor S X of the exciton state is required to reproduce strong vibrational progressions in the absorption spectrum, while a moderate difference in S X and S P and other model parameters in an appropriate parameter regime are required to reproduce the splitting of 2D peaks observed in experiments. We note that the splitting of 2D peaks is robust against the variation of model parameters in simulations, due to the strong electronic and vibronic couplings of the present system (cf. Supplementary Table 2). We assume in this scenario that the electronic coupling between exciton and polaron pair is in the weak coupling limit. In this model, we first deduce a maximum Huang-Rhys factor of 0.05 for the 730 cm -1 mode by analyzing Fourier transforms of pump-probe dynamics such as those shown in Supplementary Figure 3. We then calculate 2DES maps and observe that these show the characteristic peaks of the incoherent exciton model (Supplementary Figure 11a). A characteristic peak splitting, such as the one shown in Fig. 1 of the manuscript, is not observed ( Supplementary Figure 11a-b). In order for the 730 cm -1 mode to induce the peak splittings observed in our experiment, the Huang-Rhys factor should exceed 1 (Supplementary Figure 11c). In this case, however, the dynamics of all peaks would show large and long-lived beatings along the waiting time (Supplementary Figure 11d). These strong beatings are clearly not detected in our measured 2DES experiment (Fig. 2). Hence, we conclude that the Huang-Rhys factor of the 730 cm -1 vibrational mode is too weak to account for the experimentally observed peak splittings. + P X D( S P )(!J + !ω S XP (a † + a))D † ( S X ) + P * P * (!Ω P* + !ωa † a).

Quantification of mixing between exciton, polaron pair and C=C stretch mode
Here the vibrational energy term ωa † a describes the quantized vibrational levels of the C=C stretch mode in each electronic state manifold. The time evolution of the vibronic state ! ρ(t) = Uρ(t)U † in the displaced basis is governed by a master equation that is also of the Lindblad form. We can therefore write where the electronic and vibrational relaxation terms are given by (20) In the displaced basis, the transition dipole operators describing the optical transition between ground state and exciton can be expressed as follows where the displacement operator D( S X ) describes the Franck-Condon (FC) factors between vibrational levels in ground state and those in exciton state manifold, denoted by G,n and X,n , respectively: here n = 0,1, 2, , and the vibrational levels in electronic state k are defined by k k ⊗ ωa † a ( ) k,n = ωn k,n . For instance, the transition dipole strength between G,0 and X,n , the so-called 0 -n transition, is described by a coefficient n D( S X ) 0 .

Bright vibronic eigenstates of the Hamiltonian.
Due to the high frequency of the C=C stretch mode, which is approximately seven times higher than the thermal energy at room temperature, the initial state of the vibronic system is well approximated by G,0 . In the displaced basis, the vibronic eigenstates of the Hamiltonian are formally expressed as XP k = (ψ X,n (k) X,n + ψ P,n (k) P,n ) n=0 ∞ ∑ defined by ! H XP k = "Φ k XP k with the associated energy level of Φ k . The energy-level difference between XP k and G,0 , given by Φ k , determines the location of the absorptive or 2D peak originating from the optical transition from G,0 to XP k .
For the model parameters used in Fig. 3 of the manuscript (cf. Supplementary Table 2), Supplementary Table 5 Fig. 3f) and homogeneously broadened 2D peaks (cf. Fig. 3b, c), which are close to the location of experimentally observed splitting of 2D spectra.
In order to quantify the degree of coherent mixing of exciton and polaron pair states in each vibronic eigenstate XP k , we introduce an operator σ z = P P -X X . The expectation value of σ z shows the relative amplitudes in exciton and polaron pair states, which ranges from -1 (for a pure exciton state) to 1 (for a pure polaron pair state). An expectation value σ z close to 0 indicates strong mixing of exciton and polaron pair. Here we quantify the mixing of exciton and polaron pair by using the variance of σ z , defined by ∑ j denotes the reduced electronic state of the vibronic eigenstate XP k . In Supplementary Table 5, the degree of exciton-polaron-pair mixing of the vibronic eigenstates is displayed.
In Supplementary Table 5, the first two dominant amplitudes of each vibronic state in the displaced basis are marked in red, where the first two lowest energy vibronic states are approximately given by XP k=1,2 ≈ ψ X,0 (k) X,0 + ψ P,1 (k) P,1 , while the higher energy vibronic states are approximately given by XP k=3,4 ≈ ψ X,1 (k) X,1 + ψ P,2 (k) P,2 , because the electronic energy of the polaron pair state is lower than the exciton energy (see Supplementary Table 2 -4ρ XX (k) +1) . This implies that all the vibronic states cannot be written as product states in the form of (ψ X X + ψ P P ) ⊗ φ with a vibrational state φ , indicating that the vibronic states are the results of the coherent mixing among exciton, polaron pair and vibrational degrees of freedom due to the coherent electronic and vibronic couplings.
We note that in the absence of coherent electronic and vibronic couplings, the coherent mixing between exciton and polaron pair states is absent, for which the eigenstates of the Hamiltonian consist of bright exciton states (ψ P,n (k) = 0) and dark polaron pair states (ψ X,n (k) = 0) , which do not lead to the splitting of absorptive and 2D peaks in the simulations (cf. Supplementary Figure 10).
The mixing of electronic and vibrational degrees of freedom essentially results from the fact that the electronic and vibronic couplings in the Hamiltonian are diagonalized in different bases. In order to see this in more detail, we consider the electronic Hamiltonian H e = !(Ω X + ωS X ) X X + !(Ω P + ωS P ) P P (22) + !J( X P + P X ) + !(Ω P* + ωS P* ) P * P * .

The role of non-equilibrium vibrational motion in polaron pair formation.
With the model parameters estimated based on experimentally measured absorption and 2D spectra, we now investigate the time evolution of the vibronic system and study the role of non-equilibrium vibrational motion of the C=C stretch mode in the polaron pair dynamics.
Here we take the initial state of the vibronic system to be X ⊗ 0 G where 0 G is the vibrational ground state in the electronic ground state manifold. The vertical optical transition from ground state to the exciton state leads to such an initial state.
In Supplementary Figure 12a, we consider the case in which the polaron pair dynamics is governed by coherent couplings in the absence of relaxation (cf. J ≠ 0, S XP ≠ 0, γ rxn = 0 , see the model parameters summarized in Supplementary Table 3). Here the population dynamics of the polaron pair state (shown in black) shows oscillatory behavior with a period of ~50 fs, which is close to the frequency of the low-frequency oscillations ( < 1000 cm -1 ) shown in Supplementary Figure 9d, e. The electronic coherence between exciton and polaron pair, which is quantified by X ρ e (t) P + P ρ e (t) X for the reduced electronic state ρ e (t) , is shown in blue, while the entanglement between electronic and vibrational degrees of freedom, quantified by the logarithmic negativity 15 , is displayed in red. We note that due to the pure dephasing noise, the initial vibronic state becomes a statistical mixture of pure states, called a mixed state, for which the logarithmic negativity is the most suitable and easy to compute measure that quantifies the degree of entanglement. While it also applies to the pure state case, there it is customary to employ the von Neumann entropy due to its explicit operational meaning (cf. Supplementary Table 5). These results show the coherent features in the polaron pair rise dynamics. In particular, the non-zero entanglement implies that vibrational motion is quantum-mechanically correlated with electronic motion.
In Supplementary Figure 12b In Supplementary Figure 12c and d, we consider the parameters used in Supplementary   Figure 12a and b, respectively, but inhomogeneous broadening of the electronic energy difference between exciton and polaron pair is now included by means of considering a Gaussian distribution with a FWHM of 200 meV. It is notable that even in the presence of a large disorder in the electronic energy difference, the overall system dynamics are not changed significantly due to the strong electronic and vibronic couplings. Supplementary

Acceleration of polaron pair formation due to vibronic coupling
In the absence of vibronic coupling, the dynamics of exciton X and polaron pair P states is equivalent to that of two coupled electronic states with the vibrational mode being a simple spectator. Thus, if we assume that the excitation pulse initially prepares the system in the X state, we will observe population oscillations between X and P with a period T = π Ω R , where !Ω R = 1 2 4J 2 + ΔE 2 is the Rabi energy defined by the electronic coupling strength J and the energetic detuning E Δ between X and P states. The dynamics of the polaron pair population can be described as From this equation we can see that the maximum oscillation amplitude is reached when the bare uncoupled states X and P are perfectly resonant, i.e. 0 E Δ = . In this case, the population oscillates completely between X and P and the transfer rate is maximal.
Conversely, any energetic detuning between the bare states will reduce the amplitude of these oscillations while slightly increasing the oscillation frequency.
In the polymer samples, rapid electronic dephasing will quickly dampen these coherent population oscillations. To estimate the dynamics of the polaron pair formation in the absence of vibronic coupling ( Δ XP = 0 ), we have calculated the population dynamics from the time evolution of the density operator and then estimated the rise time of these dynamics for the polaron pair state by fitting a simple exponential rise function to the data. We assume an electronic dephasing time of 50 fs for both X and P states.
We observe that, in presence of dissipation, the populations oscillations are additionally exponentially damped with a rate defined by the electronic dephasing time γ .
As can be seen in Supplementary Figure 13, this population rise time varies from less than 10 fs when the states are resonant (zero detuning) to more than 200 fs for detuning of a few hundred meV. Hence, in this model, the energetic detuning between X and P has a large effect on the polaron pair formation.
In the polymer samples, the position of exciton and polaron pair states, and thus their relative energetic detuning, is primarily determined by the respective binding energy of the electronhole pair. DFT studies have reported typical values of detuning of a few hundred meV for conjugated polymers 18,19 . Additionally, inherent disorder in the samples induces local fluctuations of the state energy which further increases the energetic detuning between X and P states. Thus, in the absence of vibronic coupling, one might expect that polaron pair formation occurs on much longer time scales of several hundred fs. Another important point to notice is that, in the absence of vibronic coupling, the transfer integral between X,ν X and P,ν P states with different vibrational quantum numberν is always zero and is instead unity only for Δν XP = ν X −ν P = 0 . This essentially means that vibrations do not play a role in the dynamics and act as spectator modes. This is fundamentally different in the presence of vibronic coupling, i.e., in the presence of a finite displacement Δ XP between X and P states. Hereby, the transfer integrals become nonzero also for states with Δν XP ≠ 0 and thus population transfer will be allowed also between states with different vibrational quanta. The larger Δ XP , the larger is the overall transfer integral between the states and the faster the transfer rate. Physically, this is readily understood. Due to the ladder of vibrational states in the X and P states, the energetic displacement between a certain donor mode X,ν X and the nearest acceptor mode P,ν P can never be more than one half of a vibrational quantum. Hence, the energetic detuning between donor and acceptor has only a very moderate effect on the transfer dynamics. One may also say that, in the coupled (eigenstate) basis, mixing of electronic X and P states with the vibrational mode, opens up new paths for coherent population transfer through delocalized vibronic states. Thus, in the presence of strong vibronic coupling, rather fast population transfer is possible even for large energetic detuning between the electronic states.
In Supplementary Figure 13, we show the polaron pair rise time simulated for various displacements Δ XP . This suggests that strong vibronic coupling indeed speeds up polaron pair formation dynamics even in presence of relatively large disorder. Supplementary Figure 13 shows quite clearly that the vibronic coupling makes the population transfer dynamics insensitive to the precise value of the detuning and hence to the precise local realization of the disorder potential in the polymer sample. Our results suggest that this effect not only occurs in the samples that we have studied in this manuscript but is a general feature of vibronically enhanced charge transfer phenomena. As such, we believe that the efficient initial charge separation observed on ultrafast time scales is relevant not only to this particular conjugated polymer, but is a general result for systems with large electron-phonon couplings.