Quantum algorithms for quantum dynamics: A performance study on the spin-boson model

Quantum algorithms for quantum dynamics simulations are traditionally based on implementing a Trotter-approximation of the time-evolution operator. This approach typically relies on deep circuits and is therefore hampered by the substantial limitations of available noisy and near-term quantum hardware. On the other hand, variational quantum algorithms have become an indispensable alternative, enabling small-scale simulations on present-day hardware. However, despite the recent development of variational quantum algorithms for quantum dynamics, a detailed assessment of their efficiency and scalability is yet to be presented. To fill this gap, we applied a variational quantum algorithm based on McLachlan's principle to simulate the dynamics of a spin-boson model subject to varying levels of realistic hardware noise as well as in different physical regimes, and discuss the algorithm's accuracy and scaling behavior as a function of system size. We observe a good performance of the variational approach used in combination with a general, physically motivated wavefunction ansatz, and compare it to the conventional first-order Trotter-evolution. Finally, based on this, we make scaling predictions for the simulation of a classically intractable system. We show that, despite providing a clear reduction of quantum gate cost, the variational method in its current implementation is unlikely to lead to a quantum advantage for the solution of time-dependent problems.


I. INTRODUCTION
The simulation of quantum systems is one of the most promising applications of quantum computing [1], aiming to overcome the limits of classical computers when it comes to storing and manipulating exponentially large quantum states. However, many of the conceived quantum algorithms, claiming to offer exponential speed-up over classical counterparts, are too resource-intensive for available hardware and will only become practicable once fault-tolerance is reached. In turn, since today's noisy near-term quantum technology is characterised by low qubit counts (< 1000), short decoherence times (∼ 100 µs) and two-qubit gate errors (∼ 10 −3 ) [2,3], error-correction schemes cannot yet be implemented [4].
This has sparked the development of hybrid quantumclassical algorithms, or VQAs [5], that split the workload between a quantum and a classical processor. Most prominently, the variational quantum eigensolver (VQE) has become the standard-tool for eigenvalue problems [6][7][8]. With efficient encodings of variational states, VQE requires only shallow circuits and has enabled small-scale simulations of up to a few atoms already on present-day hardware [9][10][11].
However, for VQAs to be meaningful for near-term applications in the simulation of quantum dynamics, it is necessary to carefully evaluate their performance, including their stability under noisy hardware conditions. Furthermore, their versatility with different systems has to be assessed. Particularly, they rely on choosing a variational ansatz that is both compact and flexible enough to accurately represent the studied system during the entire dynamics. Finding such a variational form is itself highly non-trivial as already addressed in the literature [15], which is why often, a so-called heuristic, or hardwareefficient ansatz, is chosen. Such an ansatz is agnostic to the problem at hand and its underlying symmetries, resulting in high numbers of variational parameters which could potentially jeopardize desired quantum advantage. Hence, in order to better characterize these VQAs, their application to non-trivial systems [30] is essential.
In this work we propose a detailed study of the performance of Li's VQA [12] for solving the dynamics of a spin-boson model. Moreover, based on our results, we make predictions on scalability and possible quantum advantage with a particular focus on the comparison with Trotter-evolution. The spin-boson model presents itself as an ideal testbed due to its rich dynamics and high relevance for various areas of research, resulting in a multitude of theoretical [31][32][33] and experimental studies [34][35][36]. The generic model of a two-level system coupled to a bath of harmonic oscillators is of great importance in the study of light-matter interaction and, particularly so, in the description of optical cavities and superconducting circuits [37]. On the other hand, it may also be seen as an idealized model for the study of the non-adiabatic dynamics of molecules, where, in this case, the fermionic two-level system describes two molecular potential energy surfaces [26,38]. Recent efforts in the context of digital quantum computing have explored both the spin-boson model's stationary as well as dynamical properties [39][40][41].
In this work, we start by constructing a physically motivated time-dependent variational form. We then focus on the numerical stability of the algorithm in different physical regimes, as well as the effects of introducing realistic experimental noise. In the last section, we finally present a careful study on the scaling of the computational resources as a function of the system size, comparing the variational approach and Trotter-evolution. In particular, we present predictions for system sizes far out of reach for classical simulation and conclude on the possibility to reach quantum advantage using near-term and fault-tolerant quantum algorithms for quantum dynamics.

