# Quantum propagation using Hagedorn wave packets: generalised scheme

Rabiou Issa, Kokou Mawulonmi Robert Afansounoudji, Komi Sodoga, David Lauvergnat

► To cite this version:

Rabiou Issa, Kokou Mawulonmi Robert Afansounoudji, Komi Sodoga, David Lauvergnat. Quantum propagation using Hagedorn wave packets: generalised scheme. Molecular Physics, 2025, 8976, pp.e2554271. 10.1080/00268976.2025.2554271. hal-05338292

HAL Id: hal-05338292

https://cnrs.hal.science/hal-05338292v1

Submitted on 30 Oct 2025

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
# Quantum propagation using Hagedorn wave packets: generalized scheme

Rabiou Issa<sup>a</sup>, Kokou Mawulonmi Robert Afansounoudji<sup>a</sup>, Komi Sodoga<sup>a</sup>, and David Lauvergnat<sup>b</sup>

<sup>a</sup>Laboratoire de Physique des Matériaux et des Composants à Semi-conducteurs (LPMCS) and Centre d'Excellence Régionale pour la Maîtrise de l'Electricité (CERME), University of Lomé, 01BP1515 Lomé, Togo;

<sup>b</sup>Institut de Chimie Physique, Bâtiment 349, 15 avenue Jean Perrin, CNRS UMR8000, Université Paris-Saclay, 91405 Orsay Cedex, France

## ARTICLE HISTORY

Compiled August 27, 2025

## ABSTRACT

In this study, we present a novel wave packet propagation method that generalizes the Hagedorn approach by introducing alternative primitive basis sets that are better suited to describe various physical processes. We can mix time-dependent (Hagedorn) and time-independent basis sets, such as Fourier series, particle-in-a-box, and harmonic oscillator basis sets. Furthermore, our implementation can handle models with several electronic states, enabling the study of non-adiabatic processes. Instead of the time-dependent variational principle, our propagation scheme uses a three-step procedure (standard propagation, time-dependent parameter evaluation, and projection). This procedure relies on multidimensional integrations performed numerically with Gaussian quadrature, allowing us to avoid constraints on the form of the Hamiltonian operator. We have implemented this numerical algorithm in Fortran code, and validated it by comparing it with standard propagation schemes on harmonic and anharmonic 2D-models. Finally, we have applied our method to compute the vibrational spectrum of the 6D-modified Hénon-Heiles model and we show that our scheme reproduces well the results obtained with the standard approach. As a perspective, we show that with our generalized Hagedorn wave packet method, we are able to study the non-adiabatic dynamics of the cis-trans retinal photoisomerization with a reduced 2D-model.

## KEYWORDS

Quantum dynamics, Non-Adiabatic dynamics, Hagedorn wave packet, retinal, Time-Dependent Schrödinger Equation, Numerical Propagation

## 1. Introduction

Understanding the behavior of molecular systems requires quantum mechanics, which is essential for describing both their electronic structures and their nuclear motions. For example, the simulation of the dynamics of certain fundamental processes, such as

CONTACT R. Issa Email: issagoudouya@gmail.com

K. M. R. Afansounoudji Email: afankokourobert@gmail.com

K. Sodoga Email: ksodoga@univ-lome.tg

Corresponding author D. Lauvergnat Email: David.Lauvergnat@universite-paris-saclay.fr
multiphoton excitations [1, 2], photodissociation [3–5], photo-isomerization [6–10] and charge transfer [11–13] requires to solve numerically the time-dependent Schrödinger equation (TDSE). However, studying these processes can be very challenging (see references [14, 15]). There are several possible approaches to study these processes, such as, semiclassical[16–18], Mixed-Quantum Classical[19–23], Gaussian Wave Packets[24–27] and Wave Packet[10, 28–31] methods.

In the present study, we focus on quantum mechanical approaches to solve the TDSE using wave packet propagation techniques, so that restricting the studies to a temperature of 0 K. Additionally, the studies will be limited to systems of a finite-size, without infinite environments. Therefore, to solve the TDSE, we need to determine the time evolution of a wave packet, $Ψ$. Usually, the propagation is divided into small constant time intervals, $Δt$. Hence, the wave packet, $|Ψ(t + Δt)⟩$, at $t + Δt$ is obtained from the wave packet at $t$, $|Ψ(t)⟩$, as follows:

$$|Ψ(t + Δt)⟩ = \exp [-i/ħ \cdot \hat{H}Δt] |Ψ(t)⟩. \qquad (1)$$

where $\exp [-i/ħ \cdot \hat{H}Δt]$ is the evolution operator when the Hamiltonian is time-independent. Without being exhaustive, to solve the TDSE numerically, one can use standard methods [32]: (i) the split operator method [33], (ii) decomposition of the evolution operator into polynomials (usually Chebychev) [34], (iii) Taylor expansion of the evolution operator [35]. However, some of the above-mentioned methods[32–35] have drawbacks related to the number of nuclear degrees of freedom in the system. Indeed, when a direct-product is used, the size of the basis (or grid) scales exponentially with the number of degrees of freedom. In addition, the basis must accurately describe the evolution of the wave packet over time. Thus, if the wave packet delocalizes, the basis describing it should be very large. Similarly, if the center of the wave packet evolves strongly while remaining localized, the basis must also be very large. Furthermore, non-adiabatic processes present significant challenges to simulations. These transitions between electronic states often involve a topological feature known as a conical intersection. At this intersection, nuclear and electronic motions become intertwined by so-called non-adiabatic coupling, where quantum effects play a crucial role in the dynamics of the nuclear motion.

To solve the problem when the center of the wave packet evolves strongly, it is possible to represent the wave packet in a basis set (and grid) that evolves over time. More precisely, some basis set parameters change variationally during the propagation. Several wave packet propagation techniques have been developed, based on this idea of a time-evolving basis set (or grid). Without being exhaustive, one can mention: (i) The semiclassical approach introduced by Heller[36–39] in which a single Gaussian evolves over time. (ii) The approaches in which the wave packet is described by a sum of Gaussians evolving over time [24, 25, 27]. (iii) In the references [40, 41], the wave packet is represented on a grid that changes over time according to the evolution of the wave packet. (iv) A widely adopted solution in this context is the multi-configuration time-dependent Hartree (MCTDH) method [28, 42–44]. This approach represents the wavefunction as a small direct product of time-dependent basis sets. This method allows the propagation of quantum wave packets with a few dozen of degrees of freedom, or even larger with the multilayer version of MCTDH [29–31, 45]. However, this approach has its limitations. Both the kinetic and the potential components of the Hamiltonian must be expressed as sums of products of one-dimensional terms, which is not always
an easy task. (iv) George Hagedorn introduced the approach now known as Hagedorn wave packets[46–51], i.e. products of complex Gaussians with polynomials, also known as a prefactor polynomial [52–54], that form an orthonormal basis in $L^2$ space. Certain parameters (center, width and momentum) vary variationally during the propagation. Detailed presentations of these functions can be found in the literature [53, 55–61]. However, the standard Hagedorn scheme is not general. Indeed, it does not allow the use of other types of primitive basis sets, such as the particle-in-a-box or the Fourier basis sets, which are often used to model torsional motions in molecules. Furthermore, the standard Hagedorn scheme is restricted to a single electronic state.

This work aims to present a generalized Hagedorn wave packet propagation scheme that balances flexibility and numerical efficiency. Unlike the standard Hagedorn method, which relies on a specific basis set tied to a single electronic state, our approach: (i) combines time-dependent (Hagedorn) and time-independent basis sets (e.g., Fourier and particle-in-a-box) to better model diverse processes, such as torsional motions. (ii) handles multi-state non-adiabatic dynamics via a diabatic representation, critical for studying processes involving conical intersections. (iii) avoids variational constraints through a three-step algorithm: standard propagation, parameter update, and projection. This enables broader Hamiltonian compatibility.

More precisely, we provide a modern Fortran code (freely available on Github [62]) implementing this three-step propagation technique that generalizes the multidimensional Hagedorn basis set and works with other types of primitive basis sets.

The remainder of the paper is structured as follows. Section 2 provides a detailed presentation of the relevant theoretical methods. In Section 3, we present the results of tests conducted on 2D models with a single electronic surface. Section 4 focuses on applications, including the 6D Hénon-Heiles model. Finally, Section 5 summarizes the key points and discusses the results of this work. Furthermore, as perspective to non-adiabatic dynamics, the isomerization of retinal is considered as a 2D model with two electronic surfaces.

## 2. Methods

As mentioned in the introduction, the aim of the present study is to propose a propagation strategy with a time-dependent basis set, closely related to the Hagedorn wave packet propagation scheme. However, we will also compare this new approach against a standard approach with a time-independent basis set.

For a molecular system composed of $nc$ coordinates, $Q(Q = \{q_1, q_2, \dots, q_{nc}\})$, the total wave packet, $|\Psi(t)\rangle$, is expanded, as usual, on $ne$ electronic diabatic states, $|e\rangle$, as follows:

$$|\Psi(t)\rangle = \sum_{e=1}^{ne} |\psi^{(e)}(t)\rangle|e\rangle \qquad (2)$$

where, $|\psi^{(e)}(t)\rangle$, represents the nuclear wave packet, $\psi^{(e)}(Q, t)$, on the $e^{th}$ electronic diabatic state. This electronic component can be expanded on a multidimensional basis and written as:

$$\Psi^{(e)}(Q, t) = \sum_{I}^{NB} \tilde{C}_{I}^{t} \cdot B_{I}(Q, \Lambda), \qquad (3)$$
where $B_I(\mathbf{Q}, \mathbf{\Lambda})$ is the $I^{th}$ time-independent basis function and $\mathbf{\Lambda} = \{\lambda_1, \dots, \lambda_k, \dots, \lambda_{nc}\}$ are time-independent basis parameters (such as, the center and the scaling factor of the Harmonic Oscillator, HO, basis set). $\tilde{C}_I^t$ is the $I^{th}$ time-dependent coefficient, and $NB$ the size of the basis set.

In the present study, the standard propagation is performed with a simple but efficient Taylor expansion[35, 63] of the wave packet around $t$. The wave packet at $t + \Delta t$ is obtained as follows:

$$|\Psi(t + \Delta t)\rangle = \sum_{\ell=0}^{\ell_{max}} \frac{\Delta t^{\ell}}{\ell!} |\Psi^{(\ell)}(t)\rangle, \quad (4)$$

where in the equation 4, $|\Psi^{(\ell)}(t)\rangle$ is the $\ell^{th}$ derivatives of the wave packet at $t$. It can be obtained recursively from the time-dependent Schrödinger equation:

$$|\Psi^{(\ell+1)}(t)\rangle = -i/\hbar \cdot \hat{H} |\Psi^{(\ell)}(t)\rangle. \quad (5)$$

Remarks: $|\Psi^{(0)}(t)\rangle$ ($\ell = 0$) is simply the wave packet at time $t$. The previous relation (Eq. 5) is correct only when the Hamiltonian is time-independent. When it is not the case, time derivatives of the Hamiltonian must be taken into account[35, 63]. Furthermore, the order, $\ell_{max}$, of the Taylor series is obtained when the norm of the $\ell_{max}^{th}$ term, $\frac{\Delta t^{\ell_{max}}}{\ell_{max}!} |\Psi^{(\ell_{max})}(t)\rangle$, is smaller than a threshold (here $10^{-20}$). This guarantees that the wave packet is propagated with high accuracy, so that its norm and energy are numerically conserved over time.

## 2.1. Wave packet with time-dependent basis set

When the basis set is time-dependent, the wave packet expansion is similar to Eq. 3, except that the basis parameters, $\mathbf{\Lambda}_t$, are time-dependent. Then this expansion becomes:

$$\psi^{(e)}(\mathbf{Q}, t) = \sum_{I=1}^{NB} C_I^t \cdot B_I(\mathbf{Q}, \mathbf{\Lambda}_t). \quad (6)$$

In the previous equation, (Eq. (6)), $B_I(\mathbf{Q}, \mathbf{\Lambda}_t)$ is the $I^{th}$ time-dependent multidimensional basis function, the $C_I^t$ is its corresponding coefficient. The $\mathbf{\Lambda}_t$ are the time-dependent parameters of the basis set common to all electronic states, where $\mathbf{\Lambda}_t = \{\lambda_1^{(t)}, \dots, \lambda_k^{(t)}, \dots, \lambda_{nc}^{(t)}\}$. Each basis function is expressed as a product of 1D-primitive basis functions, $b_{i_k}^k(q_k, \lambda_k^{(t)})$, as follows:

$$B_I(\mathbf{Q}, \mathbf{\Lambda}_t) = \prod_{k=1}^{nc} b_{i_k}^k(q_k, \lambda_k^{(t)}). \quad (7)$$

In the previous equation (Eq. 7), $I$ is corresponding to a multidimensional index, $\{i_1, \dots, i_k, \dots, i_{nc}\}$, associated with the 1D-primitive basis functions. The total number of functions, $NB$, is expressed as $\prod_{k=1}^{nc} nb_k$ for a direct-product basis set, where $nb_k$ is the number of functions of the $k^{th}$ primitive basis set. The $\lambda_k^{(t)}$ are the time-dependent parameters of the $k^{th}$ primitive basis set and are particularized according to its type: (i)
for the Fourier basis set there are no time-dependent parameters (ii) for the sine or the particle-in-box basis set, the parameters (the domain range), $\lambda_k = \{A_k, B_k\}$, are time-independent. Although, they can also be implemented as time-dependent parameters. (iii) for the Hagedorn primitive basis set, the time-dependent parameters are particu-larized in $\lambda_k^{(t)} = \{\alpha_k^{(t)}, q_k^{(t)}, p_k^{(t)}\}$ where $\alpha_k^{(t)}$ is related to the basis set width and $q_k^{(t)}, p_k^{(t)}$ represent the center, and the momentum of the Hagedorn basis set, respectively. More precisely, the expression of the Hagedorn primitive basis function, $b_{i_k}^k(q_k, \lambda_k^{(t)})$ (or more concisely $b_{i_k}^k(q_k, t)$), is given by:

$$ \begin{aligned} b_{i_k}^k (q_k, t) &\equiv b_{i_k}^k (q_k, \lambda_k^{(t)}) \\ &= N (\lambda_k^{(t)}) \cdot H_{i_k-1} (x_k) \times \exp \left( -\frac{\alpha_k^{(t)}}{2} \Delta q_k^2 + i \cdot p_k^{(t)} \Delta q_k \right) \end{aligned} \quad (8) $$

where $\alpha_k^{(t)}$ is a complex number ($\alpha_k^{(t)} = a_k^{(t)} + i \cdot b_k^{(t)}$) where $a_k^{(t)}$ and $b_k^{(t)}$ are its real and imaginary parts, respectively. In Eq (8), $\Delta q_k = q_k - q_k^{(t)}$ and $x_k = \sqrt{a_k^{(t)}} \Delta q_k$. The square root of $a_k^{(t)}$ is a scaling factor and it is related to the width of the Gaussian part of the basis set. Therefore, $a_k^{(t)}$ is positive. Furthermore, $N(\lambda_k^{(t)})$ is a normalization parameter that depends only on $a_k^{(t)}$ and $H_{i_k}(x_k)$ is the usual normalized Hermite polynomial of degree $i_k$. It is worth noting that: (i) The Hagedorn and Harmonic Oscillator (HO) basis sets are similar, except that the imaginary contributions of the Hagedorn basis set are absent in the HO ones. (ii) Some parameters can be excluded for the time-dependent basis. For instance, the imaginary part $b_k^{(t)}$ of $\alpha_k^{(t)}$ or $p_k^{(t)}$. This does not prevent an exact propagation, since the missing contributions are transferred to the $C_I^t$ coefficients if the basis set is large enough. (iii) The Heller approach [36, 64, 65] can be cast in the Hagedorn scheme where only one basis set is present ($nb_k = 1$ and hence, $NB = 1$). (iv) In the standard Hagedorn or Heller schemes, the time-dependent parameter $\alpha^{(t)}$ is a full complex $nc \times nc$ matrix. In our scheme, this matrix is diagonal and the diagonal elements are the $\alpha_k^{(t)}$. This approximation does not prevent an exact propagation (the missing coupling contributions are transferred to the $C_I^t$ coefficients) and was made to facilitate the standard grid-to-base transformation required in our scheme. (v) With our scheme, we can mix time-dependent and time-independent basis sets, and it also works with several electronic states.

## 2.2. Propagation strategy

The key point of our Hagedorn approach lies in the way we propagate the time-dependent parameters $\{C^t, \Lambda_t\}$, where $\Lambda_t$ are the basis set parameters and $C^t = \{C_1^t, \dots, C_{NB}^t\}$ are the wave packet coefficients. Indeed, the usual way to propagate the Hagedorn wave packet or other approaches (Heller [37, 38, 66], Variational Multi-Configuration Gaussian[67–76], MCTDH [77, 78]) from time $t_1$ to time $t_2$ is to use the time-dependent variational principle. However, we use a different strategy with three steps (fig. 1):

(1) Standard propagation: In this step (See fig. 1), we perform a standard propagation technique (here, the Taylor propagation scheme, but other schemes could be used[32]) assuming a time-independent basis, $B(\Lambda_{t_1})$, to obtain an intermediate
wave packet at $t_2$:

$$ \Psi(Q, t_2) = \sum_{I=1}^{NB} \tilde{C}_I^{t_2} \cdot B_I(Q, \Lambda_{t_1}). \qquad (9) $$

We obtain a new set of $\{\tilde{C}^{t_2}, \Lambda_{t_1}\}$. It is important to note that $\Lambda_{t_1}$ is still at time $t_1$ and the wave packet coefficients, $\tilde{C}^{t_2}$, are not the final ones. Furthermore, and depending on the propagation technique used, the time step ($\Delta t = t_2 - t_1$) can be relatively large.

(2) New basis set: Then we build the new basis set, $\mathbf{B}(\Lambda_{t_2})$, by computing the parameters $\Lambda_{t_2}$, at time $t_2$ (See fig. 1). More precisely, the parameters, $q_k^{(t)}$, $p_k^{(t)}$ and $\alpha_k^{(t)}$ (for $k = 1, \dots, nc$), are computed from the expectation values of operators assuming the wave packet has a Gaussian shape (Heller-like wave packet):

$$ q_k^{(t)} = \langle \Psi(t) | q_k | \Psi(t) \rangle, \qquad (10) $$

$$ p_k^{(t)} = -i \langle \Psi(t) | \partial q_k | \Psi(t) \rangle \qquad (11) $$

$$ a_k^{(t)} = \left[ 2 \left( \langle \Psi(t) | q_k^2 | \Psi(t) \rangle - (q_k^{(t)})^2 \right) \right]^{-1}, \qquad (12) $$

$$ b_k^{(t)} = a_k^{(t)} \left[ 2 p_k^{(t)} q_k^{(t)} + i \left( 1 + 2 \left\langle \Psi(t) \middle| q_k \frac{\partial}{\partial q_k} \middle| \Psi(t) \right\rangle \right) \right]. \qquad (13) $$

In the previous equations, $q_k^{(t)}$ and $p_k^{(t)}$ are, respectively, the expectation values of the position operator, $q_k$, (Eq. 10) and momentum operator, $\hat{p}_k$, (Eq. 11). As mentioned above, the parameter $a_k^{(t)}$ (Eq. 12) is related to the width or the variance of the wave packet and can be computed from the expectation values of the $q_k^2$ and $q_k$ operators. There are several ways to compute $b_k^{(t)}$, but some of them are not numerically stable. Therefore, in the present implementation (Eq. 13), we used the expectation value of the $q_k \frac{\partial}{\partial q_k}$ operator and the computed expressions of $q_k^{(t)}$ and $p_k^{(t)}$. It is worth mentioning, these time-dependent parameters are not optimal (except for the Heller wave packet), but this does not prevent accurate propagation.

(3) Projection: After, the standard propagation (the first step), the wave packet is known at time $t_2$ and as a projection onto the basis set at time $t_1$ (the old basis set, $\mathbf{B}(\Lambda_{t_1})$). The wave packet expansion is characterized by the coefficients, $\tilde{C}^{t_2}$, and the time-dependent basis parameters, $\Lambda_{t_1}$ (Eq. 9). However, this wave packet must be projected onto the newly constructed basis set $\mathbf{B}(\Lambda_{t_2})$ at time $t_2$ obtained in the second step:

$$ \Psi(Q, t_2) = \sum_{I=1}^{NB} C_I^{t_2} \cdot B_I(Q, \Lambda_{t_2}) \qquad (14) $$

Then, this wave packet expansion is characterized by the coefficients, $\tilde{C}^{t_2}$, and the time-dependent basis parameters, $\Lambda_{t_2}$. It is then assumed that the two expansions,
Eq. 9 and Eq. 14, are equal:

$$ \sum_{J=1}^{NB} C_J^{t_2} \cdot B_J(\mathbf{Q}, \mathbf{\Lambda}_{t_2}) = \sum_{I=1}^{NB} \tilde{C}_I^{t_2} \cdot B_I(\mathbf{Q}, \mathbf{\Lambda}_{t_1}). \quad (15) $$

Then, to obtain the $J^{th}$ new coefficient, $C_J^{t_2}$, the previous equation (Eq 15) is projected onto the new basis function, $B_J(\mathbf{Q}, \mathbf{\Lambda}_{t_2})$:

$$ C_J^{t_2} = \sum_{I=1}^{NB} \langle B_J(\mathbf{Q}, \mathbf{\Lambda}_{t_2}) | B_I(\mathbf{Q}, \mathbf{\Lambda}_{t_1}) \rangle \tilde{C}_I^{t_2}. \quad (16) $$

It is important to note that, the equality assumption (Eq. 15) only holds true when the basis set size is infinite in the second expansion (Eq. 14). However when the basis set size is equal to NB, the numerical values of the new coefficients, $C_J^{t_2}$, obtained from the projection (Eq. 16), are still numerically exact. Furthermore, if the number of terms in the sum is be too small (Eq. 14), the wave packet norm may not be conserved (see solutions below to overcome this difficulty). This equation (Eq. 16) shows that the projection is performed on a multidimensional basis. However, numerically, it can be achieved sequentially using partial projections (for one 1D-basis set, see the Appendix A for more details). It is important to note that this sequential transformation is possible because $\alpha_k^{(t_2)}$ is considered to be a diagonal matrix. When, it is not the case, one can use an algebraic method as proposed by Vaníček et al. [53].

Some comments are in order: in the first step, the wave packet norm and energy can be conserved with some high numerical accuracy, and with some propagation techniques (such as Lanczos) they are rigorously conserved.[32] However, during the third step, the wave packet norm may not be conserved if the basis set size, NB, is not large enough or if the time step, $\Delta t$, is too large. Several strategies can be used to overcome this difficulty: (i) increase the basis set size, (ii) reduce the time step, (iii) renormalize the wave packet with the norm after the first step. Furthermore, in our implementation, we want to test several options regarding the number of time-dependent basis parameters. More precisely, $b_k^{(t)}$ (the imaginary part of $\alpha_k^{(t)}$) and the term with $p_k^{(t)}$ can be present or absent in the time-dependent basis parameters. For the figures and the tables, the following notations are used:

