Spectroscopy from quantum dynamics: a mixed wave function/analytical line shape functions approach

Quantum dynamics is the natural framework in which accurate simulation of spectroscopy of nonadiabatically coupled molecular systems can be obtained. Even if efficient quantum dynamics approaches have been developed, the number of degrees of freedom that need to be considered in realistic systems is typically too high to explicitly account for all of them. Moreover, in open-quantum systems, a quasi-continuum of low-frequency environment modes need to be included to get a proper description of the spectral bands. Here, we describe an approach to account for a large number of modes, based on their partitioning into two sets: a set of dynamically relevant modes (so-called active modes) that are treated explicitly in quantum dynamics, and a set of modes that are only spectroscopically relevant (so-called spectator modes), treated via analytical line shape functions. Linear and nonlinear spectroscopy for a realistic model system is simulated, providing a clear framework and domain of applicability in which the introduced approach is exact, and assessing the error introduced when such a partitioning is only approximate.


Introduction
Simulation of linear and nonlinear spectroscopy of electronically excited molecular systems is of paramount importance to understand complex experimental outcomes, untangle spectral contributions, and unravel the microscopic origin of the observed spectral features.Such simulations require: a) accurate computations of the microscopic (quantum chemical) inputs that characterize the molecular system of interest; b) an approach to describe the (generally nonadiabatic) system dynamics after photoexcitation (as, e.g., semiclassical or purely quantum dynamics (QD) approaches); c) a mean to describe the light-matter interaction (either via direct inclusion of the laser fields in the Hamiltonian or via perturbative methods, such as the Mukamelian response function approach [1,2]).Some of the present authors have recently presented a QD/spectroscopy approach [3] that, by recasting the firstand third-order response functions in terms of wave function overlaps, produced expressions for linear and nonlinear spectroscopy in a QD framework.A linear vibronic coupling (LVC) model Hamiltonian of the pyrene molecule, parametrized via ab-initio CASSCF/CASPT2 inputs, was employed as testbed case, where up to 15 nuclear DOFs (out of 49) 1 were accounted for in the quantum dynamics [4].The QD approach employed in Ref. [3] was multi-configuration time-dependent Hartree (MCTDH), an efficient method that can take into account tens of degrees of freedom (DOFs).
108 Page 2 of 17 Even if approaches capable of treating up to a hundred of DOFs exist (such as the multi-layer version of MCTDH (ML-MCTDH) [5][6][7], or tensor train methods [8]), one realizes that for increasing system size, not all the modes that are relevant for spectroscopy simulation can be included in the quantum dynamics.Interestingly, not all of them need to be explicitly included.The ones that must be accounted for are coupling and tuning modes that significantly couple different electronic states and/or change their energy gap.But there might be modes that are neither coupling nor tuning modes, but that can still contribute in shaping the spectral bands.In this contribution, we start from a partition of these two sets of modes (the dynamically relevant ones, and the ones that, while not relevant for the dynamics, are still relevant for the evaluation of spectroscopy) to show that the latter can actually be treated analytically.
Practically, one parametrizes a reduced model Hamiltonian based on a restricted number of modes that govern the dynamics of the photoexcited states, including the left out modes that can contribute to spectroscopy only in the computation of the response function.Let us consider a simple example, transient spectroscopy: here, a first pulse (the socalled pump pulse) creates a linear combination of states over which the nonadiabatic dynamics occurs [a manifold of states that we label with E , see Fig. 1a)]; a second pulse, the probe pulse, might further excite the system in another manifold of electronic states [labeled as F , Fig. 1a)].There might be modes for which the states in the E manifold do not have a (significant) relative displacement and/or coupling, but along with some of the states in the F manifold display a significantly different displacement [see Fig. 1c)]: these modes would indeed contribute to the vibronic broadening of E → F transitions, while being irrelevant in the E mani- fold quantum dynamics.
In what follows, we first define the framework in which such an approach that mixes explicit numerical QD propagation along some modes, and analytical line shape functions evaluation of other modes is rigorously justified.In passing, we note that approaches that couple numerical solutions of the Schrödinger equation along some modes with analytical treatment of other coordinates exist: as an example, in Ref. [9], the authors considered (one-dimensional) dissociative coordinates via explicit wave packet propagation in time, and the remaining modes via the Franck-Condon overlaps (in the frequency domain).Still, the approach derived here provides a clear and general theoretical framework for rigorously partitioning the modes.
The protocol is quite general and can be employed in tandem with whatever QD method capable of producing the wave packet overlap between different electronic states.Since we point to work out a model for large systems with many degrees of freedom, the MCTDH approach was chosen.The dynamics was computed employing the package Quantics [10,11].We derive response function expressions for linear and nonlinear spectroscopy simulation within the mixed QD/analytical approach and compare them with the exact full QD expressions.Finally, we design a simple model system that we employ to test the derived equations against the simulation of spectroscopy, and we also test the robustness of the approach in cases where the partitioning scheme is not strictly exact (i.e., when some of the modes left out from the Hamiltonian that drives the system evolution would have contributed to the dynamics if included in the QD propagation).

