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

Abstract

In some molecular systems, such as nucleobases, polyenes or the active ingredients of sunscreens, substantial amounts of photo-excitation energy are dissipated on a sub-picosecond time scale, raising questions such as: where does this energy go or among which degrees of freedom it is being distributed at such early times? Here we use transient absorption spectroscopy to track excitation energy dispersing from the optically accessible vibronic subsystem into the remaining vibrational subsystem of the solute and solvent. Monitoring the flow of energy during vibrational redistribution enables quantification of local molecular heating. Subsequent heat dissipation away from the solute molecule is characterized by classical thermodynamics and molecular dynamics simulations. Hence, we present a holistic approach that tracks the internal temperature and vibronic distribution from the act of photo-excitation to the restoration of the global equilibrium. Within this framework internal vibrational redistribution and vibrational cooling are emergent phenomena. We demonstrate the validity of the framework by examining a highly controversial example, carotenoids. We show that correctly accounting for the local temperature unambiguously explains their energetically and temporally congested spectral dynamics without the ad hoc postulation of additional ‘dark’ states. An immediate further application of this approach would be to monitor the excitation and thermal dynamics of pigment–protein systems.

---

1 Introduction

Photoexcitation of molecules in solution immediately launches a cascade of relaxation processes: solvent relaxation, internal conversion (IC), vibrational energy relaxation and any possible photo-chemical reactions, such as photo-isomerisation, to name but a few.^1,2 Even restricting to the simplest and seemingly trivial scenario of sequential electronic relaxation, determining and quantifying the energy relaxation pathways for large molecules is a daunting task due to the sheer number of participating degrees of freedom. One common simplification based on the separation of time scales is the division of vibrational relaxation into two distinct processes: internal vibrational redistribution (IVR) and vibrational cooling (VC). The former process represents the relaxation of a small sub-set of optically active vibrations following their initial excitation. Meanwhile, VC represents the flow of energy from the solute to the solvent molecules. It is typically assumed that the two processes are sequential, which underlies the separation itself. Namely, IVR leads to the redistribution of electronic excitation energy over the whole vibrational sub-system and effectively establishes a quasi-equilibrium corresponding to an elevated local temperature. Subsequently, the cooling of the thermal vibrations, VC, takes place as heat dissipation into the immediate environment. However, this separation is sometimes questionable and the interplay of these processes is important in numerous areas of chemical physics, ranging from heat transfer in proteins and nano-structures³ and thermal conduction though molecular junctions⁴ to chemical reaction mechanisms.⁵

Understanding the vibrational relaxation rates and pathways becomes important when they fall into the temporal domain of other processes,¹ such as bond breaking and formation during chemical reactions⁶ or energy transfer in pigment–protein complexes.⁷ There have been multiple studies aimed at both IVR and VC in either the first excited state, S₁, or the ground state, S₀. Earlier studies explored the principles of cooling on the relatively long-lived S₁ state in the laser dyes based on coumarins,^8,9 stilbenes,^10–12 malachite green,¹³ oxazine and rhodamine.¹⁴ On the other hand, in some systems the lifetime of S₁ state is remarkably short, the prime examples of such systems being nucleobases,¹⁵ solar screen molecules,¹⁶ polyenes¹⁷ and carotenoids.¹⁸ At the same time, the absorption spectra of these systems are lying in the visible or UV, which means that tremendous amounts of excitation energy are being disposed of on breathtakingly short time scales (as much as 2.1 eV in carotenoids¹⁹ or 4 eV in nucleobases²⁰ in times as short as 200 fs). This raises questions such as: where exactly does the disposed energy go? Or what are the effects that such ultra-fast disposal has on the molecules and their associated observables? Fundamental understanding of these questions is important because of the biological relevance of the mentioned molecules: from the role of nucleoside resilience to photodamage in early life evolution²¹ to carotenoid photo-protection of light-harvesting proteins²² to biomedical application of photo-protection by sunscreens.²³

Even though IVR and VC of molecules in solution have been subjects of extensive research,^24–27 the results to date are not general and conclusive. The relative simplicity of small molecules (2–3 atoms) allows for truly detailed mapping of internal energy relaxation.^27,28 Intermediate size molecules with a number of atoms between 10–30, such as various dyes, are comparatively well studied as well. In particular, several prominent studies of azulene^29,30 and trans-stilbene¹² report various mechanisms of IVR and VC, presenting a coarse-grained description in which the solvent is partitioned into the first solvation shell³¹ (FSS) and the bulk, Fig. 1. Within such a picture, VC is a composite process consisting of heat transfer from the solute to FSS, and subsequently from FSS to the bulk solvent. Lastly, there have been several analyses of relaxation mechanisms in truly large molecules, such as porphyrins³² or hemes in proteins,^13,33 however, a universal theory for IVR and VC in intermediate and large molecules in liquids is still lacking.¹ An appealing alternative is the non-equilibrium molecular dynamics (MD) approach.^34–36 It offers a route that avoids complicated and possibly excessive parametrization of specific rates. However, an additional task remains to connect such studies to particular experimental observables rather than just the characteristic time scales.

Fig. 1 Schematic representation of IVR and VC processes in solution. The event of photoexcitation instantaneously promotes the molecule to a non-equilibrium state and launches the IVR cascade within the high-frequency optically active modes. Subsequently, energy exchange between the molecule and the solvent takes place. In a coarse-grained picture, this exchange can be interpreted as energy transfer from the solute to the first solvation shell (FSS; shown in red) and heat dissipation from the FSS to the bulk solvent. The latter processes constitute VC and are characterized by time scales τ_ss (“solute-to-solvent”) and τ_SC (“shell cooling”), respectively. The inset summarizes the conceptual scheme of energy relaxation via successive layers.

The goal of this study is to formulate a general dynamical framework of optical excitation energy conversion into thermal energy residing both within the molecule and flowing out into the solvent. While the distinction between IVR and VC is traditionally based on the separation of their respective time scales, we do not resort to this ad hoc classification. Instead we quantum-mechanically model IVR as internal solute dynamics starting from a vibronic Hamiltonian. Conversely, VC is modeled classically by introducing the concept of local temperature. The solvent is partitioned into FSS and the bulk, Fig. 1, as suggested earlier, whereas an MD simulation is employed to determine the FSS and consequently estimate the time scale of heat diffusion from FSS to the bulk. Our unified treatment ensures that the system and local temperature evolutions are modeled self-consistently. Ultimately, we model the optical response of the solute in terms of transient absorption (TA) spectra, which elucidates the internal heating effects and identifies signatures of the IVR and VC contributions. We learn that IVR-related effects are present on much later timescales than one would expect, up to tens of picoseconds. As a particularly relevant example, the photodynamics of carotenoids are analyzed. Being members of the polyene family, carotenoids posses a complicated structure of short-lived electronic states.¹⁸ This has led to multiple conjectures regarding the observed spectroscopic signals.³⁷ Particularly, an absorptive feature termed S* is sometimes referred to as an optically dark electronic state.³⁷ Instead we shown how this feature comes about as a result of peculiar vibronic and temperature dynamics, and reproduce the evolution of carotenoid TA spectra from hundreds of femtoseconds to tens of picoseconds. This demonstrates that the presented treatment enables modeling of large polyatomic molecules in a way that disentangles the concomitant relaxations processes.

This manuscript is organized as follows: in Section 2 we present the theoretical foundations of the manuscript (following the summary shown in the inset of Fig. 1). We start with the description of relaxation within a vibronic subsystem (Subsection 2.1). We then discuss the energy redistribution within the intramolecular normal modes and subsequent transfer to the FSS and still further to the bulk solvent (Subsection 2.2). To quantify these processes we introduce the concept of effective local temperature. Next, we provide a prescription for monitoring the system dynamics by TA spectroscopy (Subsection 2.3), and introduce particular model and procedures for the implementation of the presented methodology (Subsection 2.4). In Section 3 we present the results, first focusing on an abstract model system to demonstrate the key features of the methodology (Subsection 3.1), then modeling the full dynamics of two carotenoid solutions as a particular application of interest (Subsection 3.2). Finally, in Section 4 Conclusions and outlook are presented.

