Quantum system-bath dynamics with quantum superposition sampling and coupled generalized coherent states

Previously, we introduced two versions of the Multiconfigurational Ehrenfest (MCE) approach to high dimensional quantum dynamics. It has been shown that the first version, MCEv1, converges well to the existing benchmarks for high dimensional model systems. At the same time, it was found that the second version, MCEv2, had more difficulty converging in some regimes. As MCEv2 is particularly suited for direct dynamics, it is important to facilitate its convergence. This paper investigates an efficient method of basis set sampling, called Quantum Superposition Sampling (QSS), which dramatically improves the performance of the MCEv2 approach. QSS is tested on the spin-boson model, often used for modeling of open quantum systems. It is also shown that the quantum subsystem in the spin-boson model can be conveniently treated with the help of two level system coherent states. Generalized coherent states, which combine two level system coherent states for the description of the system and Gaussian coherent states for description of the bath, are introduced. Various forms of quantum equations of motion in the basis of generalized coherent states can be developed by analogy with known quantum dynamics equations in the basis of Gaussian coherent states; in par-ticular, the multiconfigurational Ehrenfest method becomes a version of coupled generalized coherent and can generalization of sampling method known for the existing coupled coherent states method uses Gaussian coherent results have been obtained with MCEv1 for the spin-boson model a broad range of parameters and for a system of qubits coupled with environment. shown classical” cases of spin-boson model.


I. INTRODUCTION
The most straightforward way of simulating open systems would be to model them by a small quantum subsystem coupled with a large number of bath modes describing the environment and to solve the time dependent Schrödinger equation for all degrees of freedom, both those of the main system and the bath. The spinboson model, comprised of a two level system (or perhaps several coupled two level systems) and a bath of harmonic modes, is a popular example of a system-bath problem. Although much progress has been made in the development of methods of high dimensional quantum mechanics, treating large numbers of bath modes exactly still remains a difficult task. For example, hierarchical methods, such as Hierarchy of Equations of Motion (HEOM) [1][2][3] and Hierarchy of Pure States (HOPS), 4,5 allow simulations of system bath problems, but they are efficient only for certain spectral densities of the bath. Time evolving density matrix using Orthogonal Polynomials algorithm (TEDOPA) 6,7 is efficient in the low temperature/short time regime, 8 although some efforts have been made to extend it to longer times. 8 Multiconfigurational time dependent Hartree (MCTDH) 9 is a basis set technique which can treat quantum system-bath problems with a high level of accuracy. In MCTDH, the bath is discretized and a large number of bath modes are combined in groups so that a small number of very flexible time dependent basis functions per group are sufficient for adequate description of the dynamics. This mode combination is crucial for MCTDH, and hence, if the mode combination is not right, MCTDH can be hard to converge. 10 In recent years, we have developed the Multiconfigurational Ehrenfest (MCE) approach, [11][12][13] which uses trajectory guided basis sets of Gaussian coherent states to perform formally exact multidimensional quantum simulations. Several other techniques have been developed, [14][15][16][17][18] all using basis sets of Gaussian coherent states which follow the quantum wave function, thus economizing the basis set size, with these methods only differing in the way they guide the trajectories of the CS basis set. Among these techniques are variational Multiconfigurational Gaussians (vMCG) 15,16 and Davydov Ansatz 17,18 (DA), both of which rely on accurate but complicated and often numerically unstable variational trajectories. The number of degrees of freedom, which can be treated by vMCG and DA, is limited, and only the cases of spin-boson model for which the bath can be discretized with a limited number of modes can be treated. Multiple Spawning 14 (MS) in which the trajectories of Gaussians are purely classical and therefore are not always well suited for quantum dynamics, is another member of the same family of methods. MCE sits between these techniques using Ehrenfest trajectories to guide the basis. Ehrenfest trajectories guided by quantum averaged Hamiltonian are nonclassical enough to incorporate some quantum effects and to reach classically forbidden regions, but at the same time they are simple and stable. Technicalities of the MCE method, however, such as basis set initial conditions sampling, turn out to be extremely important for the convergence of MCE calculations.
Two versions of MCE, referred to as MCEv1 and MCEv2, have been developed. Both methods represent the propagating wave function as a superposition of Ehrenfest configurations guided by Ehrenfest trajectories, but MCEv1 and MCEv2 differ in the way quantum coupling between the configurations works. As a result of this difference, the trajectories which guide Ehrenfest configurations differ in the two versions of MCE. In MCEv1, the trajectories are not independent and they "push" each other so that the whole ensemble of such trajectories has to be run simultaneously. In MCEv2, the trajectories are independent and can be run one by one if desired, which makes MCEv2 suitable for direct dynamics. Both MCEv1 and MCEv2 have been tested on the spin-boson model. MCEv1 was found to converge well for a broad range of parameters of the model; 11 at the same time MCEv2 was found to be much harder to converge for some cases of the spin-boson model and special methods improving the convergence had to be developed. 19 Given the use of MCEv2 in direct dynamics simulations of nonadiabatic photochemical reactions (see, for example, our recent works 20 and 21), it is important to find new ways to facilitate the convergence of MCEv2.
In this paper, we show that a very simple sampling technique, termed Quantum Superposition Sampling (QSS), greatly improves convergence for MCEv2. We also show how this technique arises naturally if the quantum subsystem is treated with a basis of so called SU(2) coherent states so that MCEv2 can be viewed as a version of the recently introduced Coupled Generalized Coherent States (CGCS) approach. 22 In the supplementary material, we also derive various forms of equations for treating quantum system-bath models with CGCS.

