Variational dynamics of the sub-Ohmic spin-boson model on the basis of multiple Davydov D1 states

Efficient simulation of open quantum systems coupled to a reservoir through multiple channels

It is challenging to simulate open quantum systems that are connected to a reservoir through multiple channels. For example, vibrations may induce fluctuations in both energy gaps and electronic couplings, which represent two independent channels of system-bath couplings. Systems of this kind are ubiquitous in the processes of excited state radiationless decay. Combined with density matrix renormalization group (DMRG) and matrix product states (MPS) methods, we develop an interaction-picture chain mapping strategy for vibrational reservoirs to simulate the dynamics of these open systems, resulting in time-dependent spatially local system-bath couplings in the chain-mapped Hamiltonian. This transformation causes the entanglement generated by the system-bath interactions to be restricted within a narrow frequency window of vibrational modes, enabling efficient DMRG/MPS dynamical simulations. We demonstrate the utility of this approach by simulating singlet fission dynamics using a generalized spin-boson Hamiltonian with both diagonal and off-diagonal system-bath couplings. This approach generalizes an earlier interaction-picture chain mapping scheme, allowing for efficient and exact simulation of systems with multi-channel system-bath couplings using matrix product states, which may further our understanding of nonlocal exciton-phonon couplings in exciton transport and the non-Condon effect in energy and electron transfer.

Unifying methods for optimal control in non-Markovian quantum systems via process tensors

The large dimensionality of environments is the limiting factor in applying optimal control to open quantum systems beyond the Markovian approximation. Various methods exist to simulate non-Markovian systems, which effectively reduce the environment to a number of active degrees of freedom. Here, we show that several of these methods can be expressed in terms of a process tensor in the form of a matrix-product-operator, which serves as a unifying framework to show how they can be used in optimal control and to compare their performance. The matrix-product-operator form provides a general scheme for computing gradients using back propagation and allows the efficiency of the different methods to be compared via the bond dimensions of their respective process tensors.

Finite temperature dynamics in a polarized sub-Ohmic heat bath: A hierarchical equations of motion-tensor train study

The dynamics of the sub-Ohmic spin-boson model under polarized initial conditions at finite temperatures is investigated by employing both analytical tools and the numerically accurate hierarchical equations of motion-tensor train method. By analyzing the features of nonequilibrium dynamics, we discovered a bifurcation phenomenon, which separates two regimes of the dynamics. It is found that before the bifurcation time, increasing temperature slows down the population dynamics, while the opposite effect occurs after the bifurcation time. The dynamics is highly sensitive to both initial preparation of the bath and thermal effects.

Fast Quantum State Preparation and Bath Dynamics Using Non-Gaussian Variational Ansatz and Quantum Optimal Control

Fast preparation of quantum many-body states is essential for myriad quantum algorithms and metrological applications. Here, we develop a new pathway for fast, nonadiabatic preparation of quantum many-body states that combines quantum optimal control with a variational Ansatz based on non-Gaussian states. We demonstrate our approach on the spin-boson model, a single spin interacting with the bath. We use a multipolaron Ansatz to prepare near-critical ground states. For one mode, we achieve a reduction in infidelity of up to ≈60 (≈10) times compared to linear (optimized local adiabatic) ramps; for many modes, we achieve a reduction in infidelity of up to ≈5 times compared to nonadiabatic linear ramps. Further, we show that the typical control quantity, the leakage from the variational manifold, provides only a loose bound on the state's fidelity. Instead, in analogy to the bond dimension of matrix product states, we suggest a controlled convergence criterion based on the number of polarons. Finally, motivated by the possibility of realizations in trapped ions, we study the dynamics of a system with bath properties going beyond the paradigm of (sub- and/or super-) Ohmic couplings.

Dynamics of the spin-boson model: The effect of bath initial conditions

The dynamics of the (sub-)Ohmic spin-boson model under various bath initial conditions is investigated by employing the Dirac-Frenkel time-dependent variational approach with the multiple Davydov D₁ Ansatz in the interaction picture. The validity of our approach is carefully checked by comparing the results with those of the hierarchy equations of motion method. By analyzing the features of nonequilibrium dynamics, we identify the phase diagrams for different bath initial conditions. We find that for the spectral exponent s < s_c, there exists a transition from coherent to quasicoherent dynamics with increasing coupling strengths. For s_c < s ≤ 1, the coherent to incoherent crossover occurs at a certain coupling strength and the quasicoherent dynamics emerges at much larger couplings. The initial preparation of the bath has a considerable influence on the dynamics.