* For the Hagedorn propagation scheme, the keyword HAG is used, then the number of 1D-basis function is added (assuming this number is the same for all dimensions), then two times a letter to specify when $b_k^{(t)}$ and $p_k^{(t)}$ are present, T (for true) or not F (for false), respectively. Finally another letter is added to indicate whether the re-normalization at the end of step 3 is used (T) or not (F). For instance, HAG5TTF means: Hagedorn propagation scheme with 5 basis functions for each dimension and keeping the time dependent parameters, $b_k^{(t)}$ and $p_k^{(t)}$. In addition, the re-normalization is not used.

* For propagation with a time-independent basis set, the keyword, STD is used and then the number of 1D-basis functions is added. STD10 means: Standard propagation scheme with 10 basis functions for each dimension.
<figure><img src="image_2.png" /><figcaption>Figure 1.: Propagation algorithm from $t_1$ to $t_2$ in three steps. Step 1: Standard propagation with a time-independent basis set (parameters, $\Lambda_{t_1}$ at $t_1$). Step 2: Evaluation of the basis parameters, $\Lambda_{t_2}$, and basis set construction at $t_2$. Step 3: Projection of the wave packet at $t_2$ expanded on the basis set in $t_1$ into the basis set in $t_2$. re-normalization if required.</figcaption></figure>

### 3. Computational tests

In this section, the accuracy of our approach is tested on two 2D-models (Eq. 17) by comparing a wave packet obtained by the Hagedorn propagation scheme with another one converged with high accuracy obtained by a standard propagation method. More specifically, the three steps of our approach are tested. For the 2D-model ($nc = 2$), the Hamiltonian model is given by:

$$ \hat{H}(\mathbf{Q}) = -\frac{\hbar^2}{2m} \sum_{k=1}^{nc} \frac{\partial^2}{\partial q_k^2} + V(\mathbf{Q}) \quad (17) $$

where $m$ is the mass (here $m = 1$) and $V(\mathbf{Q})$ is the potential. The initial wave packet is assumed to be a Gaussian and is defined as follows:

$$ \psi_0^{(e_0)}(\mathbf{Q}) = \prod_{k=1}^{nc} \exp \left( -\frac{a_k^0}{2} \Delta_k^2 + i \cdot p_k^0 \Delta_k \right), \quad \Delta_k = q_k - q_k^0, \quad (18) $$

where, $e_0$ is the initial electronic state (in this model, $e_0 = 1$), $p_k^0$ and $a_k^0$ are, respectively, the initial momentum and a parameter related to the width of the Gaussian along the $k^{th}$ coordinates ($q_k$).
### 3.1. 2D-Harmonic Hamiltonian

In this section, we will use a simple 2D-quadratic potential defined as follows:

$$V(\mathbf{Q}) = \frac{1}{2} (q_1^2 + q_2^2). \qquad (19)$$

With this 2D-harmonic model and a Gaussian initial wave packet (Eq.18), the wave packet must remain Gaussian during propagation, which is equivalent to the Heller propagation scheme [36, 65]. This model enables us to check the influence of the time step and the basis set size, $nb_k$ on the propagation accuracy. For this model, the Hagedorn basis set parameters and the initial wave packet parameters are given in the table (1) and have been chosen so that the parameters associated with the first coordinate ($k = 1$) evolve during the propagation, while those associated with the second coordinate ($k = 2$) remain constant during the propagation. Furthermore, the number of grid points is determined by adding 5 to the number of basis functions, such that $nq_k = nb_k + 5$. This choice ensures high accuracy during numerical integration using Gaussian quadrature.

<table><caption>Table 1.: Hagedorn and HO basis set parameter values: $q_k^{(0)}$, $p_k^{(0)}$, and $\alpha_k^{(0)}$ (columns 2-4). Initial Gaussian wave packet parameters (Eq. 18, with $e_0 = 1$): $q_k^{(0)}$, $p_k^{(0)}$, and $a_k^{(0)}$ (columns 5-7)</caption><thead><tr><th rowspan="2">k</th><th colspan="3">Hagedorn and HO basis sets</th><th colspan="3">Initial wave packet</th></tr><tr><th>q<sub>k</sub><sup>(0)</sup></th><th>p<sub>k</sub><sup>(0)</sup></th><th>α<sub>k</sub><sup>(0)</sup></th><th>q<sub>k</sub><sup>0</sup></th><th>p<sub>k</sub><sup>0</sup></th><th>a<sub>k</sub><sup>0</sup></th></tr></thead><tbody><tr><td>1</td><td>2.0</td><td>0.0</td><td>(1.2, 0.0)</td><td>2.0</td><td>0.0</td><td>1.2</td></tr><tr><td>2</td><td>0.0</td><td>0.0</td><td>(1.0, 0.0)</td><td>0.0</td><td>0.0</td><td>1.0</td></tr></tbody></table>

First, we perform the overall analysis of the wave packet along time to check the numerical stability of our algorithm. As mentioned above, and since the wave packet remains Gaussian along time with this harmonic model, the wave packet must have a basis expansion where only the first coefficient, $C_1^t$, is non-zero and the other coefficients are zero ($C_I^t = 0, (I > 1)$). Since the wave packet is normalized, $|C_1^t| = 1$. Therefore, the residual coefficient error $R_C(t) = \sum_{I=2}^{NB} |C_I^t|^2$ and the norm conservation error, $\Delta N_3(t) = |N_3(0) - N_3(t)|$, allow to measure the wave packet deviation from its Gaussian shape. $N_3(t)$ denotes the norm of the wave packet at time $t$ after the three-step propagation scheme. The figure (2) shows the time evolution of the norm conservation error, $\Delta N_3(t)$, as a function of the number of basis functions, $nb_k$, in each dimension with a fixed time step of 0.25 au (left panel) and as a function of the time step, $\Delta t$ with $nb_k = 5$. As expected, the norm conservation error becomes smaller as $nb_k$ increases. For $nb_k = 10$ or larger, the error is around $10^{-9}$. Furthermore, with only five basis functions in each dimension, the norm error tends to zero as the time step decreases. For $\Delta t = 0.01$ au, the error is approximately $4.10^{-8}$. The residual coefficient error, $R_C(t)$, shows similar trends as $nb_k$ increases or as $\Delta t$ decreases (see Fig. 3). The fact that this error tends to zero shows that the wave packet remains Gaussian with time, since only the first coefficient, $C_1^t$, is non-zero. Furthermore, its absolute value must remain one, since the wave packet norm is one (Fig. 2). The energy conservation error (not shown),
<figure><img src="image_3.png" /><figcaption>Figure 2.: Time evolution of the wave packet norm deviation $\Delta N_3(t)$. The left panel shows $\Delta N_3(t)$ as a function of the basis size $nb_k$, for a fixed time step $\Delta t = 0.25$ au. The right panel displays $\Delta N_3(t)$ as a function of the time step $\Delta t$, with a fixed basis size $nb_k = 5$. The quantity $\Delta N_3(t) = |N_3(t) - N_3(t_0)|$ measures the deviation between the instantaneous norm $N_3(t)$ and the initial wave packet norm $N_3(t_0)$.</figcaption></figure>

$\Delta E_3(t) = |\langle\Psi(0)|\hat{H}|\Psi(0)\rangle - \langle\Psi(0)|\hat{H}|\Psi(t)\rangle|$, also presents a similar trend and with $\Delta t = 0.01$ au and $nb_k = 5$, this error is about $1.4 \times 10^{-12}$ au.

<figure><img src="image_4.png" /><figcaption>Figure 3.: Time evolution of the wave function residual coefficient $R_C(t)$. The left panel presents $R_C(t)$ as a function of the basis size $nb_k$ for a fixed time step $\Delta t = 0.25$ au, while the right panel shows $R_C(t)$ as a function of the time step $\Delta t$ with a fixed basis size $nb_k = 5$. The residual coefficient $R_C(t)$ quantifies the wave function deviation from the Gaussian shape at time $t$.</figcaption></figure>

Let us now analyze the three steps of our numerical scheme:

(i) The first step corresponds to a standard propagation using a wave packet Taylor expansion, for which the Taylor order is adjusted until the expansion converges (the norm of last term must be smaller than $10^{-20}$). With $\Delta t = 0.25$ au and $nb_k = 10$, the order is about 30 and the wave packet norm ($N_1$) and energy are conserved with accuracies of $10^{-16}$ and $10^{-10}$ au, respectively. This guarantees that the numerical error of this step is practically zero. Although the size of the
basis, NB, must be large enough to obtain an accurate wave packet.

(ii) In the second step, the time-dependent basis parameters $\Lambda_t$ are computed by numerical integration with Gauss-Hermite quadrature grids, and since the number of grid points is larger than the number of basis functions, these parameters can be considered as numerically exact. For our test, we chose $nq_k = nb_k+5$, although a smaller value should be sufficient. The time evolution of these parameters (with $\Delta t = 0.25$ au and $nb_k = 10$) is shown in the following figures (Figs. 4 and 5) and they show a smooth evolution over time. Furthermore, as expected, all parameters remain constant for coordinates 2 ($k=2$) along the propagation.

<figure><img src="image_5.png" /><figcaption>Figure 4.: Time evolution of the Gaussian parameters $q_k^{(t)}$ (left) and $p_k^{(t)}$ (right) in the $k^{th}$ primitive basis ($\Delta t = 0.25$ au, $nb_k = 10$). The expectation values $q_k^{(t)} \equiv \langle \Psi(t) | q_k | \Psi(t) \rangle$ (coordinate) and $p_k^{(t)} \equiv -i \langle \Psi(t) | \frac{\partial}{\partial q_k} | \Psi(t) \rangle$ (momentum) are shown for the first ($k=1$, purple) and second ($k=2$, green) degrees of freedom.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5.: Time evolution of the complex Gaussian parameter $\alpha_k^{(t)} = a_k^{(t)} + ib_k^{(t)}$ components in the $k^{th}$ primitive basis. The left panel displays the real part $a_k^{(t)}$ while the right panel shows the imaginary part $b_k^{(t)}$, computed with a time step $\Delta t = 0.25$ au and a basis size $nb_k = 10$. The purple and green curves represent the first ($k=1$) and second ($k=2$) degrees of freedom, respectively.</figcaption></figure>
(iii) The final step constitutes the most critical part of the procedure. The wave packet $\Psi(Q, t_2)$ (Eq. 9), computed at time $t_2 = t_1 + \Delta t$ within the basis set parametrized by $\Lambda_{t_1}$ (defined at $t_1$), must be subsequently projected onto the time-evolved basis set determined by $\Lambda_{t_2}$ (evaluated at $t_2$). The accuracy of this projection depends on the size of the basis set, $NB$ (or the $nb_k$) and the time step, $\Delta t$. The figure(Fig. 6) shows the norm conservation, $\Delta N_{1,3}(t) = |N_1(t) - N_3(t)|$ where $N_1(t)$ and $N_3(t)$ denote the norms of the wave packet at time $t$ after the first and third steps of the three-step propagation scheme, respectively. In fact, for a given time step (here $\Delta t = 0.25$ au), the projection is only exact if the basis set is large enough. Therefore, one can see that the error decreases to zero as $nb_k$ increases (left panel). Furthermore, for a given basis set size ($nb_k = 5$), the error decreases as $\Delta t$ decreases. Indeed, as $\Delta t$ decreases, the wave packet modification between each time step becomes less important, so the new basis set at $t_2$ will be closer to the one at $t_1$. Therefore, the projection will be more accurate.