II. THEORY
A. Two versions of the multiconfigurational Ehrenfest method: MCEv1 and MCEv2 Two versions of the multiconfigurational Ehrenfest method have been developed, dubbed MCEv1 [3] and MCEv2 [4] . The main building block of both versions is the Ehrenfest configuration which for a two level system takes the form where the Gaussian coherent state |z bth (t)⟩ = |p bth (t), q bth (t)⟩ describes the bath coupled with a superposition of the system states |1⟩ = |ϕ sstm 1 ⟩ and |0⟩ = |ϕ sstm 0 ⟩. These system states could denote electron donor and acceptor states or electronic excited and ground states, for example, and the bath often represents the vibrations of nuclei or phonons.
For simplicity, we assume that the quantum subsystem is described by a single two-level system, but generalization to a more generic case of many coupled two level systems or systems with more than two levels is not difficult. The multidimensional Gaussian coherent state |z(t)⟩ = |p(t), q(t)⟩ is a product of 1D coherent states, Each 1D CS is a Gaussian wavepacket and the center of the CS wavepacket (3) is described by a single com- relating to a point (q, p) in phase space of the bath. Both MCEv1 and MCEv2 use similar wave function Ansätze with a combination of Ehrenfest configurations, taking the form for MCEv1 and for MCEv2. The difference is in the additional coefficient A (k) (t) in the MCEv2 Ansatz (5), which is used to provide the coupling between the Ehrenfest configurations |ψ (k) (t)⟩. In the MCEv2 Ansatz (5), configurations are normalized, to guide the bath basis set, which are equivalent to Hamilton's equations of motion in a coherent states context. In Eq. (7), is the Ehrenfest average Hamiltonian for an individual configuration. In MCEv1, the amplitudes a 1 (k) (t) and a 0 (k) (t) are all coupled with each other [see Eq. (28) in Ref. 11 and supplementary material to this article]. This is in contrast to MCEv2 where Ehrenfest configurations are independent from each other, and hence for MCEv2, the amplitudes a 1 (k) (t) and a 0 (k) (t) of quantum states |1⟩ and |0⟩ are coupled only within an individual configuration, Equations (7)-(9) are simply those for the standard Ehrenfest trajectory, which yields the best wave function described by a single Ehrenfest configuration (see Sec. 3 of the supplementary material).
To provide the coupling between configurations in MCEv2, there is additional set of coupled equations for the amplitudes A (k) (t) derived from the time dependent Schrödinger equation (see Ref. 12 and Sec. 4 of the supplementary material for more details). The are simply those for the amplitudes of time dependent basis functions |ψ (k) (t)⟩ in the wave function Ansatz (5). The time dependence of |ψ (k) ⟩ is given by Eqs. (7)- (9). The coupling between amplitudes in MCEv2 is illustrated in Fig. 1 where it is also compared with coupling between the amplitudes in the MCEv1 approach. As one can see from the diagram as well as from Eq. (9) in MCEv2, the amplitudes a 1 (k) (t) and a 0 (k) (t) are only coupled within configurations, but not across them.
MCEv2 was developed for use with direct dynamics, where running independent trajectories is more convenient. An enhancement of MCEv2, which models wavepacket splitting better, referred to as Ab initio Multiple Cloning (AIMC), has been successfully used in direct dynamics calculations to simulate ultrafast photodynamics 20,21 and energy transfer in conjugated molecules "on the fly." 23,24 In principle, both MCEv1 and MCEv2 are exact quantum techniques, and must converge to the exact quantum result with a sufficiently large number of basis set configurations, but the convergence properties are different for the two techniques. Well converged   11 with a broad range of parameters and for a system of qubits coupled with environment. 25 MCEv2 was shown to converge well for short time dynamics and for "more classical" cases of spin-boson model. However, in some regimes of the spin-boson model, MCEv2 had difficulty converging. 19 Facilitating convergence of MCEv2 is important because it is used in our AIMC direct on the fly dynamics. 23 In AIMC, electronic structure theory is used to calculate potential energy surfaces ab initio and the trajectories of Ehrenfest configurations are run independently from each other. Running Ehrenfest trajectories on the fly is very expensive computationally and may take years of CPU time, but eventually they are saved on the hard drive, providing the time dependent basis set for MCEv2. Then, coupled equations (10) for the amplitudes A (k) (t) are solved as "postprocessing." To improve the convergence of MCEv2, special basis set sampling methods were developed. They were tested on the spinboson problem and used in AIMC. 19,23 These sampling methods, however, are not trivial, and finding new and better techniques to improve convergence of MCEv2 is of great importance for the future use in AIMC direct dynamics. The spin-boson model provides a very convenient testing ground, and in this paper, we use it to test a simple, yet very effective, sampling technique which greatly improves convergence on MCEv2.