The hierarchy of Davydov's Ansätze: From guesswork to numerically "exact" many-body wave functions

This Perspective presents an overview of the development of the hierarchy of Davydov's Ansätze and a few of their applications in many-body problems in computational chemical physics. Davydov's solitons originated in the investigation of vibrational energy transport in proteins in the 1970s. Momentum-space projection of these solitary waves turned up to be accurate variational ground-state wave functions for the extended Holstein molecular crystal model, lending unambiguous evidence to the absence of formal quantum phase transitions in Holstein systems. The multiple Davydov Ansätze have been proposed, with increasing Ansatz multiplicity, as incremental improvements of their single-Ansatz parents. For a given Hamiltonian, the time-dependent variational formalism is utilized to extract accurate dynamic and spectroscopic properties using Davydov's Ansätze as its trial states. A quantity proven to disappear for large multiplicities, the Ansatz relative deviation is introduced to quantify how closely the Schrödinger equation is obeyed. Three finite-temperature extensions to the time-dependent variation scheme are elaborated, i.e., the Monte Carlo importance sampling, the method of thermofield dynamics, and the method of displaced number states. To demonstrate the versatility of the methodology, this Perspective provides applications of Davydov's Ansätze to the generalized Holstein Hamiltonian, variants of the spin-boson model, and systems of cavity-assisted singlet fission, where accurate dynamic and spectroscopic properties of the many-body systems are given by the Davydov trial states.

Inspecting molecular aggregate quadratic vibronic coupling effects using squeezed coherent states

We present a systematic comparison of three quantum mechanical approaches describing excitation dynamics in molecular complexes using the time-dependent variational principle (TDVP) with increasing sophistication trial wavefunctions (ansatze): Davydov D₂, squeezed D₂ (sqD₂) and a numerically exact multiple D₂ (mD₂) ansatz in order to characterize validity of the sqD₂ ansatz. Numerical simulations of molecular aggregate absorption and fluorescence spectra with intra- and intermolecular vibrational modes, including quadratic electronic-vibrational (vibronic) coupling term, which is due to vibrational frequency shift upon pigment excitation are presented. Simulated absorption and fluorescence spectra of a J type molecular dimer with high frequency intramolecular vibrational modes obtained with D₂ and sqD₂ ansatze match the spectra of mD₂ ansatz only in the single pigment model without quadratic vibronic coupling. In general, the use of mD₂ ansatz is required to model an accurate dimer and larger aggregate's spectra. For a J dimer aggregate coupled to a low frequency intermolecular phonon bath, absorption and fluorescence spectra are qualitatively similar using all three ansatze. The quadratic vibronic coupling term in both absorption and fluorescence spectra manifests itself as a lineshape peak amplitude redistribution, static frequency shift and an additional shift, which is temperature dependent. Overall the squeezed D₂ model does not result in a considerable improvement of the simulation results compared to the simplest Davydov D₂ approach.

Non-Markovian Vibrational Relaxation Dynamics at Surfaces

Vibrational dynamics of adsorbates near surfaces plays both an important role for applied surface science and as model lab for studying fundamental problems of open quantum systems. We employ a previously developed model for the relaxation of a D-Si-Si bending mode at a D:Si(100)-(2$\times$1) surface, induced by a ``bath' of more than $2000$ phonon modes [U. Lorenz, P. Saalfrank, Chem. Phys. {\bf 482}, 69 (2017)], to extend previous work along various directions. First, we use a Hierarchical Effective Mode (HEM) model [E.W. Fischer, F. Bouakline, M. Werther, P. Saalfrank, J. Chem. Phys. {\bf 153}, 064704 (2020)] to study relaxation of higher excited vibrational states than hitherto done, by solving a high-dimensional system-bath time-dependent Schr\"odinger equation (TDSE). In the HEM approach, (many) real bath modes are replaced by (much less) effective bath modes. Accordingly, we are able to examine scaling laws for vibrational relaxation lifetimes for a realistic surface science problem. Second, we compare the performance of the multilayer multiconfigurational time-dependent Hartree (ML-MCTDH) approach with the recently developed coherent-state based multi-Davydov D2 {\it ansatz} [N. Zhou, Z. Huang, J. Zhu, V. Chernyak, Y. Zhao, {J. Chem. Phys.} {\bf 143}, 014113 (2015)]. Both approaches work well, with some computational advantages for the latter in the presented context. Third, we apply open-system density matrix theory in comparison with basically ``exact' solutions of the multi-mode TDSEs. Specifically, we use an open-system Liouville-von Neumann (LvN) equation treating vibration-phonon coupling as Markovian dissipation in Lindblad form to quantify effects beyond the Born-Markov approximation.

