Non-Kasha fluorescence of pyrene emerges from a dynamic equilibrium between excited states

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₁, 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.⁹ After a low-energy excitation into the S₂ state (3.86 eV; see Fig. 1), pyrene fluorescence has a single band at 3.26 eV (380 nm) assigned to S₁ 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₂ emission. Similar features and assignments were also reported for pyrene in a supersonic jet¹⁰ 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² argued that the energy gap between S₂ and S₁ is too small to trap pyrene long enough to emit from S₂. 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₁ emission.

Therefore, does pyrene emit from S₂? Furthermore, if it does, why does the S₂ 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,⁹ 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₂. 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₁ continuously pumping a small amount of population into S₂. If such conditions are not satisfied, the S₁ 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 S₁ and S₂ 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¹⁹ 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)²⁰ 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₇ and 20 in the S₈ state, while in the low-energy window, we started 20 trajectories in the S₂ 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²¹ via the local diabatization algorithm,^22,23 with the simplified decay of mixing decoherence correction (α = 0.1 a.u.).²⁴ 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.²⁵ Spectrum simulations, initial condition generation, and dynamics were carried out with Newton-X classical series (CS)²⁶ interfaced to TURBOMOLE.

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

where Γ_J→0 is the dimensionless differential emission rate from S_J to S₀.²⁰ 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²⁰ 

This section describes two approaches for computing these quantities. The first is the nuclear ensemble approach (NEA),²⁰ 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)²⁷—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²⁰ 

where c is the speed of light, ε₀ is the vacuum permittivity, and e and m are the electron charge and mass. $H(E−hva)$ 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₀. 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

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

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₁ minimum (R₁). This approximation assumes that the vibrational degrees of the molecule in S₁ compose a thermal bath at temperature T_vib, which can populate higher excited states. In the collisionless regime, the vibrational temperature of the S₁ state can be written as²⁸ 

where E₁ is the S₁ state energy at the S₁ 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₁ (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.¹⁰ 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,⁹ 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.³¹ 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αvαL$⁠, Eq. (15) simplifies to

where $κLJ=1−γ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₈), 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,³⁸ 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.³⁹ 

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₁ to S₂ 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₁→ S₂ 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.

Vertical excitation energies and oscillator strengths for pyrene at the S₀, S₁, and S₂ state minima are shown in Table I. They were calculated with ADC(2) using the SV(P) basis set. At the ground-state (S₀) minimum geometry, the bright states are S₂, S₄, and S₉. The S₁ state at its minimum geometry is dark with a near-zero oscillator strength. The S₂ 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₂, S_4, and S₉) 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.

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.⁴⁵ The dark ¹L_b (in Platt’s notation⁴⁶) and bright ¹L_a states are the two lowest-lying singlet transitions. These states with distinct natures (¹L_b is mainly ionic, and ¹L_a is essentially covalent) are well described with ADC(2), which correctly predicts the ¹L_a/¹L_b order.⁴¹ In contrast, time-dependent density functional theory tends to deliver wrong state ordering^41,47,48 unless double-hybrid density functionals are used.⁴³ 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.¹⁸ 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₂ at 4.12 eV, S₄ at 5.00 eV, and S₉ 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.

Starting in the high-energy excitation window (5.7 ± 0.1 eV), corresponding to excitations into S₇ and S₈, 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₃–S₈) decays exponentially with an 80-fs time constant. S₂ initially receives this population and immediately transfers it to S₁. The S₁ 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₂. Starting in the low-energy window (4.0 ± 0.1 eV), only S₂ is initially excited. It transfers the population to S₁ within 47 fs [Fig. 4 (bottom)]. Once more, after the transient phase, a small population remains in S₂.

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.²⁴ 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₁ 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₂ had the complementary occupation, 3.4% and 1.6% in the high- and low-energy windows, respectively.