<figure><img src="image_7.png" /><figcaption>Figure 6.: Time evolution of the wave packet norm deviation $\Delta N_{1,3}(t) = |N_1(t) - N_3(t)|$ in the three-step propagation scheme. The left panel shows the dependence on basis size $nb_k$ at fixed time step $\Delta t = 0.25$ au, while the right panel displays the dependence on time step $\Delta t$ for fixed basis size $nb_k = 5$. $N_1(t)$ and $N_3(t)$ denote the norms of the wave packet at time $t$ after the first and third steps of the three-step propagation scheme, respectively.</figcaption></figure>

## 3.2. 2D-anharmonic Hamiltonian

In this section, we use a modified Hénon-Heiles potential with two degrees of freedom. Indeed, in the original Hénon-Heiles potential, the presence of cubic terms leads to holes in the potential. Therefore, to avoid the wave packet delocalization around these holes, we introduced a modified version of the potential (Eq. 20), in which some $q_k$'s are replaced by $\tanh(q_k)$. This modification avoids the holes without changing the first, second and third order derivatives around the main minimum.

$$V(Q) = \frac{mw^2}{2} \sum_{k=1}^{nc} q_k^2 + \lambda \sum_{k=1}^{nc-1} \left( q_k^2 \cdot \tanh(q_{k+1}) - \frac{1}{3} \tanh(q_k)^3 \right) \quad (20)$$

The initial wave packet is a Gaussian shifted along the first coordinate (the parameters are given in Table 2) and the propagation time is $\Delta T = 60$ au with a time step
of $dt = 0.1 \text{ au}$. The parameters of the potential $V(Q)$ of the system are $m = 1.0$, $\omega = 1.0$ and $\lambda = 0.111803$ [79]. For the standard techniques, both coordinates are described using harmonic oscillator (HO) basis sets, while for the Hagedorn scheme, both coordinates are described using Hagedorn basis sets. The parameters of the basis sets for both techniques are summarized in Table 2. For the calculations, we use the number of functions in each primitive basis as $nb_k = 5, 20, 70$ and for the number of grid points, we add 5 grid points to the number of basis functions so that $nq_k = nb_k + 5$, i. e., $nq_k = 10, 25, 75$. This choice of using more grid points than basis functions is made to ensure high accuracy of numerical integration.

<table><caption>Table 2.: Harmonic Oscillator (HO) and Hagedorn basis set parameter values: $q_k^{(0)}$, $p_k^{(0)}$, and $\alpha_k^{(0)}$ (columns 2–6). Initial Gaussian wave packet parameters (Eq. 18, with $e_0 = 1$): $q_k^0$, $p_k^0$, and $a_k^0$ (columns 7–9).</caption><thead><tr><th rowspan="2">k</th><th colspan="2">HO basis set</th><th colspan="3">Hagedorn basis set</th><th colspan="3">Initial wave packet</th></tr><tr><th>q<sub>k</sub><sup>(0)</sup></th><th>a<sub>k</sub><sup>(0)</sup></th><th>q<sub>k</sub><sup>(0)</sup></th><th>p<sub>k</sub><sup>(0)</sup></th><th>α<sub>k</sub><sup>(0)</sup></th><th>q<sub>k</sub><sup>0</sup></th><th>p<sub>k</sub><sup>0</sup></th><th>a<sub>k</sub><sup>0</sup></th></tr></thead><tbody><tr><td colspan="9" style="text-align:center;">2D</td></tr><tr><td>1</td><td>2.0</td><td>1.2</td><td>2.0</td><td>0.0</td><td>(1.2, 0.0)</td><td>2.0</td><td>0.0</td><td>1.2</td></tr><tr><td>2</td><td>0.0</td><td>1.0</td><td>0.0</td><td>0.0</td><td>(1.0, 0.0)</td><td>0.0</td><td>0.0</td><td>1.0</td></tr></tbody></table>

With this model, the wave packet does not remain Gaussian, therefore, the accuracy of our procedure is tested with the norm of the difference of two wave packets, $N_{\Delta\Psi}(t)$, namely, the propagation error. The first one, $|\Psi^{\text{std}}(t)\rangle$, is obtained using a standard propagation technique (here Taylor expansion), and the second one, $|\Psi^{\text{ProjHag}}(t)\rangle$, is obtained after projecting the Hagedorn wave packet onto the basis set used for the standard propagation. This projection is similar to the one in third step of the Hagedorn procedure. More precisely, to obtain $N_{\Delta\Psi}(t)$, we used the coefficients associated with these two wave packets, $\tilde{C}^t$ (standard propagation) and $C^t_{\text{ProjHag}}$ (projected Hagedorn propagation) and $N_{\Delta\Psi}(t) = \|\tilde{C}^t - C^t_{\text{ProjHag}}\|$. Before checking the accuracy of the wave packet obtained with the Hagedorn scheme, we need to check that the standard wave packet is converged with respect to the basis set size. Our tests show that with a basis set with 70 basis functions ($nb_k = 70$) for each dimension, it is possible to converge the wave packet with an error (with an expression similar to $N_{\Delta\Psi}(t)$) smaller than $10^{-9}$ with respect to a standard propagation with $nb_k = 90$. With this converged wave packet obtained with the standard propagation scheme, we are able to check the accuracy of the Hagedorn scheme (with $\Delta t = 0.1 \text{ au}$) using several options, but without re-normalization after Hagedorn step 3. In figure 7, we present the propagation error ($N_{\Delta\Psi}(t)$) along time for different options (i) Full Hagedorn scheme, i. e. including the variation of $b_k^{(t)}$ (the imaginary part of $\alpha_k^{(t)}$) and the one of $p_k^{(t)}$. This option is labeled (TTF). (ii) A scheme, labeled (FTF), without $b_k^{(t)}$ and with the variation of $p_k^{(t)}$. (iii) A scheme without $b_k^{(t)}$ and without $p_k^{(t)}$, labeled (FFF). These three options are shown in the figure 7 and one can see that for a large basis set, $nb_k \ge 70$, the propagation errors become small for all options ($< 5 \cdot 10^{-6}$ for the TTF and FTF options and $< 10^{-7}$ for the FFF option). This shows the accuracy of our schemes. Furthermore, with this large basis set ($nb_k = 70$) and for all options, the projection error in step 3 (not shown), $\Delta N_{1,3}(t)$, is also very small ($< 10^{-11}$). Therefore, the re-normalization after step 3 is not necessary for large
basis sets. For a small value of nb (nb = 5) and with or without re-normalization after

<figure><img src="image_8.png" /><figcaption>Figure 7.: Time evolution of $N_{\Delta\Psi}(t)$ as a function of nb with $\Delta t = 0.1$ au. (i) Left panel: with $b_k^{(t)}$ and $p_k^{(t)}$, (TTF option). (ii) Middle panel: without $b_k^{(t)}$ and with $p_k^{(t)}$, (FTF option). (iii) Right panel: without $b_k^{(t)}$ and without $p_k^{(t)}$, (FFF option).</figcaption></figure>

the projection, several points need to be discussed: (i) It is important to note that the Hagedorn scheme when considering the time variation of $p_k^{(t)}$ and with or without the variation of $b_k^{(t)}$ (options: (TTF) or (FTF)) is able to propagate up to 60 au. The other Hagedorn scheme option (FFF) stops after about 1. au. Indeed, with this option, the wave packet becomes unphysical and the wave packet norm and the total energy are not conserved. The energy conservation error rises by 2.7 Hartree after a propagation of 1. au (see Fig. (8)). Therefore, in the first step of the Hagedorn scheme, the standard propagation with the Taylor series cannot converge (reducing the time step by a factor of 10 does not solve the problem). This problem is due to the projection (step 3 of the Hagedorn scheme), which is less accurate when $p_k^{(t)}$ is not included in the Hagedorn scheme. In other words, more basis functions are needed to reach a given accuracy. It is worth noting, that for the (FTF) and (TTF) options, the total energy is conserved up to $2 \times 10^{-1}$ and $7 \times 10^{-2}$ Hartree, respectively. (ii) For the (TTF) or (FTF) options, the propagation error can be very large up to one (Fig. 7), which clearly indicates that the basis set is not large enough without preventing to perform a relatively long propagation. However, we need to check if this will be problematic to compute a spectrum.

### 3.3. Comparative Analysis of Energy Spectrum

In this section, we also use the modified Hénon-Heiles potential with two degrees of freedom. The basis and initial wave packet parameters are given in the Table 3.
<figure><img src="image_9.png" /><figcaption>Figure 8.: Energy as a function of time for a propagation with nb<sub>k</sub> = 5: (i) with b<sub>k</sub><sup>(t)</sup> and p<sub>k</sub><sup>(t)</sup>, (TTF option). (ii) without b<sub>k</sub><sup>(t)</sup> and with p<sub>k</sub><sup>(t)</sup>, (FTF option). (iii) without b<sub>k</sub><sup>(t)</sup> and without p<sub>k</sub><sup>(t)</sup>, (FFF option).</figcaption></figure>

<table><caption>Table 3.: Harmonic Oscillator (HO) and Hagedorn basis set parameter values: q<sub>k</sub><sup>(0)</sup>, p<sub>k</sub><sup>(0)</sup>, and α<sub>k</sub><sup>(0)</sup> (columns 2–6). Initial Gaussian wave packet parameters (Eq. 18, with e<sub>0</sub> = 1): q<sub>k</sub><sup>0</sup>, p<sub>k</sub><sup>0</sup>, and a<sub>k</sub><sup>0</sup> (columns 7–9)</caption><thead><tr><th rowspan="2">k</th><th colspan="2">HO basis set</th><th colspan="3">Hagedorn basis set</th><th colspan="3">Initial wave packet</th></tr><tr><th>q<sub>k</sub><sup>(0)</sup></th><th>a<sub>k</sub><sup>(0)</sup></th><th>q<sub>k</sub><sup>(0)</sup></th><th>p<sub>k</sub><sup>(0)</sup></th><th>α<sub>k</sub><sup>(0)</sup></th><th>q<sub>k</sub><sup>0</sup></th><th>p<sub>k</sub><sup>0</sup></th><th>a<sub>k</sub><sup>0</sup></th></tr></thead><tbody><tr><td>1</td><td>0.0</td><td>1.0</td><td>2.0</td><td>0.0</td><td>(1.2, 0.0)</td><td>2.0</td><td>0.0</td><td>1.2</td></tr><tr><td>2</td><td>0.0</td><td>1.0</td><td>0.0</td><td>0.0</td><td>(1.0, 0.0)</td><td>0.0</td><td>0.0</td><td>1.0</td></tr></tbody></table>

However, the aim is the computation of the spectrum with the Hagedorn scheme and compare it with a reference calculation obtained from a standard calculation. For both approaches, the spectra are computed from the Fourier transform of the autocorrelation function, a(t):

$$a(t) = \langle \Psi(t_0) | \Psi(t) \rangle \qquad (21)$$

