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

Gabriel Braun, Itamar Borges Jr., Adelia Aquino, Hans Lischka, Felix Plasser, Silmar Andrade Do Monte, Elizete Ventura, Saikat Mukherjee, Mario Barbatti

► To cite this version:

Gabriel Braun, Itamar Borges Jr., Adelia Aquino, Hans Lischka, Felix Plasser, et al. Non-Kasha fluorescence of pyrene emerges from a dynamic equilibrium between excited states. The Journal of Chemical Physics, 2022, 10.1063/5.0113908. [https://arxiv.org/abs/2203.01390](https://arxiv.org/abs/2203.01390)

HAL Id: hal-03812836

https://hal.science/hal-03812836

Submitted on 13 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Non-Kasha fluorescence of pyrene emerges from a dynamic
equilibrium between excited states

Gabriel Braun,<sup>a,1</sup> Itamar Borges Jr,<sup>a,2*</sup> Adélia J. A. Aquino,<sup>b,3</sup> Hans Lischka,<sup>c,4*</sup> Felix Plasser,<sup>d,5</sup> Silmar A. do Monte,<sup>e,6</sup> Elizete Ventura,<sup>e,7*</sup> Saikat Mukherjee,<sup>f,8</sup> Mario Barbatti<sup>f,g,9*</sup>

<sup>a</sup> Departamento de Química, Instituto Militar de Engenharia (IME), Rio de Janeiro, RJ, Brazil

<sup>b</sup> Department of Mechanical Engineering, Texas Tech University, Lubbock, TX, 79409, USA

<sup>c</sup> Department of Chemistry and Biochemistry, Texas Tech University Lubbock, TX 79409-1061, USA

<sup>d</sup> Department of Chemistry, Loughborough University, Loughborough, UK

<sup>e</sup> Universidade Federal da Paraíba, 58059-900, João Pessoa-PB, Brazil

<sup>f</sup> Aix Marseille University, CNRS, ICR, Marseille, France

<sup>g</sup> Institut Universitaire de France, 75231 Paris, France

Corresponding authors.

IB: E-mail: itamar@ime.eb.br

HL: E-mail: hans.lischka@univie.ac.at

EV: E-mail: elizete@quimica.ufpb.br

MB: E-mail: mario.barbatti@univ-amu.fr; website: www.barbatti.org

<sup>1</sup> ORCID: 0000-0001-9363-5600

<sup>3</sup> ORCID: 0000-0003-4891-6512

<sup>5</sup> ORCID: 0000-0003-0751-148X5

<sup>7</sup> ORCID: 0000-0002-1015-7824

<sup>9</sup> ORCID: 0000-0001-9336-6607

<sup>2</sup> ORCID: 0000-0002-8492-1223

<sup>4</sup> ORCID: 0000-0002-5656-3975

<sup>6</sup> ORCID: 0000-0002-5878-1984

<sup>8</sup> ORCID: 0000-0002-0025-4735

ABSTRACT

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<sub>2</sub>) has been recently questioned. To elucidate this issue, we simulated the fluorescence of pyrene using two different theoretical approaches based on the vertical convolution and nonadiabatic dynamics with nuclear ensemble approaches. 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<sub>2</sub> emission. We show that the non-Kasha behavior is a dynamic-equilibrium effect, not caused by a metastable S<sub>2</sub> 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<sub>2</sub> emission was not observed in some experiments.
# I. INTRODUCTION

Kasha's empirical rule¹, ² 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³⁻⁵ and is the basis of diverse applications exploring features such as tunable emission colors and improved luminescent quantum efficiency. ⁶, ⁷

Pyrene (see Figure 1) is a prototypical chromophore with non-Kasha emission from a higher state.⁸⁻¹² 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 Figure 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 low-pressure vapor.⁸, ¹²

<figure><img src="image_2.png" /><figcaption>Figure 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.</figcaption></figure>
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 eV 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 Section II. Alternative approaches to simulate non-Kasha systems have been recently discussed in Ref. ¹⁴, ¹⁵ based on rate theory and Ref. ¹⁶ 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.

## II. THEORETICAL METHODS

### II.1. 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.¹⁷, ¹⁸ 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 10 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 Section II.3).

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,²²,²³ with the simplified decay of mixing decoherence correction ($α = 0.1$ au).²⁴ 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 Section II.4 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 Section II.2.

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.

## II.2. 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

$$ \Gamma_{TOT}(E) = \sum_{J} \rho_{1J} \Gamma_{J \to 0}(E), \qquad (1) $$

where $Γ_{J→0}$ is the dimensionless differential emission rate from $S_J$ to $S_0$.²⁰ 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²⁰

