Pyrene fluorescence after a high-energy electronic excitation exhibits a prominent band shoulder not present after excitation at low energies. The standard assignment of this shoulder as a non-Kasha emission from the second-excited state (S2) has been recently questioned. To elucidate this issue, we simulated the fluorescence of pyrene using two different theoretical approaches based on vertical convolution and nonadiabatic dynamics with nuclear ensembles. To conduct the necessary nonadiabatic dynamics simulations with high-lying electronic states and deal with fluorescence timescales of about 100 ns of this large molecule, we developed new computational protocols. The results from both approaches confirm that the band shoulder is, in fact, due to S2 emission. We show that the non-Kasha behavior is a dynamic-equilibrium effect not caused by a metastable S2 minimum. However, it requires considerable vibrational energy, which can only be achieved in collisionless regimes after transitions into highly excited states. This strict condition explains why the S2 emission was not observed in some experiments.

Kasha’s empirical rule1,2 states that “the emitting electronic level of a given multiplicity is the lowest excited level of that multiplicity.” This lack of emission from higher-lying states is simple to rationalize: internal conversion between excited states is much faster than photon emission. In other words, most chromophores with a singlet ground state end up in the minimum of the first excited state, S1, before they fluoresce. Nevertheless, fluorescence from higher-lying states, known as non-Kasha (or anti-Kasha) behavior, is well documented for several systems3–5 and is the basis of diverse applications exploring features such as tunable emission colors and improved luminescent quantum efficiency.6,7

Pyrene (see Fig. 1) is a prototypical chromophore with non-Kasha emission from a higher state.8–12 Consider, for instance, the results for low-pressure pyrene vapor reported by Baba et al.9 After a low-energy excitation into the S2 state (3.86 eV; see Fig. 1), pyrene fluorescence has a single band at 3.26 eV (380 nm) assigned to S1 emission. However, if a higher excited state is promoted with absorption at 5.32 eV, a shoulder at 3.65 eV (340 nm) appears in the emission spectrum. This new feature is commonly assigned to a non-Kasha S2 emission. Similar features and assignments were also reported for pyrene in a supersonic jet10 and other low-pressure vapor experiments.8,12

FIG. 1.

Top: absorption spectrum of pyrene in the gas phase (150 °C) from Ref. 13. Bottom: fluorescence of pyrene in the low-pressure vapor after excitation at 3.86 eV (left), 4.66 eV (middle), and 5.32 eV (right) from Ref. 9.

FIG. 1.

Top: absorption spectrum of pyrene in the gas phase (150 °C) from Ref. 13. Bottom: fluorescence of pyrene in the low-pressure vapor after excitation at 3.86 eV (left), 4.66 eV (middle), and 5.32 eV (right) from Ref. 9.

Close modal

This assignment of a non-Kasha emission, however, is not beyond controversy. del Valle and Catalán2 argued that the energy gap between S2 and S1 is too small to trap pyrene long enough to emit from S2. In fact, their experimental emission spectra of pyrene in 2-methylbutane (298 K) after excitations between 3.53 and 5.37 eV only show a single band corresponding to S1 emission.

Therefore, does pyrene emit from S2? Furthermore, if it does, why does the S2 feature not show up in the del Valle and Catalán’s results? We aim to explain the striking difference between these experimental findings using ab initio quantum chemical simulations, including nonadiabatic dynamics. However, to face the challenge of dealing with time scales of 100 ns of pyrene fluorescence,9 we needed to develop entirely new computational protocols, which are discussed in Sec. II. Alternative approaches to simulate non-Kasha systems have been recently discussed in Refs. 14 and 15 based on rate theory and Ref. 16 based on semiclassical simulations of non-Condon spectra.

We shall show that, yes, pyrene can emit from S2. Nevertheless, this non-Kasha behavior requires high-energy excitations and a collisionless regime so that the vibrational temperature of the molecule is high enough to keep S1 continuously pumping a small amount of population into S2. If such conditions are not satisfied, the S1 emission dominates the fluorescence. Thus, the non-Kasha behavior of pyrene results from a dynamic equilibrium between excited-state populations.

Excitation energies of pyrene were determined with the algebraic diagrammatic construction to the second-order [ADC(2)] using the SV(P) basis sets.17,18 The ground state was described with the Møller–Plesset perturbation theory to the second-order (MP2) with the same basis set. The S1 and S2 states were optimized at the ADC(2)/SV(P) level. Up to ten excited states were computed. No symmetry restrictions were imposed.

One thousand random nuclear geometries and velocities were sampled from a quantum harmonic oscillator Wigner distribution19 around the ground-state geometry to simulate the absorption spectrum and prepare the initial conditions for dynamics. The absorption cross section was computed with the nuclear ensemble approach (NEA)20 based on the ADC(2)/SV(P) excited-state energies and oscillator strengths of all these geometries. Dynamics simulations were initiated from two excitation windows, 5.7 ± 0.1 and 4.0 ± 0.1 eV. The ensemble points in these windows were stochastically screened to select those with the highest oscillator strength. In the high-energy window, we started 32 trajectories in the S7 and 20 in the S8 state, while in the low-energy window, we started 20 trajectories in the S2 state. Unlike conventional surface hopping simulations, in the present case, we do not require many trajectories but instead an extensive cumulative time summed over all trajectories (see Sec. II C).

Each trajectory was propagated for up to 2 ps with a 0.5 fs timestep employing ADC(2)/SV(P). Nonadiabatic effects between excited states were included with the fewest switches surface hopping21 via the local diabatization algorithm,22,23 with the simplified decay of mixing decoherence correction (α = 0.1 a.u.).24 After each hopping, the nuclear velocities were adjusted along the linear momentum direction to ensure total energy conservation. In the case of a frustrated hopping, the momentum direction was not modified. We employed the reduced kinetic energy reservoir described in Sec. II D to account for the non-extensivity of the velocity adjustment along the momentum direction.

Fluorescence spectra were simulated using two methods: (i) a vertical convolution of excitation energies and oscillator strengths at the excited-state minima and (ii) a nuclear ensemble built from dynamics results. Both are explained in Sec. II B.

All ADC(2) and MP2 calculations were done with the TURBOMOLE 7.3 program.25 Spectrum simulations, initial condition generation, and dynamics were carried out with Newton-X classical series (CS)26 interfaced to TURBOMOLE.

If we assume a distribution ρ1J of the excited population among the electronic states J (relative to S1), the total fluorescence spectrum is

ΓTOTE=Jρ1JΓJ0E,
(1)

where ΓJ→0 is the dimensionless differential emission rate from SJ to S0.20 Usually, only J = 1 matters to simulate fluorescence, but systems with non-Kasha behavior need more states for a complete description. The fluorescence intensity in terms of ΓTOT does not reflect the emission yield but rather the rate. The fluorescence lifetime can be estimated from the integral of the differential emission rate over the energy through20 

τFl=ΓTOTEdE1.
(2)

