UltraFast Relaxation, Decoherence and Localization of Photoexcited States in Conjugated Polymers: A TEBD Study
Abstract
The exciton relaxation dynamics of photoexcited electronic states in poly(phenylenevinylene) (PPV) are theoretically investigated within a coarsegrained model, in which both the exciton and nuclear degrees of freedom are treated quantum mechanically. The FrenkelHolstein Hamiltonian is used to describe the strong excitonphonon coupling present in the system, while external damping of the internal nuclear degrees of freedom are accounted for by a Lindblad master equation. Numerically, the dynamics are computed using the time evolving block decimation (TEBD) and quantum jump trajectory techniques. The values of the model parameters physically relevant to polymer systems naturally lead to a separation of time scales, with the ultrafast dynamics corresponding to energy transfer from the exciton to the internal phonon modes (i.e., the CC bond oscillations), while the longer time dynamics correspond to damping of these phonon modes by the external dissipation. Associated with these time scales, we investigate the following processes that are indicative of the system relaxing onto the emissive chromophores of the polymer: 1) Excitonpolaron formation occurs on an ultrafast time scale, with the associated excitonphonon correlations present within half a vibrational time period of the CC bond oscillations. 2) Exciton decoherence is driven by the decay in the vibrational overlaps associated with excitonpolaron formation, occurring on the same time scale. 3) Exciton density localization is driven by the external dissipation, arising from ‘wavefunction collapse’ occurring as a result of the systemenvironment interactions. Finally, we show how fluorescence anisotropy measurements can be used to investigate the exciton decoherence process during the relaxation dynamics.
I Introduction
Upon photoexcitation of a conjugated polymer with a pulse of electromagnetic radiation, the exciton state formed is no longer stabilised by the ground state nuclear geometry, leading to coupled excitonphonon dynamics as the system relaxes back to equilibrium. These dynamics have been investigated using a wide array of timeresolved spectroscopic techniques for the polymer poly(phenylenevinylene) (PPV), including fluorescence depolarization,Ruseckas et al. (2005); Grage et al. (2003); Dykstra et al. (2009) threepulse photonecho Sperling et al. (2008); Dykstra et al. (2005); Yang, Dykstra, and Scholes (2005) and coherent electronic twodimensional spectroscopy,Consani et al. (2015) with multiple time scales being identified. Based on experimental and theoretical studies,Hwang and Scholes (2011); Barford, Bittner, and Ward (2012); Grage et al. (2003) the largest time scale seen in these experiments of ps is widely accepted to correspond to Förster type exciton energy transfer between chromophores. The two shorter time scales of fs and fs, however, have been studied in less detail. While it has been suggested that these two time scales correspond to the dynamic localization of the exciton as it relaxes to the low energy chromophores,Ruseckas et al. (2005); Grage et al. (2003); Dykstra et al. (2009); Sperling et al. (2008); Dykstra et al. (2005); Yang, Dykstra, and Scholes (2005); Consani et al. (2015) to our knowledge this has yet to be confirmed by theoretical simulations.
One reason for this lack of theoretical insight is the inherent challenge in simulating the exciton relaxation dynamics numerically, due to the strong excitonphonon interactions synonymous with polymer systems. One common way to treat excitonphonon interactions within these systems is to use the one dimensional FrenkelHolstein model,Holstein (1959); *holst2; Barford and Marcus (2014) which consists of a linear array of sites through which the exciton can propagate (where each site corresponds to a moiety in the polymer chain), and where the nuclear degrees of freedom are represented by a single harmonic oscillator per site, which couples locally to the exciton. The size of the Hilbert space associated with this model grows exponentially with the number of sites in the chain, which means that computationally solving the timedependent Schrödinger equation using standard differential equation solvers becomes unfeasible for even modest chain lengths.
One way to deal with the exponentially increasing Hilbert space is to treat the nuclear degrees of freedom classically, within the socalled Ehrenfest approximation. The Hilbert space now only contains the exciton degrees of freedom and therefore grows linearly with the number of moieties in the polymer chain, making it computationally feasible to perform exciton dynamics calculations for experimentally relevant polymer chain lengths. Indeed, this technique has been applied to study the exciton relaxation dynamics of photoexcited electronic states in PPV.Tozer and Barford (2012); Sterpone and Rossky (2008); Tretiak et al. (2002) However, the Ehrenfest approximation fails because of the absence of exciton decoherence occurring during the time evolution induced by the excitonphonon coupling,Truhlar (2007); Bittner and Rossky (1995); Hwang and Rossky (2004); BedardHearn, Larsen, and Schwartz (2005); Bittner and Rossky (1997) as well as its inability to correctly describe the splitting of a nuclear wave packet when passing through a conical intersection or avoided crossing.Bittner and Rossky (1997); Prezhdo and Rossky (1997) In fact, the latter is the cause of the unphysical bifurcation of the exciton density onto separate chromophores, found to occur during the relaxation dynamics of high energy photoexcited states modeled previously within the Ehrenfest approximation.Tozer and Barford (2012) A correct description of the coupled excitonphonon dynamics therefore requires a full quantum mechanical treatment of the system, as described in this paper.
For one dimensional strongly correlated systems, the density matrix renormalization group technique has been successful in obtaining accurate numerical results for properties associated with the ground state of the system.Schollwöck (2005, 2011) Using the same Hilbert space truncation scheme, the time evolving block decimation (TEBD) techniqueDaley et al. (2004); Vidal (2003, 2004) is an efficient and accurate method for computing the short time dynamics of such systems. One appealing feature of this method is that the associated computational cost scales linearly with the number of ‘sites’ in the system,Daley et al. (2004); Vidal (2004) allowing in principle the quantum dynamics for large systems to be simulated. Indeed, the technique has already been applied to obtain accurate results for the charge dynamics within the one dimensional Holstein model.Dorfner et al. (2015); Brockt et al. (2015)
Exciton relaxation dynamics for photoexcited states of PPV are modeled in this paper using the FrenkelHolstein model and TEBD method. The phonon modes to which the exciton is coupled are assumed to be those associated with the high frequency CC bond oscillations. The low frequency torsional modes are neglected (although their possible role in relaxation and decoherence is discussed in Sec. V). We model environmental and conformational static disorder through the exciton onsite energies and hopping integrals, which leads to the exciton density associated with the low energy eigenstates of this model being spatially Anderson localized,Anderson (1958) thus defining the emissive chromophores of the polymer. External damping of the internal phonon modes is included within a Lindblad master equation approach, which allows the system to dissipate energy to the environment and thus relax into the low energy eigenstates.
This paper is structured as follows. In Sec. II, we show how a polymer chain can be represented by a coarsegrained model, generating a reduced dimensional Frenkel exciton Hilbert space in which to model the associated dynamics. In Sec. III, the FrenkelHolstein Hamiltonian is introduced for this coarsegrained model, which describes the important interactions associated with the exciton and nuclear degrees of freedom. In addition, the Lindblad master equation is introduced along with the numerical techniques used to solve the underlying equations. We then present our results in Sec. IV, including a discussion of the various processes associated with the short time relaxation dynamics, such as excitonpolaron formation, exciton decoherence and exciton density localization. We also make a connection between these processes and fluorescence depolarization measurements. Finally, we conclude in Sec. V, with a particular emphasis on comparing our results to experiment.
Ii Polymer Conformations and CoarseGraining
Polymers exist in many conformations, arising from fluctuations in the torsional angles around the single carboncarbon bonds, as well as defects in the polymer chains themselves. While in the solid state these torsional angle fluctuations are quasistatic, in solution they give rise to dynamical disorder. However, as the time scale for rotation around these single carboncarbon bonds is much longer than the associated time scale for exciton relaxation, the observables associated with the exciton relaxation dynamics can be calculated by averaging over many different static polymer chain conformations for both phases.
Figure 1 illustrates how each polymer conformation of PPV is generated. The polymer chain is built up from the starting moiety, a phenylene unit, using the parameter values given in Table 1, along with the following rules:

The torsional angles, , are taken as Gaussian random variables, with a mean and standard deviation .

With equal probability, the torsional angles, , and the vinylene bond angles, , can be positive or negative, corresponding to an anticlockwise or clockwise rotation respectively.

For a pure polymer chain, all vinylene units take the lowest energy trans geometry. However, transcis defects are introduced randomly along the chain, with probability .
In general, polymers contain many thousands of atoms and therefore performing a full atomistic simulation of exciton dynamics for each polymer conformation is unfeasible. Our model can be coarsegrained, as has been done in previous work,Tozer and Barford (2012) by representing a polymer chain as a linear array of sites, where each site corresponds to one of the moieties of the polymer. This mapping is illustrated in Fig. 1.
In conjugated polymer systems, such as PPV, the large electronelectron interactions lead to the electronhole pair of the exciton being tightly bound, typically occupying the same moiety of the chain.Barford and Paiboonvorachat (2008) This is known as a Frenkel exciton. Thus in our coarsegrained model, we only need to consider how the center of mass of the Frenkel exciton propagates along the array of sites. A single exciton state is kept per site, corresponding to the lowest energy excitation of each moiety.
Parameter  Value  Parameter  Value 

0.08 
Iii Models and Numerical Techniques
In this section we introduce the model Hamiltonians, as well as outlining the various numerical techniques used to solve the associated exciton dynamics.
iii.1 The Frenkel Model
In previous work, the initial exciton state of a PPV polymer, generated after photoexcitation, has been described by the Frenkel Hamiltonian:Barford (2013a)
(1) 
where () is the exciton creation (destruction) operator, which creates (destroys) a Frenkel exciton on moiety of the polymer chain. For PPV, the odd and even sites, , correspond to phenylene and vinylene moieties respectively. Within this Hamiltonian, the onsite energy is given by:
(2) 
where is the average moiety excitation energy, is the difference in the excitation energy between phenylene and vinylene moieties and represents the diagonal disorder, with standard deviation . This diagonal disorder arises physically from density fluctuations in the environment around the polymer chain, due to the inhomogeneity of the material.
Additionally, the nearest neighbor exciton hopping integrals present in the Frenkel Hamiltonian are given by:Barford and Trembath (2009); Barford (2013b)
(3) 
where are the torsional angles between moieties along the polymer chain. From Eq. (3), we see that there are two contributions to the nearest neighbor hopping integral. The first contribution arises from a through space dipoledipole interaction, which is incorporated through the term . While in principle this interaction term is non zero between all moieties in the polymer chain, for this analysis we only keep the nearest neighbor interactions, which are the dominant terms. The second contribution to the hopping integral arises from a through bond super exchange interaction, in which the system passes through an intermediate charge transfer exciton state. The strength of this interaction depends on the torsional angle between the moieties (, where gives the overlap of the orbitals on each moiety). Thus torsional fluctuations in the polymer chain give rise to disordered hopping integrals in the Frenkel Hamiltonian.
The disorder present in the model leads to the Frenkel Hamiltonian having two forms of eigenstates. The low energy states of the Hamiltonian are Anderson localized, nonoverlapping and nodeless states, called local exciton ground states (LEGSs).Malyshev and Malyshev (2001); Makhov and Barford (2010) A LEGS is quantified by a signedvalue parameter, , defined as:
(4) 
where is the probability amplitude of the Frenkel exciton being on site . In accord with previous work, we define a LEGS as satisfying .Malyshev and Malyshev (2001); Makhov and Barford (2010) The exciton density for three LEGSs corresponding to a particular conformation of a 99 moiety PPV polymer chain are given by the dotted curves in Fig. 2. It is the width of these states that define the size of the polymer’s absorbing chromophores.Barford and Marcus (2017) In addition, the Frenkel Hamiltonian has higher energy eigenstates that are delocalized over several chromophores, called quasiextended exciton states (QEES).Barford (2013a) The exciton density for a QEES on the same 99 moiety PPV polymer chain is given in Fig. 3.
iii.2 The FrenkelHolstein Model
Relaxation of a LEGS or QEES is physically induced by excitonnuclear coupling, which allows the nuclear geometry of the polymer to distort so that it best stabilises the newly formed excited electronic state. Most vibrational normal modes in a polymer such as PPV either couple weakly to the exciton or have a low vibrational frequency, such that they do not strongly influence the exciton dynamics on the ultrafast time scale of interest and therefore can be neglected in our model. As in previous work, we thus consider the Frenkel exciton coupling to a single high frequency normal mode per moiety (the benzenoidquinoid distortion in the phenylene unit and the CC double bond stretch in the vinylene unit) modeled using the FrenkelHolstein Hamiltonian:Holstein (1959); *holst2; Barford and Marcus (2014)
(5) 
where is the dimensionless excitonnuclear coupling parameter, which determines the strength of the excitonphonon coupling, and is the phonon energy associated with the vibrational normal modes. Within the FrenkelHolstein model, the vibrational normal mode on moiety/site is represented by a harmonic oscillator, where () is the associated dimensionless harmonic oscillator creation (destruction) operator.
One effect of the excitonnuclear coupling term in this Hamiltonian is to further localize the exciton density of the low energy LEGSs, to create socalled vibrationally relaxed states (VRSs).Barford, Marcus, and Tozer (2016) These VRSs are the low energy eigenstates of the FrenkelHolstein Hamiltonian,Tozer and Barford (2014) with their associated exciton densities given by the solid curves in Fig. 2 for a 99 moiety PPV polymer chain. Comparing these VRSs with their associated LEGSs (given by the dotted curves) shows that the degree of selflocalization arising from coupling to the high frequency vibrational modes is small, leading to the spatial extent of the VRSs still being largely determined by the disorder in the Frenkel Hamiltonian.Tozer and Barford (2014); Barford, Marcus, and Tozer (2016) It is the width of these VRSs that define the size of the polymer’s emissive chromophores.Tozer and Barford (2012)
The values for the FrenkelHolstein Hamiltonian parameters that we used to model PPVTozer and Barford (2012) are given in Table 2.
Parameter  Value  Parameter  Value 

