System-Environment Correlations in Qubit Initialization and Control

We employ numerically exact methods to study a qubit coupled to a quantum-mechanical bath, focusing on important correlation effects such as the Lamb shift and entanglement that are typically neglected in the modeling of contemporary quantum computers. We find a fundamental trade-off between speed and accuracy in qubit initialization which is important in the optimization of future reset protocols. Namely, we show analytically that at low temperatures, the deviation from the bare qubit ground state is proportional to the qubit decay rate, caused by the unavoidable entanglement between the qubit and the bath. We also find that the qubit decay is superexponential at short time scales, contrary to the result from standard Markovian approaches. However, the fidelities of quantum logic gates can be sufficiently accurately predicted by the Markovian methods. Our results can be used to develop quantum devices with engineered environments and to guide future experiments in probing the properties of the reservoirs.

Introduction-Precise control and preparation of pure quantum states are pivotal in quantum technological applications of practical interest [1,2]. For example, fast and high-fidelity initialization of a qubit to its ground state is required to realize a large-scale gate-based quantum computer since implementations of quantum error correction codes [3][4][5] call for pure ancillary qubits at each error correction cycle. However, satisfactory qubit reset still remains a technological challenge. In the most promising approaches, the qubit is steered towards the desired state by coherent driving [6][7][8][9] or by using a specifically tailored dissipative environment [10,11]. The latter has the benefit that its theoretical modeling does not rely on rotating frames which are often used in the case of time-dependent driving and can cause inaccuracies in the predicted figures of merit.
Fast initialization inherently calls for relatively strong environmental coupling, whereas coherent operations such as quantum logic are error free only in the limit of isolated quantum systems. This apparent conflict can be resolved with a dissipative low-temperature environment and temporal control over the coupling strength [10][11][12], providing very weak coupling for the coherent control and strong coupling for fast initialization to the ground state. Recently, superconducting-circuit realizations of the quantum-circuit refrigerator [13,14] and the tunable heat sink [15,16] have demonstrated that with the current technology, one can indeed control the coupling strength between a quantum system and the engineered bath over several orders of magnitude with a minimal effect on the system frequency. Such components can be conveniently integrated on the same chip with qubits, allowing scalable fabrication and low circuit complexity.
Estimates of the speed and the fidelity of the ini-tialization protocols based on qubit decay have been made in the weak-coupling, i.e., Born-Markov, approximation [8,10,11]. Stationary states then arise from a detailed balance condition of the Born-Markov rates and, therefore, they appear independent of the coupling strength. However, experiments with engineered quantum systems are entering a regime of high accuracy [17,18], where higher-order corrections need to be included. One such correction is the modification of equilibrium populations through a Lamb shift [19], which can be sizable in the case of a broadband environment [14,20,21]. System-environment entanglement is another higher-order effect detrimental to the performance of reset protocols. Any realistic analysis and optimization of the speed and the fidelity of the envisioned protocols thus calls for an exact analysis of dissipation that goes beyond the conventional weak-coupling formalism.
In this Letter, we examine open quantum dynamics of a qubit non-perturbatively and demonstrate that the effects of entanglement and Lamb shift lead to a decreasing ground-state occupation in the steady state with increasing bath coupling. This indicates a potential need to make a compromise between speed and accuracy in qubit initialization protocols. We observe a further departure from the behavior predicted by Born-Markov master equations in the transient dynamics of an initially decoupled qubit, displayed as a rapid initial decoherence into a mixture of pointer states [22,23]. For moderate and strong coupling, the initial transient dynamics has a Gaussian-shaped time dependence which is independent of the qubit frequency, also referred to as universal decoherence. In addition, we find qubit-reservoir entanglement to be the dominant source of initialization error at low temperatures, while strong-coupling effects are minor for quantum gates at experimentally relevant parameter values. Our findings can be used to improve qubit schemes involving reservoir engineering.
Qubit reset dynamics in tunable environments-As a generic situation for the qubit reset through a dissipative environment we consider, as shown in Fig. 1(a), a superconducting qubit with bare angular frequency ω q capacitively coupled to a tunable resistor at temperature T . The latter is realized using either a quantum-circuit refrigerator (QCR) or a tunable heat sink (THS). A typical power spectral density S(ω) of such an environment subject to the qubit is also depicted in Fig. 1(a) with a maximum around a cutoff frequency ω c and a zero frequency limit lim ω→0 S(ω) = κ/( βω q ) where β = 1/(k B T ) is the inverse temperature and the coupling parameter κ is identical to the zero-temperature qubit relaxation rate in the Born-Markov approximation.
These features can be conveniently modeled in terms of a spin-boson-like system with the dissipative environment bilinearly coupled to the qubit and consisting of an infinite set of harmonic oscillators [25], i.e., whereσ x = |g e| + |e g| andσ z = |g g| − |e e| are Pauli operators, |g and |e are the qubit ground and excited states, respectively, andâ k is the annihilation operator of oscillator mode k. Within this model, the power spectrum is obtained as S(ω) = J(ω)[n β (ω) + 1] with the Bose occupation of the bath modes n β (ω) = 1/[exp( βω) − 1] and the mode spectral density J(ω) = 2π k g 2 k δ(ω − ω k ), which becomes a smooth function in the limit of a large reservoir. According to Fig. 1(a), a second-order Drude model with J(ω) = (κ/ω q ) ω/[1 + (ω/ω c ) 2 ] 2 for the tunable resistor captures the relevant physics. Note that the above definitions imply that the ratio κ/ω q is independent of ω q for ω q ω c . Commonly, the qubit dynamics within this setting is described with the reduced density operatorρ the time evolution of which is assumed to follow from weakcoupling Redfield or Lindblad-type master equations (LE). However, the subtle quantum correlations between qubit and environmental degrees of freedom require a more sophisticated theoretical treatment to provide predictions which match the experimentally achievable accuracy. Suitable methods, originally developed in a condensed matter context [25], have been reviewed and found relevant in the context of quantum information [26]. Here, the Feynman-Vernon path integral formalism [27] underlying these methods will be replaced by equivalent stochastic Liouville-von Neumann equation (SLN) [28,29], unless analytic results exist.
The SLN provides an exact non-perturbative treatment of open quantum systems. It consists of the known  (dashed orange) interacting with a bath (dashed blue). We show the circuit diagram (left) of a superconducting qubit coupled to a tunable resistor which may be realized using a quantum circuit refrigerator or a tunable heat sink, together with the corresponding spectral densities (right). The bare qubit angular frequency is denoted by ωq and the Lamb-shifted quantity by ωLS. The dynamic properties of the bath are characterized by the power spectrum S(ω) of the bath fluctuations and the related mode spectral density J(ω) and the bosonic occupation n β (ω). The interaction strength between the qubit and the bath is parameterized through the zero-temperature Born-Markov decay rate κ. (b) Decay dynamics of the qubit excitation probability demonstrating initial universal decoherence, and the Lamb shift and entanglement at long times for the SLN and SLED methods (solid lines). The Redfield solution (dashed line) fails to capture these effects. Inset: The short-time behavior of the excited-state probability. We show the full universal decoherence (black dots) [24], and the short-time approximation f (t) ≈ 1 2 ω 2 c t 2 (dashed yellow) given in the text. The pink dashed line denotes the combined effect of the thermal and asymptotic results in Eq. (2). The vertical lines at ωqt = 0.1 and at ωct = 1 define the regions of validity for the indicated approximations. Here β ωq = 5, κ/ωq = 0.2, and ωc/ωq = 50. deterministic part describing the unitary evolution of the bare Hamiltonian and of a stochastic dissipative part where two complex-valued noise terms are matched with the free quantum fluctuations of the bath [30]. The physical reduced density operator is obtained by averaging over many realizations of the noise. In the parameter regimes relevant to quantum information technology, the SLN is superior compared to the path integral formulation with respect to flexibility and numerical performance. For a large cutoff frequency ω c ω q , the SLN equation can be reduced to involve only a single realvalued noise (stochastic Liouville equation with dissipation, SLED) [28,31] which is beneficial in certain ranges of parameter space. In particular, the SLN and SLED provide an excellent platform to accurately explore the qubit reset dynamics as they deliver reliable results also for strong qubit-reservoir interactions at low temperatures and for time-dependent parameters.
Relaxation dynamics-In Fig. 1(c), we monitor the decay of an excited qubit state as, in the presence of a dissipative environment, it relaxes at low temperatures towards thermal equilibrium. We observe that SLN and SLED results substantially differ from the predictions of the LE during the entire dynamics. Whereas the relaxation follows an exponential decay according to LE, the exact qubit dynamics exhibits various time domains of peculiar behavior. Note that we use in Fig. 1(b) a relatively strong environmental coupling, κ = 0.2 × ω q , as realized in recent protocols for engineered environments [13][14][15]. This is outside the strict applicability of the LE.
For short times, t 1/ω q , the bare qubit dynamics remains frozen and the qubit is only affected by the high-frequency reservoir modes [23]. As a consequence, an excited qubit state decays as ρ e (t) = e|ρ(t)|e = {1+exp[−f (t)κ/(πω q )]}/2, where both f (t) and κ/(πω q ) are system independent quantities, determined only by the reservoir [24]. Such decay, referred to as universal decoherence, can more concisely be described as dephasing in the pointer state basis ofσ x [22,23]. This behavior is depicted in the inset of Fig. 1(b), where we observe a good agreement between the analytical prediction and the numerically exact solution if ω q t 0.1. In particular, explicit expressions for f (t) for Drude-type dissipation can be found in limiting regimes, namely, f (t) ≈ 1 2 ω 2 c t 2 for ultra-short times ω c t < 1 and for times with min(t, β) 1/ω c . In this expression γ ≈ 0.577 denotes the Euler constant. These results indicate that the decay of a qubit excitation is superexponential at the time scale set by ω −1 c . Later, there is an algebraic decay at low T , especially for experimentally relevant cases of qubit initialization with ( β) −1 ω q ω c . The superexponential and asymptotic decays found above are shown in Fig. 1(b) and they agree well with the exact solution in their regimes of validity. Note that within the domain of universal decoherence, Eq. (2) predicts that ρ e → 1 2 for sufficiently large κ. This implies that a fully mixed state is reached before coherent qubit dynamics sets in. We emphasize that the phenomenon of universal decoherence is lost in the coarse-graining procedure underlying the derivation of the LE.
Steady state properties-Accurate predictions for the qubit steady state are crucial for the fidelity of initialization protocols. Figure 1(b) reveals that the bath coupling strength κ may affect the steady-state occupation of the qubit quite substantially. In fact, the ideal Boltzmann distribution at the bare qubit transition frequency ω q is obtained only in the limit κ → 0. This deviation can be attributed to both the downward Lamb shift of the qubit transition frequency which leads to excess thermal population and the entanglement of the qubit with the residual degrees of freedom.
In Fig. 2, we study in more detail the impact of qubitreservoir quantum correlations and the role of entanglement as the system approaches the steady state. We begin in Figs. 2(a) and 2(c) with the dependence of the steady-state probability on κ at an elevated temperature, βω q = 1. The excited-state occupation in the steady state ρ ∞ e = ρ e (t → ∞) increases with κ and the qubit approaches a fully mixed state already for κ/ω q > π/2. These numerical findings can further be substantiated by calculating the partition function based on a diagrammatic approach for κ/ω q 1 and large cutoff [24], and projecting it on the excited state population. This yields with the renormalized qubit frequency , with Γ(x) being the gamma function. This analytic result agrees well with the exact solution up to κ ≈ 0.2 × ω q . The inset in Fig. 2(a) shows that the decay of the Larmor precession of the pointer state with σ x = 1 also occurs at this renormalized angular frequency Ω < ω q , which is a signature of the reservoir-induced Lamb shift. More specifically, we find that for κ = 0.1 × ω q the renormalized angular frequency Ω ≈ 0.9 × ω q and has the relative error of two percent compared to the Lamb-shifted Larmor frequency obtained with an exponentially decaying sinusoidal fit for the exact result in the inset. This frequency renormalization reduces to the usual Lamb shift given in the literature only in the ultra-weak coupling limit [24,32]. Although one might expect that, at least for very small κ, the impact of the reservoir on the qubit can be captured solely by a renormalized frequency and the Markovian decay rate, this not the case. Namely, the zero temperature limit of Eq. (3) reads demonstrating the leading-order correction to the equilibrium state which can be shown to originate from the system-reservoir entanglement [24]. This is in full agreement with our numerical results at a low temperature, βω q = 10, cf. Figs. 2(b) and 2(d) [24]. Note that this result may exceed the corresponding exponentially small Boltzmann factor by orders of magnitude even for κ small enough to leave the transient dynamics virtually unaffected by dissipation.
Indeed, our findings provide a powerful tool to estimate experimentally achievable qubit initialization fidelities. For example, for an efficient implementation of the surface code one needs an error below 10 −3 [5] which, according to our SLN and SLED results, can only be reached at sufficiently low temperatures and only with coupling strengths κ/ω q 10 −3 . This suggests that there exists a trade-off between the speed for reservoirinduced reset and fidelity in qubit initialization. In fact, from the low-temperature relaxation rate for weak coupling κ, the above restriction determines a minimal reset time ω q t min ≈ 10 4 . Moreover, with a typical angular frequency of superconducting qubits ω q = 2π × 8 GHz and ω c < 2π×400 GHz, one may use a quantum circuit refrigerator to tune to κ = 10 −3 ×ω q = 1/(20 ns), and hence to reset the qubit to ρ e < 5 × 10 −4 in less than 200 ns. This would manifest a clear improvement to the current stateof-the-art experiments [8]. Note that, according to our findings, the fidelity in this example cannot be improved by simply lowering the reservoir temperature.
Impact on gate operations-Given the above subtle qubit-reservoir correlations, the question arises if they may influence also qubit operations other than the reset protocol, for example, high-fidelity gate operations. We study this in Fig. 3 for a π rotation aboutσ x , for various bath coupling strengths compared against the Rabi angular frequency g of the gate. For κ/g = 4 × 10 −5 , the conventional LE provides very accurate predictions to describe the fidelity, whereas for κ/g = 0.1, the fidelity is slightly affected by the qubit-reservoir quantum correlations. The gate error as a function of the environmental coupling strength is depicted in Fig. 3(b). Clearly, the error increases with κ as the impact of the reservoir leads to a mixing of the qubit state during the pulse. This effect is maximized at κ/g = 1 2 since we consider only the excited state as the initial state, and hence strong dissipation leads to fast increase of the desired ground state population. For the range of parameters used here, the discrepancies between LE and SLED become relevant for κ ∼ g, a value beyond the practical domain for the implementation of high-fidelity quantum gates.
Conclusions-We have shown that the impressive progress in fabricating and controlling architectures for quantum information processing calls for nonperturbative approaches beyond conventional weakcoupling master equations to reliably predict the impact of qubit-reservoir correlations for dissipative qubit initialization. In steady state, this includes quantification of the Lamb shift and bath entanglement effects, and the consequent trade-off between initialization speed and accuracy. Fortunately, our results indicate that this tradeoff does not seem to pose a fundamental problem on the route to scalable quantum computers if taken properly into account in the initialization protocol. Furthermore, we have demonstrated that the so-called universal decoherence describes early qubit dynamics up to times of the order of the inverse of the environment correlation time. This phenomenon is challenging to observe in typical Rabi-driven qubits, but may be visible in cases where the qubit Hamiltonian can be quickly controlled in the laboratory frame [33]. Finally, we have observed that the exact dynamics and that given by the Lindblad equation yield matching gate fidelities in parameter ranges of actual implementations. Our findings contain both caveats and suggestions for novel and intriguing uses of system-reservoir quantum correlations that lie beyond weak-coupling models and may direct the design of a new generation of state-of the-art experiments.
We give here a detailed derivation for the analytic results in the main text, including the steady state and the short-time evolution of the reduced density operator of a qubit interacting with a harmonic bath. We also show the equations for the SLN and SLED methods which were used to obtain the numerically exact results reported in the main text. Finally, we outline how the numerical implementation of the stochastic equations was made, and show benchmark data for our simulations.

