The full dynamics of energy relaxation in large organic molecules: from photo-excitation to solvent heating

In some molecular systems, such as nucleobases, polyenes or sunscreens, substantial amounts of photo-excitation energy are dissipated on a sub-picosecond time scale. Where does this energy go or among which degrees of freedom it is being distributed at such early times?


Equations of Motion
In this section we provide the equations of motion for the system described by the vibronic Hamiltonian as obtained from the initial Hamiltonian given in the main text: Here, the vibronic basis of states |i a = |i |a 1 i . . . |a N i is used, where |a α i denotes the wave-function of the α-th harmonic oscillator on the i-th electronic state, and the entries in the tuple a = (a 1 , . . . , a α , . . . , a N ) indicate the number of quanta of the α mode with energy spacing ω α within the associated electronic state. We defined an auxiliary vibrational energy matrix ε a = ω 1 a 1 + . . . + ω N a N , and (a α + 1) on the second line denotes a tuple that differs from a by the α-th entry, which is increased by one. On the third line F α ia, ja = i a α |a α j denotes the Franck-Condon factors. Broadly, the first term represents the static energies of the vibronic basis states, the second term generates transitions within the manifold of vibrational levels associated with the same electronic state (purely vibrational relaxation), while the third term generates transitions between vibrational manifolds belonging to different electronic states (internal conversion, IC). The final term represents the modes of the bath. The equations of motion for the populations of the vibronic states |i a , n i a , are derived from the quantum Liouville equation, by the second-order perturbation theory with respect to the system-bath coupling, as described in Ref. 1 . Essentially, they are Markovian rate equations of the following form: where the first term describes the internal conversion dynamics and the second one describes the vibrational relaxation: a ,a n i a + k ii+1 a ,a n i+1 Here, we have defined two kinds of rate constants: those between the vibrational levels of different electronic manifolds, k ji , and those between the neighboring vibrational levels of the same electronic manifold only, k ± . The latter are explicitly shown as downhill rates k + (ε) or the complementary uphill rates k − (ε), corresponding to the energy gap ε. These rates are obtained from the Fourier transform of the bath correlation function 1 , which is related to the corresponding spectral densities C c i (ε) and C f i j (ε) as where β = (k B T ) −1 is the inverse temperature. The downhill/uphill rate is then k ± (ω α ) = C c (±|ω α |). The rates k i j can be given by firstly defining the energy gap between the vibronic states: ∆ i j a,a = ω (i j) + ω aa , where ω (i j) = ε i − ε j is the purely electronic energy gap and ω aa = ε a − ε a denotes the energy gap between the vibronic states on different electronic manifolds. We then have the rates as k i> j a,a = C f (∆ i j a,a ) and k i< j a,a = C f (−∆ ji a,a ). We note that by the virtue of Eq. 5 the presented model has readily inbuilt thermodynamics. In order to describe the dynamics of the given vibronic system after photo-excitation we assume that the system is initially in the thermally populated vibronic ground state, i.e., n 0 and append to the given equations the following pumping terms: where n x y is the population of resonantly selected vibronic state(s) and Γ(t, τ p ) is the Gaussian function (τ p is the pump pulse duration).