Hamiltonian, wave function, approximations
Here, we start from a LVC nonadiabatic molecular Hamiltonian written as (1) where three manifold of electronic states are considered: the ground state manifold G (that only contains the GS, g), the E manifold in which the nonadiabatic dynamics occurs, and the F manifold of the final states probed by the probe pulse (see Fig. 1); E (ad) a is the adiabatic energy of the electronic state a.The terms Ĥa,A∕S and Ûee ′ ,A contain the dependence over the nuclear coordinates, and the vibrational modes are partitioned into a set of active modes (subscript A), explicitly considered in the QD, and a set of spectator modes (subscript S) that can be left out when selecting the relevant modes for the QD.The various terms (in dimensionless coordinates) read where the index m runs over the active modes, while s runs over the spectator modes.Moreover, the index a denote a generic state of the three manifolds, while indexes g, e∕e � , f denotes states that belong to the G , E and F manifolds, respectively.Tl is the kinetic energy operator along mode l, 2 and the potential term is simply given by a displaced harmonic potential, whose displacement (for the given electronic state a and along the given mode l) is d a,l .With- out loss of generality, we can assume d g,m = d g,s = 0 and E (ad) In what follows, we will consider projections of the Hamiltonian in Eq. ( 1) onto the various electronic states manifolds, i.e., We assume that electronic states belonging to different manifolds are coupled only via dipole interaction (i.e., only the external electromagnetic field can promote a change of manifold), 3 which is given (in the Condon approximation) by: (2) and it is here useful to define some projections of the dipole moment onto specific manifolds, i.e., Indeed, the term ∑ e≠e � ∈E Ûee � ,A �e⟩⟨e � � in Eq. ( 2) is capable of promoting population transfer only between the electronic states in the E manifold; this term is assumed to be linear in Q, and the strength of the coupling between these states is controlled by the constants ee ′ ,m (which are assumed to be not all null).Therefore, E is the only manifold in which we allow the nonadiabatic dynamics to take place.
Equation 1 prescribes that the spectator modes cannot promote (electronic state) population transfer: Indeed, we are assuming that in the construction of the LVC Hamiltonian for the quantum dynamics, all the relevant coupling modes have been included in the set of the active modes.Moreover, we would also like to include all relevant tuning modes in the active modes set, as they do influence the system dynamics.This last choice corresponds to assume that for the spectator modes, Ve,s = Ve � ,s ∀ e, e � ∈ E and ∀ s ∈ S , i.e., d e,s = d e � ,s ∀ e, e � ∈ E, ∀ s ∈ S .We will therefore write Ĥe,s = Ĥē,s ∀ e ∈ E and ∀ s ∈ S .The subscript ē indicates that the Hamiltonian does not depend on the specific state of the E manifold (i.e., it is identical for all of them).This implies that and the term PE = ∑ e∈E �e⟩⟨e� is simply the projector onto the E manifold.Under this assumptions, the E sector of the Hamiltonian, reported in Eq. ( 3), can be rewritten as where the subscript el, A and S indicate the electronic, active-vibrational and spectator-vibrational sectors of the Hamiltonian, and ( ) ĤE,A , ĤE,S = 0 2 In dimensionless coordinates the kinetic operator is given by: Tl = 1  2  l P2 l , Pl being the momentum operator along the lth coordinate. 3If, on the timescales of interest and for a given molecular system, the process that brings an excited electronic state back to the ground electronic state is relevant, one may add a term in the Hamiltonian that couples these two manifolds.i.e., the A and S sectors of the Hamiltonian commute (which would have not been the case if the various E manifold PESs were differently displaced along the spectator modes).In passing, we note that this commutation relation between the terms that depend on A and on S holds for all the three manifolds.
The assumptions employed so far are summarized hereafter (see also Fig. 1): • LVC model: The potential energy surfaces are harmonic and coupled only via linear terms (in the normal modes coordinates); furthermore, we have assumed that: -ab,l = 0 ∀ l ∈ A∕S , if a and b belongs to different manifolds: the three manifold of electronic states G , E and F are not coupled in the molecular Hamilto- nian (they are indeed only dipole coupled via the field-matter interaction Hamiltonian); -ff � ,l = 0 ∀ f , f � ∈ F, ∀ l ∈ A∕S : no population transfer can occur in the F manifold.Note that this assump- tion may look very strong, especially for high-lying states that are typically packed in a dense manifold with many potential crossings along many coordinates.Precisely because of that the E − F coherences (responsible for the signals in, e.g., transient absorption) will have a very short lifetime, i.e., a very broad and unstructured spectrum that makes the explicit evaluation of their nonadiabatic dynamics immaterial [3].-ee � ,s = 0 ∀ e, e � ∈ E, ∀ s ∈ S : the spectator modes do not couple the electronic states of the E manifold; -d e,s = d e � ,s ∀ e, e � ∈ E, ∀ s ∈ S : the electronic states in the E manifold are identically displaced along each of the spectator modes s ∈ S; • Concerning the dipole coupling with the external electromagnetic field, only transitions from G → E and from E → F are allowed (Eq.4).Overall, it should be clear that the separation of the electronic states in G , E and F manifolds is determined by the characteristics of the pump and probe pulses (as, e.g., their central frequency and bandwidth)4 ; • The Condon approximation is employed (i.e., ab (Q) = ab ), which is a sensible approximation in the diabatic (electronic states) framework here adopted.This completes the description of the framework in which we are going to derive the mixed active/spectator modes equations for spectroscopy.Note that the assumption d e,s = d ē,s ∀ e ∈ E, ∀ s ∈ S does not tell anything about the magnitude of such d ē,s , or about the magnitudes of d f ,s along the various modes and for the various f states.This is exactly why these spectator modes, even if not relevant for the nonadiabatic dynamics, may be relevant for spectroscopic observables.In what follows, we assess the impact of these modes in the response functions, and thus in linear and nonlinear spectra.