B. Enhancing initial sampling with quantum superposition sampling (QSS)
Seemingly, the most natural way of choosing basis Ehrenfest configurations is to define the amplitudes as being 1 for the populated electronic state and zero for the unoccupied electronic state |ψ (k) (t)⟩=(1|1⟩ + 0|0⟩)|z (k) (t)⟩, so that the initial MCEv2 wave function becomes (11) where for all configurations, the initial Ehrenfest amplitudes a 1 (k) (0) and a 0 (k) (0) are the same. In most of our previous calculations, the initial bath CSs z k (0) were selected randomly using sampling techniques suggested in Ref. 26. It is, however, not necessary to set the initial configuration amplitudes in (11) entirely on one state; only the wavefunction as a whole must have initial population wholly on the initial state |1⟩. QSS differs in its approach to selecting initial conditions of the Ehrenfest configurations by taking randomly not only z (k) (0) but also a 1 (k) (0) and a 0 (k) (0) so that where zinit is the initial condition for the propagating wavepacket, which in Eq. (12) is assumed to be a coherent state itself. The wave function (12) can then be propagated in exactly the same way as before. This is somewhat counterintuitive as each basis Ehrenfest configuration includes a contribution of the unpopulated state.
In practice, the construction of this wavefunction can be carried out by generating the Gaussian CSs |z k (0)⟩ from a compressed normal distribution around the initial wavefunction coherent state |zinit⟩ as was done in previous applications of MCE, but then instead of defining the amplitudes a 1 (k) (0) and a 0 (k) (0) to be equal to 1 and 0, a pair of random numbers, φ and Δφ are generated between the limits of 0 and 2π. These random numbers can then be used to generate the amplitudes such that Randomly selecting the single configuration amplitudes has a necessary effect on the calculations to find the initial values for the interconfiguration amplitude A (k) , and so this process must be modified from the process used previously for MCEv2. With QSS, the calculation for the initial A (k) amplitudes uses the identity operator for a basis set of nonorthogonal basis functions which includes the elements Ω (ψψ)−1 kl of the matrix Ω (ψψ )−1 inverse of the overlap matrix Ω (ψψ ) with the elements Through this identity operator, the amplitudes A (k) can then be given as where |Ψ⟩ = (1|1⟩ + 0|0⟩)|zinit⟩. In practice, a set of linear equations C = Ω (ψψ ) A for the vector of the amplitudes A (k) is solved instead of inverting the matrix Ω (ψψ ) . Here, C is the vector with the components Through this, the amplitudes A (k) can be calculated in such a way that the population of the resulting wave function will be entirely on the state |1⟩ without need to force the single configuration amplitudes to be wholly on this state. This sampling, while seemingly counterintuitive, is similar in spirit to techniques that have been used previously in the Coupled Coherent States (CCS) method. 24,26 It is also well known that convergence of the Gaussian CSs based methods of nonadiabatic dynamics improves if some basis Gaussian coherent states in (1) are put on the unpopulated electronic state with zero amplitude A (see Refs. 12 and 13, for instance). The idea of randomness has also been used in recent MCEv2 calculations carried out by the group of Burghardt. 27 Here, we will look at random sampling of the quantum subsystem more systematically. An important consideration is that using this sampling technique does not add any significant computational cost to the simulation but, as we will show below, quantum superposition sampling is very useful and significantly facilitates convergence of MCEv2.