Experimental details
Rhodoxanthin (Rdx) and canthaxanthin (Ctx) were purchased from Carotenature (Basel, Switzerland). Absorption and TA experiments were performed at room temperature. Temperature dependence. Temperature dependent absorption spectra of Rdx were measured in a Shimadzu UV-2600 spectrophotometer equipped with a temperature controlled cell holder TCC-240A (Shimadzu, Japan). The cell holder is an electronically controlled Peltier unit, temperature is regulated for both the sample and the reference cell. The instrument slit width was set to 2 nm. During the experiment, temperature was allowed to stabilize in the cell holder, and a 10-minute equilibration time before each measurement was used to ensure that the sample temperature was well stabilized prior to measurement of absorption spectrum at a particular temperature. Absorption spectra were measured in a 1 cm path length quartz cell.
Femtosecond transient absorption spectroscopy. Femtosecond TA data were recorded by a femtosecond spectrometer based on Integra-i (Quantronix) amplified laser system producing~130 fs pulses with 795 nm central wavelength at a repetition rate of 1 kHz. The amplifier output was divided into two paths. The first one was directed to the optical parametric amplifier TOPAS (Light Conversion) to generate tuneable excitation pulses. The second one to produce white-light continuum probe pulses in a 2 mm sapphire plate. The probe and the reference beams were brought to the slit of a spectrograph where it was dispersed onto a double photodiode array allowing measurements of transient spectra in a spectral window of~240 nm. In all measurements the mutual polarization of pump and probe beams was set to the magic angle (54.7 o ) by placing a polarization rotator in the pump beam. All measurements were carried out in a rotational cuvette consisting of two 1 mm quartz windows separated by a 1 mm Teflon spacer.

Temperature dependence of Rdx absorption spectrum
The absorption spectra shown in the main text were fitted first obtaining the values of parameters ω (02) , ∆ω (02) , d 1,2 as given in the  Table S1. The temperature dependent spectra of Rdx are shown in Fig. S1 (left panel). The spectra are normalized to the maximum peak, which eliminates the amplitude variation due to the change of the properties of the solvent, as discussed by Lenzer et al. 2

Figure S1
Temperature dependence of Rdx absorption spectrum. Panels left to right: spectral fits; thermal peak-shift; width dependence on temperature (green circles -results from spectral fits, black line -global fit of the dependence). spectra were then fitted using the same procedure varying only the 0-0 frequency, which shifts due to temperature dependence of the refractive index of the solute 2 , and the width. The fitting results are shown by dashed lines in the left panel of Fig. S1. For the sake of completeness, the extracted peak shift is given in the central panel of Fig. S1. Width dependence on temperature is shown in the rightmost panel: the width fitting from the experimental data results are shown by the circles; the line represents the global fitting of the results by the model The determined parameter values are ∆ω (i j) (T ∞ ) = 170 cm −1 and η = 29.76 cm −1 K −1/2 yields the parametrization closest to the experiment in the given temperature window.

Carotenoid parameters from TA
Modeling and fitting carotenoid TA spectra requires determining the remaining molecular parameters (in addition to the readily-known values of ω (02) , ∆ω (02) , d 1,2 from the absorption fitting). They can be grouped into "line-shape" parameters and "dynamic" parameters, as emphasized by black and blue notation in Table S1, respectively. The former group consists of the ω (1n) , ∆ω (1n) , d 1,2 to describe the IA from the |S 1 state, as well as the ratio of the GSB/IA amplitudes, |µ 02 | 2 /|µ 1n | 2 ; a normalization condition |µ 1n | 2 = 1 is used. The Stokes shift for |S 2 -to-|S 0 transition, δ ω 1,2 comprise the second group of parameters. While the reorganization energies λ i j IC mostly govern only the decay of the corresponding signal components, strictly speaking, λ i IVR and the displacements also influence the overall line-shape (as opposed to just individual vibronic transition σ (ω, ∆ω (i j) )). Namely, λ i IVR govern the narrowing of the spectra due to vibrational relaxation. The extent of such narrowing in turn depends on the "dynamic" displacements d 1,2 (conf., the "static" displacements d 1,2 ). They influence the global line-shape by determining which vibronic levels are the primary recipients of the excitation during IC (the larger the displacement -the higher vibronic levels get populated upon the IC first; the higher vibronic levels involved in an optical transition -the broader overall IA signal from that vibronic manifold). Some of the dynamic parameters can be verified experimentally, e.g. ω (12) , d 1,2 can be determined from the TA in the near infrared 3 . Here, we only resort to the typical values, also assuming ω (01) + ω (12) = ω (02) , and adjust the reorganization energies accordingly. The long-time fitting results are shown in the main text. Here, we additionally provide the early/intermediate times, Fig. S2, to demonstrate that the method uniformly captures all times from tens of femtoseconds (to tens of picoseconds in the main text).