This section describes two approaches for computing these quantities. The first is the nuclear ensemble approach (NEA),20 in which the emission band is given by a sum of transitions from many different nuclear geometries, representing the nuclear wavefunction density in the source states. (The source states are those the molecule initially populates before the emission.) The second is the vertical convolution approach (VCA)27—a widely used spectrum simulation technique based on the vertical transitions at the equilibrium geometries convoluted with arbitrary band shape functions.

Within NEA, the differential emission rate is20 

ΓJ0NEAE=e22πmc3ε01HEhva1NpJlNpJΔE0J2RlJ×f0JRlJgLorentzEΔE0JRlJ,δJ,
(3)

where c is the speed of light, ε0 is the vacuum permittivity, and e and m are the electron charge and mass. H(Ehva) is the Heaviside function centered at the absorption energy a. It ensures that emission always occurs with energy gaps smaller than the excitation energy. ΔE0J is the energy gap, and f0J is the oscillator strength; both are calculated between SJ and S0. They are computed for each of the NpJ geometries in the nuclear ensemble for state SJ. The sum runs over all these geometries, too. gLorentz is the normalized Lorentzian line shape,

gLorentzε,2δ=1πδε2+δ2,
(4)

with an arbitrary line half-width δJ = δ/2.

The population distribution in NEA is defined by the number of points in the ensemble,

ρ1JE=NpJNT,
(5)

where

NT=JNpJ.
(6)

Inserting Eqs. (3) and (5) into Eq. (1) gives the NEA spectrum

ΓTOTNEAE=Jρ1JEΓJ0NEAE.
(7)

Alternatively, we can use VCA to simulate the fluorescence spectrum. In this case, we suppose that the sum of the many narrow Lorentzian peaks in NEA yields a broad Gaussian band. Thus, we apply the following approximation to Eq. (3):

1NpJlNpJΔE0J2RlJf0JRlJgLorentzEΔE0JRlJ,δJΔE0J2RJf0JRJgGaussianEΔE0JRJ,σJ,
(8)

where, in the right-side term, ΔE0J and f0J are computed at the minimum geometry RJ of SJ. gGaussian is the normalized Gaussian function,

gGaussianε,σ=12π1/2σe12ε2σ2.
(9)

The differential emission rate simplifies to

ΓJ0V CAE=e22π3/2mc3ε01HEhva×ΔE0J2RJf0JRJσJe12EΔE0JRJσJ2.
(10)

The band width σJ is an open parameter, which can be empirically adjusted.

The total emission in Eq. (1) requires the distribution ρ1J of excited-state populations. This quantity can be determined via rate theory (as discussed in Ref. 14) or via dynamics (as we discuss in Sec. II C). Nevertheless, these approaches are too involved to be used in a fast routine tool like VCA. Another possibility is to assume that the excited population is distributed according to a Boltzmann distribution,

ρ1JB=expΔE0JR1ΔE01R1kBTvib/×LexpΔE0LR1ΔE01R1kBTvib,
(11)

where the Boltzmann factors are computed at the S1 minimum (R1). This approximation assumes that the vibrational degrees of the molecule in S1 compose a thermal bath at temperature Tvib, which can populate higher excited states. In the collisionless regime, the vibrational temperature of the S1 state can be written as28 

Tvib=lnhνa+EZP0E1+EZP1hνa+EZP0E1EZP112EZP1NDFkB,
(12)

where E1 is the S1 state energy at the S1 minimum, EZPJ is the zero-point energy of state SJ, and NDF is the number of degrees of freedom of the molecule.

With the previous equations, the vertical convolution approximation for the emission is

ΓTOTV CAE=Jρ1JBΓJ0V CAE.
(13)

The NEA [Eq. (7)] and VCA [Eq. (13)] spectra can be directly compared. They do not require normalization. The main advantage of the VCA spectrum is that its computation is straightforward as it is just a Gaussian convolution of lines centered at the minimum of the source states. Nevertheless, the VCA spectrum does not contain information about vibronic couplings between states, and the bandwidth σJ is arbitrary.

The microcanonical formulation adopted here imposes another limitation on VCA. The population distribution in Eq. (11) is computed at the minimum of S1 (vertical transition), not at the minimum of the occupied state (adiabatic transition). This choice was made because a microcanonical system (which, by definition, does not count on a thermal reservoir) has different vibrational temperatures in each minimum. Thus, if we considered adiabatic transitions, a Boltzmann distribution would not be valid anymore. Later, we shall see that the Boltzmann distribution based on the vertical transitions, although it yields a useful first approximation for the fluorescence spectrum, does not quantitatively match the state distribution from dynamics simulations.

NEA spectrum simulations are significantly more involved because they require a nuclear ensemble construction and electronic-structure calculations for each nuclear geometry in the ensemble. (We will discuss how to build the ensemble for fluorescence spectra in Sec. II C.) Nonetheless, NEA spectra give absolute band shapes. Although the line half-width δJ is arbitrary, its value is chosen to be much narrower than the band width σJ. Thus, the band emerging from the NEA simulations reflects the actual band shape. Moreover, because the oscillator strengths are computed for every point in the ensemble, this approach includes information beyond the Condon approximation, recovering vibronic couplings.

As explained in Sec. II B, the NEA simulations require an ensemble of nuclear geometries representing the distribution of configurations in the source state. When simulating absorption spectra, sampling a Wigner distribution is convenient for building such an ensemble. However, for fluorescence simulations, the Wigner sampling would only be adequate for a regime in which the entire photon energy excess was dissipated into the environment, leaving the molecule at its zero-point energy level. This intense dissipation is not expected to occur in pyrene in a collisionless regime.10 Therefore, dynamics simulations are more adequate for building the ensemble in this case.

To build such an ensemble, we start the dynamics at the excited states, after excitation with energy a, and let the molecule relax through the state manifold. A microcanonical propagation ensures the collisionless regime, and we adopted trajectory surface hopping to guide the nonadiabatic relaxation. Nevertheless, there is one main concern we must deal with: fluorescence is a long-term process. In pyrene, for instance, it should take place within 100 ns,9 which is much longer than the state-of-the-art simulation capabilities. We can overcome this problem by noticing that nonadiabatic relaxation is ultrafast, occurring on the sub-picosecond scale. Then, supposing that the molecular motion is ergodic after this initial relaxation, we can propagate many independent, relatively short trajectories, discard the initial relaxation in each of them, and sample points from the remaining dynamics to compose the nuclear ensemble. Discarding the initial dynamics should erase the main memories from the Franck–Condon region. Propagating independent trajectories (instead of a single long one) should reduce auto-correlation in the ensemble. The technical details of this procedure are discussed in Sec. III B.

We used surface hopping for the nonadiabatic dynamics to build the nuclear ensembles. An issue not directly related to the spectrum simulations but crucial to the results’ accuracy is the velocity adjustment after hopping between states. The conventional way of making such adjustments yields artifacts when applied to large molecules,29,30 such as pyrene. We discuss in this subsection how this problem was solved.

