Optical spectra of oligofurans. A theoretical approach to the transition energies, reorganization energies, and the vibronic activity

: There is experimental evidence of high vibronic activity that accompanies the strongly allowed transition between the ground state and the lowest electronic singlet excited state of oligofurans that contain 2,3, and 4 furan rings. The absorption and emission spectra of the three lowest oligofurans measured in liquid nitrogen temperature show distinct fine structures that are reproduced using the projection-based model of vibronic coupling (with Dushinsky rotation included) parameterized utilizing either DFT (with several different exchange-correlation functionals) or ab initio (CC2) quantum chemistry calculations. Using as reference the experimental data concerning the electronic absorption and fluorescence for the 8 lowest oligofurans we first analyze the performance of the exchange-correlation functionals for the electronic transition energies and the reorganization energies. Subsequently, we use the best functionals alongside with the CC2 method to explore how the reorganization energies are distributed among the totally symmetric vibrations, identify the normal modes that dominate in the fine structures present in the absorption and emission bands, and trace their evolution with the increasing number of rings in the oligofuran series. Confrontation of the simulated spectra with the experiment allows for verification of performance of the selected DFT functionals and the CC2 method.


Introduction
Oligofurans belong to the class of organic chromophores that, owing to accessibility of the low-energy excited states and their high processability, are viewed as potential building blocks of the organic electronic devices such as light emitting diodes or field effect transistors. Nowadays, the chromophores that are most commonly used for this purpose are based on oligoacenes, oligothiophenes, and their derivatives [2][3][4][5][6][7][8][9][10][11][12][13][14]. Recent experimental evidence suggests, however, that oligofurans may also be promising materials, owing to their higher rigidity and better solubility than the thiophene-based compounds, their high fluorescence quantum yield, as well as their biodegradable properties and availability of biorenewable sources [15][16][17][18][19]. It is believed that the replacing the highly polarizable sulfur atoms by oxygen is less advantageous for charge mobilities and explains the lower conductivity of oligofurans with respect to oligothiophenes [20]. It has been shown by Gidron at al. [21], however, that the tighter packing of the oligofuran films can more than compensate for the smaller size and less diffuse nature of the oxygen orbitals, leading to the charge mobilities that can surpass those observed for the oligothiophenes, which has opened the way to the development of new, environmentfriendly organic semiconductors.
Despite these encouraging reports, the number of both experimental and theoretical studies concerning excited states of oligofurans is still rather limited in comparison with the thoroughly studied oligothiophenes. One of the prominent experimental papers that gives an account of the photophysical properties of the lowest excited singlet state for the three smallest members of the oligofuran series (bifuran -2O, tertfuran -3O, and quaterfuran -4O) comes from the group of Becker[1] and provides detailed information on energetics and intensities of the transitions between the S0 and S1 states of the studied oligofurans, as well as the absorption and emission spectra measured in the liquid nitrogen temperatures (77K). The spectra display distinct fine structures for both absorption and emission, which indicates considerable vibronic coupling associated with the S0-S1 electronic transitions. Experimental UV/VIS spectra of higher oligothiophenes (5O-9O) were also reported [17,22], but owing to their low-resolution none but for the crudest vibronic features can be discerned. Nevertheless, they provide information on the vertical transition energies that we will use as a reference to select the optimum approximate version of the DFT exchange-correlation functional for studying the excited states of oligofurans in general, and for modeling the spectral vibronic activities in particular.
The vibronic activity is brought about by coupling of the nuclear and electronic degrees of freedom that occurs due to geometry changes induced by the electronic transitions. The primary effect -linear in the Taylor expansion -is the shift of the minimum of the potential energy surface (PES) for the nuclear motion. Quadratic terms: changes in vibrational frequencies (force constant modifications) or in normal mode identities (Dushinsky rotations) are usually of secondary importance and may pass unnoticed unless the spectra display sufficiently high resolution (like those obtained using the supersonic jet expansion or recorded in the Shpolsky matrices cooled in liquid helium). In this study, we aim to calculate the global reorganization energies associated with the investigated electronic transitions, reproduce the available experimental absorption and emission spectra of oligofurans, and identify the normal modes that contribute to the observed vibronic structures. To achieve these goals we employ the DFT/TDDFT methodology with several representative exchange-correlation functionals, as well as the ab initio approach using the approximate coupled clusters singles and doubles model (CC2) [23][24][25] that has already proven successful in similar applications [26][27][28][29][30][31] and currently can be regarded as the method of choice for theoretical analysis of spectroscopic properties of medium size (containing several tens of heavy atoms) organic systems. In order to select the best XC functionals for our purposes, as well as to further test the performance of CC2, we first calculate the electronic transition energies and the related reorganization energies for the oligofurans containing from 2 to 9 rings and compare them with the available experimental data [1,17]. The best performing functionals and the CC2 method are then employed in theoretical analysis of the vibronic spectral activities for the three lowest oligofurans. It is worth noting that oligofurans are planar in their ground states, in contrast to oligothiophenes, which adopt trans-gauche conformations between the adjacent rings. This is caused by the relatively low aromaticity of furan [32], owing to which the π-electrons, being less confined within the rings, are considerably delocalized over the whole oligomer chain (which increases the coupling between the rings and hinders the rotations along the inter-ring bonds) [18]. This structural regularity and considerable electron delocalization make oligofurans an ideal, even if demanding (owing to their size), playground for testing the performance of theoretical methods. Such tests are especially important for Kohn-Sham DFT, as it suffers from the arbitrariness of choosing the appropriate version of the exchange-correlation functional from the ever growing number of approximations. It is worth noting that selecting the functionals that can correctly reproduce the reorganization energies and vibronic activities is essential from the standpoint of organic electronics, as the nature and efficiency of the energy transport in the molecular crystals (and aggregates) are governed by the reorganization energy (relative to the charge/excitation transfer integral) and modulated by the same Franck-Condon overlap integrals (in the case of the coherent transfer) that determine the relative intensities of the vibronic features observed in the optical spectra [33][34][35].
One of the major sources of errors in DFT calculations, crucial for description of the extended, unsaturated systems, is the many-electron self-interaction error (SIE) [36,37], which manifests itself in the incorrect (exponential rather than r -1 ) long-distance behavior of the DFT exchange-correlation potentials that results in seriously underestimated values of Rydberg and charge-transfer states, overshot polarizabilities, and the tendency to overestimate electron delocalization shown by pure DFT functionals [38]. Since the HF method overly favors the localized structures, the so-called hybrid functionals that contain some portion of the non-local (or exact) exchange (calculated using the occupied Kohn-Sham orbitals in the same way as in the HF method) were designed as a remedy to these problems. Unfortunately, replacing the DFT exchange with the non-local one leaves unbalanced the SIE of the correlation part, which in pure DFT functionals is at least partly compensated by the SIE of the DFT exchange. Interestingly, the optimum fraction of the non-local exchange in the functional seems to depend not only on the particular property that a functional is dedicated to but also on the size of the molecule, especially if it contains extended systems of conjugated double (or triple) bonds.
Such behavior can be readily observed for the linear carbon chains (polyynes), which are one of the simplest extended systems of delocalized π-electrons. In this case, the calculated differences between the adjacent CC bonds (BLDs) were found to be closely related to the admixture of the non-local exchange in the DFT functional: the pure GGA functionals yielded strongly underestimated BLDs, while the hybrid functionals performed better, with the BHLYP (50% of the exact exchange) results being the closest to the experimental values [39]. The frequency of the dominant raman active mode (the socalled α mode, or the Я mode in the Effective Conjugated Coordinate Theory [40]), obtained using the DFT method were found to decrease with the growing chain length in accordance with the experimental observations[oligoeny]. The frequencies themselves, however, as well as the degree of their changes are functional-dependent, since they are functions of the BLDs [41,42]. An analogous situation was found for the HOMO-LUMO gap, which was strongly underestimated by the GGA functionals and overshot by those with large contents of the non-local exchange (e.g. BHLYP) [39]. The errors were not constant but depended on both the functional and the chain length. [43]. The quality of the BLDs obtained with the ab initio methods increased predictably with improvement of the model of the N-electron wavefunction, with the coupled cluster singles and doubles with perturbative correction from the triple amplitudes (CCSD(T)) yielding nearly accurate values for several short poliynes. [39,[44][45][46] Organic chromophores are another example of such challenging systems, and since their size in most cases prevents doing the correlated ab initio calculations, DFT is then the method of choice. Owing to the potential sensitivity of the DFT results to the choice of the approximate exchange-correlation functional, we initially study the behavior of the representative XC functionals for oligofurans with increasing chain length (up to 9O, containing 36 conjugated double bonds) and subsequently select the most suitable ones that will be used to model the experimental spectra of 2O, 3O, and 4O.

Computational details
Preparation of the theoretical spectra presented in this paper require molecular geometries optimized for the ground and excited electronic states, as well the vibrational analyses performed at both stationary points. Differences between the molecular geometries in the ground and the S1 electronic states projected onto the normal coordinates of either of these states yield the Franck-Condon (FC) factors that govern the vibronic activity accompanying the electronic transitions. The FC factors are calculated using the Dushin program of Reimers [47]. The theoretical spectra are simulated using the ezSpectrum program [48] that, apart from the spectra based on the FC factors only (assuming that the normal modes in both electronic states are virtually the same -which is the so-called parallel approximation) allows also for the inclusion of the normal mode mixing (Dushinsky rotations) in the spectra. Since it did not increase computational time, we have decided to take the Dushinsky effect into account even if the impact of the normal mode mixing is likely to be negligible for the spectra recorded in 77K, owing to their limited resolution.
The quantum chemistry calculations have been performed employing predominantly DFT and TDDFT methods. Owing to the sensitivity of DFT results to the composition of the exchange-correlation (XC) functional (especially strong for extended unsaturated systems), we employed several functionals with different admixtures of the nonlocal exchange: BP(86) -pure generalized gradient approximation (GGA) functional (0% of the non-local exchange) [49,50], TPSSh meta-GGA functional containing 10% of the nonlocal exchange [51,52], B3LYP -hybrid functional with 21% of the non-local exchange [53][54][55], PBE0 -a hybrid functional with 25% of the non-local exchange [56], MN15 -the most recent of the Minnesota meta-GGA hybrid functionals of Thrular et al., containing 44% of the non-local exchange [57], BHLYP -a hybrid functional with half of the exchange calculated in the HF-like manner [58], and three range-corrected functionals, in which the contribution of the non-local exchange varies with the distance from the reference electron: CAM-B3LYP (with 19-65% of the non-local exchange) [59], tuned-CAM-B3LYP -a version of CAM-B3LYP originally dedicated to description of the excited states of diarylethene derivatives (8-100% of the non-local exchange) [60], and ω-B97X (16-100% of the non-local exchange) [61]. To compare DFT with the ab initio calculations we have selected the coupled clusters singles and doubles model, in which the amplitudes for the singly excited configurations are computed rigorously, whereas the amplitudes for doubles are added perturbatively [23]. Such a simplification of the parent CCSD model reduces the computational cost by an order of magnitude (from N6 to N5), which allows for tackling systems as large as hexaphenyl [62] with a reasonable basis sets. Interestingly, owing to a fortuitous cancellation of errors the CC2 calculations often yield results that surpass those of the more complete (CCSD) model.
All the calculations have been performed with the def2-TZVPP basis set, which provides accuracy close to the complete basis set limit (especially for the DFT methods) at still reasonable computational cost. Preliminary calculations with the cc pVDZ basis set showed that reduction of the basis set brought minor changes to the calculated energies. It is thus reasonable to expect that further increase of the basis set (e.g. to QZVPP) would not lead to significant improvement of the results.
All the quantum chemistry calculations were done using Turbomole 7.5 package [63], apart from those that employed the meta-GGA hybrid functionals (TPSSh and MN15), for which we used Gaussian 16, [64] as Turbomole lacks the current density response for meta-GGA functionals. Without it some dependence of the results on the orbital rotation may occur. It is not likely to be strong, but we did observe some instability of the transition energies in the exploratory calculations. The errors were negligible as far as the transition energies were concerned, but could be quite significant for the energy differences, or, in the present context -for the reorganization energies.

Transition energies and reorganization energies
The adiabatic transition energies (Ead) between the ground state and the lowest excited singlet state are calculated as differences between the electronic energies corresponding to optimized geometries in both electronic states. The vertical absorption energy (E v.abs. ) is the difference between the energies of the S1 state and the ground state when both are calculated in the optimized geometry of the ground state, while the vertical emission Energy (E v.em. ) is the same difference but between the energies computed for the optimized geometry of the excited state. As follows, the excited state reorganization energy λ exc is given as E v.abs. -E ad. , whereas its ground state counterpart λ gr equals E ad. -E v.em. . Thus defined energies are schematically depicted in figure 1. E ad. corresponds to the experimentally discernible/observable 0-0 energy, provided that the zero-point energies in both electronic states are close to each other. E v.abs. / E v.em. is commonly attributed to the position (in the energy scale) of the maximum of the absorption/emission band. Strictly speaking, the vertical energies correspond to the centers of gravity of the respective transition bands. In the case of the room-temperature absorption spectra, the maxima are clearly visible, whereas estimation of the centers of gravity is difficult owing to the higher energy absorption bands that overlap with the signal of the So-S1 excitation. Therefore, for our study we chose to use directly the energies of the maxima given in the experimental papers. For the emission the situation is quite different, as the spectral bands are not 'contaminated' by spurious signals, but even at room temperature they display a distinct, vibronic structure, without a single dominating maximum. Therefore, the vertical emission energies for the studied oligofurans have been determined by calculating the centres of gravity of the respective fluorescence bands (which directly correspond to the vertical emission energies) [65]. Note that all the experimental spectra have been transferred from the wavelength to the wavenumber domain using the following conversion factors: The experimental values of the vertical absorption and emission energies are collected in tables 1 and 2, respectively, for comparison with the transition energies obtained from the quantum chemistry calculations. One should bear in mind, however, that the spectra were recorded in solutions (acetonitrile or dioxan) and are slightly red-shifted with respect to measurements in vacuum. Therefore, to serve as the reference for the quantum chemistry results, the experimental transition energies have been increased by 0.15 eV, which is a typical value of gas-to-solvent shift for organic chromophores in non-polar solvents. The calculated values of vertical absorption and emission energies are collected in tables 1 and 2, respectively. The electronic transition energies monotonously decrease as the number of furan rings in the oligomer grows, or -as the system of the delocalized π-electrons becomes more extended. This dependence is to a good approximation linear with respect to the reciprocal of the number of rings (1/n), but saturation has been observed observed for larger n (deviations from the linear trend begin at n>12 [66]) for several unsaturated oligomer series studied at various levels of theory [67][68][69][70]. Our present calculations for oligomers up to nonafuran confirm the linearity of the relation E(1/n) (see Figures 2 and 3), even if some deviations can be noticed for the pure DFT functional (BP), for which the dependence is slightly concave, and also for the ωB97X, for which it becomes convex for longer oligomers. Note, that also experimental energies seem to closely follow the linear trends. The figures 2 and 3 clearly show that not only the calculated energy values but also the rates of their changes with respect to 1/n strongly depend on the choice of the approximate exchange-correlation functionals. All of the functionals with the fixed contributions of the non-local exchange (from BP to BHLYP) generally underestimate the transition energies, and the deviations from the experimental values increase with the oligomer length. Only the BHLYP functional slightly overshoots the absorption energies for the three shortest oligofurans (2O-4O), yet for the longer ones the errors become negative, even if relatively small (not exceeding 0.23 eV for 9O). This agreement, however, is more accidental than systematic, as the slope of the E(1/n) dependence for this functional is considerably steeper than for the experimentally based trend (cf . table  3), so for longer oligomers the errors would most certainly increase. The reoptimized version of CAM-B3LYP (tuned-CAM-B3LYP) yields larger errors than the original version of this functional: the transition energies (both for absorption and emission) are underestimated for all oligofurans, with the errors growing for longer oligomers. It may seem surprising at first, because the long-range contribution of the nonlocal exchange reaches 100% in the tuned version of CAM-B3LYP, yet it behaves more like the classic hybrid functionals, such as MN15 or BHLYP, than the parent CAM-B3YP functional or ωB97X, which also has 100% contribution of the exact exchange in the long-range limit. Tuned CAM-B3LYP, however, contains considerably less non-local exchange at short distances (from the reference electron) -8% with respect to ωB97X -16%, or CAM-B3LYP -19%. Also the long-range limit is reached at considerably larger distances 1 . It thus appears that the short-range content of the non-local exchange and the size of the transition region are equally important for correct reproduction of the transition energies as the long-range composition of the exchange part of the functional. Interplay of these three factors may offer sufficient flexibility for a successful optimization of the range-corrected functional (e.g. CAM-B3LYP) for the specific case of the extended unsaturated systems (oligomers and other organic chromophores) and their excited states.
The transition energies obtained using the CC2 method have relatively large errors for both the absorption and emission energies. The energies also decrease more steeply with 1/n as compared to the experiment (for both absorption and emission). Such discrepancies have not been expected from the method that was used previously in several successful attempts of theoretical analyses of the optical spectra of unsaturated organic systems [26][27][28][29][30]. A closer inspection, however, shows that the respective errors of the CC2based absorption and emission energies are close to each other, the behavior that is mimicked only by two of the tested DFT functionals: MN15 and tuned-CAM-B3LYP (cf. figures S1 and S2 of the electronic supplement). Moreover, the slopes of E(1/n) functions fitted to the CC2 absorption and emission energies differ by the same amount (0.3) from one another as the slopes of the functions fitted to the experimental energies. Such performance, unique among the tested methods, is essential for modeling the optical spectra as it leads to the correct reproduction of the reorganization energies, which determine the width, shape, and -distributed among the totally symmetric normal modes -the fine structure of the spectral band. To obtain the reorganization energies based on the experimental data one needs to estimate the vertical absorption and emission energies, as well as locate the adiabatic energies. The former can be determined from the first moment of the spectral bands (or their "centers of gravity") [65], or even approximated as the position of their maxima if the bands lack the fine structure and are not overly asymmetric. Finding the adiabatic energies is more difficult, especially for the spectra recorded in room temperatures because the onset of the band does not necessarily correspond to the 0-0 line, which is often poorly resolved or overlapped by the hot bands. Bearing in mind that even minute errors of the adiabatic energies may lead to quite large errors of the reorganization energies (as they are given as differences between the two similar, and relatively large, values) we have decided to eliminate the adiabatic energies from the picture by taking as the experimental reference the halved differences between the vertical absorption and emission energies. Thus defined values correspond to the average of the excited state and the ground state reorganization energies. For the planar, organic molecules, such as oligofurans, for which the S0-S1 transitions are of ππ* origin (dominated by the HOMO-LUMO electronic excitations) these two reorganization energies can be expected not to differ much from one another. Browsing through the calculated values of the reorganization energies for all of the studied oligofurans (cf. table S2 of the electronic supplement) shows that this is indeed true.
The reorganization energies do not follow the linear trend with respect to the reciprocal number of rings so closely as the transition energies do. This, however, is quite understandable as the reorganization energies are strongly sensitive to even small inaccuracies of the transition energies (either calculated or estimated from the experiment). Also, for long oligomer chains, the calculated reorganization energies can be expected to saturate with the chain length, because the electronic transition-induced geometry changes tend to localize in the central region of the oligomer chain, leaving the terminal rings virtually intact. The extent of this region, however, strongly depends on the choice of XC functional. This localization does not seem to occur whatsoever for pure DFT functionals, while it is the most pronounced for the exchange-correlation functionals with large contributions of the exact exchange, which can be rationalized in terms of the inherent preference of the Hartree-Fock method for localized electronic structures. Analogous behavior was reported in the case of the oligothiophene cations containing up to 50 thiophene rings for the functionals with high large content of the non-local exchange [71]. For the oligofurans analyzed here, the chains are too short to appreciate the full extent of this effect, but the tendency can already be observed, as demonstrated in figure 3, which shows the striking difference between the pictures obtained using BP (pure GGA functional) and the ωB97X functional (in which the contribution of the exact exchange reaches 100% for long distances from the reference electron). The differences between the respective CC bonds in the excited and the ground electronic states are concentrated in the center of the chain for the ωB97X-based geometries, while for BP the differences are delocalized over the whole molecule (the central part is even somewhat less affected). There is no direct experimental evidence that might confirm the picture obtained with either of the two functionals, yet the overall improvement of the transition energies and the reorganization energies observed for the hybrid functionals suggests that it is these functionals that can grasp (at least qualitatively) the genuine physical effect. This conjecture is indirectly corroborated by an analogous mechanism operational e.g. in molecular crystals, in which the geometry distortions localized within the range of several nearest neighbors in the crystal lattice accompany the electronic excitation, leading to self-trapping of the excitons and to the formation of the polarons.
Given the above rationale, however, one might expect that the calculated reorganization energies would start to saturate with the number of rings first for the functionals that contain high contributions of the non-local exchange (BHLYP, ωB97X). Figure 4 shows, however, that the most pronounced saturation occurs for the pure DFT functional (BP) and those with small admixtures of the exact exchange (TPSSh, B3LYP). The dependence close to linear (no visible saturation) yielded PBE0, MN15, and BHLYP. The range corrected functionals, on the other hand, behave unexpectedly, as for them the calculated reorganization energies seem to decrease faster (instead of showing saturation) for long oligofurans. A deeper investigation into this effect falls beyond the scope of this study. One may speculate, however, that the varying contribution of the non-local exchange in the range-corrected functionals is at least partly responsible -it may have a beneficial influence on the transition energies, but it may also introduce subtle errors, which are revealed only when the difference values (such as the reorganization energies) are calculated.
The reorganization energies calculated in our study using the DFT functionals with the fixed contribution of the non-local exchange are strongly underestimated for all the studied oligofurans by the functionals, for which the contribution of the non-local exchange does not exceed 25% (figure 4). The BP functional is completely inadequate for computing the reorganization energies (and for further modeling of the spectra), as the errors reach nearly 70% for longer oligofurans. The TPSSh does perform better, but the 10% admixture of the exact exchange is insufficient to reach even semi-quantitative agreement with the experiment. B3LYP and PBE0 can be used 2O, or 3O, but even for the shortest oligofurans the errors are quite substantial. Further increase of the content of the nonlocal exchange in the functional leads to larger values of the reorganization energies, thus reducing the deviations from the experimental estimates. The reorganization energies obtained with MN015, containing 44 of the exact exchange turn out to be surprisingly accurate, with the errors less than 4% for all oligofurans except 2O, for which the calculated value is overestimated by 14%). BHLYP-based reorganization energies are considerably higher than those obtained with MN15, being thus overshot by about 20% with respect to the experiment. Apparently, the fixed 50% contribution of the non-local exchange is already too large for modeling of the optical spectra of oligofurans. From among the range-

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 19 October 2021
corrected functionals, the tuned-CAM-B3LYP has yielded the smallest errors, underestimating the reorganization energies by less than 15%. Using its general-purpose counterpart (CAM-B3LYP) has led to the reorganization energies that are by 13-23% larger (the errors decrease for longer oligomers) than the experimental estimates. The reorganization energies obtained with the ωB97X functional bear relatively high, positive errors for all the oligofurans under study.
The values of reorganization energies obtained with all DFT functionals considered here decrease too rapidly with respect to the experiment, with the slope that depends on the content of the non-local exchange. This is clearly visible in figure 5, as well as in table  4, which contains the slope parameters of the linear functions λ(1/n) fitted to the calculated and experimental values of the reorganization energies. The reorganization energies computed using CC2 on the other hand, seem to follow the trend that is nearly parallel to the experimental one. Hence, the deviations of the CC2 reorganization energies from experiment are not only small (1-12%), but also similar for all of the oligomers (up to 7O, as for 8O and 9O the calculations could not be converged for technical reasons). This feature of the CC2 method is likely to be the foundation of its good performance in modeling the optical spectra of organic chromophores. This method will therefore be the basis of our theoretical analysis of the vibronic structures in the optical spectra of 2O, 3O, and 4O. For comparison we will also demonstrate the theoretical spectra obtained using the results of the MN15 and the tuned CAM-B3LYP calculations, as those two functionals yielded the smallest errors of the reorganization energies.

