Herzberg-Teller Approximation in Vibronic Spectroscopy

• Herzberg-Teller approximation is a first-order expansion of the electronic transition dipole moment with respect to nuclear coordinates, enhancing the analysis of non-Condon effects.
• It systematically enables intensity borrowing and activates vibronic sidebands by incorporating linear dipole gradients that relax Franck–Condon selection rules.
• The method is pivotal in time-dependent spectral formulations, wavepacket dynamics, and open-system models to accurately simulate vibrationally resolved spectra.

The Herzberg–Teller approximation is a first-order vibronic expansion of the electronic transition dipole moment with respect to nuclear coordinates. In place of the Condon approximation, which assumes a coordinate-independent dipole, it retains the leading linear dependence of the dipole on normal modes and thereby provides a systematic description of non-Condon effects, intensity borrowing, and vibronic activation of electronically weak or symmetry-forbidden transitions. In the cited treatments, it is central to vibrationally resolved absorption and emission spectroscopy, coherent multidimensional spectroscopy, and open-system descriptions of excitonic dynamics (Prlj et al., 2020, Patoz et al., 2018, Troiani, 14 Jan 2026).

1. Formal definition and relation to the Condon approximation

In the Condon approximation, the electronic transition dipole moment is taken to be independent of nuclear geometry, for example

$\mu(Q)\approx \mu_0$

or, in the notation of one of the semiclassical formulations,

$\vec{\mu}(q)\approx \vec{\mu}(q_0).$

Under this assumption, vibronic intensities are controlled solely by Franck–Condon overlaps between vibrational states on the initial and final electronic surfaces (Prlj et al., 2020, Patoz et al., 2018).

The Herzberg–Teller approximation relaxes that constraint by expanding the transition dipole in a Taylor series around a reference geometry and retaining the linear term. In normal-mode coordinates,

$$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$

and, equivalently, in generalized coordinates,

$$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$

The first derivative term is the Herzberg–Teller, or first-order non-Condon, contribution. Higher-order terms may also be written, but the cited works consistently emphasize the linear truncation as the operative HT approximation (Prlj et al., 2020, Patoz et al., 2018, Zhang et al., 2016).

The physical content of the linear term is that vibronic motion modulates the electronic transition moment. This enables intensity borrowing from bright states, activates vibronic sidebands, and can make electronically forbidden transitions vibronically allowed when the zeroth-order dipole vanishes but its gradient does not (Patoz et al., 2018, Zhang et al., 2016).

2. Vibronic transition moments, intensity borrowing, and selection rules

For a transition between vibrational states on two electronic surfaces, the vibronic transition moment in the first-order HT approximation is

$$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$

When the amplitude is squared, the intensity separates into a Franck–Condon term, a pure Herzberg–Teller term, and an FC–HT cross term,

$$|\mathcal{M}|^2 \approx |\mu_0\langle\chi_{v'}|\chi_{v''}\rangle|^2 + \left|\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right|^2 + 2\,\mathrm{Re}\!\left[ \mu_0\langle\chi_{v'}|\chi_{v''}\rangle \left(\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right)^{\!*} \right].$$

This decomposition makes explicit that HT terms both add new intensity channels and interfere with the FC contribution (Prlj et al., 2020).

Because the coordinate operator connects different vibrational states, HT terms relax FC selection rules. In the harmonic-mode picture they activate transitions involving odd vibronic quanta and, more generally, transitions with nonzero matrix elements such as $\langle \chi_e|Q_i|\chi_g\rangle$. The cited works therefore describe HT coupling as a route by which electronically weak or symmetry-forbidden bands gain measurable intensity (Fang et al., 2022, Troiani, 14 Jan 2026).

One exact identity used to clarify intensity borrowing relates the gradient of the transition dipole to nonadiabatic couplings and dipole couplings: $$\frac{\partial}{\partial q_i}\vec{\bm{\mu}} = [\vec{\bm{\mu}},\mathbf{F}_i] + \left(e\sum_{j=1}^N Z_j \frac{\partial \vec{R}_j}{\partial q_i}\right)\mathbf{1}.$$ For off-diagonal transition elements near the ground-state equilibrium geometry, this reduces to the statement that a nonzero dipole gradient can arise from vibronic coupling to intermediate bright states. In that sense, the HT gradient is a direct formal expression of intensity borrowing (Patoz et al., 2018, Begušić et al., 2019).