Surface hopping propagates the molecular dynamics on a single potential energy surface. At each time step, one calculates the hopping probability between the currently occupied state and the others, and a stochastic process determines whether the molecule should remain in the same state or hop to another one. After a successful surface hopping, energy conservation is enforced by changing the kinetic energy to compensate for the potential energy variation due to the hopping. Thus, an equivalent amount of nuclear kinetic energy is added or removed depending on whether the system hops down to a lower electronic state or up to a higher state (the latter is called back hopping). If no kinetic energy change can ensure energy conservation, the hopping is frustrated, and the molecule remains in the same state.31 Following the semiclassical prescription for nonadiabatic forces,32–35 the kinetic energy change should be made by rescaling the nuclear velocities in the direction of the nonadiabatic coupling vector between the current and target electronic state.36,37 However, when this vector is unavailable [as in our simulations using ADC(2)], it is usual to rescale the velocities in the linear momentum direction.

In general, if a hopping event occurs from a state L to a state J, through a potential energy gap ΔELJ = EJEL, conservation of the total energy is ensured by constraining that variation of kinetic energy ΔΚLJ = KJKL equals minus the variation of potential energy,

ΔKLJ=ΔELJ.
(14)

In practice, this constraining is imposed by changing the velocity of the nucleus α as

vαJ=vαL+γLJuαMα,
(15)

where uα is a vector indicating the direction along which the change should be applied, and Mα is the nuclear mass. γLJ is determined as explained in Ref. 36. In the particular case that uα is the momentum direction MαvαL, Eq. (15) simplifies to

vαJ=κLJvαL,
(16)

where κLJ=1γLJ satisfies Eq. (14) if

κLJ=1ΔELJKL.
(17)

The hopping is allowed if κLJ is real, which means

KLΔELJ0.
(18)

Therefore, the kinetic energy of the entire system is considered when evaluating whether a back hopping is allowed or not, which may violate causality because energy from far away regions can be instantaneously used to allow the hopping.

This is not the case when rescaling is made in the nonadiabatic vector direction. To see it, note that for an arbitrary rescaling direction, the hopping is allowed if γLJ in Eq. (15) is real. This condition implies that it is allowed when the following inequality holds (see Ref. 36 for details):

12αuα2MααvαLuα2ΔELJ0.
(19)

Thus, inserting Eq. (14) in equality (19), the absolute amount of kinetic energy removed from the molecule after hopping from a lower to a higher state (L < J) must be

ΔKLJ12αuα2MααvαLuα2.
(20)

This amount is controlled by the orientation between v(L) and u. The adjustment in the nonadiabatic coupling vector direction implies that the kinetic energy reservoir used to ensure energy conservation is

ΔKLJ12αhLJ,α2MααvαLhLJ,α2,
(21)

which depends on the projection of v(L) on hLJ. Since hLJ is limited to the molecular region causing the state crossing, only a subset of atoms contributes to the kinetic energy reservoir. On the other hand, when the adjustment is made in the direction of the momentum, Inequality (20) simplifies to

ΔKLJKL,
(22)

and the kinetic energy of all atoms contributes to the reservoir, as we already know from Inequality (18).

In the present case of pyrene, direct evaluation of the nonadiabatic coupling vector is not possible due to the chosen ADC(2) electronic structure theory, and we have carried out the velocity adjustment in the direction of the nuclear momentum at hoppings. Since pyrene is a relatively large polyatomic molecule (26 atoms) and prepared with high-energy initial conditions (trajectories are initiated from excited states as high as S8), there is always a considerable amount of kinetic energy available in the reservoir to facilitate the back hoppings. Consequently, surface hopping propagation with velocity adjustment in the momentum direction over-allows back hoppings, thus, biasing the nuclear ensemble toward the excited states.

To overcome this artificial back-hopping excess, we propose reducing the available kinetic energy reservoir by dividing the nuclear kinetic energy by the number of degrees of freedom (NDF) when evaluating whether a back hopping can occur. This means replacing Inequality (18) by

KLNDFΔELJ0.
(23)

The rationale for this choice is the following. As mentioned, the semiclassical prescription32–35 implies that the correct kinetic energy reservoir should be restricted to atoms contributing to the nonadiabatic coupling vector [Inequality (21)]. Because this vector corresponds to one degree of freedom among NDF degrees,38 dividing the available kinetic by NDF should (statistically speaking) deliver a sensible reservoir, as it also refers only to one degree of freedom. In the case of pyrene with NDF = 72, this prescription significantly reduces the available kinetic energy, decreasing the number of back hoppings. It is triggered only to determine whether a back hopping will occur or will be frustrated; all the other features of the surface hopping propagation remain the same. An analogous reduction of the kinetic energy reservoir is typically used in quantum-mechanics/molecular mechanics (QM/MM) surface hopping dynamics, where only the kinetic energy of the QM atoms is considered.39 

Figure 2 illustrates the impact of the kinetic energy reservoir on a typical pyrene trajectory starting with the same initial conditions. In both cases, the molecule quickly relaxes to lower excited states within less than 100 fs. However, as expected, the conventional adjustment with the total kinetic energy [Fig. 2(a)] has many more back hoppings than the adjustment with the reduced kinetic energy [Fig. 2(b)]. For example, the trajectory in Fig. 2(a) attempted a back hopping from S1 to S2 at 99.5 fs. The potential energy gap was 0.22 eV. With the conventional kinetic energy reservoir, 3.68 eV of kinetic energy was available, and back hopping was allowed. With the reduced kinetic energy reservoir, the available kinetic energy would be only 3.68/72 = 0.05 eV, and the back hopping would be frustrated. Most back hoppings behave in this way. This same trajectory allowed 21 S1→ S2 back hoppings, with a mean potential energy gap of 0.30 eV and mean available kinetic energy of 3.44 eV. With the reduced reservoir, only four of these hoppings would be allowed. This analysis demonstrates that the amount of kinetic energy available for back hopping with the conventional reservoir would unreasonably allow most back hopping events, while the reduced kinetic energy reservoir delivers a more balanced description of the nonadiabatic events.

FIG. 2.

State occupation of a single surface hopping trajectory at ADC(2)/SV(P) level with (a) conventional velocity adjustment after hopping and (b) reduced kinetic energy adjustment.

FIG. 2.

State occupation of a single surface hopping trajectory at ADC(2)/SV(P) level with (a) conventional velocity adjustment after hopping and (b) reduced kinetic energy adjustment.

Close modal

Vertical excitation energies and oscillator strengths for pyrene at the S0, S1, and S2 state minima are shown in Table I. They were calculated with ADC(2) using the SV(P) basis set. At the ground-state (S0) minimum geometry, the bright states are S2, S4, and S9. The S1 state at its minimum geometry is dark with a near-zero oscillator strength. The S2 state, at its minimum, has a strong oscillator strength of 0.544. The available experimental absorption and emission spectra are redshifted compared to the theory, but we should be cautious with the comparison. All three bright absorption bands have well-resolved vibrational progressions. In each case, we report the band maximum, but there is no guarantee that it would correspond to the vertical excitation. The accuracy of ADC(2)/SV(P) for describing pyrene will be better addressed when discussing the absorption and emission spectrum simulations later in this section. For comparing these results to other electronic structure methods (for the first three excited states), see Ref. 40. Vertical excitation energies and oscillator strengths of the bright states (S2, S4, and S9) computed with several single and multi-configurational methods are discussed in Ref. 41. For recent results for the first seven excited states obtained with restricted active space perturbation theory to the second order (RASPT2) with several active spaces, see Ref. 42.