In addition, and to attenuate the oscillations of each band due to the Gibbs phenomenon, a filter ($\cos^2(\pi(t-t_0)/t_f)$, where $t_f$ is the final propagation time) is added to the Fourier transform. For the standard approach, the autocorrelation function calculation is just a scalar product between the wave packet coefficients, $\tilde{C}_I$, at t and $t_0$. However, for the Hagedorn scheme, the wave packets at t and $t_0$ do not share the same basis set. Therefore, before performing the scalar product, the wave packet at t developed on the basis set $\mathbf{B}(\mathbf{\Lambda}_t)$, is projected onto the basis set at $t_0$, $\mathbf{B}(\mathbf{\Lambda}_{t_0})$ using the technique of the third step of the Hagedorn scheme. In figure 9, three spectra are shown, one from the standard propagation with 20 HO basis functions in each dimension, labeled (STD20) and two from the Hagedorn scheme with 5 basis functions in each dimension with different options. For both propagation approaches, the time step is 0.1 au and the propagation time is 60 au. The first spectrum (blue solid line, positive spectrum, STD20) is our reference spectrum and it shows 6-7 peaks more or less equally spaced by 1.0 au.
This reflects the small anharmonicity of the potential. Furthermore, the width of each band is too large to distinguish its components. The second (labeled HAG5TTF) and third (labeled HAG5TTT) spectra are obtained using the Hagedorn scheme with $b_k^{(t)}$ and $p_k^{(t)}$ and without or with re-normalization, respectively. The figure 9 shows that the spectra obtained using both Hagedorn approaches accurately reproduce the spectrum computed with the help of the standard method.

<figure><img src="image_10.png" /><figcaption>Figure 9.: Energy spectrum obtained from a propagation of 60 au using: (i) the standard propagation approach with $nb_k = 20$, labeled STD20 (blue solid line, positive spectrum). (ii) the Hagedorn scheme with $nb_k = 5$, without re-normalization (orange solid line, with orange dot, negative spectrum, labeled HAG5TTF). (iii) the Hagedorn scheme with $nb_k = 5$, with re-normalization (green dashed line, negative spectrum, labeled HAG5TTT).</figcaption></figure>

Finally, to further analyze the Hagedorn scheme, we perform propagations with a larger final propagation time of 600 au instead of 60 au for the previous calculations. Since, we have obtained a good spectrum without re-normalization for short propagation of 60 au, we have tried the long propagation under the same conditions (HAG5TTF). Unfortunately, the propagation have stopped before reaching the final time (the total energy has increased by 0.6 Hartree). Therefore, we have performed the propagation with re-normalization (HAG5TTT). The corresponding spectrum (orange solid line, labeled HAG5TTT) compares very well with the STD20 one (see Fig. 10). In addition, the figure 11 shows the spectrums around 5 au for the standard approach and for the Hagedorn scheme with several basis set sizes. The spectrum obtained with the standard approach shows a main peak and a small side band at higher energy due to the five components around this energy. However, the first Hagedorn spectrum (with $nb_k = 5$, HAG5TTT, left panel) shows a peak on the left that is not present in the reference spectrum. The other two peak positions are well reproduced, although the intensity of the central peak is too small. When the basis set size increases, this additional peak disappears (Fig. 11, central and right panels). The last figure (Fig. 12) of this section shows spectra around 6. au obtained with only 5 basis functions in each dimension for both the standard (STD5) and the Hagedorn (HAG5TTT) approaches. The position of the Hagedorn main peak (marked with an orange vertical line) is in excellent agreement with that of the standard approach obtained with $nb_k = 20$ (blue vertical line). However, the main peak position of the STD5 calculation is strongly overestimated. This is to be expected, since the basis set is too small to reproduce the 6 components of the energy levels around 6 au.
<figure><img src="image_11.png" /><figcaption>Figure 10.: Energy spectrum obtained from a propagation of 600 au using: (i) the standard propagation approach with nb$_{k}$ = 20, labeled STD20 (blue solid line). (ii) the Hagedorn scheme with nb$_{k}$ = 5 and with re-normalization (orange solid line, labeled HAG5TTT).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11.: Energy spectrum around 5 au obtained from a propagation of 600 au using: (i) the standard propagation approach with nb$_{k}$ = 20, labeled STD20 (blue solid line, with blue dot). (ii) the Hagedorn scheme with re-normalization (orange solid line). Left panel with nb$_{k}$ = 5 (HAG5TTT), middle panel with nb$_{k}$ = 10 (HAG10TTT) and right panel with nb$_{k}$ = 20 (HAG20TTT).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12.: Energy spectrum around 6.0 au obtained from a propagation of 600 au using: (i) the standard propagation approach with nb<sub>k</sub> = 5, labeled STD5 (green solid line, with green dot). (ii) the Hagedorn scheme with re-normalization (orange solid line) with nb<sub>k</sub> = 5 (HAG5TTT). In addition, vertical lines show the main peak positions of STD5 (green dashed line), STD20 (blue dashed line) and HAG5TTT (orange dashed line), respectively.</figcaption></figure>

4. Application on the 6D Hénon-Heiles model

In the previous section, the propagation was done with models with two degrees of freedom. In this section, we compare the Hagedorn scheme with the standard approach using the modified Hénon-Heiles model (Eq. 20) with six degrees of freedom [79]. For these calculations, the propagation time is 60 *au* and the time step is 0.1 *au*. The basis set and the initial wave packet parameters are given in the Table 4 and differ slightly from the 2D ones.

<table><caption>Table 4.: Harmonic Oscillator (HO) and Hagedorn basis set parameter values: $q_k^{(0)}$, $p_k^{(0)}$, and $\alpha_k^{(0)}$ (columns 2–6). Initial Gaussian wave packet parameters (Eq. 18, with $e_0 = 1$): $q_k^0$, $p_k^0$, and $a_k^0$ (columns 7–9).</caption><thead><tr><th rowspan="2">k</th><th colspan="2">HO basis set</th><th colspan="3">Hagedorn basis set</th><th colspan="3">Initial wave packet</th></tr><tr><th>q<sub>k</sub><sup>(0)</sup></th><th>a<sub>k</sub><sup>(0)</sup></th><th>q<sub>k</sub><sup>(0)</sup></th><th>p<sub>k</sub><sup>(0)</sup></th><th>α<sub>k</sub><sup>(0)</sup></th><th>q<sub>k</sub><sup>0</sup></th><th>p<sub>k</sub><sup>0</sup></th><th>a<sub>k</sub><sup>0</sup></th></tr></thead><tbody><tr><td colspan="9" style="text-align: center;">6D</td></tr><tr><td>1</td><td>0.0</td><td>1.0</td><td>2.0</td><td>0.0</td><td>(1.0, 0.0)</td><td>2.0</td><td>0.0</td><td>1.0</td></tr><tr><td>2-6</td><td>0.0</td><td>1.0</td><td>0.0</td><td>0.0</td><td>(1.0, 0.0)</td><td>0.0</td><td>0.0</td><td>1.0</td></tr></tbody></table>

The figure 13 shows the spectrum (blue solid line, SDT10) for the standard propagation with 10 basis functions in each dimension, as well as, the spectra for the Hagedorn scheme with nb<sub>k</sub> = 5 (green dashed line, HAG5TTT and orange solid line, with orange dot, HAGTTF). As with the 2D-model (Fig. 9), and with this short propagation time, the agreement is almost perfect. However, two non-physical peaks appear at an energy below the first physical level (around 3. *au*) in the spectra obtained with re-
normalization (HAG5TTT). It is important to note that these two peaks are absent from the spectra obtained without re-normalization (HAG5TTF). We believe, these two peaks are due to the deformation of the autocorrelation function resulting from the wave packet re-normalization. To confirm this deformation, we have computed the error of the Hagedorn autocorrelation function with respect to the standard autocorrelation function as follows:

$$Err = \frac{\int_{0}^{t_f} |a^{\text{STD}}(t) - a^{\text{HAG}}(t)| dt}{t_f}, \qquad (22)$$

where, $a^{\text{STD}}(t)$ and $a^{\text{HAG}}(t)$ denote the auto-correlation functions obtained from the standard and Hagedorn calculations, respectively. The errors associated to the Hagedorn calculation with re-normalization (HAG5TTT) is about twice larger than the one without re-normalization (HAG5TTF). Furthermore, the intensities of these non-physical peaks are reduced when the basis set sizes increase. Figure 14 illustrates this reduction as the basis set sizes increase from $nb_k = 5$ to $nb_k = 7$. This feature is expected, as a larger basis set reduces the loss of wavepacket norm after the third step (see this effect for the Harmonic model in Fig 6).

Finally, the last figure (Fig. 15) shows spectra around 7. and 8. Hartree obtained with only 5 basis functions in each dimension for both the standard (STD5) and the Hagedorn (HAG5TTF) approaches and as well as the reference spectrum (STD10). As for the 2D-model, the position of the Hagedorn main peak (marked with an orange vertical line) is in excellent agreement with that of the standard approach obtained with $nb_k = 10$ (blue vertical line). However, the main peak positions of the STD5 calculation are slightly overestimated at 7. au and strongly overestimated at 8. au. For this 6D-anharmonic model, the Hagedorn scheme allows the use of 1D-basis sets, which are approximately two times smaller than those used in the standard approach. Therefore, the 6D-Hagedorn basis set is approximately 64 times smaller than the one used in the standard scheme. However, in the Hagedorn scheme, the potential representation must to be recomputed at each time step because the basis set (and thus the grid) changes. Nevertheless, the total Hagedorn propagation is approximately three times faster than the standard propagation. Regarding the accuracy of the peak positions obtained from the HAG5TTF calculation compared to the STD10 calculation, the largest error is approximately $3 \times 10^{-3}$ au (a relative error of 0.03 %) for the peak around 8. au.

## 5. Conclusion and perspectives

In this study, we present a numerical implementation of the Hagedorn wave packet with the possibility to use a combination of time-dependent basis sets (Hagedorn ones) and time-independent basis sets (standard ones) with one or several diabatic electronic states. Some remarks are in order with respect to the actual numerical implementation:

* The propagation is performed in three steps: (i) standard propagation (ii) parameter computation of the time-dependent basis set, assuming the wave packet remains Gaussian (iii) projection of the wave packet onto the new basis set. On the one hand, this strategy does not rely on the time-dependent variational principle, so we do not need to compute an inverse of the overlap matrix, which may be singular. On the other hand, this scheme is not optimal because the time-dependent parameters are computed assuming that the wave packet remains Gaussian, which
<figure><img src="image_14.png" /><figcaption>Figure 13.: Energy spectrum obtained from a propagation of 60 au using: (i) the standard propagation approach with nb$_{k}$ = 10, labeled STD10 (blue solid line, positive spectrum). (ii) the Hagedorn scheme with nb$_{k}$ = 5, without re-normalization (orange solid line, with orange dot, negative spectrum, labeled HAG5TTF). (iii) the Hagedorn scheme with nb$_{k}$ = 5, with re-normalization (green dashed line, negative spectrum, labeled HAG5TTT).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14.: Energy spectrum obtained from a propagation of 60 au for energies below 2.5 au with the Hagedorn scheme with re-normalization and with nb$_{k}$ = 5 (blue solid line, labeled HAG5TTT) or with nb$_{k}$ = 7 (orange dashed line, labeled HAG7TTT).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15.: Energy spectra around 7.0 au (left panel) and around 8.0 au (right panel) obtained from a propagation of 60 au with: (i) the standard propagation approach with $nb_k = 5$, labeled STD5 (green solid line), (ii) the standard propagation approach with $nb_k = 10$, labeled STD10 (blue solid line). (iii) the Hagedorn scheme with re-normalization (orange solid line) with $nb_k = 5$ (HAG5TTF). In addition, vertical lines show the main peak positions of the STD5 (green dashed line), the STD10 (blue dashed line) and the HAG5TTF (orange dashed line), respectively.</figcaption></figure>

is not the case for realistic systems. Therefore, we are planning to implement the time-dependent variational principle [80] to propagate the time-dependent parameters.