eV  eV  
eV  eV  
meV  
eV  0.033 
iii.3 The Lindblad Master Equation
To describe the time evolution of a Frenkel exciton (initially in a high energy QEES or LEGS) relaxing to the low energy VRSs, the model must allow the system to lose excess energy to the environment. Dissipation of energy arising from systemenvironment coupling is commonly accounted for by describing the time evolution of the quantum system by a Lindblad master equation:Breuer and Petruccione (2002)
(6) 
where is the dimensionless time. Such an equation describes the time dependence of the system density matrix, , from which any observable associated with the system of interest can be obtained. From Eq. (6), we see that the time dependence of the system density matrix depends on two terms. The first term on the right of Eq. (6) corresponds to the evolution of the system density matrix under the Hamiltonian, , given by:
(7) 
where is the dimensionless displacement operator and is the dimensionless momentum operator, for the oscillator on site , and the last term in Eq. (7) accounts for the damping correction to the harmonic oscillator frequency induced by the external dissipation.Barchielli and Vacchini (2015) The second term in Eq. (6) accounts for energy dissipation of the system induced by the systemenvironment coupling. The effect of the systemenvironment coupling is controlled by the dimensionless dissipation parameter , as well as the Lindblad operator for site , , which for the system considered in this paper becomes the dimensionless harmonic oscillator destruction operator.
While the Lindblad master equation only approximately describes dissipation effects due to systemenvironment coupling, the perturbative nature of its derivation means that we expect the master equation to describe the dynamics of systems in which the systemenvironment coupling is weak (i.e., when is small). Such a situation is indeed physically relevant to the exciton relaxation dynamics in polymer chains, in which the excitonphonon coupling strength is much larger that any coupling terms to the environment.
It is not necessarily obvious why the harmonic oscillator destruction operator is a good choice for the Lindblad operator in our master equation in Eq. (6). To show that it is, we first note that the expectation value of some operator corresponding to a system observable, , can be calculated from the system density matrix as follows:
(8) 
where corresponds to taking the trace of the quantity inside the brackets. Using Eq. (6), we can then derive the following equations of motion for the dimensionless nuclear displacement and momentum:
(9) 
which have the form of Newton’s equations of motion for a damped harmonic oscillator Breuer and Petruccione (2002); Barchielli and Vacchini (2015) and are the identical equations of motion used previously for modeling exciton relaxation dynamics for classical nuclei.Tozer and Barford (2012)
iii.4 Optical Intensity
In an experiment, the initial exciton state from which the relaxation dynamics proceed is formed by photoexcitation (typically with a broadband pulse) from the electronic ground state. On application of this pulse, the nuclear degrees of freedom remain static (i.e., the harmonic oscillators remain in their ground state), while the exciton is created in a state that corresponds to a linear combination of the eigenstates of the Frenkel Hamiltonian. The probability that an eigenstate of the Frenkel Hamiltonian with energy contributes to this initial exciton state is given by the spectral function:
(10) 
where is the energy eigenvalue associated with eigenstate of the Frenkel Hamiltonian. In addition, is the associated oscillator strength for this eigenstate, given by:
(11) 
where corresponds to the electronic ground state of the system. On calculating the spectral function using Eq. (10), we average our result over the various instances of the disorder in the Frenkel Hamiltonian, as well as over all possible conformations of the polymer chain. For a particular conformation, the transition dipole moment operator, , takes the form:
(12) 
where is the magnitude of the transition dipole moment for a Frenkel exciton localized on a single moiety of the polymer and () is the exciton creation (destruction) operator for site/moiety . Additionally, is a unit vector that points along the polymer axis at moiety , which is illustrated in Fig. 1.
Figure 4 shows the calculated spectral function for PPV, using the parameter values given in Tables 1 and 2. In this work, we chose to perform our relaxation dynamics for initial Frenkel exciton energies of 2.68 eV and 2.88 eV, shown by the dotted lines in Fig. 4. Both of these Frenkel exciton energies have appreciable formation probabilities and also allow us to investigate the differing relaxation dynamics of LEGSs, low energy QEESs and high energy QEESs.
iii.5 Numerical Techniques
Section III.3 outlines the necessary theoretical framework for modeling the exciton relaxation dynamics. However, evolving the system density matrix numerically using the Lindblad master equation given in Eq. (6) is not straight forward for the following reasons. The FrenkelHolstein model contains excitonphonon interactions, which lead to its associated Hilbert space growing exponentially with the number of sites/moieties in the polymer chain. Thus for chains containing a physically realistic number of moieties/sites, this Hilbert space is too large to manipulate on a computer. This issue is compounded further, as the size of the system density matrix is the Hilbert space squared, meaning the use of standard differential equation solvers to compute the system density matrix from Eq. (6) is intractable. In this work, these issues are overcome by using the time evolving block decimation (TEBD)tnt (); Vidal (2003, 2004) and the quantum jump trajectory methods,Daley (2014); Dalibard, Castin, and Molmer (1992); Molmer, Castin, and Dalibard (1993); Dum, Zoller, and Ritsch (1992) respectively. The details of these techniques are not discussed here, but since these are unfamiliar to the chemical physics community, an overview of both can be found in the Appendices.
Before a discussion of our results, we briefly outline the basic constituents of the quantum jump trajectory method, which will prove useful in interpreting the role of the external dissipation in our calculated dynamics. The quantum jump trajectory method obtains system observables associated with the numerical solution of a Lindblad master equation, such as Eq. (6), by averaging over socalled quantum trajectories. A single trajectory is represented by a state in the system Hilbert space, with its time evolution determined as follows. Typically, the trajectory evolves under the effective Hamiltonian:Daley (2014)
(13) 
where is the system Hamiltonian given in Eq. (7). However, at probabilistically determined times, the trajectory undergoes a ‘quantum jump’, which corresponds to application of one of the Lindblad operators onto the state. Averaging over a large number of these stochastically determined trajectories is known to accurately reproduce the timedependent observables associated with the Lindblad master equation.Daley (2014)
Iv Results
The dynamics must fulfil certain criteria if the time evolution determined by our theoretical methodology is physically correct. One such requirement is that the total energy of the system must decrease during the dynamics, approaching the energies of the VRSs, which define the emissive chromophores of the polymer. We indeed find that the system energy does decrease in our dynamics model, as shown by the time evolution of the total energy of the system for two different initial Frenkel exciton energies given in Figure 5. This confirms that the form of the external dissipation in our model is having the desired effect. Also shown in this figure are the energies associated with the VRSs, for various instances of the disorder, given by the dotted curves. The fact that the total system energy at fs is greater than the energy of these VRSs for both initial Frenkel exciton energies, as well as being non zero at this time, suggests that the dynamics are far from equilibrium. This is not surprising, since the value of the exponential dissipation time scale, fs, means that the system will only reach equilibrium at times much longer than computationally practical simulations of fs.
We observe a separation of time scales within our computed dynamics. This is illustrated in Fig. 6, which shows the timedependent phonon energy for two different initial Frenkel exciton energies, with and without external dissipation. For very short times, the evolution of the system observables depend on the FrenkelHolstein parameters and are largely independent of the external dissipation, . From Fig. 6, we see that the phonon energy initially increases during the dynamics, suggesting that this ultrafast time scale is characterised by energy transfer from the exciton to the internal nuclear degrees of freedom (i.e., the phonon degrees of freedom act as a heat bath for the exciton). Eventually, the energy associated with the nuclear degrees of freedom saturates, with the time evolution of the observables now dependent on the longer external dissipation time scale. Figure 6 shows that for , the phonon energy decreases at long times, suggesting that the long time scale dynamics are characterised by dissipation of the phonon energy to the environment (i.e., the environment acts as a heat bath for the internal nuclear degrees of freedom).
Even though the system is far from equilibrium for the times considered in our relaxation dynamics simulation, we should expect the time evolution of observables to still exhibit features associated with the steady states. Indeed, if the system does relax to the low energy VRSs as physically expected, processes such as excitonpolaron formation, exciton decoherence and exciton density localization should be present in the short time dynamics. We now consider each of these processes.
iv.1 ExcitonPolaron Formation
In systems with strong electronphonon interactions, the low energy excitations of the system are described by polarons, which are quasiparticles consisting of an electron or exciton ‘dressed’ or ‘selftrapped’ by local displacements of the nuclei.Landau (1933) (To distinguish dressed electrons and excitons, we use the term ‘excitonpolaron’ in this paper for the dressed exciton.) The effective size of the excitonpolaron can be described by the excitonphonon correlation function:Hoffmann and Soos (2002)
(14) 
where the factor normalizes the function for eigenstates of the FrenkelHolstein Hamiltonian (i.e., ).Romero, Brown, and Lindenberg (2000) Physically, this correlation function corresponds to the average nuclear displacement sites away from the exciton. From previous work, it is known that for the VRSs of the FrenkelHolstein model, the excitonphonon correlation function has an exponential form, with a correlation length that is a function of the parameters , and .Tozer and Barford (2014) If the system does evolve to the VRSs, then we would expect the excitonphonon correlation function to acquire this exponential behavior during the exciton relaxation dynamics.
Figure 7 shows the time evolution of the excitonphonon correlation function for an initial high energy QEES given in Fig. 3. This correlation function acquires an exponential form during the dynamics, confirming that excitonpolaron formation is occurring within our model. Excitonpolaron formation occurs on an ultrafast time scale ( fs), with the dynamics being largely independent of the external dissipation. Indeed, we find that the time scale associated with excitonpolaron formation solely depends on the parameter , with the local excitonphonon correlations present within half a vibrational time period. In addition, persisting oscillations occur in the excitonphonon correlation function with a time period ( fs) that matches the vibrational time period of the harmonic oscillators in the model. These oscillations thus arise from the time evolution of the nuclear displacements , which physically correspond to the CC bond oscillations.
We now consider whether the correlation length associated with the dynamical excitonphonon correlation function given in Fig. 7 is the same as for the VRSs. The solid curves in Fig. 8, which correspond to the excitonphonon correlation functions at fs for two initial Frenkel exciton energies ( eV and eV), as well as several instances of the disorder, all coincide. This is to be expected, as the parameters , and have the same values for all these curves, giving rise to an identical excitonpolaron formation time scale and excitonphonon correlation length. Also shown in Fig. 8 is the average excitonphonon correlation function associated with the VRSs, given by the dashed curve. The form of the dynamical excitonphonon correlation functions at fs shows good agreement with the same function for the VRSs, which is consistent with the system evolving to these states during the dynamics. However the agreement is not perfect, with deviations arising due to the persistent oscillations in the dynamical correlation function, which are symptomatic of the time evolution having not yet reached the steady state.
iv.2 Exciton Decoherence
An exciton confined to a one dimensional polymer chain in general can have long range quantum coherences between the moieties, which physically give rise to interference effects in the calculated observables. While the range of these coherences are often limited by the presence of disorder in the system, for QEESs the exciton coherences can still persist over a distance of several chromophores. If the exciton localizes onto a single chromophore during the relaxation dynamics, then these long range coherences will decay.
One way to quantify the magnitude of the exciton coherences is from the offdiagonal elements of the exciton reduced density matrix, :
(15) 
which is obtained by taking the trace of the system density matrix over the nuclear degrees of freedom (characterised by the quantum number ). The following exciton coherence correlation function can then be defined, which gives the average magnitude of the exciton coherences between moieties/sites a distance apart:Kuhn and Sundstrom (1997); Smyth, Fassioli, and Scholes (2012)
(16) 
where is the ket corresponding to the exciton on site .
Figure 9 shows the time dependence of the exciton coherence correlation function for an initial high energy QEES given in Fig. 3. The correlation function rapidly localizes, showing that within our model, decoherence of the exciton occurs on an ultrafast time scale. This time scale is more evident from the time dependence of the exciton coherence number, :
(17) 
which corresponds to the average number of moieties over which exciton coherences persist and is given in the inset of Fig. 9. Indeed, we find that the coherence number reaches its equilibrium value within fs. As for excitonpolaron formation, the short time scale associated with this decoherence leads to the associated dynamics being largely independent of the external dissipation and depending solely on the FrenkelHolstein parameters.
The exciton decoherence mechanism present within our model can be elucidated by considering the following state representation for the FrenkelHolstein model:
(18) 
Within this representation, corresponds to the probability amplitude for the Frenkel exciton residing on site/moitey , given by , while the ket corresponds to the state of the harmonic oscillators associated with the exciton on site . Using this state representation, the exciton coherence correlation function takes the form:
(19) 
This expression contains two terms. The first term, , corresponds to the exciton wavefunction overlap between the two sites/moieties. As shown in Sec. IV.3, this term remains essentially stationary over the time scale of exciton decoherence. The second term, , corresponds to the overlap of the vibrational states associated with the exciton being on sites/moieties and . In Sec. IV.1, we saw that during the exciton relaxation dynamics, excitonpolaron formation occurs on an ultrafast time scale, leading to the nuclear displacements of the polymer becoming locally correlated to the exciton. This means that during the dynamics, the vibrational state will only have nuclear displacements spatially close to moiety/site (where the exciton resides), thus causing the vibrational overlap to decrease with increasing site separation . It is this local nature of the vibrational overlap that accounts for the local nature of the exciton coherence correlation function in Fig. 9. Excitonpolaron formation is therefore the mechanism by which the exciton decoheres, which is consistent with decoherence mechanisms found in other electronnuclear coupled systems.Bittner and Rossky (1995); Hwang and Rossky (2004); BedardHearn, Larsen, and Schwartz (2005)
Knowledge of the exciton decoherence mechanism now gives insight into the dependence of the associated time scale on the FrenkelHolstein parameters. Based on this mechanism, we would expect the exciton decoherence time scale to be identical to that for excitonpolaron formation, depending solely on and with decoherence occurring within half a vibrational time period. While this is consistent with our results, we find that the exciton decoherence time scale additionally depends on the excitonphonon coupling strength, , with the time scale decreasing with increasing . This arises, because increasing leads to an increase in the nuclear displacements associated with the excitonpolaron, which results in a more rapid decay in the vibrational overlaps that lead to decoherence.
iv.3 Exciton Density Localization
In Sec. IV.2, we saw that rapid exciton decoherence occurs during the relaxation dynamics, driven by the decay in the overlap of the vibrational states associated with excitonpolaron formation. As exciton decoherence arises from the vibrational configuration of the system, the exciton density can in principle be delocalized along the polymer chain, even if the exciton coherences between different sites/moieties are short ranged. Therefore, if the system does relax to the low energy VRSs, the exciton density must also localize during the time evolution.
In Fig. 10, the solid curve corresponds to the exciton density at time fs for an initial high energy QEES given in Fig. 3. Also shown are the scaled exciton densities for the three VRSs of the polymer chain, given by the dotted curves. Three of the peaks in the exciton density at fs match those corresponding to the VRSs, suggesting that the system is relaxing into these low energy states during the time evolution. The presence of additional peaks in the time evolved exciton density are probably a result of the system having not reached the steady state at fs, but could also signify that the steady state solutions of our Lindblad master equation do not correspond to the VRSs. As we are primarily interested in understanding the short time dynamics in this paper, we will pursue this potential issue further in a subsequent paper, where the steady state solutions of the Lindblad master equation will be investigated.
While the time evolved exciton density contains peaks corresponding to the VRSs of the polymer chain, the solid curve in Fig. 10 seems to show that the exciton density remains quasidelocalized during the dynamics. However, as the system is described by a mixed state density matrix, the associated observables physically correspond to an ensemble average over many different environment configurations (or alternatively, over many quantum trajectories). This means that it is not a priori obvious whether the nonlocal nature of the ensemble averaged exciton density is symptomatic of an actual absence of exciton density localization in our model or whether it is simply a consequence of ensemble averaging over exciton density that has localized onto different chromophores. To distinguish between these two possibilities requires a correlation function that measures the spatial extent of the exciton for a single environment configuration, such as:Spano et al. (2008); Tempelaar et al. (2014)
(20) 
where the scaling factor is chosen so that . In this definition, the ket corresponds to the exciton being on site/moiety , while all the harmonic oscillators are in their ground state, signified by the ‘0’ index.
The physical interpretation of this exciton localization correlation function can be understood by considering the state representation for the FrenkelHolstein model, given in Eq. (18). Using this state representation, the correlation function takes the form:
(21) 
where is the probability amplitude for the exciton being on moiety/site , corresponds to the state of the nuclear degrees of freedom when the exciton resides on site and corresponds to the ground state of the harmonic oscillators in the system. To a good approximation, the quantity will be independent of the exciton site index . This is because the spatial distribution of the oscillator displacements around the exciton depend on the excitonphonon correlation length, which will be largely independent of the exciton site index. Applying this approximation to Eq. (21) leads to the following simplified form for the exciton localization correlation function:
(22) 
This expression shows that the correlation function can be regarded as giving the average magnitude of the exciton wavefunction overlap between moieties/sites a distance apart and therefore gives a measure of the spatial extent of the exciton on a polymer chain.
Figure 11 shows the time dependence of the exciton localization correlation function for an initial high energy QEES given in Fig. 3. The main figure corresponds to the time evolution with the standard external dissipation parameter of , while the upper inset corresponds to the time evolution without external dissipation to the environment. When external dissipation is present, this correlation function does localize, confirming that exciton density localization is present within our model of the relaxation dynamics. However, the time scale associated with this process appears substantially different from the ultrafast excitonpolaron formation and exciton decoherence time scales studied previously. Indeed, the exciton localization correlation function remains essentially static for initial times on the order of the ultrafast time scale ( fs) with the function only starting to localize at longer times. This time scale is more evident from the time dependence of the exciton localization number, :Tempelaar et al. (2014)
(23) 
which corresponds to the average number of moieties over which the exciton wavefunction overlap remains non zero and is given in the lower inset of Fig. 11. The exciton localization number continues to decrease for fs, illustrating the much longer time scale associated with this process. This suggests that the exciton density localization is driven by the external dissipation to the environment. Further confirmation of this is given in the upper inset of Fig. 11, which shows that without external dissipation, the correlation function remains delocalized during the time evolution.
The mechanism by which external dissipation drives this exciton density localization can be understood by investigating the time evolution of individual trajectories within the quantum jump trajectory method. The quantum jump trajectory method is a way of simulating the time evolution described by a Lindblad master equation and is outlined briefly in Sec. III.5, with a more detailed summary given in Appendix A. Figure 12 shows the time dependence of the exciton density for a typical trajectory associated with this technique, where the system is initially in a high energy QEES given in Fig. 3. For fs, the dynamics correspond to time evolution under the effective Hamiltonian, given in Eq. (13), where the figure shows that during this evolution, the exciton density remains quasidelocalized along the polymer chain. The discontinuity in the exciton density at fs corresponds to a ‘quantum jump’, in which one of the harmonic oscillator destruction operators, , is randomly applied to the system. The effect of this ‘quantum jump’ process is to cause an immediate localization of the exciton density, as seen in Fig. 12. ^{1}^{1}1Therefore, it is the lack of these ‘quantum jump’ processes within the Ehrenfest dynamics approximation that lead to the unphysical bifurcation of the exciton density found in previous work.Tozer and Barford (2012) This arises because only states of the harmonic oscillator on site with a non zero displacement remain after application of the operator, . As these nuclear displacements are locally correlated to the exciton through excitonpolaron formation, application of therefore leads to a state with the exciton localized around site . Physically this process corresponds to ‘wavefunction collapse’, where local environment interactions with the internal phonon degrees of freedom act as a quantum measurement. Hence, it is the position of these local environment interactions that determines the particular chromophore that the exciton density relaxes onto during the time evolution.
iv.4 LEGSs Dynamics
In the preceding analysis, the main focus has been on understanding the short time exciton relaxation dynamics of initial QEESs. However, the features of the dynamics, as well as their associated time scales, are largely identical for initial LEGSs, with a few important differences. For QEESs, we saw that the system evolves to become a mixed state combination of the various VRSs on the polymer chain, as seen in the time evolved exciton density. In contrast, we find that an initial LEGS adiabatically relaxes almost entirely onto a single VRS, with the exciton density remaining on the same chromophore throughout the time evolution. This finding is in agreement with previous work on the relaxation dynamics of LEGSs.Barford, Boczarow, and Wharram (2011); Tozer and Barford (2012) Finally, as the spatial extent of LEGSs and VRSs are very similar (see Fig. 2) the extent of exciton decoherence and exciton density localization during the time evolution is much less pronounced compared to that for an initial QEES, resulting in a smaller time dependence of the associated correlation functions.
iv.5 Time Resolved Fluorescence Anisotropy
So far, we have outlined many features present in our model of the exciton relaxation dynamics in polymer systems, such as excitonpolaron formation, exciton decoherence and exciton density localization. We now consider whether such features can be seen in experimental observables used to study the dynamics.
One such experimental technique is the measurement of the timedependent decay in the fluorescence anisotropy after photoexcitation. It is known experimentally that as the exciton relaxes, the polarization axis associated with the fluorescence rotates from that of the incident radiation, caused by a rotation of the transition dipole moment of the polymer. This can be quantified using the fluorescence anisotropy, :Lakowicz (1983)
(24) 
where and are the intensities of the fluorescence radiation polarised parallel and perpendicular to the incident radiation, respectively.
For an arbitrary state of a quantum system, , the fluorescence intensity polarised along the axis is related to the component of the transition dipole operator, , by:
(25) 
where corresponds to the system in the ground electronic state, with the nuclear degrees of freedom in the state characterised by the quantum number . In principle, the expression for also has an energy dependent term. However, as long as the variance of the energy associated with the system remains small during the time evolution, this term just becomes a multiplicative constant that can be neglected in Eq. (25). The expression given in Eq. (25) can be generalised for a mixed state density matrix, , as follows:Rahman, Knox, and Kenkre (1979)
(26) 
where the final equation is obtained by using the expression for the transition dipole moment given in Eq. (12). In this equation, corresponds to the exciton reduced density matrix, defined in Eq. (15), while is the ket corresponding to the exciton residing on site . In addition, is the component of the unit vector that points along the polymer axis at site , illustrated in Fig. 1. Inserting the required intensity components, given in Eq. (26), into Eq. (24) allows the fluorescence anisotropy to be obtained for a specific polymer conformation. The fluorescence anisotropy can then be averaged over several different polymer conformations and initial exciton states using the following expression:Lakowicz (1983)
(27) 
where is the total fluorescence intensity and is the fluorescence anisotropy, both associated with conformation at time . The factor of 0.4 is included on the assumption that the polymers are oriented uniformly in the bulk material.Lakowicz (1983)
Figure 13 shows the calculated time dependence of the fluorescence anisotropy for two different initial Frenkel exciton energies. The mechanism for the decay of the fluorescence anisotropy can be understood by comparing the expression for the components of the fluorescence intensity given in Eq. (26) with the expression for the exciton coherence correlation function given in Eq. (16). While the fluorescence intensity components cannot be directly written in terms of the exciton coherence correlation function, since both quantities depend on the elements of the exciton reduced density matrix, the time dependence of the two quantities are related. Indeed, we find numerically that the time dependence of the fluorescence intensity components are dominated by the decay of the off diagonal elements of the exciton reduced density matrix, which also give rise to the localization of the exciton coherence correlation function. The decay in the fluorescence anisotropy can therefore be attributed to the decoherence of the exciton during the relaxation dynamics, induced by excitonpolaron formation. This is also consistent with the observations that the decay of the fluorescence anisotropy is independent of and also occurs on a time scale equal to the exciton decoherence time scale from Fig. 9.
The difference in the magnitude of the fluorescence anisotropy decay for the two initial Frenkel exciton energies can now be explained in terms of this exciton decoherence mechanism. For the lower initial Frenkel exciton energy of eV, Fig. 4 shows that a large fraction of the initial Frenkel exciton states are LEGSs. As the exciton is already localized onto a single chromophore in a LEGS, the extent of exciton decoherence during the relaxation dynamics is much reduced for an initial LEGS compared to an initial QEES. This therefore explains why the magnitude of the fluorescence anisotropy decay is smaller for the lower initial Frenkel exciton energy in Fig. 13, consistent with experimental observations.Ruseckas et al. (2005)
V Discussion and Conclusions
In this paper, we have introduced a model for exciton relaxation dynamics in polymer systems. This model is based on describing the system in terms of a parameterized FrenkelHolstein Hamiltonian, which explicitly includes excitonphonon coupling. While this Hamiltonian has been previously used to describe exciton relaxation dynamics in polymer systems,Tozer and Barford (2012) our approach differs as we treat the nuclear degrees of freedom quantum mechanically, with external dissipation of these modes included in a Lindblad master equation formalism. Our dynamics simulations are then performed numerically using the quantum jump trajectory and time evolving block decimation techniques.
We find that our initial results for the time evolution are physically sensible, with the total system energy decreasing as a function of time, confirming that the dissipation in our model is having the desired effect. The values of the parameters relevant to polymer systems naturally lead to a separation of time scales within the dynamics. The ultrafast time scale depends on the FrenkelHolstein parameters and physically corresponds to energy transfer between the exciton and nuclear degrees of freedom, while the longer time dynamics are driven by the external dissipation parameter, , corresponding to damping of these phonon modes by the environment. These time scales are also found to be relevant to several features of our dynamics model, such as excitonpolaron formation, exciton decoherence and exciton density localization, which are symptomatic of the system relaxing to the socalled vibrational relaxed states (VRSs). In general we find:

Excitonpolaron formation occurs within half a vibrational time period of the internal phonon modes (i.e., the CC bond oscillations), namely within 10 fs. We also showed that the timeevolved excitonphonon correlation function is in good agreement with that for the VRSs, which is consistent with the dynamics relaxing to these states.

Exciton decoherence is driven by the decay in the vibrational overlaps associated with excitonpolaron formation. The time scale associated with this process is similar to that of excitonpolaron formation.

The exciton density of initially prepared quasiextended states evolves to have peaks corresponding to the more localized VRSs, which is consistent with the system relaxing onto emissive chromophores. Exciton density localization, whose ultimate cause is disorder, is driven by the external dissipation through ‘wavefunction collapse’. Therefore, this process has a much longer associated time scale than for excitonpolaron formation and exciton decoherence.
We emphasize that one key result of this work is that exciton decoherence and exciton density localization, which are terms that are often used interchangeably in the literature, are shown to be different physical processes that occur on different time scales. While exciton decoherence occurs on an intrinsic ultrafast time scale of fs as a result of the vibrational overlaps associated with excitonpolaron formation, exciton density localization (or more precisely, localization of the excitonpolaron itself) Tozer and Barford (2014); Barford, Marcus, and Tozer (2016) occurs on a much longer extrinsic time scale that is driven by the dissipation to the environment through a ‘wavefunction collapse’ process.
We now attempt to interpret experimental observations in the light of these simulations. As stated in the introduction, experimental studies of the exciton relaxation dynamics in PPV observe two sub picosecond time scales ( fs and fs). Fluorescence depolarization experiments on PPV observe that the initial anisotropy value is , suggesting that an ultrafast relaxation process is occurring with a time scale shorter than the time resolution of the experiment. This is consistent with the anisotropy decay time scale in our model and therefore we assign this ultrafast process to the exciton decoherence arising from coupling to the high frequency CC bond oscillations.
We also noted, however, that our calculated fluorescence anisotropy decay is essentially independent of the dissipation time scale, , and therefore our simulation cannot account for the time scale of fs observed in fluorescence depolarization experiments. As suggested in previous work,Dykstra et al. (2009) we predict that this further decay in the anisotropy arises from coupling of the exciton to the low frequency torsional modes in the polymer, which have an oscillation period of ca. 100 fs.Cornil et al. (1997); Tretiak et al. (2002) Indeed, as these torsional modes behave classically, they would be expected to form selflocalized Landau polaronsTozer and Barford (2014); Barford, Marcus, and Tozer (2016) and therefore cause further exciton decoherence in addition to that arising from coupling to the high frequency CC bond oscillations. However, because classical normal modes give rise to a diverging excitonphonon correlation length, the exciton decoherence and exciton density localization processes will occur via different mechanisms than for the high frequency modes. As most spectroscopic quantities depend on the exciton reduced density matrix, it is likely that the two time scales seen in these experiments arise from exciton decoherence through coupling of the exciton to these two different normal modes.
The role of the low frequency torsional modes on exciton decoherence and spectroscopic measurements will be the subject of future work. Also outstanding is the question as to what are the experimental observables associated with exciton density localization.^{2}^{2}2Our simulations indicate that the fluorescence anisotropy defined by emission into the zeroth vibronic peak is qualitatively affected by exciton density localization. We intend to compute the coherent electronic 2D spectra and threepulse photonecho spectra in an attempt to address this question. Finally, we are in the process of preparing a paper that describes the steady state solutions of the Lindblad master equation introduced here, to ascertain whether our dissipation model does cause the system to relax into the low energy VRSs, as physically expected.
Vi Acknowledgements
This work was performed using the Tensor Network Theory Library, Beta Version 1.2.1 (2016), S. AlAssam, S. R. Clark, D. Jaksch and TNT Development team, www.tensornetworktheory.org. JRM would like to thank the following for financial support: the EPSRC Centre for Doctoral Training, Theory and Modelling in Chemical Sciences, under grant EP/L015722/1, as well as University College Oxford, through the Radcliffe Scholarship. The authors would like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work (http://dx.doi.org/10.5281/zenodo.22558). We also thank S. R. Clark for helpful discussions in the early stages of this work.
References
 Ruseckas et al. (2005) A. Ruseckas, P. Wood, I. D. W. Samuel, G. R. Webster, W. J. Mitchell, P. L. Burn, and V. Sundstrom, Phys. Rev. B 72, 115214 (2005).
 Grage et al. (2003) M. M. L. Grage, Y. Zaushitsyn, A. Yartsev, M. Chachisvilis, V. Sundström, and T. Pullerits, Phys. Rev. B 67, 205207 (2003).
 Dykstra et al. (2009) T. E. Dykstra, E. Hennebicq, D. Beljonne, J. Gierschner, G. Claudio, E. R. Bittner, J. Knoester, and G. D. Scholes, J. Phys. Chem. B 113, 656 (2009).
 Sperling et al. (2008) J. Sperling, A. Nemeth, P. Baum, F. Sanda, E. Riedle, H. F. Kauffmann, S. Mukamel, and F. Milota, Chem. Phys. 349, 244 (2008).
 Dykstra et al. (2005) T. E. Dykstra, V. Kovalevskij, X. Yang, and G. D. Scholes, Chem. Phys. 318, 21 (2005).
 Yang, Dykstra, and Scholes (2005) X. Yang, T. E. Dykstra, and G. D. Scholes, Phys. Rev. B 71, 045203 (2005).
 Consani et al. (2015) C. Consani, F. Koch, F. Panzer, T. Unger, A. Köhler, and T. Brixner, J. Chem. Phys. 142, 212429 (2015).
 Hwang and Scholes (2011) I. Hwang and G. D. Scholes, Chem. Mater. 23, 610 (2011).
 Barford, Bittner, and Ward (2012) W. Barford, E. R. Bittner, and A. Ward, J. Phys. Chem. A 116, 10319 (2012).
 Holstein (1959) T. Holstein, Ann. Phys. (N. Y.) 8, 325 (1959).
 (11) T. Holstein, Ann. Phys. (N. Y.) 8, 343.
 Barford and Marcus (2014) W. Barford and M. Marcus, J. Chem. Phys. 141, 164101 (2014).
 Tozer and Barford (2012) O. R. Tozer and W. Barford, J. Phys. Chem. A 116, 10310 (2012).
 Sterpone and Rossky (2008) F. Sterpone and P. J. Rossky, J. Phys. Chem. B 112, 4983 (2008).
 Tretiak et al. (2002) S. Tretiak, A. Saxena, R. L. Martin, and A. R. Bishop, Phys. Rev. Lett. 89, 097402 (2002).
 Truhlar (2007) D. G. Truhlar, “Decoherence in combined quantum mechanical and classical mechanical methods for dynamics as illustrated for non–BornâOppenheimer trajectories,” in Quantum Dynamics of Complex Molecular Systems, edited by D. A. Micha and I. Burghardt (Springer Berlin Heidelberg, 2007) pp. 227–243.
 Bittner and Rossky (1995) E. R. Bittner and P. J. Rossky, J. Chem. Phys. 103, 8130 (1995).
 Hwang and Rossky (2004) H. Hwang and P. J. Rossky, J. Phys. Chem. B 108, 6723 (2004).
 BedardHearn, Larsen, and Schwartz (2005) M. J. BedardHearn, R. E. Larsen, and B. J. Schwartz, J. Chem. Phys. 123, 234106 (2005).
 Bittner and Rossky (1997) E. R. Bittner and P. J. Rossky, J. Chem. Phys. 107, 8611 (1997).
 Prezhdo and Rossky (1997) O. V. Prezhdo and P. J. Rossky, J. Chem. Phys. 107, 825 (1997).
 Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
 Schollwöck (2011) U. Schollwöck, Ann. Phys. (N.Y.) 326, 96 (2011).
 Daley et al. (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech.: Theory Exp. , P04005 (2004).
 Vidal (2003) G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
 Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
 Dorfner et al. (2015) F. Dorfner, L. Vidmar, C. Brockt, E. Jeckelmann, and F. HeidrichMeisner, Phys. Rev. B 91, 104302 (2015).
 Brockt et al. (2015) C. Brockt, F. Dorfner, L. Vidmar, F. HeidrichMeisner, and E. Jeckelmann, Phys. Rev. B 92, 241106 (2015).
 Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
 Barford and Paiboonvorachat (2008) W. Barford and N. Paiboonvorachat, J. Chem. Phys. 129, 164716 (2008).
 Barford (2013a) W. Barford, Electronic and optical properties of conjugated polymers (OUP, 2013).
 Barford and Trembath (2009) W. Barford and D. Trembath, Phys. Rev. B 80, 165418 (2009).
 Barford (2013b) W. Barford, J. Phys. Chem. A 117, 2665 (2013b).
 Malyshev and Malyshev (2001) A. V. Malyshev and V. A. Malyshev, Phys. Rev. B 63, 195111 (2001).
 Makhov and Barford (2010) D. V. Makhov and W. Barford, Phys. Rev. B 81, 165201 (2010).
 Barford and Marcus (2017) W. Barford and M. Marcus, J. Chem. Phys. 146, 130902 (2017).
 Barford, Marcus, and Tozer (2016) W. Barford, M. Marcus, and O. R. Tozer, J. Phys. Chem. A 120, 615 (2016).
 Tozer and Barford (2014) O. R. Tozer and W. Barford, Phys. Rev. B 89, 155434 (2014).
 Breuer and Petruccione (2002) H.P. Breuer and F. Petruccione, The theory of open quantum systems (OUP, 2002).
 Barchielli and Vacchini (2015) A. Barchielli and B. Vacchini, New J. Phys. 17, 083004 (2015).
 (41) S. AlAssam, S. Clark, D. Jaksch, and TNT Development Team AlAssam, Clark, and Jaksch (2017).
 Daley (2014) A. J. Daley, Adv. Phys. 63, 77 (2014).
 Dalibard, Castin, and Molmer (1992) J. Dalibard, Y. Castin, and K. Molmer, Phys. Rev. Lett. 68, 580 (1992).
 Molmer, Castin, and Dalibard (1993) K. Molmer, Y. Castin, and J. Dalibard, J. Opt. Soc. Am. B 10, 524 (1993).
 Dum, Zoller, and Ritsch (1992) R. Dum, P. Zoller, and H. Ritsch, Phys. Rev. A 45, 4879 (1992).
 Landau (1933) L. D. Landau, Z. Phys 3, 664 (1933).
 Hoffmann and Soos (2002) M. Hoffmann and Z. G. Soos, Phys. Rev. B 66, 024305 (2002).
 Romero, Brown, and Lindenberg (2000) A. H. Romero, D. W. Brown, and K. Lindenberg, Phys. Lett. A 266, 414 (2000).
 Kuhn and Sundstrom (1997) O. Kuhn and V. Sundstrom, J. Chem. Phys. 107, 4154 (1997).
 Smyth, Fassioli, and Scholes (2012) C. Smyth, F. Fassioli, and G. D. Scholes, Philos. Trans. R. Soc. A. 370, 3728 (2012).
 Spano et al. (2008) F. C. Spano, S. C. J. Meskers, E. Hennebicq, and D. Beljonne, J. Chem. Phys. 129, 024704 (2008).
 Tempelaar et al. (2014) R. Tempelaar, F. C. Spano, J. Knoester, and T. L. C. Jansen, J. Phys. Chem. Lett. 5, 1505 (2014).
 (53) Therefore, it is the lack of these ‘quantum jump’ processes within the Ehrenfest dynamics approximation that lead to the unphysical bifurcation of the exciton density found in previous work.Tozer and Barford (2012).
 Barford, Boczarow, and Wharram (2011) W. Barford, I. Boczarow, and T. Wharram, J. Phys. Chem. A 115, 9111â (2011).
 Lakowicz (1983) J. Lakowicz, Principles of Fluorescence Spectroscopy (Plenum Press, New York, 1983).
 Rahman, Knox, and Kenkre (1979) T. S. Rahman, R. S. Knox, and V. M. Kenkre, Chem. Phys. 44, 197 (1979).
 Cornil et al. (1997) J. Cornil, D. Beljonne, C. M. Heller, I. H. Campbell, B. K. Laurich, D. L. Smith, D. D. C. Bradley, K. Mullen, and J. L. Bredas, Chem. Phys. Lett. 278, 139 (1997).
 (58) Our simulations indicate that the fluorescence anisotropy defined by emission into the zeroth vibronic peak is qualitatively affected by exciton density localization.
 AlAssam, Clark, and Jaksch (2017) S. AlAssam, S. R. Clark, and D. Jaksch, J. Stat. Mech. Theor. Exp. , 093102 (2017).
 Orus (2014) R. Orus, Ann. Phys. (N.Y.) 349, 117 (2014).
 White and Feiguin (2004) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
Appendix A The Quantum Jump Trajectory Method
The quantum jump trajectory method is a technique that allows the numerical solution of a Lindblad master equation, such as Eq. (6), to be obtained using timedependent Schrödinger equation based methods. More specifically, observables associated with the timedependent density matrix solution of a Lindblad master equation are obtained by averaging over many socalled quantum jump trajectories. The time evolution of a single quantum jump trajectory is obtained using the following procedure:

Step 1: We first spilt up the total time evolution of a single trajectory into several time steps of length . We next define an effective Hamiltonian for the system as follows:Daley (2014)
(28) where is the system Hamiltonian, is the Lindblad master equation dissipation parameter and is the Lindblad operator for site . Evolving this single quantum trajectory, initially in state , under the effective Hamiltonian gives us a trial state at time :
(29) 
Step 2: As the effective Hamiltonian in Eq. (28) is nonHermitian, the time evolution in Eq. (29) does not conserve the quantum state’s norm:
(30) where is the amount the norm has decayed over the time step . Using Eq. (30), the state of the single quantum trajectory at time is determined probabilistically as follows:

With probability , the state at time is:
(31) 
With probability , the state at time is:
(32)
where is the Lindblad operator for site . The site at which the Lindblad operator is applied in Eq. (32) is again determined probabilistically, where the probability associated with applying the Lindblad operator at site is given by:
(33) 
Hence, we now need a technique by which to apply the time step evolution operator, , for the system of interest. This is achieved for the FrenkelHolstein model in this work using the TEBD technique, which is introduced in Appendix B.
Appendix B Time Evolving Block Decimation
b.1 Matrix Product States
It has been found that a practical quantum state representation for onedimensional systems involves using matrix product states (MPS). The MPS representation of a generic quantum state is as follows:Schollwöck (2011)
(34) 
where are known as the physical indices, which correspond to the quantum numbers of the site basis, , and are known as the internal indices. For systems made up of sites, with a physical index dimension , an arbitrary quantum state can be exactly represented by a MPS with an internal index dimension of . However, most states can be accurately represented with a much smaller internal index dimension, where increased internal index dimensions are required to represent quantum states with a greater amount of entanglement between sites.
MPSs can be helpfully represented and manipulated using tensor network diagrams.Orus (2014) Figure 14 shows the tensor network diagram for a three site MPS. In these diagrams, the circles represent the tensors in the MPS, , while the lines correspond to their indices. If a line is ‘closed’, then the index it represents is implicitly summed over.
b.2 Trotter Decomposition
In order to obtain the quantum dynamics for a system, the evolution operator for a time step, , must first be computed. For systems that have a large Hilbert space, this is not possible. Within the TEBD technique, this issue is circumvented by using the Trotter decomposition to expand the evolution operator into several smaller and more manageable terms.
For Hamiltonians that contain solely onsite and nearest neighbor terms, it can be shown that the full Hamiltonian can be written as the following sum:
(35) 
where is the Hamiltonian for the “bond” linking sites and and is the total number of sites in the system. Using Eq. (35), the Trotter decomposition for the evolution operator is given by:White and Feiguin (2004)
(36) 
The Trotter decomposition given in Eq. (36) is not exact, as the Hamiltonian terms in general do not commute with each other. However, the error associated with this Trotter decomposition is of , which becomes negligible if the time steps are chosen to be small.
b.3 The Procedure
The procedure for applying a single Trotter decomposition evolution operator onto a MPS is illustrated in Fig. 15, with the first tensor network diagram in this figure corresponding to the object . This figure illustrates the following steps:

Step 2: A singular value decomposition is applied on the tensor to obtain three new tensors, as well as a new internal index :Daley et al. (2004)
(39) If the size of the original internal index was , then the new internal index has grown to size on completion of the procedure outlined above, where is the size of the site basis. To keep the internal indices constant in size, we truncate the index by discarding the ‘states’ associated with the smallest values of .Daley et al. (2004) The error associated with this truncation procedure can be quantified by computing the sum of the discarded singular values, , which gives a bound on the error associated with the time evolved quantum state. It is this truncation procedure that keeps the total Hilbert space of the system finite and manageable throughout the dynamics.

Step 3: The resulting tensor network diagram is transformed back into MPS form by creating the new tensor :
(40)
Hence to complete a dynamics time step within the TEBD algorithm, the same three step procedure outlined above is performed for every evolution operator in the Trotter decomposition, given in Eq. (36).
Outlined above is the necessary framework for numerically obtaining the dynamics associated with the Lindblad master equation given in Eq. (6). For our simulations of the exciton relaxation dynamics for PPV polymer chains containing 99 moieties/sites, the parameters used in the TEBD algorithm are given in Table 3. These were obtained as follows. The maximum number of phonons per site, , was chosen as it gives physically sensible results for the dynamics, while keeping the size of the site Hilbert space manageable. The simulations were also performed for , where we found that the results did not deviate significantly from the calculations. Using dimensionless time, , the value of the time step, , was obtained by comparing the TEBD results for a four site FrenkelHolstein model with the numerically exact dynamics obtained from exact diagonalization. As for the case of determining , we chose to use as a balance between obtaining accurate results for the quantum dynamics, as well as making the calculation not too computationally intensive. The error associated with this Hilbert space truncation (quantified from the sum of the discarded singular values, , obtained during the singular value decompositions within the TEBD algorithm) was found to be small for each time step during our dynamics simulations, suggesting that we can be confident that the Hilbert space truncation does not lead to a large error in our results for these parameter values.
Parameter  Value  Parameter  Value 

2 