2 Theory and model

2.1 Excitation relaxation in a vibronic system

In many large molecules several vibrational modes are optically active. Therefore we will consider the evolution of excitation in a vibronic subsystem, where the vibronic states are quantum states that require both the electronic and vibrational indices for their characterization. There are several approaches to model vibrational energy relaxation with explicit treatment of the vibronic states, ranging from rigorously quantum mechanical^38,39 to semi-phenomenological.^40,41 Here, we employ an intermediate scheme suitable for modeling optical response of strongly-pronounced vibronic transitions, termed the vibrational energy relaxation approach.⁴² It treats population dynamics via quantum mechanical master equations and models the TA spectra in a semi-phenomenological fashion, and it has been shown to reproduce numerous TA features of carotenoids in particular with high accuracy.^43,44

Let us outline this scheme by considering a manifold of electronic states of a molecule linearly coupled to its vibrational subsystem. The following approach rests upon the observation that typically only a limited set of high-frequency vibrational modes strongly couple to the optical transitions in molecules. These high-frequency modes require explicit treatment and corresponding optical transitions should be considered in terms of vibronic states. The remaining modes can be treated either classically or implicitly as a thermal bath. A good criterion for such a separation is the temperature: the set of high-frequency under-damped modes (ħω ≫ kT) can then be quantized as harmonic oscillators, while all the low-frequency (ħω < kT) modes constitute the thermal bath. To describe the relaxation processes within the given molecular system, in this paper we consider the following Hamiltonian of electronic states |i〉 with energies ε_i coupled to a finite set of harmonic oscillators, α = 1, …K, via mutually shifted potential energy surfaces:

Here, q_α and p_α denote the coordinate and momentum of the α-th vibrational mode (oscillator of mass m_α and frequency ω_α), and d_α^(ij) denotes the mutual displacement of this mode's potentials between states |i〉 and |j〉. The modes are coupled with strength c_iκ to an infinite set of harmonic oscillators of the bath with coordinates x_κ and momenta p_κ. Several approximations are made here. Firstly, we assume that the same normal modes are associated with each electronic state, which is not strictly true due to vibronic interactions and differences in potential energy surfaces (conf., Duschinsky rotations⁴⁵). To make the model simpler, the coupling strengths c_iκ are set identical for all modes associated with the same electronic state. Secondly, the off-diagonal elements of the electronic subspace are coupled to the oscillators of the bath with coupling strengths f_ijκ, which means that such coupling is vanishing at equilibrium and affects the system dynamics only via fluctuations. Such a treatment disregards the fact that some electronic states can be coupled non-adiabatically (i.e., conical intersections), these effects could be introduced separately.

By quantizing the high-frequency system modes, we recast the Hamiltonian into the vibronic basis spanned by the states |i_a〉 = |i〉|a₁〉ⁱ…|a_K〉ⁱ, where |a_α〉ⁱ are the wave-functions of the high-frequency modes on the i-th electronic state. We employ a tuple a = (a₁, …, a_α, …, a_K) as a bookkeeping devise, where the entries indicate the number of quanta of the α-th mode with energy spacing ω_α within the associated electronic state. We also refer to all the vibronic states (combinations of a) associated with the electronic state |i〉 as a vibronic manifold. The Hamiltonian eqn (1) recast into the vibronic basis is given in the ESI.† The dynamics of the vibronic states |i_a〉 are described by the equations of motion for the corresponding populations n_aⁱ(t). The equations are derived from the quantum Liouville equation under the secular and Markov approximations^42,46 and are given in the ESI.† Essentially, they are of the second order with respect to the system–bath couplings, c_iκ/f_ijκ, which are characterized by the spectral densities and . A particular form of the spectral density should be chosen to reflect the specific physical situation under consideration. In principle, the spectral densities include both the intramolecular degrees of freedom (low-frequency modes) and the solvent phonons that couple to the electronic or vibronic degrees of freedom. This effectively yields a master equation, with Fermi golden rule-type transition rates between the vibronic states. Namely, the rates are proportional to the spectral densities [capital script C](ω) probed at the energy gaps of the transitions. Additionally, the rates within a vibronic manifold are proportional to the respective vibrational quantum numbers, while the inter-manifold transitions are weighted by corresponding Franck–Condon factors, F_ia,ja′^α = ⁱ〈a_α|a_α′〉^j. We note that by virtue of detailed balance in the bath correlation functions, the presented model has readily inbuilt thermodynamics.

2.2 Local molecular heating and cooling

The majority of formal excitation energy relaxation studies consider a molecule or an assembly of molecules to be in contact with an infinite and infinitely fast heat reservoir. Indeed such a depiction is reasonable insofar as the IC and the associated vibrational relaxation are the limiting steps when compared to the inter-molecular energy redistribution. However, in some situations this assumption can be questionable. E.g., when ultrafast ICs are involved, the energy dissipation to the medium might not appear as instantaneous and the finite vibrational energy transfer between molecules should be considered. Additionally, a quasi-equilibrium of intra-molecular vibrations can be established in a situation where molecule-to-medium relaxation is the rate limiting step. This quasi-equilibrium gives rise to the effective local temperature.⁴⁷ Let us consider the vibrational subsystem of a molecule. This subsystem can be considered the immediate recipient of the energy given away by the electronic subsystem during IC: after all, a large molecule is its own immediate bath. In such a case, to quantify the redistribution of energy we turn to the internal energy of the molecule, which is given for a collection of harmonic oscillators in thermal equilibrium with the bath of temperature T by the well-known relationship:here, ω_n is the frequency of the n-th normal mode out of N in total (we set ħ = 1 throughout this paper), β is the inverse temperature and [scr E, script letter E]₀ formally denotes the zero-point energy. Eqn (2) provides the link between the energy received by the intra-molecular vibrations from the electronic subsystem and the internal molecular temperature. Namely, let us consider the energy of the electronic subsystem only, ε_el, and the vibrational subsystem of energy [scr E, script letter E](T_∞), where T_∞ is the constant homogeneous temperature of the bulk solvent. The electronic relaxation will eventually deliver energy ε_el to the vibrational subsystem, where it is shared among the normal modes as if these modes were in quasi-equilibrium with a thermostat of temperature T_loc: [scr E, script letter E](T_∞) + ε_el = [scr E, script letter E](T_loc). This effectively describes the local molecular heating. We note that in real systems some energy from the electronic subsystem likely relaxes directly to the solvent phonon modes. So realistically, under such “leaky” conditions, [scr E, script letter E](T_∞) + ε_el ≥ [scr E, script letter E](T_loc) should be expected. In this work such leaks are disregarded partially due to technical complications of evaluating the size of the effect, partially expecting that in larger molecules such contribution would be diminishing.

In the case of a vibronic system, the energy released during IVR also contributes to the local heating. We assume that the amount of energy dissipated by the vibronic subsystem gets instantaneously spread among the remaining modes thus establishing the Boltzmann quasi-equilibrium and the local temperature.⁴⁷ By this we also implicitly assume that the thermalization of the non-optical intramolecular modes is not hindered by any localization effects.^4,48 The evolution of the local temperature is governed both by IC and IVR and by the heat exchange between the molecule and the solvent, namely VC. The latter process leads to the equilibration with the bulk solvent. In the following we will review the respective in-coming and out-going terms of the equation for the local temperature. We note that typically the build-up of the local temperature is ignored,^21,29 and it is assumed that the starting point of evolution is [scr E, script letter E](T₀) = [scr E, script letter E](T_∞) + ε, where ε is either equal to the energy difference between energy delivered by the pumping laser, ε_pump, and the closest excited state or even ε ≈ ε_pump. Such a condition, however, implies instantaneous IC and IVR. Here, we seek to bridge the gap between the infinitesimal and finite IC + IVR in the calculation of temperature evolution.