Wave function and Response functions
The natural place where this machinery should be applied is transient spectroscopy: The first pulse prepares the system in the E manifold, where the nonadiabatic dynamics takes place, while the probe pulse promotes the time-evolved wave packet either in the GS or in the F manifold.The cen- tral frequency of the two pulses determines the states that should be placed in the various manifolds.In UV/Vis transient spectroscopy, both E and F manifolds are made of (neutral) valence excited states; if the probe pulse is higher in energy, the F manifold can include ionized states (time- resolved photoelectron spectroscopy); if the probe lies in the X-ray domain, core-excited or core-ionized states will be considered (time-resolved XANES and time-resolved XPS, respectively).
Nonetheless, we also consider linear absorption, where the LVC modeling might include all the relevant coupling and tuning modes, and still have left out some relevant spectator modes, which have a negligible relative displacement along the E PESs, but a non-negligible displacement with respect to the GS PES minimum.
By following the steps of Ref. [3], the molecular wave function at a given time t is defined as: where � a (t, Q)⟩ represents the nuclear wave packet evolving on the a-th electronic state � a (q;Q)⟩ which can be repre- sented in a time-dependent or -independent basis [12]. 5q and Q denote the electronic and nuclear degrees of freedom, respectively.By employing a diabatic formulation, the electronic wave function becomes (ideally) independent on the nuclear coordinates, i.e., � a (q;Q)⟩ = � a (q)⟩ .The condition ⟨Ψ(Q, q, t)�Ψ(Q, q, t)⟩ = 1 holds for each time t.The elec- tronic wave function � a (q)⟩ is conveniently rewritten as �g⟩ , ( 9) �e⟩ and �f ⟩ for states in the three manifolds; in what follows, we will drop the explicit coordinate dependence, so that: The complete WP on each electronic state a is simply the Hartree product of the WP along the two sets of modes (A and S), i.e., so that the wave function can be written as: . As shown in Ref. [3], the response functions can be expressed in terms of the wave packets obtained in a QD simulation; this is the starting point of our derivation.Note that, in what follows, we work in atomic units, where ℏ = 1 .

Linear absorption
The first-order response function reads [3] where ⟨ g (t)� e→e � (t)⟩ denotes the time-dependent over- lap between the ground state WP (bra side) and a timeevolved WP on the E manifold (ket side).In particular, the symbol � e→e � (t)⟩ indicates that fraction of the nuclear WP that was initially prepared in e at t = 0 , and that is now found in e ′ at time t.This means that, at each time t, e eg � e→e � (t)⟩ , i.e., the total nuclear WP at time t on the electronic state e ′ is simply given by the (dipole- weighted) sum over all the terms that, initially prepared in the various electronic states e at time t = 0 , contribute to ( 10) shape the WP of state e ′ at time t (N is a normalization fac- tor that is easily obtained from the set of eg dipoles, where eg = ⋅ eg , being the pump electric field polarization).By taking advantage of the projector-like nature of the dipole moments and of the block diagonal (manifolds driven) definition of the total Hamiltonian, we can conveniently rewrite Eq. ( 13) as Aimed at partitioning the response function into an active modes expression, and a spectator modes expression, the Hamiltonian operators are conveniently expressed in terms of the active and spectator modes.Due to the null commutation between the A and S sectors of the Hamiltonian, one has that Moreover, we observe that: so that the obtained terms Ĥg,S and Ĥē,S do not only commute with the A projected Hamiltonians, but also with the transition dipole moment operators (as Ĥg,S and Ĥē,S are operators that only act on the nuclear space).One eventually gets i.e., we have factored all the S-dependent operators on one side of the operator expression.
We now rewrite Eq. ( 14) (dropping the i ℏ (t) prefac- tor) taking advantage of the expression obtained in Eq. ( 17) and employing the wave function definition given in Eq. ( 12), obtaining ( 14) R (1) where in the last step we have defined �Ψ a,A (t)⟩ as the pro- jection of the wave function onto the A modes subspace.Very interestingly, one notices that the first term of the final expression is the response function for a quantum dynamics simulation performed on the sole A modes.
The second term of Eq. ( 18) is the time-dependent nuclear WP overlap between the � g,S (0)⟩ WPs propagated on the ground and the excited state potential energy surfaces.Recasting all the pieces together, it becomes apparent that the first-order response function from a quantum dynamics simulation of the full A + S system can be exactly recast in a A-dependent response function term (denoted as QD(A)), and a S-dependent WP overlap, i.e., where we have rewritten ⟨ g,S (0)�e i Ĥg,S t e −i Ĥē,S t � g,S (0)⟩ as ⟨ g,S (t)� ē,S (t)⟩ , by leveraging on the fact that � g,S (0)⟩ = � ē,S (0)⟩ due to the impulsive nature of the laser field that creates identical copies of the WP from the starting electronic state to the arrival electronic state.
We thus arrived at the core of the present contribution.The two obtained terms can be evaluated separately: the one that involves the active modes, via a proper QD simulation performed on the sole A modes (that will take care of the nonadiabatic dynamics driven by the active modes), while the one that depends on the spectator modes, can be solved analytically.
Let us focus on the term that pertains the spectator modes and show how to rewrite it by means of an analytical expression.In the zero temperature limit, one could consider � g,S (0)⟩ being the Hartree product of the lowest energy vibrational states of each s ∈ S .Still, a more gen- eral derivation is possible that also accounts for temperature effects: In particular, one could consider an incoherent sum over the vibrational eigenstates, weighted by Boltzmann contributions, i.e., rewrite this expression in terms of the equilibrium GS density matrix ρg,S of the spectator modes.The equilibrium GS density matrix is defined as where = 1∕(K B T) ( K B being the Boltzmann constant, and T being the temperature), n s run over the vibrational states of the mode s, and the overlap of Eq. ( 19) can now be generalized to account for thermal effects, obtaining: Taking Ĥg,S as reference, the time evolution of such operator can be evaluated in the Heisenberg picture as: Making use of this expression, the e −i Ĥē,S t propagator can be rewritten in terms of the energy gap fluctuation as: where the time-ordered exponential of the energy gap fluctuation is defined as: Now plugging Eq. ( 24) into equation ( 21), one gets rid of the e −i Ĥg,S t propagator obtaining the following expression for the WP overlap along the spectator modes: The last term in Eq. ( 26) can be calculated with the cumulant expansion method proposed by Prof. Ryogo Kubo [13], by leveraging on formal properties of statistical distributions.Notably, the coupling with a bath of harmonic vibrations (like in the LVC model employed here) induces Gaussian fluctuations of the δgē,S energy gap, for which a second-order truncation of the cumulant expansion is exact.Therefore, the generalized (ensemble) WP overlap along the spectator modes can be rewritten simply as: where g(t) is the so-called line Finally, the advantage of this formulation became apparent since both Eqs.(28) and (29) have an analytical expression which exploits the relation between the so-called spectral density J( ) and the two-time correlation function of the energy gap: where the spectral density for a collection of displaced harmonic oscillator (with dimensionless displacements) reads: The a and b indexes refer to different electronic states.In the specific case of the linear absorption investigated here, however, the g ēē (t) entirely encodes the response of the system along the spectator modes, and thus it can be used to rewrite the expression for the linear response function, obtaining: We remark the fact that, under the assumptions detailed in the previous section, this expression is exact, i.e., it gives the same results of Eq. ( 13) (where both A and S modes are included in the quantum dynamics).