* Since we are using grid, we are not formally restricted in the choice of the Hamiltonian form. Indeed, we can use simple analytical models (as we did in the present study), but also more complex ones for which the potentials do not have simple expressions. Furthermore, we can easily extend our implementation so that complex curvilinear coordinates can be used with analytical kinematic energy operator (KEO) expressions [81–88] or with a KEO obtained numerically [89–95]. Indeed, the first step of our propagation scheme can be used with any Hamiltonian form (as it is done in discrete variable representation or pseudo-spectral approaches [96–101]) and the other two steps are independent of the Hamiltonian.

This implementation has been tested on several models. In particular, it was able to reproduce the Heller scheme for the quadratic model. This test was important to check the numerical stability of our three-step scheme. The tests on anharmonic systems (2D and 6D) show that it was able to reproduce a standard propagation scheme with high accuracy. Furthermore, the spectra, computed from the Fourier transform of the autocorrelation function, are also well reproduced, but with a smaller basis set than the one required for the standard propagation scheme. In addition, with the direct product expansion of the wave packet, the number of operations required to compute the Hamiltonian action on a wave packet grows as $nc \cdot n_k^{nc+1}$ [102]. Therefore and assuming a 1D-Hagedorn basis set two times smaller than the standard one, we can expect a speedup of approximately of $2^{nc+1}$ for the Hagedorn scheme with respect of the standard one. However, this has to be reduced since the potential must be recalculated with each time step in the Hagedorn scheme.
As mentioned in section 2, our implementation allows the study of models with several electronic states. To show this feature, we test our implementation with the 2D-model of retinal photoisomerization developed by Susanne Hahn and Gerhard Stock [7] (see the appendix B). This model has played a key role in understanding ultrafast photochemical processes in biological systems and was later extended to a 2+23D model that includes additional vibrational modes, providing a more comprehensive description of the system [6]. Over the years, this model has been extensively studied and refined, contributing significantly to the field of photoinduced molecular dynamics [10, 103–106].

This model is described by a 2 × 2 Hamiltonian in the diabatic basis with two coordinates, $\Phi$ associated with the torsion (it describes the cis-trans isomerization) and $Q$ associated with the active vibronic coupling mode. For this system, we use a Fourier basis set for the $\Phi$ coordinate (with $nb_{\Phi} = nq_{\Phi} = 256$) for both the standard and Hagedorn propagation approaches. For the coordinate $Q$, we use the HO basis set and the Hagedorn one for the standard and Hagedorn propagation schemes, respectively. The size of the basis set is denoted by $nb_Q$ and the corresponding grid size, $nq_Q$, is $nb_Q + 5$. Initially, the system is assumed to have already absorbed light from the cis isomer. Therefore, the initial wave packet is expressed as a 2D-Gaussian on the second diabatic electronic state ($e_0 = 2$). The basis set and initial wave packet parameters are given in the appendix B.

For the comparison between the standard and the Hagedorn propagation scheme, we evaluate the populations, $p_e$, along time for each diabatic state $e$ ($p_e = \langle\psi^{(e)}(t) | \psi^{(e)}(t)\rangle$). For a propagation of 10000. au (with a time step of 1. au), the numerical convergence is reached with $nb_Q = 20$ and $nb_Q = 50$ for the standard and Hagedorn propagation schemes, respectively (see the first diabatic state population in the Fig. 16). The fact that more basis functions are needed for the Hagedorn scheme seems surprising. However, since we have a common basis set for both diabatic states (equivalent to the single-set formalism in MCTDH), the Hagedorn basis parameters are better adapted to the state with a large population (here the second one). To overcome this difficulty, we are planing to describe the wave packet with different basis sets for the different electronic states (equivalent to the multi-set formalism in MCTDH). Furthermore, the spawning technique, which was originally developed for a wave packet described by a sum of multidimensional Gaussian wave function [27, 107], can be adapted to the Hagedorn wave packet approach, as shown by Gradinaru et al. [108].

In the present study, the initial wave packet is described as a product of 1D-Gaussian along the coordinates. Although, this approximation is widely used (in particular for non-adiabatic dynamics [14]), a relaxation approach can easily be implemented to obtain the vibrational or vibronic ground state with the Hagedorn approach. Furthermore, the photoinduced processes can be modelized through a propagation: (i) for which the initial wave packet is obtained by the action of the transition dipole moment operator on the ground vibronic state. (ii) with a time-dependent Hamiltonian[35] including the time-dependent electric field explicitly with a semiclassical field/matter description. As mentioned in the introduction, the aim of this study is the development of a Hagedorn wave packet technique and therefore at a temperature of 0 K. Developing such Hagedorn technique to solve Liouville–von Neumann equation is certainly doable but beyond the scope of this study. Although the Monte Carlo Wave-Function (MCWF)[109] or the Variational Wave Packet (VWP)[110] approaches have been developed to efficiently solve the Liouville–von Neumann equation.
<figure><img src="image_17.png" /><figcaption>Figure 16.: Time evolution of the ground state population during retinal isomerization obtained from two computational approaches: (i) the standard method (STD20, blue dashed line) and (ii) the Hagedorn approach with two variants (HAG20TTT, orange solid line with solid dot; HAG50TTT, green solid line with squares).</figcaption></figure>

Furthermore, the actual implementation relies on a direct-product of the basis sets and grids. Therefore, it is limited to small system sizes since the total number of basis functions, NB, or grid points, NQ, grows exponentially with the number of degrees of freedom, nc, although this is somehow attenuated since a small basis set size could be used when spectra are computed. To overcome this curse of dimensionality, a basis function selection scheme has already been implemented for Hagedorn wave packets using algebraic approaches [53]. However, in our actual implementation, a standard approach with a direct-product multidimensional basis set is used. A pruned basis set or a basis set with a selection of basis functions is insufficient since a multidimensional grid is required. Therefore, a Smolyak scheme [111] combined with standard propagation techniques[10, 50, 112, 113] could be well-adapted, as it avoids exponential scaling for both the basis set and the grid. The Smolyak technique was first applied to vibrational spectroscopy by Avila and Carrington, followed by several other groups [114–119]. Additionally, the Smolyak scheme has been applied to wave packet propagation, both with and without non-adiabatic couplings [10, 50, 112, 113]. Therefore, we are also planning to combine the Smolyak scheme with the Hagedorn propagation technique. With this improvement, we should be able to study larger systems, such as spin-boson models, although with a finite bath size. The Smolyak approach has already been applied to system/bath models. [10, 118].

## Acknowledgments

K. S., R. I. and K. M. R. A. thank the Institut de Chimie Physique (ICP) of the Université Paris Saclay and the CNRS, France, for computing facilities. R. I. and K. M. R. A. gratefully acknowledge the ICP for hosting our research stays and also the ThéoSim
group for their scientific support and hospitality. This work was supported through mobility grants from the Service de Coopération et d'Action Culturelle (SCAC) of the French Embassy in Togo. We additionally thank Campus France for their administrative assistance during our stays in France.

## Disclosure statement

No potential conflict of interest was reported by the author(s).

## Data availability statement

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## Appendix A. Appendix: Projection

Once the new basis is constructed, we can proceed with the projection. From Eqs 9,14,16, we recall that, the $C_I^{t_2}$ are given by:

$$C_J^{t_2} = \langle B_J(\mathbf{\Lambda}_{t_2}) | \Psi(t_2) \rangle = \sum_{I=1}^{NB} \tilde{C}_I^{t_2} \langle B_J(\mathbf{\Lambda}_{t_2}) | B_I(\mathbf{\Lambda}_{t_1}) \rangle, \quad (A1)$$

where, $\tilde{C}_I^{t_2}$ and $C_J^{t_2}$ are the wave packet coefficients at time $t_2$ associated, respectively, to the $I^{th}$ basis function at time $t_1$ and $J^{th}$ basis function at time $t_2$. This equation shows that the projection is performed on a multidimensional basis. The key elements are the overlaps between the basis functions at different times, $t_1$ and $t_2$. However, numerically, it can be achieved sequentially using partial projections (one coordinate after another). To perform the sequential transformations, the coefficients $\tilde{C}_I^{t_2}$ and $C_J^{t_2}$ notations need to be re-expressed as :

$$\begin{aligned} \tilde{C}_I^{t_2} &= C_{i_1, \dots, i_k, \dots, i_{nc}}^{b_{t_1}^1 \cdots b_{t_1}^k \cdots b_{t_2}^{nc}}(t_2), \\ C_I^{t_2} &= C_{i_1, \dots, i_k, \dots, i_{nc}}^{b_{t_2}^1 \cdots b_{t_2}^k \cdots b_{t_2}^{nc}}(t_2). \end{aligned} \quad (A2)$$

This notation enables to show that the coefficients come from primitive basis sets at different times $t_1$ or $t_2$. Furthermore, we recall that:

$$\begin{aligned} B_I(\mathbf{Q}, \mathbf{\Lambda}_{t_1}) &= b_{i_1}^1(q_1, t_1) \cdots b_{i_k}^k(q_k, t_1) \cdots b_{i_{nc}}^{nc}(q_{nc}, t_1), \\ B_I(\mathbf{Q}, \mathbf{\Lambda}_{t_2}) &= b_{i_1}^1(q_1, t_2) \cdots b_{i_k}^k(q_k, t_2) \cdots b_{i_{nc}}^{nc}(q_{nc}, t_2). \end{aligned} \quad (A3)$$

The functions $b_{i_k}^k(q_k, t_1) = b_{i_k}^k(q_k, \lambda_k^{(t_1)})$ and $b_{i_k}^k(q_k, t_2) = b_{i_k}^k(q_k, \lambda_k^{(t_2)})$ represent the $i_k^{th}$ function of the primitive basis $k$ at time $t_1$ and time $t_2$, respectively. With these notations, the only difference between the $\tilde{C}_I^{t_2}$ and $C_I^{t_2}$ or between the $B_I(\mathbf{Q}, \mathbf{\Lambda}_{t_1})$ and $B_I(\mathbf{Q}, \mathbf{\Lambda}_{t_2})$ are the time indices ($t_1$ or $t_2$) of the $\lambda_k^{th}$ time dependent basis parameters ($k$).
The partial projection for the coordinate k reads:

$$C_{j_1, \dots, j_k, t_{k+1}, \dots, i_{nc}}^{b_{t_2}^1 \dots b_{t_2}^k b_{t_1}^{k+1} \dots b_{t_1}^{nc}}(t_2) = \sum_{i_k} \langle b_{j_k}^k(t_2) | b_{i_k}^k(t_1) \rangle \\ \times C_{j_1, \dots, j_{k-1}, i_k, \dots, i_{nc}}^{b_{t_2}^1 \dots b_{t_2}^{k-1} b_{t_1}^k \dots b_{t_1}^{nc}}(t_2) \quad (A4)$$

For this projection, some remarks are in order:

(1) The partial projection starts with $k = 1$ and ends at $k = nc$.

(2) For a given $k$, the partial projection has to be applied for all values of the indices $\{j_1, \dots, j_k, i_{k+1}, \dots, i_{nc}\}$.

(3) The overlap elements obtained by nD-multidimensional integrals, $\langle B_J(\Lambda_{t_2}) | B_I(\Lambda_{t_1}) \rangle$, are substituted by several 1D-integrals, $\langle b_{j_k}^k(t_2) | b_{i_k}^k(t_1) \rangle \equiv \langle b_{j_k}^k(\lambda_k^{(t_2)}) | b_{i_k}^k(\lambda_k^{(t_1)}) \rangle$.

(4) These 1D-integrals are just the overlaps, $S^k(j_k, i_k)$, between the $k^{th}$ 1D-basis sets at different time.

The calculation of the $S^k(j_k, i_k)$ overlap matrix elements are performed numerically using Gauss-Hermite quadratures:

$$S^k(j_k, i_k) = \langle b_{j_k}^k(t_2) | b_{i_k}^k(t_1) \rangle \\ \approx \sum_{u=1}^{nq_k} b_{j_k}^k(q_k^u, t_2) \cdot b_{i_k}^k(q_k^u, t_1) \cdot w_k^u \quad (A5)$$

In Eq. A5, the $q_k^u$ and the $w_k^u$ ($u \in \{1 \cdots nq_k\}$) represent the Gauss quadrature, $g^k$, for the shifted and scaled Harmonic Oscillator basis set. In term of standard Harmonic Oscillator grid ($x_u$ and $w_u$), $q_k^u$ and $w_k^u$ are expressed as follows:

$$q_k^u = q_k^{(t)} + x_u / \sqrt{a_k^{(t)}}, \quad w_k^u = w_u / \sqrt{a_k^{(t)}} \quad (A6)$$

It is important to note that the grid used to compute $S^k(j_k, i_k)$ is neither the grid associated with time $t_1$ nor the one associated with time $t_2$. Indeed, since the product of two Gaussians is also a Gaussian, this grid is an intermediate grid constructed from the time-dependent basis set parameters at times $t_1$ and $t_2$. More precisely, it can be shown that the intermediate parameters, $a_k^{Int}$ and $q_k^{Int}$ are:

$$a_k^{Int} = \frac{a_k^{(t_1)} + a_k^{(t_2)}}{2} \quad (A7)$$

$$q_k^{Int} = \frac{a_k^{(t_1)} \cdot q_k^{(t_1)} + a_k^{(t_2)} \cdot q_k^{(t_2)}}{a_k^{(t_1)} + a_k^{(t_2)}} \quad (A8)$$

Remarks: (i) we do not transform the basis sets $b^k(q_k, \lambda_k^{(t_1)})$ and $b^k(q_k, \lambda_k^{(t_2)})$ to a new basis set with the intermediate parameters, but we just compute $b_{j_k}^k(q_k^u, \lambda_k^{(t_2)})$ and $b_{i_k}^k(q_k^u, \lambda_k^{(t_1)})$ on the new grid points. (ii) When, $b_k^{(t)}$ (imaginary part of $\alpha_k^{(t)}$) and $p_k^{(t)}$ are not taken into account, the computation of $S^k(j_k, i_k)$ are numerically exact as soon as $nq_k = nb_k$. (iii) To reach a good accuracy, we chose $nq_k = nb_k + 5$.
Appendix B. Appendix: 2D-retinal model

The 2D-retinal model is described by a 2 × 2 Hamiltonian, $\hat{H}^d$, in the diabatic basis, expressed as:

$$ \hat{H}^d(\boldsymbol{Q}) = \hat{T}(\boldsymbol{Q}) + \begin{bmatrix} V_{11}^d(\Phi, Q) & V_{12}^d(Q) \\ V_{21}^d(Q) & V_{22}^d(\Phi, Q) \end{bmatrix}, \qquad (B1) $$

where the potentials for the two diabatic electronic states and their coupling are defined by:

$$ V_{11}^d(\Phi, Q) = \frac{1}{2}W_1(1 - \cos \Phi) + \frac{1}{2}\omega Q^2, \qquad (B2) $$

$$ V_{22}^d(\Phi, Q) = E - \frac{1}{2}W_2(1 - \cos \Phi) + \frac{1}{2}\omega Q^2 + \kappa Q, \qquad (B3) $$

$$ V_{12}^d(Q) = V_{21}(Q) = \lambda Q. \qquad (B4) $$

Here, $m$ represents the effective mass of the reaction coordinate while $\omega$ denotes the frequency of the coupling mode. The parameter $\lambda$ corresponds to the inter-state coupling, and $W_0$, $W_1$, and $E$ are adjusted to reproduce the properties of the electronic states $S_1$ and $S_2$. The coordinate $\Phi$ represents the mode along which the molecule undergoes isomerization, and $Q$ corresponds to the active vibronic mode coupling the two electronic states. The kinetic energy operator, assumed to be the same for both diabatic electronic states, is given by:

$$ \hat{T}(\boldsymbol{Q}) = -\frac{1}{2m} \frac{\partial^2}{\partial \Phi^2} - \frac{\omega}{2} \frac{\partial^2}{\partial Q^2}. \qquad (B5) $$

The values of the model parameters, determined through fitting to experimental data and theoretical calculations, are:

* the inverse of the effective mass, $m^{-1} = 4.84 \times 10^{-4}$ eV

* the torsional barriers $W_0 = 3.6$ eV and $W_1 = 1.09$ eV

* the linear coupling $\kappa = 0.1$ eV

* the vibronic frequency $\omega = 0.19$ eV

* the inter-state coupling $\lambda = 0.19$ eV,

* the offset energy $E = 2.48$ eV.

The basis set and initial wave packet parameters are given in the table B1.

<table><caption>Table B1.: Hagedorn basis set parameter values: $q_k^{(0)}$, $p_k^{(0)}$, and $\alpha_k^{(0)}$ (columns 2–4). Initial Gaussian wave packet parameters (Eq. 18, with $e_0 = 2$): $q_k^0$, $p_k^0$, and $\Re(a_k^0)$ (columns 5–7).</caption><thead><tr><th rowspan="2">k (coord.)</th><th colspan="3">Hagedorn and HO basis set</th><th colspan="3">Initial wave packet</th></tr><tr><th>q<sub>k</sub><sup>(0)</sup></th><th>p<sub>k</sub><sup>(0)</sup></th><th>α<sub>k</sub><sup>(0)</sup></th><th>q<sub>k</sub><sup>0</sup></th><th>p<sub>k</sub><sup>0</sup></th><th>Re(a<sub>k</sub><sup>0</sup>)</th></tr></thead><tbody><tr><td>1 (Φ)</td><td>/</td><td>/</td><td>/</td><td>0.0</td><td>0.0</td><td>60.9836267</td></tr><tr><td>2 (Q)</td><td>0.0</td><td>0.0</td><td>0.9202033</td><td>0.0</td><td>0.0</td><td>0.9202033</td></tr></tbody></table>
References

[1] B. Walker, B. Sheehy, L. Dimauro, P. Agostini, K. Schafer and K. Kulander, Physical Review Letters **73** (9), 1227–1230 (1994).

[2] A. Rudenko, K. Zrost, B. Feuerstein, V. DeJesus, C. Schröter, R. Moshammer and J. Ullrich, Physical Review Letters **93** (25), 253001 (2004).

[3] E. Engel, N. Doss, G. Harris and J. Tennyson, Monthly Notices of the Royal Astronomical Society **357** (2), 471–477 (2005).

[4] H. Pedersen, S. Altevogt, B. Jordon-Thaden, O. Heber, M. Rappaport, D. Schwalm, J. Ullrich, D. Zajfman, R. Treusch, N. Guerassimova et al., Physical Review Letters **98** (22), 223202 (2007).

[5] K. Sodoga, J. Loreau, D. Lauvergnat, Y. Justum, N. Vaeck and M. Desouter-Lecomte, Physical Review A **80** (3), 033417 (2009).

[6] S. Hahn and G. Stock, Chemical Physics **259** (2-3), 297–312 (2000).

[7] S. Hahn and G. Stock, the Journal of Physical Chemistry B **104** (6), 1146–1149 (2000).

[8] E. Marsili, M. Olivucci, D. Lauvergnat and F. Agostini, Journal of Chemical Theory and Computation **16** (10), 6032–6048 (2020).

[9] M. Manathunga, X. Yang, Y. Orozco-Gonzalez and M. Olivucci, The Journal of Physical Chemistry Letters **8** (20), 5222–5227 (2017).

[10] A. Pereira, J. Knapik, A. Chen, D. Lauvergnat and F. Agostini, The European Physical Journal Special Topics **232** (12), 1917–1933 (2023).

[11] M. Kimura and C.D. Lin, Physical Review A **34** (1), 176 (1986).

[12] A. Jain, C.D. Lin and W. Fritsch, Physical Review A **36** (5), 2041 (1987).

[13] J. Loreau, K. Sodoga, D. Lauvergnat, M. Desouter-Lecomte and N. Vaeck, Physical Review A **82** (1), 012708 (2010).

[14] L.L.E. Cigrang, B.F.E. Curchod, R.A. Ingle, A. Kelly, J.R. Mannouch, D. Accomasso, A. Alijah, M. Barbatti, W. Chebbi, N. Došlić, E.C. Eklund, S. Fernandez-Alberti, A. Freibert, L. González, G. Granucci, F.J. Hernández, J. Hernández-Rodríguez, A. Jain, J. Janoš, I. Kassal, A. Kirrander, Z. Lan, H.R. Larsson, D. Lauvergnat, B. Le Dé, Y. Lee, N.T. Maitra, S.K. Min, D. Peláez, D. Picconi, Z. Qiu, U. Raucci, P. Robertson, E. Sangiogo Gil, M. Sapunar, P. Schürger, P. Sinnott, S. Tretiak, A. Tikku, P. Vindel-Zandbergen, G.A. Worth, F. Agostini, S. Gómez, L.M. Ibele and A. Prlj, The Journal of Physical Chemistry A (2025).

[15] B. Wu, B. Li, X. He, X. Cheng, J. Ren and J. Liu, Journal of Chemical Theory and Computation **21** (8), 3775–3813 (2025).

[16] W.H. Miller, The Journal of Physical Chemistry A **105** (13), 2942–2955 (2001).

[17] J. Liu and W.H. Miller, The Journal of Chemical Physics **126** (23), 234110 (2007).

[18] P. Frantsuzov, A. Neumaier and V.A. Mandelshtam, Chemical Physics Letters **381** (1-2), 117–122 (2003).

[19] J.C. Tully, The Journal of Chemical Physics **93** (2), 1061–1071 (1990).

[20] D.F. Coker and L. Xiao, The Journal of Chemical Physics **102** (1), 496–510 (1995).

[21] R. Kapral, Chemical Physics **481**, 77–83 (2016).

[22] F. Agostini, S.K. Min, A. Abedi and E.K.U. Gross, Journal of Chemical Theory and Computation **12** (5), 2127–2143 (2016).

[23] J.E. Runeson and D.E. Manolopoulos, The Journal of Chemical Physics **159** (9), 104111 (2023).

[24] T.J. Martinez, M. Ben-Nun and R.D. Levine, The Journal of Physical Chemistry **100** (19), 7884–7895 (1996).

[25] G.A. Worth, M.A. Robb and I. Burghardt, Faraday discussions **127**, 307–323 (2004).

[26] J. Liu and W.H. Miller, The Journal of Chemical Physics **134** (10), 104102 (2011).
[27] B.F.E. Curchod and T.J. Martinez, Chemical Reviews **118** (7), 3305–3336 (2018).

[28] H.D. Meyer, U. Manthe and L.S. Cederbaum, Chemical Physics Letters **165** (1), 73–78 (1990).

[29] H. Wang and M. Thoss, The Journal of Chemical Physics **119** (3), 1289–1299 (2003).

[30] U. Manthe, The Journal of Chemical Physics **128** (16), 164116 (2008).

[31] O. Vendrell and H.D. Meyer, The Journal of Chemical Physics **134** (4), 044135 (2011).

[32] C. Leforestier, R.H. Bisseling, C. Cerjan, M.D. Feit, R. Friesner, A. Gulderberg, A. Hammerich, G. Jolicard, W. Karrlein, H.D. Meyer et al., Journal of Computational Physics **94** (1), 59–80 (1991).

[33] M.D. Feit, J.A. Fleck and A. Steiger, Journal of Computational Physics **47** (3), 412–433 (1982).