GSB narrowing of Ctx: model vs experiment.
Although quantitative agreement between the model and the experiment is poorer for Ctx, it is possible to show extremely good qualitative agreement in the following way. Normalizing the experimental TA spectra to the maximum of the IA from the |S 1 state reveals two trends as demonstrated in the top panel of Fig. S3. Firstly, an increase of positive signal around 18000 cm −1 is observed in time; this is a good technique to demonstrate the existence of S * signal when it is not well separated from the |S 1 lifetime (which is the case for, e.g., Rdx). Secondly, the associated narrowing of GSB is clearly taking place. It has been observed earlier that such GSB narrowing accompanies the presence of S * signal. In this work we provide the underlying details of such relationship. This can be seen from the bottom panel demonstrating the corresponding modeling results.

Quantum-chemical calculation details
The normal modes of carotenoids for the distribution function were calculated after vacuum optimization using the Gaussian 09 package 4 with a B3LYP exchange correlation functional with the 6-311g(2d,p) Pople basis set. Normal mode frequencies for Rdx and Ctx   are given in Table S2.

Molecular dynamics
Molecular dynamics (MD) simulations have been performed on Ctx, Rdx in benzene and, additionally, zeaxanthin (Zea) in THF in order to explore the structure of typical solvation shells of carotenoids.

Simulation details
Sample preparation. Ctx, Rdx and Zea were optimized in vacuum using the Gaussian 09 package 4 with a B3LYP/6-311(2d,p) basis set. MD were conducted using AMBER 16 package 5 . A 70 x 70 x 70 Å 3 cubic periodic solvent boxes were created around the carotenoids, where the solvent could be benzene (for Ctx/Rdx) or THF (for Zea), using xleap program embedded in AMBER. The total vdw box size was around 156 x 132 x 132 Å 3 , depending on the particular solute and solvent combination. For all simulations General Amber Force Field (GAFF) 6 was applied for both carotenoids and solvents. The Coulomb cutoff was set to 8. No restraints were applied to the systems.
Equilibration/Relaxation. All simulations were minimized in energy with the steepest descend method for 500 steps followed by same steps of conjugate gradient method. The system was then relaxed and equilibrated in the constant density and temperature (NVT) ensemble and the constant pressure (NPT) ensemble with a time step of 2 femtosecond for 20 ps each, using Langevin thermostat with a target temperature of 300 K.
Production trajectories. The production runs were all performed with a 2 fs time step and a constant temperature of 300 K at an NVT ensemble to make analysis easier. To make sure the system is fully equilibrated, an additional 1 ns NPT run was conducted before the actual NVT production run.

Radial distribution function and solvation shells
In the three panels of Fig. S4 we show the RDFs for different carotenoid/solvent combinations: Ctx in benzene (left panel); Rdx in benzene (middle panel), zeaxanthin (Zea) in THF (right panel). Apparently, the results are similar for all segments for all carotenoidsolvent combination. Head groups tend to have larger radius in general, and RDF for the head group in Rdx is not satisfactorily resolved, yet the general trend is clearly observed.
To give a better impression of the shell structure, in Fig. S5 we show the distribution of the number of neighboring solvents within the given distance from the carotenoid (at this point, in order to avoid double counting, the whole molecule is considered rather than just segments). The counts are over different snapshots from the MD trajectory. The distributions of neighboring solvents within 5 and 6 angstroms are shown in green and orange, respectively. FIg. S6 shows the modeled Rdx TA spectra (thin full/dashed/dotted lines) for various parameter values around the determined ones (shown by thin full lines); thick faint lines show the experimental results for reference. From left to right τ SC , τ ss and N S are varied respectively.