Transient absorption
The third-order response for a pump-probe transient absorption experiment reads [3]: where t 2 denotes the time interval between the pump and probe pulse, while t 3 indicates the time interval between the probe pulse and the detection of the signal.In writing Eq. (33), we have removed the common (i∕ℏ) 3 (t 2 ) (t 3 ) prefactor, and we have expressed both Hamiltonian and dipole operators as acting in/between specific manifolds.Equation (33) is derived in the Condon approximation and in the impulsive limit (extremely short duration of the laser pulses). 6  (29) e −g ēē (t) (33) The three terms of Eq. ( 33) are, in order: the ground state bleaching (GSB), the stimulated emission (SE), and the excited state absorption (ESA) [2,3].We label the three contribution to the total third-order response function as , respectively.Let us analyze these three terms one by one and follow similar steps that brought us, in the linear absorption case, from Eqs. ( 13)-(32).
GSB: This term closely resembles the one analyzed for the linear absorption.After performing the A/S partitioning, the S-dependent overlap term reads: where we again exploit the fact that Hamiltonian operators defined in different manifolds commute.Thus, one eventually gets: As before, the term that depends on the S can be generalized to consider a thermal distribution of vibrational states and return the same analytical expression we found for the linear absorption.

SE:
For the SE, one similarly has that the S-dependent overlap term in the SE response function contribution is given by which is the overlap between the bra side WP (that, after a t 2 propagation in the E manifold, is then projected back to the ground state and there evolved along t 3 ) and the ket side WP (which evolved on the E manifold for the full t 2 + t 3 time interval).Following the same procedure discussed for the linear absorption, and generalizing the overlap term as an ensemble average, this expression can be reformulated inserting the energy gap fluctuation δgē,S (t) and employing the second-order cumulant expansion, so that: (34) ⟨ g,S (0)�e +i Ĥg,S (t 2 +t 3 ) e −i Ĥē,S t 3 e −i Ĥg,S t 2 � g,S (0)⟩ = =⟨ g,S (0)�e +i Ĥg,S t 2 e −i Ĥē,S t 3 � g,S (0 (36) ⟨ g,S (0)�e +i Ĥē,S t 2 e +i Ĥg,S t 3 e −i Ĥē,S( t 2 +t 3 ) � g,S (0)⟩ = =⟨ ē;g,S (t 2 ;t 3 )� ē,S (t 2 + t 3 )⟩ 108 Page 8 of 17 In particular, the function  ēgē (t 2 , t 3 ) is a combination of line shape functions g ij (t) obtained integrating different two-time correlation function of the energy gap fluctuation, i.e., where * is denotes the complex conjugate of the line shape function, and many terms that were functions of the GS displacement (that we set to zero, without loss of generality) simply vanish.Thus, we have obtained that ESA: A similar procedure leads to the corresponding S-dependent ESA term, given by: which depends on both the ē displacement, and the f states displacements.Following similar steps to those reported for the SE contribution, it is possible to write the four-point line shape function as a sum of two-point g(t) functions, i.e., Note that such an expression is much more complex than the one found for the SE, as it involves displacements of e and f states, that, in principle, are both non-vanishing.Therefore, the final expression for the transient absorption in the mixed QD/analytical approach becomes: and is exact in the presented framework.Fully worked out expressions for the various R (3)  QD(A) (t 3 , t 2 ) contributions are reported in the supporting information.
LA and TA spectra are obtained by Fourier transformation of the first-and third-order response functions (with respect to t and t 3 , respectively). (37)

Results
In this section, we design a simple model system and use it to compare the results obtained by a) including the S modes explicitly in the QD (labeling this approach QD(A + S) ), b) including them via the analytical line shape functions ( QD(A) * An(S) approach), and c) just ignoring them (QD(A) approach, which is what would happen when considering only a limited number of modes in the system model for QD, and ignoring the effect of all the left out modes).
A few assumptions are made to simplify the problem: we consider an oriented sample (i.e., no orientational averaging is considered here); in transient absorption, the pump and probe pulses are assumed to be extremely narrow in time (impulsive limit); finally, in the evaluation of the analytical line shape functions, we took the zero temperature limit (i.e., we consider the → ∞ ); this allows one to directly compare the QD(A) * An(S) and the full QD(A + S) approach, which is indeed run in the T = 0 limit.In the evaluation of the response functions Fourier transform, an additional damping exponential term exp (−t∕(2 )) is multiplied to the expres- sions (regardless the level of theory); this term summarizes the effect of the many environment (or bath) modes that are not explicitly considered in the dynamics.Note that such bath modes can be explicitly considered (together with their temperature dependence) in the QD(A) * An(S) approach.This will be discussed in the Discussion section.Note that both orientational averaging and realistic pulses can be straightforwardly accounted for [14].
After assessing the identity between the full QD approach and the QD/analytical approach, we explore a region of parameters for which d e,s ≠ d e ′ ,s for different states in the E manifold, i.e., we relax the requirement for them to be identical.Still, we require the commutation relation to be close to zero (even if not strictly zero), i.e., ĤE,A , ĤE,S ∼ 0 .In other words, we test the robustness of the QD/analytical equations considering a small (but non-vanishing) commutator between the two sets of A and S modes.
Additional computational details are reported in the SI.

