24.1 Theoretical Background

24.1.1 Vibronic spectra at zero temperature

The accurate prediction of absorption and emission spectra often requires a quantum mechanical treatment of nuclear vibration[405]. In addition, the geometrical differences between ground and excited state structures, and the differences in vibrational spectra and normal modes lead to mode mixing, making the prediction of the vibrational structure in electronic spectra non-trivial. Under the neglect of anharmonicity, these effects can be described by the Duschinsky rotation[406] \[\begin{equation} \mathbf{Q}_{\text{i}}=\mathbf{D}+\mathbf{J}\mathbf{Q}_{\text{f}}\, , \end{equation}\](24.1) where \(\mathbf{Q}_{\text{i}}\) and \(\mathbf{Q}_{\text{f}}\) denote initial and final state vibrational coordinates, respectively. \(\mathbf{D}\) is the geometric displacement between ground and excited state vibrational structures and \(\mathbf{J}\) is the Duschinsky rotation matrix. In case of absorption, the initial state is the electronic ground state and the final state is the electronically excited state; in case of emission, the opposite order applies.

Within the harmonic oscillator approximation, the absorption and emission spectra are given by: \[\begin{align} \sigma_{abs}(\omega)&= \frac{4\pi^2\omega}{3c}|\mu_{if}|^2 \sum_{{\bf v}_{\text{f}}}|\langle \theta_{0}(\mathbf{Q}_{\text{i}})|\theta_{\bf v_{\text{f}}} (\mathbf{Q}_{\text{f}})\rangle|^2\delta(E_{\bf v_{\text{f}}}-E_0-\omega) \\ \end{align}\](24.2) and \[\begin{align} \sigma_{em}(\omega)&= \frac{4\omega^3}{3c^3}|\mu_{if}|^2 \sum_{{\bf v}_{\text{f}}}|\langle \theta_{0}(\mathbf{Q}_{\text{i}})|\theta_{\bf v_{\text{f}}} (\mathbf{Q}_{\text{f}})\rangle|^2\delta(E_0-E_{\bf v_{\text{f}}}-\omega), \\ \end{align}\](24.3) respectively. Here, \(E_0\) denotes the absolute energy of the initial state where all quantum numbers are zero; \(E_{\bf v_{\text{f}}}\) denotes the final state with quantum numbers \({\bf v_{\text{f}}}\). Both cases assume that absorption and emission occurs from the lowest vibrational level of the initial state (zero temperature approximation). The absorption spectrum (\(\sigma_{abs}\)) is given as the absorption cross section in atomic units (Bohr\(^2\)); the emission spectrum (\(\sigma_{em}\)) is given as the emission rate in inverse atomic time units.