II. THEORY
A. Quantum dynamics with product formulas As eluded to in the introduction, the most widely used method for time-evolution in the context of quantum computing remains the approximation of the unitary time evolution operator with a Trotter-Suzuki formula. At first order and with H = with an error that scales with O(N 2 h t 2 /d). It can be shown, however, that for Hamiltonians which can be mapped to a qubit-lattice and split into even and odd parts, as is the case for the spin-boson Hamiltonian introduced below, this scaling reduces to linear in the number of Hamiltonian terms [25,42], The drawback of this method is that it typically requires long circuits due to the error scaling quadratically with the simulation time.

B. Variational quantum algorithm for real time evolution
Alternatively, variational time-evolution algorithms for quantum dynamics aim to drastically reduce the circuit depth. A time-dependent variational ansatz |Φ(θ) , with θ = θ(t), seeks to approximate the true state |Ψ(t) , obtained as a solution to the time-dependent Schrödinger equation (TDSE) i d|Ψ dt = H |Ψ . The parameter's timedependence will be left implicit in the following and we set = 1.
On a quantum computer, a variational ansatz is prepared by acting upon a reference qubit-state |φ with a parameterized unitary operator, the quantum circuit, |Φ(θ) = U (θ) |φ , where θ = (θ 1 , θ 2 , . . .) ∈ R N θ is a set of real parameters. Although variational parameters can generally be complex, they are, in fact, required to be real in the setting of quantum computation since they will be encoded as angles of rotational quantum gates. As outlined in [12,43], such a time-dependent varational ansatz can be employed in a hybrid quantum-classical algorithm.
One of three variational principles (VPs) [44][45][46], the Dirac-Frenkel variational principle (DFVP) [47,48], the McLachlan variational principle (MVP) [49], and the time-dependent variational principle (TDVP) [50], may then be used to derive a set of equations of motion (EOMs) dictating the parameter evolution. In fact, as is intelligibly shown in [46], all three principles are equivalent under the condition that the variational manifold M is such that |δΦ and i |δΦ are both elements of the same tangent space. This is typically satisfied for a complex parameterization but not for purely real parameters [44], as is the case here. In fact, while parameters have to be made real artificially with the DFVP, both the MVP and the TDVP naturally maintain a real parameterization [43,44]. Due to known instabilities in the integration of the EOMs resulting from the TDVP, we will make use of MVP, where variation is with respect to |Θ = |Φ . Assuming the evolution of |Φ to be governed by the same TDSE as that of |Ψ , this means to minimize the distance between the projection H |Φ and the variational tangent vector d |Φ /dt. Eq. (3) results in the condition δΦ|i∂ t − H|Φ = 0. With all time-dependence residing in the parameters θ and accounting for a potential global phase mismatch between exact and approximate state, i.e. taking |Φ → e iα(t) |Φ , we obtain as EOMs [12,43,44] with the matrix elements and the vector components where the respective second term results from the inclusion of a global phase. Eq. (4) may then be solved by any numerical ODE-solver, e.g., a Runge-Kutta method. . . .
. . . Figure 1. Schematic representation of the qubit-mapped spin-boson model. The spin's state is captured by a single qubit, while under the direct mapping for bosonic modes, each energy level n k corresponds to a qubit.
We highlight that the MVP and the previous derivation is not immanent to quantum computation but may be used for any classical variational ansatz. What is distinct in the quantum setting is the preparation of the ansatz and the evaluation of individual terms by means of quantum circuits [12,51,52].
Recently, several other VQAs for quantum dynamics were proposed, relying either on propagating parameters by means of an EOM like Eq. (4) but differing in the way the ansatz is constructed [16,17,19], or by carrying out an optimization at each timestep [13-15, 18, 20]. In this last case, one can minimize for instance the distance between a variational state and the outcome of a small Trotter-step, avoiding the measurement-intensive construction of the matrix elements required in Eq. (4) as well as its inversion, which is a potential source of numerical instabilities. Concerning the optimization of variational quantum circuits, although it was shown in Ref. 53 that such optimization is in general NP-hard due to unresolvable local minima, approximate solutions suffice and can be found efficiently in practical simulations (cf. Solovay-Kitaev theorem [54]). Herein, we will make use of the original variational approach in Eq. (4), which solely relies on the integration of an EOM and does not involve any parameter optimization.