TABLE I.

Pyrene ADC(2)/SV(P) vertical excitation energies and oscillator strengths at the S0, S1, and S2 minima. Experimental results are given in parentheses.

S0 minimumS1 minimumS2 minimum
StateE (eV)Osc. st.E (eV)Osc. st.E (eV)Osc. st.
S0 0.00 ⋯ 0.11 ⋯ 0.16 ⋯ 
S1 3.78 (3.46a0.001 3.55 (3.37b0.001 3.61 0.001 
S2 4.12 (3.87c0.443 3.94 0.428 3.79 0.544 
S3 4.74 0.000 4.65 0.000 4.42 0.000 
S4 5.00 (4.68c0.450 4.79 0.473 4.77 0.000 
S5 5.10 0.000 4.99 0.000 4.83 0.439 
S6 5.23 0.000 5.07 0.000 4.93 0.000 
S7 5.64 0.000 5.48 0.000 5.49 0.000 
S8 5.75 0.000 5.60 1.250 5.62 0.000 
S9 5.79 (5.39c1.185 5.61 0.000 5.86 1.030 
S10 6.21 0.000 6.11 0.000 5.98 0.000 
S0 minimumS1 minimumS2 minimum
StateE (eV)Osc. st.E (eV)Osc. st.E (eV)Osc. st.
S0 0.00 ⋯ 0.11 ⋯ 0.16 ⋯ 
S1 3.78 (3.46a0.001 3.55 (3.37b0.001 3.61 0.001 
S2 4.12 (3.87c0.443 3.94 0.428 3.79 0.544 
S3 4.74 0.000 4.65 0.000 4.42 0.000 
S4 5.00 (4.68c0.450 4.79 0.473 4.77 0.000 
S5 5.10 0.000 4.99 0.000 4.83 0.439 
S6 5.23 0.000 5.07 0.000 4.93 0.000 
S7 5.64 0.000 5.48 0.000 5.49 0.000 
S8 5.75 0.000 5.60 1.250 5.62 0.000 
S9 5.79 (5.39c1.185 5.61 0.000 5.86 1.030 
S10 6.21 0.000 6.11 0.000 5.98 0.000 
a

Experimental absorption maximum of pyrene in a supersonic jet.49 

b

Experimental S1 band origin of pyrene in a supersonic jet.49 

c

Experimental absorption maximum of pyrene vapor.13 

Describing the electronic structure of the first two singlet excited states of pyrene (and, more generally, of polycyclic aromatic hydrocarbons) is still a tremendous theoretical challenge,41–44 playing a decisive role in the fluorescence, and absorption spectra.45 The dark 1Lb (in Platt’s notation46) and bright 1La states are the two lowest-lying singlet transitions. These states with distinct natures (1Lb is mainly ionic, and 1La is essentially covalent) are well described with ADC(2), which correctly predicts the 1La/1Lb order.41 In contrast, time-dependent density functional theory tends to deliver wrong state ordering41,47,48 unless double-hybrid density functionals are used.43 The results from Ref. 41 also make clear that the basis set size—SV(P) compared with def2-TVZP—has a minor effect on the absorption properties.

ADC(2) may not describe well the ultrafast relaxation of pyrene after high energy excitation due to the low accuracy of the doubly excited states with this method.18 As discussed in Ref. 42, RASPT2 calculations show such a state near the second bright state. Nevertheless, this limitation of ADC(2) does not affect our simulations, which are concerned only with dynamics at later times, after the initial relaxation.

Pyrene’s absorption cross sections are shown in Fig. 3. The theoretical results are from a nuclear ensemble simulation based on 1000 points with excitations computed with ADC(2)/SV(P) (δJ = 0.05 eV for all states). The experimental spectrum is the vapor data at 150 °C from Ref. 13. Both sets of results are given in absolute units. Below 6 eV, pyrene has three prominent absorption bands corresponding to S2 at 4.12 eV, S4 at 5.00 eV, and S9 at 5.79 eV. The corresponding experimental maxima are at 3.87, 4.68, and 5.39 eV, respectively. Experimental and theoretical envelop band intensities nicely agree, but the experimental bands show vibrational progressions that the nuclear ensemble approach based on ensembles built exclusively in the source state cannot simulate. See Ref. 50 for a discussion about how the method can be generalized to describe vibrational progressions.

FIG. 3.

Pyrene gas-phase photoabsorption cross-section as a function of the excitation energy. Theoretical results using the nuclear ensemble approach with ADC(2)/SV(P). The shaded regions indicate the low- and high-energy windows excited in the dynamics simulations. The experimental results are from Ref. 13.

FIG. 3.

Pyrene gas-phase photoabsorption cross-section as a function of the excitation energy. Theoretical results using the nuclear ensemble approach with ADC(2)/SV(P). The shaded regions indicate the low- and high-energy windows excited in the dynamics simulations. The experimental results are from Ref. 13.

Close modal

ADC(2)/SV(P) predicts well the energy of the low-energy band, despite the 0.25-eV shift between the vertical excitation and the band maximum reported in Table I. Nevertheless, the middle- and high-energy bands are blue-shifted by about 0.20 and 0.25 eV, respectively.

Starting in the high-energy excitation window (5.7 ± 0.1 eV), corresponding to excitations into S7 and S8, pyrene quickly populates the manifold of excited states, as we can see in the population evolution shown in Fig. 4 (top). The population of the high-excited state (S3–S8) decays exponentially with an 80-fs time constant. S2 initially receives this population and immediately transfers it to S1. The S1 population grows with a 112-fs time constant. After this transient phase, most of the population is in S1, with a minor fraction in S2. Starting in the low-energy window (4.0 ± 0.1 eV), only S2 is initially excited. It transfers the population to S1 within 47 fs [Fig. 4 (bottom)]. Once more, after the transient phase, a small population remains in S2.

FIG. 4.

Excited-state populations (dotted lines) and occupations (solid lines) of pyrene in the first 500 fs from the nonadiabatic dynamics at the ADC(2)/SV(P) level. Top: dynamics starting from S7 and S8 after 5.7 eV excitation. Bottom: dynamics starting from S2 after 4 eV excitation.

FIG. 4.

Excited-state populations (dotted lines) and occupations (solid lines) of pyrene in the first 500 fs from the nonadiabatic dynamics at the ADC(2)/SV(P) level. Top: dynamics starting from S7 and S8 after 5.7 eV excitation. Bottom: dynamics starting from S2 after 4 eV excitation.

Close modal

Figure 4 shows the excited-state population (mean value of the time-dependent Schrödinger equation coefficients for each state) and occupation (fraction of trajectories in each state). Usually, the two quantities match in decoherence-corrected surface hopping.24 However, the reduced kinetic-energy adjustment (see Sec. II D) causes a slight deviation of about 5% between them.

We considered that trajectories after 500 fs were already in the stationary regime for both excitation windows. Thus, discarding the first 500 fs of each trajectory, we were left with a cumulative time of 21 ps in the high-energy window and 22 ps in the low-energy window (Table II). S1 was occupied in 96.6% of time steps during the stationary regime in the high-energy window and 98.4% in the low one. S2 had the complementary occupation, 3.4% and 1.6% in the high- and low-energy windows, respectively.

TABLE II.

Statistics over the stationary regime (>500 fs) of the trajectories in the high- and low-energy windows. Standard deviations are given in parentheses.

Time stepsOccupation (%)Mean ΔE (eV)Mean f
High-energy excitation (5.7 ± 0.1 eV) 
S1 40 871 96.6 3.31 (0.18) 0.027 (0.073) 
S2 1460 3.4 3.57 (0.20) 0.319 (0.153) 
Cumulative time (fs) 21 165.5    
Low-energy excitation (4.0 ± 0.1 eV) 
S1 45 070 98.4 3.39 (0.15) 0.017 (0.053) 
S2 752 1.6 3.60 (0.15) 0.339 (0.183) 
Cumulative time (fs) 22 911    
Time stepsOccupation (%)Mean ΔE (eV)Mean f
High-energy excitation (5.7 ± 0.1 eV) 
S1 40 871 96.6 3.31 (0.18) 0.027 (0.073) 
S2 1460 3.4 3.57 (0.20) 0.319 (0.153) 
Cumulative time (fs) 21 165.5    
Low-energy excitation (4.0 ± 0.1 eV) 
S1 45 070 98.4 3.39 (0.15) 0.017 (0.053) 
S2 752 1.6 3.60 (0.15) 0.339 (0.183) 
Cumulative time (fs) 22 911    

In the high energy window, the mean value of the ΔE10 energy gap in the stationary regime was 3.31 eV (Table II), which is 0.24 eV redshifted from the vertical excitation at the S1 minimum (Table I). In the low-energy window, this redshift is 0.16 eV. The mean ΔE20 energy gap is redshifted from the vertical value at the S2 minimum by 0.22 eV in the high- and 0.19 eV in the low-energy windows.

The mean S1 → S0 oscillator strength grows 27 times in the high-energy window (17 times in the low one) compared with its value at the S1 minimum. This growth is due to the vibronic coupling with the S2 state when pyrene loses its planarity during dynamics.51 The mean oscillator strength of S2 drops from 0.544 at the S2 minimum to about 0.3 in both excitation windows.

Recently, Aleotti et al.42 reported nonadiabatic wave packet dynamics simulations for pyrene after initial excitation into the first (S2) and second (S5) excited bright states. Their results are based on a 49-dimensional linear vibronic coupling (LVC) Hamiltonian parameterized with RASPT2 data and multi-layer multiconfigurational time-dependent Hartree (ML-MCTDH) propagation. They focus their simulations on the initial relaxation within 2 ps, which does not provide enough information for analyzing fluorescence. Nevertheless, they observed that the internal conversion between the S1 and S2 states is in the 100-fs scale even though pyrene does not ever reach the crossing between these states. Based on this result, they raise the question of “whether semi-classical trajectory-based approaches, which treat nuclei classically, are capable of capturing the ultrafast nature of the internal conversion.” As discussed in Sec. III B (see Fig. 4), the answer is positive.

Figure 5(a) shows the VCA emission spectrum computed with Eq. (13) for pyrene after photoexcitation at 5.7 eV at ADC(2)/SV(P) level. Using Eq. (12), the vibrational temperature is 1001 K. The S2 emission dominates the spectrum, while the S1 emission appears as a shoulder. The Boltzmann populations are 98.9% for S1 and 1.1% for S2. Although the S2 population is tiny, this state’s much bigger oscillator strength (f02/f01 = 544) explains the dominance of S2 emission.

FIG. 5.

Pyrene emission spectra. S1 and S2 emission peaks are shown with dotted and dashed lines, respectively. Computed with the vertical convolution approach [Eq. (13)] using σJ = 0.15 eV for all states at ADC(2)/SV(P) level. (a) Excitation at 5.7 eV and S1 oscillator strength as computed at the S1 minimum. (b) Excitation at 5.7 eV and S1 oscillator strength is assumed to be 0.015 due to vibronic couplings. (c) Excitation at 4 eV and S1 oscillator strength as computed at the S1 minimum. Note the shift in the vertical scale. (d) Excitation at 4 eV and S1 oscillator strength is assumed to be 0.015 due to vibronic couplings. The intensities indicate emission rate, not yield.

FIG. 5.

Pyrene emission spectra. S1 and S2 emission peaks are shown with dotted and dashed lines, respectively. Computed with the vertical convolution approach [Eq. (13)] using σJ = 0.15 eV for all states at ADC(2)/SV(P) level. (a) Excitation at 5.7 eV and S1 oscillator strength as computed at the S1 minimum. (b) Excitation at 5.7 eV and S1 oscillator strength is assumed to be 0.015 due to vibronic couplings. (c) Excitation at 4 eV and S1 oscillator strength as computed at the S1 minimum. Note the shift in the vertical scale. (d) Excitation at 4 eV and S1 oscillator strength is assumed to be 0.015 due to vibronic couplings. The intensities indicate emission rate, not yield.

Close modal

Nevertheless, this S2-emission dominance is not experimentally observed.9,10 The disagreement is due to vibronic couplings induced by the symmetry breaking during the vibrational motion,49,51 causing the oscillator strength of S1 to be much stronger. VCA cannot predict vibronic couplings, but if we assume that the mean oscillator strength of S1 is 0.015, the emission spectrum changes substantially, as shown in Fig. 5 (b). (We have seen in Sec. III B that the mean value of the S1 oscillator strength is 0.027 for the high excitation window. However, we want to keep the discussion about the VCA spectrum entirely independent of the dynamics. For this reason, we will proceed with this hypothetical 0.015 value, corresponding to 15 times the oscillator strength at the minimum of the state.) Considering vibronic coupling, S1 dominates the emission, and the S2 signal becomes a shoulder, as experimentally observed. Note yet that the S2 emission does not change significantly between both simulations for a 5.7-eV excitation. Its maximum intensity remains at about 0.6 × 10−8. However, the S1 intensity increases by a factor of 15.

When pyrene is excited at 4 eV, the vibrational temperature is 582 K, and the Boltzmann populations of S1 and S2 are 99.96% and 0.04%, respectively. The strong S2 → S0 oscillator strength cannot compensate for this near-zero S2 population, and the S1 emission dominates the spectrum [Fig. 5(c)]. Nevertheless, the differential emission rate is a factor of ten smaller than in the high-energy excitation cases. Assuming a factor 15 increase of the S1 → S0 oscillator strength due to vibronic coupling enhances the S1 dominance, and the S2 shoulder disappears [Fig. 5(d)], as also observed experimentally.9,10

These VCA simulations show two main relevant variables affecting pyrene’s fluorescence: the relative population and the vibronic coupling between S1 and S2. The relative population controls the relative amount of emission from S2, whereas the vibronic coupling controls the emission intensity from S1. Nevertheless, the vibronic coupling was included on a crude ad hoc basis. Section IV B discusses this feature on more solid grounds based on vibronic couplings estimates within the nuclear ensemble approach.

All 42 331 time steps accumulated in the stationary regime of the dynamics in the high-energy window (Table II) were used to compute the nuclear ensemble fluorescence spectrum with Eq. (7) using δJ = 0.05 eV for all states. Figure 6(a) shows that the emission band peaks at 3.38 eV and features a prominent shoulder at 3.7 eV. The spectrum decomposition reveals the contributions of S1 and S2 emissions. S1 dominates the band, and the shoulder is due to S2 emission.

FIG. 6.

Fluorescence spectrum of pyrene. Nuclear ensemble (NEA) simulations after an (a) high- and (b) low-energy excitation. The contributions of S1 and S2 emissions are shown as well. Panels (c) and (d) show experimental (Ref. 9) and theoretical [same as in (a) and (b)] fluorescence after a (c) high- and (d) low-energy excitations. The experimental intensity was rescaled to match the theoretical maximum. The theoretical emission energies were redshifted by 0.26 eV in (c) and 0.27 eV in (d) to match the experiments. The intensities indicate emission rate, not yield.

FIG. 6.

Fluorescence spectrum of pyrene. Nuclear ensemble (NEA) simulations after an (a) high- and (b) low-energy excitation. The contributions of S1 and S2 emissions are shown as well. Panels (c) and (d) show experimental (Ref. 9) and theoretical [same as in (a) and (b)] fluorescence after a (c) high- and (d) low-energy excitations. The experimental intensity was rescaled to match the theoretical maximum. The theoretical emission energies were redshifted by 0.26 eV in (c) and 0.27 eV in (d) to match the experiments. The intensities indicate emission rate, not yield.

Close modal

Similarly, the 45 822 time steps accumulated in the stationary regime of low-energy-window dynamics were used to compute the fluorescence spectrum shown in Fig. 6(b). The emission band peaks at 3.50 eV. Once more, S1 dominates the spectrum. Although the shoulder is not visible, S2 still significantly contributes to the emission.

The band shapes in both excitation windows qualitatively agree with the experimental results from Ref. 9, as shown in Figs. 6(c) and 6(d). The theoretical energies are blueshifted by about 0.26 eV due to the ADC(2)/SV(P) level (see discussion in Sec. III A). The experimental shoulder reaches larger emission energies than the theoretical one. The reason is that such large emission energies correspond to the transitions from high vibrational levels in S2 to low vibrational levels in S0. This feature is not modeled by NEA, which only considers the potential energy gap when computing the transition lines.

The fluorescence lifetime [Eq. (2)] is reported in Table III together with the experimental results from Ref. 9. NEA predicts a fluorescence lifetime of 55 ns for the high-energy window, which is in excellent agreement with the experimental result, 62 ns. For the low-energy window, the theoretical prediction, 89 ns, is 2.3 times shorter than the experimental measurement. The most likely reason for this difference is that while the experimental results are for excitation at the S2 band origin, our simulations consider excitations at 4.0 ± 0.1, 0.45 eV above the computed S2 band origin (see Table I). This larger energy excess causes transition energies, oscillator strengths, and S2 populations to be bigger than those at the band origin, yielding shorter fluorescence lifetimes.

TABLE III.

Experimental (Ref. 9) and theoretical (nuclear ensemble) fluorescence lifetimes after high- and low-energy excitations.

Fluorescence lifetime (ns)
Experimental (Ref. 9)Theoretical (NEA)
High energy 62 55 
Low energy 210 89 
Fluorescence lifetime (ns)
Experimental (Ref. 9)Theoretical (NEA)
High energy 62 55 
Low energy 210 89 

Non-Kasha behavior can be achieved either by trapping part of the population in the S2 state minimum (static equilibrium) or continuously repopulating S2 from S1 (dynamic equilibrium).4 Both scenarios are schematically illustrated in Fig. 7. In Ref. 4, Itoh distinguished between three mechanisms yielding non-Kasha behavior. The static-equilibrium scenario corresponds to Itoh’s mechanism C. The dynamic-equilibrium scenario corresponds to Itoh’s mechanism A if the molecule relaxes to the S1 minimum and mechanism B if it does not. These non-Kasha mechanisms have also been characterized through transition rate theory by Veys and Escudero.15 The static-equilibrium scenario is equivalent to their type I, while the dynamic equilibrium is type II. They distinguish yet a type III mechanism where S1 is populated via electron-energy transfer.

FIG. 7.

Schematic S1 and S2 occupations. Each horizontal line represents a single molecule’s state occupation as a function of time (trajectory). (a) Most trajectories populate S1, but a few are trapped in the S2 minimum (static equilibrium). (b) All trajectories populate S1, but they can spend short periods in S2 (dynamic equilibrium). The relative population of both states is the same in both scenarios.

FIG. 7.

Schematic S1 and S2 occupations. Each horizontal line represents a single molecule’s state occupation as a function of time (trajectory). (a) Most trajectories populate S1, but a few are trapped in the S2 minimum (static equilibrium). (b) All trajectories populate S1, but they can spend short periods in S2 (dynamic equilibrium). The relative population of both states is the same in both scenarios.

Close modal

Our VCA and NEA simulations show that pyrene’s non-Kasha behavior emerges from the dynamic equilibrium between S1 and S2. Pyrene does not get trapped in S2. It populates the S1 minimum, but due to the vibrational energy excess, it can briefly access S2. An example of such behavior is illustrated in Fig. 2(b). Statistically, this dynamic equilibrium yields a steady population on S2, which can spontaneously fluoresce from there. The duration of pyrene’s stay in S2 is exponentially distributed with a 30-fs mean value. This time is too short for complete relaxation on the S2 surface, meaning that pyrene’s non-Kasha fluorescence tends to occur from unrelaxed vibronic levels (corresponding to mechanism B in Itoh’s classification4).

The distinction between static and dynamic equilibrium may have practical implications for the design of materials with non-Kasha behavior. Static-equilibrium non-Kasha emission should be favored by large energy gaps between S1 and S2, reducing the internal conversion rate between these states.16 We may recall that del Valle and Catalán2 invoked the small S1/S2 gap in pyrene to dismiss the non-Kasha assignments. Nevertheless, dynamic-equilibrium non-Kasha emission may be enhanced by small S2/S1 gaps, which facilitates populating S2.

Both VCA and NEA agree that pyrene’s complex fluorescence results from the dynamic equilibrium between the S1 and S2 populations and the vibronic coupling between S1 and S2. There are, however, some disagreements on the details of their description. As discussed in Sec. II B, the vertical transition approximation adopted in VCA for a microcanonical [Eq. (11)] system should reduce the quality of its population prediction. After high-energy excitation, VCA predicted only 1.1% of the population on S2, significantly lower than the 3.4% predicted by dynamics. Consequently, the S2 shoulder in the VCA description [Fig. 5(b)] is weak compared to the experiments. After low-energy excitation, VCA predicted that the S2 population should be near-zero (0.04%), while NEA predicted a much larger population of 1.6%. Thus, according to NEA (but not VCA), the S2 emission is still present after low-energy excitation. It is only hidden under the broader S1 band. Because NEA is a better approximation than VCA, its picture of pyrene fluorescence should likely be more accurate.

Alternatively to the NEA sampling based on dynamics, we could consider using a high-temperature Wigner distribution, which would be much less time-consuming. If we did so, the distribution would populate highly excited vibrational states of S1, but it would not add any population to S2. We could try compensating it by a second Wigner sampling around the S2 minimum. Nevertheless, we would need to determine the relative populations between the two states, likely through a Boltzmann distribution. We would have to assume that the non-Kasha emission followed a dynamic equilibrium to adopt such a procedure as we did in the VCA approach. On the other hand, we do not need such assumptions with the dynamics sampling, letting the simulations reveal whether the molecule follows static or dynamic equilibrium. Therefore, NEA with dynamics sampling is a superior approach despite the computational costs.

We simulated the fluorescence spectrum of pyrene in the gas phase employing two different theoretical procedures, the nuclear ensemble (NEA) and vertical convolution (VCA) approaches. Both predicted that the fluorescence after a high-energy excitation should contain an S1-emission band along an S2-emission shoulder. Furthermore, both approaches predicted that this shoulder should not be present in fluorescence after a low-energy excitation. These results agree with the standard assignment of the fluorescence shoulder to non-Kasha S2 emission.8–10 

Our simulations also help elucidate why the spectra reported in Ref. 2 do not show any explicit sign of S2 emission. This emission depends on the population of the S2 state which, in turn, depends on how much vibrational energy is available in the excited molecule. In Ref. 2, the emission spectra of pyrene were measured in 2-methylbutane at 298 K. Under this condition, the photon energy should dissipate into the environment, cooling down the vibrational modes and reducing the S2 population. Indeed, this effect is corroborated by the experimental results of Baba et al.,9 who showed that the band shoulder disappeared when the pressure of cyclohexane added as inert gas is increased.

Both sets of simulations agree that the non-Kasha fluorescence of pyrene results from a dynamic equilibrium, where the molecule in S1 may access S2 briefly due to the vibrational energy excess, creating a steady S2 population that may fluoresce. The alternative, a static equilibrium where a part of the population is trapped in S2 minimum, was not observed. This distinction between dynamic and static equilibrium has significant implications for designing novel non-Kasha chromophores because they impose opposite requirements regarding S1/S2 energy gaps.

Fluorescence simulations of non-Kasha systems are still challenging for computational chemistry. The research protocols developed in this work—including the fluorescence spectrum simulation approaches, the ergodic approximation to build nuclear ensembles from dynamics, and the size-extensivity correction for velocity adjustment in surface hopping—should help close this methodological gap, providing new insights into these crucial photophysical phenomena.

S.M. and M.B. thank the European Research Council (ERC) advanced grant SubNano (Grant Agreement No. 832237) and the Centre de Calcul Intensif d’Aix-Marseille. E.V. thanks the Brazilian agencies CAPES for funding the grant through the Capes/PrInt project (Grant No. 88887.467063/2019-00). E.V. and S.A.d.M. thank the Brazilian agency CNPq (Grant Nos. 423112/2018-0, 308371/2021-6, and 310123/2020-8) for support. I.B. thanks the Brazilian agencies CNPq (Grant Nos. 304148/2018-0 and 409447/2018-8) and Faperj (Grant No. E26/201.197/2021) for support. This work was supported by the Center for Integrated Nanotechnologies (Project No. 2021BC0076), an Office of Science User Facility operated for the U.S. Department of Energy Office of Science by Los Alamos National Laboratory (Contract No. 89233218CNA000001) and Sandia National Laboratories (Contract No. DE-NA-0003525).

The authors have no conflicts to disclose.

Gabriel Braun: Investigation (equal); Visualization (equal). Itamar Borges, Jr.: Funding acquisition (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Adélia J. A. Aquino: Investigation (equal); Writing – review & editing (equal). Hans Lischka: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Felix Plasser: Methodology (supporting); Writing – review & editing (equal). Silmar A. do Monte: Investigation (equal); Writing – review & editing (equal). Elizete Ventura: Investigation (equal); Writing – review & editing (equal). Saikat Mukherjee: Formal analysis (equal); Investigation (equal); Software (lead); Writing – review & editing (equal). Mario Barbatti: Conceptualization (equal); Funding acquisition (equal); Methodology (lead); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are openly available in Zenodo.52 

1.
M.
Kasha
,
Discuss. Faraday Soc.
9
,
14
(
1950
).
2.
J. C.
del Valle
and
J.
Catalán
,
Phys. Chem. Chem. Phys.
21
,
10061
(
2019
).
3.
T.
Itoh
,
T.
Takemura
, and
H.
Baba
,
Chem. Phys. Lett.
40
,
481
(
1976
).
4.
T.
Itoh
,
Chem. Rev.
112
,
4541
(
2012
).
5.
A. P.
Demchenko
,
V. I.
Tomin
, and
P.-T.
Chou
,
Chem. Rev.
117
,
13353
(
2017
).
6.
D.
Malpicci
,
E.
Lucenti
,
C.
Giannini
,
A.
Forni
,
C.
Botta
, and
E.
Cariati
,
Molecules
26
,
6999
(
2021
).
7.
S. K.
Behera
,
S. Y.
Park
, and
J.
Gierschner
,
Angew. Chem., Int. Ed.
60
,
22624
(
2021
).
8.
P. A.
Geldof
,
R. P. H.
Rettschnick
, and
G. J.
Hoytink
,
Chem. Phys. Lett.
4
,
59
(
1969
).
9.
H.
Baba
,
A.
Nakajima
,
M.
Aoi
, and
K.
Chihara
,
J. Chem. Phys.
55
,
2433
(
1971
).
10.
Y.
Numata
,
Y.
Suzuki
, and
I.
Suzuka
,
J. Photochem. Photobiol., A
237
,
49
(
2012
).
11.
A.
Nakajima
and
H.
Baba
,
Bull. Chem. Soc. Jpn.
43
,
967
(
1970
).
12.
T.
Deinum
,
C. J.
Werkhoven
,
J.
Langelaar
,
R. P. H.
Rettschnick
, and
J. D. W.
van Voorst
,
Chem. Phys. Lett.
12
,
189
(
1971
).
13.
A.
Thöny
and
M. J.
Rossi
,
J. Photochem. Photobiol., A
104
,
25
(
1997
).
14.
K.
Veys
and
D.
Escudero
,
J. Phys. Chem. A
124
,
7228
(
2020
).
15.
K.
Veys
and
D.
Escudero
,
Acc. Chem. Res.
55
,
2698
(
2022
).
16.
A.
Prlj
,
T.
Begušić
,
Z. T.
Zhang
,
G. C.
Fish
,
M.
Wehrle
,
T.
Zimmermann
,
S.
Choi
,
J.
Roulet
,
J.-E.
Moser
, and
J.
Vaníček
,
J. Chem. Theory Comput.
16
,
2617
(
2020
).
17.
C.
Hättig
, in
Advances in Quantum Chemistry
, edited by
H. J. Å.
Jensen
(
Academic Press
,
2005
), p.
37
.
18.
A.
Dreuw
and
M.
Wormit
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
5
,
82
(
2015
).
19.
M.
Hillery
,
R. F.
O’Connell
,
M. O.
Scully
, and
E. P.
Wigner
,
Phys. Rep.
106
,
121
(
1984
).
20.
R.
Crespo-Otero
and
M.
Barbatti
,
Theor. Chem. Acc.
131
,
1237
(
2012
).
21.
J. C.
Tully
,
J. Chem. Phys.
93
,
1061
(
1990
).
22.
G.
Granucci
,
M.
Persico
, and
A.
Toniolo
,
J. Chem. Phys.
114
,
10608
(
2001
).
23.
F.
Plasser
,
G.
Granucci
,
J.
Pittner
,
M.
Barbatti
,
M.
Persico
, and
H.
Lischka
,
J. Chem. Phys.
137
,
22A514
(
2012
).
24.
G.
Granucci
and
M.
Persico
,
J. Chem. Phys.
126
,
134114
(
2007
).
25.
S. G.
Balasubramani
,
G. P.
Chen
,
S.
Coriani
,
M.
Diedenhofen
,
M. S.
Frank
,
Y. J.
Franzke
,
F.
Furche
,
R.
Grotjahn
,
M. E.
Harding
,
C.
Hättig
,
A.
Hellweg
,
B.
Helmich-Paris
,
C.
Holzer
,
U.
Huniar
,
M.
Kaupp
,
A. M.
Khah
,
S. K.
Khani
,
T.
Müller
,
F.
Mack
,
B. D.
Nguyen
,
S. M.
Parker
,
E.
Perlt
,
D.
Rappoport
,
K.
Reiter
,
S.
Roy
,
M.
Rückert
,
G.
Schmitz
,
M.
Sierka
,
E.
Tapavicza
,
D. P.
Tew
,
C.
van Wüllen
,
V. K.
Voora
,
F.
Weigend
,
A.
Wodyński
, and
J. M.
Yu
,
J. Chem. Phys.
152
,
184107
(
2020
).
26.
M.
Barbatti
,
M.
Ruckenbauer
,
F.
Plasser
,
J.
Pittner
,
G.
Granucci
,
M.
Persico
, and
H.
Lischka
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
26
(
2014
).
27.
S.
Bai
,
R.
Mansour
,
L.
Stojanović
,
J. M.
Toldo
, and
M.
Barbatti
,
J. Mol. Model.
26
,
107
(
2020
).
28.
M.
Barbatti
,
J. Chem. Phys.
156
,
204304
(
2022
).
29.
F.
Plasser
,
S.
Mai
,
M.
Fumanal
,
E.
Gindensperger
,
C.
Daniel
, and
L.
González
,
J. Chem. Theory Comput.
15
,
5031
(
2019
).
30.
A.
Carof
,
S.
Giannini
, and
J.
Blumberger
,
J. Chem. Phys.
147
,
214113
(
2017
).
31.
A. W.
Jasper
,
S. N.
Stechmann
, and
D. G.
Truhlar
,
J. Chem. Phys.
116
,
5424
(
2002
).
32.
P.
Pechukas
,
Phys. Rev.
181
,
174
(
1969
).
33.
M. F.
Herman
,
J. Chem. Phys.
81
,
754
(
1984
).
34.
J. C.
Tully
,
Int. J. Quantum Chem.
40
,
299
(
1991
).
35.
F.
Webster
,
P. J.
Rossky
, and
R. A.
Friesner
,
Comput. Phys. Commun.
63
,
494
(
1991
).
36.
M.
Barbatti
,
J. Chem. Theory Comput.
17
,
3010
(
2021
).
37.
E.
Fabiano
,
T. W.
Keal
, and
W.
Thiel
,
Chem. Phys.
349
,
334
(
2008
).
38.
G. A.
Worth
,
M. A.
Robb
, and
B.
Lasorne
,
Mol. Phys.
106
,
2077
(
2008
).
39.
M.
Ruckenbauer
,
M.
Barbatti
,
T.
Müller
, and
H.
Lischka
,
J. Phys. Chem. A
117
,
2790
(
2013
).
40.
I. S. K.
Kerkines
,
I. D.
Petsalakis
,
G.
Theodorakopoulos
, and
W.
Klopper
,
J. Chem. Phys.
131
,
224315
(
2009
).
41.
B.
Shi
,
D.
Nachtigallová
,
A. J. A.
Aquino
,
F. B. C.
Machado
, and
H.
Lischka
,
J. Chem. Phys.
150
,
124302
(
2019
).
42.
F.
Aleotti
,
D.
Aranda
,
M.
Yaghoubi Jouybari
,
M.
Garavelli
,
A.
Nenov
, and
F.
Santoro
,
J. Chem. Phys.
154
,
104106
(
2021
).
43.
L.
Goerigk
and
S.
Grimme
,
J. Chem. Theory Comput.
7
,
3272
(
2011
).
44.
A.
Prlj
,
M. E.
Sandoval-Salinas
,
D.
Casanova
,
D.
Jacquemin
, and
C.
Corminboeuf
,
J. Chem. Theory Comput.
12
,
2652
(
2016
).
45.
S.
Shirai
and
S.
Inagaki
,
RSC Adv.
10
,
12988
(
2020
).
46.
J. R.
Platt
,
J. Chem. Phys.
17
,
484
(
1949
).
47.
S.
Grimme
and
M.
Parac
,
ChemPhysChem
4
,
292
(
2003
).
48.
M.
Parac
and
S.
Grimme
,
Chem. Phys.
292
,
11
(
2003
).
49.
N.
Ohta
,
H.
Baba
, and
G.
Marconi
,
Chem. Phys. Lett.
133
,
222
(
1987
).
50.
A. S.
Petit
and
J. E.
Subotnik
,
J. Chem. Phys.
141
,
014107
(
2014
).
51.
A. Y.
Freidzon
,
R. R.
Valiev
, and
A. A.
Berezhnoy
,
RSC Adv.
4
,
42054
(
2014
).
52.
G.
Braun
,
I.
Borges
, Jr.
,
A. J. A.
Aquino
,
H.
Lischka
,
F.
Plasser
,
S. A.
do Monte
,
E.
Ventura
,
S.
Mukherjee
, and
M.
Barbatti
(2022). “
Pyrene fluorescence data obtained by nonadiabatic dynamics calculations
,” Zenodo.