Ultrafast Excited-State Dynamics of Carotenoids and the Role of the SX State.

Carotenoids are natural pigments with multiple roles in photosynthesis. They act as accessory pigments by absorbing light where chlorophyll absorption is low, and they quench the excitation energy of neighboring chlorophylls under high-light conditions. The function of carotenoids depends on their polyene-like structure, which controls their excited-state properties. After light absorption to their bright S2 state, carotenoids rapidly decay to the optically dark S1 state. However, ultrafast spectroscopy experiments have shown the signatures of another dark state, termed SX. Here we shed light on the ultrafast photophysics of lutein, a xanthophyll carotenoid, by explicitly simulating its nonadiabatic excited-state dynamics in solution. Our simulations confirm the involvement of SX in the relaxation toward S1 and reveal that it is formed through a change in the nature of the S2 state driven by the decrease in the bond length alternation coordinate of the carotenoid conjugated chain.

In photosynthetic organisms, carotenoids act as accessory pigments in light-harvesting (LH) complexes by absorbing in regions of the visible spectrum where the absorption of the main pigments (chlorophylls or bacteriochlorophylls) is not efficient. Moreover, carotenoids play a crucial role in photoprotective processes, and they allow a proper assembly of the LH complexes by stabilizing their structure.[1−3] The multiple roles that carotenoids play in LH complexes are made possible by their polyene-like structure and the specificity of their electronically excited states.[4−7] The light absorption of carotenoids involves the second singlet excited state (S2) because the transition to the lowest (S1) state is forbidden. By considering the approximate C2 symmetry of the polyene chains, S2 has B+ symmetry and is dominated by a single electron promotion from the HOMO (gerade, g) to the LUMO (ungerade, u) orbital, and S1 (∼A–) is mainly a HOMO → LUMO double excitation.[4,8] Upon photoexcitation, the generated S2 state decays in a few hundred femtoseconds to S1, which then relaxes back to S0 on a picosecond time scale.[8]

During the years, this two-state decay model has been enriched with additional (dark) states, which have been proposed to explain spectroscopic observations, and should lie in between the S2 and S1 states.[9−18] In particular, a state, commonly indicated as S, was initially detected as an intermediate in the S2–S1 decay of carotenoids in organic solvents, but then its signature was also observed in LH complexes.[16,19] Its energy position with respect to S2 and S1 was shown to be dependent on the number (N) of conjugated CC double bonds of the carotenoid, and the S features were assigned to an electronic state of B– symmetry. Further studies have also shown that the activation of S as a decay intermediate is also affected by the environment.[15]

The unique properties of polyenes have always attracted the interest of quantum chemical investigations, and many different levels of theory have been used to clarify the nature of their electronic states and simulate their photophysical properties.[20−35] On the contrary, much less has been done for what concerns the explicit simulation of the excited-state relaxation following the photoexcitation to the S2 state.[36−39] This lack of dynamic investigations has prevented a clear theoretical confirmation of the involvement of additional dark states in the excited-state decay.

Here, we focus on lutein (N = 10), which belongs to the xanthophyll class of carotenoids because it contains one hydroxyl group in each of its terminal rings. Lutein is among the most abundant carotenoids in nature because it is present in the major LH complex (LHCII) of photosystem II in plants.[40] Lutein is required for proper LHCII folding and organization.[41] In addition to extending the absorption window of the LHCII complex, lutein also plays a fundamental role in the so-called nonphotochemical quenching (NPQ), the photoprotective strategy used under excess light to dissipate the energy absorbed by chlorophylls into harmless heat.[3,42−45]

To unveil the ultrafast decay pathway of the S2 state of lutein and characterize the intermediate S state, we use nonadiabatic excited-state dynamics based on the mixed quantum-classical surface hopping (SH) method.[46] To the best of our knowledge, this kind of computational investigation for carotenoids is still missing in the literature. In fact, the molecular dimensions of the carotenoids, the number of correlated electrons, and the presence of states with different character make ab initio descriptions extremely costly and unsuitable for nonadiabatic dynamics.