Theoretical modeling of the optical spectra
Vibronic activity observed in the optical spectra of organic chromophores is often characterized in terms of Huang-Rhys (HR) factors Si=bi 2 /2. The Franck-Condon (FC) parameters, bi, can be evaluated from the formula: in which xA and xB are the 3N-dimensional vectors of the equilibrium Cartesian coordinates in the A-th and B-th states, M is the 3Nx3N diagonal matrix of the atomic masses, and is the 3N vector of the mass-weighted Cartesian displacements constituting the ith normal coordinate calculated for the A-th (final) electronic state. For symmetry reasons, the HR factors are non-zero only for the totally symmetric normal modes. The sum of the HR factors multiplied by the respective oscillation quanta is equal (within the accuracy of the calculated harmonic frequencies) to the reorganization energy for the final electronic state in the studied transition. The HR factors for 2O, 3O, and 4O are collected in table 5. For the sake of clarity the normal modes, whose contribution to the total reorganization energy is less than 1% are omitted (the complete sets of the FC parameters, HR factors, and the fraction of the reorganization energies for the three smallest oligofurans, calculated using all of the methods studied in this work can be found in table 3S of the electronic supplement). We have chosen the HR factors instead of the FC parameters since the former better discriminate the progression-forming modes. The HR factors are also direct measures of the relative intensity the fundamental replicas in the FC progressions with respect to the 0-0 line.
The largest values of the HR factors for all three oligofurans are associated with the normal mode around 1700 cm-1, in which the adjacent single and double CC bonds of the carbon skeleton of the molecule vibrate in antiphase. This is the same mode that is termed the Я mode in the framework of the Effective Conjugation Coordinate Theory [], according to which for the π-conjugated systems with strong electron-phonon coupling such a normal mode must manifest itself in the Raman spectrum as a single, very strong scattering at about 1500-1700 cm -1 . The scattering frequency decreases as the extent of the electron delocalization gets larger. The ECC theory was used to estimate the effective conjugation length in oligofurans, yielding the estimate of 14.5 repeat units [72]. It turns out that the Я mode dominates the vibronic activity for the analyzed oligofurans, being responsible for 55% of the total reorganization energy, and together with the satellite vibrations from the spectral region above 1425 cm-1 they carry from 59% (2O, emission) to 72% (4O, absorption) of the total reorganization energy. The other type of vibrations that have a noticeable impact on the vibronic activity are the ring breathing modes, located between 815 and 1035 cm 1, and contributing from 20% (4O, absorption) to 33% (2O, emission) to the total reorganization energy. The remaining normal modes cover for less than 10% of the reorganization energy and cannot be resolved in the spectra recorded in the liquid nitrogen temperature. Nonetheless, those vibrations cannot be neglected as their contributio is important for the overall width of the model spectra.
The theoretical spectra are formed by the overlapping vibronic bands, the energy of which (with respect to the 0-0 line) is given as: Where denote the numbers of vibrational quanta of the totally symmetric normal modes that are excited during the electronic transition. The vibronic bands are represented as a Lorenz-type shape function with the width of 250 cm -1 for the absorption spectra, and half of that for the fluorescence. The relative intensities of the individual vibronic bands are determined by the FC overlap integrals, available from the well-known recursion formulas [73] using the computed HR parameters. Table 5. The frequencies [in cm -1 ] and HR factors (italicized) calculated at the CC2/def2-TZVPP level of theory for three lowest members of the oligofuran series. The normal modes that contribute less than 1% to the total reorganization energy are omitted for clarity.