$$ \tau_{Fl} = \hbar \left[ \int \Gamma_{TOT}(E) dE \right]^{-1}. \qquad (2) $$

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²⁰

$$ \Gamma_{J \to 0}^{\text{NEA}}(E) = \frac{e^2}{2\pi\hbar mc^3\epsilon_0} [1 - H(E - h\nu_a)] \times \frac{1}{N_p^{(J)}} \sum_l^{N_p^{(J)}} \Delta E_{0J}^2 (\mathbf{R}_l^{(J)}) f_{0J}(\mathbf{R}_l^{(J)}) g_{\text{Lorentz}}(E - \Delta E_{0J}(\mathbf{R}_l^{(J)}), \delta_J), \quad (3) $$

where c is the speed of light, $\epsilon_0$ is the vacuum permittivity, and e and m are the electron charge and mass. $H(E - h\nu_a)$ is the Heaviside function centered at the absorption energy $h\nu_a$. It ensures that emission always occurs with energy gaps smaller than the excitation energy. $\Delta E_{0J}$ is the energy gap, and $f_{0J}$ is the oscillator strength, and both are calculated between $S_J$ and $S_0$. They are computed for each of the $N_p^{(J)}$ geometries in the nuclear ensemble for state $S_J$. The sum runs over all these geometries, too. $g_{\text{Lorentz}}$ is the normalized Lorentzian line shape

$$ g_{\text{Lorentz}}(\epsilon, 2\delta) = \frac{1}{\pi} \frac{\delta}{\epsilon^2 + \delta^2}, \quad (4) $$

with an arbitrary line half-width $\delta_J = \delta/2$.

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

$$ \rho_{1J}^{E} = \frac{N_{p}^{(J)}}{N_{T}}, \quad (5) $$

where

$$ N_T = \sum_J N_p^{(J)}. \quad (6) $$

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

$$ \Gamma_{\text{TOT}}^{\text{NEA}}(E) = \sum_{J} \rho_{1J}^{E} \Gamma_{J \to 0}^{\text{NEA}}(E). \quad (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):
$$ \frac{1}{N_p^{(J)}} \sum_l^{N_p^{(J)}} \Delta E_{0J}^2(\mathbf{R}_l^{(J)}) f_{0J}(\mathbf{R}_l^{(J)}) g_{\text{Lorentz}}(E - \Delta E_{0J}(\mathbf{R}_l^{(J)}), \delta_J) \\ \approx \Delta E_{0J}^2(\mathbf{R}_J) f_{0J}(\mathbf{R}_J) g_{\text{Gaussian}}(E - \Delta E_{0J}(\mathbf{R}_J), \sigma_J), \quad (8) $$

where, in the right-side term, $\Delta E_{0J}$ and $f_{0J}$ are computed at the minimum geometry $\mathbf{R}_J$ of S$_J$. $g_{\text{Gaussian}}$ is the normalized Gaussian function

$$ g_{\text{Gaussian}}(\varepsilon, \sigma) = \frac{1}{(2\pi)^{1/2} \sigma} e^{-\frac{1}{2} \frac{\varepsilon^2}{\sigma^2}}. \quad (9) $$

The differential emission rate simplifies to

$$ \Gamma_{J \to 0}^{\text{VCA}}(E) = \frac{e^2}{(2\pi)^{3/2} \hbar mc^3 \varepsilon_0} \left[1 - H(E - h\nu_a)\right] \frac{\Delta E_{0J}^2(\mathbf{R}_J) f_{0J}(\mathbf{R}_J)}{\sigma_J} e^{-\frac{1}{2} \left(\frac{E - \Delta E_{0J}(\mathbf{R}_J)}{\sigma_J}\right)^2}. \quad (10) $$

The band width $\sigma_J$ is an open parameter, which can be empirically adjusted.

The total emission in Eq. (1) requires the distribution $\rho_{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 Section II.3). 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,

$$ \rho_{1J}^{B} = \frac{\exp\left(-\frac{(\Delta E_{0J}(\mathbf{R}_1) - \Delta E_{01}(\mathbf{R}_1))}{k_B T_{vib}}\right)}{\sum_L \exp\left(-\frac{(\Delta E_{0L}(\mathbf{R}_1) - \Delta E_{01}(\mathbf{R}_1))}{k_B T_{vib}}\right)}, \quad (11) $$

