Pedro J Castro

Aurelie Perveaux

David Lauvergnat

Mar Reguero

Ultrafast internal conversion in 4aminobenzonitrile occurs sequentially along the seam

Upon UV light absorption, 4-aminobenzonitrile undergoes an ultrafast radiationless decay process from a charge-transfer state to a locallyexcited state. This pathway proceeds through an extended seam of conical intersections between the second-and first-excited singlet electronic states.

Quantum dynamics simulations show that planar geometries dominate the earlier times (< 20 fs), after which the whole seam becomes explored, as the wavepacket spreads and breathes along both sides of the unstable ridge with respect to the amino-bending coordinate.

Introduction

4-aminobenzonitrile (ABN), and its N-dimethyl derivative, DMABN [1] are well-known prototype systems for studying intramolecular charge transfer (ICT) in electron donor/acceptor molecules. Despite their chemical similarity, their luminescent patterns are quite different according to the environment: ABN shows single fluorescence in any solvent (normal band), whereas DMABN exhibits normal fluorescence in non-polar solvents but dual fluorescence in polar ones (normal and anomalous bands) [ 2 , 3 ]. The sensitivity of ABN derivatives to their environment makes them good candidates as fluorescent probes, chemical sensors, molecular switches, or electro-optical switches [4,[START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF][START_REF] Bedoya | [END_REF]7,8].

It is well-established that the normal band is originated from a locallyexcited (LE) state, while the anomalous band is due to a charge-transfer (CT) state, with a large dipole moment that is stabilized in polar solvent environments. The exact structure of the species responsible for the anomalous band and the mechanism that populates them are still controversial due to contradictory arguments, based on both experimental observations and theoretical calculations, which support different models and hypotheses [1,3,9,10,11].

Three received models have been proposed regarding the structure of the luminescent CT species: the planar ICT model (PICT) [12], the twisted one (TICT) [10,13], and the partially twisted or pretwisted one (pTICT) [14]. It seems sufficiently proven that the PICT and TICT species actually correspond to two minimum-energy points in the potential-energy surface (PES) of the excited CT state in both ABN and DMABN, but their respective roles in the fluorescent spectra remain to be clarified. A fourth model, the rehybridized ICT (RICT) [15], is thought to be a stable species corresponding to a πσ* excited state but is nowadays discarded as a luminescent species [16,17,18,19].

The ICT mechanism is triggered upon UV light absorption, which populates first the S2 state of CT character [20]. The subsequent sequence of steps leading to the S1 PES (with both LE and CT regions) until the luminescent species are formed is still under discussion. Gómez et al. proposed the following mechanism in 2005 [20]. After excitation to the S2 state, the system relaxes to the PICT minimum. The system follows then a barrierless pathway along the amino-bending coordinate to the optimized S2/S1 conical intersection (denoted CI-Cs). Because of the extended conical intersection seam, ultrafast radiationless decay from S2 to S1 can take place at various torsion angles of the amino group leading either to the LE or the TICT geometries on S1. Experimental works suggest initial production of the LE species on S1 through internal conversion from S2 and later equilibration with a CT species on S1 [9]. Park et al. reported recently an experimental study based on highly timeresolved fluorescence spectra (TRFS) over the whole emission domain of DMABN in acetonitrile [18]. Such accurate measurements provide valuable information about the dynamics of the ICT process free from the interferences of the solvent reorganization and vibrational relaxation dynamics that occur on the same time scale. The experimental techniques used in that work also give access to an analysis of the ultrafast events that occur within the first few femtoseconds. From their observations it is concluded that, after photoexcitation to the S2 state (CT), there is significant internal conversion to the S1 state in less than 30 fs for geometries that are still quasiplanar, thus in the LE region. A similar ultrafast mechanism was also proposed in a previous work by the McCamant group [21,22].

Computational studies carried out on ABN and DMABN have shown that the topographies of the S2 and S1 PESs of both systems are similar, but small energetic changes are suspected to induce a different interplay between the LE and CT species and lead to either normal or dual fluorescence [3,20]. The initial steps of the deactivation mechanism, though, are expected to be analogous in both systems, so the experimental conclusions, derived by Park et al. [18] for DMABN, should hold in a qualitative way for ABN as well.

In a previous work [23] on ABN, we have shown that ultrafast production of the LE species in S1 can occur efficiently via S2/S1 internal conversion involving a newly-investigated planar conical intersection (denoted CI-C2v), both in the gas phase and a polar solvent. Although this conical intersection is not the lowest-energy one within the S2/S1 intersection seam, it is easily accessible from the Franck-Condon point both geometrically and energetically. In other words, the system is able to deactivate first upon involving coordinates associated to such in-plane skeletal deformations in a first stage, while bending and twisting motions may occur on longer time scales once vibrational energy is redistributed. Curchod et al. did recently a study on DMABN [START_REF] Curchod | Ab Initio Multiple Spawning Photochemical Dynamics of DMABN Using GPUs[END_REF] using time-dependent DFT combined with Multiple-Spawning dynamics. While their molecular system (DMABN) is slightly different from ours (ABN), they have reached similar conclusions: ultrafast transfer from S2 to S1 can occur without torsion of the amino group.

In the present work, we investigate in more detail the competition between the deactivation pathways through a planar conical intersection, CI-C2v, and through the low-lying Cs one, CI-Cs, so as to provide further information about the ultrafast process that populates the LE state almost immediately after the initial photoexcitation of the CT state. Our objective is to determine whether deactivation occurs specifically around CI-C2v, if the system rather bifurcates early on towards CI-Cs, as suggested by Gómez et al. [20], or if the situation is intermediate, with significant delocalization along the aminobending coordinate so that the wavepacket will reach the whole portion of the S1/S2 seam connecting the bent conical intersection to the planar one at similar times.

Computational Details

The topographies of both coupled PESs involved in the initial steps of the process was first determined from a CASSCF (12,11)/cc-pVDZ study using the Gaussian 09 package [START_REF] Frisch | Gaussian 09[END_REF]. Energies in these regions were further refined at the CASPT2 level with the Molcas 7 package [START_REF] Aquilante | Software News and Update. MOLCAS 7: The Next Generation[END_REF] so as to include the effect of dynamic correlation in the calculations. In our previous study, we showed that the effect of a polar solvent on ABN only enhances the efficiency of the radiationless decay but does not change the mechanism of the ultrafast process under study (< 30fs) [3,23]. This was not accounted for in the present study, as calculations in vacuo seem sufficient to provide reliable mechanistic information.