2O
3O 4O By virtue of the ezSpectrum software of Krylov et al. we have included the Duschinsky mixing of the normal modes induced by the electronic transition, even though this effect has turned out to be quite small for the spectra of medium resolution that we model in this work (figure S4 of the electronic supplement shows comparison of the selected spectra obtained with and without the Dushinsky mixing). The intensity throughout thus calculated emission bands is multiplied by the factor to account for the fact that the intensity distribution in the emission spectrum is proportional to the third power of the energy of the emitted light, according to the relation [74]: in which A is the Einstein coefficient for the spontaneous emission, is the permittivity of the vacuum, and μ represents the electric transition dipole moment. Analogous correction for the absorption spectra requires multiplication by , as the oscillator strength is proportional to the first power of the energy of the absorbed radiation: Here, is the mass of an electron, e -its charge, and c is the speed of light. The transition dipole moments were computed together with the transition energies, yet we have not used them as the absolute intensities of the experimental spectra were not reported in the study of Seixas de Melo at al [1]. Therefore we chose instead to normalize the theoretical spectra so that they match the total (integral) intensity of the experimental ones. Note that the vibrational frequencies have been multiplied by 0.95, which is an appropriate scaling factor for all of the three methods selected for this analysis [75,76].
The experimental absorption and emission spectra of 2O and their theoretical counterparts based on the vibronic parameters calculated with the MN15, tuned-CAM-B3LYP exchange-correlation functionals, as well as with the CC2 method, are presented in figure 6. The reorganization energies obtained with either of the selected theoretical methods The integral intensities of the theoretical spectra have been normalized to match the experiment. To ease the comparisons, the onsets of the computed spectra have been aligned to the experimental one. The stick spectrum has been obtained using the MN15 functional. are close to the experimental estimate: MN15 has yielded the value that is overshot by 13%, the error of the CC2 result is 8%, while that of tuned-CAM-B3LYP is -3%. A good match between the theoretical and experimental spectra can therefore be expected. The actual picture shows that while the frequencies of the vibronic features are reproduced almost quantitatively (after the appropriate scaling), some minor discrepancies between theory and experiment can be observed in the intensity domain. For the theoretical absorption spectra, the first two maxima are too low with respect to the third one, which might suggest that the HR parameters for the dominant normal mode (1776 cm -1 ) are somewhat overestimated. This deviation is more pronounced in the spectra based on the MN15 and CC2 results, while the best agreement with experimental absorption spectrum has been reached indeed for the tuned-CAM-B3LYP functional (which is consistent with the smallest error for the reorganization energy). The relative intensities of the first four vibronic features are here reproduced nearly perfectly (the third maximum is very slightly overshot). Some intensity deficiency visible at higher energies (4.5-4.7 eV) is probably an experimental artifact, owing to the considerable background intensity in that spectral region brought about by the long-wavelength tail of some higher energy absorption band, clearly visible as the steeply rising intensity above 4.8 eV.
Theoretical fluorescence spectra follow closely the experimental one in the low energy (tail) region, which indicates that the overall reorganization energy for the ground state (λgr) -governing the width of the spectral band has been indeed correctly reproduced. Deviations from the experiment are mainly observed in the higher energy part of the band (close to the onset): the intensity of the first maximum (corresponding to the 0-0 transition) is noticeably overshot, while the next two vibronic features have slightly insufficient intensities (especially for the spectrum based on the tuned-CAM-B3LYP results). It might seem that the HR parameters for the dominant vibrations (1687 cm -1 and 930 cm -1 ) are underestimated, thus overly accentuating the 0-0 line and underestimating the fundamental vibronic replicas. A closer look on the values of the calculated reorganization energies (table S2 of the electronic supplement) shows that λexc are larger than λgr by 7-9% for all the three theoretical methods selected in the previous section, which may affect the HR parameters for the progression forming modes and account for at least part of the observed discrepancies. The situation, however, is probably more complex, as the spectral region of the second and third maxima is shaped not only by the FC replicas of two dominant vibrations but also by numerous combination bands involving other normal modes, contributing to the background intensity (cf the stick spectrum in figure 6). Moreover, larger values of the HR factors for the two abovementioned normal modes would have to result in the elevated intensity of the low-energy tail of the calculated spectra and, in consequence, in deteriorated agreement with the experiment. It appears that some non-Condon mechanism influences the intensity pattern in the experimental emission spectrum. One of the possibilities that come to mind is the frequently occurring self-absorption that consumes some of the intensity of the first maximum, modifying the intensity ratio of the first and two following peaks. Whether this conjecture is true, however, would have to be verified in an independent experiment.
The calculated spectra of 3O are displayed in figure 7. The agreement with the experiment is nearly perfect for the frequencies, while the intensities show discrepancies that are of the same character as those for 2O. Again, the 0-0 line in the absorption spectrum is slightly underestimated, mostly for the MN15 and CC2 based spectra, while the high-energy tails of the theoretical spectra follow closely the experimental one. The tuned-CAM-B3LYP based spectrum best reproduces the experimental intensity ratio of the third maximum and the 0-0 line. The second maximum, however, is slightly overshot, while the high-energy tail carries insufficient intensity. The latter makes the band appear too narrow, which is in accordance with the reorganization energy underestimated by this functional by 11%. With the respective errors of 4% (MN15), or 1% (CC2), the intensity of the theoretical spectra based on those two methods are in nearly perfect agreement with the experimental one at high energies (perfectly reproduced overall width of the spectral band). The intensity of the 0-0 line, however, is again slightly too low. On the other hand, it is overestimated in the theoretical fluorescence spectra. Since the calculated values of λexc are larger than λgr by 7% for tuned-CAM-B3LYP), 8% for MN15, and 12% for CC2 (cf table S of the electronic supplement), it means that λexc are actually overestimated, and λgr -underestimated, by a few percent (especially when the average reorganization energy is practically exact). It might then explain the slightly underestimated intensity of the 0-0 line with respect to the 2nd, or the 3rd line in the absorption spectra as if the related HR factors were a little too large. This conjecture is corroborated by the second vibronic replica in the FC progression of the mode 1759 cm -1 (the 5 th maximum in the spectrum) being also a little overshot with respect to the experiment. The overestimated intensity of the 0-0 line in the fluorescence spectra and too low intensity of the third maximum seem to be consistent with the negative errors of the λgr. Whether these errors account for all of the discrepancies, or self-absorption was also operational as postulated for the case of 2O, cannot be resolved at the moment without additional experiments. Apart from the minor deviations from the experiment discussed above, the overall quality of the theoretical spectra can be regarded as quite satisfactory.
The theoretical spectra calculated for 4O show the best agreement with the experiment for both absorption and fluorescence. Both frequencies and intensities of all the vibronic features in the spectra are nearly quantitatively reproduced. The discrepancies observed for 2O and 3O in the intensity domain are strongly reduced or inexistent. For the tuned-CAM-B3LYP-based spectra, slight intensity deficiencies can be observed far from the onset, which reflects the fact that the reorganization energies are underestimated by this functional by 13%. It is worth mentioning, that the MN15 and CC2 reproduced the reorganization energy for 4O with quantitative accuracy (the errors of 2% and 1%, respectively), which shows in the excellent quality of the relevant theoretical spectra. It would require the experimental spectra recorded in the liquid helium temperature to further verify the performance of these methods, concerning the fine details of reorganization energy distribution among the normal modes, which are not resolved in the available spectra recorded at 77K. Preparations for such measurements are underway in our laboratory. At the moment, however, one must conclude that increasing the chain length 2O to 4O brings about little changes to the overall pattern of vibronic activity, formed by the two dominant types of molecular motions -the ring breathing modes and the CC stretches along with the carbon skeleton of the molecules. The maxima in the absorption spectra are more regularly spaced than in the fluorescence ones, owing to the frequency shifts of the progression forming modes due to the electronic excitation (the averaged frequencies of the ring breathing modes is down-shifted by 30-50 cm -1 , while its counterpart for the CC stretches increases by about 50 cm -1 ). As far as the evolution of vibronic activity is concerned, the two dominating types of normal modes behave differently. The HR factor for the CC stretching mode of about 1750 cm -1 remains virtually unchanged, while its counterparts for the ring breathing modes (800 1000 cm -1 ) undergo significant reduction: the total HR factor for these modes drops from 0.9 for 2O to 0.5 for 4O. The impact of this decrease on the actual spectra can be observed as a gradual reduction of the intensity of the second maximum with respect to both the third maximum and the 0-0 line. .