In order to quantify the increase of the local temperature, we need to keep track of the energy stored in the vibronic subsystem at any given time, ε_el(t), as well as the energy delivered by the pumping laser by that time, ε_pump(t). The difference between the pumped-in energy and the stored energy is transferred to the vibrational subsystem. Thus, we assume that at every time-step δt the change of intramolecular vibrational energy equals

δ[scr E, script letter E] = [ε_pump(t + δt) − ε_el(t + δt)] − [ε_pump(t) − ε_el(t)]. (3)

Alternatively, δ[scr E, script letter E] can be viewed as the difference of the energy received (ε_pump(t + δt) − ε_pump(t)); vanishes outside the pumping pulse) and lost (ε_el(t + δt) − ε_el(t)) by the vibronic system during time-step δt. Eqn (2) can then be used to obtain the inverse dependence, T([scr E, script letter E]), and to calculate the increase of local temperature δ_{[scr E, script letter E]}T per time-step δt. Again, we note that if we included the solvent phonon modes, possibly contributing to the spectral densities of the previous subsection, the local temperature would be lower due to larger number of participating degrees of freedom.

The populations n_aⁱ(t) can be used to evaluate the δ[scr E, script letter E] in eqn (3) as follows. The total energy of the vibronic subsystem ε_el(t) is given by

where ε_a = ω₁a₁ + … + ω_Ka_K. We can also keep track of the energy actually delivered to the vibronic subsystem by the pumping pulse, which typically excites a limited number of vibronic levels a′ on manifold(s) i′. This amounts tohere, “pump only” signifies that populations n_a′^i′(t) are the solutions to the equations of motion with only the pumping term and no subsequent relaxation present. Thus, ε_pump(t) starts from zero and acquires a stationary value corresponding to the pumping frequency within the duration of the pumping time.

Should the molecule be isolated from the solvent, its local temperature would rise on the time scale corresponding to its excitation lifetime to the maximal value, T_max, determined by the laser pulse and the distribution function eqn (2). However, in real situations the molecule is in contact either with the solvent molecules or some kind of embedding matrix. In this paper we are considering the former situation, although the latter is of considerable biological interest in the case of pigments in proteins. The molecular cooling of solute molecules has been investigated to a different degree in several studies.^12,29 Here, we employ an approach analogous to the one given by Kovalenko et al., where the solvent is partitioned into the FSS and the bulk.¹² The local temperature, from hereon denoted as T, decreases due to the heat transfer to the FSS of the temperature T_S, which in turn cools down by heat diffusion to the remainder of the solvent of temperature T_∞. Such a process can be described, along with the earlier prescription for heating, by coupled differential equations:

Here, N_S is the number of solvent molecules within the FSS of radius R, γ is an empirical coupling strength, c_S/c is the heat capacity of the solvent/solute, and χ is the diffusivity of the solvent. The second term on the r.h.s. of eqn (7) describes the diffusive cooling of a sphere of radius R, and it is derived starting from the Fourier equation.¹² The parameters entering eqn (6) are of very different nature: while γ is purely phenomenological and depends on the specific solvent–solute interactions, the diffusivity can be evaluated from the other well known solvent parameters, as χ = κM/Cρ, where κ, M, C and ρ are the thermal conductivity, molecular mass, heat capacity and density, respectively. We denote the solute-to-solvent heat transfer time as τ_ss = (N_Sγ)⁻¹ and the FSS cooling, or “shell cooling”, time as τ_SC = R²/3χ. Taken together with the equations governing the vibronic population evolution, eqn (6) and (7) provide a method for the description of local heating and cooling.

2.3 Probing dynamics by transient absorption spectroscopy

We inspect the dynamics of the photoexcited system by means of TA spectroscopy.⁴⁹ The time- and frequency-resolved TA spectrum A(ω, t) can be given as the following decomposition:

A(ω, t) = A_A(ω, t) − A_E(ω, t) − A_A(ω, ∞). (8)

Here, A_A(ω, t) and A_E(ω, t) denote all the absorptive (i.e., both ground state and induced absorption; IA) and emissive (stimulated emission; SE) contributions. When the delays between the pumping and probing pulses are larger than the pulse durations, the A_A(ω, t) and A_E(ω, t) terms can be given by factorizing them into the temporal and frequency parts (coherent effects are neglected as well).⁴² Particularly, the A_A(ω, t) is given by:

Here, the response in the time domain is given by the populations of the vibronic levels, n_aⁱ(t), while the response in frequency domain is described by the spectral shape term I_aⁱ(ω):

where μ_ij is the electronic transition dipole moment of the i → j transition, and σ(ω, Δω^(ij)) determines the line-shape of the i → j transition (Gaussian or Lorentzian of FWHM Δω^(ij)). ω^(ij) = ε_i − ε_j is the 0–0 (purely electronic) energy gap and ω_aa′ = ε_a − ε_a′ denotes the energy gap between the vibronic states on different manifolds. Similarly, the SE contribution, A_E(ω, t), is given by the same eqn (9) and (10), except that now 0 < j < i (transitions to vibronic levels of lower electronic manifolds).

It is noteworthy that in principle the internal temperature manifests in TA spectra in several ways. Firstly, it affects the relaxation rates and thereby also governs the quasi-equilibrium population distribution, once it is established. Secondly, it changes the line-widths of the spectra due to increased participation of the low frequency modes.⁵⁰ The latter effect is incorporated via a phenomenological relation:^12,51

(Δω^(ij)(T))² = (Δω^(ij)(T_∞))² + η²(T − T_∞), (11)

which is based on the fact that the line-width of the transition can be partitioned into the temperature-dependent part and the part describing temperature-independent broadening mechanisms (such as static disorder⁵²). The broadening parameter η is to be obtained from an experiment for a particular molecule and solvent combination.

2.4 Protocol for modeling multiple relaxation layers

In the previous sections we described a general model of a molecule coupled to the environment. Here, we introduce a particular model considered further in the paper and particular procedures used to implement the presented framework. For the ease of discussion, in the following we will refer to the manifold of vibronic levels |i_a〉 as a state |S_i〉. We will consider three such states, |S₀〉, |S₁〉 and |S₂〉, and their vibronic manifolds of two high-frequency modes, hence, a = (a₁, a₂). Such modes associated with carbon–carbon double- and single-bond stretching (frequencies at around 1150 and 1520 cm⁻¹, respectively) are typical for polyenes and aromatics.⁵¹ The three states are shown schematically in Fig. 2a along with the explicit manifold of the vibronic levels in the electronic ground state, |S₀〉.

Fig. 2 (a) Schematic representation of a three electronic level system with the associated vibronic manifolds. Ground state vibronic manifold of the two high-frequency carbon–carbon stretching modes is explicitly shown on the right. (b) Connection between the parameters governing the evolution of the joint molecule-and-solvent system and the possible experimental observables. The bottom tier summarizes the model parameters: molecule/solvent parameters in the solid black frames are mostly model-independent; parameters in the dashed frame depend on a specific model, but can be verified independently; parameters in the red frames are specific to the presented model and can only be inferred in this specific set-up.

As described earlier, there are two kinds of relaxation processes in the considered system: the IC transitions between states |S_i〉 and vibrational relaxation within these states, namely IVR. In order to parametrize the relaxation rates, we assume the overdamped Brownian oscillator spectral density⁵² of the form [capital script C]_OBO(ω) = 2λωΛ/(ω² + Λ²), where λ is the reorganization energy and Λ is the damping rate. The latter parameter determines the overall shape of the spectral density. For simplicity it is assumed identical for all the processes and assigned a value of (163 fs)⁻¹ based on earlier modeling.⁴² The reorganization energies govern the strength of coupling to the bath and are given individually for a specific process in the model. The |S_i〉 → |S_j〉 IC rates are governed by reorganization energies λ_IC^ij, while vibrational relaxation within each of the three manifolds is governed by λ_IVRⁱ (identical for both modes as stated earlier). Together with more directly experimentally accessible parameters, such as energy gaps ω^(ij), frequencies ω_α and curve displacements d_α^(ij), they determine the time scales of particular vibronic transitions. Let τ_IC^ij and τ_IVRⁱ denote the sets of inverse rates for IC transitions and IVR within the manifolds, accordingly. Depending on the separation of these time scales, they may of may not directly correspond to the experimentally observed lifetimes of the states and vibrations, e.g. τ_S1 or τ_C=C. The hierarchy structure of model parameters and observables is illustrated in Fig. 2b: the bottom tier represents the parameters dependent on a particular model, such as given here; the middle tier corresponds to the internal parameters of the system that could in principle be obtained from an alternative model; the top tier represents the typical experimental observables, which correspond to the parameters of the middle tier either directly or indirectly, depending, e.g., on the separation of time scales.