Model parameters
A simple model system in which all the ingredients described above are taken into account was designed.It comprises: • Five electronic states: the ground state g, two states in the E manifold ( e 1 and e 2 ), and two states in the F manifold ( f 1 dipole coupled to e 1 , f 2 to e 2 ).We assume that only the g → e 2 transition is bright from the GS; • Three modes: one coupling mode ( m 1 ), one tuning mode ( m 2 ), and one spectator mode ( s 1 ).
Page 9 of 17 108 The model parameters are reported in Table 1.Such parameters were chosen in order to have a relatively fast population dynamics, obtained by tuning the coupling between the two E manifold states; a significant role of the spectator mode in the spectra (including or excluding the spectator mode shall significantly affect the spectroscopic observables), obtained imposing a large displacement d ē,s 1 along the s 1 mode that will translate in a pronounced vibronic progression in the spectra.Note that, for this exploratory work, we assumed the transition dipole moments (for the bright transitions) to have the same orientation and module equal to 1.

Population dynamics
Figure 2 reports the first 300 fs of the population dynamics that follows the e 2 state photoexcitation, for both the QD(A + S) and the QD(A) cases.
Under the assumptions we made (form of the Hamiltonian prescribed by Eqs. ( 1) and ( 2), where the E manifold states have equal displacements of along the spectator modes (Eq.( 6)) and the nonadiabatic coupling between states in the E manifold and the GS is set to zero), the QD(A + S) and QD(A) population dynamics are identical.This confirms that including or excluding the spectator modes from the QD does not affect the dynamics itself.In the supporting information, we also reported a comparison between QD(A + S) results obtained with the MCTDH dynamics and with a numerically exact propagation method (Figure s1).

Spectra
We first analyze the linear absorption spectrum of the model system reported in Table 1.Transient spectra are also shown in the following Section.

Linear absorption
Figure 3 shows the time-dependent WP overlap obtained in different scenarios: QD(A + S) model (considering both active and spectator modes in the QD simulation), QD(A) * An(S) model [considering the active modes in the QD and the spectator mode with the analytical expression given in Eq. (30)], and QD(A) model (where the spectator modes are ignored).We also reported the analytical spectator mode(s) contribution, labeled as An(S).Note that, at variance with the population dynamics, the QD(A + S) and QD(A) are different.The full modes model QD(A + S) is instead identical to the QD(A) * An(S) , as predicted in the derivation of Eq. (32) reported in the previous section.
Figure 4 shows the comparison between QD(A + S) , QD(A) * An(S) and QD(A) linear absorption spectra.QD(A + S) and QD(A) spectra were obtained employing (the Fourier transform of) Eq. ( 13), while the QD(A) * An(S) spectrum was obtained via Fourier transform of Eq. (32). 7ince the LA is essentially the Fourier transform of the overlap (centered at the vertical energy gap between the GS and the bright e 2 state), what we get here mirrors the results of Fig. 3: QD(A + S) and QD(A) * An(S) are identi- cal, while QD(A) is significantly different.In particular, the role of the spectator mode s 1 is that of introducing addi- tional broadening, in the form of a vibronic progression that Dipole coupling between pair of states.Note that e 1 is not dipole coupled to g.All the non-null transition dipoles are assumed to have the same direction and module 1 Fig. 2 Population dynamics of the photoexcited state e 2 (red curve) and of the dark state e 1 (blue curve), within the model system reported in Table 1, obtained with a QD (MCTDH) simulation that comprises both active and spectator modes, referred to as QD(A + S) , compared with the QD(A) population dynamics of the same states (black dashed curves).The curves at the two levels of theory coincide 108 Page 10 of 17 originates from the significant shift of the E manifold states with respect to the GS along such mode.Figure s3 of the supporting information shows the Fourier transform of the analytical contribution of the S mode(s) alone, which indeed presents well resolved vibronic bands with ℏ s 1 spacing and peaks height determined by the (common E manifold state) displacement d ē,s 1 .
In the supporting information, Figure s2, we also report the comparison of the QD(A + S) linear spectrum obtained via MCTDH quantum dynamics and via numerically exact propagation where the wave packet is directly expanded on the time-independent primitive basis set, and the timedependent coefficients are obtained with a short iterative Lanczos propagation scheme [10,15].The two curves are identical.