Variational Squeezed Davydov Ansatz for Realistic Chemical Systems with Nonlinear Vibronic Coupling

Chemical systems normally possess strong nonlinear vibronic couplings at both zero and finite temperature. For the lowest-order quadratic couplings, here, we introduce a squeezing operator into a variational coherent-state-based method, Davydov ansatz, to simulate the quantum dynamics and the respective spectroscopy. Two molecular systems, pyrazine and the 2-pyridone dimer, are taken as calculated model systems, both of which involve nontrivial quadratic vibronic couplings in high- and low-frequency regions, respectively. Upon a comparison with the benchmarks, the method manifests its advantage for nonlinear couplings. The squeezed bases are also proven to be applicable for the finite temperature by adapting with the thermofield dynamics.

Simulation of Ab Initio Optical Absorption Spectrum of β-Carotene with Fully Resolved S0 and S2 Vibrational Normal Modes

The electronic absorption spectrum of β-carotene (β-Car) is studied using quantum chemistry and quantum dynamics simulations. Vibrational normal modes were computed in optimized geometries of the electronic ground state S₀ and the optically bright excited S₂ state using the time-dependent density functional theory. By expressing the S₂-state normal modes in terms of the ground-state modes, we find that no one-to-one correspondence between the ground- and excited-state vibrational modes exists. Using the ab initio results, we simulated the β-Car absorption spectrum with all 282 vibrational modes in a model solvent at 300 K using the time-dependent Dirac-Frenkel variational principle and are able to qualitatively reproduce the full absorption line shape. By comparing the 282-mode model with the prominent 2-mode model, widely used to interpret carotenoid experiments, we find that the full 282-mode model better describes the high-frequency progression of carotenoid absorption spectra; hence, vibrational modes become highly mixed during the S₀ → S₂ optical excitation. The obtained results suggest that electronic energy dissipation is mediated by numerous vibrational modes.

Modeling Molecular J and H Aggregates using Multiple-Davydov D2 Ansatz

The linear absorption spectrum of J and H molecular aggregates is studied using the time-dependent Dirac-Frenkel variational principle (TDVP) with the multi-Davydov D2 (mD2) trial wavefunction (Ansatz). Both the electronic...

Lamb Shift and the Vacuum Rabi Splitting in a Strongly Dissipative Environment

We study the vacuum Rabi splitting of a qubit ultrastrongly coupled to a high-Q cavity mode and a radiation reservoir. Three methods are employed: a numerically exact variational approach with a multiple Davydov ansatz, the rotating-wave approximation (RWA), and the transformed RWA. Agreement between the variational results and the transformed RWA results is found in the regime of validity of the latter, where the RWA breaks down completely. We illustrate that the Lamb shift plays an essential role in modifying the vacuum Rabi splitting in the ultrastrong coupling regime, leading to off-resonant qubit-cavity coupling even though the cavity frequency equals the bare transition frequency of the qubit. Specifically, the emission spectrum exhibits one broad low-frequency peak and one narrow high-frequency peak in the presence of relatively weak cavity-qubit coupling. As the cavity-qubit coupling increases, the low-frequency peak narrows while the high-frequency peak broadens until they have similar widths.

Vibration-mediated energy transport in bacterial reaction center: Simulation study

Exciton energy relaxation in a bacterial Reaction Center (bRC) pigment-protein aggregate presumably involves emission of high energy vibrational quanta to cover wide energy gaps between excitons. Here, we assess this hypothesis utilizing vibronic two-particle theory in modeling of the excitation relaxation process in bRC. Specific high frequency molecular vibrational modes are included explicitly one at a time in order to check which high frequency vibrations are involved in the excitation relaxation process. The low frequency bath modes are treated perturbatively within Redfield relaxation theory. The analysis of the population relaxation rate data indicates energy flow pathways in bRC and suggests that specific vibrations may be responsible for the excitation relaxation process.

Singlet fission dynamics and optical spectra of pentacene and its derivatives