C. The spin-boson model
We consider a two-level system coupled to a bath of M bosons. The two-level system may represent an atom with two energy levels, a spin-1 2 particle, or any artificial system such as, for instance, a superconducting qubit. For brevity, we will refer to it simply as 'the spin'. Such a system is described by the spin-boson Hamiltonian [37,40], The bosonic operators a † k (a k ) create (annihilate) harmonic basis states with eigenfrequencies ω k , Pauli matrices σ i , i ∈ {x, y, z} act on the state of the spin with eigenfrequency and tunneling rate ∆. The coupling be-tween spin and bosons is via σ x with coupling constants g k .
Simulation on a quantum device requires to encode states in qubit registers and map operators to quantum gates, e.g., to strings of Pauli operators. Note that, in the following, the notation will be largely adapted from [40]. The excitation space of the k-th bosonic state will be truncated at a maximum occupation number n max k , leaving n max k + 1 possible occupations per mode k, including the ground state. Under the direct qubit-mapping [55], the occupation number vector (ONV) is then mapped to a qubit-register of size n max k Requiring to maintain correct spin-statistics, the corresponing mapping of bosonic creation and annihilation operators follows immediately as a †

D. The variational ansatz
The so-called polaron transformation (PT) enjoys popularity in the classical simulation of spin-boson models [33,56] and has successfully been used for VQE ground state calculations recently [40]. However, it proved to be insufficient for the use with variational time evolution and the Hamiltonian Eq. (7). Instead of the PT, here we employ a Variational Hamiltonian Ansatz (VHA) [57]. Inspired by the unitary time evolution operator, the timeparameter is simply replaced with a variational parameter that is distinct for each term in H, yielding Note that all Hamiltonian parameters are absorbed into variational parameters.
Translating U H into a sum of Pauli strings is now straight-forward. Employing the above operator mapping, we find with I the identity. Such identity terms contribute nothing but a global phase upon exponentiation and can thus be neglected in the variational ansatz.
Similarly, for the interaction term, one obtains Since this expression consists of mutually non-commuting terms, the summation over n k is split into even and odd parts, X e k := −i n k even √ n k + 1(σ x n k σ x n k +1 + σ y n k σ y n k +1 )/2, and analogously for the odd part X o k , such that Notably, we have [X e k , X o k ] = 0 while all terms within X e,o k commute. The resulting exponential is approximated with a Trotter-series of depth d, yielding an ansatz suitable for implementation in terms of quantum gates. Exponentials of Pauli terms may directly be written as rotational gates thereafter, e.g., . Although the final variational ansatz appears bulky, it can be compactly expressed as a series of one-and two-qubit gates and may be looked up in Appendix A.