Similarly to the molecular evolution, the principles of quantifying the solute–solvent heat exchange are summarized on the right hand side of Fig. 2b scheme. Two important parameters, the mean solvation number N_S, and FSS radius R, can be estimated from the MD simulations. This is done in the following way. From an MD trajectory of a single solute molecule in the solvent we calculate the radial distribution function³¹ (RDF), the first minimum of which gives the radius R. In the case of large molecules which are spatially distributed to the extent comparable or even larger than R, certain approximations need to be invoked. Firstly, the molecule is partitioned into segments. The segmentation is rather arbitrary, but the obvious criterion is to reduce the volume of interest below the sphere of radius R. Then both the segments and the solvent molecules are contracted onto their centers of masses. Subsequently, the RDFs are calculated for each individual segment. Lastly, after estimating each individual segment-radius the “solvation volume” of the whole molecule is constructed. The global FSS radius R, required in eqn (7), is approximated as the radius of a sphere covering the volume equivalent to the determined solvation volume. Note, that the accuracy of the value of τ_SC is limited by the deviation of the solvent volume from a sphere for which the Fourier equation was solved. A particular instance of such a procedure is described in the following section. The remaining bulk properties of the solvent, such as the heat diffusivity χ, are readily available in chemical datasets. By contrast, the detailed heat exchange between the system and its FSS falls outside the current level of theory, and the corresponding parameter γ is determined purely by fitting, as emphasized by the red rectangle in the scheme. From the basic bottom tier parameters we can estimate the solute-to-solvent heat transfer and shell cooling times, τ_ss and τ_SC, accordingly. However, the real observables of the top tier, such as the induced line-shape changes, also depend on the left part of the scheme, therefore both sides need to be treated self-consistently.

3 Results and discussion

3.1 Energy relaxation in a vibronic system and its tracing via transient absorption

We demonstrate the key features of the developed methodology by applying it to an arbitrary artificial system introduced in the previous subsection. The value of 10 000 cm⁻¹ is assumed for both energy gaps ω⁽¹⁰⁾ and ω⁽²¹⁾. We consider pumping the system by a 70 fs pulse centered at 20 000 cm⁻¹, which resonantly excites the lowest vibronic level of |S₂〉, i.e. |2₀₀〉, creating population n₀²(t). The reorganization energies λ²¹, λ¹, and curve displacements d⁽²¹⁾ are tuned in such a way as to yield both the |S₂〉 → |S₁〉 conversion and the vibrational relaxation within the |S₁〉 manifold on sub-picosecond time-scale. The choice of λ⁰ is dictated by similar values in real systems, as shown in next subsection. Reorganization energy λ¹⁰ provides an effective measure of the |S₁〉 decay rate. The specific parameter values underlying the studied regimes are given in the ESI.† Without loss of generality, the [scr E, script letter E](T) relationship, eqn (2), is taken from a particular molecule analyzed in the following subsection. The local heating parameters γ, N_S are also close to the values of particular molecules reported further, but otherwise completely arbitrary in order to highlight the model features. Namely, the shell consisting of 10 molecules (heat capacity identical to the model system, c = c_S) is set to cool on a time scale of τ_SC = 10 ps, while we set τ_ss = 2 ps.

Fig. 3a shows the dependence of the system evolution on the lifetime of |S₁〉 state, τ_S1. This parameter is evaluated by an exponential fit of the sum of vibronic level populations: . Evolutions of the lowest vibronic levels of C=C and C–C modes in the ground state, |0₁₀〉 and |0₀₁〉 respectively, are also shown as populations n⁰₁₀ and n⁰₀₁. The two panels of Fig. 3a correspond to the cases of τ_S1 = 10 and 1 ps (top and bottom); the dashed lines show the population equilibrium values at T = T_∞ = 300 K. For the time scale comparison, the vibrational relaxation rates based on the chosen λ⁰ are (2.34 ps)⁻¹ and (1.77 ps)⁻¹ for C=C and C–C modes, respectively. Let us firstly look at the τ_S1 = 10 ps case (top panel). The vibronic ground state levels are initially depopulated by the pumping pulse (note the logarithmic scale), then obtain some transient populations, which slowly relax sustained by the decaying |S₁〉. The relaxation of the |0₀₁〉 state is clearly not mono-exponential and reaches equilibrium within 50 ps. This suggests that starting with sub-10 ps |S₁〉 lifetime, the IC ceases to be the rate limiting step and the effects of populated ground state vibronic levels could be observed. Conversely, with extremely short-lived |S₁〉 state, the vibronic states both obtain considerable transient populations and outlast the excited state. They can clearly be seen biexponentially approaching the long-time equilibrium values, the fast exponent corresponding to the IVR times τ⁰_IVR for each mode, while the slow exponent matches τ_SC.

Fig. 3 Evolution of the artificial system after optical pumping. (a) Population evolutions. Evolutions of the lowest vibronic levels of C=C and C–C modes in the ground state (i.e., populations of states |0₁₀〉 and |0₀₁〉) are shown in red and blue, accordingly. (b) Molecular temperature evolution for different τ_ss values (τ_S1 = 1 ps).

The corresponding evolution of the molecular temperature for different values of the time scale τ_ss when τ_S1 = 1 ps are shown in Fig. 3b; the thickened line corresponds to the τ_ss value from panel a. For smaller τ_ss values it is evidently a biexponential process. As can be expected for a solution of two coupled rate equations, the maximal values of the transient species rise/fall with an increase/decrease of τ_ss. Also, with longer τ_ss the biexponential character of the temperature evolution diminishes. The initial part and the first exponent correspond to the establishment of equilibrium with the FSS, which heats up during this process (not shown) reaching its maximal temperature. The subsequent decay and the second exponent correspond to τ_SC (the solute and FSS cooling in tandem). The time scale of the initial part and the ratio of maximal temperatures are largely governed by the γ parameter, N_S and the heat capacity ratio. It is noteworthy that according to the distribution function used, the full conversion of 20 000 cm⁻¹ excitation energy into internal energy would result in molecular temperature of ∼555 K, while the value of 440 K, shown by the dashed line, corresponds to the conversion of 10 000 cm⁻¹, i.e. the energy gap covered by the IVR on the |S₁〉 state. The fact that for shorter τ_ss the maximal temperatures do not reach the latter value signifies that most of the temperature build-up is from the |S₂〉 → |S₁〉 IC, which is then efficiently dissipated towards the FSS. The fast time scale of this process underlines that a clear distinction between IVR and VC based on time scale alone may not be feasible. In most cases of molecular cooling analysis it is assumed that the formation of local temperature is nearly instantaneous, which is based on the idea that the intramolecular energy redistribution mediated by nonlinear interactions inside the molecule is much faster than the energy transfer from the molecule to its surroundings. Consequently, the temperature equivalent of the deposited excitation energy is considered to be the starting value of local temperature.^21,29,53,54 Here, we see that due to the details of the actual build-up of the temperature, its maximal values can be expected to be much smaller than the mentioned predictions based on instantaneous heating.