A multimode Brownian oscillator model is employed to investigate the absorption spectra of pentacene and its derivatives in solution and thin films. Excellent agreement has been obtained between simulated and measured absorption spectra. Furthermore, using parameters obtained from fitting the absorption spectra of these pentacene derivatives, the singlet fission dynamics and two-dimensional electronic spectra of an ab initio Hamiltonian are investigated by Dirac-Frenkel time-dependent variation with multiple Davydov trial states. It is found that the periodic wave packet motion induced in the displaced excited state, and the accompanying vibrational relaxation, can be visualized by two-dimensional electronic spectra at short times.

Schrödinger-Cat States in Landau-Zener-Stückelberg-Majorana Interferometry: A Multiple Davydov Ansatz Approach

Employing the time-dependent variational principle combined with the multiple Davydov D₂ Ansatz, we investigate Landau-Zener (LZ) transitions in a qubit coupled to a photon mode with various initial photon states at zero temperature. Thanks to the multiple Davydov trial states, exact photonic dynamics taking place in the course of the LZ transition is also studied efficiently. With the qubit driven by a linear external field and the photon mode initialized with Schrödinger-cat states, asymptotic behavior of the transition probability beyond the rotating-wave approximation is uncovered for a variety of initial states. Using a sinusoidal external driving field, we also explore the photon-assisted dynamics of Landau-Zener-Stückelberg-Majorana interferometry. Transition pathways involving multiple energy levels are unveiled by analyzing the photon dynamics.

Photon-assisted Landau-Zener transitions in a periodically driven Rabi dimer coupled to a dissipative mode

We investigate multiple photon-assisted Landau-Zener (LZ) transitions in a hybrid circuit quantum electrodynamics device in which each of two interacting transmission-line resonators is coupled to a qubit, and the qubits are driven by periodic driving fields and also coupled to a common phonon mode. The quantum state of the entire composite system is modeled using the multi-D₂Ansatz in combination with the time-dependent Dirac-Frenkel variational principle. Applying a sinusoidal driving field to one of the qubits, this device is an ideal platform to study the photon-assisted LZ transitions by comparing the dynamics of the two qubits. A series of interfering photon-assisted LZ transitions takes place if the photon frequency is much smaller than the driving amplitude. Once the two energy scales are comparable, independent LZ transitions arise and a transition pathway is revealed using an energy diagram. It is found that both adiabatic and nonadiabatic transitions are involved in the dynamics. Used to model environmental effects on the LZ transitions, the common phonon mode coupled to the qubits allows for more available states to facilitate the LZ transitions. An analytical formula is obtained to estimate the short time phonon population and produces results in reasonable agreement with numerical calculations. Equipped with the knowledge of the photon-assisted LZ transitions in the system, we can precisely manipulate the qubit state and successfully generate the qubit dynamics with a square-wave pattern by applying driving fields to both qubits, opening up new venues to manipulate the states of qubits and photons in quantum information devices and quantum computers.

A hierarchical effective mode approach to phonon-driven multilevel vibrational relaxation dynamics at surfaces

We discuss an efficient Hierarchical Effective Mode (HEM) representation of a high-dimensional harmonic oscillator bath, which describes phonon-driven vibrational relaxation of an adsorbate-surface system, namely, deuterium adsorbed on Si(100). Starting from the original Hamiltonian of the adsorbate-surface system, the HEM representation is constructed via iterative orthogonal transformations, which are efficiently implemented with Householder matrices. The detailed description of the HEM representation and its construction are given in the second quantization representation. The hierarchical nature of this representation allows access to the exact quantum dynamics of the adsorbate-surface system over finite time intervals, controllable via the truncation order of the hierarchy. To study the convergence properties of the effective mode representation, we solve the time-dependent Schrödinger equation of the truncated system-bath HEM Hamiltonian, with the help of the multilayer extension of the Multiconfigurational Time-Dependent Hartree (ML-MCTDH) method. The results of the HEM representation are compared with those obtained with a quantum-mechanical tier-model. The convergence of the HEM representation with respect to the truncation order of the hierarchy is discussed for different initial conditions of the adsorbate-surface system. The combination of the HEM representation with the ML-MCTDH method provides information on the time evolution of the system (adsorbate) and multiple effective modes of the bath (surface). This permits insight into mechanisms of vibration-phonon coupling of the adsorbate-surface system, as well as inter-mode couplings of the effective bath.