III. TRAINS AND CLONING
Previously, we found that two sampling tricks improved the performance of MCEv2. First, it was shown that trains or time displaced basis sets of CSs of Ehrenfest configurations, which follow the same guiding trajectory but with time delay, enhance the convergence. A second very important sampling technique in MCEv2 and related methods is cloning, which allows resampling of the basis set during propagation. 19 The idea of cloning is to create adaptive basis sets by splitting Ehrenfest configurations into two after passing through regions of nonadiabatic coupling so as to avoid the possibility of configurations propagating along an average potential energy surface which is unphysical in nature. Currently, in MCEv2 an Ehrenfest configuration can be "cloned" to the states |1⟩ and |0⟩ as follows: Cloning of a configuration is performed when two criteria are met: first when the population of both states is not small, which is regulated by the cloning threshold |a 1 (k) a 0 (k) |, and second when the gradients of their potential energies differ significantly (see Ref. 19). This way of cloning is largely inspired by spawning in the multiple spawning approach 14,28 to quantum nonadiabatic dynamics. In Sec. IV, we will show that using together with QSS trains and cloning improves the performance of MCEv2 even further.

IV. TEST OF THE QSS
Following the work of Symonds et al., 19 here we also use the spin-boson model to provide the first test of the QSS sampling techniques in MCEv2. In the introductory work for the MCE method, 11 MCEv1 was found to converge well in all major regimes of the spin-boson model, for which MCTDH benchmark data existed that time. At the same time, MCEv2 had difficulty converging in more "quantum" cases. In the previous studies of MCE with the spinboson model, 11,19 an Ohmic spectral density with an exponential cutoff was used, and this spectral density is also used in this paper. The discretization scheme was exactly the same as in Refs. 11 and 19 where more details can be found. Here, we will only mention that the parameters of the model ωc, α k , Δ, ε and β = 1 T describe the shape of the spectral density of the bath (ωc, α k ), the coupling between two states of system (Δ), energy shift between two quantum states of the system (ε), and the temperature. Following the previous work, 11,19,29 the units are chosen such that the energy is scaled to the parameter Δ. In addition, similarly to the previous work, a thermally averaged population difference between the two states of the system was calculated. A set Nrpt of bath initial conditions |zinit⟩ is generated from quantum Boltzmann distribution and initial wave function |Ψ⟩ = |1⟩|zinit⟩ is propagated Nrpt times using the MCEv2 Ansatz (5) with quantum superposition sampling of the initial coefficients a 1 (k) (0) and a 0 (k) (0) until the average population difference converges.
We do not consider the cases of the spin-boson model for which MCEv2 converged previously, and instead focus on the cases where MCEv2 had more difficulty converging in order to demonstrate the efficacy of the QSS method. The initial test case considered for QSS was that of a pair of asymmetric wells, drawing on the previous work 14 where a large disparity between the results of MCEv2 and benchmark has been demonstrated for this case. 14 At the same time, MCEv1 achieved a very good agreement with MCTDH. 11 Figure 2 shows the original MCEv2 result for N bf = 200 basis Ehrenfest configurations which without QSS was quite far away from MCTDH benchmark. Increasing the basis set size in MCEv2 to 1000 led to almost no improvement. A direct comparison between the population difference calculated using MCEv2 with the original sampling methodology and MCEv2 with quantum superposition sampling immediately demonstrates the improved convergence of this enhancement. Figure 2 shows that with a basis set as small as 100 basis functions, there is a noticeable improvement in convergence towards the benchmark MCTDH 29 or MCEv1 11 results. When the number of basis functions is increased to 1000, describing a larger area of phase space, the agreement is almost perfect within the first few oscillations and significantly closer afterward. Consequently, increasing the number of basis functions alone is clearly not sufficient for complete convergence. Following on the success of MCEv2 with cloning, 14 cloning was applied to the simulation of MCEv2 with QSS. Because of the relatively good agreement between the MCTDH benchmark and MCEv2 using QSS and a relatively small basis sets, it was now possible to increase the maximum The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp number of clones from four as used by Symonds et al. 19 to 14 without risking inflating the basis set to an unmanageable size. It can be seen in Fig. 3 that with a relatively small initial basis set and a large number of cloning events, convergence with the benchmark data is greatly improved from when no cloning is used, with a deviation between the two population differences visible only after the fourth oscillation.
The result for a model decoherent system is shown in Fig. 4. When used on its own, MCEv2 gives oscillations of the decaying population difference, which are much too large in magnitude in comparison to the benchmark result. Again, an increase of the basis set size does not produce any visible improvement of MCEv2 without QSS. The use of QSS dampens these oscillations so they align more closely with the MCTDH benchmark data, which can be seen in Fig. 4.
We have also considered the case of localized wavefunctions, where the wavefunction remains in the initial electronic state. Unlike MCEv1, which was able to reproduce benchmark successfully, localization proved to be a very difficult case for the MCEv2 method in the previous simulations. 14 These wavefunctions are characterized by very high dimensionality (i.e., very many bath modes are needed to discretize the bath properly), high cut-off frequencies of the bath spectrum, strong system/bath coupling, and very low temperatures, making them highly nontrivial to simulate. The MCEv2 result displays high frequency oscillations that are not seen in the benchmark MCTDH 30 and MCEv1 11 data. While the wavefunction remains localized and wholly on the initial electronic state, the average population difference is higher than is seen in the benchmark data. As can be seen in Fig. 5, however, quantum superposition sampling vastly improves the agreement with the MCTDH benchmark and lowers the average population difference to closer to the correct level although the unwanted oscillations still exist, albeit in a much dampened form.
To improve the performance of MCEv2 with QSS further, we employed train basis sets, with encouraging result shown in Fig. 6.
Increasing the number of trains and their length and at the same time keeping the distance between the "carriages" of the train as small as possible systematically improves the agreement with the benchmark. The case of tunneling between symmetric wells at low temperature with reasonably strong system/bath coupling also converged FIG. 7. Comparison of the population differences for tunneling between a pair of low temperature symmetric wells with fairly strong system/bath coupling using MCEv2 with and without QSS, and these are then compared to a MCTDH benchmark calculation (dashed-dotted). 29 The parameters of the spin-boson model are ωc/∆ = 7.5, α K = 0.6, β∆ = 5.0, and ε/∆ = 0 with M = 60 frequencies used to discretize the bath. For the simulation with no QSS used N bf = 200 basis functions (dashed), and the simulations with QSS used N bf = 200 (short dash) and N bf = 1000 (solid) basis functions. The number of repetitions is N rpt = 128 for all simulations.