Conclusions
We have tested a representative set of exchange-correlation functional with different contributions of the non-local exchange (from 0 to 100%) for the reproduction of the vibronic activities that accompany electronic excitations of organic chromophores. Specifically, we focused on the vertical electronic transition energies that led in consequence to the reorganization energies, the key parameters that govern the degree of vibronic coupling, and determine the shape and width of the spectral bands.
We have observed that for oligofurans containing from 2 to 9 rings the vertical transition energies depend linearly on the reciprocal number of rings in the molecule. The slopes of the functions E(1/n) fitted to the calculated values appear to decrease with growing contribution of the exact exchange in the XC functional. The functional that performed best both for absorption and emission is CAM-B3LYP, for which the errors have been the smallest. Nonetheless, they were also size-dependent and for even longer oligofurans they could become substantial. Only the ωB97X functional yielded the absorption energies that were uniformly overestimated with the constant error of approximately 0.3 eV. The error of the emission energies reproduced by this functional (0.14 eV on average), however, slowly increased with the growing chain length. The other functional for which the longrange contribution of the non-local exchange reaches 100% performed more like the similarly to the classic hybrid-functionals with larger content of the nexact exchange (MN15, BHLYP). Apparently, the exact-exchange dominated long-range limit of the range-corrected functional is not the only factor that affects the excitation energies of highly delocalized systems. Others (like the short-range content of the exact-exchange, or the transition rate between the limiting regimes) appear to be equally important. There seems to be room for improvement -simultaneous optimization of these factors may lead to functionals that are better suited for reproducing accurately both the values of the electronic transition energies of the unsaturated oligomers (and more generally, of the organic chromophores) and of the trends of their changes with the size of the system.
The hybrid functionals with fixed amounts of the exact exchange have not performed so well as the range-corrected ones, as they have managed to yield small errors for only one or two oligomers at best, while the deviations for the experiment have been considerably larger. The CC2 method seems to work in a similar manner as the hybrid functionals with moderate admixtures of the non-local exchange. It is, however, the only method, which has yielded reorganization energies with relatively small (<12%) and approximately uniform deviations from the experimental estimates. For all of the DFT methods used in our study, these errors of the reorganization energies are distinctly size-dependent. They are also larger than for CC2, with the exception of the MN15 functional, which has performed remarkably well, yielding errors below 4% for all of the oligomers but for the smallest one (for 2O the calculated reorganization energy has been overestimated by 14%). MN15 is closely followed by tuned-CAM-B3LYP, for which the errors are slightly larger (3-15%) and negative. Nevertheless, for the reorganization energies, which as difference values is very sensitive to inaccuracies of the vertical transition energies (both theoretical and experimental), the performance of the CC2 method and the two abovementioned DFT functionals has to be regarded as very good. This assessment is confirmed by the high quality of the theoretical spectra obtained using the HR factors based on the quantum chemistry results. Besides nearly quantitative agreement with the experiment of the frequencies of the vibronic features observed in the absorption and emission bands, the overall bandwidths are also reproduced with excellent accuracy, which is consistent with the small errors of the calculated reorganization energies. Some discrepancies of the intensities, observed for the 0-0 line and the vibronic fundamental replicas of the progression forming modes are tentatively put down to a slight imbalance of the calculated reorganization energies for the excited states, which are larger than their counterparts for the ground state by at most 12% (depending on the oligomer and theoretical method involved). This conjecture, however, requires further investigation that would employ an alternative theoretical approach, in which the excited state is not accessed by the response theory, as it is done for CC2 and TDDFT. CASPT2 might be a good choice as it treats the ground state and the excited states on equal terms. At the same time, an alternative experiment would be beneficial, especially concerning the fluorescence of oligofurans, in order to assess (or exclude) the impact of self-absorption on the intensity of the 0-0 line. New spectra, if measured at the liquid helium temperature, would shed more light on details of the vibronic activity of the remaining totally symmetric normal modes (apart from the CC stretching and the ring breathing ones), providing new reference data for further testing of the theoretical methods on such excellent model systems as the oligofurans.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Table S1: Optimized cartesian coordinates for the studied oligofurans calculated using the selected theoretical methods (exchange-correlation DFT functionals) and CC2 with the def2-TZVPP basis set., Table S2: The excited state and the ground state reorganization energies calculated using the studied exchange correlation functionals with the def2-TZVPP basis set. Table S3: Vibrational frequencies, FC parameters, HR factors, and the related reorganization energies for the totally symmetric normal modes of 2O, 3O, and 4O, calculated using the studied exchange correlation functionals with the def2-TZVPP basis set.