The symmetry analysis of pyrrole provides an explicit selection-rule example. For the $$1^1A_2(\pi\sigma^*) \leftarrow \tilde{X}^1A_1(\pi\pi)$$ transition in $C_{2v}$, the direct dipole moment vanishes by symmetry, but first-order HT terms activate the transition when the dipole symmetry and vibrational symmetry combine to contain $\vec{\mu}(q)\approx \vec{\mu}(q_0).$0. In the notation given there, $\vec{\mu}(q)\approx \vec{\mu}(q_0).$1 vibrational distortions activate $\vec{\mu}(q)\approx \vec{\mu}(q_0).$2, $\vec{\mu}(q)\approx \vec{\mu}(q_0).$3 distortions activate $\vec{\mu}(q)\approx \vec{\mu}(q_0).$4, and $\vec{\mu}(q)\approx \vec{\mu}(q_0).$5 ring modes activate $\vec{\mu}(q)\approx \vec{\mu}(q_0).$6 (Picconi et al., 2017).

3. Time-dependent spectral formulations

A major virtue of the HT approximation is that it enters naturally in time-domain correlation-function formulations of spectra. For zero-temperature absorption, one cited formulation writes

$\vec{\mu}(q)\approx \vec{\mu}(q_0).$7

with

$\vec{\mu}(q)\approx \vec{\mu}(q_0).$8

For emission from $\vec{\mu}(q)\approx \vec{\mu}(q_0).$9, an analogous expression propagates the initial dipole-prepared state on the ground-state surface (Prlj et al., 2020).

In the Condon limit, $$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$0, so the autocorrelation reflects only nuclear dynamics on the final surface. With HT coupling,

$$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$1

and the correlation function acquires pure HT and FC–HT cross contributions involving coordinate operators at $$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$2 and during propagation. This is the basic time-domain origin of non-Condon structure in vibronic spectra (Prlj et al., 2020).

The same viewpoint is used in an explicitly tensorial orientational average beyond the Condon approximation. For coordinate-dependent dipoles, the isotropic average is written as

$$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$3

Under Condon conditions this reduces to the textbook $$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$4 rule based on $$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$5, but the tensor formulation remains valid in the HT case (Patoz et al., 2018).

At finite temperature, the correlation function becomes

$$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$6

and can be mapped, via thermo-field dynamics, to an autocorrelation in a doubled space,

$$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$7

with $$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$8. In the cited implementation, the HT dipole acts only on the physical coordinates, while the thermal purification is carried by the doubled Gaussian initial state (Begušić et al., 2020).

In coherent multidimensional spectroscopy, the same non-Condon mechanism appears as multiplicative oscillatory factors in the time-domain response. For the displaced harmonic oscillator model with linear HT coupling,

$$\mu(Q)\approx \mu_0 + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 Q_k + \text{higher orders},$$9

the response functions factor into an electronic term, a Franck–Condon vibrational factor, and finite sums of oscillatory HT corrections. Upon Fourier transformation, these oscillations generate replicas of FC multidimensional spectral features shifted by integer multiples of the vibrational frequency (Troiani, 14 Jan 2026).

4. Semiclassical, wavepacket, and ab initio realizations

Several cited works implement HT effects through wavepacket propagation. In the thawed Gaussian approximation, the nuclear wavepacket is propagated in a local harmonic approximation,

$$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$0

with center, width, and phase evolving according to the local gradient and Hessian of the potential (Patoz et al., 2018, Begušić et al., 2019).

The extended thawed Gaussian approximation propagates a Gaussian multiplied by a linear polynomial,

$$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$1

or, in the notation used for non-Condon spectra of azulene,

$$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$2

Within the local harmonic potential this functional form is preserved, so a single Gaussian dynamics yields both FC and HT contributions analytically and efficiently (Prlj et al., 2020, Patoz et al., 2018).

To reduce cost further, one implementation replaces the time-dependent Hessian in the width equations by a fixed adiabatic reference Hessian evaluated once at a final-state optimized geometry. This “single-Hessian” construction keeps the Gaussian center on the true force field while using constant curvature for widths. In the cited azulene study, the HT propagation formula remains valid under this approximation, and the method treats anharmonicity and mode mixing efficiently through a single on-the-fly ab initio trajectory (Prlj et al., 2020).

A related variant, the three thawed Gaussians approximation, represents the initial HT contribution not as a linear prefactor on one Gaussian but as an antisymmetric superposition of two displaced Gaussians, with a third Gaussian carrying the Condon term. The motivation stated there is that the HT packet has a node at the center and can be viewed as split into lobes in anharmonic potentials. Each Gaussian is propagated independently, and the spectra are assembled from their coherent superposition (Begušić et al., 2019).

For finite temperature, the doubled-space ETGA carries over with nearly no additional computational cost because the ground-state branch is typically stationary and the formal doubling does not require new on-the-fly electronic-structure calculations. The cited benzene application uses this framework to separate finite-temperature, anharmonicity, and HT effects in a symmetry-forbidden spectrum (Begušić et al., 2020).