Quantum Dynamics Simulation of Doublet Excitation and Magnetic Field Effect in Neutral Radical Materials

The photon absorption and the relevant magnetic field effect of a doublet in neutral radical materials are investigated by combining the density functional theory and the variational quantum dynamics with Davydov ansatz. The doublet is a perfect model system to study the full-quantum dynamics of a two-level system coupling to a realistic molecular vibrational environment. In this work, we simulate the optical absorption spectroscopy of the neutral radical material, (4-N-carbazolyl-2,6-dichlorophenyl)bis(2,4,6-trichlorophenyl)-methyl, and find a good agreement with experiments for both highest occupied molecular orbital-singly occupied molecular orbital (SOMO) and SOMO-lowest unoccupied molecular orbital transitions. The nontrivial role of the intramolecular vibronic couplings is comprehensively discussed with separate spectroscopy and population dynamics, suggesting different contributions and the long time scale effect stemming from the vibrations, according to different symmetries. On the basis of the model, an applied magnetic field is taken into account to qualitatively investigate its magnetic properties in a dynamics manner, leading to a result which can be described by a sum of Lorentzian functions.

Modeling irreversible molecular internal conversion using the time-dependent variational approach with sD2 ansatz

A model of irreversible molecular internal conversion dynamics due to molecular thermal energy dissipation to the bath is presented.

Chebyshev hierarchical equations of motion for systems with arbitrary spectral densities and temperatures

The time evolution in open quantum systems, such as a molecular aggregate in contact with a thermal bath, still poses a complex and challenging problem. The influence of the thermal noise can be treated using a plethora of schemes, several of which decompose the corresponding correlation functions in terms of weighted sums of exponential functions. One such scheme is based on the hierarchical equations of motion (HEOM), which is built using only certain forms of bath correlation functions. In the case where the environment is described by a complex spectral density or is at a very low temperature, approaches utilizing the exponential decomposition become very inefficient. Here, we utilize an alternative decomposition scheme for the bath correlation function based on Chebyshev polynomials and Bessel functions to derive a HEOM approach up to an arbitrary order in the environmental coupling. These hierarchical equations are similar in structure to the popular exponential HEOM scheme, but are formulated using the derivatives of the Bessel functions. The proposed scheme is tested up to the fourth order in perturbation theory for a two-level system and compared to benchmark calculations for the case of zero-temperature quantum Ohmic and super-Ohmic noise. Furthermore, the benefits and shortcomings of the present Chebyshev-based hierarchical equations are discussed.

Davydov-Ansatz for Landau-Zener-Stueckelberg-Majorana transitions in an environment: Tuning the survival probability via number state excitation

We theoretically investigate transitions in a two-level system, which are induced by a sweep through an avoided crossing in the presence of coupling to a single, excited bosonic mode. This allows us to propose an initial number-state bosonic excitation as a new possible control parameter for the survival probability at long times. The expansion of number states in terms of coherent states centered around points on a circle in phase space makes a multi-Davydov-Ansatz the method of choice to perform the required numerical calculations. It is revealed that the starting time of the transition greatly affects the final transition probabilities. In addition, we found that the mixing angle, which is tuning between the diagonal and off-diagonal coupling, is decisive for the ability to control the transition via number state excitation. For a mixing angle of π/4, we found the maximal effect of number state excitation on the transition probability.

Dissipative dynamics in a tunable Rabi dimer with periodic harmonic driving

Recent progress on qubit manipulation allows application of periodic driving signals on qubits. In this study, a harmonic driving field is added to a Rabi dimer to engineer photon and qubit dynamics in a circuit quantum electrodynamics device. To model environmental effects, qubits in the Rabi dimer are coupled to a phonon bath with a sub-Ohmic spectral density. A nonperturbative treatment, the Dirac-Frenkel time-dependent variational principle together with the multiple Davydov D₂ ansatz, is employed to explore the dynamical behavior of the tunable Rabi dimer. In the absence of the phonon bath, the amplitude damping of the photon number oscillation is greatly suppressed by the driving field, and photons can be created, thanks to the resonance between the periodic driving field and the photon frequency. In the presence of the phonon bath, one can still change the photon numbers in two resonators and indirectly alter the photon imbalance in the Rabi dimer by directly varying the driving signal in one qubit. It is shown that qubit states can be manipulated directly by the harmonic driving. The environment is found to strengthen the interqubit asymmetry induced by the external driving, opening up a new venue to engineer the qubit states.