Writing the delta function in Eqs. 24.2 and 24.3 as a Fourier transform and applying Mehler’s formula[407, 408], the infinite sum in equations 24.2 and 24.3 can be eliminated and written in terms of a generating function \(G(t)\) in the time domain, which for absorption reads \[\begin{equation} \sum_{{\bf v}_{\text{f}}}|\langle \theta_{0}(\mathbf{Q}_{\text{i}})|\theta_{\bf v_{\text{f}}} (\mathbf{Q}_{\text{f}})\rangle|^2\delta(E_{\bf v_{\text{f}}}-E_0-\omega) = \int_{-\infty}^{\infty} dt \exp\big[-i t(\Delta E_{if}-\frac{1}{2}\sum_{j}\omega^f_j-\omega)\big] G(t), \end{equation}\](24.4) with \(\Delta E_{if}\) being the adiabatic excitation energy. The generating function is given by \[\begin{align} G(t) & = 2^{\frac{N}{2}} \left(\frac{\det({\bf S}^{-1} \mathbf{\Omega}_{\text{i}}\mathbf{\Omega}_{\text{f}})}{\det({\bf L})\det({\bf M})}\right)^\frac{1}{2} \exp\left(\mathbf{D}^{T}(\mathbf{\Omega}_{\text{f}}{\bf BJM}^{-1}\mathbf{J}^{T}\mathbf{\Omega}_{\text{f}}\mathbf{B}-\mathbf{\Omega}_{\text{f}}{\bf B}){\bf D}\right), \\ \end{align}\](24.5) where \(\mathbf{\Omega}_{\text{i}}\), \(\mathbf{\Omega}_{\text{f}}\), \(\bf S\), \(\bf B\) are diagonal matrices with \((\Omega_{\text{i}})_{kk}=\omega^{\text{i}}_{k}\), \((\Omega_{\text{f}})_{kk}=\omega^{\text{f}}_{k}\), \(S_{kk}=\sinh(i\omega^{\text{f}}_k t)\), and \(B_{kk}=\tanh(i\omega^{\text{f}}_k t/2)\). \(^T\) denotes the transpose of the matrix. \(\omega^{\text{i}}_k\) and \(\omega^{\text{f}}_k\) denote vibrational frequencies of initial and final state, respectively. Matrices \(\bf L\) and \(\bf M\) are obtained as \({\bf M}= {\bf J}^{T}\mathbf{\Omega}_{\text{f}}{\bf BJ}+\mathbf{\Omega}_{\text{i}}\) and \({\bf L}= {\bf J}^{T} \mathbf{\Omega}_{\text{f}}{\bf B}^{-1}{\bf J}+\mathbf{\Omega}_{\text{i}}\). Hence, knowledge of ground and excited state structures and their vibrational spectra allows the construction of all matrices appearing in Eq. 24.5 and \(G(t)\) can be propagated in time. In practice, \(G(t)\) has to be truncated after a maximum time \(t_{max}\). In addition, \(G(t)\) is multiplied by a damping function \(\exp(-t/\tau)\) with lifetime \(\tau\), which leads to a Lorentzian broadening of the spectral lines. More details can be found in references [405, 409].

24.1.2 Single Vibronic Level Spectra

Radless allows the calculation of spectra originating from vibrationally singly excited initial states, i.e. single vibronic level (SVL) emission spectra and Vibrationally Promoted Electronic Resonance (VIPER) spectra[409]. In the frequency domain the SVL emission spectrum for a vibronic initial state \(|\theta^1_k\rangle\), where mode \(k\) is singly-excited, reads \[\begin{align} {\sigma^1}_{em,k}(\omega)&= \frac{4\omega^3}{3 c^3}|\mu_{if}|^2 \sum_{{\bf v}_\text{f}}|\langle \theta^1_k|\theta_{\bf v_\text{f}} \rangle|^2\delta(\Delta E_{if}+E_{{\bf 0}_\text{i}}+\omega^{\text{i}}_{k}-E_{\bf v_\text{f}}-\omega). \\ \end{align}\](24.6) Analogously, the absorption spectrum from a singly excited vibrational state reads \[\begin{align} {\sigma^1}_{abs,k}(\omega)&= \frac{4\pi^2\omega}{3c} |\mu_{if}|^2 \sum_{{\bf v}_\text{f}}|\langle \theta^1_k|\theta_{\bf v_\text{f}} \rangle|^2\delta(\Delta E_{if}-E_{{\bf 0}_\text{i}}-\omega^{\text{i}}_{k}+E_{\bf v_\text{f}}-\omega). \\ \end{align}\](24.7) Using Mehler’s formula and recursive harmonic oscillator recursive relationships, Eq. (24.6) can be formulated in time-domain for emission \[\begin{align} {\sigma^1}_{em,k}(\omega)& = \frac{4\omega^3}{3 c^3} |\mu_{if}|^2 \frac{1}{2\pi} \int_{-\infty}^{\infty} dt \exp\big[it \left(\Delta E_{if}+E_{{\bf 0}_\text{i}}+\omega^{\text{i}}_{k}-\omega\right)\big] G^1_k(t), \end{align}\](24.8) with the generating function \[\begin{align} G^1_k(t) & = G(t) \times \omega^i_k \left[2 ({ M}^{-1})_{kk}+4{\bf D}^{\dagger}\mathbf{\Omega}_{\text{f}}{\bf BJR}^k{\bf J}^{\dagger}\mathbf{\Omega}_{\text{f}}{\bf BD}-2 ({ L}^{-1})_{kk}\right], \\ \end{align}\](24.9) with \(R^k_{ij}=(M^{-1})_{ik}(M^{-1})_{kj}\). For absorption, the time-domain expression reads \[\begin{align} {\sigma^1}_{abs,k}(\omega)& = \frac{4\pi^2\omega}{3c} |\mu_{if}|^2 \frac{1}{2\pi} \int_{-\infty}^{\infty} dt \exp\left(-it \left(\Delta E_{if}-E_{{\bf 0}_\text{i}}-\omega^{\text{i}}_{k}-\omega\right)\right) G^1_k(t). \end{align}\](24.10) Here, \(E_{{\bf 0}\text{i}}\) denotes the zero point vibrational energy of the initial state. The generating function \(G^1_k(t)\) in Eq. (24.10) is identical in structure to Eq. (24.9), but initial state index \(i\) refers to the electronic ground state and state index \(f\) refers to the electronically excited state.