E. Resource estimates and scaling
In this section, we discuss the scaling of the different computational resources for computing the dynamics of the spin-boson model with both the variational and the Trotter approach, Eq. (1) and Eq. (4), respectively. First off, the classical cost per timestep of the variational algorithm is determined by the number of variational parameters, which, for our ansatz (Eq. (8)), is given by The quantum cost is determined by qubit-and gatecounts, as well as the number of circuit evaluations. In the variational case, the total number of qubits is where an extra qubit was added to account for the possibility of evaluating gradients by means of an ancilla qubit [52]. Trotter evolution does not require any ancilla, hence requiring one qubit less, N q − 1.
The number of CNOT gates in the quantum circuit is ansatz-dependent and, for Eq. (8), can be estimated as Assuming the worst qubit-connectivity, i.e., a linear chain, we included a factor N q to account for swap gates that enter the circuit upon transpilation. This means, in the worst-case scenario, one needs to swap over the entire qubit register to execute a CNOT gate. This is true for both variational and Trotter simulation and, since our ansatz and the Trotter circuit differ only in the gate angles (which are variational parameters in the case of variational simulation), N cx is the same for both Trotter and variational simulation. Note, however, that the Trotter depth, d, differs in the two approaches; In the case of the Trotter algorithm, the circuit depth increases quadratically with the simulation time (cf. Eq. (2)), while for the variational approach, the depth (and the corresponding number of variational parameters) determines the size and nature of the sub-manifold governing the dynamics. The number of circuit evaluations per timesteps to evaluate the elements of M and V in Eq. (4) is determined by the number of Hamiltonian terms N h , the number of circuits necessary to evaluate all gradients ∂ i |Φ , which we denote N dθ , and the number of samples per circuit, In the variational ansatz Eq. (8), we have N θ = 2d(M n max + 1) variational parameters, a total of N dθ = d(5M n max + 2) gradient circuits, and N h = 7M n max + 2 Hamiltonian terms. The number of gradient terms differs from the number of parameters, as some parameters are repeated in the circuit.
Finally, we estimate the number of timesteps taken by the ODE solver to reach the final time T . Throughout this work, we will use an adaptive solver, however, adaptively choosing a step size is highly system dependent. Therefore, to simplify the estimate, we base it on the local error of a non-adaptive Runge-Kutta solver. For an order-p Runge-Kutta solver and a fixed timestep τ , the local error scales as ε local = O(τ p+1 ). For a desired final accuray of ε T , we thus estimate timesteps. With N t = T /τ , this means the timestep must satisfy τ p = O(ε T /T ). We emphasize that this is indeed a very rough estimate as scaling coefficients of the local error may heavily depend on the system under study. Especially so when considering an adaptive timestep. In that case, the number of function calls is what determines the number of circuit evaluations and thus also the cost of the algorithm, regardless of the number of accepted or rejected steps.