Dynamical scaling in the Ohmic spin-boson model studied by extended hierarchical equations of motion

Through a decomposition of the bath correlation function, the hierarchical equations of motion are extended to the Ohmic spin-boson model at zero temperature. For two typical cutoff functions of the bath spectral density, the rate kernel of spin dynamics is numerically extracted by a time-convolution equation of the average magnetic moment. A characteristic time is defined accordingly as the inverse of the zeroth-order moment of the rate kernel. For a given Kondo parameter in the incoherent regime, the time evolution of average magnetic moments gradually collapses onto a master curve after rescaling the time variable with the characteristic time. The rescaled spin dynamics is nearly independent of the cutoff frequency and the form of cutoff functions. For a given cutoff frequency, the characteristic time with the change of the Kondo parameter is fitted excellently as a function of the renormalized tunneling amplitude. Despite a significant difference in definition, our result is in good agreement with the characteristic time of the noninteracting blip approximation.

Finite-temperature time-dependent variation with multiple Davydov states

The Dirac-Frenkel time-dependent variational approach with Davydov Ansätze is a sophisticated, yet efficient technique to obtain an accurate solution to many-body Schrödinger equations for energy and charge transfer dynamics in molecular aggregates and light-harvesting complexes. We extend this variational approach to finite temperature dynamics of the spin-boson model by adopting a Monte Carlo importance sampling method. In order to demonstrate the applicability of this approach, we compare calculated real-time quantum dynamics of the spin-boson model with that from numerically exact iterative quasiadiabatic propagator path integral (QUAPI) technique. The comparison shows that our variational approach with the single Davydov Ansätze is in excellent agreement with the QUAPI method at high temperatures, while the two differ at low temperatures. Accuracy in dynamics calculations employing a multitude of Davydov trial states is found to improve substantially over the single Davydov Ansatz, especially at low temperatures. At a moderate computational cost, our variational approach with the multiple Davydov Ansatz is shown to provide accurate spin-boson dynamics over a wide range of temperatures and bath spectral densities.

Tracing feed-back driven exciton dynamics in molecular aggregates

Excitation, exciton transport, dephasing and energy relaxation, and finally detection processes shift molecular systems into a specific superposition of quantum states causing localization, local heating and finally excitonic polaronic effects.

Polaron dynamics in oligoacene stacks

The dynamical properties of polarons in organic molecular crystals are numerically studied in the framework of an one-dimensional Holstein-Peierls approach that includes lattice relaxation. Particularly, the present study is aimed at designing a tight-binding Hamiltonian that can address the charge transport mechanism in model oligoacene stacks. Our findings show that the definition of a particular oligoacene system depends strictly on the employed set of parameters. The usefulness of this methodology is highlighted by analyzing the polaron's saturation velocity and, consequently, its stability in the presence of a damping term and substantially high electric field strengths. Importantly, these results may be useful for the designing of novel materials to be employed in the field of molecular electronics.

Dynamics of the two-spin spin-boson model with a common bath

Dynamics of the two-spin spin-boson model in the presence of Ohmic and sub-Ohmic baths is investigated by employing a multitude of the Davydov D1 trial states, also known as the multi-D1 Ansatz. Its accuracy in dynamics simulations of the two-spin SBM is improved significantly over the single D1 Ansatz, especially in the weak to moderately strong coupling regime. Validity of the multi-D1 Ansatz for various coupling strengths is also systematically examined by making use of the deviation vector which quantifies how faithfully the trial state obeys the Schrödinger equation. The time evolution of population difference and entanglement has been studied for various initial conditions and coupling strengths. Careful comparisons are carried out between our approach and three other methods, i.e., the time-dependent numerical renormalization group (TD-NRG) approach, the Bloch-Redfield theory, and a method based on a variational master equation. For strong coupling, the multi-D1 trial state yields consistent results as the TD-NRG approach in the Ohmic regime while the two disagree in the sub-Ohmic regime, where the multi-D1 trial state is shown to be more accurate. For weak coupling, the multi-D1 trial state agrees with the two master-equation methods in the presence of both Ohmic and sub-Ohmic baths, but shows considerable differences with the TD-NRG approach in the presence of a sub-Ohmic bath, calling into question the validity of the TD-NRG results at long times in the literature.