[34] H. Tal-Ezer and R. Kosloff, The Journal of Chemical Physics **81** (9), 3967–3971 (1984).

[35] D. Lauvergnat, S. Blasco, X. Chapuisat and A. Nauts, The Journal of chemical physics **126** (20), 204103 (2007).

[36] E.J. Heller, The Journal of Chemical Physics **64** (1), 63–73 (1976).

[37] D. Huber and E.J. Heller, The Journal of chemical physics **87** (9), 5302–5311 (1987).

[38] G. Drolshagen and E.J. Heller, The Journal of chemical physics **79** (4), 2072–2082 (1983).

[39] J. Vaníček and T. Begušić, Chapter 6 - Ab Initio Semiclassical Evaluation of vibrationally Resolved Electronic Spectra With Thawed Gaussians in Molecular Spectroscopy and Quantum Dynamics Elsevier 2021.

[40] J. Sielk, H.F.V. Horsten, F. Krüger, R. Schneider and B. Hartke, Physical Chemistry Chemical Physics **11** (3), 463–475 (2009).

[41] S. Choi and J. Vaníček, The Journal of Chemical Physics **151** (23) (2019).

[42] A.D. Hammerich, R. Kosloff and M.A. Ratner, Chemical Physics Letters **171** (1-2), 97–108 (1990).

[43] H.D. Meyer and G.A. Worth, Theoretical Chemistry Accounts **109**, 251–267 (2003).

[44] D. Peláez and H.D. Meyer, Chemical Physics **482**, 100–105 (2017).

[45] L. Shi, M. Schröder, H.D. Meyer, D. Peláez, A.M. Wodtke, K. Golibrzuch, A.M. Schönemann, A. Kandratsenka and F. Gatti, The Journal of Chemical Physics **159** (19), 194102 (2023).

[46] G. Hagedorn, Annals of Physics **135** (1), 58–70 (1981).

[47] G.A. Hagedorn, V. Rousse and S.W. Jilcott, Annales Henri Poincaré **7**, 1065–1083 (2006).

[48] G. Hagedorn, Annals of Physics **269** (1), 77–104 (1998).

[49] E. Faou, V. Gradinaru and C. Lubich, SIAM Journal on Scientific Computing **31** (4), 3027–3041 (2009).

[50] C. Lasser and C. Lubich, Acta Numerica **29**, 229–401 (2020).

[51] Z.T. Zhang and J.J. Vaníček, The Journal of Chemical Physics **161** (11), 111101 (2024).

[52] H. Dietert, J. Keller and S. Troppmann, Journal of Mathematical Analysis and Applications **450** (2), 1317–1332 (2017).

[53] J.J.L. Vaníček and Z.T. Zhang, Journal of Physics A: Mathematical and Theoretical **58** (8), 085303 (2025).

[54] T. Ohsawa, Journal of Fourier Analysis and Applications **25**, 1513–1552 (2019).

[55] G.A. Hagedorn, Spectral Analysis, Differential Equations and Mathematical Physics: A Festschrift in Honor of Fritz Gesztesy’s 60th Birthday **87**, 183–190 (2013).

[56] T. Ohsawa and M. Leok, Journal of Physics A: Mathematical and Theoretical **46** (40), 405201 (2013).

[57] V. Gradinaru and G.A. Hagedorn, Numerische Mathematik **126** (1), 53–73 (2014).

[58] C. Lasser and S. Troppmann, Journal of Fourier Analysis and Applications **20**, 679–714 (2014).

[59] X. Li and A. Xiao, International Journal of Modeling, Simulation, and Scientific Computing **5** (04), 1450013 (2014).

[60] T. Ohsawa, Journal of Mathematical Physics **56** (3) (2015).

[61] P. Punoševac and S.L. Robinson, Journal of Mathematical Physics **57** (9), 092102 (2016).

[62] R. Issa, Wave packet propagation Fortran code (standard and Hagedorn schemes) https:
//github.com/asodoga/TD_Schrod_Rabiou 2025, Accessed: May 02, 2025.

[63] T. Mercouris, Y. Komninos, S. Dionissopoulou and C.A. Nicolaides, Physical Review A 50, 4109 (1994).

[64] E.J. Heller, Accounts of Chemical Research **14** (12), 368–375 (1981).

[65] E.J. Heller, The Journal of Chemical Physics **62** (4), 1544–1555 (1975).

[66] R. Moghaddasi Fereidani and J.J. Vaníček **159** (9), 094114 (2023).

[67] B. Lasorne, G.A. Worth and M.A. Robb, Wiley Interdisciplinary Reviews: Computational Molecular Science **1** (3), 460–475 (2011).

[68] B. Lasorne, G.A. Worth and M.A. Robb, in *Molecular Quantum Dynamics: From Theory to Applications*, edited by Fabien Gatti (Springer Berlin Heidelberg, Berlin, Heidelberg, 2014), pp. 181–211.

[69] B. Lasorne, M.A. Robb and G.A. Worth, Chemical Physics Letters **432** (4-6), 604–609 (2006).

[70] B. Lasorne, M.A. Robb and G.A. Worth, Physical Chemistry Chemical Physics **9** (25), 3210–3227 (2007).

[71] D. Mendive-Tapia, B. Lasorne, G.A. Worth, M.A. Robb and M.J. Bearpark, The Journal of Chemical Physics **137** (22), 22A513 (2012).

[72] C.S. Allan, B. Lasorne, G.A. Worth and M.A. Robb, The Journal of Physical Chemistry A **114** (33), 8713–8729 (2010).

[73] M. Araújo, B. Lasorne, A.L. Magalhães, G.A. Worth, M.J. Bearpark and M.A. Robb, The Journal of Chemical Physics **131** (14), 144301 (2009).

[74] B. Lasorne, M.J. Bearpark, M.A. Robb and G.A. Worth, The Journal of Physical Chemistry A **112** (50), 13017–13027 (2008).

[75] M. Araujo, B. Lasorne, A.L. Magalhães, M.J. Bearpark and M.A. Robb, The Journal of Physical Chemistry A **114** (45), 12016–12020 (2010).

[76] G.W. Richings, I. Polyak, K.E. Spinlove, G.A. Worth, I. Burghardt and B. Lasorne, International Reviews in Physical Chemistry **34** (2), 269–308 (2015).

[77] M.H. Beck, A. Jäckle, G.A. Worth and H.D. Meyer, Physics Reports **324** (1), 1–105 (2000).

[78] H.D. Meyer, U. Manthe and L. Cederbaum, The Multiconfiguration Time-Dependent Hartree (MCTDH) Method. https://www.pci.uni-heidelberg.de/tc/mctdh.html.

[79] M. Nest and H.D. Meyer, The Journal of Chemical Physics **117** (23), 10499 (2002).

[80] A. McLachlan and M. Ball, Reviews of Modern Physics **36** (3), 844 (1964).

[81] E.B. Wilson, J.C. Decius and P.C. Cross, *Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra*, mcgraw-hil ed. New York.

[82] R. Meyer and H.H. Günthard, The Journal of Chemical Physics **49** (4), 1510–1520 (1968).

[83] X. Chapuisat, A. Nauts and J.P. Brunet, Molecular Physics **72** (1), 1–31 (1991).

[84] X. Chapuisat and C. Iung, Physical Review A **45** (9), 6217–6235 (1992).

[85] X.G. Wang and T. Carrington, The Journal of Chemical Physics **113** (17), 7097–7101 (2000).

[86] F. Gatti, C. Muñoz and C. Iung, Journal of Chemical Physics **114** (19), 8275–8282 (2001).

[87] M. Ndong, L. Joubert-Doriol, H.D. Meyer, A. Nauts, F. Gatti and D. Lauvergnat, The Journal of Chemical Physics **136** (3), 034107 (2012).

[88] M. Ndong, A. Nauts, L. Joubert-Doriol, H.D. Meyer, F. Gatti and D. Lauvergnat, The Journal of Chemical Physics **139** (20), 204107 (2013).

[89] R. Meyer, Journal of Molecular Spectroscopy **76**, 266–300 (1979).

[90] M. Harthcock and J. Laane, Journal of Molecular Spectroscopy **324**, 300–324 (1982).

[91] M.L. Senent, D.C. Moule and Y.G. Smeyers, Journal of Molecular Spectroscopy **372**, 257–266 (1995).

[92] D. Luckhaus, The Journal of Chemical Physics **113** (4), 1329–1347 (2000).

[93] D. Lauvergnat and A. Nauts, The Journal of Chemical Physics **116** (19), 8560–8570 (2002).
[94] E. Mátyus, G. Czakó and A.G. Császár, The Journal of chemical physics **130** (13), 134112 (2009).

[95] E. Marsili, F. Agostini, A. Nauts and D. Lauvergnat, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences **380** (2223), 20200388 (2022).

[96] D. Gottlieb and S.A. Orszag, Numerical analysis of spectral methods: theory and applications SIAM 1977.

[97] J.C. Light and T. Carrington, Advances in Chemical Physics **114**, 263–310 (2000).

[98] R.A. Friesner, The Journal of Chemical Physics **85** (3), 1462–1468 (1986).

[99] V. Szalay and P. Ádám, The Journal of Chemical Physics **137** (6), 064118 (2012).

[100] H.G. Yu, The Journal of Chemical Physics **122** (16), 164107 (2005).

[101] F. Gatti, C. Iung, C. Leforestier and X. Chapuisat, The Journal of chemical physics **111** (16), 7236–7243 (1999).

[102] J.C. Light and T. Carrington Jr, Advances in Chemical Physics **114**, 263–310 (2000).

[103] P.J.M. Johnson, M.H. Farag, A. Halpin, T. Morizumi, V.I. Prokhorenko, J. Knoester, T.L.C. Jansen, O.P. Ernst and R.J.D. Miller, The Journal of Physical Chemistry B **121** (16), 4040–4047 (2017).

[104] L. Seidner and W. Domcke, Chemical Physics **186** (1), 27–40 (1994).

[105] W. Domcke and G. Stock, Advances in Chemical Physics **100**, 1–169 (1997).

[106] B. Balzer, S. Hahn and G. Stock, Chemical Physics Letters **379** (3-4), 351–358 (2003).

[107] M. Ben-Nun and T.J. Martínez, The Journal of Chemical Physics **108** (17), 7244–7257 (1998).

[108] V. Gradinaru and O. Rietmann, Journal of Computational Physics **509** (April), 113029 (2024).

[109] J. Dalibard, Y. Castin and K. Mølmer, Physical Review Letters **68** (5), 580–583 (1992).

[110] T. Gerdts and U. Manthe, Journal of Chemical Physics **106** (8), 3017–3023 (1997).

[111] S.A. Smolyak, Soviet Mathematics Doklady **4**, 240 (1963).

[112] V. Gradinaru, SIAM Journal on Numerical Analysis **46** (1), 103–123 (2008).

[113] Y. Xie, Y. Yang, X. Zhu, A. Chen and B. Gu, Journal of Chemical Theory and Computation **20** (21), 9512–9521 (2024).

[114] G. Avila and T. Carrington, The Journal of chemical physics **131** (17), 174103 (2009).

[115] G. Avila and T. Carrington, Chemical Physics **482**, 3–8 (2017).

[116] D. Lauvergnat and A. Nauts, Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy **119**, 18–25 (2014).

[117] G. Avila and E. Mátyus, The Journal of Chemical Physics **150** (17), 174107 (2019).

[118] A. Chen, D.M. Benoit, Y. Scribano, A. Nauts and D. Lauvergnat, Journal of Chemical Theory and Computation **18** (7), 4366–4372 (2022).

[119] D. Lauvergnat and A. Nauts, ChemPhysChem **24** (21), e202300501 (2023).