A. Noisy variational quantum simulation
In the following, we study the spin-boson model with the Hamiltonian in Eq. (7) for various Hamiltonian parameters and system sizes. In particular, we consider here the resonant case ω k ≡ ω, g k ≡ g and the regime of ultrastrong coupling (USC), where g/ω ∈ [0.1, 1]. Note that, from here on, we will take H/ω such that all Hamiltonian parameters are expressed in terms of bosonic eigenfrequencies. We begin with the simplest case of a single bosonic mode with an excitation number cutoff at n max = 1, resulting in three qubits under the direct mapping. The coupling strength is fixed at g/ω = 0.5 and we distinguish ( , ∆) ∈ {(0, 0), (−1, 0), (0, 1)}. Furthermore, we prepare the initial state in the non-interacting ground state |01 b |0 s = |010 and monitor its evolution through the orientation of the spin, P z = σ z + 1 /2. Note that we use the reverse qubit-ordering notation as conventional in Qiskit [58].
We aim to investigate how much variational simulations are affected by varying levels of noise. To this end, we differentiate four regimes; one with statistical (shot) noise only, one with full hardware noise mimicking IBM's ibmq_santiago device [2], which belongs to IBM's 5-qubit Falcon processors, as well as two intermediate regimes.
The two intermediate regimes are achieved by mimicking a device through Qiskit's noise model feature and the possibility to isolate and manipulate specific noise components. This allows for a detailed study of the influence of current hardware noise. Particularly, for the intermediate noise regimes, we decrease the average one-and two-qubit gate errors, e 1qg and e 2qg , respectively, as well as readout errors of the device, e read , while simultaneously increasing average relaxation and dephasing times, T 1 and T 2 , respectively, by a factor η, e 1qg = e dev 1qg /η , e 2qg = e dev 2qg /η , e read = e dev read /η , To simulate this setup, we employ Qiskit's shot-based Qasm-simulator with 8192 shots per circuit evaluation, and SciPy's adaptive Runge-Kutta solver of order 5(4) [59]. Fig. 2 shows the results of these simulations with η ∈ {1, 2, 10, ∞}, where η = 1 denotes full hardware noise and statistical noise, whereas η = ∞ means no hardware noise, i.e., only statistical noise. We employed complete readout error mitigation [58] via 2 Nq calibration circuits where N q is the number of qubits. All evolutions in Fig. 2 were obtained with a variational circuit of Trotter depth d = 1, containing N θ = 4 variational parameters. The top row of panels a)-c) displays the evolution of P z (t), while the respective bottom row gives the infidelities ∆Φ(t) = 1 − | Φ(t)|Ψ(t) |, where Φ(t) is the propagated (noisy) variational state at time t, while Ψ(t) is the corresponding exact solution obtained by exponentiation of the Hamiltonian matrix. It is evident that mere statistical noise (η = ∞) yields high accuracy throughout the entire simulation time, with a final infidelity of O(10 −4 ) to O(10 −3 ). We note that this is achieved with O(10 3 ) integration steps and a total of O(10 7 ) shots throughout one simulation. While this is highly model dependent, we achieve a final accuracy at a fixed number of samples several orders of magnitude better than those estimated in [18].
Moreover, although the accuracy decreases with the introduction of hardware noise, the variational algorithm (using the proposed variational ansatz in Eq. (8)) achieves a final infidelity of O(10 −2 ) to O(10 −1 ), even with full hardware noise (η = 1). Importantly, despite a deviation of the variational state from the true state trajectory over time, basic physical properties of the system's evolution, such as the oscillation frequency, are reproduced at least qualitatively.
Systematically reducing the noise (η = 2 and η = 10) as in Eq. (17) gradually increases accuracy. Remarkably, for η = 10, the simulation accuracy is comparable to that without hardware noise (η = ∞). This becomes even more clear in Fig. 2 d), where we plot the mean error ∆Φ for all η and the respective system from a)-c).

B. Scaling up -simulating larger spin-boson systems
Next, we enlarge the system to five qubits, which will enable us to make better scaling predictions for classically intractable systems in Section IV. Based on the direct qubit-mapping, five qubits may represent two different systems -a spin coupled to one bosonic mode (M = 1) with an excitation cutoff at n max = 3, or a spin coupled to two bosonic modes (M = 2) with an excitation cutoff at n max = 1 each. The time-evolution of these two systems within the same setup as before are displayed in Fig. 3 a) and b), respectively. Note that, in this and the following subsection, we consider statevector simulations only. Top rows a.1), b.1) and bottom rows a.2), b.2) represent system parameters ( , ∆) = (−1, 0), (0, 1), respectively.
We observe that Trotter depth d = 1 does not offer enough variational flexibility in all cases anymore and a depth of d = 2 is necessary to account for the correct dynamics. This could be anticipated from a simple dimensional analysis of the Hilbert space. However, with a total of N θ = 16 (a) and N θ = 12 (b) real variational parameters at d = 2, the dimensionality of the variational state remains well below the exponential size of the full 5-qubit wavefunction, which would require 31 complex or 62 real parameters for full parameterization (two of the 2×2 5 real parameters may be fixed with norm and global phase).