24.1.3 Thawed Gaussian Approximation

The thawed Gaussian approximation (TGA) [410] accounts for anharmonicity of the potential energy surface by propagating a Gaussian wavepacket along a classical trajectory. The method implemented in radless is a variant called single-Hessian TGA, which requires a single MD trajectory in the final electronic state (ground state for emission, excited state for absorption), along with the information needed for computing vibronic spectra within the harmonic approximation. Below we give a brief outline of the method; the reader is referred to Refs. [411, 412] for more details.

Within this wavepacket-based method, the generating function is computed as \[\begin{equation} G(t) = |\mu_{if}|^2 \langle \psi_0 | \psi_t \rangle, \end{equation}\](24.11) where \[\begin{equation} \psi_t(q) = \frac{1}{(\pi \hbar)^{D/4} (\det Q_t)^{1/2}} \exp \left\{ \frac{i}{\hbar} \left[ \frac{1}{2} (q - q_t)^T \cdot P_t \cdot Q_t^{-1} \cdot (q-q_t) + p_t^T \cdot (q - q_t) + S_t \right] \right\} \end{equation}\](24.12) is the vibrational Gaussian wavepacket defined by position and momentum vectors \(q_t\) and \(p_t\), classical action \(S_t\), and matrices \(P_t\) and \(Q_t\). Initially, the parameters of the wavepacket are determined by the vibrational structure of the initial electronic state, namely, \[\begin{align} &q_0 = 0, p_0 = 0 \\ &Q_0 = (\Omega_i/\hbar)^{-1/2}, P_0 = i (\Omega_i / \hbar)^{1/2}, \end{align}\](24.13–24.14) and \(S_0 = 0\), while the parameters evolve according to \[\begin{align} \dot{q}_t = p_t, &\qquad \dot{p}_t = - V_f^{\prime}(q_t) \\ \dot{Q}_t = P_t, &\qquad \dot{P}_t = - K \cdot Q_t, \\ \dot{S}_t = &\frac{1}{2}p_t^2 - V(q_t), \end{align}\](24.15–24.17) Mass-scaled normal-mode coordinates corresponding to the initial electronic state are used throughout. Here, \(V(q_t)\) and \(V^{\prime}(q_t)\) are the energy and gradient in the final electronic state evaluated at position \(q_t\), while \(K\) is the force constant evaluated at the equilibrium geometry of the final electronic state. Equations (24.15) are simply classical dynamics in the final electronic state, starting from the equilibrium geometry of the initial electronic state. The generating function is computed as prescribed in Sec. 2.4 of Ref. [412].