The wavepacket perspective also underlies the multiconfigurational time-dependent Hartree treatment of pyrrole. There, coordinate-dependent transition dipole functions constructed from a linear HT expansion prepare the initial wavepacket, and the subsequent propagation on a locally diabatic Hamiltonian yields frequency-resolved absorption spectra and dissociation time scales for an optically dark $$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$3 state (Picconi et al., 2017).

5. Open-system and excitonic formulations

In excitonic open-system models, HT coupling is often represented by promoting the transition dipole to an operator that depends on bath coordinates. One cited DEOM formulation writes

$$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$4

where $$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$5 is the collective bath operator of an underdamped Brownian oscillator mode. In the Franck–Condon limit $$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$6; with $$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$7, the field couples directly to the vibrational coordinate during electronic excitation (Zhang et al., 2016).

For linear absorption of a two-level monomer in that framework, the cited analytical result is

$$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$8

The cross term and quadratic term explicitly redistribute intensity: increasing $$\vec{\mu}(q)\approx \vec{\mu}(q_0)+\mathrm{grad}_q\vec{\mu}|_{q_0}^{\top}\cdot (q-q_0).$$9 suppresses the zero-phonon peak and enhances vibronic sidebands, especially on the blue side. The same study emphasizes that non-Condon coupling intensifies dynamical electronic–vibrational energy transfer, enhances total system-and-bath quantum coherence, and strengthens long-lived vibronic oscillations in 2D spectroscopy (Zhang et al., 2016).

A closely related open-system construction in CODDE and DEOM uses a total dipole

$$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$0

which is the linear HT form expressed in system–bath language. In that formalism, non-Condon response functions are obtained through one-dissipaton auxiliary variables and correlated driving kernels, and the linear absorption correlation function contains both the reduced density operator and first-tier hybrid bath objects (Fang et al., 2022).

HEOM treatments extend the same idea to both coordinate-dependent transition dipoles and coordinate-dependent interstate couplings. The cited dimer Hamiltonian includes first- and second-order HT terms in the resonance coupling,

$$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$1

$$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$2

as well as a first-order non-Condon dipole

$$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$3

In HEOM space these terms introduce additional couplings between adjacent auxiliary density operators, and for sufficiently strong excitonic coupling the resulting linear absorption spectra can be captured very well by effective Huang–Rhys factors (Seibt et al., 2018).

6. Representative molecular applications and spectroscopic signatures

The cited applications show that HT effects are highly system dependent. For the phenyl radical, on-the-fly semiclassical spectra with and without HT are very similar and the HT contribution only slightly broadens peaks; the conclusion stated there is that anharmonicity outweighs HT corrections. By contrast, for benzene’s symmetry-forbidden $$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$4 transition, the Condon spectrum is zero because $$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$5, whereas the HT spectrum agrees well with experiment and reproduces intensities of the main progression within about $$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$6 of experimental values (Patoz et al., 2018).

Azulene provides a different non-Condon regime. In the cited study, $$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$7 absorption is FC-dominated and HT does not improve the band appreciably, while $$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$8 absorption and $$\mathcal{M}_{v'v''} = \langle \chi_{v'}(e')| \mu(Q) |\chi_{v''}(e'')\rangle \approx \mu_0 \langle \chi_{v'}|\chi_{v''}\rangle + \sum_k \left(\frac{\partial \mu}{\partial Q_k}\right)_0 \langle \chi_{v'}|Q_k|\chi_{v''}\rangle.$$9 emission require HT to reproduce the experimental vibronic structure and intensity. The same work states that the breaking of absorption–emission mirror symmetry for $$|\mathcal{M}|^2 \approx |\mu_0\langle\chi_{v'}|\chi_{v''}\rangle|^2 + \left|\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right|^2 + 2\,\mathrm{Re}\!\left[ \mu_0\langle\chi_{v'}|\chi_{v''}\rangle \left(\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right)^{\!*} \right].$$0 arises predominantly from HT terms interacting with Duschinsky mode mixing, and that turning off either HT or mode mixing largely restores symmetry (Prlj et al., 2020).

Pyrrole illustrates HT-enabled access to dissociative dark states. The $$|\mathcal{M}|^2 \approx |\mu_0\langle\chi_{v'}|\chi_{v''}\rangle|^2 + \left|\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right|^2 + 2\,\mathrm{Re}\!\left[ \mu_0\langle\chi_{v'}|\chi_{v''}\rangle \left(\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right)^{\!*} \right].$$1 transition is electric-dipole forbidden at the $$|\mathcal{M}|^2 \approx |\mu_0\langle\chi_{v'}|\chi_{v''}\rangle|^2 + \left|\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right|^2 + 2\,\mathrm{Re}\!\left[ \mu_0\langle\chi_{v'}|\chi_{v''}\rangle \left(\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right)^{\!*} \right].$$2 equilibrium geometry, but becomes weakly allowed through vibronic coupling. In the 15-dimensional treatment cited there, the total absorption maximum is at $$|\mathcal{M}|^2 \approx |\mu_0\langle\chi_{v'}|\chi_{v''}\rangle|^2 + \left|\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right|^2 + 2\,\mathrm{Re}\!\left[ \mu_0\langle\chi_{v'}|\chi_{v''}\rangle \left(\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right)^{\!*} \right].$$3, the envelope width is $$|\mathcal{M}|^2 \approx |\mu_0\langle\chi_{v'}|\chi_{v''}\rangle|^2 + \left|\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right|^2 + 2\,\mathrm{Re}\!\left[ \mu_0\langle\chi_{v'}|\chi_{v''}\rangle \left(\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right)^{\!*} \right].$$4, and the polarization-resolved intensities satisfy $$|\mathcal{M}|^2 \approx |\mu_0\langle\chi_{v'}|\chi_{v''}\rangle|^2 + \left|\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right|^2 + 2\,\mathrm{Re}\!\left[ \mu_0\langle\chi_{v'}|\chi_{v''}\rangle \left(\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right)^{\!*} \right].$$5. The spectrum is diffuse, structured by ring-mode progressions together with obligatory H-bend excitation, and the associated photodissociation is predominantly direct with short survival times (Picconi et al., 2017).