In the high energy window, the mean value of the ΔE₁₀ 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₁ minimum (Table I). In the low-energy window, this redshift is 0.16 eV. The mean ΔE₂₀ energy gap is redshifted from the vertical value at the S₂ minimum by 0.22 eV in the high- and 0.19 eV in the low-energy windows.

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

Recently, Aleotti et al.⁴² reported nonadiabatic wave packet dynamics simulations for pyrene after initial excitation into the first (S₂) and second (S₅) 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₁ and S₂ 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 S₂ emission dominates the spectrum, while the S₁ emission appears as a shoulder. The Boltzmann populations are 98.9% for S₁ and 1.1% for S₂. Although the S₂ population is tiny, this state’s much bigger oscillator strength (f₀₂/f₀₁ = 544) explains the dominance of S₂ emission.

Nevertheless, this S₂-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₁ to be much stronger. VCA cannot predict vibronic couplings, but if we assume that the mean oscillator strength of S₁ 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₁ 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₁ dominates the emission, and the S₂ signal becomes a shoulder, as experimentally observed. Note yet that the S₂ emission does not change significantly between both simulations for a 5.7-eV excitation. Its maximum intensity remains at about 0.6 × 10⁻⁸. However, the S₁ 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₁ and S₂ are 99.96% and 0.04%, respectively. The strong S₂ → S₀ oscillator strength cannot compensate for this near-zero S₂ population, and the S₁ 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₁ → S₀ oscillator strength due to vibronic coupling enhances the S₁ dominance, and the S₂ 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₁ and S₂. The relative population controls the relative amount of emission from S₂, whereas the vibronic coupling controls the emission intensity from S₁. 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 S₁ and S₂ emissions. S₁ dominates the band, and the shoulder is due to S₂ 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₁ dominates the spectrum. Although the shoulder is not visible, S₂ 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₂ to low vibrational levels in S₀. 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₂ band origin, our simulations consider excitations at 4.0 ± 0.1, 0.45 eV above the computed S₂ band origin (see Table I). This larger energy excess causes transition energies, oscillator strengths, and S₂ populations to be bigger than those at the band origin, yielding shorter fluorescence lifetimes.

Non-Kasha behavior can be achieved either by trapping part of the population in the S₂ state minimum (static equilibrium) or continuously repopulating S₂ from S₁ (dynamic equilibrium).⁴ 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₁ minimum and mechanism B if it does not. These non-Kasha mechanisms have also been characterized through transition rate theory by Veys and Escudero.¹⁵ 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₁ 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₁ and S₂. Pyrene does not get trapped in S₂. It populates the S₁ minimum, but due to the vibrational energy excess, it can briefly access S₂. An example of such behavior is illustrated in Fig. 2(b). Statistically, this dynamic equilibrium yields a steady population on S₂, which can spontaneously fluoresce from there. The duration of pyrene’s stay in S₂ is exponentially distributed with a 30-fs mean value. This time is too short for complete relaxation on the S₂ surface, meaning that pyrene’s non-Kasha fluorescence tends to occur from unrelaxed vibronic levels (corresponding to mechanism B in Itoh’s classification⁴).

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₁ and S₂, reducing the internal conversion rate between these states.¹⁶ We may recall that del Valle and Catalán² invoked the small S₁/S₂ gap in pyrene to dismiss the non-Kasha assignments. Nevertheless, dynamic-equilibrium non-Kasha emission may be enhanced by small S₂/S₁ gaps, which facilitates populating S₂.

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

Our simulations also help elucidate why the spectra reported in Ref. 2 do not show any explicit sign of S₂ emission. This emission depends on the population of the S₂ 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₂ population. Indeed, this effect is corroborated by the experimental results of Baba et al.,⁹ 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₁ may access S₂ briefly due to the vibrational energy excess, creating a steady S₂ population that may fluoresce. The alternative, a static equilibrium where a part of the population is trapped in S₂ 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₁/S₂ 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.⁵²