Transient absorption
In this section, we consider transient absorption.Here, the system interacts with two laser pulses, with a controlled time delay, t 2 , that is varied in the interval t 2 ∈ (0, 200) fs.The pump pulse promotes the WP in the (bright states of the) E manifold, where the various PESs are potentially displaced with respect to the GS, and coupled among each other.This means that, the WP, that originally resided in the GS, after photoexcitation is out of equilibrium in the E PESs, therefore evolves in time driven by the purely intra-state nuclear dynamics (along both the A and the S modes), and by the inter-state nonadiabatic (linear) coupling.After a t 2 time interval, the probe pulse promotes the t 2 evolved wave packet either back to the GS (giving rise to SE signals), or to some higher lying states of the F manifold (giving rise to ESA signals).The spectrum is then obtained by Fourier transformation of the t 3 time-dependent overlaps between the ( t 2 evolved) wave packet residing in different states.Very importantly, along the spectator modes, the analytical line shape functions presented in Eqs. ( 38) and (41) are capable of capturing both the t 2 and t 3 -dependent dynamics in an exact way [1,2].
We still compare the three scenarios as we did for the linear absorption: the full quantum dynamics model, QD(A + S) , the mixed QD and analytical approach, QD(A) * An(S) , and the approach in which the spectator modes are simply neglected, QD(A).The QD(A + S) and QD(A) spectra are obtained via Fourier transformation (along t 3 ) of Eqs.s1-s3 of the supporting information, while for the QD(A) * An(S) spectrum we also considered the addi- tional contributions specified in Eq. ( 42). 8 The model system parameters are reported in Table 1.
Figure 5 shows the comparison of the three levels of theory (Figure s5 of the supporting information shows the separate TA contributions).As a first general observation, the spectra of the full QD approach and of the mixed QD/ Fig. 3 WP overlaps: full QD treatment, QD(A + S) (only real part, red curve, bottom); active modes in the quantum dynamics and spectator modes treated analytically, QD(A) * An(S) (only real part, black dashed curve, bottom); only active modes considered in the quantum dynamics, QD(A) (real and imaginary parts, blue and cyan curves respectively, middle); analytical line shape function expression ( e −g ēē (t) ) for the spectator mode, An(S) (real and imaginary parts, dark and light green curves respectively, top).Some curves have been vertically shifted to avoid congestion of the data.The mixed QD/analytical curve coincides with the full QD curve.The E (ad) e 2 energy gap was here down-shifted to 0.8 eV (from 4.02 eV) as this makes it easier to visually compare the coherences (avoiding excessively high-frequency oscillations) analytical approach are, as expected, identical.A few signals can be distinguished in all three panels: in the center of the spectrum, around 4 eV, one recognizes the (positive) GSB contribution, which is constant along t 2 (at least in the ideal case of impulsive interaction between the laser and the molecules of the sample).This signal mirrors the LA spectrum, and is therefore much narrower in the QD(A) case, while the vibronic progressions due to the s 1 mode are only present in the first two panels.On the top part of the spectrum, the ESA signal from the e 2 state is clearly visible as a blue feature that extends from around 5 eV to 7 eV.This contribution quickly decreases in intensity in the first 50 fs (as dictated by the decreasing e 2 population), and then stabilizes.Note that signatures of WP bifurcation within the e 2 well are apparent in all three spectra in the time window between 100 and 150 fs.In the bottom part of the spectrum, ESA from the e 1 state can be seen, which, at variance with the previous signals, are null at t 2 = 0 and quickly increase due to the e 2 → e 1 popula- tion transfer.Note that, for both ESA signals, the QD(A + S) and QD(A) * An(S) spectra show pronounced oscillations that modulates the emission frequency along the t 2 time, that can be assigned to the relative displacement between the E and F states along the spectator mode (such modulation has indeed a period of ∼ 25 fs, consistent with the s 1 mode frequency of s 1 = 0.16 eV).Similar oscillations are also recorded in the e 2 SE signal, which is partially covered by the e 1 ESA, and can be seen more clearly in Figure s5 of the supporting information.
The spectra here reported further validates the QD(A) * An(S) approach as an inexpensive and exact (within the derived framework) alternative to the full QD(A + S) quantum dynamics.