We now focus on how the described dynamics are reflected in terms of experimental observables. In optical experiments, molecular cooling is typically associated with the changes in the red edge of the TA spectra,^12,30 integrated peak intensity^12,55 or fluorescence narrowing.⁵⁶ Here, we focus on the signatures of molecular cooling in the ground state, as reflected in the spectral changes of the ground state bleach (GSB) contribution of the TA spectrum. TA spectra of the artificial system are shown in Fig. 4a–c. Due to the ultrafast (100–200 fs) decay of the |S₂〉 state and the energetic position of the |S₁〉 state, the frequency window around 20 000 cm⁻¹ at times longer than 1 ps is virtually free of contributions from the vibronic levels |2₀₀〉 and |1_a〉. This leaves only the GSB component and potentially some IA from higher ground state vibronic levels. Fig. 4a shows the simplest case corresponding to τ_S1 = 100 ps. It features a monotonic decrease (on the time scale of |S₁〉 lifetime) of the GSB signal, which is just the negative absorption spectrum (conf., shaded area in Fig. 5). By contrast, the case of τ_S1 = 10 ps having qualitatively nearly identical system dynamics, demonstrates entirely new features in the TA spectrum, Fig. 4b. At short times, t < τ_S1, there appears a significant positive feature on the red side of the GSB signal, while at long times, t > τ_S1, alternating positive and negative peaks, reminiscent of the absorption structure, are observed. For the τ_S1 = 1 ps case, Fig. 4c, the positive feature at ∼19 000 cm⁻¹ can now be clearly seen to rise up until 2 ps and then fall. At the same time, the spectrum above 20 000 cm⁻¹ is monotonically decaying until the purely negative signal “spills over” into positive values at around 4 ps; see especially the trough at ∼20 800 cm⁻¹. To better understand this behavior, let us look at only the IA from the ground vibronic levels, , shown in Fig. 5 along with the absorption spectrum—which technically is (shaded contour). The broad and relatively flat IA at 1 ps (black) originates mostly from transitions |0_a1>1a2>1〉 → |2_a〉. Around 2 ps (green) both vibrational modes acquire their largest transient populations in the first vibronic level, and the IA signal is mostly a composition of |0₁₀〉, |0₀₁〉, |0₀₀〉 → |2_a〉 transitions. Between 5 and 10 ps the IVR is over and quasi-equilibrium is established, n⁰_a (10 ps) ≈ n⁰_a (∞), as can be seen from both the population evolution and the faint red lines of Fig. 5. During the remaining time the only tangible change is due to the spectral line narrowing, which represents the remaining vibrations undergoing VC. The difference between the broad IA and the narrow absorption yields the typical oscillatory structure, which has been observed by Lenzer et al. as a difference spectrum of β-carotene solution at different ambient temperatures.⁵⁴

Fig. 4 TA spectra of the artificial system. Panels a to c correspond to the |S₁〉 state lifetimes of 100, 10 and 1 ps, accordingly. All the spectra are normalized to the maximum level of the spectrum at 1 ps for the τ_S1 = 100 ps parameter set (panel a). The faint gray lines correspond to spectra separated by 1 ps time step, while the faint red lines correspond to 10 ps time step, starting from 10 ps.

Fig. 5 IA component, A_A(ω, t), from the ground state of the artificial system. All parameters are the same as in Fig. 4c. The shaded contour shows A_A(ω, ∞), which is the absorption spectrum of the artificial system. The faint gray lines correspond to spectra separated by 1 ps time step, while the faint red lines correspond to 10 ps time step, starting from 10 ps.

3.2 Case study: carotenoids

In this subsection we apply the described strategy and model to a pair of real physical systems. Particularly relevant candidates for such a study are carotenoids, a family of ubiquitous natural pigments participating in vision, photosynthesis and photoprotection.⁵⁷ Their lowest excited state is both short-lived and optically dark, therefore the information about its properties is mostly accessible via TA.¹⁸ Spectral dynamics are typically complex with a number of signals not easily explained in terms of the quantum-chemical electronic structure. In the past decades various signals were being attributed to additional dark states,³⁷ which in some cases are now tenuously implicated in key biological functions.^58,59 One such example is the so-called S* signal. It has been reported in TA spectra of carotenoids already in 1995 and was interpreted as a hot ground state signal,⁶⁰ but later “rediscovered” in bacterial light-harvesting pigment–protein complexes and rather interpreted as a separate electronic state.⁶¹ Although this signal has been argued to represent vibronic transitions from S₁ and/or hot ground state,^43,54 this interpretation has not yet been unanimously accepted.^62,63 We demonstrate how the presented methodology exhaustively describes the carotenoid TA peculiarities, and elucidates the origin of the S* signal in carotenoids in a conclusive manner.

Particularly, we examine two carotenoids: canthaxanthin (Ctx) and rhodoxanthin (Rdx). Their absorption spectra in benzene along with the molecular structures are shown in Fig. 6a. The measured spectra are shown by thick faint lines, thin lines correspond to the model absorption spectra, I⁰_a(ω)n⁰_a(∞), determined by fitting. For carotenoids μ₀₁ ≈ 0 because of the inversion and so-called particle-hole/alternancy symmetries,⁶⁴ hence the linear absorption corresponds to the |S₀〉 → |S₂〉 transition only. The typical vibronic structure is well pronounced for Rdx, while for Ctx it is obscured by larger line-width. The fitting reproduces the spectra up to the second vibronic peak very well, but then falls below the experimental values. The latter feature is inherent for the harmonic model regardless of the specific implementation, either the “dressed” stick spectrum, as employed here, or the line-shape function via cumulant expansion.⁶⁵ The determined parameters are given in full in the ESI.†

Fig. 6 (a) Absorption spectra of Ctx and Rdx in benzene. Thick-faint lines correspond to the measurement, thin lines represent fitting. (b) Internal energy as function of temperature. (c) Schematic representation of the FSS as a cylindrical mesh. Highlighted are the segments of a carotenoid chosen for RDF calculation. The inset shows RDF of an end group of Ctx as obtained from the MD simulation. Red line corresponds to the Savitzky–Golay smoothing for the extraction of FSS radius.

In order to asses the energy distribution over the internal degrees of freedom in the carotenoids we perform the quantum chemical calculations of the normal modes. The details of the calculations along with the determined normal modes and their frequencies are provided in the ESI.† The subsequent dependence of the internal energy on the temperature [scr E, script letter E](T), eqn (2), is given in Fig. 6b. Evidently, the dependencies are nearly identical for both carotenoids. At 300 K they both posses [scr E, script letter E](T) ≈ 10 000 cm⁻¹, while the laser excitation into their |2₀₀〉 levels provides the excess energy of ∼20 000 cm⁻¹ (exact numbers are used in calculations). Should this excess energy be instantaneously converted into the equivalent internal vibrational energy, it would yield the maximal achievable temperature of T_max = 550 K. This temperature is sometimes assumed as a starting point of temperature evolution,⁵⁴ however, as demonstrated on the previous subsection (we used Ctx [scr E, script letter E](T) for the artificial system), this maximal temperature is unlikely to be achieved due to the dissipation channels which we will consider next. We parametrize the [scr E, script letter E](T) dependencies by analytical functions to obtain the inverted dependencies, T([scr E, script letter E]), for the dynamics calculations.

To parametrize eqn (6) and (7) from the statistics of FSS dynamics, we performed MD simulations of the carotenoids in benzene. The details of MD simulations are given in the ESI.† As mentioned earlier, for molecules as large as Ctx and Rdx (94 and 90 atoms, respectively) some coarse-graining prior to RDF calculations is necessary. We chose to divide the molecules into 6 segments shown in Fig. 6c and contract them onto their mass-weighted centers. In turn, each benzene was coarse-grained onto one particle at the center of mass. The RDFs for all the segments of both carotenoids are largely similar (conf., Fig. S4 in the ESI†), hence, only the RDF for the end group of Ctx is shown in Fig. 6c. Having applied some data smoothing, we can identify the minimum corresponding to the radius of the FSS. Subsequently we assume that the FFS forms a hemisphere-ended cylinder of length equal to the distance between the outermost coarse-grained segments, shown as the mesh in Fig. 6c. We then approximate the radius R in eqn (7) as the radius of a sphere equal in volume to the given cylinder (in which case , where r is the radius determined from RDF and l is the distance between the head groups). Taken with the diffusivity of benzene²⁹χ = 0.95 (10⁻⁷ m² s⁻¹), this yields shell cooling times of 6.0 ps and 6.45 ps for Rdx and Ctx, respectively.