Because the main features of the ultrafast dynamics of lutein have been measured to be similar in protein and in solution,[13,17,19] here we model a methanol solution. To account for the effects of the solvent molecules, we use a quantum mechanics/molecular mechanics (QM/MM) description (Figure ). The electronic structure of lutein is described by a semiempirical configuration interaction technique,[47,48] which allows a balanced description of all excited states on the same grounds. The semiempirical parameters are specifically optimized for the system under study. (See Methods for more details.)

Figure 1

(a) Quantum mechanics/molecular mechanics (QM/MM) system used in the simulations. The QM lutein is shown in a licorice representation, and the MM methanol molecules are represented by thin lines. (b) Molecular structure of lutein.

As a preliminary analysis, we computed the three lowest excited states (S1–S3) for an isolated lutein in its ground-state (S0) minimum geometry and compared it to the analogous but symmetric C20H22 polyene (Table ). In both molecules, the lowest excited state presents double excitation character (2A– for the symmetric system), and the second is mostly a HOMO–LUMO excitation (1B+). Notably, the third excited state (1B–) is only 0.22 eV higher in energy than the second. Although 1B– was not explicitly considered in the parametrization of our semiempirical method, this result is in agreement with previous DFT/MRCI calculations.[13] Moreover, ab initio DMRG/MRPT2 calculations have suggested that the adiabatic 1B– energy is lower than the 1B+ energy for carotenoids with eight or more double bonds.[35] A similar result is suggested by experimental Raman excitation profiles.[49] Other excited states are significantly higher in energy.

Table 1

Main Electronic Configurations in the Ground and the Three Low-Lying Singlet States of Lutein and the All-trans Linear Polyene with 10 Double Bonds (C20H22) at their S0 Minimum Geometry (Point Groups C1 for Lutein and C2 for C20H22), Computed with the R-AM1/FOMO-CASCI(6,6) Methoda

		weight		energy (eV)		oscilator strength
state	configuration	lutein	C20H22	lutein	C20H22	lutein	C20H22
S0(1Ag–)	···h, h	0.845	0.857	0.00	0.00
S1(2Ag–)	h, h → l, l	0.360	0.367	2.19	2.26	0.000	0.000
	h → l + 1	0.171	0.174
	h – 1 → l	0.161	0.158
S2(1Bu+)	h → l	0.826	0.826	2.62	2.60	2.865	3.247
S3(1Bu–)	h, h → l, l + 1	0.232	0.234	2.84	2.97	0.002	0.000
	h → l + 2	0.214	0.230
	h – 1, h → l, l	0.209	0.208
	h – 2 → l	0.173	0.170

The states are labeled according to their (pseudo)symmetry. For each configuration, the corresponding weight in the electronic wavefunction is reported. The complete list of CI coefficients, in the determinant basis, for each state of lutein is provided in Table S10. State energies relative to the ground state and the oscillator strengths are also reported. In the definition of the electronic configurations, we used h to indicate HOMO and l to indicate LUMO.

To confirm the robustness of these findings with respect to geometrical distortions, we sampled the Franck–Condon region of the solvated lutein through QM/MM S0 dynamics at room temperature. The vertical excitation energies calculated along the trajectory (Figure S5) confirm a close separation of the S2 and S3 states, whereas higher excited states remain well separated. While S2 is generally the brightest state, S3 acquires dipole strength (Figure S7), which suggests a mixing of these two states.

Starting from the S0 initial conditions, we propagated 200 SH trajectories by populating the excited states according to their radiative transition probability within an energy interval (see Methods for more details). Initially, S2 is the most populated state, but S3 is also significantly populated. In Figure a, we report the adiabatic state populations as functions of time obtained from the QM/MM surface-hopping simulations.

Figure 2

Adiabatic (panel a) and diabatic (panel b) state populations as functions of time obtained from the simulations of excited-state dynamics for lutein in methanol solution. The reported results are obtained by averaging over all trajectories and time intervals of 1 fs. In panel b, the fitting functions for the diabatic populations are also shown (dashed lines, see eqs S7–S9 in Section S3.3). The extracted time constants are τ2 = 21.8 fs and τ = 132.3 fs.

In our simulations, the photogenerated population of S2 decays to S1 within 200 fs. The main relaxation pathway toward S1 involves a direct transition (Table ). However, the population decay of S2 is clearly nonexponential, and strong population oscillations occur between S2 and S3 within the first 50 fs of the photoexcitation (Figure a). In all cases, the S1 state is directly populated by transitions from S2 (Tables and S7).

