Finite-temperature vibronic spectra from the split-operator coherence thermofield dynamics

We present a numerically exact approach for evaluating vibrationally resolved electronic spectra at finite temperatures using the coherence thermofield dynamics. In this method, which avoids implementing an algorithm for solving the von Neumann equation for coherence, the thermal vibrational ensemble is first mapped to a pure-state wavepacket in an augmented space, and this wavepacket is then propagated by solving the standard, zero-temperature Schr¨odinger equation with the split-operator Fourier method. We show that the finite-temperature spectra obtained with the coherence thermofield dynamics in a Morse potential agree exactly with those computed by Boltzmann-averaging the spectra of individual vibrational levels. Because the split-operator thermofield dynamics on a full tensor-product grid is restricted to low-dimensional systems, we briefly discuss how the accessible dimensionality can be increased by various techniques developed for the zero-temperature split-operator Fourier method.


I. INTRODUCTION
Vibrationally resolved electronic spectroscopy has proven to be a powerful tool for probing molecular systems and understanding photochemical processes.6][7][8] These thermal effects may become more pronounced in two-dimensional spectroscopy experiments. 9][11][12][13][14][15][16] Given the large number of vibrational levels that may be occupied, time-independent sum-over-state expressions become impractical at higher temperatures for more than a few vibrational degrees of freedom.Timedependent approaches offer a promising alternative, in which the spectrum is evaluated as the Fourier transform of an appropriate autocorrelation function without the need to select relevant transitions and compute their Franck-Condon factors. 11,17However, the most widely used quantum mechanical [18][19][20][21][22] and semiclassical 17,[23][24][25][26][27][28] methods were developed only for pure states and cannot directly handle a thermal ensemble of states.
Standard approaches to include temperature effects either directly evolve the thermal density matrix or employ statistical averaging over pure states.An alternative way to represent quantum dynamics at finite temperature, called thermofield dynamics, 29,30 maps the thermal ensemble to a pure state in an augmented space, which can be propagated by conventional wavepacket a) Electronic mail: jiri.vanicek@epfl.chmethods for solving the time-dependent Schrödinger equation.0][41][42][43] In a similar spirit, Grossmann performed exact calculations of the thermal dipole time correlation function for vibrational spectra by propagating two thermal wavepackets and evaluating the matrix element of the dipole moment between them. 44n vibronic spectroscopy applications, the evolution of the thermofield wavepacket replaces the dynamics of the coherence between two electronic states.Due to the doubling of dimensions, the coherence thermofield dynamics obviously benefits from a combination with semiclassical or other approximate propagation methods that scale favorably with the number of degrees of freedom.Combining the thermofield dynamics with the extended thawed Gaussian approximation [45][46][47] has enabled the simultaneous description of finite-temperature, non-Condon, 40 and even anharmonicity effects. 9,12,41Yet, one would sometimes like to avoid dynamical approximations and obtain the exact quantum solution.
For small systems with a known potential energy surface, the dynamical Fourier method may be used to solve the time-dependent Schrödinger equation on a grid with the split-operator algorithm. 18,22,26,48,49Although limited to low dimensions, the split-operator Fourier method offers a numerically exact solution, which can validate semiclassical and other more approximate approaches.
Here, we combine the coherence thermofield dynamics with the split-operator Fourier method in order to simulate vibronic spectra at finite temperature.This approach is validated on a one-dimensional anharmonic Morse model by comparing the spectra computed with the split-operator coherence thermofield dynamics to the results of the Boltzmann averaging of spectral contributions from different initial vibrational levels.We also examine the thermofield wavepacket's evolution in the augmented configuration space consisting of both the physi-cal and fictitious degrees of freedom.