Additionally, we approximate N_S as the number of solvent molecules in the vicinity of the solute. The average number of neighboring benzene molecules at distances from 0.5 nm to 0.6 nm increase from 5 to 10 for Rdx and from 6 to 14 for Ctx (see Fig. S5 in the ESI†). We picked the latter values for simulations. We note that while the radius R is combined with the fixed parameter χ to uniquely yield the time scale τ_SC, the number N_S is combined with the free parameter γ to yield the time scale τ_ss, hence τ_ss is still a fitting parameter and N_S merely acts as a constraint to extracting γ. The time scale of relaxation to the FSS was taken τ_ss = (N_Sγ)⁻¹ = 0.5 ps based on similar values determined for stilbene.¹² As this value could not be taken from an independent measurement or modeling, we attempted to verify its value indirectly. Namely, the extent to which the local temperature manifests in spectra depends on its magnitude, governed by τ_ss and τ_SC, and the parameter η viaeqn (11). Measuring the temperature dependence of Rdx absorption spectrum revealed that the value of η = 29.76 cm⁻¹ K^−1/2 is a reasonable estimate, which in turn used in the TA calculations together with τ_ss = 0.5 ps give a good fit. Details of the temperature-dependent measurement are given in the ESI.†

Now we apply the developed model and the determined parameters to reconstruct the TA spectra of the two carotenoids. Fig. 7 shows the spectra at late times, relevant for the discussion of the molecular cooling, while spectra at early times are provided in Fig. S2 in the ESI† (along with all the parameter values), to demonstrate the consistency of the methodology. Let us consider Ctx, Fig. 7a, as the first example. The measured spectra are shown by thick faint lines, the model corresponds to thin lines. Even though the spectral profile between 16 500 cm⁻¹ and 19 000 cm⁻¹ shows a mismatch at earlier times, the spectra are clearly decaying on at least two time scales. The IA from |S₁〉 at 16 000 cm⁻¹ decays on 4.9 ps time scale, whereas the decay of a signal around 18 000 cm⁻¹, termed S*, is visibly slower. An important feature accompanying the emergence of the separate S* signal is the narrowing of GSB.³⁷ This narrowing is qualitatively captured by the model very well (conf., Fig. S3 in the ESI†) and is a direct signature of the spectral line-width change. On the one hand, the effect of internal temperature on the line-width is essential to account for the extent of the process, as the presence of “non-thermal” ground state vibronic populations alone is insufficient. On the other hand, evidently these ground state vibronic levels decay at least on the same or slower timescale as the VC, because the signal in the GSB frequency window never crosses the zero line and is purely negative. The molecular temperature is found to decay on the time scale of 10–12 ps (depending on the fitting window). The extracted lifetime of |0₁₀〉 population is 8–10 ps. An even longer time constant of |0₁₀〉 decay in Ctx (15–20 ps) has been reported from the anti-Stokes Raman measurements.⁶⁶ These results provide a clear example of a system in which IVR is at least of the same or possibly even longer time scale as VC.

Fig. 7 (a, b) TA of Ctx and Rdx, accordingly. Thick-faint lines correspond to measurements, thin lines represent the model fitting. (c) TA of Rdx without the line-width dependence on temperature (η = 0).

By contrast, in the Rdx case, both |S₁〉 relaxation and IVR are considerably faster: τ_S1 = 0.9 ps, |0₁₀〉 lifetime is 3.5 ps. At 20 ps (violet line, Fig. 7b) the oscillatory structure is clearly observable and is identical to the red structure depicted in Fig. 4c. This shows that in Rdx both |S₁〉 → |S₀〉 IC and subsequent IVR are rapid enough to reveal the alternating structure typical for VC alone. This also reveals the dual nature of the S* signal at 17 500 cm⁻¹: at 4 ps (red line) this red-detuned satellite of GSB is mostly due to IA from higher vibronic levels on the ground state, at 8 ps (green line) it is a mixed result of higher vibronic levels and elevated molecular temperature, while at the latest times, e.g. 20 ps, this positive feature is mostly due to the line-shape differences between the locally “hot” (after pumping) and “cold” (no pumping) molecule. Such interpretation is further strongly supported by comparison of the spectra in Fig. 7b with the modeled TA assuming temperature-independent line-width (by setting η = 0 in eqn (11)) shown in Fig. 7c. The most striking difference is the lack of GSB narrowing, demonstrating that the temperature dependence of the line-width is a crucial feature of the underlying model rather than a mere correction. By contrast, the temperature dependence has little effect to the S* signal at times earlier than 20 ps. Even though the inclusion of temperature dependence visibly improves the overall line-shape of the peak, IVR alone is sufficient to explain its presence at these times.

Regarding the long-time evolution, of particular relevance and interest is the FSS radius, R, which is the scale parameter crucial for heat conduction, as for any diffusion process. This has been pointed out in the earliest attempts to describe VC.²⁹ Regardless, whether the FSS is modeled as a sphere,^12,29 point-source¹³ or a box,¹⁰ the scale parameter will be decisive in determining the order of τ_SC. The question of the model accuracy beyond that is rather open. For instance, we have disregarded the stability of FSS, which in principle is possible to asses from MD trajectories by estimating the mean dwelling time of a solvent molecule within the FSS. However, the interchange of molecules between FSS and the bulk would clearly reduce τ_SC by introducing convection as an additional relaxation channel. Hence, we can only expect the model to provide us the upper value of τ_SC. Interestingly, we find that the radii of FSS for various carotenoids/solvents do not differ much (conf., also results for zeaxanthin in THF in Fig. S4c in the ESI†), nor do the diffusivities of various solvents (they vary mostly within 0.8–1.2 × 10⁻⁷ m² s⁻¹).^10,29,55 This means that τ_SC is seemingly rather similar for all carotenoids. Moreover, in various studies of VC the correlation of cooling times and solvent diffusivities has been reported both present^10,55 and absent.⁵⁶ This shows that the decomposition of VC into solute-to-solvent and FSS-to-bulk steps is not straightforward. On the other hand, the γ parameter (or, conversely, τ_ss) for a given combination of solute and solvent is variable, yet there is no apparent single method to establish an exact microscopic dependence. For some large organic molecules the transfer to solute time is estimated to be of the order of 10 ps.^8,36 Taking all this into account, it could very well be that the reported variation of S* signal lifetimes from picoseconds¹⁹ to tens of picoseconds⁶⁷ can be due to the latter process.

As mentioned earlier, the current γ parameter is unknown from any direct complementary measurement or model. Therefore we complemented the TA experiment with temperature dependence of the ground state absorption spectra, as pointed out by Lenzer et al. The amplitude of the oscillatory structure provides a rather well established connection between γ and η parameters. The robustness of the model with respect to the variation of τ_SC, τ_ss and N_S is demonstrated in Fig. S6 in the ESI,† where we show that the determined values are optimal for the overall fitness of the model. The possible obstacles for relating the current η to the change of line-width of absorption spectrum with varied temperature Δω^(ij)(T), apart from the physical limitations (the thermal degradation of the sample), are two-fold. Firstly, uniformly heating the solution also changes some purely solvent properties, such as the refractive index and density of the solvent, leading to the peak shift and change in amplitude. Reversing these effects by data analysis might as well affect the retrieved Δω^(ij)(T) dependence. Secondly, the premises of the two experiments are rather different, and the question whether the temperatures of hot and cold bulk solvent yield identical Δω^(ij)(T) with respect to the molecular temperature alone remains open. Therefore we consider the uniform heating experiments only as complimentary means to constrain η and, subsequently, γ parameters. It would essentially be intriguing to obtain direct estimates of γ from non-equilibrium MD simulations.³⁶