C. Comparison with Trotter-evolution
With circuit depth and two-qubit gate-count being the main limiting factors in noisy near-term quantum simulation due to short coherence times, it is worthwhile to compare the variational results from Figs. 2 and 3 to the more resource-intensive Trotter-evolution. Importantly, in Trotter-evolution, the depth increases with simulation time, while in variational simulation, the ansatz-depth remains constant throughout the simulation. Henceforth, the question is how the ansatz-depth required by the variational approach scales with system size compared with Trotter-simulation.
To address this question, we use the first-order formula of Eq. (1) from Section II A and aim at finding the minimal Trotter depth d to achieve a final accuracy of ε thresh . For this, we compute the infidelity ∆Φ(t) after each Trotter step. Beginning with a single circuit layer, d = 1, we append a layer to the circuit every time the infidelity increases above the threshold and repeat the step until ∆Φ(t) < ε thresh again.
The findings shown in Fig. 4 emphasize the resourceefficiency of the variational approach when used in combination with a well-chosen ansatz. Here we plot the final Trotter depth necessary to keep ∆Φ(t) < ε thresh throughout a fixed simulation time of T = 10 with ε thresh ∈ {10 −2 , 10 −3 , 10 −4 } and for several system sizes indicated by the number of qubits, N q ∈ [3, . . . , 11]. These are compared to the smallest Trotter depth of the variational ansatz that achieves ∆Φ(t) ≤ 10 −4 in the variational simulations. Note that, with a growing number of qubits, the number of distinct spin-boson systems that can be mapped to N q increases. For example, systems with M = 2, n max = 1 and M = 1, n max = 3 both result in 5 qubits; systems with M = 2, n max = 4 and M = 5, n max = 1 both result in 11 qubits. Data points in Fig. 4 represent simulation results averaged over all possible systems with the same number of qubits and n max k ≡ n max . Detailed numbers may be looked up in Appendix C, Tables II and III. In the variational simulations of up to 11 qubits, d ≤ 4 Trotter steps sufficed to maintain a target infidelity ∆Φ(t) ≤ 10 −4 . On the other hand, using Trotter-evolution, the circuit size grows significantly faster with system size.
To underline the different scaling behaviors, we linearly fit the depths in Fig. 4 according to Eq. (2). Note that d ∝ N h ∝ N q since the number of terms in the qubitmapped Hamiltonian is N h = 7M n max +2. The fit results are represented by the lines on different scales (linear and log scale in the top and bottom row, respectively, as well as different system size regimes, left and right) and the parameters may be looked up in Appendix C, Table IV. Although the depth scales linearly with system size both with Trotter-evolution as well as with the variational approach, it becomes visible from the fitting lines that the scaling coefficients vastly differ in both cases (cf . Table IV). In fact, despite the obvious savings of the variational method in terms of circuit depth, the extrapolated number of circuit layers for a 120-qubit system is merely around two orders of magnitude smaller than that estimated for Trotter-evolution with a target accuracy of ε thresh = 10 −4 . While this may seem like a large resource saving, it has to be put in relation with additional computational costs associated with the variational scheme, as will become more clear in the following section.    spin-boson systems with 10 degrees of freedom per mode (M = 12, n max = 9). For the simulation of such a system as described in the previous sections, we would need to control N q = 121 qubits (122 with an ancilla in variational simulation). In the following, we want to estimate the computational effort necessary to simulate such a system with both the variational and the Trotter-based approach, making use of the scaling laws presented in Section II E and the fits reported in Fig. 4. Throughout this section, we will consider a maximal target error of ε ≤ 10 −4 during the entire simulation.
From the results in Fig. 4, we estimate that the propagation of such a system using the first-order product formula (Eq. (1)) will require a circuit depth d ≈ 3400 in order to achieve the desired accuracy. This value is obtained by taking the average prediction for both Hamiltonian regimes in Fig. 4, leading to N cx ≈ 10 7 twoqubit gates. Importantly, this is the only computational cost associated with Trotter-based simulation in order to reach the fixed final time of T = 10, as no classical data processing is required.
On the other hand, in variational simulations the cost is split into three main components: the one associated to the circuit length (gate counts), the number of circuit evaluations to compute the different matrix elements, and the classical data processing to obtain the parameter update. The same extrapolation based on Fig. 4 predicts a variational form with a depth d ≈ 31 to simulate N q = 122 qubits. In this case the number of 2-qubit gates amounts to N cx ≈ 10 5 . With a final time of T = 10 and a desired accuracy of ε ≤ 10 −4 , the total number of timesteps N t = ε/ε local (cf. Eq. (16)), needed to integrate the EOM, amounts to N t ≈ 100. Note that in the numerical studies presented above, we found good agreement with this scaling even with an adaptive timestep when carefully choosing absolute and relative error tolerances for acceptance criteria. We report exact numbers of function calls for statevector as well as noisy simulations in Appendix C, Table V. The total number of circuit evaluations for the variational case becomes Now, assuming a two-qubit gate length of 100 ns, executing the Trotter-circuit takes approximately 1 second. A single evaluation of the variational circuit, on the other hand, would last roughly 0.01 seconds. We neglect the additional time required to measure and reset qubits after each circuit evaluation and the speed-up from possibly parallelizing circuit evaluations in the variational approach, since these two effect counter each other. Under these assumptions, a variational simulation will take approximately 0.01 s × 10 18 = 10 15 s ≈ 3 × 10 7 yr.
This example illustrates that, although the variational procedure allows for shallower circuits and hence opens avenues for performing simulations of small systems on Final depth, i.e., number of Trotter steps, required to achieve an accuracy ∆Φ below ε thresh , comparing Trotter and variational simulation for different system sizes (data points). Top and bottom plots show the same data with a linear and log scale, respectively. Legend entries denote different values of ε thresh and variational simulation, respectively. Variational simulation achieves a final infidelity below 10 −4 throughout all numerical examples and is therefore to be compared to the ε thresh = 10 −4 Trotter curve. The scaling of the depth is estimated with a linear fit (see main text), which is used to extrapolate to system sizes of N qubits = 120, the number of qubits necessary for a simulation comparable to state of the art classical spin-boson simulations.
near-term quantum computers, it is unlikely that it will lead to quantum advantage for the simulation of spinboson models. In fact, the number of circuit evaluations quickly becomes prohibitive in this case. This issue was also raised by Barison and coworkers [18] who proposed a variational algorithm which reduces the scaling of N circ from quadratic to linear in the number of parameters. Although this step goes in the right direction, it is by itself not enough to make the variational approach feasible for the applications described here. It is worth mentioning that the scaling could be further reduced by improving the sampling procedure. However, for an optimal number of shots, N shots = O(log(1/ε)), and a linear scaling of N circ with the number of parameters, the number of circuit evaluation would still amount to N tot circ ≈ 10 10 , taking roughly 10 7 s ≈ 0.3 yr. At the same time, the scaling of product formulas is sub-optimal and novel algorithms exhibiting reduced complexity have been proposed recently. Most notably, qubitization [29] achieves a gate complexity linear and additive in time, O(t + log(1/ )), which is provably optimal. A rough estimate based on the asymptotic bounds presented in Ref. 29 suggests that qubitization requires a two-qubit gate count of N cx = O(10 5 ) for the above example, taking O(0.01 s) to run, and additional O(10) ancilla qubits to implement the needed oracles. Despite its potential, the implementation of this algorithm poses important challenges to near-term quantum computing, which cannot be addressed in this work.