II. THEORY
Let us consider a molecular system with D nuclear degrees of freedom and two uncoupled adiabatic electronic states: the ground state |g⟩ and the excited state |e⟩.The absorption spectrum for vibronic transitions from |g⟩ to |e⟩ is obtained as the Fourier transform, of the dipole time correlation function, where μ represents the transition dipole moment operator, Ĥg and Ĥe are the vibrational Hamiltonians associated with the ground and excited electronic states, ρ is the vibrational density operator of the ground electronic state, and const = 2πω/(ℏc). 17,22As we are interested only in the dependence of the spectrum on frequency and temperature, in the following we set const = 1/(2 * π).
At a finite temperature T , the density of occupied vibrational levels in the initial population follows the Boltzmann distribution with the inverse temperature β = 1/(k B T ).We can then rewrite C(t) as p n e iωn,gt ⟨n, g|μ † e −i Ĥet/ℏ μ|n, g⟩, where p n = e −βℏωn,g / k e −βℏω k,g is the Boltzmann probability, |n, g⟩ denotes the n-th vibrational eigenfunction of Ĥg , and ℏω n,g is the corresponding vibrational energy. 11The thermal spectrum can be obtained from by propagating the initial vibrational wavepackets |φ n (0)⟩ = μ|n, g⟩ individually as |φ n (t)⟩ = e −i Ĥet/ℏ |φ n (0)⟩.][52] Using the coherence thermofield dynamics, the thermal ensemble is instead mapped to a pure state in a doubled configuration space.The thermal dipole time correlation function (2) can be rewritten exactly 12 (see the Appendix) as a wavepacket autocorrelation function, of an initial thermofield wavepacket, represented in the augmented thermofield space (denoted by a bar) as where q = (q, q ′ ) T is a 2D-dimensional coordinate vector in the augmented configuration space and |q ′ ⟩ represents the position states in the "fictitious" Hilbert space (denoted by a prime).The wavepacket is propagated according to the time-dependent Schrödinger equation with the "augmented" Hamiltonian Ĥ expressed in the position representation as with a 2D × 2D mass matrix, replacing the D × D mass matrix m.In contrast to standard thermofield dynamics, where the two Hamiltonians in the right-hand side of Eq. ( 8) are identical, the thermofield Hamiltonian (8) in vibronic spectroscopy depends on the two Hamiltonians associated with the two electronic states involved in the transition.In essence, the thermofield wavepacket encodes the coherence between the dynamics on the two potential energy surfaces.
If the ground surface is harmonic, the initial thermofield wavepacket 12 ψ(q) = e (i/ℏ)[(q−q0) T • Ā0•(q−q0)/2+ pT is a Gaussian centered at in the phase space.This Gaussian has a temperaturedependent width matrix, 12 where A and B are the with serves to ensure normalization.The spectrum at a finite temperature can now be obtained by applying standard techniques for solving the time-dependent Schrödinger equation to the thermofield wavepacket with a doubled number of coordinates.Although the increase in dimensionality implies that an exact thermofield solution does not save computational time compared to the direct propagation of the density matrix, the wavepacket approach permits the use of zerotemperature numerical methods and avoids implementing methods for solving the von Neumann equation.
Among various numerical methods for propagating wavepackets, we choose to combine the thermofield dynamics with the second-order split-operator algorithm, 18,22,26,48,49 a popular choice for separable Hamiltonians, because it is explicit, numerically stable, and easy to implement and because it preserves several geometric properties of the exact evolution operator, such as unitarity, symplecticity, and time reversibility.

III. NUMERICAL EXAMPLES
In the following, we use natural units with m = ℏ = k B = μ = 1.As in Ref. 12, the excited-state surface is a one-dimensional Morse potential, with the equilibrium position q ref = 1.5, the minimum energy V e (q ref ) = 10, the fundamental frequency ω e = 0.9 at q ref , and the anharmonicity measure χ = 0.02.The initial, ground-state potential energy surface is harmonic and given by Eq. ( 10) with q ′ ref = 0 and κ = 1.The spectra are computed either by Boltzmann averaging [Eq.( 4)] or with thermofield dynamics [Eq.( 6)] at scaled temperatures T ω = k B T /(ℏω g ) = 0, 0.5, 1, and 2, where ω g = κ/m is the ground-state vibrational frequency.Here, as a demonstration, we apply the standard second-order TVT algorithm on a grid with 256 points for each degree of freedom.The wavepacket is propagated for a total time of 1000 with a time step of 0.1.A Gaussian damping function with a half-width at halfmaximum of 15 is applied to the autocorrelation functions before the spectra are evaluated with Eq. ( 1).
As shown in Fig. 1, the thermofield split-operator Fourier method gives the same results as the Boltzmannweighted average.As the temperature increases, the spectral range broadens and additional peaks (the "hot bands") appear.In Ref. 12, a comparison was made between the thermal spectra obtained using the thermofield thawed Gaussian approximation and the "exact" spectra obtained by computing individual Franck-Condon factors by numerical integration.With the time-dependent splitoperator approach, we avoid the need to preselect the relevant transitions and can give numerically exact results for arbitrary global potentials, even if their vibrational eigenfunctions do not have analytical expressions.This approach thus has the potential to validate semiclassical methods in a wider range of situations, including those with a fitted global potential energy surface.FIG. 3. Initial thermofield wavepackets at different scaled temperatures Tω; q is the physical degree of freedom, whereas q ′ is fictitious.
In Fig. 2, we show the decomposition of the thermal spectrum at T ω = 1 into the contributions from different initial vibrational levels.Only spectra from levels n = 0, 1, 2, and 3 shown since the contributions from n > 3 are negligible here (the Boltzmann weight is 0.0002 for n = 4).As expected, the higher vibrational levels are the source of hot bands and expand the spectral range.Whereas the Boltzmann averaging requires propagating each initial vibrational level separately, a single thermofield wavepacket trajectory takes into account all thermal contributions.
With a grid-based method, we can easily visualize the wavepacket and its evolution.Figure 3 shows the initial thermofield wavepacket at different temperatures in the augmented coordinate space.As the temperature increases, the off-diagonal block B of the initial width matrix (13) increases and the initial thermofield wavepacket rotates and stretches as a result.
In Fig. 4, we show the evolution of the wavepacket at T ω = 1.Whereas the extent of the wavepacket over the fictitious degree of freedom q ′ remains approximately constant, the initial excitation of the coordinate q on the anharmonic excited-state surface leads to a significant evolution of the wavepacket in the physical dimension.

IV. DISCUSSION AND CONCLUSION
We have described a very simple method for computing vibronic spectra at finite temperatures.This method applies the Fourier split-operator propagation algorithm to the thermofield wavepacket representation of the electronic coherence, avoiding both Boltzmann averaging and solving the von Neumann equation.Instead, we solve in the augmented configuration space at different times t; q is the physical degree of freedom, whereas q ′ is fictitious.
the zero-temperature Schrödinger equation in an augmented space.The method yields quantum mechanically exact spectra of arbitrary global potentials, making it a valuable tool for validating semiclassical and other approximate techniques for finite-temperature spectra calculations.In contrast to Boltzmann averaging spectral contributions from different initial vibrational levels, the thermofield dynamics approach obtains the exact result for a given temperature from a single propagation.At high temperatures, as the number of contributing states increases, the thermofield approach becomes more favorable.][55][56] The method presented here is based on the full tensorproduct grid and, therefore, is limited to systems with very few degrees of freedom.To extend its applicability to higher-dimensional systems, one can combine the split-operator thermofield dynamics with adaptive moving grids, [57][58][59] matching-pursuit coherent-state basis, 60 tensor-train representation, 61 or any other technique designed for the zero-temperature case.In fact, tensor-train representations have already been used extensively 32,36 for the thermofield dynamics, although not in conjunction with the split-operator Fourier method.Here, we have intentionally restrained ourselves to the simplest, full-grid method, because we did not want to obscure the main idea by combining it with various strategies for increasing the accessible dimensionality.We believe that even the full-grid version can be interesting for researchers who would like to quickly check the thermal effects on their zero-temperature spectra of lowdimensional systems without implementing new codes.
The method is not restricted to linear spectra.In the same way as the Gaussian, 9 multiconfigurational Ehrenfest, 42 or optimized mean-trajectory 43 semiclassical thermofield dynamics, the method can be applied to more sophisticated experimental techniques, such as pump-probe and two-dimensional spectroscopies.Beyond vibronic spectra, this approach may also find applications in other finite-temperature dynamics problems, such as the calculation of the rates of internal conversion. 27,37,62,63

FIG. 4 .
FIG.4.Evolution of the thermofield wavepacket (for Tω = 1) in the augmented configuration space at different times t; q is the physical degree of freedom, whereas q ′ is fictitious.