4 Conclusions and outlook

In this study the explicit treatment of IVR is used to bridge the initial act of photoexcitation of a molecule and the eventual VC mediated by a solvation shell. We specifically analyzed the interplay of the characteristic time scales: τ_S1 (lifetime of the lowest excited state of the solute), τ_C=C (lifetimes of particular vibrational modes), τ_ss and τ_SC (solute-to-FSS and FSS-to-bulk heat transfer times, respectively). Particular combinations of these time scales determine which effects “ring out” in the observed TA spectra. When τ_S1 ceases being the limiting rate, the GSB signal is attributed to the hot ground state. Yet, one can evidently discern two types of “hot” state: associated with either the non-equilibrium distribution of vibronic populations or with an elevated local temperature (quasi-equilibrium of all the molecular vibrations). Both effects are sources of change in the GSB signal: IVR leads to the overall decay of the signal and its narrowing due to the vanishing vibronic IA contributions (Fig. 5), while the local temperature manifests via the line-width Δω^(ij)(T). If τ_C=C and τ_ss/τ_SC are comparable, then both the vibrational populations and the local temperature determine the appearance of GSB. Specifically, a positive shoulder, red-detuned with respect to the ground state absorption, readily stems from non-equilibrium vibronic populations (conf., Fig. 5 and 7c), however, it is enhanced by the thermal dependence of the line-width. On the other hand, given a further separation of timescales, τ_C=C < τ_ss/τ_SC, we have demonstrated—both in experiment and in the model—the emergence of GSB as alternating positive and negative values (conf., Fig. 4c and 7b). This oscillatory structure is a telltale sign of the elevated local temperature and immediately informs that IVR is mostly finished. Finally, our model shows that the initial build-up of the local temperature due to time dependence of IVR cannot be disregarded and the conjecture of ultra-fast IVR should be treated with caution.⁵⁶ The details of the initial heating determine the maximal achievable temperature (conf., Fig. 3b), which in turn is crucial in calibrating any internal “thermometer”, a physical observable indicative of the local temperature. The presented model could be nicely complemented (and tested) by measurements of the solvent temperature, and indeed such type of measurements for certain systems are already reported.^68,69

The applicability of the theoretical analysis was demonstrate by inspecting photo- and temperature dynamics in carotenoids, Ctx and Rdx.⁷⁰ The well resolved spectra and a very short-lived S₁ state (∼1 ps) of Rdx lead to the observation of interplay between the vibronic and thermal effects. After S₁ signal is gone, substantial populations of levels |0_a1>0a2>0〉 are clearly present, as the signal in the GSB region is purely negative. By contrast, at the latest times, the oscillatory structure emerges, which is associated with elevated local temperature of the molecule relaxed to the lowest ground state vibronic level |0₀₀〉. Conversely, Ctx never demonstrates the oscillatory spectral features. The main reason for this is the slower IVR, thus, Ctx is an example of a system where IVR and VC are concomitant rather than sequential. The obtained results also illustrate the dual origin of the S* signal as absorption from the hot ground state: before IVR is over, S* is a “hot” state in a dual sense of both non-equilibrium vibronic populations and increased local temperature; at times exceeding IVR, the S* signal indicates increased local temperature only. The complicated nature of the multiple relaxation processes of similar time scales explains the multi-exponential decay of S* in longer-chain carotenoids.¹⁹ We conclude that the presented research ultimately corroborates the interpretation of the S* signal in long carotenoids in terms of a hot ground state.^{42,43,54,60,71}

In summary, our model unifies energy relaxation processes within a molecule and into the surrounding medium, and reveals specific features of TA spectra, characteristic to IVR and VC processes. Most of the model parameters can be obtained by complementary measurements, therefore the model is highly consistent. Its predictive power is demonstrated by disentangling various relaxation pathways in molecules as large and complicated as carotenoids, where the coincidental signal lifetimes lead to the appearance of the absorptive features of the hot molecule in disguise as an excited electronic state. Hence, carotenoids are a suitable candidate “thermometer” for testing specific theories of solute–solvent interactions and intermolecular vibrational energy transfer. Such studies are important in understanding the cooling of pigment–protein complexes¹³ or the possible roles of hot states in the biological protein function. For example, a recent suggestion that the S* state triggers the photoactivation of the orange carotenoid protein (OCP) is a striking example.⁵⁹ The reported lifetime of the S* state in OCP is much longer than that of the S₁ state and is in line with our assignment to the hot ground state. Thus, the energy of a hot ground state of the carotenoid bound to OCP could be exploited to break the hydrogen bonds and initiate the OCP photocycle. Moreover, the role of solute-to-solvent energy transfer and local heating effects in the chemical reactivity and selectivity has been brought into question,³⁶ as it is becoming evident that in some cases the inclusion of non-instantaneous thermalization might be necessary to adequately capture the chemistry of reactive intermediates.^5,6

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank Radek Litvin for help with measurements of temperature-dependent absorption spectra. VB thanks prof. Mikas Vengris for inspiring discussions. This research utilized Queen Mary's Apocrita HPC facility, supported by QMUL Research-IT. VB and CD acknowledge the support from the Leverhulme Trust RPG-2015-337. TP thanks the Czech Science Foundation (16-10417S) for financial support. JH acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC-2089 as well as by the Austrian Science Fund (FWF): START project Y631-N27.

Notes and references