Exploring the validity regime of the QD/analytical approach
In this section, we explore the robustness of the QD/analytical methods when the null commutation relation of Eq. ( 8) is not strictly respected: this means that the displacements of the E manifold PESs along the S modes are allowed to be slightly different.For simplicity, we restrict the analysis to two electronic states, e 1 and e 2 , in the E manifold, so that: therefore, to have a quasi-vanishing commutator, we require that By inserting the appropriate expressions for each of the three terms in Eq. ( 44), and developing the square of the harmonic PESs expressions, one gets: By introducing the compact notation Δd In the model system considered here, in which we only have one coupling mode ( m 1 ) and one spectator mode ( s 1 ), the above expression simplifies to The first part of the expression collects all the coupling coefficients, while the part in square brackets depends upon the explicit form of the WP (via the operators Q ).Therefore, one notes that there are three main factors that control the magnitude of this element: the strength of the inter-state coupling, the mode frequency, and the relative displacements between the wells.Moreover, the term in square brackets prescribes that the closer the WP centroid along the s 1 coordinate is to the average e 2 − e 1 displacement, the smaller this term will be.
Even if an exhaustive search of the parameter space is outside of the scope of the present paper, hereafter we assess the impact of neglecting this term in two cases: • A high-frequency spectator mode, with s 1 = 0.16 eV, de 2 e 1 ,s 1 = 1.25 (reorganization energy 0.125 eV, and Huang Rhys factor S = 0.78 ), i.e., the same values con- sidered in the model system above, and for which the relative displacement was varied with Δd steps of magni- tude 0.1 in the interval (−0.5, 0.5) (so that, for the largest displacements |Δd| = 0.5 , the well shifted toward the GS minimum has = 0.08 eV and S = 0.5 , while the other well has = 0.18 eV and S = 1.13).
• A low-frequency spectator mode, with s 1 = 0.05 eV, de 2 e 1 ,s 1 = 2.0 (reorganization energy = 0.1 eV, and Huang Rhys factor S = 2 ), and for which the relative displacement was varied with Δd steps of magnitude 0.35 in the interval (−1.75, 1.75) (so that, for the largest displacements |Δd| = 1.75 , the well shifted toward the GS minimum has = 0.032 eV and S = 0.63 , while the other well has = 0.207 eV and S = 4.13 ).The rationale with which we choose the de 2 e 1 ,s 1 value was again to have the QD(A) and QD(A + S) spectra significantly different.
In figure s6 of the supporting information we draw the PESs for Δd = 0 and for the largest Δd displacements for both the high-and low-frequency modes considered.We are not going to change the strength of the inter-state coupling e 1 e 2 ,m 1 and the average displacement de 2 e 1 ,s 1 (i.e., we will keep these quantities fixed, with e 1 e 2 ,m 1 = 0.15 eV as reported in Table 1).The accuracy of the QD/analytical expressions under the approximation of quasi-vanishing commutation are explored against a full QD approach by comparing both population dynamics and linear absorption spectra.In particular, in Figs. 6 and 7, we show the difference between the exact (full QD) population dynamics and LA spectra at various relative displacements Δd against the QD/analytical approach (which is exact only at Δd = 0 ), in the two cases of an high-fre- quency (Fig. 6) and a low-frequency (Fig. 7) spectator mode.Figure 8 reports a schematic representation of the change in the system PESs according to positive b) and negative c) Δd displacements.Note that positive and negative dis- placement will impact the spectra differently: if one focuses on the bright e 2 state, for positive displacements (and since de 2 e 1 ,s 1 is also positive), its PES will be pushed further away from the GS equilibrium position.Upon excitation from the GS, the spectrum will show a decrease in intensity of the low energy vibrational bands, and an increase in intensity of high-energy vibrational bands (for which the Franck-Condon overlap with the GS WP will increase for increasing displacement of the e 2 well).For negative displacements, instead, the e 2 PES will be pushed closer to the GS equilib- rium position, and the spectrum will change in the opposite direction (increase in low-frequency vibronic bands, and decrease in high-frequency ones).
The results of Figs. 6 and 7 shows that the QD/analytical approach is robust enough against significant relative displacements between the E manifold PESs along the spec- tator modes.Moreover, it is worth noticing that even if an error is introduced in the QD(A) * An(S) equations when compared to the full QD(A + S) approach, still the results is significantly improved with respect to the bare QD(A) approach that simply avoid considering the spectator modes altogether.
Here, we considered a simple case of one coupling and one spectator mode (Eq.47).We anticipate that it would be extremely interesting to analyze the more general case of Eq. ( 46) with multiple coupling and multiple spectator modes.In fact, even if such an analysis is outside of the scope of this exploratory work, an interesting feature can be noticed: by looking at the summation of contributions present in Eq. ( 46) one would guess that the more the modes, the highest the error (and therefore the worst the effect of not including such modes explicitly in the QD).Nonetheless, in the limit of many modes, along with the coupled electronic states can have a random distribution of both positive and negative relative displacements, it is possible that the net effect of such a multitude of disordered modes might quench this term.This can be the case of a solid in a disordered phase (for which one would additionally have a random distribution of energies that will further make the impact of the neglected commutator irrelevant).In this case, the exclusion of many potentially relevant modes might have a very small impact on the dynamics.

Discussion
In the following, we report some additional discussion on the general philosophy of the approach we propose and we briefly compare it with alternative approaches reported in the literature to reduce the number of modes in LVC models.When selecting a restricted number of modes that should enter in the QD simulation, one typically includes modes according to their (intra-or inter-state) coupling to the various electronic excited states of interest: since some spectator modes might be (similarly) strongly coupled to some electronic states, these would be included in the QD, and therefore contribute to the spectroscopy.Here, we showed an approach that can safely exclude them from the dynamics, making the QD computationally more efficient and/or allowing to include additional (dynamically relevant) modes in the model, while still accounting for the removed spectator modes in the spectra via inexpensive evaluation of analytical functions.The presented approach shows a practical strategy to let the QD focus only on the dynamically relevant modes, while still capturing the spectroscopic impact the left out modes exactly.
It should be mentioned that in the literature a rigorous and elegant approach has been proposed and exploited in several cases to reduce the dimensionality problem in the computation of linear absorption spectra of nonadiabatic systems driven by LVC Hamiltonians.Such scheme is based on a hierarchical transformation of the normal modes in effective modes, which are divided in blocks with the property that more blocks are included in the dynamics more momenta of the spectrum are exactly reproduced [16][17][18].This method has been also employed in combination with the introduction of spectral densities [19][20][21].Additionally, it has been generalized to quadratic vibronic coupling (QVC) Hamiltonians, [22] and adopted in combination with hybrid quantum/classical solutions of the dynamics [23,24].This approach revealed extremely effective for linear spectroscopy but, to the best of our knowledge, it has never been tested for nonlinear spectroscopy.As a matter of fact there are indications that the number of modes that need to be included to accurately reproduce the full system dynamics may rapidly increase with time, and for nonlinear spectroscopy, t 2 can be of sev- eral hundreds of femtoseconds [24].Moreover, the effectiveness of this approach decreases for a large number n of coupled states, since the first block of the hierarchy includes in principle n × (n + 1)∕2 modes).Our approach is less general and is not designed to reduce the number of couplings modes.Nonetheless, it is conceptually very simple, being grounded on an hypothesis that is simple to physically grasp and verify: the existence of displaced modes that do not remarkably dephase during the t 2 dynamics in the E manifold.Moreover, it sets the proper framework to understand to what extent one can straightforwardly generalize to nonadiabatic systems the adoption of vibronic analytical bandshape approaches, which are popularly employed in the computation of electronic spectra in systems with negligible inter-state couplings.It will be interesting in the future to perform some well-designed model study to compare the two approaches.