Table 2

Main Relaxation Pathways in Both the Adiabatic Basis and the Diabatic One Obtained in the Excited-State Simulations of Lutein in Methanol Solutiona

adiabatic basis	no. traj.	%
S2 → S1	70	35.0
S2 → S3 → S2 → S1	31	15.5
S2 → S3 → S2 → S3 → S2 → S1	14	7.0
S3 → S2 → S1	12	6.0
S2 → S1 → S0	11	5.5

For the diabatic basis, the pathways towards the 2A– and 1A– states are grouped together. For each identified pathway, the total number of surface-hopping trajectories and the corresponding percentages are also reported. Only the pathways with percentages above 5.0% are reported. The complete lists of pathways are reported in Tables S7 and S8.

To better characterize the nature of the electronic states involved in the dynamics and their changes along the relaxation, we introduced a diabatic description. The diabatic states are defined so that they maximally resemble the lowest adiabatic singlet states of lutein (S0–S3) at the S0 minimum geometry and are indicated using their (pseudo)symmetry: 1A–, 2A–, 1B+, and 1B– (see Table and Methods for more details on the diabatization).

The computed populations of diabatic states along the relaxation are reported in Figure b. Diabatic populations reveal a strikingly different picture with respect to adiabatic ones: here, in fact, the initially most populated 1B+ bright state rapidly transfers its population to the 1B– dark state. The 1B– state becomes lower in energy and can be identified with S2 in most of the trajectories. At variance with the adiabatic populations, analyzing diabatic states allows a clear identification of 1B– as intermediate in the relaxation. The 1B+ → 1B– transfer of population is ultrafast and occurs in the first ∼50 fs of our simulations. At longer times, 1B– starts transferring its population to the lower-lying 2A– and 1A– states, which are the most populated states at the end of the simulations (t = 200 fs). We see in Table that the main relaxation mechanism is 1B+ → 1B– → 2A– + 1A–.

We also identified several minor pathways involving multiple transitions back and forth between 1B+ and 1B– (Tables and S8). This indicates the oscillation of population between these two states as their crossing seam is traversed multiple times by the nuclear trajectory (see below). The 2A– + 1A– states are mainly populated by transitions from 1B–, and the 1B+ → 2A– + 1A– direct transfer of population represents a negligible route: the 1B+ → 2A– transition was observed in only 3 trajectories (out of 200), and no direct transition from 1B+ to 1A– was identified (Tables S5 and S8).

In Figure , we report the energies of the three low-lying diabatic excited states and the bond-length alternation (BLA) coordinate, averaged over all of the SH trajectories, as a function of time (see Figure S9 for the adiabatic energies). Here, the BLA is defined as the average difference between single- and double-bond lengths in the π-conjugated system. Throughout the ultrafast dynamics, the energies of the excited states undergo high-frequency fluctuations, which are generated by oscillations of the BLA. Notably, the 1B– energy is more sensitive than the 1B+ energy to the BLA change, which causes a swap between the two states at the beginning of the dynamics. These observations suggests that the entire relaxation to 2A– occurs along the BLA coordinate and is driven by the high-frequency C=C and C–C stretching modes of the polyene chain of lutein. The oscillations show a period of ∼20 fs, which corresponds to ∼1700 cm–1, compatible with the C=C frequencies, although slightly overestimated by the semiempirical method. Overall, the BLA coordinate undergoes a significant decrease during the ultrafast dynamics, indicating the shortening of the C–C bonds and the corresponding elongation of the conjugated C=C bonds.

Figure 3

Energies (eV) of the three low-lying diabatic excited states relative to the diabatic ground state (1A–) and bond-length alternation (BLA, Å) as functions of time. The reported results are obtained by averaging over all trajectories and time intervals of 1 fs.

To get a more detailed description of the relaxation pathways, we calculated the average energy gaps, electronic couplings, and values of BLA both at the Franck–Condon point (t = 0) and at the transitions between pairs of diabatic states. All of the values are reported in Table S4 of the SI.