FIG. 8.
Comparison of the population differences for tunneling between a pair of low temperature symmetric wells with fairly strong system/bath coupling for MCEv2 (dashed) and MCEv2 using QSS and swarm clones (solid), compared to a MCTDH benchmark calculation (dashed-dotted). 29 The parameters of the spin-boson model are ωc/∆ = 7.5, α K = 0.6, β∆ = 5.0, and ε/∆ = 0 with M = 60 frequencies used to discretize the bath. Both simulations used or initially had N bf = 200 basis functions, with the QSS simulation allowing N cln = 14 cloning events per basis function with cloning threshold |a 1k a 0k | > 0.1. N rpt = 128 repetitions for both simulations.
with MCEv1 11 but proved to be a limiting case for MCEv2, 19 with a "blind cloning" process needed to place empty configurations on the acceptor state on the other side of the potential barrier. QSS prevents the unwanted localization in the initial electronic state that occurs with MCEv2 without QSS, by randomly placing some basis functions in the other electronic state. This can be seen in Fig. 7, where quantum superposition sampling dramatically decreases the disagreement with the benchmark data. As expected, increasing the size of the basis set improves convergence, as shown in Fig. 7. It is clearly necessary, however, to use other techniques to encourage greater agreement to the benchmark. Figure 8 demonstrates how cloning can be used to fit the simulation closer to the MCTDH, the result being a simulation that decays to the correct magnitude even if the rate of the decay has not fully converged. Like in the localized case, QSS places basis functions in the classically disallowed electronic state, which then allows the simulation to better model the tunneling between the two states.