1. A. Rosspeintner, B. Lang and E. Vauthey, Annu. Rev. Phys. Chem., 2013, 64, 247–271.
2. T. Kumpulainen, B. Lang, A. Rosspeintner and E. Vauthey, Chem. Rev., 2017, 117, 10826–10939.
3. Proteins: Energy, Heat and Signal Flow, ed. D. M. Leitner and J. E. Straub, CRC Press, 2009.
4. D. M. Leitner, in Energy Relaxation and Thermal Transport in Molecules, ed. W. Andreoni and S. Yip, Springer International Publishing, 2018, pp. 1–22.
5. B. K. Carpenter, Chem. Rev., 2013, 113, 7265–7286.
6. B. K. Carpenter, J. N. Harvey and A. J. Orr-Ewing, J. Am. Chem. Soc., 2016, 138, 4695–4705.
7. M. Jakučionis, V. Chorosajev and D. Abramavicius, Chem. Phys., 2018, 515, 193–202.
8. T. Dahinten, J. Baier and A. Seilmeier, Chem. Phys., 1998, 232, 239–245.
9. K. Ohta, T. J. Kang, K. Tominaga and K. Yoshihara, Chem. Phys., 1999, 242, 103–114.
10. K. Iwata and H. Hamaguchi, J. Phys. Chem. A, 1997, 101, 632–637.
11. H. Okamoto, T. Nakabayashi and M. Tasumi, J. Raman Spectrosc., 2000, 31, 305–309.
12. S. A. Kovalenko, R. Schanz, H. Hennig and N. P. Ernsting, J. Chem. Phys., 2001, 15, 3256–3273.
13. T. Lian, B. Locke, Y. Kholodenko and R. M. Hochstrasser, J. Phys. Chem., 1994, 98, 11648–11656.
14. J.-Y. Liu, W.-H. Fan, K.-L. Han, W.-Q. Deng, D.-L. Xu and N.-Q. Lou, J. Phys. Chem. A, 2003, 107, 10857–10861.
15. C. E. Crespo-Hernández, B. Cohen, P. M. Hare and B. Kohler, Chem. Rev., 2004, 104, 1977–2019.
16. L. A. Baker, S. E. Greenough and V. G. Stavros, J. Phys. Chem. Lett., 2016, 7, 4655–4665.
17. N. A. Anderson, S. H. Pullen, L. A. Walker II, J. Shiang and R. J. Sension, J. Phys. Chem. A, 1998, 102, 10588–10598.
18. T. Polívka and V. Sundström, Chem. Rev., 2004, 104, 2021–2071.
19. H. Staleva, M. Zeeshan, P. Chábera, V. Partali, H.-R. Sliwka and T. Polívka, J. Phys. Chem. A, 2015, 119, 11304–11312.
20. C. T. Middleton, K. de La Harpe, C. Su, Y. K. Law, C. E. Crespo-Hernández and B. Kohler, Annu. Rev. Phys. Chem., 2009, 60, 217–239.
21. J.-M. L. Pecourt, J. Peon and B. Kohler, J. Am. Chem. Soc., 2001, 123, 10370–10378.
22. A. V. Ruban, M. P. Johnson and C. D. Duffy, Biochim. Biophys. Acta, 2012, 1817, 167–181.
23. N. D. N. Rodrigues, M. Staniforth and V. G. Stavros, Proc. R. Soc. A, 2016, 472.
24. T. Elsaesser and W. Kaiser, Annu. Rev. Phys. Chem., 1991, 42, 83–107.
25. J. C. Owrutsky, D. Raftery and R. M. Hochstrasser, Annu. Rev. Phys. Chem., 1994, 45, 519–555.
26. J. Assmann, M. Kling and B. Abel, Angew. Chem., Int. Ed., 2003, 42, 2226–2246.
27. J. L. Skinner, Theor. Chem. Acc., 2011, 128, 1.
28. M. P. Grubb, P. M. Coulter, H. J. B. Marroux, A. J. Orr-Ewing and M. N. R. Ashfold, Chem. Sci., 2017, 8, 3062–3069.
29. U. Sukowski, A. Seilmeier, T. Elsaesser and S. F. Fischer, J. Chem. Phys., 1990, 93, 4094–4101.
30. D. Schwarzer, J. Troe, M. Votsmeier and M. Zerezke, J. Chem. Phys., 1996, 105, 3121–3131.
31. D. Chandler, Introduction to Modern Statistical Mechanics, Oxford University Press, New York, 1987.
32. Y. Mizutani and T. Kitagawa, J. Mol. Liq., 2001, 90, 233–242.
33. M. Kondoh, M. Mizuno and Y. Mizutani, J. Phys. Chem. Lett., 2016, 7, 1950–1954.
34. G. Stock, Phys. Rev. Lett., 2009, 102, 118301.
35. S.-M. Park, P. H. Nguyen and G. Stock, J. Chem. Phys., 2009, 131, 184503.
36. S. Essafi and J. N. Harvey, J. Phys. Chem. A, 2018, 122, 3535–3540.
37. T. Polívka and V. Sundström, Chem. Phys. Lett., 2009, 477, 1–11.
38. S. H. Lin, J. Chem. Phys., 1974, 61, 3810–3820.
39. O. Kühn, V. May and M. Schreiber, J. Chem. Phys., 1994, 101, 10404–10415.
40. J.-P. Zhang, C.-H. Chen, Y. Koyama and H. Nagae, J. Phys. Chem. B, 1998, 102, 1632–1640.
41. T. M. Kardaś, B. Ratajska-Gadomska, A. Lapini, E. Ragnoni, R. Righini, M. Di Donato, P. Foggi and W. Gadomski, J. Chem. Phys., 2014, 140, 204312.
42. V. Balevičius Jr, A. Galestian Pour, J. Savolainen, C. N. Lincoln, V. Lukeš, E. Riedle, L. Valkunas, D. Abramavicius and J. Hauer, Phys. Chem. Chem. Phys., 2015, 17, 19491–19499.
43. V. Balevičius Jr, D. Abramavicius, T. Polívka, A. Galestian Pour and J. Hauer, J. Phys. Chem. Lett., 2016, 7, 3347–3352.
44. V. Balevičius Jr, C. N. Lincoln, D. Viola, G. Cerullo, J. Hauer and D. Abramavicius, Photosynth. Res., 2017, 135, 55–64.
45. Q. Peng, Y. Yi and Z. Shuai, J. Chem. Phys., 2007, 126, 114302.
46. L. Valkunas, D. Abramavicius and T. Mančal, Molecular Excitation Dynamics and Relaxation, Wiley-VCH, Weinheim, 2013.
47. D. Abramavicius, V. Chorosajev and L. Valkunas, Phys. Chem. Chem. Phys., 2018, 20, 21225–21240.
48. D. M. Leitner, Adv. Phys., 2015, 64, 445–517.
49. W. Demtröder, Laser Spectroscopy 2: Experimental Techniques, Springer-Verlag, Berlin, 5th edn, 2015.
50. J. J. Markham, Rev. Mod. Phys., 1959, 31, 956–989.
51. O. Braem, T. J. Penfold, A. Cannizzo and M. Chergui, Phys. Chem. Chem. Phys., 2012, 14, 3513–3519.
52. S. Mukamel, Principles of Nonlinear Optical Spectroscopy, Oxford University Press, New York, 1995.
53. S. Velate, E. O. Danilov and M. J. Rodgers, J. Phys. Chem. A, 2005, 109, 8969–8975.
54. T. Lenzer, F. Ehlers, M. Scholz, R. Oswald and K. Oum, Phys. Chem. Chem. Phys., 2010, 12, 8832–8883.
55. X. Tan, T. L. Gustafson, C. Lefumeux, G. Burdzinski, G. Buntinx and O. Poizat, J. Phys. Chem. A, 2002, 106, 3593–3598.
56. A. Pigliucci, G. Duvanel, L. M. Lawson Daku and E. Vauthey, J. Phys. Chem. A, 2007, 111, 6135–6145.
57. The Photochemistry of Carotenoids, ed. H. Frank, A. Young, G. Britton and R. Cogdell, Springer, Netherlands, 1999, vol. 8.
58. N. Liguori, P. Xu, I. H. M. van Stokkum, B. van Oort, Y. Lu, D. Karcher, R. Bock and R. Croce, Nat. Commun., 2017, 8, 1–9.
59. P. E. Konold, I. H. M. van Stokkum, F. Muzzopappa, A. Wilson, M. L. Groot, D. Kirilovsky and J. T. M. Kennis, J. Am. Chem. Soc., 2019, 141, 520–530.
60. P. O. Andersson and T. Gillbro, J. Chem. Phys., 1995, 103, 2509–2519.
61. C. C. Gradinaru, J. T. M. Kennis, E. Papagiannakis, I. H. M. van Stokkum, R. J. Cogdell, G. R. Fleming, R. A. Niederman and R. van Grondelle, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 2364–2369.
62. H. Hashimoto, C. Uragami, N. Yukihira, A. T. Gardiner and R. J. Cogdell, J. R. Soc., Interface, 2018, 15, 20180026.
63. B. Robert and H. A. Frank, Carotenoids: electronic states and biological functions, in Light Harvesting in Photosynthesis, ed. R. Croce, R. van Grondelle, H. van Amerongen and I. van Stokkum, CRC Press, 2018, ch., pp. 51–70.
64. T. Ritz, A. Damjanović and K. Schulten, ChemPhysChem, 2002, 3, 243–248.
65. C. Uragami, K. Saito, M. Yoshizawa, P. Molnár and H. Hashimoto, Arch. Biochem. Biophys., 2018, 650, 49–58.
66. T. Nakabayashi, H. Okamoto and M. Tasumi, J. Phys. Chem. A, 1997, 101, 3494–3500.
67. A. E. Jailaubekov, M. Vengris, S.-H. Song, T. Kusumoto, H. Hashimoto and D. S. Larsen, J. Phys. Chem. A, 2011, 115, 3905–3916.
68. S. C. Nguyen, J. P. Lomont, B. W. Caplins and C. B. Harris, J. Phys. Chem. Lett., 2014, 5, 2974–2978.
69. B. Dereka and E. Vauthey, Chem. Sci., 2017, 8, 5057–5066.
70. P. Chábera, M. Fuciman, P. Hřibek and T. Polívka, Phys. Chem. Chem. Phys., 2009, 11, 8795–8803.
71. F. Ehlers, M. Scholz, J. Schimpfhauser, J. Bienert, K. Oum and T. Lenzer, Phys. Chem. Chem. Phys., 2015, 17, 10478–17488.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9sc00410f

---