Notably, the 1B+ ↔ 1B– transitions occur at geometries where the energy gap between the two states is nearly zero. This quasi-degeneracy condition is reached through a significant decrease in the BLA, namely, from ∼0.05 Å at the initial geometries to ∼0.025 Å at the transition. Such transition geometries are characterized by large values of the S2/S3 nonadiabatic coupling (Figure S10). On the other hand, the average electronic couplings between 1B+ and 1B– at their transition geometries are small and similar to the values at the Franck–Condon point. Contrary to the 1B+ ↔ 1B– transitions, the jumps from 1B– to 2A– occur at geometries where the energy gap between states is quite large, i.e., ∼0.64 eV on average, a value very close to the one calculated at the initial geometries. However, now the transitions occur with a larger electronic coupling: ∼100 meV compared to around 60 meV at the Franck–Condon point. There is also a further decrease in the BLA, which averages to 0.002 Å at the 1B– → 2A– transition geometries.

Comparing Figure a,b, one can notice that the 2A– and 1A– populations differ significantly from the populations of the S1 and S0 states, respectively, especially toward the end of the simulations. In addition, the two A– populations seem to follow the same time evolution. This is caused by a strong mixing of the 2A– and 1A– states when lutein approaches the S1 minimum. There, the S1 wave function, while dominated by the 2A– diabatic state, acquires an important contribution from 1A–. We quantified the electronic coupling between 2A– and 1A– at the geometries where S1 is the active state of the SH trajectories. The mean absolute value is ∼550 meV, much larger than the ∼80 meV at the Franck–Condon point (Table S6). The two A– states are also much closer in energy, with an average energy gap of ∼0.5 eV. The strong mixing makes the adiabatic S0 and S1 PESs much different from their diabatic counterparts. Given this mixing, in the diabatic analysis we have considered the two A– states together.

We quantified the characteristic times for the excited-state transitions by fitting the diabatic state populations in Figure b through a simple kinetic model. We assumed that the state populations obey the kinetic scheme 1B+ → 1B– → 2A– + 1A– with fixed rates. (See Section S3.3 for more details.) The fitting functions (eqs S7–S9) are shown in Figure b by dashed lines and well fit the state populations, with the exception of the oscillations at around 20 fs. We obtain the time constants τ2 = 21.8 fs for the 1B+ → 1B– transition and τ = 132.3 fs for the 1B– → 2A– + 1A– transition. These time constants also agree with the mean transition times reported in Table S4, confirming that the kinetic model well describes the population evolution.

Before summarizing the ultrafast evolution of lutein, we compare the results obtained so far with the excited-state dynamics of gas-phase lutein. (See Section S5 for the computational details.) Figure S12 shows the adiabatic and diabatic state populations as functions of time. The evolution is qualitatively very similar to the solution one shown in Figure . However, in the gas phase we observe more pronounced oscillations of the S2 and S3 populations during the first ∼50 fs and a slower excited-state decay to S1. The ultrafast relaxation mechanism identified in methanol is preserved in our gas-phase simulations. Specifically, in the adiabatic picture the S2 state exchanges population with S3 before decaying to S1, and in the diabatic basis, the 1B– state mediates the ultrafast decay 1B+ → 2A– + 1A–. The extracted lifetimes for the 1B+ and 1B– states in vacuum are τ2 = 49 fs and τ = 171 fs, respectively, which are slightly longer than the ones in methanol. The larger internal conversion rate in methanol, compared to the gas-phase simulations, can be mainly attributed to the faster 1B+ → 1B– transfer of population, which in turn is due to the larger 1B+/1B– average coupling, both at the starting geometries and at the 1B+ → 1B– transitions (compare Tables S4 and S13). Because the gas phase and methanol can be taken as two limits of low and high polarity, this result suggests a decrease in ultrafast relaxation times with increasing solvent polarity.

We can now characterize the detailed mechanism for the ultrafast excited-state dynamics of lutein. To this end, we make use of the potential energy curves calculated for the four low-lying states as functions of BLA (Figure ). These curves have been obtained by fitting the diabatic PESs along a relaxed scan on the S0 of lutein (Figure S8).

Figure 4

