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 (S_{2}) 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 S_{2} emission. We show that the non-Kasha behavior is a dynamic-equilibrium effect not caused by a metastable S_{2} 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 S_{2} emission was not observed in some experiments.

## I. INTRODUCTION

Kasha’s empirical rule^{1,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, S_{1}, before they fluoresce. Nevertheless, fluorescence from higher-lying states, known as non-Kasha (or anti-Kasha) behavior, is well documented for several systems^{3–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 S_{2} state (3.86 eV; see Fig. 1), pyrene fluorescence has a single band at 3.26 eV (380 nm) assigned to S_{1} 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 S_{2} emission. Similar features and assignments were also reported for pyrene in a supersonic jet^{10} and other low-pressure vapor experiments.^{8,12}

This assignment of a non-Kasha emission, however, is not beyond controversy. del Valle and Catalán^{2} argued that the energy gap between S_{2} and S_{1} is too small to trap pyrene long enough to emit from S_{2}. 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 S_{1} emission.

Therefore, does pyrene emit from S_{2}? Furthermore, if it does, why does the S_{2} 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 S_{2}. 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 S_{1} continuously pumping a small amount of population into S_{2}. If such conditions are not satisfied, the S_{1} emission dominates the fluorescence. Thus, the non-Kasha behavior of pyrene results from a dynamic equilibrium between excited-state populations.

## II. THEORETICAL METHODS

### A. Computational details

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 S_{1} and S_{2} 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 distribution^{19} 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 S_{7} and 20 in the S_{8} state, while in the low-energy window, we started 20 trajectories in the S_{2} 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 hopping^{21} 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.

### B. Fluorescence spectrum simulations

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

where Γ_{J→0} is the dimensionless differential emission rate from S_{J} to S_{0}.^{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 through^{20}

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 is^{20}

where *c* is the speed of light, *ε*_{0} is the vacuum permittivity, and *e* and *m* are the electron charge and mass. $H(E\u2212hva)$ is the Heaviside function centered at the absorption energy *hν*_{a}. It ensures that emission always occurs with energy gaps smaller than the excitation energy. Δ*E*_{0J} is the energy gap, and *f*_{0J} is the oscillator strength; both are calculated between S_{J} and S_{0}. They are computed for each of the $NpJ$ geometries in the nuclear ensemble for state S_{J}. The sum runs over all these geometries, too. *g*_{Lorentz} is the normalized Lorentzian line shape,

with an arbitrary line half-width *δ*_{J} = *δ*/2.

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

where

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):

where, in the right-side term, Δ*E*_{0J} and *f*_{0J} are computed at the minimum geometry **R**_{J} of S_{J}. *g*_{Gaussian} is the normalized Gaussian function,

The differential emission rate simplifies to

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,

where the Boltzmann factors are computed at the S_{1} minimum (**R**_{1}). This approximation assumes that the vibrational degrees of the molecule in S_{1} compose a thermal bath at temperature *T*_{vib}, which can populate higher excited states. In the collisionless regime, the vibrational temperature of the S_{1} state can be written as^{28}

where *E*_{1} is the S_{1} state energy at the S_{1} minimum, *E*_{ZPJ} is the zero-point energy of state S_{J}, and *N*_{DF} is the number of degrees of freedom of the molecule.

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

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 S_{1} (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.

### C. Nuclear ensemble construction

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 *hν*_{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.

### D. Velocity adjustment after hopping

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 Δ*E*_{LJ} = *E*_{J} − *E*_{L}, conservation of the total energy is ensured by constraining that variation of kinetic energy Δ*Κ*_{LJ} = *K*_{J} − *K*_{L} equals minus the variation of potential energy,

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

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\alpha v\alpha L$, Eq. (15) simplifies to

where $\kappa LJ=1\u2212\gamma LJ$ satisfies Eq. (14) if

The hopping is allowed if *κ*_{LJ} is real, which means

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):

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

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

which depends on the projection of **v**^{(L)} on **h**_{LJ}. Since **h**_{LJ} 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

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 S_{8}), 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 (*N*_{DF}) when evaluating whether a back hopping can occur. This means replacing Inequality (18) by

The rationale for this choice is the following. As mentioned, the semiclassical prescription^{32–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 *N*_{DF} degrees,^{38} dividing the available kinetic by *N*_{DF} should (statistically speaking) deliver a sensible reservoir, as it also refers only to one degree of freedom. In the case of pyrene with *N*_{DF} = 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 S_{1} to S_{2} 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 S_{1}→ S_{2} 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.

## III. THE EXCITED STATES OF PYRENE

### A. Electronic structure and absorption spectrum

Vertical excitation energies and oscillator strengths for pyrene at the S_{0}, S_{1}, and S_{2} state minima are shown in Table I. They were calculated with ADC(2) using the SV(P) basis set. At the ground-state (S_{0}) minimum geometry, the bright states are S_{2}, S_{4}, and S_{9}. The S_{1} state at its minimum geometry is dark with a near-zero oscillator strength. The S_{2} 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 (S_{2}, S_{4,} and S_{9}) 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.

. | S_{0} minimum
. | S_{1} minimum
. | S_{2} minimum
. | |||
---|---|---|---|---|---|---|

State . | E (eV)
. | Osc. st. . | E (eV)
. | Osc. st. . | E (eV)
. | Osc. st. . |

S_{0} | 0.00 | ⋯ | 0.11 | ⋯ | 0.16 | ⋯ |

S_{1} | 3.78 (3.46^{a}) | 0.001 | 3.55 (3.37^{b}) | 0.001 | 3.61 | 0.001 |

S_{2} | 4.12 (3.87^{c}) | 0.443 | 3.94 | 0.428 | 3.79 | 0.544 |

S_{3} | 4.74 | 0.000 | 4.65 | 0.000 | 4.42 | 0.000 |

S_{4} | 5.00 (4.68^{c}) | 0.450 | 4.79 | 0.473 | 4.77 | 0.000 |

S_{5} | 5.10 | 0.000 | 4.99 | 0.000 | 4.83 | 0.439 |

S_{6} | 5.23 | 0.000 | 5.07 | 0.000 | 4.93 | 0.000 |

S_{7} | 5.64 | 0.000 | 5.48 | 0.000 | 5.49 | 0.000 |

S_{8} | 5.75 | 0.000 | 5.60 | 1.250 | 5.62 | 0.000 |

S_{9} | 5.79 (5.39^{c}) | 1.185 | 5.61 | 0.000 | 5.86 | 1.030 |

S_{10} | 6.21 | 0.000 | 6.11 | 0.000 | 5.98 | 0.000 |

. | S_{0} minimum
. | S_{1} minimum
. | S_{2} minimum
. | |||
---|---|---|---|---|---|---|

State . | E (eV)
. | Osc. st. . | E (eV)
. | Osc. st. . | E (eV)
. | Osc. st. . |

S_{0} | 0.00 | ⋯ | 0.11 | ⋯ | 0.16 | ⋯ |

S_{1} | 3.78 (3.46^{a}) | 0.001 | 3.55 (3.37^{b}) | 0.001 | 3.61 | 0.001 |

S_{2} | 4.12 (3.87^{c}) | 0.443 | 3.94 | 0.428 | 3.79 | 0.544 |

S_{3} | 4.74 | 0.000 | 4.65 | 0.000 | 4.42 | 0.000 |

S_{4} | 5.00 (4.68^{c}) | 0.450 | 4.79 | 0.473 | 4.77 | 0.000 |

S_{5} | 5.10 | 0.000 | 4.99 | 0.000 | 4.83 | 0.439 |

S_{6} | 5.23 | 0.000 | 5.07 | 0.000 | 4.93 | 0.000 |

S_{7} | 5.64 | 0.000 | 5.48 | 0.000 | 5.49 | 0.000 |

S_{8} | 5.75 | 0.000 | 5.60 | 1.250 | 5.62 | 0.000 |

S_{9} | 5.79 (5.39^{c}) | 1.185 | 5.61 | 0.000 | 5.86 | 1.030 |

S_{10} | 6.21 | 0.000 | 6.11 | 0.000 | 5.98 | 0.000 |

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 ^{1}L_{b} (in Platt’s notation^{46}) and bright ^{1}L_{a} states are the two lowest-lying singlet transitions. These states with distinct natures (^{1}L_{b} is mainly ionic, and ^{1}L_{a} is essentially covalent) are well described with ADC(2), which correctly predicts the ^{1}L_{a}/^{1}L_{b} order.^{41} In contrast, time-dependent density functional theory tends to deliver wrong state ordering^{41,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 S_{2} at 4.12 eV, S_{4} at 5.00 eV, and S_{9} 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.

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.

### B. Nonadiabatic dynamics

Starting in the high-energy excitation window (5.7 ± 0.1 eV), corresponding to excitations into S_{7} and S_{8}, 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 (S_{3}–S_{8}) decays exponentially with an 80-fs time constant. S_{2} initially receives this population and immediately transfers it to S_{1}. The S_{1} population grows with a 112-fs time constant. After this transient phase, most of the population is in S_{1,} with a minor fraction in S_{2}. Starting in the low-energy window (4.0 ± 0.1 eV), only S_{2} is initially excited. It transfers the population to S_{1} within 47 fs [Fig. 4 (bottom)]. Once more, after the transient phase, a small population remains in S_{2}.

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). S_{1} was occupied in 96.6% of time steps during the stationary regime in the high-energy window and 98.4% in the low one. S_{2} had the complementary occupation, 3.4% and 1.6% in the high- and low-energy windows, respectively.

. | Time steps . | Occupation (%) . | Mean ΔE (eV)
. | Mean f
. |
---|---|---|---|---|

High-energy excitation (5.7 ± 0.1 eV) | ||||

S_{1} | 40 871 | 96.6 | 3.31 (0.18) | 0.027 (0.073) |

S_{2} | 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) | ||||

S_{1} | 45 070 | 98.4 | 3.39 (0.15) | 0.017 (0.053) |

S_{2} | 752 | 1.6 | 3.60 (0.15) | 0.339 (0.183) |

Cumulative time (fs) | 22 911 |

. | Time steps . | Occupation (%) . | Mean ΔE (eV)
. | Mean f
. |
---|---|---|---|---|

High-energy excitation (5.7 ± 0.1 eV) | ||||

S_{1} | 40 871 | 96.6 | 3.31 (0.18) | 0.027 (0.073) |

S_{2} | 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) | ||||

S_{1} | 45 070 | 98.4 | 3.39 (0.15) | 0.017 (0.053) |

S_{2} | 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 Δ*E*_{10} energy gap in the stationary regime was 3.31 eV (Table II), which is 0.24 eV redshifted from the vertical excitation at the S_{1} minimum (Table I). In the low-energy window, this redshift is 0.16 eV. The mean Δ*E*_{20} energy gap is redshifted from the vertical value at the S_{2} minimum by 0.22 eV in the high- and 0.19 eV in the low-energy windows.

The mean S_{1} → S_{0} oscillator strength grows 27 times in the high-energy window (17 times in the low one) compared with its value at the S_{1} minimum. This growth is due to the vibronic coupling with the S_{2} state when pyrene loses its planarity during dynamics.^{51} The mean oscillator strength of S_{2} drops from 0.544 at the S_{2} 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 (S_{2}) and second (S_{5}) 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 S_{1} and S_{2} 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.

## IV. THE FLUORESCENCE SPECTRUM OF PYRENE

### A. Pyrene fluorescence with vertical convolutions

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 S_{2} emission dominates the spectrum, while the S_{1} emission appears as a shoulder. The Boltzmann populations are 98.9% for S_{1} and 1.1% for S_{2}. Although the S_{2} population is tiny, this state’s much bigger oscillator strength (*f*_{02}/*f*_{01} = 544) explains the dominance of S_{2} emission.

Nevertheless, this S_{2}-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 S_{1} to be much stronger. VCA cannot predict vibronic couplings, but if we assume that the mean oscillator strength of S_{1} 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 S_{1} 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, S_{1} dominates the emission, and the S_{2} signal becomes a shoulder, as experimentally observed. Note yet that the S_{2} 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 S_{1} 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 S_{1} and S_{2} are 99.96% and 0.04%, respectively. The strong S_{2} → S_{0} oscillator strength cannot compensate for this near-zero S_{2} population, and the S_{1} 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 S_{1} → S_{0} oscillator strength due to vibronic coupling enhances the S_{1} dominance, and the S_{2} 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 S_{1} and S_{2}. The relative population controls the relative amount of emission from S_{2}, whereas the vibronic coupling controls the emission intensity from S_{1}. 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.

### B. Pyrene fluorescence spectrum with nuclear ensembles

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 S_{1} and S_{2} emissions. S_{1} dominates the band, and the shoulder is due to S_{2} emission.

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, S_{1} dominates the spectrum. Although the shoulder is not visible, S_{2} 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 S_{2} to low vibrational levels in S_{0}. 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 S_{2} band origin, our simulations consider excitations at 4.0 ± 0.1, 0.45 eV above the computed S_{2} band origin (see Table I). This larger energy excess causes transition energies, oscillator strengths, and S_{2} populations to be bigger than those at the band origin, yielding shorter fluorescence lifetimes.

## V. DYNAMIC EQUILIBRIUM IN PYRENE FLUORESCENCE

Non-Kasha behavior can be achieved either by trapping part of the population in the S_{2} state minimum (*static equilibrium*) or continuously repopulating S_{2} from S_{1} (*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 S_{1} 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 S_{1} is populated via electron-energy transfer.

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

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 S_{1} and S_{2}, reducing the internal conversion rate between these states.^{16} We may recall that del Valle and Catalán^{2} invoked the small S_{1}/S_{2} gap in pyrene to dismiss the non-Kasha assignments. Nevertheless, dynamic-equilibrium non-Kasha emission may be enhanced by small S_{2}/S_{1} gaps, which facilitates populating S_{2}.

Both VCA and NEA agree that pyrene’s complex fluorescence results from the dynamic equilibrium between the S_{1} and S_{2} populations and the vibronic coupling between S_{1} and S_{2}. 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 S_{2}, significantly lower than the 3.4% predicted by dynamics. Consequently, the S_{2} shoulder in the VCA description [Fig. 5(b)] is weak compared to the experiments. After low-energy excitation, VCA predicted that the S_{2} 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 S_{2} emission is still present after low-energy excitation. It is only hidden under the broader S_{1} 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 S_{1}, but it would not add any population to S_{2}. We could try compensating it by a second Wigner sampling around the S_{2} 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.

## VI. CONCLUSIONS

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 S_{1}-emission band along an S_{2}-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 S_{2} emission.^{8–10}

Our simulations also help elucidate why the spectra reported in Ref. 2 do not show any explicit sign of S_{2} emission. This emission depends on the population of the S_{2} 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 S_{2} 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 S_{1} may access S_{2} briefly due to the vibrational energy excess, creating a steady S_{2} population that may fluoresce. The alternative, a static equilibrium where a part of the population is trapped in S_{2} 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 S_{1}/S_{2} 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.

## ACKNOWLEDGMENTS

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).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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