Conclusion
In this contribution, we have shown an effective way to simulate linear and nonlinear spectroscopy of a nonadiabatic system via a mixed quantum dynamics/analytical approach based upon partitioning of the molecular modes into two groups: a set of dynamically relevant modes, that should be included in the Hamiltonian that drives the (nonadiabatic) quantum dynamics; and a set of dynamically not relevant but spectroscopically relevant modes that do not need to be included in the dynamics but can still contribute to the spectral line shape, and that can be accounted for in an exact way via analytical line shape functions.We applied the protocol by means of the MCTDH quantum dynamics method, employing the package Quantics [10,11].The derived expressions are nonetheless completely general and whatever QD method (capable of producing the wave packet overlap) and software might be employed.
We have designed a simple model system, with two nonadiabatically coupled electronic states and three modes within a linear vibronic coupling framework: a coupling mode, a tuning mode (both included in the QD) and a spectator mode, characterized by a potential energy surface equally displaced for both the electronic states in the E manifold (i.e., the manifold of states explored by the quantum dynamics after the initial photoexcitation).We have shown that including or excluding the spectator mode in the QD does not change the population dynamics, while it significantly affects the linear absorption spectrum.We have also shown that instead of explicitly including it in the QD, it is possible to recover the same exact spectrum accounting for such a mode via analytical line shape functions.Such a result is obtained for both linear absorption and transient absorption spectra, with the latter being considerably more complex and expensive to simulate (which makes the reduced cost of the mixed QD/analytical approach even more attractive).
Even if the proposed QD/analytical approach is exact only when the nonadiabatically coupled electronic states are identically displaced along the spectator modes, we have shown that upon relaxation of such an hypothesis, the obtained approximate spectra are still very close to the exact spectra for a significant variation of the relative well displacements.Moreover, such spectra outperform the results obtained when such spectator modes are simply not accounted for.In general, one can foresee that this approach will be more accurate when the number of coupled states is small and, as it happens in many cases, show some similarities (e.g., those prepared through the excitation of an electron from the same occupied orbital) since, in this case, it is more likely that many modes shows similar displacements with respect to the ground state.
Before concluding, we highlight that a few noticeable and readily applicable extensions of the proposed approach can be envisaged: a) the possibility of considering temperature effects in bath modes was clearly demonstrated in the derivation of the line shape functions.This will allow to go beyond the simple Lorentzian or Gaussian broadening of the spectra, and consider a quasi-continuum of bath modes via arbitrarily specified spectral densities (as, e.g., the widely employed Drude-Lorentz spectral density).The only requirement (to have exact expressions) is again that of having all nonadiabatically coupled states identically coupled to the bath modes (even if the method proved to be robust under slight coupling differences that can therefore still be considered); b) orientational averaging of the spectra can be performed, by considering a small number of initial dipoles orientation (or field polarizations).The spectra would be obtained by averaging the results of the multiple QD propagation with different initial conditions (see ref. [14,25]; c) the derived equations are not limited to LVC Hamiltonians: they can be easily generalized to quadratic vibronic coupling (QVC) models, capable of including frequency changes and Duschinsky mixings also along the spectator modes [26].The requirement would be that of having a Duschinsky matrix separable in blocks (for active and spectator modes).In general, if a rigorous partitioning of the modes in active and spectator is doable, one could consider even more complex potential energy surface shapes and perform separate QD on the active and the spectator modes, considerably reducing the scaling of the calculation while being still able to get accurate spectra; d) exploring the extent of the error introduced in the limit of a large number of randomly displaced modes (as, e.g., in disordered solids); e) the present approach can be extended to the ML-MCTDH framework, in which even a larger number of active modes can be considered, making it applicable to very large molecular systems.

Fig. 1 a
Fig. 1 a Manifold of states and allowed dipole transitions; b PES displacements along a generic active mode m; c PES displacements along a generic spectator mode s (note that the E PESs are identically displaced) �e⟩⟨f � + h.c.

Fig. 4
Fig. 4 Comparison of a QD(A + S) (red curve) and QD(A) * An(S) (black dashed curve) linear absorption spectra with b the QD(A) spectrum (blue curve).In all the curves, a Lorentzian broadening was employed, with = 15 fs

e 2 e 1 ,0 2 −Fig. 5
Fig. 5 Comparison of transient absorption spectra for a full quantum approach QD(A + S) , b mixed QD/analytical approach QD(A) * An(S) , and c) QD approach that simply neglects the spectator modes QD(A).The (negative) ESA contributions are reported in

Fig. 6 a
Fig. 6 a, b (Bright state) Population dynamics and c, d linear absorption spectra for positive (left plots) and negative (right plots) Δd displacements, in the case of an high-frequency spectator mode with s = 0.16 eV, an average displacement between the two wells

Fig. 7 a
Fig. 7 a, b (Bright state) Population dynamics and c, d linear absorption spectra for positive (left plots) and negative (right plots) Δd displacements, in the case of an low-frequency spectator mode with

Fig. 8 and e 2 PESs: a zero relative displacement (both PESs are centered in de 2 e 1 along
Fig. 8 Schematic of the relative displacements between the e 1 and e 2 PESs: a zero relative displacement (both PESs are centered in de 2 e 1 along the mode s); b positive relative displacement (keeping the aver- brackets ⟨⋯⟩ g indicate the ensemble aver- age over the GS thermal distribution of vibrational states.It is now convenient to rewrite the operator that appears within the angular brackets of Eq. (21) in terms of the energy gap fluctuation between the g and the ē state.Let us thus define the energy gap operator δgē,S as Ĥē,S t ρg,S ⟩ = � e +i Ĥg,S t e −i Ĥē,S t � g where the angular

Table 1
Summary of the system parameters: modes frequency (in eV); intra-state coupling (in eV); reorganization energy (in eV); mode-dependent inter-state coupling in the E manifold ee ′ ,m (in eV); adiabatic energies E (ad) i and vertical energies E i (in eV) for the QD(A + S) and QD(A) models; and pair of dipole coupled states ij