Energies (eV) of the four low-lying diabatic (solid lines) and adiabatic (dashed lines) states as functions of BLA (Å) obtained by fitting the diabatic PESs along a relaxed scan on the S0 state of lutein in vacuum (Figure S8). A schematic representation of the main relaxation pathway identified in the QM/MM surface hopping simulations is also shown. Note that the BLA axis is inverted for an easier interpretation of the mechanism.

Upon photoexcitation, the system, initially in the 1B+ bright state, is rapidly driven toward the crossing between the 1B+ and 1B– PESs and changes its electronic character while remaining in the S2 surface. Then, the nuclear wavepacket is divided among these two states, and the electronic population is split between the upper (S3) and lower (S2) surfaces. Multiple crossings can occur at this point between the two surfaces, causing the electronic populations to oscillate as observed in Figure . Finally, once the 1B– state is populated, the system relaxes to the S1/2A– state.

Comparing the diabatic to the adiabatic PESs (dashed lines in Figure ) also reveals the large electronic mixing between the two A– states, as the adiabatic S0 and S1 PESs strongly deviate from the diabatic counterparts as soon as the BLA decreases. This indicates sizable vibronic mixing along the BLA coordinate, which generates a widely avoided crossing between the two surfaces and alters the curvature of S0 and S1 PESs. Specifically, the S0 curvature is decreased, whereas the S1 curvature is increased. This effect causes a downshift in the vibrational frequencies of the C=C stretching mode in the S0 state and a frequency upshift in the S1 state.[50,51] We have monitored the frequency upshift by comparing the BLA power spectrum calculated with the S0 Born–Oppenheimer dynamics with the one calculated on those SH trajectories that remain in the S1 state for at least 175 fs (Figure S11). In the S1 state, the BLA power spectrum features a peak at substantially higher frequencies compared to the S0 dynamics.

Our results show that the 1B– state indeed mediates the S2 → S1 internal conversion by changing the electronic character of the S2 surface. We can then identify 1B– with the spectroscopically detected S dark state. This assignment also supports the main interpretations of the ultrafast spectroscopy of lutein in solution.[13,17]

In particular, Miki et al.[17] investigated the vibrational dynamics in the electronic excited states by pump-degenerate four-wave mixing and observed a frequency downshift of the C=C and C–C stretching modes in the first 200 fs after the photoexcitation. This effect was ascribed to the strong diabatic mixing between 1B+ and 1B– states. The reported values for the decay time of the S2 state are between 20 and 50 fs, depending on the solvent, while the extracted lifetime of S is ≲100 fs.[17] These decay times agree well with our simulations, suggesting that the observed ultrafast decay of lutein is indeed mediated by the 1B– state. In the combined experimental and theoretical investigation by Ostroumov et al.,[13] the authors observed large-amplitude oscillations in the transient absorption signals, which were attributed to electronic quantum beats caused by the coherent excitation of the strongly coupled 1B+ and 1B– states. The damping of the quantum beats was observed with decoherence times well below 100 fs (which can be related to the lifetime of 1B+), and the rise time for the absorption signal of S1 (∼2A–) was ∼600 fs.

Our extracted lifetimes for the S2 (1B+) and S (1B+) states of lutein, i.e., τ2 = 22 fs and τ = 132 fs, respectively, are also very close to the corresponding experimental times determined for other carotenoids with similar conjugation lengths.[9,10,12,14,15] In particular, the involvement of an intermediate S state in S2 → S1 was already observed by Polli et al.[12] for neurosporene (N = 9) and by Maiuri et al.[15] for spheroidene (N = 10). For clarity, all of the aforementioned experimental relaxation times, together with the ones determined in the present work for lutein, are reported in Table S9 of the SI.

Our simulations finally showed that high-frequency BLA modes are fundamental in driving the transition to the 1B–/S state. By analyzing the transient grating signal of β-carotene, Ghosh et al. have suggested that the S state is reached through a twist of one C=C double bond.[52] We could not observe any excited-state torsion around any double C=C bond that could drive the excited-state relaxation, at least in the S2 → S → S1 ultrafast transition studied here. Importantly, with our atomistic surface-hopping dynamics we have not selected a priori the relevant coordinates, which are instead a result of the simulations.