In order to run quantum dynamics calculation with the multilayer version [ 27 ] of the multiconfiguration time-dependent Hartree (ML-MCTDH) method implemented in the Heidelberg MCTDH package [28,29], analytical expressions of quasidiabatic coupled PESs are required, often in the form of what is called a vibronic-coupling Hamiltonian (VCH) model [30,31,32]. Its entries are expressed usually as linear (linear vibronic-coupling model -LVC) or quadratic (quadratic vibronic-coupling model -QVC) functions of normal Cartesian coordinates originated from the ground-state equilibrium geometry (Franck-Condon point). This is the type of approach that we used in the present work. However, the main originality of our strategy is that we explicitly used analytical relationships between adiabatic data and diabatic parameters so as to obtain the latter automatically. Consistently with our previous study [23], three diabatic states are involved (corresponding diagonal entries of the electronic-energy matrix: 𝐻 !! (𝐐) , 𝐻 "" (𝐐) , and 𝐻 ## (𝐐)). In addition, along some directions or coordinates, the model was adjusted to obtain a better description of the conical intersections. We thus circumvent numerical fitting procedures that, in some cases, can occur to be time consuming and tedious tasks from a technical perspective. Our fulldimensional models of coupled PESs (39D in the present study) are based on CASSCF energies and geometries and further refined to match CASPT2 data in the gas phase.

The choice of an adequate coordinate system depends on the process under study. In particular, for molecular systems with large-amplitude motions, normal-mode coordinates are not adequate to describe motions leading far from the equilibrium position [33]. Thus, it is often beneficial to describe the molecular system with curvilinear coordinates, i.e., distances and angles, since they describe large-amplitude motions such as for example torsions in a more natural way; in other words, they will provide a simpler expression of the PES. Unfortunately, the use of curvilinear coordinates can lead to very complicated expressions of the kinetic-energy operator (KEO) [ 34 , 35 ], which can be expressed numerically (but exactly) [36,37,38,39,40] or analytically [41,42,43,44,45]. An analytical approach is more practical, as there is no need to compute the numerical the KEO on a grid and then fit the results or make further approximations (for example upon considering Taylor expansions). However, an analytical expression of the KEO is not always compatible with an "MCTDH format" (see below), where operators must be written as sums of products of low-dimensional functions. Some specific types of coordinates allow this condition to be fulfilled, in particular so-called polyspherical coordinates, denoted QPoly [41,43,44], which were indirectly used in this study. However, the actual coordinates, Q, used here to define the model and to run the dynamics were obtained after several transformations of the polyspherical coordinates.

Sets of Coordinates

Fig. 1 depicts the polyspherical vectors and subsystems used to define relative frames (red: S1; green: S1,1; blue: S2,1; pink: S1,2,1). Subsystems are defined so as to avoid numerical issues related to poles due to collinear vectors (such as 𝐑 ! (!) and 𝐑 " (!,!) ) around the reference geometry [45]. Each vector is expressed in terms of spherical coordinates. These coordinates, Qpoly, were carefully selected so as to correspond to the most relevant motions involved in the photoprocess. In particular, 𝑅 ! (!) describes

N 15 C 14 C 4 C 5 C 6 C 1 C 2 C 3 N 11 H 12 H 13 H 9 H 8 H 7 H 10 G 2 G 1 paper% S 1 S BF R 1 (1) R 1 (1,1) R 2 (1,1) R 3 (1,1) R 4 (1,1) R 5 (1,1) R 6 (1,1) S 1,1 S 2,1 R 4 (2,1) S 1,2,1 R 2 (2,1) R 3 (2,1) R 1 (2,1) R 1 (1,2,1) R 2 (1,2,1) R 3 (1,2,1)
the quinoidal stretching, 𝑅 " (",!) the symmetric stretching of the allyl-like fragment on the amino side, 𝛾 (",!) the envelope ring-puckering on the amino side, 𝛼 (!,",!) the pyramidalization of the C-atom attached to the amino group, 𝛾 (!,",!) and 𝜑 # (!,",!) the torsion and pyramidalization of the amino group.

Then, two transformations are required: (i) the polyspherical coordinates, Qpoly, are symmetrized with respect to the C2v point group (see

𝐑 !! LE = 𝐐 $%& '()' !" -𝐐 $%& *+ %𝐐 $%& '()' !" -𝐐 $%& *+ % , 𝐑 ,, PICT = 𝐐 $%& '()' !" -𝐐 $%& 1('2 %𝐐 $%& '()' !" -𝐐 $%& 1('2 % ,
Eq. 1

where 𝐐 -./ 01 and 𝐐 -./ 2345 are the PICT or the LE minima, respectively, and

𝐐 -./ 4364 !"
is the geometry of C2v conical intersection. Note that our procedure can define the gradient of the electronic coupling as any linear combination of the two vectors of the branching space rotated together such that 𝐻 "# (𝐐 < ) = 0 where 𝐐 < is some chosen geometry where the states are not degenerate. This point is used as a reference for setting the value of the arbitrary mixing angle between both degenerate states at the conical intersection. The extra condition thus implies that the adiabatic and quasidiabatic representations coincide at 𝐐 < . Here, symmetry guarantees that 𝐻 "# (𝐐) = 0 at all C2v/Cs points (i.e., as long as there is no distortion along B2/A" modes). For example, 𝐻 "# (𝐐 2345 ) = 0. In practice, the B2/A" derivative coupling is identified to the gradient of the electronic coupling in the present situation.

As a starting approximation, each diagonal entry was expanded as a quadratic form, 𝐻 == (𝐐) ≈ 𝐻 == quad (𝐐) with 𝑖 = 1,2,3 , around a reference

Gradient Difference Derivative Coupling

Gradient Difference Derivative Coupling geometry, 𝐐 BCD,== , corresponding to the relevant C2V stationary point: 𝐐 EF = 𝐐 BCD,!! , 𝐐 01 = 𝐐 BCD,"" , and 𝐐 2345 = 𝐐 BCD,## . The quasidiabatic curvatures of the diagonal entries, 𝐻 == quad (𝐐), were obtained from the adiabatic ones at those points through a second-order Jahn-Teller procedure (whereby the difference between the adiabatic and quasidiabatic second derivatives is considered as a second-order effect involving the non-adiabatic coupling).

In a second stage, the curvatures along the two directions 𝐐 01 -𝐐 CI-C 2v (proportional to 𝐑 "" LE , see Eq. 1) and 𝐐 2345 -𝐐 CI-C 2v (proportional to 𝐑 ## PICT , see Eq. 1) were slightly adjusted to enforce that the two quasidiabatic PESs 𝐻 "" (𝐐) and 𝐻 ## (𝐐) crossed exactly at 𝐐 CI-C 2v and did so with the right energy. Fig. 3 shows the agreement between the ab-initio energies and those of our VCH model (based on CASPT2 data in the gas phase) for a cut along the direction 𝐐 2345 -𝐐 CI-C 2v . In addition, accessing the bent conical intersection (CI-Cs) requires breaking the C2v symmetry of the molecule along the non-totally-symmetric bending motion of the amino group, 𝛼 (!,",!) . This coordinate, the first Euler angle of the sub-system (1,2,1), describes the motion of the amino group (see Fig. 1).

In other words, there is a pair of equivalent bent conical intersections, denoted CI-C 8 G and CI-C 8 6 , where the amino group is either up or down with respect to the molecular plane, as shown on Fig. 4. To account for the instability of the C2v pathway (see Fig. 5), we refined the

𝐹 switch (𝑥) = 1 + tanh T𝑥 - 𝑥 CI-4 % & 2 U 2 + 1 -tanh T𝑥 + 𝑥 CI-4 % & 2 U 2 .
Eq. 3 Fig. 5 shows the agreement between the ab-initio energies and those of our VCH model (based on CASPT2 data in the gas phase) along the aminobending direction. The largest difference between the ab-initio data and the model is 0.04 eV, which is within the range of error for high-level ab-initio methods. P) , where both the coefficients and the basis functions are optimized (as in an MCSCF electronic wavefunction),

𝜓 ABCDEFG 5Q 1 , ⋯ Q 𝑀 , 𝑡7 = 8 … 8 𝐴 H # ,…,H $ (𝑡) = 𝜑 H % (J) 5Q 𝐾 , 𝑡7 K JL: M $ H $ L: M # H # L:
. Eq. 4

The SPFs are themselves expanded in terms of primitive basis functions,

𝜑 H % (J) 5Q 𝐾 , 𝑡7 = 8 𝐶 H % ;O % (J) (𝑡) P % O % L: 𝜒 O % (J) 5Q 𝐾 7.
Eq. 5

Therefore, MCTDH can be seen as a two-layer scheme with time-dependent coefficients: 𝐴 O ) ,…,O * (𝑡) at the top layer, and sets of second-layer timedependent coefficients 𝐶 O ( ;7 ( (P) (𝑡) for each degree of freedom.

The size of the SPF basis can be further reduced by combining the physical

coordinates [Q ! ⋯ Q N ] into logical coordinates (also called combined modes) [𝐐 ! ! ⋯ 𝐐 V ! ],
such that each logical coordinate, 𝐐 P ! , comprises one or several of the physical coordinates, as

𝐐 P ! = 1Q ! ( ⋯ Q W ( 2.
The superscript 1 in the notation represents the layer number of the combined modes (notation introduced to make explicit the multilayer formulation).

The MCTDH nuclear wavefunction with combined modes reads

𝜓 ABCDEFG A𝐐 1 1 , ⋯ 𝐐 𝑃 1 , 𝑡B = 8 … 8 𝐴 :;H # ,…,H & : (𝑡) = 𝜑 H % (:;J) A𝐐 𝐾 1 , 𝑡B Q JL: M & H & L: M # H # L:
. Eq. 6

The time-dependent basis functions 𝜑 O ( (!;P) are now multidimensional.

Introducing mode combination implies that the computational effort is transferred from the propagation of a large vector of

𝐴 !;O ) ,…,O + ! (𝑡) coefficients
with one-dimensional SPFs, to a shorter vector of coefficients but multidimensional SPFs. Some experience and knowledge of the system under study is required to find an efficient mode-combination scheme for the study.

For example, combining modes with similar frequencies is a possible strategy, as shown by O. Vendrell et al. [27].

The mode-combined SPFs expressed in the primitive basis are given by

𝜑 H % (:;J) A𝐐 𝐾 1 , 𝑡B = 8 … 8 𝐶 H % ;O # % …O ' % (!;J) (𝑡) = 𝜒 O ( % (J,R % ) (Q 𝑖 𝐾 ) S % RL: P ' % O ' % P # % O # % .
Eq. 7

In a high-dimensional system one should combine groups of degrees of freedom into high-dimensional SPFs in order to make the size of the vector of coefficients in Eq. 6 manageable (i.e., to get a wavefunction propagation that is reasonable in terms of computation time). However, the combined SPFs are too large to be efficiently propagated. The ML-MCTDH layering scheme is a very flexible way of dealing with this issue. One treats the combined mode as a "sub-configuration" involving smaller groups of logical coordinates. This introduces a new layer of coefficients, the size of which is manageable. The procedure can be repeated over and over until the primitive degrees of freedom are reached.

According to standard use, we may refer to the one-layer scheme as the standard method (primitive basis), to the two-layer scheme simply as MCTDH, and to deeper layering schemes as ML-MCTDH.

Our system, consisting in 39 nuclear coordinates, was described with up to a five-layer wavefunction, and for each layer eight SPFs were used except for the last layer. For the primitive basis sets (the so-called last layer) we used, as usual, the discrete variable representation (DVR) with 40 grid points corresponding to harmonic-oscillator basis sets. The vibrational ground state on S0 was obtained with the relaxation procedure implemented for the ML-MCTDH method.

As shown in the previous section, to adjust the linear vibronic-coupling model, the coordinates used in the dynamics are defined as linear combinations of polyspherical coordinates. Therefore, the analytical KEO is not simply expressed as a sum of products and it cannot be used with MCTDH directly. To overcome this difficulty and with the help of TNUM [38], we generated a Taylor expansion of the metric tensor, G(Q), around the reference geometry of S0. In this study, we only used a zero-order approximation. In other words, the G-matrix will be considered constant all over the coordinate grid.

Results and Discussion

Fig. 6 shows a comparison of the evolution of the quasidiabatic populations obtained with: (i) our previous model [23], essentially focused on the description of the deactivation pathway through the C2v conical intersection;

(ii) the present model, which also describes correctly the C2v-breaking motion along the amino-bending direction that leads to the Cs conical intersection.

It shall be noted that the generalised VHC model and the quantum simulations in the present study enable to describe only the short-time (subpicosecond) internal-conversion process. In addition, in order for the system to return to the low-lying electronic state and thus reach thermal equilibrium, other processes must be taken into account, such as other internal-conversion mechanisms, internal vibrational redistribution, fluorescence, solvent effect... One can notice that the transfer of quasidiabatic population between the CT and LE electronic states is faster when both deactivation pathways are described on the same footing (plain line) but it becomes similar on longer time scales (after 50-100 fs).

However, this result is global and does not discriminate the relative efficiencies and roles of the two different conical intersections (i.e., CI-Cs and CI-C2v) with respect to the deactivation process. To this end, we defined different regions along the amino-bending coordinate, centred on each conical intersection and with borders located midway between CI-C2v and CI-Cs +/-. The partial populations integrated within these regions will be called local quasidiabatic populations (for each specific conical intersection) and their time evolutions are presented in Fig. 7 for the LE state. The evolution of the local quasidiabatic populations within the inner (around CI-C2v) and outer (around CI-Cs +/-) regions shows two sequences. The first sequence concerns times less than 20fs: the transfer of population occurs mainly in the vicinity of CI-C2v (sharp peak corresponding to a maximum of 24% of population transferred to LE around CI-C2v whereas population transferred around both CI-Cs + and CI-Cs -increases more smoothly, up to 16% within the same period of time; see Fig. 7). After 20 fs, the local quasidiabatic population within the outer region becomes larger and evolves as a mirror image of the local quasidiabatic population within the inner region. One can notice that the global quasidiabatic population no longer evolves after about 25 fs (green plain line in Fig. 6). This indicates that the width of the wavepacket has now increased enough for both regions to be populated significantly and only exhibits a slight oscillatory breathing along the amino-bending direction. In summary, the non-adiabatic transfer of population occurs sequentially along the seam. Before 20 fs, the wavepacket travels on a ridge along the planar deactivation channel with no significant spreading. After 20 fs, delocalization along the amino-bending coordinate implies that the whole seam is reached. Comparison with the results our previous simulations shows that accounting for this spreading makes the transfer of population a bit faster but not more efficient overall.

Conclusions

In this study, we have performed quantum dynamics simulations on ABN in the gas phase in order to describe the competition between two deactivation pathways from the initially-excited S2 state to the luminescent S1 state and provide further information about the ultrafast process that populates the LE state almost immediately after the initial photoexcitation.

We have developed a full-dimensional (39D) quasidiabatic model of coupled PESs based on CASSCF (12,11)/cc-pVDZ data further refined to match to CASPT2 data in the gas phase. This model, which describes the whole portion of the S2/S1 deactivation seam along the amino-bending coordinate, was used to run quantum dynamics simulations with ML-MCTDH.

Our quantum dynamics calculations clearly indicate that deactivation is first localized around the planar conical intersection (CI-C2v) and becomes delocalized after about 20 fs along the S2/S1 seam connecting the pair of lower-energy amino-bent conical intersections (CI-Cs + and CI-Cs -) to the planar one (centre of the seam). The PICT species on S2 is populated right after photoexcitation where it undergoes internal conversion to the LE species on S1 very early, mainly through the planar conical intersection. The amino-bending motion is not crucial for the overall efficiency of the internal conversion. However, delocalization of the wavepacket along the aminobending direction implies that the wavepacket explores a larger portion of the seam, which makes the deactivation process faster. Such computational results are consistent with the recent theoretical work of Curchod et al. on the parent system DMABN [START_REF] Curchod | Ab Initio Multiple Spawning Photochemical Dynamics of DMABN Using GPUs[END_REF] and the recent experimental observation of Park et al. [18] in DMABN as well where the first deactivation step is attributed to quasiplanar motions inducing ultrafast internal conversion and direct formation of the LE species in less than 30 fs.

As shown in our previous study [23], this feature is expected to be even more efficient in polar solvents than in the gas phase due to the strong stabilization of the CT state with respect to the LE state, which results in making the conical intersections closer to the Franck-Condon region.

Fig. 1 .

 1 Fig. 1. Polyspherical vectors and subsystems.

2 2. 2

 22 !,:) + ∆𝜃 ! (:,!,:) Quasidiabatic Model of Coupled PESsThe quasidiabatic PESs and couplings are represented as a 3×3 electronic Hamiltonian matrix based on a generalization of the VCH model[30]. It consists in a real symmetric matrix, 𝐇(𝐐), made of three diagonal potentialenergy functions, 𝐻 !! (𝐐), 𝐻 "" (𝐐), and 𝐻 ## (𝐐) , and three off-diagonal electronic couplings, 𝐻 !" (𝐐), 𝐻 !# (𝐐), and 𝐻 "# (𝐐) . In the Franck-Condon region the three diabatic states are chosen to correlate with the adiabatic ones: state 1 (S0/GS), state 2 (S1/LE), and state 3 (S2/CT).The non-adiabatic coupling terms between S0 and the other two states can be neglected. In other words, 𝐻 !! (𝐐) is identified to the S0 PES and the electronic couplings 𝐻 !" (𝐐) and 𝐻 !# (𝐐) are set to zero. The remaining coupling, 𝐻 "# (𝐐), is expanded linearly around 𝐐 CI-C 2v , the geometry of the C2v conical intersection between S1 and S2, along the derivative coupling calculated at this point. The derivative couplings at both the CI-C2v and the CI-Cs points are very similar and correspond essentially to a B2/A" Kekulélike motion that breaks the left/right symmetry of the molecule, and which does not change much all over the seam as illustrated on Fig.2.

Fig. 2 .

 2 Fig. 2. Gradient differences and derivative couplings at CI-C2V (upper panel) and CI-Cs (lower panel) calculated at the CASSCF(12,11)/cc-pVDZ level of theory with 0.5:0.5 stateaveraging on S1 and S2.

Fig. 3 .

 3 Fig. 3. C2v pathway along a linear-interpolation coordinate Q linking the PICT point (Q = 0) to the CI-C2v point (Q = 13). Energy differences are given with respect to the groundstate minimum. Dashed lines: ab initio; plain lines: model. Ab-initio level of theory: CASPT2/cc-pVDZ and MS2-CASPT2 at the CI-C2v point.

Fig. 4 .

 4 Fig. 4. Optimized geometry at the CASSCF(12,11)/cc-pVDZ level of theory of the pair of S2/S1 conical intersections denoted CI-𝐂 𝒔 " (a) and CI-𝐂 𝒔 # (b).

Fig. 5 .

 5 Fig.5. Cs pathway along a linear-interpolation coordinate Q linking the PICT point (Q = 0) to the pair of CI-Cs + and CI-Cs -points (Q = +8 and -8). Energy differences are given with respect to the ground-state minimum. Dashed lines: ab initio; plan lines: model. Ab-initio level of theory: CASPT2/cc-pVDZ and MS2-CASPT2 at the CI-Cs +/-points.

Fig. 6 .

 6 Fig. 6. Time evolution of the quasidiabatic populations of the LE and CT electronic states in the gas phase. Coupled PES model based on CASPT2 data. Plain line: model describing both bent and planar deactivation pathways; dashed line: model describing the planar deactivation pathway [23].

Fig. 7 .

 7 Fig. 7. Time evolution of the local quasidiabatic populations transferred to the LE state in the gas phase. Coupled PES model based on CASPT2 data. Plain line: local population

Table 1

 1

); such symmetrized coordinates, Qsym, make it possible to enforce C2V symmetry into ab-initio data (Hessians, gradients) obtained numerically; (ii) effective coordinates are defined from a linear transformation of Qsym such that a single coordinate is defined according to each conical-intersection direction (Eq. 1). The remaining linear combinations belong to the orthogonal complement. This coordinate change is performed automatically with the TNUM program that uses the vectors 𝐑 "" LE and 𝐑 ## PICT expressed in terms of the symmetrized coordinates,

Table 1 .

 1 Variations of the original set of polyspherical coordinates (∆𝐐 𝐏𝐨𝐥𝐲 = 𝐐 𝐏𝐨𝐥𝐲 -𝐐 𝐏𝐨𝐥𝐲 𝟎 ) symmetrized with respect to a planar C2v reference geometry (𝐐 𝐏𝐨𝐥𝐲 = 𝐐 𝐏𝐨𝐥𝐲 𝟎 ). In this, we used a convention where (zy) is the molecular plane (||) and (zx) is the perpendicular plane (⊥). 𝑨 𝟐 deformations are torsions, 𝑩 𝟐 deformations are in-plane left-right-breaking motions, and 𝑩 𝟏 deformations are out-of-plane up-down-breaking motions.

	𝐴 ! [𝐶 "7 ⇒ 𝐶 "7 ]	𝐴 " [𝐶 "7 ⇒ 𝐶 " ]	𝐵 ! [𝐶 "7 ⇒ 𝐶 8 9 ]	𝐵 " 1𝐶 "7 ⇒ 𝐶 8 ∥ 2
	∆𝑅 : (:)

  expression of 𝐻 == (𝐐) (see Eq. 2) for 𝑖 = 2,3 upon using a switch-type function denoted 𝐹 The value of 𝑓 C == is obtained at 𝐐 CI-4 % & and thus at 𝑥 CI-4 % & (or 𝐐 CI-4 % ' and 𝑥 CI-4 % ' )

			(a)	(b)
	intersection (CI-C F G and CI-C F 6 ) energies, read:
	𝐻 == (𝐐) = 𝐻 == quad (𝐐)		Eq. 2
	+	1 2	𝑓 C == G𝐹 switch (𝑥) -𝐹 switch (0)I J	tanh (0.5 • 𝑥) 0.5 P
	upon imposing 𝐻 == Q𝐐 CI-4 % & R = 𝑉 Q𝐐 CI-4 % & R in Eq. 2, where 𝑉 Q𝐐 CI-4 % & R is the
	ab-initio energy of the conical intersection at 𝐐 CI-4 % & ,

switch (see Eq. 3). Its purpose is to modify symmetrically the curvature with respect to the amino-bending coordinate, x (x is one of the coordinates Q), and enforce that 𝐻 "" (𝐐) and 𝐻 ## (𝐐) intersect exactly at 𝐐 CI-4 % & and 𝐐 CI-4 % ' with the right energy. The switch-modified diagonal quasidiabatic PESs and 𝑓 C == , the new curvatures that controls the conical " .

5 Acknowledgements P. J. C. and M. R. thank the financial support provided by the Spanish Administration (CTQ2011-23140), the Generalitat de Catalunya (2014SGR199 and Xarxa d'R+D+I en Química Teòrica i Computacional, XRQTC) and the European Union (COST Action CODECS CM1002). A. P., D. L. and B. L. acknowledge the financial support from the French National Research Agency (ANR) through the CoConicS Project (No. ANR-13-BS08-0013-03).


Robert Binder

Matteo Bonfanti

David Lauvergnat

Irene Burghardt

First-principles description of intra-chain exciton migration in an oligo(para-phenylene vinylene) chain. I. Generalized Frenkel-Holstein Hamiltonian

A generalized Frenkel-Holstein Hamiltonian is constructed to describe exciton migration in oligo(paraphenylene vinylene) chains, based upon excited-state electronic structure data for an oligomer comprising 20 monomer units (OPV-20). Time-Dependent Density Functional Theory (TDDFT) calculations using the !B97XD hybrid functional are employed in conjunction with a transition density analysis (TDA) to study the low-lying singlet excitations and demonstrate that these can be characterized to a good approximation as a Frenkel exciton manifold. Based on these findings, we employ the analytic mapping procedure of [Binder et al., J. Chem. Phys. 141, 014101 (2014)] to translate one-dimensional (1D) and two-dimensional (2D) potential energy surface (PES) scans to a fully anharmonic, generalized Frenkel-Holstein Hamiltonian. A 1D PES scan is carried out for intra-ring quinoid distortion modes, while 2D PES scans are performed for the anharmonically coupled inter-monomer torsional and vinylene bridge bond-length alternation modes. The kinetic energy is constructed in curvilinear coordinates, using the TNUM procedure. As a result, a fully molecular-based, generalized Frenkel-Holstein Hamiltonian is obtained which is employed for quantum dynamical exciton dynamics simulations, as shown in a companion paper [Binder, Burghardt, J. Chem. Phys., submitted].

I. INTRODUCTION

Exciton migration in semiconducting polymers is a multi-scale process, [START_REF] Köhler | Electronic Processes in Organic Semiconductors[END_REF][START_REF] Barford | Electronic and Optical Properties of Conjugated Polymers[END_REF][START_REF] Barford | [END_REF] which drives photogenerated excitations towards donor-acceptor interfaces where charge separation -i.e., exciton dissociation -occurs. The fastest dynamical components fall into the femtosecond to picosecond regime, as is also the case for the subsequent charge separation step in various donor-acceptor systems. [4][5][6] On this time scale, we expect that molecular details will play an important role in the dynamics of the electronically excited system, necessitating a full molecular-level analysis. Indeed, ultrafast timeresolved experiments reveal that coherent features arise in the femtosecond to picosecond scale dynamics, 4,5,7,8 requiring a more detailed analysis than the conventional Förster rate theory [START_REF] Köhler | Electronic Processes in Organic Semiconductors[END_REF][START_REF] Barford | [END_REF]9 description.

The theoretical description of exciton migration is challenging, both on account of the excited-state electronic structure level of treatment and on account of the dynamical description which needs to include electronphonon coupling as a prominent feature. 4,6,7,10 Full quantum-dynamical simulations are usually restricted to small system sizes and precomputed potential energy surfaces (PES). 6,[11][12][13] By contrast, Ab Initio Molecular Dynamics (AIMD) type approaches, 14 which are for the most part combined with mean-field Ehrenfest a) Electronic mail: burghardt@chemie.uni-frankfurt.de (Author to whom correspondence should be addressed.) dynamics 15,16 or surface hopping dynamics, 17,18 often do not lead to a qualitatively correct result, notably due to the failure to describe detailed balance (in the case of Ehrenfest simulations [19][20][21] ) and due to major approximations in the description of multiple curve crossing events (in the case of essentially all standard quantumclassical methods 22 ). Finally, model Hamiltonianse.g., the Frenkel-Holstein (FH) Hamiltonian [START_REF] Barford | [END_REF]23 and the Haken-Strobl-Reineker (HSR) model [START_REF] Kenkre | Exciton Dynamics in Molecular Crystals and Aggregates[END_REF] -can be most useful to capture the key traits of the dynamics, but necessitate a realistic parametrization. In the context of time-dependent simulations, these model Hamiltonians have been combined with various quantum-classical propagation schemes [START_REF] Tozer | [END_REF]26 as well as with master equation approaches. 11,[27][28][29] More recently, we introduced a first-principles approach that combines the formal framework of the FH Hamiltonian with the explicit computation of PES of oligomer species. Specifically, an analytic mapping procedure is employed which systematically translates PES information to a generalized FH Hamiltonian. 30 Based on this Hamiltonian, high-dimensional quantum dynamical simulations were carried out employing the multiconfiguration time-dependent Hartree (MCTDH) method 31 and its multi-layer (ML-MCTDH) variant. [32][33][34][35] Using this set-up, we showed for oligothiophene (OT) systems that torsional fluctuations induce coherent exciton motion, and that the active dynamical species should be thought of as exciton-polarons generated by trapping due to high-frequency CC stretch modes. 12,13 This scenario accounts for exciton relaxation accompanied by torsional planarization at low temperatures, 12 and temperature-dependent exciton migration under thermal conditions, induced by torsional fluctuations. 13 These studies highlight that the initially photogenerated species -chromophores or spectroscopic units 36,37 -undergo a characteristic two-time scale dynamics where ultrafast trapping due to high-frequency modes is followed by coherent exciton-polaron dynamics driven by fluctuations. 12 Besides accurate multiconfigurational quantum techniques, we also employed approximate semiclassical and quantum-classical approaches, notably the symmetrical quasi-classical Meyer-Miller (SQC/MM) approach and mean-field Ehrenfest dynamics, 38 showing that these approaches are able to produce good results if the zero-point energy of the high-frequency modes is disregarded.

In the present work, we aim to carry out an analogous analysis for another paradigmatic semiconducting polymer system, poly(para-phenylene vinylene) (PPV) and the corresponding OPV-n oligomer species (where n relates to the number of repeat units). PPV has been extensively employed in the substituted form of poly-[2-methoxy-5-((2-ethylhexyl)oxy)phenylene vinylene] (MEH-PPV) in device applications in organic electronics. [START_REF] Köhler | Electronic Processes in Organic Semiconductors[END_REF]39 From previous electronic structure work based on highlevel methods 12,40 -notably the second-order algebraic diagrammatic construction (ADC( 2)) method 41 -we conclude that OPV-n and OT-n oligomers share basic features of the excited singlet manifold, and both can be well described as Frenkel J aggregates (i.e., head-to-tail aligned molecular monomer units), 42 with possible admixtures of charge transfer states. However, di↵erences between OPV and OT type systems may emerge as far as the e↵ect of vibronic couplings and nonadiabaticity are concerned. Specifically, we seek to investigate whether trapping e↵ects due to high-frequency modes are as pronounced in OPV-n as in their OT-n analogs, and whether the torsion-induced dynamics are similar.

Hence, the first part of the present analysis is concerned with the construction of a generalized FH Hamiltonian for OPV-n systems, whereas the second partpresented in a companion paper [START_REF] Binder | First-principles description of intrachain exciton migration in an oligo(para-phenylene vinylene) chain. II. ML-MCTDH simulations of exciton dynamics at a torsional defect[END_REF] -addresses the dynamical description based upon this Hamiltonian, using ecient multiconfigurational quantum dynamics methodology (specifically, the ML-MCTDH method as mentioned above), similarly to Refs. [12 and 13]. By analogy with our previous analysis for OT-n systems, we focus on the combined e↵ect of high-frequency modes and "soft" torsional modes, which have an important e↵ect on the electronic structure properties, by creating and removing conjugation breaks. As an additional feature, the present analysis also explores the e↵ect of anharmonic correlations between torsional modes and bond-length alternation (BLA) modes. To this end, two-dimensional PES scans are carried out, and it is demonstrated that the mapping procedure of Ref. [30] is capable of handling this more general situation. In order to carry out a large number of single-point calculations, we employ an approximate electronic structure methodology, i.e., Time-Dependent Density Functional Theory (TDDFT) with the range-separated !B97XD hybrid functional. 44,45 For a comparatively large system, OPV-20, this approach yields a very good description of the Frenkel manifold over a wide range of geometries.

The remainder of the manuscript is organized as follows. Sec. II summarizes electronic structure properties and addresses the vibrational modes that are most relevant to the electronic excitation. Sec. III addresses the construction of the generalized Frenkel-Holstein Hamiltonian based upon the analytic mapping procedure of Ref. [30] and presents the explicit potentials obtained for OPV-n, based upon one-dimensional and two-dimensional PES scans. Finally, Sec. IV concludes. An Appendix addresses details of the construction scheme of the FH Hamiltonian.

II. OPV-20 OLIGOMER: ELECTRONIC STRUCTURE CALCULATIONS AND RELEVANT MODES

In this section, we describe the lowest singlet states of the OPV-20 oligomer species under study and ascertain that these states are compatible with the Frenkel exciton model (Sec. II A). To this end, a transition density analysis (TDA) is carried out for representative geometries (Sec. II B). Furthermore, the last part of this section (Sec. II C) is concerned with identifying the vibrational modes which are most strongly coupled to the electronic excitations.

A. Lowest singlet states of OPV-20

The first five singlet excited states of an OPV-20 oligomer -capped with phenylene units at both ends, such that (PV-19)P is a more precise designationwere computed using TDDFT with the range-separated !B97XD hybrid functional 44,45 and the def2-SVP basis set, using the Gaussian 09 program package. [START_REF] Frisch | Gaussian 09, revision d.01[END_REF] Vertical excitation energies, oscillator strengths and dominant orbital contributions of the S n , n = 1, . . . , 5, states calculated with a C 2h symmetry constraint are collected in Table I. Subsequently, various PES cuts were carried out, i.e., 1D cuts along an intra-phenylene quinoid distortion coordinate and 2D cuts along the ring-torsional and BLA modes at the vinylene junction (see Sec. II C for the definition of the relevant coordinates and Sec. III C for details about the PES scans). Overall, excited-state calculations were carried out for more than 500 molecular geometries.

At the DFT level of treatment, a ground state equilibrium value of ✓ eq = 0 for the ring-torsional modes is obtained. This is in contrast to earlier calculations for smaller systems -notably the (PV-7)P oligomer -using second-order Møller-Plesset (MP2) theory 40 which provides evidence for a ⇠30 twisted geometry and a shallow barrier at the planar geometry. Similarly, the OT-FIG. 1. Electron-hole maps ("omega matrices", see Eq. ( 1)) obtained by TDA for OPV-20, using the wavefunction analysis tool of Ref. [47]. Omega matrices are shown for the lowest five adiabatic excited states (S1, . . . , S5), for the planar geometry (✓10,11 = 0 ) and a twisted geometry (✓10,11 = 45 ) that will serve as the initial condition for our quantum dynamics simulations in the companion paper. 8 systems we investigated 12,30 exhibit a slightly twisted ground-state geometry. Due to the shallow barrier in both cases, of the order of 10 meV or less, the vibrational wavefunction in the electronic ground state is delocalized 13,[START_REF] Tukachev | [END_REF] such that the ground-state structures are on average planar, i.e., broadly distributed around the planar geometry.

In line with previous studies of shorter OPV-n species 40,49 notably using the ADC(2) method as mentioned above, 40 the first singlet states shown in Table I can be assigned to ⇡ ⇡ ⇤ transitions whose nodal structure is well described by a Frenkel exciton model. 40,49 Within the Frenkel exciton picture, the oligomer chain corresponds to a J aggregate, 42 composed of head-to-tail aligned monomer units. From Table I, we can infer that the nodeless S 1 ( 1 B u ) state carries the dominant oscillator strength (f = 17.28); further, we note that transitions from the electronic ground state to even-parity (gerade) excited states are optically forbidden, whereas transitions to higher odd-parity (ungerade) states are allowed, albeit with a sharply decreasing transition probability. A very similar agreement with the J aggregate model is found for other oligomer systems, notably oligothiophenes. 12 Given that TDDFT correctly describes the excitonic character of the low-lying singlet states, and the computed PESs exhibit the correct topology (see Sec. III.C), the TDDFT methodology is a viable approach in the present context.

B. Transition density analysis

The electronic character of the low-lying singlet states can be best captured by the TDA 40,[START_REF] Plasser | TheoDORE -A package for Theoretical Density, Orbital Relaxation and Exciton analysis[END_REF][49][50][51] which yields an electron-hole (e-h) representation as illustrated in Fig. 1. Here, a characteristic charge transfer index ⌦ (k) ⌫µ is computed for pairs of molecular fragments (⌫, µ) 51

⌦ (k) ⌫µ = 1 2 X n2⌫ X m2µ (⇢ (k) S) nm (S⇢ (k) ) nm (1)
where ⇢ (k) is the single-particle transition density matrix for the kth excited state (which relates to the creation of an e-h pair, from the perspective of many-body perturbation theory [START_REF]many-body perturbation theory, the simplest treatment relating to the e-h state representation is in terms of the Tamm-Danco↵ approximation. 61,62 Here[END_REF] ), and S is the overlap matrix. Both ⇢ (k) and S are specified in the atomic orbital (AO) basis with indices (n, m), and the summations run over AOs belonging to the ⌫th and µth molecular fragments. If ⌦ (k) ⌫µ was strictly diagonal (⌦

(k) ⌫µ = ⌫µ ⌦ (k)
⌫µ ), the corresponding state would be of pure Frenkel type, with localized e-h pairs. O↵-diagonal contributions ⌫ 6 = µ, indicate charge transfer character. The ⌦ denotes the position of the electron and µ denotes the position of the hole. [START_REF]many-body perturbation theory, the simplest treatment relating to the e-h state representation is in terms of the Tamm-Danco↵ approximation. 61,62 Here[END_REF] Fig. 2 illustrates the molecular fragments chosen in the TDA. A symmetric partitioning is used, where each phenylene unit and the left and right "half" vinylene units are assigned to a "monomer" fragment. Obviously, alternative partitionings are possible, including the standard chemical partitioning where a pair of neighboring phenylene and vinylene units are assigned to a repeat unit, or a partitioning where each phenylene and vinylene moiety is addressed as an individual unit. The present scheme is chosen because (i) it is immediately compatible with the procedure of Ref. [30] where homo-oligomer species are addressed, and (ii) this partitioning corresponds to the natural symmetry of the relevant vibrations, as detailed in Sec. II C. Also, a direct comparison can be made with our closely related analysis of OT-n systems. 12 In Fig. 1, the "omega matrices" ⌦ (k) are shown for the first five singlet states for two geometries, i.e., the planar geometry (✓ 10,11 = 0 ) and a twisted geometry (✓ 10,11 = 45 ), where ✓ 10,11 denotes the ring torsion at the center of the lattice, connecting the monomer units 10 and 11 (see Fig. 2 and Fig. 4). The diagonal dominance of the TDA maps indicates that the Frenkel model employed in the present study, which is restricted to the center-of-mass exciton (⌫ = µ), is a good approximation. While small nearest-neighbor charge transfer contributions exist, these are less pronounced and the sum of charge transfer indices adds up to zero, i.e., no net charge transfer between the fragments occurs; hence, the relevant states correspond to local excitations. 53 Further, a significant di↵erence is observed for the two geometries: For the planar geometry (✓ 10,11 = 0 ), the S n , n = 1, . . . , 5 states exhibit (n 1) nodes along the diagonal (⌫ = µ) and feature the characteristic particlein-the-box type eigenstate structure. 40,49 For the twisted geometry (✓ 10,11 = 45 ), the e-h probabilities decrease to zero at the center of the lattice, indicating a conjugation break. Indeed, the 45 angle involves an orthogonal position of neighboring phenyl rings, as further explained below (see Sec. II.B). This observation is similar to related studies of OT oligomer chain segments. 12 The Supp. Mat. provides further details about the evolution of the e-h map as a function of the ✓ 10,11 angle and the forma-tion of a conjugation break (see Fig. S1).

Overall, we conclude that (i) the Frenkel J aggregate model is a valid description of the electronic structure of the lowest OPV singlet states, supporting the description in terms of a FH Hamiltonian, and (ii) the torsional modes will play an important role in modifying the electronic structure properties and inducing conjugation breaks, very similarly to our previous analysis for OT type systems. 12

C. Relevant vibrational degrees of freedom

We now identify a subset of local vibrations that drive the exciton dynamics and will be included in the FH Hamiltonian (noting that these coordinates do not represent a complete set of relevant modes). From the properties of the HOMO-LUMO transition, depicted in Fig. 3, we can infer which vibrational modes are most strongly coupled to the optically bright S 1 S 0 transition. Upon excitation to the Frenkel excitonic manifold, the formal CC single bonds connecting the vinylene moieties to the neighboring phenylene rings will contract acquiring partial double-bond character, while the formal CC double bonds of the vinylene bridges will elongate leading to a BLA mode. The associated inter-monomer torsional degrees of freedom are expected to sti↵en due to the loss of "pure" single-bond character within the vinylene bridges. Hence, it can be expected that the intermonomer BLA and torsion modes are highly correlated. Furthermore, the partial inversion of bond orders leads to a quinoid distortion of the phenylene rings, see the Lewis structure of panel (b) (bottom) of Fig. 3.

To assign these modes, we refer to the monomerbased partitioning of Fig. 2. As illustrated in Fig. 4, three primary vibrational coordinates are identified: sitecorrelated inter-monomer torsion modes ✓ n,n+1 and vinylene bridge BLA modes y n,n+1 as well as site-local intramonomer quinoid-distortion modes x n . 40,[START_REF] Tozer | [END_REF] Here, we deliberately choose the same mode labels as in our previous studies of OT systems, 12,30 as both sets of modes can be considered as structurally similar -and they also turn out similar as far as their role in the dynamics is concerned. [START_REF] Binder | First-principles description of intrachain exciton migration in an oligo(para-phenylene vinylene) chain. II. ML-MCTDH simulations of exciton dynamics at a torsional defect[END_REF] In order to define these vibrational degrees of freedom unambiguously, we adopt the atom numbering shown in Fig. 2 for an OPV dimer fragment. Here, a symmetric definition of the monomer unit has been chosen, for which each vinylene bridge is partially assigned to two monomer units. Due to the more complex composition of the OPV monomer unit in comparison to OT, the definition of the site-correlated modes is also somewhat more involved.

In order to define the active inter-monomer torsion mode ✓ n,n+1 , we first introduce the two torsional angles ✓ n,L and ✓ n,R which describe the relative angles between the phenylene ring of site n and the neighboring, planar vinylene bridges linking it to sites n 1 and n + 1:

✓ n,L = # ⇣ C (n) 5 -C (n) 4 -C (n) 8 -C (n 1) 7 ⌘ ✓ n,R = # ⇣ C (n) 2 -C (n) 1 -C (n) 7 -C (n+1) 8 ⌘ (2)
where # denotes a dihedral angle. We now construct the positive linear combination of these torsional angles, 40

✓ n,n+1 ⌘ ✓ sum n,n+1 = 1 2 (✓ n,R + ✓ n+1,L ) (3)
Adopting this convention, ✓ n,n+1 = 0 refers to a fully planar structure, while ✓ n,n+1 = 45 refers to a conformation in which the two molecular fragments separated by ✓ n,n+1 are perpendicular to each other. (The negative superposition ✓ di↵ n,n+1 = 1 2 (✓ n,R ✓ n+1,L ) corresponds to a buckling mode [START_REF] Tukachev | [END_REF] which does not modify the relative orientation of the phenylene rings; this type of mode is not included here but could also play a role in the photochemistry.)

In order to define the active inter-monomer BLA mode y n,n+1 , we introduce the two formal single bond lengths s n,L and s n,R as well as the formal double bond length d n,n+1 of the vinylene bridge, 40

s n,L = r ⇣ C (n) 4 -C (n) 8 ⌘ s n,R = r ⇣ C (n) 1 -C (n) 7 ⌘ d n,n+1 = r ⇣ C (n) 7 -C (n+1) 8 ⌘ (4)
where r denotes a bond length. Following Ref. [40], we now introduce the BLA coordinate y n,n+1 as a linear combination of s sum n,n+1 = 1 2 (s n,R + s n+1,L ) and d n,n+1 ,

y n,n+1 = 1 2 (s n,R + s n+1,L ) d n,n+1 = s sum n,n+1 d n,n+1 (5)
Adopting this convention, positive values of the BLA coordinate indicate the normal bond order of the vinylene bridge, while negative values indicate an inversion of the bond order. We further introduce the deviation y n,n+1 from the ground state optimal value y eq y n,n+1 ⌘ y n,n+1 y eq (6) Using this convention, y n,n+1 > 0 refers to an elongation of the formal single bonds and a compression of the formal double bond of the associated vinylene bridge, while y n,n+1 < 0 refers to a compression of the formal single bonds and an elongation of the formal double bond. Finally, the intra-monomer quinoid-distortion mode x n will be defined in a similar fashion to the symmetric ringbreathing mode in OTs. 12 This definition can be inferred directly from the intra-ring bond order inversion depicted in Fig. 3. Along the trajectory of the quinoid-distortion mode, the intra-ring bonds C

(n) 1 -C (n) 2 , C (n) 1 -C (n) 6 , C (n) 4 -C (n) 3 and C (n) 4 -C (n) 5
are elongated in a synchronous fashion, while bonds C

(n) 2 -C (n) 3 and C (n) 5 - C (n) 6
are synchronously contracted in such a way that all bond angles within the six-membered ring remain constant at around 120 . Hence, we define x n as

x n = 1 4 ⇣ r ⇣ C (n) 1 -C (n) 2 ⌘ + r ⇣ C (n) 1 -C (n) 6 ⌘ +r ⇣ C (n) 4 -C (n) 3 ⌘ + r ⇣ C (n) 4 -C (n) 5 ⌘⌘ (7)
We further introduce the deviation of the quinoiddistortion coordinate, x n , from its ground state optimal value x eq

x n ⌘ x n x eq (8)

Using this convention, x n > 0 refers to a quinoid distortion of the phenylene ring with respect to its ground state optimal structure, while x n < 0 refers to an antiquinoid distortion.

III. GENERALIZED FRENKEL-HOLSTEIN HAMILTONIAN AND ANALYTIC MAPPING PROCEDURE

Based on the above analysis, we now construct a generalized FH Hamiltonian as a function of the selected set of coordinates {x n , y n,n+1 , ✓ n,n+1 }. As pointed out in Sec. II B, the Frenkel picture represents to a good approximation the diagonally dominant e-h states visualized in Fig. 1, i.e., the first excited singlet states essentially feature localized e-h states. Alternatively, one can state that the center-of-mass exciton with position n = (⌫ + µ)/2 is a good approximation to the actual excitonic state.

Similarly to our previous investigations of the OT system, 12,30 the Hamiltonian is constructed by employing the analytic mapping procedure of Ref. [30]. This procedure translates PES cuts obtained from electronic structure calculations of oligomer species to an FH Hamiltonian. The procedure entails no loss of information for the adiabatic S 1 and S 2 PES, and yields an approximate representation for the remaining Frenkel manifold, in line with the structure of the Frenkel Hamiltonian.

The mapping procedure can be understood as a special type of diabatization scheme, which employs the solution to an inverse eigenvalue problem. The procedure works in a strictly local, pointwise fashion and returns values of the potential functions that enter into the FH type Hamiltonian specified below. In the appendix, a summary of the procedure is given; details can be found in Ref. [30].

By analogy with Refs. [12 and 30], PES scans are carried out to construct a generalized FH Hamiltonian as a function of the set of modes {x n , y n,n+1 , ✓ n,n+1 }. Differently from Ref. [30], we now include two-dimensional PES scans for the vinylene bridge modes {y n,n+1 , ✓ n,n+1 } which are expected to be anharmonically coupled. By contrast, the quinoid distortion modes {x n } are included in terms of a 1D scan, similarly to our previous approach.

In the following, we first present the generalized FH Hamiltonian employed in this work (Sec. III A), followed by the details of the monomer potentials appearing in this Hamiltonian (Sec. III B), and the PES scans that were used as input for the monomer-based functions of the FH Hamiltonian (Sec. III C).

A. Generalized FH Hamiltonian

The general form of the exciton Hamiltonian in the basis of Frenkel configurations |ni, i.e., single-monomer excitations in the direct-product space spanned by N monomers is as follows,

ĤFH = N X n,n 0 =1 ĤFH n,n 0 |nihn 0 | + Ĥbath 1 (9)
combining contributions for an N -site system ( Ĥn,n 0 ) with an external harmonic-oscillator bath ( Ĥbath ). Here, Ĥn,n 0 comprises diagonal kinetic ( T ) and on-site potential ( V site n ) contributions, as well as next-neighbor excitonic couplings ( V exc n,n 0 ),

ĤFH n,n 0 = n,n 0 T + n,n 0 V site n + V exc n,n 0 (10)
The kinetic energy operator is taken to be electronically diagonal, in keeping with a diabatic representation, and takes the following form,

T = 1 2 N X n=1 G xx p2 xn + N 1 X n=1 ⇣ G yy p2 yn,n+1 + G ✓✓ p2 ✓n,n+1 ⌘ ! (11)
This expression refers to curvilinear coordinates, using the TNUM procedure 55 applied to an OPV dimer. As discussed in the Supp. Mat., the metric tensor is assumed to be constant and o↵-diagonal momentum couplings are found to be negligible to a good approximation. The diagonal entries of the metric tensor G are equivalent to inverse masses, i.e., G xx = m 1

x and G yy = m 1 y , and inverse moments of inertia, i.e., G ✓✓ = I 1 ✓ , respectively. Explicit values of these quantities are collected in Table II.

where c E is the constant excitation energy of a single monomer and vG and vE are ground-state and excitedstate monomer potentials, respectively. Further, the excitonic couplings of Eq. ( 10) are given as nearest-neighbor couplings,

V exc n,n 0 ⇣ ŷn 1,n , ŷn,n+1 , ✓n 1,n , ✓n,n+1 ⌘ = ŵ ⇣ ŷn 1,n , ✓n 1,n ⌘ + ŵ ⇣ ŷn,n+1 , ✓n,n+1 ⌘ (15)
where an additional dependence on the quinoiddistortion modes xl (l = n 1, n, n + 1) has been neglected.

In the above expressions, the summations over the full physical range of N monomer units include contributions by "virtual" or "dangling" modes at the left and right ends of the chain, i.e., torsion modes ✓0,1 and ✓N,N+1 as well as BLA modes ŷ0,1 and ŷN,N+1 , respectively. These modes are included in order to counteract edge e↵ects; they are fixed at their ground-state equilibrium values and contribute constant values to the potential energy, in the absence of any kinetic energy contributions.

In the present study, we do not include couplings between the intra-monomer modes (i.e., {x n }) and "left" (n 1, 1) and "right" (n, n + 1) inter-monomer modes (i.e., {y n 1,n , ✓ n 1,n } etc.). However, we do include anharmonic couplings between the inter-monomer BLA (y n,n+1 ) and ring torsion (✓ n,n+1 ) modes. As a result, the following form of the monomer potentials is employed, vG/E ⇣ xn , ŷn 1,n , ŷn,n+1 , ✓n 1,n , ✓n,n+1

⌘ = vG/E ⇣ xn ⌘ +v G/E ⇣ ŷn 1,n , ✓n 1,n ⌘ +v G/E ⇣ ŷn,n+1 , ✓n,n+1 ⌘ (16)
In the next sections (Sec. III B and Sec. III C), we show how these potentials are constructed explicitly. Finally, the bath Hamiltonian Ĥbath of Eq. ( 9) comprises sets of harmonic oscillator modes that are bilinearly coupled to the low-frequency torsional mode(s). The bath is constructed such as to act exclusively on the torsions since these are e↵ectively damped (or overdamped) in the real polymer material. Conversely, the high-frequency modes are known to undergo underdamped vibrations on time scales extending up to a picosecond. 5,8 Hence, Ĥbath comprises sets of bath modes { bj,k } acting on the individual torsions ✓k,k+1 ,

Ĥbath ({ b}) = N 1 X k=1 N b X j=1 p2 b j,k 2m b + m b ! 2 j 2 bj,k c j m b ! 2 j ⇣ ✓k,k+1 ✓ 0k ⌘ ! 2 (17)
where N b denotes the number of modes of each sub-bath, m b is the common mass of the bath oscillators, and ! j are their frequencies. Further, ✓ 0k = h✓ (0) k,k+1 i denotes the expectation value of the relevant torsional mode at time t = 0. In order to achieve phenomenological damping of the active torsion(s), the couplings c j are chosen such as to conform to an Ohmic spectral density [START_REF] Binder | First-principles description of intrachain exciton migration in an oligo(para-phenylene vinylene) chain. II. ML-MCTDH simulations of exciton dynamics at a torsional defect[END_REF] J(!)

= I ✓ !, i.e., c j = ! j p (2/⇡) I ✓ m b !.
The number of bath oscillators needs to be large enough to guarantee e↵ectively irreversible decay behavior for a given frequency spacing !, which is in turn chosen to match the Poincaré recurrence time T P = 2⇡/ ! with the observation interval. In the present case, N b = 20 represents a minimal set of bath modes that fulfill these requirements. The details of the implementation of the discretized bath are described in the companion paper. [START_REF] Binder | First-principles description of intrachain exciton migration in an oligo(para-phenylene vinylene) chain. II. ML-MCTDH simulations of exciton dynamics at a torsional defect[END_REF]

B. Monomer potentials: Functional form

The one-and two-dimensional functions v G , v E and w of Eqs. ( 13), ( 14) and ( 15) as well as the excitation energy of a single monomer of Eq. ( 14), c E = 0.199 E h = 5.41 eV, are determined via the abovementioned mapping procedure. 30 As the mapping procedure works in a discrete, point-wise fashion, the resulting monomer potentials have to be fitted to a continuous analytical form. Specifically, the quinoid-distortion monomer potentials were fitted to a Morse oscillator functional form,

vG/E ( x) = D G/E exp ↵ G/E x x G/E,0 1 2 +E G/E,0 (18)
where x ⌘ xn refers to the deviation of xn from its ground-state optimal value. D G/E denotes the well width and ↵ G/E the well depth in the monomer ground and excited state, respectively.

x G/E,0 describes the equilibrium coordinate value (as gauged to the groundstate optimal value) and E G/E,0 is an energetic o↵set. See Tab. III for the parameters of the quinoid-distortion monomer potentials.

The functional form of the two-dimensional, nonseparable BLA/torsion monomer potentials is somewhat more involved, combining polynomial contributions of the BLA mode and cosine contributions of the torsion mode,

vG/E ( ŷ, ✓) = j max G/E X j=0 k max G/E X k=0 c (j,k) G/E ŷj cos ⇣ k ✓⌘ (19)
Likewise, the potential type coupling takes the analogous form ŵ( ŷ, ✓) =

j max w X j=0 k max w X k=0 c (j,k) w ŷj cos ⇣ k ✓⌘ ( 20
)
where ŷ = ŷn refers to the deviation of the BLA mode ŷn,n+1 from its ground state optimal value and where (j max 19)- (20), with the specified values of j max and k max , are of the order of 1 meV (more specifically, = 0.469 meV for the groundstate monomer potential v G , = 1.037 meV for the excited state monomer potential v E , and = 0.865 meV for the coupling w).

The higher orders in the description of the excitedstate monomer potential and potential type coupling reflect the more complex structure of the corresponding discrete data sets obtained from our analytic mapping procedure, see also Sec. III C. See Tabs. IV, V and VI for the BLA/torsion ground and excited-state monomer potential as well as potential type coupling parameters, respectively. As further explained in the Appendix, the mapping procedure returns a set of {v G , v E , w} values from adiabatic {S 0 , S 1 , S 2 } data at a given geometry (noting that the full Frenkel manifold is generated by the mapping procedure, see App. A 3 for details). Hence, the input data for the mapping procedure consist in a set of electronic structure data points. In the present application, we choose data from (i) one-dimensional PES scans along the intra-monomer x n mode, and (ii) two-dimensional PES scans as a function of the junction modes ✓ 10,11 and y 10,11 . As explained in App. A 3 and Ref. [30], a full parametrization of the FH Hamitonian is obtained from these selected PES scans. In the following subsections, these two types of scans are discussed and general features of the PES are highlighted.

One-dimensional quinoid-distortion PES scan

For reference, one-dimensional rigid PES scans were performed for all three types of coordinates, and the corresponding monomer potentials are illustrated in Fig. 5. While a more detailed 2D scan is subsequently undertaken for the junction modes (✓ 10,11 and y 10,11 ), the onedimensional PES scan for the quinoid distortion mode x n is employed in the final definition of the PES. Hence, we now summarize the conditions for the latter scan (panels (c) and (f) in Fig. 5).

For the quinoid distortion mode, a one-dimensional PES scan was carried out by varying the values of the central modes x 10 and x 11 in a synchronous fashion while keeping all other modes at their ground state equilibrium geometries. This procedure conserves the overall symmetry. Here, x 10 and x 11 were varied from 0.16 Å to 0.16 Å in increments of 0.01 Å for a total of 33 steps. In panel (c) of Fig. 5, the molecular adiabatic ground and excited states are visualized, while panel (f ) shows the monomer ground and excited state potentials and potential type coupling as obtained from our analytical mapping procedure. 19) for the functional form of the monomer potential.

c (0,0) G c (0,1) G c (0,2) G c (0,3) G c (0,4) G 8.25 ⇥ 10 3 8.56 ⇥ 10 3 1.67 ⇥ 10 3 2.51 ⇥ 10 3 1.16 ⇥ 10 3 c (1,0) G c (1,1) G c (1,2) G c (1,3) G c (1,4) G 9.02 ⇥ 10 3 8.86 ⇥ 10 3 1.23 ⇥ 10 3 2.63 ⇥ 10 3 1.22 ⇥ 10 3 c (2,0) G c (2,1) G c (2,2) G c (2,3) G c (2,

Two-dimensional torsion/bond-length alternation PES scan

Next, a two-dimensional PES scan along both torsion and bond-length alternation modes was carried out in order to capture correlation e↵ects between these modes. To this end, the values of the central torsion mode ✓ 10,11 and the central BLA mode y 10,11 were varied simultaneously while keeping all other modes at their ground state equilibrium geometries, thereby following the "defect scan" procedure introduced in Ref. [30]. Here, ✓ 10,11 was varied from its ground state optimal value of 0 to 90 in increments of 5 for a total of 19 steps, while y 10,11 was varied from 0.14 Å to 0.12 Å in increments of 0.01 Å for a total of 27 steps. In the latter case, the s sum 10,11 coordinate of Eq. ( 5) was incremented by 0.005 Å, while the d 10,11 coordinate of Eq. ( 4) was decremented by the same value for each step. As a result, a total number of 513 individual molecular geometries was evaluated.

One-dimensional sections of the obtained surfaces are shown in Fig. 5. In panels (a) and (b), the molecular adiabatic ground and excited states are visualized as a function of ✓ 10,11 for y 10,11 = 0 and as a function of y 10,11 for ✓ 10,11 = 0 , respectively. In panels (d) and (e), the corresponding monomer ground and excited state potentials v G and v E as well as the potential type coupling w as obtained from the analytic mapping procedure are shown. Fig. 6 shows the two-dimensional surfaces. Panels (a) through (c) show the molecular S 0 , S 1 and S 2 surfaces, respectively, while panels (d) through (f ) show the monomer ground state potential v G , monomer excited state potential v E and potential type coupling w, respectively. The correlation e↵ects between ✓ n,n+1 and y n,n+1 are subtle for the S 0 and S 2 states as well as for v G (which is directly derived from S 0 ). In contrast, for the S 1 PES, planarization of the torsion entails a quite significant shift of the BLA mode to negative values (from 0.018 a.u. at 45 to 0.057 a.u. at 0 ), which corresponds to a (partial) inversion of the bond-orders of the vinylene bridge, and vice versa. This functional behavior is amplified for v E , as expected. For w, which exhibits the strongest correlation e↵ects, the absolute value of the next-neighbor excitonic coupling is found to be maximal for planar geometries and negative values of the BLA mode (indicating maximum ⇡-overlap and a BLA mode adjusted to the planar torsion).

D. PES features

The adiabatic and monomer potentials discussed above exhibit some generic features that are very similar to the potentials we obtained previously for the related OT species. 12 First, the S 1 and S 2 potentials exhibit a strong dependence on the torsional coordinate; these two states are nearly degenerate at the twisted geometry, but S 1 is more stable by ⇠0.1 eV at the planar geometry. This is in agreement with the TDA of Fig. 1. The S 1 state exhibits a larger curvature than the S 0 ground state (i.e., the respective torsional periods are 1431 fs in S 0 , 1123 fs in S 1 , and 1424 fs in S 2 ). Excitation at the twisted geometry will lead to planarization in the excited state, and likely involves nonadiabatic e↵ects due to the proximity of S 1 and S 2 . Besides these similarities, a more detailed comparison with the OT case reveals that the di↵erence of curvatures in the S 0 and S 1 states is significantly more pronounced in the latter.

Turning to the vinylene BLA modes y n,n+1 and the intra-ring quinoid distortion modes x n , these reflect a quinoid distortion as a result of photoexcitation, as discussed in Sec. II C. As a result, the potentials of Fig. 5bc) show that a modified equilibrium geometry prevails in the excited states. Upon photoexcitation, excited-state stabilization therefore occurs, leading to a trapping effect, or exciton-polaron formation. Again, these e↵ects are similar to OT systems, and the relevant frequencies are similar as well (with e↵ective ground state frequencies of around 500 cm 1 for the BLA modes and around 1200 cm 1 for the quinoid distortion modes). However, both molecular and monomer potentials exhibit a much less anharmonic shape when compared to the OT intermonomer bond-stretch modes. Also, both the change in equilibrium value and the associated energetic stabilization when going from the molecular ground state to the excited states is much less pronounced. The same is true for the change in curvature. Hence, we expect that the contribution of both the BLA modes and quinoid distortion modes to the polaronic stabilization will be considerably weaker as compared with the analogous OT modes.

Regarding the monomer potentials and diabatic couplings, these reflect the scaling properties of oligomer systems with respect to system size and feature (i) larger energetic di↵erences than the oligomer potentials, and (ii) large excitonic couplings, of the order of w = 1.2 eV, very similarly to the OT systems we investigated. 12 Here, the negative sign of the excitonic coupling is consistent with the J aggregate property, as mentioned above. 30,42 We further note that the excitonic coupling does not feature a strong dependence on the torsional mode, differently from other model Hamiltonians where the electronic coupling tends to vanish at a conjugation break. 56 That is, the e↵ective decoupling of subunits at a conjugation break is realized by a di↵erent type of "diabatization" in our procedure 57 (see also Fig. S2 of the Supp. Mat. which illustrates the angle dependent eigenvector structure generated by our mapping scheme).

IV. CONCLUSION

To summarize, a generalized Frenkel-Holstein Hamiltonian was constructed for OPV oligomers, as a function of high-frequency quinoid distortion and bond-length alternation modes, as well as torsional low-frequency modes. The mapping procedure developed in Ref. [30] was employed to translate PES sections in one and two dimensions to the Frenkel-Holstein form. The procedure is a highly useful tool in the present context, since it produces an exact one-to-one mapping between computed electronic structure data and eigenvalues of the Frenkel-Holstein Hamiltonian for the S 1 and S 2 PES which play a key role in the dynamical evolution of these J aggregate type systems. The mapping procedure can be carried out at di↵erent levels of molecular detail, e.g, relying on one-dimensional PES scans or including anharmonic couplings between pairs of modes -here, BLA and torsional modes at the vinylene junctions. In the present work, we demonstrated for the first time that the procedure can straightforwardly handle two-dimensional PES sections. Overall, the mapping procedure rigorously connects the electronic structure perspective with a class of model Hamiltonians which are employed extensively to describe transport phenomena in semiconducting polymers. While other methods, notably machine learning approaches, [START_REF] Smith | [END_REF]59 are more generally applicable, the present method is a highly accurate way to combine FH type Hamiltonians with anharmonic potentials, including large amplitude motions of soft modes. In future work, the optimal selection of active modes should be further explored, possibly guided by more approximate excitedstate dynamics.

Besides the OT systems that we studied previously, 12,30 OPVs represent basic paradigm materials in organic photovoltaics. The molecular-level information obtained in the present study informs about the role of trapping e↵ects due to high-frequency modes, and about the role of torsional displacements. We commented in several instances on the comparison with OTs, for which we constructed a very similar extended FH Hamiltonian. Besides obvious similarities regarding the PES shape as a function of high-frequency vs. lowfrequency modes and the trend towards planarization in the excited state, some di↵erences are observed as well. Notably, it appears that the high-frequency modes in OPVs give rise to less pronounced exciton-phonon coupling, and the di↵erence in curvature of the groundstate vs. excited-state torsional potentials is also less pronounced.

In a companion paper [43], we describe quantum dynamical results obtained with the present extended Frenkel-Holstein Hamiltonian, to permit a direct comparison with our previous results for OT systems. The combination of ab initio parametrized FH Hamiltonians and quantum dynamical calculations in many dimensions provides important benchmarks for more approximate methods -both regarding the PES representation and regarding the dynamical methodologies -that are used standard-wise in the context of molecular exciton dynamics.

SUPPLEMENTARY MATERIAL

See supplementary material for (i) the explicit form of the kinetic energy derived from the TNUM procedure and additional information on the TDA and excitonic eigenvectors derived from the FH Hamiltonian, and (ii) a For-tran90 program providing the potential energy surfaces obtained as described in the text. (See the file README for instructions on how to compile and run the test code.)

  FIG. 2. Definition of monomer units and atom numbering for the oligo(para-phenylene vinylene) system under study. The boundary of the two monomers is indicated by a dashed line bisecting the vinylene double bond.

FIG. 3 .

 3 FIG. 3. (a) Frontier orbitals of OPV-20 pertaining to the S1 S0 transition, as obtained at the !B97XD/def2-SVP level; a zoom-in on the central two monomer units is shown. top: initial orbital, highest occupied molecular orbital (HOMO), bottom: final orbital, lowest unoccupied molecular orbital (LUMO). (b) Lewis structures approximating the bonding situation in: top: electronic ground state, bottom: S1 state. The excitation to the S1 entails a (partial) inversion of bond orders.

FIG. 4 .

 4 FIG. 4. Relevant vibrational coordinates of OPV systems: (a) inter-monomer torsion mode ✓n,n+1, (b) inter-monomer BLA mode yn,n+1, (c) intra-monomer quinoid-distortion mode xn; see Eqs. (2)-(8) for the relevant definitions.

G = 2 P k 2 k

 22 . The quality of the fit was assessed by computing the standard deviation = p /(N data N param ) where k are the residues, N data is the number of data points, and N param is the number of parameters of the fit function. The standard deviations for the power series of Eqs. (

4 ) G 7 . 5 c 4 c 3 c 11 ⇥ 3 c 48 ⇥

 4754311348 79 ⇥ 10 2 5.08 ⇥ 10 4 3.73 ⇥ 10 3 1.19 ⇥ 10 4 1.52 ⇥ 10 4TABLE V. Excited-state bond-length alternation (yn,n+1) and torsion (✓n,n+1) monomer potential parameters. All values given in atomic units. See Eq. (19) for the functional form of the monomer potential. All c ⇥ 10 2 1.47 ⇥ 10 4 1.54 ⇥ 10 2 1.12 ⇥ 10 4 3.10 ⇥ 10 ⇥ 10 2 1.63 ⇥ 10 3 2.45 ⇥ 10 2 4.06 ⇥ 10 4 2.74 ⇥ 10 ⇥ 10 2 9.44 ⇥ 10 3 1.48 ⇥ 10 2 8.54 ⇥ 10 4 1.87 ⇥ 10 10 3 9.60 ⇥ 10 3 6.59 ⇥ 10 3 8.26 ⇥ 10 3 3.20 ⇥ 10 10 2 6.46 ⇥ 10 2 6.30 ⇥ 10 3 1.74 ⇥ 10 2 5.23 ⇥ 10 3

FIG. 5 .

 5 FIG.5. Panels (a) through (c): molecular adiabatic ground (S0) and excited state (S1 and S2) PESs, (d) through (f ): monomer ground and excited state potentials vG and vE and potential type coupling w as obtained from analytic mapping procedure of Ref.[30]. Top to bottom: torsion modes ✓n,n+1 (one-dimensional section for yn,n+1 = 0), bond-length alternation modes yn,n+1 (one-dimensional section for ✓n,n+1 = 0 ) and quinoid-distortion modes xn. Note the di↵erent y-axis scalings and partially broken y-axes. All scans were carried out for the central junction (n= 10).

FIG. 6 . 3 . 21 ⇥ 3 c 1 . 18 ⇥ 3 c 2 c 3 c

 63213118323 FIG. 6. Panels (a) through (c): two-dimensional molecular adiabatic ground and excited state potential energy surfaces S0, S1 and S2 as dependent on torsion mode ✓n,n+1 and bond-length alternation mode yn,n+1. Panels (d) through (f ) monomer ground and excited state potentials vG and vE as well as potential type coupling w. Note the di↵erent y-axis scalings. All scans were again carried out for the central junction (n= 10).

TABLE I .

 I Excitation energies E, symmetries, oscillator strengths f and dominant orbital contributions (from initial orbital to final orbital) of the first five singlet excitations Sn, n = 1, . . . , 5, of OPV-20 as optimized in the electronic ground state with a C 2h symmetry constraint, employing the !B97XD functional and def2-SVP basis set.

	Sn E [eV] symmetry f initial orbital final orbital
	1 2.965	Bu	17.28 HOMO	LUMO
	2 3.048	Ag	0.00 HOMO	LUMO + 1
			HOMO 1 LUMO
	3 3.157	Bu	1.56 HOMO	LUMO + 2
			HOMO 2 LUMO
	4 3.284	Ag	0.00 HOMO	LUMO + 3
			HOMO 3 LUMO
	5 3.419	Bu	0.51 HOMO	LUMO + 4
			HOMO 4 LUMO

TABLE II .

 II Values of the elements of the metric tensor G Eq. (11) as well as corresponding masses m and moments of inertia I, respectively. All values given in atomic units.

							Gxx	Gyy	G ✓✓
					1.05 ⇥ 10 5 7.15 ⇥ 10 5 7.32 ⇥ 10 7
							mx	my	I ✓
					9.48 ⇥ 10 4	1.40 ⇥ 10 4	1.37 ⇥ 10 6
	The on-site potential of Eq. (10) reads as follows,
	V site n	⇣n	x, ŷ,	✓o⌘	= V0 + ˆ n ⇣n ⇣ xn , ŷn 1,n , ŷn,n+1 , ✓n 1,n , ✓n,n+1 x, ŷ, ✓o⌘	⌘
							(12)
	with the ground-state potential
	V0	⇣n	x, ŷ,	✓o⌘	=	⌘
							(13)
	and the di↵erence potential
	ˆ

N X l=1 vG ⇣ xl , ŷl 1,l , ŷl,l+1 , ✓l 1,l , ✓l,l+1 n ⇣ xn , ŷn 1,n , ŷn,n+1 , ✓n 1,n , ✓n,n+1 ⌘ = c E + vE ⇣ xn , ŷn 1,n , ŷn,n+1 , ✓n 1,n , ✓n,n+1 ⌘ vG ⇣ xn , ŷn 1,n , ŷn,n+1 , ✓n 1,n , ✓n,n+1 ⌘

TABLE III .

 III Quinoid-distortion mode (xn) monomer potential parameters. All values given in atomic units. See Eq. (18) for the Morse potential form of the monomer potentials.

	DG	↵G	xG,0	EG,0
	1.43 ⇥ 10 1	3.24 ⇥ 10 1 0.00	0.00
	DE	↵E	xE,0	EE,0
	6.74 ⇥ 10 1	1.51 ⇥ 10 1 9.94 ⇥ 10 2 1.63 ⇥ 10 2
	C. Electronic structure input: PES scans

TABLE IV .

 IV Ground state bond-length alternation (yn,n+1) and torsion (✓n,n+1) monomer potential parameters. All values given in atomic units. See Eq. (

ACKNOWLEDGMENTS

We thank Dominik Brey for his valuable help with the preparation of the figures. Funding by the Deutsche Forschungsgemeinschaft (project BU-1032-2) is gratefully acknowledged.

DATA AVAILABILITY STATEMENT

The data that supports the findings of this study are available within the article and its supplementary material.

Appendix A: Monomer-based Frenkel Hamiltonian

In this appendix, we describe the construction of the generalized Frenkel-Holstein Hamiltonian of Eqs. ( 9)- (17) that was presented in Sec. III A.

Site-based Hamiltonian

Our starting point is a general site-based Hamiltonian,

where { ˆ l } denotes a set of local modes comprising all relevant intra-monomer vibrations on monomer l, while { ⇣l 1,l } and { ⇣l,l+1 } denote sets of correlated modes, i.e., inter-monomer vibrations which are shared by neighboring monomer units. In the OPV-n system, the former modes correspond to the quinoid distortion modes {x l } whereas the latter comprise the BLA modes {ŷ l,l+1 } and torsional modes { ✓l,l+1 }. On the r.h.s. of Eq. (A1), T denotes the kinetic energy, see Eq. ( 11) of the main text. Further, Vsite l denotes the on-site potential of site l, which depends on both local modes ˆ l and correlated modes ⇣l 1,l and ⇣l,l+1 involving site l. (Besides, constant electronic o↵sets are also absorbed into Vsite l .) The excitonic coupling Hamiltonian Vexc l comprises nearest-neighbor couplings at the (l 1, l) and (l, l + 1) junctions and depends on the local modes ˆ l 1 , ˆ l and ˆ l+1 as well as the correlated modes ⇣l 1,l and ⇣l,l+1 .

For simplicity, we will refer to the full set of modes at a given site l as ⌧l ⌘

⇣

ˆ l , ⇣l 1,l , ⇣l,l+1

Assuming a two-level model for the monomer species, the explicit form of the on-site contribution is as follows,

where |G (l) i and |E (l) i denote the monomer ground and excited state of site l, respectively, and ✏ G and ✏ E denote the ground and excited state energies. Besides the latter, the site potential includes ground and excited state monomer potentials vG and vE . These are not restricted to a particular form and can be both anharmonic and non-separable. The excitonic coupling is given as

annihilates (creates) excitations on the lth monomer and simultaneously creates (annihilates) excitations on the neighboring (l ± 1)th monomer. The sign of the coupling function w(⌧ ) determines whether the system is of J-aggregate or H-aggregate type: 30,42 i.e., if w < 0, the system is a J-aggregate whose lowest eigenstate is nodefree and carries maximal oscillator strength. Conversely, w(⌧ ) > 0 corresponds to the H-aggregate case where the highest eigenstate is node-free and carries maximal oscillator strength. Again, w(⌧ ) is not restricted to a particular form and can be anharmonic and non-separable.

E↵ective Frenkel-Holstein Hamiltonian

The site-based Hamiltonian described above is now cast in a basis of Frenkel configurations {|ni}, with n = 1, . . . , N, which represent single-monomer excitations in the direct product space spanned by N monomers, |ni = Q j6 =n |G (j) i ⌦ |E (n) i. These configurations are orthogonal to each other and to the overall ground state, |0i = Q N j=1 |G (j) i. Linear combinations of these configurations yield delocalized Frenkel states,

In the Frenkel basis, the elements of the Hamiltonian of Eqs. (A1)-(A4) take the following form,

where the potential matrix elements comprise on-site and excitonic coupling contributions,

resulting in a real symmetric tridiagonal potential matrix of "Hückel" type. Here, the site potential is given as

noting that Roman letters ( V ) denote operators in the nuclear subspace, while calligraphic letters ( V) denote operators in the nuclear and electronic subspaces. In Eq. (A7), V0 denotes the ground state potential,

with the bare electronic ground state energy c G = N✏ G , and ˆ n denotes the di↵erence potential,

where ✏ = ✏ E ✏ G corresponds to the excitation energy of a single monomer. Due to the excitonic basis, each matrix element in Eq. (A7) depends on all degrees of freedom of the system {⌧ }.

For the excitonic Frenkel coupling elements one obtains,

Eq. (A7) for the diagonal site potential and Eq. (A10) for the excitonic coupling correspond to Eq. ( 12) and Eq. ( 15) of the main text.

Mapping procedure

The mapping procedure developed in Ref. [30] (and extended in Ref. [60]) is an analytic protocol that maps electronic structure based PES information onto the tridiagonal Frenkel-Holstein potential matrix of Eq. (A6). The procedure works in a pointwise fashion and leads to the construction of a coordinate-dependent nearestneighbor coupled Hamiltonian including anharmonicities and mode-mode couplings. The protocol as described in Ref. [30] is restricted to homo-oligomers composed of identical building blocks. For these systems, the parameters of the nearest-neighbor coupled Frenkel-Holstein Hamiltonian as a function of the vibrational coordinates can be found by solving an inverse eigenvalue problem for a Hückel type problem with a single "defect", e.g., the displacement of the torsional coordinate at a specific junction (n, n + 1). The procedure can equally account for the simultaneous displacement of several coordinates at the same junction, e.g., a combined displacement of a torsional and bond-length alternation mode as discussed in the present paper. Typically, a scan in one or several dimensions along these coordinate(s) is carried out.

In the following, we assume tht the active coordinates at the "defect" junction (n S , n S + 1), which are subject to a PES scan, take the values ⌧ S (where "S" stands for scan). The vibrational coordinates at all other sites are taken to be fixed at a reference geometry denoted ⌧ 0 , e.g., the equilibrium geometry. This leads to a matrix representation of the Frenkel-Holstein Hamiltonian which features two types of diagonal entries: (i) for all diagonal entries away from the defect junction, the potential matrix elements read

and (ii) for the diagonal entries at the defect junction n = n S or n = n S + 1, we obtain

Further, the following o↵-diagonal entries are obtained in the Frenkel-Holstein part,

for all entries away from the defect, and

for the couplings at the location of the defect.

Finally, the ground state potential reads

Eqs. (A11)-(A14) represent a tridiagonal Hückel type matrix with a defect at the junction (n S , n S + 1). Despite the apparent simplicity of this "extended Hückel Hamiltonian", the matrix representation of the Hamiltonian cannot be diagonalized analytically. However, the parameters {v G , v E , w} can be analytically determined from a combination of two eigenstates of the Frenkel manifold -i.e., two adiabatic S n states along with the groundstate PES values. For a J-aggregate type system where the lowest singlet state is the bright state, the S 1 and S 2 states are a natural choice.

The input required by the mapping procedure as employed in the present context are therefore, for each geometry, PES data for the adiabatic S 1 and S 2 states. along with the ground-state S 0 PES at the same geometry. These three states, S n , n = 0, 1, 2, are reproduced exactly by the mapping procedure. The remaining states, S n , n = 3, 4, . . . are predicted by the procedure, within the assumptions of a Frenkel model, and were found in good agreement with calculated PESs. 30 Once the parametrization for all geometries has been obtained for the set of available PES points, these data are fitted and the resulting functions can be used in the fully flexible Hamiltonian where all coordinates are variable.

Details of the procedure can be found in Ref. [30] which focuses on the case of correlated modes. The notation employed in the present appendix material is closer to Ref. [60] where the procedure was generalized to a twodimensional lattice. (However, the approach of Ref. [60] employs an alternative scan procedure where all sitebased coordinates are varied in a uniform scan.)


Mathilde Duny

Guillaume Jubeaux

Olivier Geffard

Aurélie Cortéjade

Laure Wiest

Mickaël Nicolas

Emmanuelle Vulliet

Keywords: 0, 12 PFOS: 10, 24 PFOA: 1, 58 PFNA: 0, 2 PFDA: 1, 73 PFUnA: 1, 12 PFDoA: 3, 54 PFTriA: 1, 09 (in ng/g) PFHxS: 0, 34 PFHpS: 0, 2 PFOS: 15, 66 PFOA: 2, 69 PFNA: 0, 6 PFDA: 1, 96 PFUnA: 0, 36 PFDoA: 1, 03 (in ng/g) Carbamazepine: 0, 36 PFHxS: 0, 47 PFHpS: 0, 15 PFOS: 7, 71 PFOA: 2, 44 PFNA: 0, 39 PFDA: 1, 77 PFUnA: 0, 33 PFDoA: 0, 98 PFTriA: 0, 21 ∑= 20

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Development of a quantitative method based on liquid chromatography coupled to tandem mass spectrometry for the quantification of micropollutants in Gammarus fossarum

Development of a quantitative method based on liquid chromatography coupled to tandem mass spectrometry for the quantification of micropollutants in Gammarus fossarum c Monitoring micropollutants is essential to protect water and the aquatic ecosystem. One of the means of assessing river water quality is the use of sentinel species such as Gammarus fossarum. This species is an amphipod crustacean found in most rivers in Western Europe. Its capacity to bioaccumulate environmental contaminants enables to get precious information on water quality. Thanks to its resistance to pollution and its wide availability in the environment, it is a good bioindicator. Actually, many regulated contaminants are monitored using liquid chromatography (LC) coupled to tandem mass spectrometry (MS/MS). However, regarding emerging compounds, data are still incomplete or insufficient for risk assessment. Furthermore, most current extraction methods are based on high sample masses, usually greater than 1 g, and are not suitable for small organisms such as gammarids.

The aim of this study is to develop an analytical method that uses low sample mass, to identify and quantify 45 emerging environmental micropollutants in G. fossarum. These contaminants are mostly Pharmaceuticals and Personal Care Products (PPCPs) (paracetamol, ketoprofen, amiodarone…), perfluorinated compounds (PFOS, PFOA…), flame retardant (hexabromocyclododecan) and organotin compounds (monobutyltin, monophenyltin, tributyltin…). In order to obtain sufficient sensitivity even when using a low matrix mass, great care will be taken to optimize the detection parameters of the mass spectrometer. A single 15 min method was developed to quantify 45 selected micropollutants by LC-MS/MS in both positive and negative ionization mode. The interface parameters (interface temperature, gas flow, etc) were optimized and the best compromise for all the compounds was chosen. For the extraction, a miniaturized QuEChERS (EN Method 15662) with only 0.125 g of matrix was optimized.

The method validation was undertaken according to ICH guide criteria. Limits of quantification were comprised between 0.04 and 313.5 ng/g which is coherent with the analysis of traces in sentinel species.

This method was finally successfully applied to real samples of Gammarus fossarum exposed in different sites in France, chosen according to the levels of contaminants reported in water samplings. This application revealed a global contamination of gammarids with perfluorinated compounds.

  , Carbamazepine, Ketoprofen, Fenofibric acid, Irbesartan, Mercaptodimethur, Boscalide, Irganox 1076, Diclofenac, Triclocarban, Quinoxyfen Amiodarone, Octocrylene, Bisphenol-A, Triclosan, Tetrabromobisphenol A, 4-tert-octylphenol Perfluorinated compounds PFBA, PFPeA, PFPeS, PFHxA, PFBS, PFHxS, PFHpA, PFHpS, PFOA, PFOS, PFNA, PFNS, PFDA, PFDS, PFUnA, PFDoA, PFDoS, PFTriA Flame retardant Hexabromocyclododecane (α, β, γ) Organotin compounds Monobutyltin, Dibutyltin, Tributyltin, Tetrabutyltin, Monophenyltin, Diphenyltin, Triphenyltin Optimization of interface parameters: -Nebulizing gas flow (1-3L/min) -Heating gas flow (3-25L/min) -Drying gas flow (3-20L/min) -Interface temperature (100-400°C)

  Mathilde Duny 1 , Guillaume Jubeaux 2 , Olivier Geffard 2 , Aurélie Cortéjade 1 , Laure Wiest 3 , Mickaël Nicolas 1 , Emmanuelle Vulliet 3 1 CARSO-LSEHL, Vénissieux, France (mduny@groupecarso.com) 2 BIOMAE, Château-Gaillard, France 3 Institut des Sciences Analytiques -ISA, UMR 5280 CNRS Université Lyon1, Villeurbanne, France

	Context and objective
	Materials and methods
							Acetonitrile
							fraction
	Grinding
	+ water						ANALYSIS
	+ Acetonitrile + Hexane Doping (IS)	Vortex	citrate buffer) + 50 mg salts (MgSO4, NaCl, 100 mg QuEChERS	Vortex	Ultrasound	salts purification 100 mg dSPE Centrifugation	Centrifugation	LCMS8060NX triple quadrupole tandem, LC20ADXR, ESI source
			ammonium acetate
	Results
	Conclusions