where the Boltzmann factors are computed at the S$_1$ minimum ($\mathbf{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}$

$$ T_{vib} = \left( \ln \left[ \frac{h\nu_a + E_{ZP0} - E_1 + E_{ZP1}}{h\nu_a + E_{ZP0} - E_1 - E_{ZP1}} \right] \right)^{-1} \frac{2E_{ZP1}}{N_{DF}k_B}, \quad (12) $$

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
$$ \Gamma_{TOT}^{VCA}(E) = \sum_{J} \rho_{1J}^{B} \Gamma_{J \to 0}^{VCA}(E). \qquad (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 $\sigma_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 Section II.3.) Nonetheless, NEA spectra give absolute band shapes. Although the line half-width $\delta_J$ is arbitrary, its value is chosen to be much narrower than the band width $\sigma_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.

## II.3. Nuclear ensemble construction

As explained in Section II.2, 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\nu_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 Section III.2.

## II.4. 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,²⁹,³⁰ 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,³²⁻³⁵ 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.³⁶, ³⁷ 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 $\Delta E_{LJ} = E_J - E_L$, conservation of the total energy is ensured by constraining that variation of kinetic energy $\Delta K_{LJ} = K_J - K_L$ equals minus the variation of potential energy:

$$ \Delta K_{LJ} = -\Delta E_{LJ}. \qquad (14) $$
In practice, this constraining is imposed by changing the velocity of the nucleus $\alpha$ as

$$ \mathbf{v}_{\alpha}^{(J)} = \mathbf{v}_{\alpha}^{(L)} + \gamma_{LJ} \frac{\mathbf{u}_{\alpha}}{M_{\alpha}}, \qquad (15) $$

where $\mathbf{u}_{\alpha}$ is a vector indicating the direction along which the change should be applied, and $M_{\alpha}$ is the nuclear mass. $\gamma_{LJ}$ is determined as explained in Ref. 36. In the particular case that $\mathbf{u}_{\alpha}$ is the momentum direction $M_{\alpha} \mathbf{v}_{\alpha}^{(L)}$, Eq. (15) simplifies to

$$ \mathbf{v}_{\alpha}^{(J)} = \kappa_{LJ} \mathbf{v}_{\alpha}^{(L)}, \qquad (16) $$

where $\kappa_{LJ} = 1 - \gamma_{LJ}$ satisfies Eq. (14) if

$$ \kappa_{LJ} = \sqrt{1 - \frac{\Delta E_{LJ}}{K_L}}. \qquad (17) $$

The hopping is allowed if $\kappa_{LJ}$ is real, which means

$$ K_L - \Delta E_{LJ} \geq 0. \qquad (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 $\gamma_{LJ}$ in Eq. (15) is real. This condition implies that it is allowed when the following inequality holds (see Ref. 36 for details)

$$ \frac{1}{2 \sum_{\alpha} \frac{u_{\alpha}^2}{M_{\alpha}}} \left( \sum_{\alpha} \mathbf{v}_{\alpha}^{(L)} \cdot \mathbf{u}_{\alpha} \right)^2 - \Delta E_{LJ} \geq 0, \qquad (19) $$

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

$$ |\Delta K_{LJ}| \leq \frac{1}{2 \sum_{\alpha} \frac{u_{\alpha}^2}{M_{\alpha}}} \left( \sum_{\alpha} \mathbf{v}_{\alpha}^{(L)} \cdot \mathbf{u}_{\alpha} \right)^2. \qquad (20) $$

This amount is controlled by the orientation between $\mathbf{v}^{(L)}$ and $\mathbf{u}$. The adjustment in the nonadiabatic coupling vector direction implies that the kinetic energy reservoir used to ensure energy conservation is
$$|\Delta K_{LJ}| \le \frac{1}{2 \sum_{\alpha} \frac{(h_{LJ,\alpha})^2}{M_{\alpha}}} \left( \sum_{\alpha} \mathbf{v}_{\alpha}^{(L)} \cdot \mathbf{h}_{LJ,\alpha} \right)^2, \qquad (21)$$

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

$$|\Delta K_{LJ}| \le K_L, \qquad (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 $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

$$\frac{K_L}{N_{DF}} - \Delta E_{LJ} \ge 0. \qquad (23)$$

The rationale for this choice is the following. As mentioned, the semiclassical prescription³²⁻³⁵ implies that the correct kinetic energy reservoir should be restricted to atoms contributing to the nonadiabatic coupling vector [Eq. (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 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 [Figure 2 (a)] has many more back hoppings than the adjustment with the reduced kinetic energy [Figure 2 (b)]. For example, the trajectory in Figure 2 (a) attempted a back hopping from S<sub>1</sub> to S<sub>2</sub> 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<sub>1</sub>→S<sub>2</sub> 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.

<figure><img src="image_3.png" /><figcaption>Figure 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.</figcaption></figure>

### III. THE EXCITED STATES OF PYRENE

#### III.1. Electronic structure and absorption spectrum

Vertical excitation energies and oscillator strengths for pyrene at the S<sub>0</sub>, S<sub>1</sub>, and S<sub>2</sub> state minima are shown in Table 1. They were calculated with ADC(2) using the SV(P) basis set. At the ground-state (S<sub>0</sub>) minimum geometry, the bright states are S<sub>2</sub>, S<sub>4</sub>, and S<sub>9</sub>. The S<sub>1</sub> state at its minimum geometry is dark with a near-zero oscillator strength. The S<sub>2</sub> 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. ⁴⁰. Vertical excitation energies and oscillator strengths of the bright states (S₂, S₄, and S₉) computed with several single and multi-configurational methods are discussed in Ref. ⁴¹. 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. ⁴².

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,⁴¹⁻⁴⁴ playing a decisive role in the fluorescence and absorption spectra.⁴⁵ The dark ¹L<sub>b</sub> (in Platt's notation⁴⁶) and bright ¹L<sub>a</sub> states are the two lowest-lying singlet transitions. These states with distinct natures (¹L<sub>b</sub> is mainly ionic, and ¹L<sub>a</sub> is essentially covalent) are well described with ADC(2), which correctly predicts the ¹L<sub>a</sub>/¹L<sub>b</sub> order.⁴¹ In contrast, time-dependent density functional theory tends to deliver wrong state ordering⁴¹, ⁴⁷, ⁴⁸ unless double-hybrid density functionals are used.⁴³ The results from Ref. ⁴¹ also make clear that the basis set size—SV(P) compared to def2-TVZP—has a minor effect on the absorption properties.

<table><caption>Table 1. Pyrene ADC(2)/SV(P) vertical excitation energies and oscillator strengths at the S₀, S₁, and S₂ minima. Experimental results are given in parentheses.</caption><thead><tr><th rowspan="2">State</th><th colspan="2">S<sub>0</sub> minimum</th><th colspan="2">S<sub>1</sub> minimum</th><th colspan="2">S<sub>2</sub> minimum</th></tr><tr><th><i>E</i> (eV)</th><th>Osc. st.</th><th><i>E</i> (eV)</th><th>Osc. st.</th><th><i>E</i> (eV)</th><th>Osc. st.</th></tr></thead><tbody><tr><td>S<sub>0</sub></td><td>0.00</td><td>-</td><td>0.11</td><td>-</td><td>0.16</td><td>-</td></tr><tr><td>S<sub>1</sub></td><td>3.78 (3.46<sup>a</sup>)</td><td>0.001</td><td>3.55 (3.37<sup>c</sup>)</td><td>0.001</td><td>3.61</td><td>0.001</td></tr><tr><td>S<sub>2</sub></td><td>4.12 (3.87<sup>b</sup>)</td><td>0.443</td><td>3.94</td><td>0.428</td><td>3.79</td><td>0.544</td></tr><tr><td>S<sub>3</sub></td><td>4.74</td><td>0.000</td><td>4.65</td><td>0.000</td><td>4.42</td><td>0.000</td></tr><tr><td>S<sub>4</sub></td><td>5.00 (4.68<sup>b</sup>)</td><td>0.450</td><td>4.79</td><td>0.473</td><td>4.77</td><td>0.000</td></tr><tr><td>S<sub>5</sub></td><td>5.10</td><td>0.000</td><td>4.99</td><td>0.000</td><td>4.83</td><td>0.439</td></tr><tr><td>S<sub>6</sub></td><td>5.23</td><td>0.000</td><td>5.07</td><td>0.000</td><td>4.93</td><td>0.000</td></tr><tr><td>S<sub>7</sub></td><td>5.64</td><td>0.000</td><td>5.48</td><td>0.000</td><td>5.49</td><td>0.000</td></tr><tr><td>S<sub>8</sub></td><td>5.75</td><td>0.000</td><td>5.60</td><td>1.250</td><td>5.62</td><td>0.000</td></tr><tr><td>S<sub>9</sub></td><td>5.79 (5.39<sup>b</sup>)</td><td>1.185</td><td>5.61</td><td>0.000</td><td>5.86</td><td>1.030</td></tr><tr><td>S<sub>10</sub></td><td>6.21</td><td>0.000</td><td>6.11</td><td>0.000</td><td>5.98</td><td>0.000</td></tr></tbody></table>

<footnote><sup>a</sup>Experimental absorption maximum of pyrene in a supersonic jet.⁴⁹ <sup>b</sup>Experimental absorption maximum of pyrene vapor.¹³ <sup>c</sup>Experimental S<sub>1</sub> band origin of pyrene in a supersonic jet.⁴⁹</footnote>
<figure><img src="image_4.png" /></figure>

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

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. ⁴², 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 Figure 3. The theoretical results are from a nuclear ensemble simulation based on 1000 points with excitations computed with ADC(2)/SV(P) ($\delta_I = 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. ⁵⁰ 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 1. Nevertheless, the middle- and high-energy bands are blue-shifted by about 0.20 and 0.25 eV, respectively.
### III.2. Nonadiabatic dynamics

Starting in the high-energy excitation window ($5.7 \pm 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 Figure 4-top. The population of the high-excited state ($S_3$ to $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 \pm 0.1$ eV), only $S_2$ is initially excited. It transfers the population to $S_1$ within 47 fs (Figure 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.²⁴ However, the reduced kinetic-energy adjustment (see Section II.4) 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 2). $S_1$ was occupied in 96.6% of time steps during the stationary period 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.

In the high energy window, the mean value of the $\Delta E_{10}$ energy gap in the stationary regime was 3.31 eV (Table 2), which is 0.24 eV redshifted from the vertical excitation at the $S_1$ minimum (Table 1). In the low-energy window, this redshift is 0.16 eV. The mean $\Delta 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.
<figure><img src="image_5.png" /><figcaption>Figure 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 S<sub>7</sub> and S<sub>8</sub> after 5.7 eV excitation. Bottom: dynamics starting from S<sub>2</sub> after 4 eV excitation.</figcaption></figure>

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

Recently, Aleotti *et al.*<sup>42</sup> reported nonadiabatic wave packet dynamics simulations for pyrene after initial excitation into the first (S<sub>2</sub>) and second (S<sub>5</sub>) 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<sub>1</sub> and S<sub>2</sub> states is in the 100-fs scale even though pyrene does not ever reach the crossing between these states, and raises 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 Section III.2 (see Figure 4), the answer is positive.
<table><caption>Table 2. Statistics over the stationary regime (> 500 fs) of the trajectories in the high- and low-energy windows. Standard deviations are given in parentheses.</caption><thead><tr><th></th><th>Time steps</th><th>Occupation (%)</th><th>Mean ΔE (eV)</th><th>Mean f</th></tr></thead><tbody><tr><td colspan="5">High-energy excitation (5.7 ± 0.1 eV)</td></tr><tr><td>S<sub>1</sub></td><td>40871</td><td>96.6</td><td>3.31 (0.18)</td><td>0.027 (0.073)</td></tr><tr><td>S<sub>2</sub></td><td>1460</td><td>3.4</td><td>3.57 (0.20)</td><td>0.319 (0.153)</td></tr><tr><td>Cumulative time (fs)</td><td>21165.5</td><td></td><td></td><td></td></tr><tr><td colspan="5">Low-energy excitation (4.0 ± 0.1 eV)</td></tr><tr><td>S<sub>1</sub></td><td>45070</td><td>98.4</td><td>3.39 (0.15)</td><td>0.017 (0.053)</td></tr><tr><td>S<sub>2</sub></td><td>752</td><td>1.6</td><td>3.60 (0.15)</td><td>0.339 (0.183)</td></tr><tr><td>Cumulative time (fs)</td><td>22911</td><td></td><td></td><td></td></tr></tbody></table>

# IV. THE FLUORESCENCE SPECTRUM OF PYRENE

## IV.1. 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₂ 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_{02}/f_{01} = 544$) explains the dominance of S₂ emission.
<figure><img src="image_6.png" /></figure>

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

Nevertheless, this S<sub>2</sub>-emission dominance is not experimentally observed.<sup>9, 10</sup> The disagreement is due to vibronic couplings induced by the symmetry breaking during the vibrational motion,<sup>49, 51</sup> causing the oscillator strength of S<sub>1</sub> to be much stronger. VCA cannot predict vibronic couplings, but if we assume that the mean oscillator strength of S<sub>1</sub> is 0.015, the emission spectrum changes substantially, as shown in Figure 5 (b). (We have seen in Section III.2 that the mean value of the S<sub>1</sub> 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.) Considering vibronic coupling, S<sub>1</sub> dominates the emission, and the S<sub>2</sub> signal becomes a shoulder, as experimentally observed. Note yet that the S<sub>2</sub> emission does not change significantly between both simulations for a 5.7-eV excitation. Its maximum intensity remains at about 0.6×10<sup>-8</sup>. However, the S<sub>1</sub> intensity increases by a factor of fifteen.

When pyrene is excited at 4 eV, the vibrational temperature is 582 K, and the Boltzmann populations of S<sub>1</sub> and S<sub>2</sub> are 99.96 and 0.04%, respectively. The strong S<sub>2</sub>→S<sub>0</sub> oscillator strength cannot compensate
for this near-zero S₂ population, and the S₁ emission dominates the spectrum [Figure 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 [Figure 5 (d)], as also observed experimentally.⁹,¹⁰

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.2 discusses this feature on more solid grounds based on vibronic couplings estimates within the nuclear ensemble approach.

## IV.2. Pyrene fluorescence spectrum with nuclear ensembles

All 42331 time steps accumulated in the stationary regime of the dynamics in the high-energy window (Table 2) were used to compute the nuclear ensemble fluorescence spectrum with Eq. (7) using $\delta_f = 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 45822 time steps accumulated in the stationary regime of low-energy-window dynamics were used to compute the fluorescence spectrum shown in Figure 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 Figure 6 (c) and (d). The theoretical energies are blueshifted by about 0.26 eV due to the ADC(2)/SV(P) level (see discussion in Section III.1). 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.
<figure><img src="image_7.png" /><figcaption>Figure 6. Fluorescence spectrum of pyrene. Nuclear ensemble (NEA) simulations after an (a) high- and (b) low-energy excitation. The contributions of S₁ and S₂ 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.</figcaption></figure>

<table><caption>Table 3. Experimental (Ref. 9) and theoretical (nuclear ensemble) fluorescence lifetimes after high- and low-energy excitations.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Fluorescence lifetime (ns)</th></tr><tr><th>Experimental (Ref. 9)</th><th>Theoretical (NEA)</th></tr></thead><tbody><tr><td>High energy</td><td>62</td><td>55</td></tr><tr><td>Low energy</td><td>210</td><td>89</td></tr></tbody></table>

The fluorescence lifetime [Eq. (2)] is reported in Table 3 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 eV, 0.45 eV above the computed S₂ band origin (see Table 1). 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.

V. DYNAMIC EQUILIBRIUM IN PYRENE FLUORESCENCE

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

<figure><img src="image_8.png" /><figcaption>Figure 7. Schematic S₁ and S₂ occupations. Each horizontal line represents a single molecule's state occupation as a function of time (trajectory). (a) Most trajectories populate S₁, but a few are trapped in the S₂ minimum (static equilibrium). (b) All trajectories populate S₁, but they can spend short periods in S₂ (dynamic equilibrium). The relative population of both states is the same in both scenarios.</figcaption></figure>

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 Figure 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. In turn, 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 Section II.2, 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 (Figure 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
produce 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.

# 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<sub>1</sub>-emission band along an S<sub>2</sub>-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<sub>2</sub> emission.<sup>8-10</sup>

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<sub>1</sub> may access S<sub>2</sub> briefly due to the vibrational energy excess, creating
a steady S<sub>2</sub> population that may fluoresce. The alternative, a static equilibrium where a part of the
population is trapped in S<sub>2</sub> 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<sub>1</sub>/S<sub>2</sub> 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.

DATA AVAILABILITY

All input and output results from dynamics are available in a public repository.⁵² The dataset also contains
a script to build the ensemble from dynamics to be used in the fluorescence spectrum simulations.

ACKNOWLEDGMENTS

S.M and M.B. thank the European Research Council (ERC) advanced grant SubNano (Grant agreement 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 number 88887.467063/2019-00). E.V. and S. A. M. thank the Brazilian agency CNPq (Grant Numbers 423112/2018-0, 308371/2021-6 and 310123/2020-
8) for support. I. B. thanks the Brazilian agencies CNPq (Grant numbers 304148/2018-0 and 409447/2018-8) and Faperj (Grant number 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 89233218CNA000001) and Sandia National Laboratories (Contract DE-NA-0003525).

# AUTHOR DECLARATIONS

## Conflict of Interest

The authors have no conflicts to disclose.

# AUTHORS CONTRIBUTIONS

Conceptualization: HL, MB; Formal Analysis: SM; Funding acquisition: IB, HL, MB; Investigation: GB, AJAA, SAM, EV, SM; Methodology: FP, MB; Project administration: HL; Software: SM; Supervision: IB, HL, MB; Visualization: GB, MB; Writing – original draft: IB, HL, MB; Writing – review & editing: IB, AJAA, HL, FP, SAM, EV, SM, MB.

# REFERENCES

¹ M. Kasha, Discuss. Faraday Soc. **9** (1950) 14.

https://doi.org/10.1039/DF9500900014

² J. C. del Valle, and J. Catalán, Phys. Chem. Chem. Phys. **21** (2019) 10061.

https://doi.org/10.1039/C9CP00739C

³ T. Itoh, T. Takemura, and H. Baba, Chem. Phys. Lett. **40** (1976) 481.

https://doi.org/10.1016/0009-2614(76)85123-8

⁴ T. Itoh, Chem. Rev. **112** (2012) 4541.

https://doi.org/10.1021/cr200166m

⁵ A. P. Demchenko, V. I. Tomin, and P.-T. Chou, Chem. Rev. **117** (2017) 13353.

https://doi.org/10.1021/acs.chemrev.7b00110

⁶ D. Malpicci, E. Lucenti, C. Giannini, A. Forni, C. Botta, and E. Cariati, Molecules **26** (2021)

https://doi.org/10.3390/molecules26226999

⁷ S. K. Behera, S. Y. Park, and J. Gierschner, Angew. Chem. Int. Ed. **60** (2021) 22624.
https://doi.org/10.1002/anie.202009789

⁸P. A. Geldof, R. P. H. Rettschnick, and G. J. Hoytink, Chem. Phys. Lett. **4** (1969) 59.

https://doi.org/10.1016/0009-2614(69)85066-9

⁹H. Baba, A. Nakajima, M. Aoi, and K. Chihara, J. Chem. Phys. **55** (1971) 2433.

https://doi.org/10.1063/1.1676429

¹⁰Y. Numata, Y. Suzuki, and I. Suzuka, J. Photochem. Photobiol. A **237** (2012) 49.

https://doi.org/10.1016/j.jphotochem.2012.04.004

¹¹A. Nakajima, and H. Baba, Bull. Chem. Soc. Jpn. **43** (1970) 967.

https://doi.org/10.1246/bcsj.43.967

¹²T. Deinum, C. J. Werkhoven, J. Langelaar, R. P. H. Rettschnick, and J. D. W. van Voorst, Chem. Phys. Lett. **12** (1971) 189.

https://doi.org/10.1016/0009-2614(71)80647-4

¹³A. Thöny, and M. J. Rossi, J. Photochem. Photobiol. A **104** (1997) 25.

https://doi.org/10.1016/S1010-6030(96)04575-3

¹⁴K. Veys, and D. Escudero, J. Phys. Chem. A **124** (2020) 7228.

https://doi.org/10.1021/acs.jpca.0c05205

¹⁵K. Veys, and D. Escudero, Acc. Chem. Res. (2022) https://doi.org/10.1021/acs.accounts.2c00453

¹⁶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** (2020) 2617.

https://doi.org/10.1021/acs.jctc.0c00079

¹⁷C. Hättig, in *Adv. Quantum Chem.*, edited by H. J. Å. Jensen (Academic Press, 2005), pp. 37.

https://doi.org/10.1016/S0065-3276(05)50003-0

¹⁸A. Dreuw, and M. Wormit, WIREs: Comp. Mol. Sci. **5** (2015) 82.

https://doi.org/10.1002/wcms.1206

¹⁹M. Hillery, R. F. O'Connell, M. O. Scully, and E. P. Wigner, Phys. Rep. **106** (1984) 121.

https://doi.org/10.1016/0370-1573(84)90160-1

²⁰R. Crespo-Otero, and M. Barbatti, Theor. Chem. Acc. **131** (2012) 1237.

https://doi.org/10.1007/s00214-012-1237-4

²¹J. C. Tully, J. Chem. Phys. **93** (1990) 1061.

https://doi.org/10.1063/1.459170

²²G. Granucci, M. Persico, and A. Toniolo, J. Chem. Phys. **114** (2001) 10608.
https://doi.org/10.1063/1.1376633

23 F. Plasser, G. Granucci, J. Pittner, M. Barbatti, M. Persico, and H. Lischka, J. Chem. Phys. **137** (2012) 22A514.

https://doi.org/10.1063/1.4738960

24 G. Granucci, and M. Persico, J Chem Phys **126** (2007) 134114.

https://doi.org/10.1063/1.2715585

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. v. Wüllen, V. K. Voora, F. Weigend, A. Wodyński, and J. M. Yu, J. Chem. Phys. **152** (2020) 184107.

https://doi.org/10.1063/5.0004635

26 M. Barbatti, M. Ruckenbauer, F. Plasser, J. Pittner, G. Granucci, M. Persico, and H. Lischka, WIREs: Comp. Mol. Sci. **4** (2014) 26.

https://doi.org/10.1002/wcms.1158

27 S. Bai, R. Mansour, L. Stojanović, J. M. Toldo, and M. Barbatti, J. Mol. Model. **26** (2020) 107.

https://doi.org/10.1007/s00894-020-04355-y

28 M. Barbatti, J. Chem. Phys. **156** (2022) 204304.

https://doi.org/10.1063/5.0090205

29 F. Plasser, S. Mai, M. Fumanal, E. Gindensperger, C. Daniel, and L. González, J. Chem. Theory Comput. **15** (2019) 5031.

https://doi.org/10.1021/acs.jctc.9b00525

30 A. Carof, S. Giannini, and J. Blumberger, J. Chem. Phys. **147** (2017) 214113.

https://doi.org/10.1063/1.5003820

31 A. W. Jasper, S. N. Stechmann, and D. G. Truhlar, J. Chem. Phys. **116** (2002) 5424.

https://doi.org/10.1063/1.1453404

32 P. Pechukas, Phys. Rev. **181** (1969) 174.

https://doi.org/10.1103/PhysRev.181.174

33 M. F. Herman, J. Chem. Phys. **81** (1984) 754.

https://doi.org/10.1063/1.447708

34 J. C. Tully, Int. J. Quantum Chem. **40** (1991) 299.
https://doi.org/10.1002/qua.560400830

<sup>35</sup> F. Webster, P. J. Rossky, and R. A. Friesner, Comput. Phys. Commun. <strong>63</strong> (1991) 494.

https://doi.org/10.1016/0010-4655(91)90272-M

<sup>36</sup> M. Barbatti, J. Chem. Theory Comput. <strong>17</strong> (2021) 3010.

https://doi.org/10.1021/acs.jctc.1c00012

<sup>37</sup> E. Fabiano, T. W. Keal, and W. Thiel, Chem. Phys. <strong>349</strong> (2008) 334.

https://doi.org/10.1016/j.chemphys.2008.01.044

<sup>38</sup> G. A. Worth, M. A. Robb, and B. Lasorne, Mol. Phys. <strong>106</strong> (2008) 2077.

https://doi.org/10.1080/00268970802172503

<sup>39</sup> M. Ruckenbauer, M. Barbatti, T. Müller, and H. Lischka, J. Phys. Chem. A <strong>117</strong> (2013) 2790.

https://doi.org/10.1021/jp400401f

<sup>40</sup> I. S. K. Kerkines, I. D. Petsalakis, G. Theodorakopoulos, and W. Klopper, J. Chem. Phys. <strong>131</strong> (2009) 224315.

https://doi.org/10.1063/1.3271347

<sup>41</sup> B. Shi, D. Nachtigallová, A. J. A. Aquino, F. B. C. Machado, and H. Lischka, J. Chem. Phys. <strong>150</strong> (2019) 124302.

https://doi.org/10.1063/1.5086760

<sup>42</sup> F. Aleotti, D. Aranda, M. Yaghoubi Jouybari, M. Garavelli, A. Nenov, and F. Santoro, J. Chem. Phys. <strong>154</strong> (2021) 104106.

https://doi.org/10.1063/5.0044693

<sup>43</sup> L. Goerigk, and S. Grimme, J. Chem. Theory Comput. <strong>7</strong> (2011) 3272.

https://doi.org/10.1021/ct200380v

<sup>44</sup> A. Prlj, M. E. Sandoval-Salinas, D. Casanova, D. Jacquemin, and C. Corminboeuf, J. Chem. Theory Comput. <strong>12</strong> (2016) 2652.

https://doi.org/10.1021/acs.jctc.6b00245

<sup>45</sup> S. Shirai, and S. Inagaki, RSC Advances <strong>10</strong> (2020) 12988.

https://doi.org/10.1039/C9RA10483F

<sup>46</sup> J. R. Platt, J. Chem. Phys. <strong>17</strong> (1949) 484.

https://doi.org/10.1063/1.1747293

<sup>47</sup> S. Grimme, and M. Parac, ChemPhysChem <strong>4</strong> (2003) 292.

https://doi.org/10.1002/cphc.200390047
⁴⁸ M. Parac, and S. Grimme, Chem. Phys. **292** (2003) 11.
https://doi.org/10.1016/s0301-0104(03)00250-7

⁴⁹ N. Ohta, H. Baba, and G. Marconi, Chem. Phys. Lett. **133** (1987) 222.
https://doi.org/10.1016/0009-2614(87)87055-0

⁵⁰ A. S. Petit, and J. E. Subotnik, J. Chem. Phys. **141** (2014) 014107.
https://doi.org/10.1063/1.4884945

⁵¹ A. Y. Freidzon, R. R. Valiev, and A. A. Berezhnoy, RSC Advances **4** (2014) 42054.
https://doi.org/10.1039/C4RA05574H

⁵² G. Braun, I. Borges Jr, A. J. A. Aquino, H. Lischka, F. Plasser, S. Andrade do Monte, E. Ventura, S. Mukherjee, and M. Barbatti, in *Zenodo DOI: 10.5281/zenodo.68076242022*).
https://doi.org/10.5281/zenodo.6807624