V. MCEv2 AS A VERSION OF COUPLED GENERALISED COHERENT STATES
The idea to use random initial conditions for the amplitudes a 1 (k) (0) and a 0 (k) (0) in the Ehrenfest configurations seems counterintuitive, but it leads to a big improvement in the convergence of MCEv2 calculations. The idea becomes very natural, however, if one thinks of the subsystem [a 1 (t)|1⟩ + a 0 (t)|0⟩] of the Ehrenfest configuration (1) as a two level system coherent state, also called SU(2) coherent state 31 which can be written as |ζ(a 1 , a 0 )⟩ = a 1 |1⟩ + a 0 |0⟩.
The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp Usually, the SU(2) coherent state is written as where one of the coefficients is assumed to be real. Having both a 1 and a 0 as complex would introduce an insignificant phase factor, but the advantage is that now the amplitudes are treated on an equal footing. A multidimensional SU(2)-CS representing M two-level systems as simply the Hartree product of several 1D CSs, In these notations, the MCEv2 Ansatz becomes where the wave function is a superposition of generalized coherent states |ζ (n) (t), z (n) (t)⟩, a mixture of SU(2) |ζ (n) (t)⟩ and Gaussian |z (n) (t)⟩ CSs. While the difference between Eq. (22) and the MCEv2 Ansatz (5) is only in the notations, the new notations help to rationalize the new QSS sampling in the context of established theory. Recently, it has been shown that the machinery developed for Gaussian CSs in the Coupled Coherent States (CCS) method 24 can be adopted for the use with generalized CSs. 22 From this prospective, MCEv2 is equivalent to Coupled Generalized Coherent States (CGCS). The details of working equations of Coupled Generalized Coherent States (CGCS) are given in the supplementary material. Sampling techniques developed for standard coupled coherent states in Refs. 24 and 26 can therefore be used in the CGCS approach. One of the ideas suggested in Ref. 26 for the CCS approach was to sample the most important and the most quantum degrees of freedom of Gaussian coherent states from a broad and random distribution and to sample less important degrees of freedom from a random but very narrow distribution, creating a "pancake" in a multidimensional phase space. 26 When applied to CGCS, this idea of a "pancake" distribution gives the grounds for quantum superposition sampling if the degrees of freedom associated with two level system are regarded as the most important ones.
Thus, the idea of coupled generalized coherent states suggested in Ref. 22 can be extended and numerous methods of quantum propagation developed for the basis sets of Gaussian coherent states can also be used for generalized CSs. In addition, various forms of quantum coupled equations suggested previously for Gaussian coherent states can be written for the generalized coherent states. In the supplementary material, we use the language of generalized coherent states to derive and to generalize the equations of MCEv1 and MCEv2 as well as fully variational generalized coherent states.

VI. SUMMARY AND CONCLUSIONS
In this paper, we demonstrated that the new quantum superposition sampling significantly improves the performance of the MCEv2 method. We have tested QSS on several cases of the spin-boson model where MCEv2 failed previously and observed a significant improvement. Even when QSS still disagrees with the benchmark, it brings the result closer to it and therefore represents a better starting point for using other sampling methods, such as cloning and trains. QSS is a natural extension of important sampling technique known in the CCS method to the coupled generalized coherent states approach. It greatly improves the efficiency of MCEv2, and now, the hope is that it should be useful in the future for simulations of system-bath quantum problems and in direct dynamics "on-the-fly" simulations.

SUPPLEMENTARY MATERIAL
The supplementary material provides more details on the definition of generalized coherent states and derivations of various types of quantum equations of motions for the wave functions in the basis of GCSs.