Our results can finally help clarify the involvement of S in light harvesting and in energy transfer to chlorophylls.[16,19] In particular, an S feature was observed by two-dimensional electronic spectroscopy in one of the lutein molecules bound to the major LHCII pigment–protein complex.[19] By analyzing the amplitude oscillations at different points on the two-dimensional map, Son et al. revealed the nonadiabatic formation of a dark state, assigned to S of the lutein bound to site L2.[19] This state is populated in less than 20 fs from S2 of the same lutein and transfers energy to the neighboring chlorophylls. Such a transfer forbids the observation of the S → S1 lifetime. However, the time scale of the nonadiabatic transition from S2 is again consistent with our results as well as with the measurements in solution. We therefore suggest that the S assignment to 1B– will also be relevant in light-harvesting complexes. Because the polarity of LHCII is intermediate between the gas phase and highly polar methanol, we expect similar excited-state dynamics to what we observed in our simulations. Steric constraints imposed by the protein may, however, influence the lutein excited states and their dynamics. Efforts are underway in our group to simulate the ultrafast dynamics of lutein in pigment–protein complexes.

In summary, the ultrafast decay pathway from bright state S2 to S1 has been simulated here for lutein in methanol solution using QM/MM nonadiabatic dynamics. These simulations have allowed the characterization of the nature of the intermediate state S, subject of controversial discussions in the literature.[53] The outcomes of our simulations indicate that the S state is formed through a change in nature of the S2 state, which switches from 1B+ to 1B– character within ∼20 fs after the photoexcitation. Then, from the intermediate state S (1B–) the system decays to S1 in about 130 fs. Unlike the S2 → S transition, which mainly occurs in the S2 potential energy surface, the S → S1 decay involves a transition between two different adiabatic PESs, namely, from the relaxed S2 state (∼1B–) to S1. The whole relaxation pathway S2 → S → S1 is driven by the decrease in the bond-length alternation coordinate of the carotenoid conjugated chain.

Methods

The electronic energies and wave functions of lutein were computed in a semiempirical AM1 framework using the floating occupation molecular orbital-configuration interaction (FOMO-CI) method.[47,48,54] The CI was of the complete active space (CAS) type, with an active space of six electrons in six molecular orbitals (MOs). A Gaussian width for floating occupation of 0.1 hartree was used. The active MOs were of the π type, located on the polyene chain of lutein. The standard AM1 parameters were reoptimized in order to reproduce the best experimental and computational data available on excitation energies, oscillator strengths, and geometrical parameters of lutein. More details about the parametrization can be found in the Supporting Information (Section S1). In the simulations of lutein photodynamics, the interaction with the solvent (methanol) was taken into account using a QM/MM approach with electrostatic embedding. In particular, the lutein molecule was placed in a spherical cluster of 920 methanol molecules, which were all treated at the MM level using the OPLS-AA force field[55] (Figure ). The starting conditions for the surface hopping (SH) nonadiabatic simulations were sampled from a ground-state thermal trajectory at 300 K, performed using the Bussi–Parrinello stochastic thermostat.[56,57] The QM/MM equilibration trajectory of lutein in methanol was run for 100 ps (Figures S3–S7 in Section S3.1). The last 80 ps of the thermalization were used to sample the initial conditions for the SH trajectories, which were selected by considering an excitation energy interval of 2.5 ± 0.1 eV and taking into account the radiative transition probability, according to the procedure outlined in ref (54). The SH calculations were performed using the local diabatization algorithm[48,58] and with a time step of 0.1 fs, as employed in the integration of both the nuclear degrees of freedom and the electronic ones. Quantum decoherence was approximately taken into account with the energy-based decoherence correction (EDC) algorithm,[59] setting the constant C to 0.1 hartree. A total of 200 SH trajectories (178 starting from S2 and 22 from S3) were propagated for 200 fs. The six low-lying singlet states were taken into account in the nonadiabatic dynamics. For each simulation time, the population of the ith adiabatic state was computed as the fraction of SH trajectories running on the ith adiabatic PES. To characterize the physical nature of the electronic states during the SH simulations, we applied a diabatization procedure previously devised in the framework of the FOMO-CI method.[60] More details about the diabatization can be found in the Supporting Information (Section S2).