In strongly coupled vibronic dimers treated by HEOM, first-order HT coupling enhances the vibrational progression in one exciton band and diminishes it in the other, while non-Condon dipoles redistribute intensities without changing vibronic spacings. The same study reports that effective Huang–Rhys factors provide a compact description of these HT-induced spectral changes when excitonic coupling is sufficiently strong (Seibt et al., 2018).

Across these applications, a recurring diagnostic pattern is that HT terms become decisive when a transition is electronically forbidden or weak at the reference geometry, when absorption–emission mirror symmetry is broken, or when vibronic sidebands and intensity patterns cannot be reproduced within FC-only models (Patoz et al., 2018, Prlj et al., 2020).

7. Validity regime, approximations, and practical use

All cited implementations emphasize that the standard HT approximation is a truncation at first order in the dipole expansion. It is accurate when the transition dipole varies moderately over the support of the nuclear wavepacket. Strongly nonlinear $$|\mathcal{M}|^2 \approx |\mu_0\langle\chi_{v'}|\chi_{v''}\rangle|^2 + \left|\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right|^2 + 2\,\mathrm{Re}\!\left[ \mu_0\langle\chi_{v'}|\chi_{v''}\rangle \left(\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right)^{\!*} \right].$$6 requires higher orders, and several works explicitly note that second-order vibronic coupling or quadratic HT terms are omitted from the first-order treatments (Prlj et al., 2020, Patoz et al., 2018, Begušić et al., 2020).

The wavepacket and semiclassical realizations impose further assumptions. The ETGA and related thawed-Gaussian methods assume a locally harmonic potential around the trajectory center, and single-trajectory versions cannot describe strong wavepacket splitting or strongly anharmonic long-time dynamics. The single-Hessian variant is most reliable in near-harmonic regions, although it still captures part of the anharmonicity through propagation of the packet center on the full ab initio surface (Prlj et al., 2020, Patoz et al., 2018).

Open-system approaches impose different restrictions. The DEOM and CODDE formulations use harmonic baths and Gaussian statistics; the HEOM implementation discussed for HT and non-Condon effects relies on exponential decompositions of bath correlation functions and requires increasing hierarchy depth for convergence, especially for undamped oscillators and second-order HT terms (Zhang et al., 2016, Fang et al., 2022, Seibt et al., 2018).

Within these limits, the practical guidance given in the cited studies is consistent. HT terms should be included when electronically weak or symmetry-forbidden bands are of interest, when mirror-image symmetry between absorption and emission is broken, or when nearby bright states and vibronic couplings suggest strong intensity borrowing. Mode mixing should be treated accurately, $$|\mathcal{M}|^2 \approx |\mu_0\langle\chi_{v'}|\chi_{v''}\rangle|^2 + \left|\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right|^2 + 2\,\mathrm{Re}\!\left[ \mu_0\langle\chi_{v'}|\chi_{v''}\rangle \left(\sum_k \mu_k' \langle\chi_{v'}|Q_k|\chi_{v''}\rangle\right)^{\!*} \right].$$7 and its gradient should be evaluated at the appropriate reference geometry, and in single-surface spectral simulations the Born–Oppenheimer assumption should be validated independently when nonadiabatic effects may be significant (Prlj et al., 2020, Patoz et al., 2018).