V. DISCUSSION AND CONCLUSIONS
In this work, we investigated the performance of a time-evolution VQA by simulating the quantum dynamics of a spin-boson Hamiltonian, a model which is widely used to describe the embedding of a two-level system in a thermal bath. Aside from assessing the VQA's numerical stability, the purpose of our investigation is to provide scaling estimates and predictions, particularly in comparison to conventional Trotter-evolution. In particular, we analyzed the performance of these time-evolution algorithms in the regimes of near-term and fault-tolerant quantum computing.
To this end, we studied the dynamics of several spinboson systems, varying in size as well as in the Hamiltonian parameter space (i.e., the Hamiltonian coefficients). Furthermore, we introduced hardware noise into the variational simulations by using the noise model of one of IBM's quantum computers, which provided a clear upper bound for the level of noise tolerated by the algorithm. Throughout all simulations, the physically motivated variational ansatz, which we constructed based on the system Hamiltonian, offered a great deal of flexibility and correctly captured various system's dynamics without the need of further tuning the variational quantum circuits. Moreover, it exhibits linear scaling of both the number of variational parameters and circuit depth.
Concerning the scaling of the two considered methods for time-evolution, namely the variational and the Trotter-based approach, we presented approximate scaling laws for the classical and the quantum computational resources required by both methods. We further performed a series of simulations for system sizes in the range N q ∈ {3, . . . , 11} to determine the required circuit depths for a fixed target error, and extrapolated these values to larger numbers of qubits using appropriate fitting models. Based on these extrapolations, we could estimate the computational cost of both methods for simulating system sizes, which are barely accessible with cutting-edge classical algorithms.
From this analysis, we can conclude that the variational approach in the current implementation is an efficient and reliable approach in the case of relatively small setups, especially in the context of current hardware limitations. However, the costs associated with the number of circuit evaluations will quickly become unaffordable, hampering its applicability to large setups. Note that this occurs despite the fact that the number of resources (variational parameters and gate count) only increases linearly with the system size. Although Trotter-evolution has a two-qubit gate count two orders of magnitude larger than the variational method, it does not suffer from the same prohibitive scaling of required measurements.
In conclusion, the variational algorithm might be a use-  Figure 5. a) Quantum circuit representation of the variational ansatz UH in Eq. (8), showing the spin and one bosonic qubit register. The blue box represents the coupling term in UH, while the red two-qubit gates represent bosonic selfinteraction terms. b) Decomposition of the coupling gate from a), with each cascade of three-qubit gates representing a term inX e,o k . d) Final decomposition of the three-qubit gate into one-and two-qubit gates, coupling the spin and the respective bosonic mode. c) Decomposition of the bosonic selfinteraction gate. Variational parameters enter through rotational gates Rz(θ). The Hadamard gate H and Y † = Rx(π/2) rotate qubits from σ z -into the σ x -and σ y -basis, respectively.

Appendix A: The variational quantum circuit
In this section, we detail the construction of the quantum circuit from the qubit-mapped ansatz U H in Section II D. A detailed representation for one bosonic mode is shown in Fig. 5, whereby the circuit for M bosonic modes is obtained by appending resepective bosonic qubit registers that are coupled to the qubit representing the spin via respective coupling gates (blue box in Fig. 5 Qubit a) and which have their own self-interaction gates (red two-qubit gates in Fig. 5 a). For all simulations in the main-text, we initialize each bosonic register in its non-interacting ground state |0 k = |0 . . . 01 , obtained with a bit-flip on the first k-mode qubit, |1 = X |0 = σ x |0 . Variational parameters enter through rotational gates as R z (θ) = exp(−iθσ z /2). Moreover, all parameters are n k -dependent, θ = θ n k , such that each bosonic selfinteraction and coupling gate (red two-and blue threequbit gates) is parameterized individually. Within each of these gates, however, for instance within the gate in c), all three R z -gates have the same parameter. While these parameters could be made independent as well if more flexibility is required of the ansatz, we found no advantage in doing so. 340 1  Table IV. Results of a linear fit f (x) = p1x + p0 to the data points in Fig. 4.