I. ANALYTIC SHORT-TIME DYNAMICS AND STEADY STATE
In the main text, we have presented analytic results for the short-time and asymptotic behavior of the reduced density operatorρ S of the qubit, the time evolution of which is determined by the spin-boson Hamiltonian where {â k } are the annihilation operators of the bath oscillators, andσ z = |g g| − |e e| andσ x = |g e| + |e g| are Pauli matrices where |g and |e are the ground and excited states of the qubit. Here, we give a detailed derivation of these results.
A. Short-time decoherence In Fig. 1(b) of the main text, we observe that the short-time evolution of the excited state occupation ρ e = e|ρ S |e displays a rapid non-exponential drop. The drop occurs at a time scale that is shorter than the characteristic time scale of the system, set by ω −1 q . In this short-time limit, one can neglect the free evolution of the system and calculate the decoherence analytically using the Hamiltonian The calculation of this so-called universal decoherence was first carried out by Braun et al. [1] by employing the phase-space representation of the density operator for the whole derivation. Here, we give an alternative derivation using the operator formalism. We assume that initially the system and the bath are statistically independent. Accordingly, the initial density operator can be written asρ We further assume that the bath oscillators are initially in the thermal state which is determined by the inverse temperature β = 1/(k B T ) and can be expressed aŝ where Z k = Tr B exp(− ω kâ † kâ k ). In the interaction picture, the von Neumann equation can be written as withQ(t) = k g k (â † k e iω k t +â k e −iω k t ). This has the formal solution where T denotes time ordering.
Here, we study the operator partρ nm =Î ⊗ n|ρ|m ⊗Î of the joint density operator in the eigenbasis of operator σ x formed by the pointer states |n which obeyσ x |n = n|n . Here,Î is an identity operator for the bath and hencê ρ nm is a density operator for the bath. As a consequence, one obtainŝ We express where δt = t/N . Iteratively applying the Baker-Campbell-Hausdorff formula with [â k ,â † l ] = δ kl and [â k , where we have denoted D k (t) = sin(ω k t) + i[1 − cos(ω k t)].
The reduced density operator is obtained by tracing over the bath degrees of freedom as where on the last line we have used the thermal initial state for the bath, defined in Eq. (S4). Here, the trace can be simplified with the unitary rotationÛ = e −i k α k (t)â † kâ k where α k (t) are chosen such that the imaginary part of D k (t)â † k + D * k (t)â k is eliminated. We also note that the thermal-state density operator is diagonal in the eigenbasis ofâ † kâ k and, thus, unaffected by the rotation. We obtain where in the last equality, we have set α k (t) = −arg[D k (t)]. As a consequence of the transformation, we can write the reduced density operator as Thus, we need to calculate e ic k (â † k +â k ) in the thermal state, where This can be carried out in the phase-space representation, where the Wigner function for the thermal density operator ρ th k reads with the thermal averages q 2 k = q 2 k,0 coth( βω k /2) and p 2 k = q 2 k 2 /(4q 2 k,0 ) and q 2 k,0 being the ground-state width of mode k.
As a consequence, one arrives withq k = q k,0 (a † k + a k ) at Finally, we can write the elements of the reduced density matrix as where and we have recalled that where the latter equality holds for an ohmic bath with a second-order Drude cutoff at ω c . We emphasize that these are identical relations to those obtained previously by Braun et al. [1].
In the main text, we show data for the excited-state occupation of the operatorσ z given by where in the last equality we have assumed that initially ρ e (0) = 1. We have denoted where the pointer states |± are the eigenstates of theσ x operator obeyingσ x |± = ±|± .

Very-short-time behaviour
For very short times (ω c t < 1), the argument in the exponential reduces into where in the last equality we have also assumed zero temperature. We thus observe that, contrary to the Fermi's golden rule, the short-time dependence of the excited-state occupation is proportional to t 2 . This implies that for times obeying ω c t < 1, we expect faster than exponential decay of the initial excited-state occupation. This feature is demonstrated in Fig. 1(c) of the main text.

Thermal part
In general, the function f (t) that determines the short-time decoherence can be expressed as a sum of the zerotemperature and finite-temperature parts as One obtains an analytic solution for the thermal part by noticing that the Drude cutoff can be neglected at the presence of the thermal cutoff. Thus, one obtains for ohmic spectral density that the finite-temperature part of f (t) is given by This was obtained by first calculating the integral for the time derivative of f β (t) and then integrating the result in time.

Asymptotic behavior
We have calculated symbolically using Maple the asymptotic behavior of the zero-temperature part. The result can be written as where γ ≈ 0.577 . . . is the Euler constant. We emphasize that for the experimentally relevant case with ( β) −1 < ω q ω c , the short-time decoherence is determined accurately by the zero-temperature part of the function f (t) because the thermal time scale is longer than that of the system. The sum of Eqs. (S37) and (S36) is equal to the Eq. (2) of the main text.

B. Partition function approach for the steady state
The qubit occupation in the steady state can be calculated using the canonical partition function Z = Tr e −βĤ (S38) of the whole qubit-bath system. For clarity, we study the results in terms of the Kondo parameter K = κ/(πω q ). If the system depends on a parameter λ, one can write and, consequently, Due to the coupling, the expectation values of the qubit in the steady state differ from those obtained with the partition function Z 0 = exp( βω qσz ) of the bare qubit. In the following, we calculate these deviations in the regime of weak coupling with K sufficiently smaller than 1, and show that they are caused by the Lamb shift and the entanglement with the bath. The equilibrium properties of an open quantum system can be described by a reduced partition function [2] Z q with the property for any parameter λ which appears only in the system Hamiltonian (∂H I /∂λ = ∂H R /∂λ = 0). By setting λ = ω q in Eq. (S40), one obtains whereσ z = |g g| − |e e| is a Pauli operator of the bare qubit. The reduced partition function of the dressed qubit can be expressed in the weak coupling regime as [2] Z q = 2 cosh( βΩ/2), where ψ(z) is the digamma function, and with Γ(x) being the gamma function. The above result is equivalent to Eq. (4) of the main text and valid for all values of β. It has been derived using exponential cutoff, but we will show later that this assumption leads into minor deviations from the results given by the Drude cutoff used in the numerical simulations. Using the chain-rule of derivation in Eq. (S42), we obtain Eq.(3) of the main text: The expectation value comprises of two factors. We demonstrate below that the first factor describes the Lamb shift due to the renormalization of the qubit frequency by the bath. The other factor, ∂Ω/∂ω q , of the expectation value is a measure of entanglement between the qubit and the bath, and can be written as ∂Ω/∂ω q = ∂Ω/∂ω eff (∂ω eff /∂ω q ) where If one neglects the factor ∂Ω/∂ω q in the expression in Eq. (S47), one obtains where the derivative is with respect to the renormalized frequency Ω instead of the bare frequency ω q as in Eq. (S42). This describes a qubit with only a frequency renormalization ω q → Ω which is what one would expect for a system experiencing only the Lamb shift and no entanglement with the bath. In the literature [3], the Lamb shift is typically derived using the second-order perturbation theory with respect to the couplings g k . As a consequence, the Lamb-shifted transition frequency of the qubit can be written as where P stands for principal value. In the second equality, we have assumed zero temperature and ω c ω q , and used the spectral density J(ω) = πKω exp(−ω/ω c ) with an exponential cutoff. For Drude cutoff, the Euler constant γ = 0.577 . . . is replaced with 1 2 . On the other hand, one can make a linear expansion in K for the renormalized qubit frequency given in Eq. (S44). Note that at zero temperature and for ω c ω q the expression (S44) reduces to ω LS only under the much stricter constraint K ln(ω c /ω q ) 1. In Fig. S1, we compare the renormalized frequency Ω in Eq. (S44) with the Lamb-shifted qubit frequency ω LS in Eq. (S51). As expected, we observe that the perturbative result follows closely the renormalized frequency in Eq. (S44) for small values of the coupling constant K. At low temperatures with β > 5, the deviations appear for K 0.1. The perturbative nature of ω LS is emphasized by the fact that it decreases without a bound and becomes negative at K = [ln(ω c /ω q ) − γ] −1 ≈ 0.3, where the latter equality has been calculated for the parameters used in Fig. S1. On the other hand, the renormalized frequency Ω approaches zero asymptotically.
However, one cannot obtain Eq. (S50) as a limiting case to Eq. (S47). This would require ∂Ω/∂ω eff → 1 and ∂ω eff /∂ω q → 1. These limits can be reached only at zero temperature and zero K, i.e., when the bath can be neglected altogether. Therefore, the steady-state occupation of a qubit is never given by the Boltzmann distribution for the renormalized qubit frequency as the entanglement with the bath generates a notable correction for all β and K. Especially, when the temperature is zero, the excited-state occupation of the qubit is given solely by the entanglement with the bath, as we will show in the following sections. We have calculated the relative error numerically for five values of κ (dots). We also show an interpolated curve (solid) that goes through the data points. We have used the parameters βωq = 1 and ωc/ωq = 50.

Comparison with numerically obtained Larmor frequency
We have compared the renormalized frequency Ω in Eq. (S44) with the numerically obtained Larmor frequency ω 0 . In Fig. S2, we show the relative error The data shows that for κ 0.1 × ω q the relative error is below two percent. This further justifies the interpretation of Ω as the renormalized transition frequency of the qubit in the weak-coupling limit. We also note that the data for κ = 0.1 × ω q corresponds to the decaying Larmor oscillations shown in the inset of Fig. 2(a) of the main text.

Zero-temperature occupation
Here, we show that the coupling to the bath gives rise to a non-vanishing excited-state occupation of the qubit in the steady state, even in the zero-temperature limit βω eff → ∞. We assume a weak coupling (K 1) and a high cutoff (ω c ω q ). With these approximations, we obtain in Eq. (S47) that tanh( βΩ/2) ≈ 1 and Ω ≈ ω eff . Thus, We are interested in the steady state occupation ρ ∞ e = −( σ z − 1)/2 in the excited state |e of theσ z operator. We obtain up to the first order in K that It should be noted that for K > 0, this result deviates from the ρ ∞ e = 0 prediction given by the Boltzmann distribution for the bare qubit at T = 0. Thus, one can interpret that the non-zero occupation of the excited state cannot be treated as a Lamb shift and, thus, has to be generated by the entanglement between the qubit and the bath.
The difference between the constant term is caused by the fact that Eq. (S55) was derived using an exponential cutoff whereas here we used a Drude cutoff similar to our numerical simulations. Thus, we observe that the ∂Ω/∂ω q factor in Eq. (S47) arises due to entanglement. As a consequence, the entanglement is the dominating source of initialization error at low temperatures.

II. SLN AND SLED METHODS
If the coupling between the qubit and the bath cannot be treated as a weak perturbation, or the environmental correlation time is long, the typical Born-Markov approximation leading to Redfield and Lindblad master equations becomes inaccurate [2]. In such situations, one has to rely on more accurate methods, such as the formally exact Feynman-Vernon path integral formalism.

A. Stochastic Liouville-von Neumann equation
For bilinear coupling and the qubit-bath system starting from a factorized initial state given in Eq. (S3), one can show [4] that the path-integral representation for the reduced density operator dynamics can be cast into the form of the so-called stochastic Liouville-von Neumann (SLN) equation The SLN equation comprises of a deterministic coherent part given by the first term on the right side of the equation, followed by the stochastic dissipative part, the dynamical properties of which are set by the complex-valued random variables ξ and ν. The correlations between the qubit and the bath are encoded into the correlations of the noise terms which follow the equations where Θ(t) is the Heaviside step function and is the autocorrelation function of the bath withζ = k g k (â † k +â k ). We note that the correlations ξ(t)ξ * (t ) , ξ(t)ν * (t ) , and ν(t)ν * (t ) are not fixed by the bath correlation function, and can be thus chosen to optimize the efficiency of the numerical realization.

B. Stochastic Liouville equation with dissipation
In the case of Ohmic dissipation with a high Drude cutoff frequency ω c ω q , the path integral formalism can be reduced into the form of so-called stochastic Liouville equation with dissipation (SLED) [5,6] The above SLED has a stochastic part characterized by a single real-valued noise term ξ(t). The remaining terms form the deterministic part of the SLED. The autocorrelation function of the noise term is given by the real part of the bath correlation function as We emphasize that we have treated separately the classical white noise part in the bath correlation function, resulting in convergence of the numerical implementation of the method which is faster than is obtained without such separation.

III. IMPLEMENTATION OF THE SLN AND SLED METHODS
The SLN equation (S66) and the SLED in Eq. (S71) can be solved for each noise realization, i.e. sample, using the conventional methods for deterministic differential equations. The samples are generated in a discretized time grid [0, h, . . . , (N − 1)h] with the finite step size h, spanning from the initial value at t = 0 over an interval lasting for several relaxation times κ −1 (N − 1)h.

A. Generation of correlated-noise samples
The numerical implementation of the correlated-noise samples is the main difference to the corresponding Lindblad algorithm. Here, we describe a numerical scheme for the generation of such samples obeying the correlation functions (S67)-(S69). We note that these correlation functions arise in the SLN description of the exact path-integral formalism where the bath correlation function L(t) is defined as in Eq. (S70). In the case of the SLED, one only needs to consider Eq. (S67) where the real part of the bath correlation function is defined in Eq. (S72). We give here the method of generating the correlated complex-valued noise terms ξ and ν for the SLN equation, and then discuss in the end how the noise generation is simplified for the case of the SLED.
We first divide the ξ-noise into two parts as ξ(t) = ξ r (t) + ξ c (t), where ξ r (t) is assumed real and and ξ c (t)ξ c (t ) = ξ r (t)ξ c (t ) = ξ r (t)ν(t ) = 0. Moreover, we denote ξ c (t) = ξ R c (t) + iξ I c (t) and ν(t) = ν R (t) + iν I (t), where the terms ξ R,I c and ν R,I are assumed real. These noises can be generated by filtering the independent Gaussian noise samples x 1 (t), x 2 (t), and x 3 (t) with appropriate window functions W 1 (t) and W 2 (t) as where we have defined the Fourier transformation of a function f (t) as The window functions are determined by the correlation functions in Eq. (S73) and Eq. (S74), and can be expressed as where the second equality in the last equation has been written for the ohmic spectral density J(ω) defined in Eq. (S64). Each noise term can be generated numerically by the following protocol: 1. Produce an array of N independent Gaussian variables {x(t )} corresponding to the grid points t = h with = 0, . . . , N − 1.
In the case of SLED, one needs to do this procedure only once, for ξ r , whereas for the SLN equation all five real-valued random variables ξ r , ξ R,I c , and ν R,I are needed.

B. Details about the numerical implementation
The individual solutions for a given sample do not have a physical interpretation but, nevertheless, the density operator can be obtained by taking an average over the solutions obtained with different samples. We have solved the SLN/SLED equations by representing the density operator using a vector notation which allows writing the stochastic equations in the formρ S = Lρ S where L is the Liouvillian superoperator of the SLN/SLED equation including both the deterministic and the stochastic parts. For a given noise sample, the 'Liouvillian' equation is solved deterministically using the Magnus integrator method up to the first order in the time step h.
As both equations can be treated using deterministic methods for each noise sample, the solution for a given sample has the same complexity as the corresponding Lindblad equation. The difference with respect to the conventional Born-Markov master equations arises from the fact that, for a given set of parameters, one has to solve the dynamical equation for many noise samples. The convergence of the averaging procedure depends heavily on the temperature of the bath as the number of needed samples increases rapidly with decreasing temperature.
For the parameters listed in the main text and Table I, our simulation runs approximately 10 4 sample points per second on a modern CPU. Because the noise samples are independent, the SLN and SLED methods are easily parallelizable. We have exploited this property at low temperatures where a large number of samples is needed. The parallelization was implemented with supercomputers at CSC -the Finnish IT Center for Science. The benchmarks for the numerical solution are listed in Table I for the data relevant for Figures 1 and 2 of the main text.