System-Environment Correlations in Qubit Initialization and Control

TL;DR

This paper investigates quantum correlations in superconducting qubits, revealing a fundamental speed-accuracy trade-off in initialization and providing insights for optimizing reset protocols and predicting gate fidelities.

Contribution

It employs numerically exact methods to analyze system-reservoir correlations in transmon qubits, highlighting intrinsic constraints in qubit reset speed and accuracy.

Findings

Identifies a speed-accuracy trade-off in qubit initialization.

Shows standard models suffice for predicting gate fidelities.

Provides a framework for experimental probing of qubit environments.

Abstract

The impressive progress in fabricating and controlling superconducting devices for quantum information processing has reached a level where reliable theoretical predictions need to account for quantum correlations that are not captured by the conventional modeling of contemporary quantum computers. This applies particularly to the qubit initialization as the process which crucially limits typical operation times. Here we employ numerically exact methods to study realistic implementations of a transmon qubit embedded in electromagnetic environments focusing on the most important system-reservoir correlation effects such as the Lamb shift and entanglement. For the qubit initialization we find a fundamental trade-off between speed and accuracy which sets intrinsic constraints in the optimization of future reset protocols. Instead, the fidelities of quantum logic gates can be sufficiently accurately predicted by standard treatments. Our results can be used to accurately predict the performance of specific set-ups and also to guide future experiments in probing low-temperature properties of qubit reservoirs.

Tables

Table 1. Table 1: Benchmarks for the numerical solution of the SLN equation and the SLED. We give examples for the data in Figs.  1 and  2 of the main text with fixed parameters κ / ω q = 0.2 𝜅 subscript 𝜔 q 0.2 \kappa/\omega_{\rm q}=0.2 and ω c / ω q = 50 subscript 𝜔 c subscript 𝜔 q 50 \omega_{\rm c}/\omega_{\rm q}=50 . The data have been obtained with the time step h = 2 − 7 × ω q − 1 ℎ superscript 2 7 superscript subscript 𝜔 q 1 h=2^{-7}\times\omega_{\rm q}^{-1} . We denote the number of time steps, i.e. the sample length, with N 𝑁 N and the number of samples with N S subscript 𝑁 S N_{\rm S} . The standard deviation of σ 𝜎 \sigma in the steady state is calculated in the interval t ∈ [ 9 ​ κ T − 1 , 10 ​ κ T − 1 ] 𝑡 9 superscript subscript 𝜅 T 1 10 superscript subscript 𝜅 T 1 t\in[9\kappa_{\rm T}^{-1},10\kappa_{\rm T}^{-1}] for all data sets.

Figure	$\hslash \beta {\omega }_q$	$N$	$N_S$	$\sigma$	method
1	5	$2^{13}$	$2\times {10}^7$	$1.0\times {10}^{-2}$	SLN
1	5	$2^{13}$	$5\times {10}^4$	$5.3\times {10}^{-4}$	SLED
2a	1	$2^{12}$	$2\times {10}^4$	$3.0\times {10}^{-4}$	SLED
2c	10	$2^{13}$	$2\times {10}^5$	$7.5\times {10}^{-4}$	SLED

Equations

$$\hat{H}=\hslash \sum _{n=0}^{N-1}{\omega }_n|n⟩⟨n|+\hslash \sum _k{\Omega }_k{\hat{a}}_k^{†}{\hat{a}}_k+\hslash \hat{q}\sum _k{g}_k({\hat{a}}_k^{†}+{\hat{a}}_k),$$

$$f(t)=2\left(\gamma -\frac{1}{2}+\ln ({\omega }_{\mathrm{c}}t)+\ln \left\{\frac{\sinh [\pi t/(\hslash \beta )]}{\pi t/(\hslash \beta )}\right\}\right),$$

$$⟨{\hat{\sigma }}_{\mathrm{z}}{⟩}^{\infty }=\tanh (\hslash \beta \Omega /2)\frac{\partial \Omega }{\partial {\omega }_{\mathrm{q}}},$$

$$\Omega ={\omega }_{\mathrm{eff}}{\left\{1+2K\left[\text{Re}\psi \left(i\frac{\hslash \beta {\omega }_{\mathrm{eff}}}{2\pi }\right)-\ln \left(\frac{\hslash \beta {\omega }_{\mathrm{eff}}}{2\pi }\right)\right]\right\}}^{1/2}.$$

$${\rho }_{\mathrm{e}}^{\infty }\approx \frac{\kappa }{2\pi {\omega }_{\mathrm{q}}}[-1-\gamma +\ln ({\omega }_{\mathrm{c}}/{\omega }_{\mathrm{q}})],$$

$${\hat{H}}_{\mathrm{S}}=4E_{\mathrm{C}}\hat{n}-E_{\mathrm{J}}\cos \hat{\phi },$$

$${\alpha }_{\mathrm{r}}=\frac{{\omega }_2-{\omega }_1}{{\omega }_1-{\omega }_0}-1\approx -\mathrm{0.04.}$$

$$\hat{H}=\hslash \sum _{n=0}^{N-1}{\omega }_n|n⟩⟨n|+\hslash \sum _k{\Omega }_k{\hat{a}}_k^{†}{\hat{a}}_k+\hslash \hat{q}\sum _k{g}_k({\hat{a}}_k^{†}+{\hat{a}}_k),$$

$$\hat{H}=-\frac{\hslash {\omega }_{\mathrm{q}}}{2}{\hat{\sigma }}_{\mathrm{z}}+\hslash \sum _k{\Omega }_k{\hat{a}}_k^{†}{\hat{a}}_k+\hslash {\hat{\sigma }}_{\mathrm{x}}\sum _k{g}_k({\hat{a}}_k^{†}+{\hat{a}}_k),$$

$$\hat{H}\approx {\hat{H}}_{\mathrm{B}}+{\hat{H}}_{\mathrm{I}}=\hslash \sum _k{\Omega }_k{\hat{a}}_k^{†}{\hat{a}}_k+\hslash {\hat{\sigma }}_{\mathrm{x}}\sum _k{g}_k({\hat{a}}_k^{†}+{\hat{a}}_k).$$

$$\hat{\rho }(0)={\hat{\rho }}_{\mathrm{S}}(0)\otimes {\hat{\rho }}_{\mathrm{B}}(0).$$

$${\hat{\rho }}_{\mathrm{B}}(0)=\underset{k}{⨂}{\hat{\rho }}_k^{\mathrm{th}}=\underset{k}{⨂}\frac{1}{Z_k}{e}^{-\hslash \beta {\Omega }_k{\hat{a}}_k^{†}{\hat{a}}_k},$$

$$\frac{\text{d}\hat{\rho }(t)}{\text{d}t}=-\frac{i}{\hslash }[{\hat{H}}_{\mathrm{I}}(t),\hat{\rho }(t)],$$

$${\hat{H}}_{\mathrm{I}}(t)=\hslash {\hat{\sigma }}_{\mathrm{x}}\hat{Q}(t),$$

$$\hat{\rho }(t)=\mathcal{T}{e}^{-\frac{i}{\hslash }{\int }_0^t\text{d}{t}^{\prime }{H}_{\mathrm{I}}(t^{\prime })}\hat{\rho }(0)\mathcal{T}{e}^{\frac{i}{\hslash }{\int }_0^t\text{d}{t}^{\prime }{H}_{\mathrm{I}}(t^{\prime })},$$

$${\hat{\rho }}_{nm}=\mathcal{T}{e}^{-in{\int }_0^t\text{d}{t}^{\prime }\hat{Q}(t^{\prime })}{\rho }_{nm}^{\mathrm{S}}(0){\hat{\rho }}_{\mathrm{B}}(0)\mathcal{T}{e}^{im{\int }_0^t\text{d}{t}^{\prime }\hat{Q}(t^{\prime })}.$$

$$\mathcal{T}{e}^{im{\int }_0^t\text{d}{t}^{\prime }\hat{Q}(t^{\prime })}=\underset{N\to \infty }{\lim }\prod _{k=N}^0{e}^{im\hat{Q}(k\delta t)\delta t},$$

$$\mathcal{T}{e}^{im{\int }_0^t\text{d}{t}^{\prime }\hat{Q}(t^{\prime })}$$

$$=$$

$${\rho }_{nm}^{\mathrm{S}}(t)$$

$$=$$

$${\hat{U}}^{†}\left[D_k(t){\hat{a}}_k^{†}+D_k^{\ast }(t){\hat{a}}_k\right]\hat{U}=D_k(t)e^{i{\alpha }_k(t)}{\hat{a}}_k^{†}+D_k^{\ast }(t)e^{-i{\alpha }_k(t)}{\hat{a}}_k=\sqrt{2[1-\cos ({\Omega }_{k}t)]}({\hat{a}}_k^{†}+{\hat{a}}_k),$$

$${\rho }_{nm}^{\mathrm{S}}(t)=\prod _k{\text{Tr}}_k\left[e^{-i(n-m)\frac{g_k}{{\Omega }_k}\sqrt{2[1-\cos ({\Omega }_{k}t)]}({\hat{a}}_k^{†}+{\hat{a}}_k)}{\hat{\rho }}_k^{\mathrm{th}}\right]e^{i(n^2-m^2){\sum }_{\ell }\frac{g_{\ell }^2}{{\omega }_{\ell }^2}[{\omega }_{\ell }t-\sin ({\omega }_{\ell }t)]}{\rho }_{nm}^{\mathrm{S}}(0).$$

$$c_k=-(n-m)\frac{g_k}{{\Omega }_k}\sqrt{2[1-\cos ({\Omega }_{k}t)]}.$$

$$W_{\beta }(q_k,p_k)=\frac{1}{2\pi \sqrt{⟨{\hat{q}}_k^2⟩⟨{\hat{p}}_k^2⟩}}{\mathrm{e}}^{-q_k^2/(2⟨{\hat{q}}_k^2⟩)-p_k^2/(2⟨{\hat{p}}_k^2⟩)},$$

$$⟨e^{i{c}_k({\hat{a}}_k^{†}+{\hat{a}}_k)}⟩=e^{-\frac{c_k^2}{2}\coth [\hslash \beta {\Omega }_k/2]}.$$

$${\rho }_{nm}^{\mathrm{S}}(t)=e^{-(n-m{)}^{2}f(t)\kappa /(4\pi {\omega }_{\mathrm{q}})+i(n^2-m^2)\phi (t)\kappa /(4\pi {\omega }_{\mathrm{q}})}{\rho }_{nm}^{\mathrm{S}}(0),$$

$$f(t)$$

$$\phi (t)$$

$$J(\omega )=2\pi \sum _k{g}_k^2\delta (\omega -{\Omega }_k)=\frac{\kappa (\omega /{\omega }_{\mathrm{q}})}{{\left(1+\frac{{\omega }^2}{{\omega }_c^2}\right)}^2},$$

$${\rho }_{\mathrm{e}}^{\mathrm{S}}(t)$$

Full text

††preprint: APS/123-QED

I Introduction

Precise control and preparation of pure quantum states are pivotal in quantum technological applications of practical interest DiVincenzo (2000); Nielsen and Chuang (2000). 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 Shor (1995); Fowler et al. (2009); Fowler (2011) 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 Valenzuela et al. (2006); Grajcar et al. (2008); Geerlings et al. (2013); Jin et al. (2015) or by using a specifically tailored dissipative environment Jones et al. (2013); Tuorila et al. (2017). 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 Jones et al. (2013); Tuorila et al. (2017); Silveri et al. (2017), providing very weak coupling during the coherent control and strong coupling during fast initialization to the ground state. Recently, superconducting-circuit realizations of the quantum-circuit refrigerator Tan et al. (2017); Silveri et al. (2019) and the tunable heat sink Partanen et al. (2018a, b) have demonstrated that with the current technology, one can indeed control the coupling strength between the 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 initialization protocols based on qubit decay have been made in the weak-coupling, i.e., Born–Markov, approximation Jones et al. (2013); Tuorila et al. (2017); Geerlings et al. (2013). Stationary states then arise from a detailed balance condition of the Born–Markov rates and, therefore, appear independent of the coupling strength. However, experiments with engineered quantum systems are entering a regime of high accuracy Barends et al. (2013); Kelly et al. (2015), where higher-order corrections need to be included. One such correction is the modification of equilibrium populations through a Lamb shift Fragner et al. (2008), which can be sizable in the case of a broadband environment Silveri et al. (2019); Gramich et al. (2011, 2014). System–environment entanglement is another higher-order effect detrimental to the performance of reset protocols. Any realistic analysis and optimization of the speed and fidelity of the envisioned protocols thus calls for an exact analysis of dissipation that goes beyond the conventional weak-coupling formalism.

Here, we examine non-perturbatively the open quantum dynamics of both an ideal two-level quantum system and a superconducting transmon qubit with $N$ energy levels Koch et al. (2007). We focus on the figures of merit important to the quantum information community and leave, e.g., more detailed studies of non-Markovianity Breuer et al. (2016); Rivas et al. (2014); Wolf et al. (2008); Luo et al. (2012); Chruściński and Maniscalco (2014); Basilewitsch et al. (2017) for future work. We 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 Zurek (2003); Braun et al. (2001). For moderate and strong coupling, the initial transient dynamics has a Gaussian temporal shape 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, whereas strong-coupling effects are minor for quantum gates at experimentally relevant parameter values. During initialization using strong coupling to an engineered bath, effects of the intrinsic qubit dissipation are small and, thus, can be neglected here (see Appendix C). Our findings can be used to improve qubit schemes involving reservoir engineering.

This paper is organized as follows. In Sec. II, we introduce a prototype system for studies of strong bath-coupling effects, consisting of a superconducting transmon qubit bilinearly coupled to a thermal bath. We also describe the numerically exact method used in our simulations. In Sec. III, we study the decay dynamics of the qubit and give a detailed description of the shortcomings of the Born–Markov master equations in terms of universal decoherence. Section IV presents an accurate calculation of the steady state. We compare the numerically exact data against the Boltzmann distribution of the bare qubit, and interpret the discrepancies analytically in terms of Lamb shift and entanglement with the bath. In Sec. V, we study the gate error arising from the weak-coupling approximation. We summarize our results in Sec. VI.

II Tunable environment for qubit initialization

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 $\omega_{\rm q}$ capacitively coupled to a tunable resistor at temperature $T$. The latter is realized using either a quantum-circuit refrigerator or a tunable heat sink. A typical power spectral density $S(\omega)$ of such an environment is also depicted in Fig. 1(a) with a maximum around a cutoff frequency $\omega_{\rm c}$ and a zero-frequency limit $\lim_{\omega\to 0}S(\omega)=\kappa/(\hbar\beta\omega_{\rm q})$ where $\beta=1/(k_{\rm B}T)$ and the coupling parameter $\kappa$ is identical to the zero-temperature qubit relaxation rate in the Born–Markov approximation.

These features can be conveniently modeled (see Appendix A) with a dissipative environment bilinearly coupled to the $N$-level transmon qubit and consisting of an infinite set of harmonic oscillators, i.e.,

[TABLE]

where $\hat{q}=\sum_{k,m}\langle k|\hat{n}|m\rangle|k\rangle\langle m|$, $\hat{n}$ is the Cooper-pair number operator of the transmon, $\omega_{n}$ and $|n\rangle$ are the eigenfrequencies and eigenstates of the transmon, respectively, and $\hat{a}_{k}$ is the annihilation operator of oscillator mode $k$. The transmon comprises a weakly nonlinear resonator with $\omega_{\rm q}=\omega_{1}-\omega_{0}$ and relative anharmonicity $\alpha_{\rm r}=(\omega_{2}-\omega_{1})/\omega_{\rm q}-1$. We restrict our discussion to $N$ lowest energy eigenstates. In the limit $N=2$ (ideal qubit), Eq. (1) reduces to the well-known spin-boson model Leggett et al. (1987); Weiss (1999) (see Appendix B) which we use in our analytic calculations and in some numerically exact simulations. We find that $N=5$ (transmon qubit) is enough for accurate studies of the low-energy dynamics at low temperatures.

Within this model, the power spectrum is obtained as $S(\omega)=J(\omega)[n_{\beta}(\omega)+1]$ with the Bose occupation of the bath modes $n_{\beta}(\omega)=1/[\exp(\hbar\beta\omega)-1]$ and the mode spectral density $J(\omega)=2\pi\sum_{k}g_{k}^{2}\delta(\omega-\omega_{k})$, which becomes a smooth function in the limit of a large reservoir. According to Fig. 1(a), a Drude model with $J(\omega)=(\kappa/\omega_{\rm q})\,\omega/[1+(\omega/\omega_{\textrm{c}})^{2}]^{2}$ for the tunable resistor captures the relevant physics. Accordingly, the power spectrum gives rise to a Markovian behavior (independent of frequency) only at high temperatures, while at low $T$ it displays a strong frequency dependence inducing non-Markovian dynamics. Note that the above definitions imply that the ratio $\kappa/\omega_{\rm q}$ is independent of $\omega_{\rm q}$ for $\omega_{\rm q}\ll\omega_{\rm c}$.

Commonly, the quantum dynamics within this setting is described with the reduced density operator $\hat{\rho}$ the time evolution of which is assumed to follow from weak-coupling Redfield- or Lindblad-type master equations (LEs). However, the subtle quantum correlations between a qubit and environment require a more sophisticated theoretical treatment that provides predictions which match the experimentally achievable accuracy. Suitable methods, originally developed in a condensed matter context Weiss (1999), have found use in quantum information Makhlin et al. (2001). Here, the Feynman–Vernon path-integral formalism Feynman and Vernon (1963), underlying these methods, is replaced by equivalent stochastic Liouville–von Neumann equation (SLN) Stockburger (1999a); Stockburger and Grabert (2002), unless analytic results exist (see Appendix D).

The SLN provides an exact non-perturbative treatment of open quantum systems. It augments the Liouville equation with two noise terms which are matched to the free quantum fluctuations of the bath 111The symmetric and antisymmetric parts of the two-sided spectral function of the bath are related by the fluctuation-dissipation theorem.. The physical reduced density operator is obtained by averaging over many realizations of the noise. For a high cutoff frequency $\omega_{\rm c}\gg\omega_{\textrm{q}}$, the SLN equation can be reduced to involve only a single real-valued noise [stochastic Liouville equation with dissipation (SLED)] Stockburger and Mak (1998); Stockburger (1999a).

III Decay dynamics

In Fig. 1(b), we monitor the decay of the first excited transmon state as 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 dynamics exhibits various time domains of peculiar behavior. Note that we use in Fig. 1(b) a relatively strong environmental coupling, $\kappa=0.2\times\omega_{\rm q}$, as realized in recent protocols for engineered environments Tan et al. (2017); Silveri et al. (2019); Partanen et al. (2018a). This is outside the strict applicability of the LE.

For early times, $t\ll 1/\omega_{\rm q}$, the ideal-qubit dynamics remains frozen and the qubit is only affected by the high-frequency reservoir modes Braun et al. (2001). As a consequence, the initialization error of the qubit decays as 1-\rho_{\rm g}(t)$$=\{1+\exp[-f(t)\kappa/(\pi\omega_{\rm q})]\}/2, where $\rho_{\rm g}=\langle 0|\hat{\rho}|0\rangle$, and both $f(t)$ and $\kappa/(\pi\omega_{\rm q})$ are system independent quantities, determined only by the reservoir (see Appendix B). Such decay, referred to as universal decoherence, can more concisely be described as dephasing in the pointer state basis of \hat{\sigma}_{\rm x}$$=|0\rangle\langle 1|+|1\rangle\langle 0| Zurek (2003); Braun et al. (2001). 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 $\omega_{\rm q}t\lesssim 0.1$. In particular, explicit expressions for $f(t)$ can be found for an ideal qubit in limiting regimes, namely, $f(t)\approx\frac{1}{2}\,\omega_{\rm c}^{2}t^{2}$ for ultrashort times $\omega_{\rm c}t<1$ and

[TABLE]

for times with $\min(t,\hbar\beta)\gg 1/\omega_{\rm c}$. Here, $\gamma\approx 0.577$ denotes the Euler constant. These results indicate that the decay of a qubit excitation is superexponential at the time scale set by $\omega_{\rm c}^{-1}$. Later, there is an algebraic decay at low $T$, especially for experimentally relevant cases of qubit initialization with $(\hbar\beta)^{-1}\ll\omega_{\rm q}\ll\omega_{\rm 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. The difference between the exact two- and multi-level dynamics during the early evolution indicates a leakage to transmon states $|n\rangle$ with $n>1$ at low temperature with $\hbar\beta\omega_{\rm q}=5$. This further validates the conclusion that the exact dynamics cannot be reconciled with a simple detailed-balance rate structure of the LE. We emphasize that the phenomenon of universal decoherence is lost in the coarse-graining procedure underlying the derivation of the LE.

IV 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 $\kappa$ may affect the steady-state occupation of the qubit substantially. In fact, the ideal Boltzmann distribution at the bare qubit transition frequency $\omega_{\rm q}$ is obtained only in the limit $\kappa\rightarrow 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 bath degrees of freedom.

In Fig. 2, we study in more detail the impact of qubit–reservoir quantum correlations and the role of entanglement as the system approaches the steady state. We show in Figs. 2(a) and 2(c) the dependence of the steady-state probability on $\kappa$ at an elevated temperature, $\hbar\beta\omega_{\rm q}=1$. The initialization error in the steady state $1-\rho_{\rm g}^{\infty}$ with $\rho_{\rm g}^{\infty}=\rho_{\rm g}(t\rightarrow\infty)$ increases with $\kappa$ and the transmon qubit approaches a fully mixed state already for $\kappa/\omega_{\rm q}>\pi/2$. These numerical findings can further be substantiated in the case of an ideal qubit by calculating the partition function based on a diagrammatic approach for $\kappa/\omega_{\rm q}\ll 1$ and large cutoff, and projecting it on the excited state population (see Appendix B). This yields  $\rho_{\rm e}^{\infty}=(1-\langle\hat{\sigma}_{\rm z}\rangle^{\infty})/2$, where $\hat{\sigma}_{\rm z}=|0\rangle\langle 0|-|1\rangle\langle 1|$ and Weiss (1999)

[TABLE]

with the renormalized qubit frequency

[TABLE]

Here, $\psi(x)$ is the digamma function, $\omega_{\rm eff}=G(\omega_{\rm q}/\omega_{\rm c})^{K/(1-K)}\omega_{\rm q}$, $K=\kappa/(\pi\omega_{\rm q})$, $G=[\Gamma(1-2K)\cos(\pi K)]^{1/[2(1-K)]}$, and $\Gamma(x)$ is the gamma function. This result agrees well with the exact solution up to $\kappa\approx 0.2\times\omega_{\rm q}$. The inset in Fig. 2(a) shows that the decay of the Larmor precession of the pointer state of an ideal qubit with $\langle\hat{\sigma}_{x}\rangle=1$ also occurs at this angular frequency $\Omega<\omega_{\rm q}$, which is a signature of the reservoir-induced Lamb shift. More specifically, we find that for $\kappa=0.1\times\omega_{\rm q}$ the renormalized angular frequency $\Omega\approx 0.9\times\omega_{\rm q}$ agrees within 2% compared to the Lamb-shifted Larmor frequency extracted using a fit for the exact result in the inset. This frequency renormalization reduces to the usual Lamb shift given in the literature only in the ultraweak-coupling limit Carmichael (1999) (see Appendix B). We also find that the inclusion of transmon states $|n\rangle$ with $n>1$ renders our system nearly harmonic and, consequently, leads to a suppression of the Lamb shift at high cutoff frequencies (data not shown, see Appendix B).

Although one might expect that, at least for very small $\kappa$, 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, taking the zero-temperature limit of Eq. (3), we obtain

[TABLE]

demonstrating the leading-order correction to the equilibrium state of an ideal qubit originating from the system–reservoir entanglement (see Appendix B). This is in full agreement with our numerical results at a low temperature, $\hbar\beta\omega_{\rm q}=10$ [cf. Figs. 2(b) and 2(d)]. Note that this result may exceed the corresponding exponentially small Boltzmann factor by orders of magnitude even for $\kappa$ small enough to leave the transient dynamics virtually unaffected by dissipation. We also show in Fig. 2(d) that the transmon states $|n\rangle$ with $n>1$ do not contribute significantly to the initialization error at low temperatures.

Our findings provide a powerful tool to estimate experimentally achievable qubit initialization fidelities. For example, efficient implementation of the surface code requires an error below $10^{-3}$  Fowler (2011) which, according to our results, can only be reached at sufficiently low temperatures and only with coupling strengths $\kappa/\omega_{\rm q}\lesssim 10^{-3}$. This suggests that there exists a trade-off between the speed and fidelity in reservoir-induced qubit initialization. In fact, from the low-temperature relaxation rate for weak coupling $\kappa$, the above restriction determines a minimal reset time $\omega_{\rm q}t_{\rm min}\approx 10^{4}$. Moreover, with a typical angular frequency of superconducting qubits $\omega_{\textrm{q}}=2\pi\times 8$ GHz and $\omega_{\textrm{c}}<2\pi\times 400$ GHz, one may use a quantum-circuit refrigerator to tune to $\kappa=10^{-3}\times\omega_{\textrm{q}}=1/(20\textrm{ ns})$, and hence to reset the qubit to 1-\rho_{\rm g}$$<5\times 10^{-4} in less than 200 ns. This would manifest a clear improvement to the current state-of-the-art experiments Geerlings et al. (2013). The fidelity in this example cannot be improved by simply lowering the reservoir temperature.

V Gate error

Given the above subtle qubit–reservoir correlations, the question arises if they may influence also other qubit protocols, such as high-fidelity gate operations. We study this for an ideal qubit in Fig. 3 for a $\pi$ rotation about $\hat{\sigma}_{\textrm{x}}$, for various bath coupling strengths compared against the Rabi angular frequency $g$ of the gate. For $\kappa/g=4\times 10^{-5}$, the LE provides very accurate predictions for the gate error, whereas for $\kappa/g=0.1$, the error is slightly affected by the qubit–reservoir quantum correlations. The gate error as a function of the bath coupling strength is depicted in Fig. 3(b). Clearly, the error increases with $\kappa$ as the impact of the reservoir leads to a mixing of the qubit state during the pulse. This effect is maximized at $\kappa/g=\frac{1}{2}$ since we consider only the excited state as the initial state, and hence strong dissipation leads to a fast increase of the desired ground-state population. For our range of parameters, the discrepancies between LE and SLED become relevant for $\kappa\sim g$, a value beyond the practical domain for the implementation of high-fidelity quantum gates.

VI Conclusions

We have shown that the impressive progress in fabricating and controlling devices for quantum information processing calls for non-perturbative approaches beyond conventional weak-coupling 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 trade-off 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 universal decoherence describes early qubit dynamics up to times of the order 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 Nakamura et al. (1999). 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 may direct the design of a new generation of state-of-the-art experiments.

Acknowledgements.

We thank Sahar Alipour, Rebecca Schmidt, and Matti Silveri for useful discussions. This work was financially supported in part by European Research Council under Grant No. 681311 (QUESS) and by Academy of Finland through its QTF Centre of Excellence program Grants No. 312298 and No. 312300. J.A. and J.T.S. thank the German Science Foundation (DFG) through AN336/12-1 and the IQST for financial support. The authors wish to acknowledge CSC – IT Center for Science, Finland, for generous computational resources.

Appendix A Transmon

In the main text, we have studied the system–bath correlations of a superconducting transmon qubit Koch et al. (2007). The transmon can be modeled with the Cooper-pair-box Hamiltonian

[TABLE]

where $E_{\rm C}$ and $E_{\rm J}$ are the charging energy of a Cooper pair and the Josephson energy, respectively, and $\hat{\varphi}$ and $\hat{n}=-i\partial/\partial\varphi$ are the superconducting phase and Cooper-pair number operators of the superconducting island. Contrary to a Cooper-pair box, the transmon is operated in the regime $E_{\rm J}\gg E_{\rm C}$ which leads to a suppression of charge noise. The charge noise depends exponentially on $-\sqrt{E_{\rm J}/E_{\rm C}}$. As a consequence, the Hamiltonian in Eq. (6) reduces to that of a harmonic oscillator with a weak anharmonicity proportional to $\hat{\varphi}^{4}$. In our exact simulations, we have used the numerically solved angular eigenfrequencies $\omega_{n}$ and the corresponding eigenvectors $|n\rangle$, with non-negative integer $n$, of the Hamiltonian in Eq. (6). We have used $E_{\rm J}/E_{\rm C}=100$, resulting in the relative anharmonicity

[TABLE]

A.1 Coupling with a harmonic bath

We assume that the transmon is bilinearly coupled through the Cooper-pair-number operator $\hat{n}$ to a dissipative environment consisting of an infinite set of harmonic oscillators. We represent the total Hamiltonian of the transmon–bath system in the eigenbasis of the transmon as

[TABLE]

where $\hat{q}=\sum_{k,m}\langle k|\hat{n}|m\rangle|k\rangle\langle m|$, and $\hat{a}_{k}$ are the annihilation operators of the bath oscillators.

Appendix B Analytic early dynamics and steady state

In the main text, we have presented analytic results for the early and asymptotic behavior of the reduced density operator $\hat{\rho}_{\rm S}$ in the case of $N=2$, i.e., an ideal qubit, the time evolution of which is determined by the spin-boson Hamiltonian [$N=2$ in Eq. (8)]

[TABLE]

where $\omega_{\rm q}=\omega_{1}-\omega_{0}$, and $\hat{\sigma}_{\rm z}=|\textrm{g}\rangle\langle\textrm{g}|-|\textrm{e}\rangle\langle\textrm{e}|$ and $\hat{\sigma}_{\rm x}=|\textrm{g}\rangle\langle\textrm{e}|+|\textrm{e}\rangle\langle\textrm{g}|$ are Pauli matrices where $|\textrm{g}\rangle=|0\rangle$ and $|\textrm{e}\rangle=|1\rangle$ are the ground and excited states of the qubit. Here, we give a detailed derivation of these results.

B.1 Early decoherence

In Fig. 1(b), we observe that the early evolution of the initialization error determined by the excited-state occupation $\rho_{\rm e}=\langle\textrm{e}|\hat{\rho}_{\rm S}|\textrm{e}\rangle$ 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 $\omega_{\rm q}^{-1}$. In this early-time limit, one can neglect the free evolution of the system and calculate the decoherence analytically using the Hamiltonian

[TABLE]

The calculation of this so-called universal decoherence was first carried out by Braun et al. Braun et al. (2001) 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

[TABLE]

We further assume that the bath oscillators are initially in the thermal state which is determined by the inverse temperature $\beta=1/(k_{\rm B}T)$ and can be expressed as

[TABLE]

where $Z_{k}=\textrm{Tr}_{\rm B}\exp(-\hbar\Omega_{k}\hat{a}^{{\dagger}}_{k}\hat{a}_{k})$. In the interaction picture, the von Neumann equation can be written as

[TABLE]

where

[TABLE]

with $\hat{Q}(t)=\sum_{k}g_{k}(\hat{a}_{k}^{{\dagger}}e^{i\Omega_{k}t}+\hat{a}_{k}e^{-i\Omega_{k}t})$. This has the formal solution

[TABLE]

where $\mathcal{T}$ denotes time ordering.

Here, we study the operator part $\hat{\rho}_{nm}=\hat{I}\otimes\langle n|\hat{\rho}|m\rangle\otimes\hat{I}$ of the joint density operator in the eigenbasis of operator $\hat{\sigma}_{\rm x}$ formed by the pointer states $|n\rangle$ which obey $\hat{\sigma}_{\rm x}|n\rangle=n|n\rangle$. Here, $\hat{I}$ is an identity operator for the bath and hence $\hat{\rho}_{nm}$ is a density operator for the bath. As a consequence, one obtains

[TABLE]

We express

[TABLE]

where $\delta t=t/N$. Iteratively applying the Baker–Campbell–Hausdorff formula with $[\hat{a}_{k},\hat{a}_{l}^{{\dagger}}]=\delta_{kl}$ and $[\hat{a}_{k},\hat{a}_{l}]=[\hat{a}_{k}^{{\dagger}},\hat{a}_{l}^{{\dagger}}]=0$, we obtain

[TABLE]

where we have denoted $D_{k}(t)=\sin(\Omega_{k}t)+i[1-\cos(\Omega_{k}t)]$.

The reduced density operator is obtained by tracing over the bath degrees of freedom as

[TABLE]

where on the last line we have used the thermal initial state for the bath, defined in Eq. (12). Here, the trace can be simplified with the unitary rotation $\hat{U}=e^{-i\sum_{k}\alpha_{k}(t)\hat{a}_{k}^{{\dagger}}\hat{a}_{k}}$ where $\alpha_{k}(t)$ are chosen such that the imaginary part of $D_{k}(t)\hat{a}_{k}^{{\dagger}}+D_{k}^{*}(t)\hat{a}_{k}$ is eliminated. We also note that the thermal-state density operator is diagonal in the eigenbasis of $\hat{a}_{k}^{{\dagger}}\hat{a}_{k}$ and, thus, unaffected by the rotation. We obtain

[TABLE]

where in the last equality, we have set $\alpha_{k}(t)=-\textrm{arg}[D_{k}(t)]$. As a consequence of the transformation, we can write the reduced density operator as

[TABLE]

Thus, we need to calculate $\langle e^{ic_{k}(\hat{a}_{k}^{{\dagger}}+\hat{a}_{k})}\rangle$ in the thermal state, where

[TABLE]

This can be carried out in the phase-space representation, where the Wigner function for the thermal density operator $\hat{\rho}_{k}^{\rm th}$ reads

[TABLE]

with the thermal averages $\langle\hat{q}_{k}^{2}\rangle=q_{k,0}^{2}{\rm coth}(\hbar\beta\Omega_{k}/2)$ and $\langle\hat{p}_{k}^{2}\rangle=\langle q_{k}^{2}\rangle\hbar^{2}/(4q_{k,0}^{2})$ and $q_{k,0}^{2}$ being the ground-state width of mode $k$.

As a consequence, one arrives with $\hat{q}_{k}=q_{k,0}(a_{k}^{\dagger}+a_{k})$ at

[TABLE]

Finally, we can write the elements of the reduced density matrix as

[TABLE]

where

[TABLE]

and we have recalled that

[TABLE]

where the latter equality holds for an ohmic bath with a second-order Drude cutoff at $\omega_{\rm c}$. We emphasize that these are identical relations to those obtained previously by Braun et al. (Braun et al., 2001).

In the main text, we show data for the excited-state occupation of the operator $\hat{\sigma}_{\rm z}$ given by

[TABLE]

where in the last equality we have assumed that initially $\rho_{\rm e}(0)=1$. We have denoted

[TABLE]

where the pointer states $|\pm\rangle$ are the eigenstates of the $\hat{\sigma}_{\rm x}$ operator obeying $\hat{\sigma}_{\rm x}|\pm\rangle=\pm|\pm\rangle$.

B.1.1 Very early behavior

For very early times ($\omega_{\rm c}t<1$), the argument in the exponential reduces into

[TABLE]

where in the last equality we have also assumed zero temperature. We thus observe that, contrary to the Fermi’s golden rule, the early dependence of the excited-state occupation is proportional to $t^{2}$. This implies that for times obeying $\omega_{\rm c}t<1$, we expect faster than exponential decay of the initial excited-state occupation. This feature is demonstrated in Fig. 1(b).

B.1.2 Thermal part

In general, the function $f(t)$ that determines the early decoherence can be expressed as a sum of the zero-temperature and finite-temperature parts as

[TABLE]

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

[TABLE]

This was obtained by first calculating the integral for the time derivative of $f_{\beta}(t)$ and then integrating the result in time.

B.1.3 Asymptotic behavior

We have calculated symbolically using maple the asymptotic behavior of the zero-temperature part. The result can be written as

[TABLE]

where $\gamma\approx 0.577\ldots$ is the Euler constant. We emphasize that for the experimentally relevant case with $(\hbar\beta)^{-1}<\omega_{\rm q}\ll\omega_{\rm c}$, the early 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. (42) and (41) is equal to Eq. (2).

B.2 Partition function approach for the steady state

The qubit occupation in the steady state can be calculated using the canonical partition function

[TABLE]

of the whole qubit-bath system. For clarity, we study the results in terms of the Kondo parameter $K=\kappa/(\pi\omega_{\rm q})$. If the system depends on a parameter $\lambda$, one can write

[TABLE]

and, consequently,

[TABLE]

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(\hbar\beta\omega_{\rm q}\hat{\sigma}_{\rm 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 Weiss (1999) $Z_{q}$ with the property

[TABLE]

for any parameter $\lambda$ which appears only in the system Hamiltonian ($\partial H_{I}/\partial\lambda=\partial H_{R}/\partial\lambda=0$). By setting $\lambda=\omega_{\rm q}$ in Eq. (45), one obtains

[TABLE]

where $\hat{\sigma}_{\rm z}=|\textrm{g}\rangle\langle\textrm{g}|-|\textrm{e}\rangle\langle\textrm{e}|$ is a Pauli operator of the bare ideal qubit. The reduced partition function of the dressed qubit can be expressed in the weak coupling regime as Weiss (1999)

[TABLE]

where

[TABLE]

$\psi(z)$ is the digamma function, and

[TABLE]

with $\Gamma(x)$ being the gamma function. The above result is equivalent to Eq. (4) and valid for all values of $\beta$. It has been derived using an 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. (47), we obtain Eq. (3):

[TABLE]

The expectation value comprises 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, $\partial\Omega/\partial\omega_{\rm q}$, of the expectation value is a measure of entanglement between the qubit and the bath, and can be written as $\partial\Omega/\partial\omega_{\rm q}=\partial\Omega/\partial\omega_{\rm eff}(\partial\omega_{\rm eff}/\partial\omega_{\rm q})$ where

[TABLE]

B.2.1 Lamb shift

If one neglects the factor $\partial\Omega/\partial\omega_{\rm q}$ in the expression in Eq. (52), one obtains

[TABLE]

where the derivative is with respect to the renormalized frequency $\Omega$ instead of the bare frequency $\omega_{\rm q}$ as in Eq. (47). This describes a qubit with only a frequency renormalization $\omega_{\rm q}\rightarrow\Omega$ which is what one would expect for a system experiencing only the Lamb shift and no entanglement with the bath.

In the literature Carmichael (1999), 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 ideal qubit can be written as

[TABLE]

where $\mathcal{P}$ stands for principal value. In the second equality, we have assumed zero temperature and $\omega_{\rm c}\gg\omega_{\rm q}$, and used the spectral density $J(\omega)=\pi K\omega\exp(-\omega/\omega_{\rm c})$ with an exponential cutoff. For Drude cutoff, the Euler constant $\gamma=0.577\ldots$ is replaced with $\frac{1}{2}$. On the other hand, one can make a linear expansion in $K$ for the renormalized qubit frequency given in Eq. (49). Note that at zero temperature and for $\omega_{\rm c}\gg\omega_{\rm q}$ the expression (49) reduces to $\omega_{\rm LS}$ only under the much stricter constraint $K\ln(\omega_{c}/\omega_{\rm q})\ll 1$. We also emphasize that for a purely harmonic system with natural frequency of $\omega_{\rm q}$, the Lamb shifted transition frequency is given by Silveri et al. (2019)

[TABLE]

Thus, contrary to the logarithmic divergence of the Lamb shift of a maximally anharmonic system in Eq. (56), the Lamb shift of a harmonic oscillator converges towards zero with increasing $\omega_{\rm c}$. Being only a weakly anharmonic system, we thus expect that the Lamb shift of the transmon is also small compared to that of a two-level system. This is also supported by our numerical data (not shown) which display a negligible shift of the Larmor frequency for $N=5$ if compared against the data for $N=2$ shown in the inset of Fig. 2(a).

In Fig. 4, we compare the renormalized frequency $\Omega$ in Eq. (49) with the Lamb-shifted qubit frequency $\omega_{\rm LS}$ in Eq. (56). As expected, we observe that the perturbative result follows closely the renormalized frequency in Eq. (49) for small values of the coupling constant $K$. At low temperatures with $\beta>5$, the deviations appear for $K\gtrsim 0.1$. The perturbative nature of $\omega_{\rm LS}$ is emphasized by the fact that it decreases without a bound and becomes negative at $K=[\ln(\omega_{\rm c}/\omega_{\rm q})-\gamma]^{-1}\approx 0.3$, where the numerical value has been calculated for the parameters used in Fig. 4. On the other hand, the renormalized frequency $\Omega$ approaches zero asymptotically.

However, one cannot obtain Eq. (55) as a limiting case to Eq. (52). This would require $\partial\Omega/\partial\omega_{\rm eff}\rightarrow 1$ and $\partial\omega_{\rm eff}/\partial\omega_{\rm q}\rightarrow 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 $\beta$ and $K$. Especially, when the temperature is zero, the excited-state occupation of the qubit in the steady state is given solely by the entanglement with the bath, as we will show in the following sections.

B.2.2 Comparison with numerically obtained Larmor frequency

We have compared the renormalized frequency $\Omega$ in Eq. (49) with the numerically obtained Larmor frequency $\omega_{0}$. In Fig. 5, we show the relative error

[TABLE]

The data shows that for $\kappa\lesssim 0.1\times\omega_{\rm q}$ the relative error is below 2%. This further justifies the interpretation of $\Omega$ as the renormalized transition frequency of the qubit in the weak-coupling limit. We also note that the data for $\kappa=0.1\times\omega_{\rm q}$ correspond to the decaying Larmor oscillations shown in the inset of Fig. 2(a).

B.2.3 Zero-temperature occupation

Here, we show that the coupling to the bath gives rise to a nonvanishing excited-state occupation of the ideal qubit in the steady state, even in the zero-temperature limit $\hbar\beta\omega_{\rm eff}\rightarrow\infty$. We assume a weak coupling ($K\ll 1$) and a high cutoff ($\omega_{\rm c}\gg\omega_{\rm q}$). With these approximations, we obtain in Eq. (52) that $\tanh(\hbar\beta\Omega/2)\approx 1$ and $\Omega\approx\omega_{\rm eff}$. Thus,

[TABLE]

We are interested in the steady-state occupation $\rho_{\rm e}^{\rm\infty}=(1-\langle\hat{\sigma}_{\rm z}\rangle)/2$ in the excited state $|\textrm{e}\rangle$ of the $\hat{\sigma}_{\rm z}$ operator. We obtain up to the first order in $K$ that

[TABLE]

It should be noted that for $K>0$, this result deviates from the $\rho_{\rm e}^{\infty}=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.

B.3 Perturbative treatment of the zero-temperature ground-state entanglement

We calculate perturbatively the reduced density operator of the qubit in the steady state in the zero-temperature limit. We assume that the system is described with the total Hamiltonian

[TABLE]

and obeys the Boltzmann distribution in the steady state. Thus at zero temperature, the whole qubit-bath system is in its ground state. Since the Jaynes–Cummings type terms $\propto\hat{\sigma}_{+}\hat{a}_{k}+\hat{\sigma}_{-}\hat{a}_{k}^{{\dagger}}$ conserve the occupation number, they do not affect the ground state and are, thus, neglected in the following. As a consequence, we can write the total Hamiltonian as

[TABLE]

where $\sigma_{-}=|\textrm{g}\rangle\langle\textrm{e}|$ is the annihilation operator of the qubit. This Hamiltonian can be approximately diagonalized up to the second-order in the couplings $g_{k}$ with the unitary transformation

[TABLE]

where $\hat{S}^{(k)}=[g_{k}/(\omega_{\rm q}+\Omega_{k})](\hat{a}_{k}^{{\dagger}}\hat{\sigma}_{+}-\hat{a}_{k}\hat{\sigma}_{-})$ diagonalize the interaction terms of the Hamiltonian in Eq. (64). As a result, we obtain the Hamiltonian

[TABLE]

We assume that in the transformed frame, the qubit-bath system is in a thermal state at zero temperature, i.e. in the ground state of the Hamiltonian in Eq. (66). Since the coupling terms in Eq. (66) conserve the occupation number, the ground state can be written as $|\textrm{g},0,0,\ldots\rangle$, where the first quantum number labels the state of the dressed qubit and the rest those of the dressed bath oscillators. The corresponding density operator of the qubit-bath system can be written as

[TABLE]

The bare qubit occupation can be obtained by transforming the ground-state density operator back to the laboratory frame as

[TABLE]

where $\hat{S}=\sum_{k}\hat{S}^{(k)}$, and the second equality holds up to the second order in the coupling constants $\{g_{k}\}$. The reduced density operator for the qubit is obtained by tracing over the bath, resulting in

[TABLE]

where

[TABLE]

is a measure of the entanglement between the qubit and the bath, i.e., hybridization of the reduced density operator of the qubit even in the zero-temperature limit. Above, we have used the definition for the mode spectral density given by

[TABLE]

Finally, let us compare this result with Eq. (62). For cutoff with $\omega_{\rm c}\gg\omega_{\rm q}$, the excited-state occupation is given by

[TABLE]

The difference between the constant terms is caused by the fact that Eq. (62) was derived using an exponential cutoff whereas here we used a Drude cutoff similar to our numerical simulations. Thus, we observe that the $\partial\Omega/\partial\omega_{\rm q}$ factor in Eq. (52) arises due to entanglement. As a consequence, the entanglement is the dominating source of initialization error at low temperatures.

Appendix C Effect of intrinsic dissipation

In addition to the low-temperature engineered environment described above, the transmon qubit is also coupled to an intrinsic and uncontrollable environment characterized by the zero-temperature dissipation rate $\gamma$ and temperature $T_{\rm i}$. During an initialization protocol, we assume that $\kappa\gg\gamma$. Consequently, the intrinsic environment has a negligible effect on the steady-state occupation of the transmon.

We show this by studying a harmonic oscillator that is coupled weakly to two thermal baths. We assume weak coupling and, as a consequence, the quantum dynamics of the reduced density operator is governed by the master equation

[TABLE]

where $N_{\rm ee}$ and $N_{\rm i}$ are the Bose–Einstein occupations of the engineered and intrinsic baths, respectively. We denote the occupation probabilities for the bare oscillator eigenstates $|n\rangle$ with $P_{n}=\langle n|\hat{\rho}|n\rangle$. Using the master equation (73), one obtains

[TABLE]

We assume that the steady state ($\dot{P}_{n}=0$) is given by the thermal occupation with

[TABLE]

where $N$ is the effective Bose–Einstein occupation for the harmonic oscillator interacting with two independent baths. After straightforward algebra, we obtain

[TABLE]

Clearly in the limit $\kappa\gg\gamma$ and for a relatively low intrinsic thermal occupation $N_{\rm i}$ of the oscillator, one obtains that $N\approx N_{\rm ee}$. In the main text, we assume that this holds also for relatively large dissipation rate $\kappa\sim 0.1\times\omega_{\rm q}$. As a consequence, we neglect the intrinsic bath in order to clarify our discussions and to emphasize our main message.

Appendix D 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 Weiss (1999). In such situations, one has to rely on more accurate methods, such as the formally exact Feynman–Vernon path integral formalism.

D.1 Stochastic Liouville–von Neumann equation

For bilinear coupling and the qubit-bath system starting from a factorized initial state given in Eq. (11), one can show Stockburger and Grabert (2002) 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

[TABLE]

The SLN equation comprises a deterministic coherent part given by the first term on the right-hand side of the equation, followed by the stochastic dissipative part, the dynamical properties of which are set by the complex-valued random variables $\xi$ and $\nu$. The correlations between the qubit and the bath are encoded into the correlations of the noise terms which follow the equations

[TABLE]

where $\Theta(t)$ is the Heaviside step function and

[TABLE]

is the autocorrelation function of the bath with $\hat{\zeta}=\sum_{k}g_{k}(\hat{a}_{k}^{{\dagger}}+\hat{a}_{k})$. We note that the correlations $\langle\xi(t)\xi^{*}(t^{\prime})\rangle$, $\langle\xi(t)\nu^{*}(t^{\prime})\rangle$, and $\langle\nu(t)\nu^{*}(t^{\prime})\rangle$ are not fixed by the bath correlation function, and can be thus chosen to optimize the efficiency of the numerical realization.

D.2 Stochastic Liouville equation with dissipation

In the case of ohmic dissipation with a high Drude cutoff frequency $\omega_{\rm c}\gg\omega_{\rm q}$, the path integral formalism can be reduced into the form of so-called stochastic Liouville equation with dissipation (SLED) Stockburger (1999b); Wiedmann et al. (2016)

[TABLE]

The above SLED has a stochastic part characterized by a single real-valued noise term $\xi(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

[TABLE]

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.

Appendix E Implementation of the SLN and SLED methods

The SLN equation (77) and the SLED in Eq. (82) 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,\ldots,(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 $\kappa^{-1}\ll(N-1)h$.

E.1 Generation of correlated-noise samples

The numerical implementation of the correlated-noise samples is the main difference from the corresponding Lindblad algorithm. Here, we describe a numerical scheme for the generation of such samples obeying the correlation functions (78)–(80). 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. (81). In the case of the SLED, one only needs to consider Eq. (78) where the real part of the bath correlation function is defined in Eq. (83). We give here the method of generating the correlated complex-valued noise terms $\xi$ and $\nu$ 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 $\xi$ noise into two parts as $\xi(t)=\xi_{\textrm{r}}(t)+\xi_{\rm c}(t)$, where $\xi_{\rm r}(t)$ is assumed real and

[TABLE]

and $\langle\xi_{\rm c}(t)\xi_{\rm c}(t^{\prime})\rangle=\langle\xi_{\rm r}(t)\xi_{\rm c}(t^{\prime})\rangle=\langle\xi_{\rm r}(t)\nu(t^{\prime})\rangle=0$. Moreover, we denote $\xi_{\rm c}(t)=\xi_{\rm c}^{\rm R}(t)+i\xi_{\rm c}^{\rm I}(t)$ and $\nu(t)=\nu^{\rm R}(t)+i\nu^{\rm I}(t)$, where the terms $\xi_{\rm c}^{\rm R,I}$ and $\nu^{\rm 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

[TABLE]

where we have defined the Fourier transformation of a function $f(t)$ as

[TABLE]

The window functions are determined by the correlation functions in Eqs. (84) and (85), and can be expressed as

[TABLE]

where $L_{\rm i}(\omega)$ is the Fourier transform of Im[$L(t)$]. For odd spectral densities, we can write

[TABLE]

where the second equality in the last equation has been written for the ohmic spectral density $J(\omega)$ defined in Eq. (71).

Each noise term can be generated numerically by the following protocol:

Produce an array of $N$ independent Gaussian variables $\{x(t_{\ell})\}$ corresponding to the grid points $t_{\ell}=\ell h$ with $\ell=0,\ldots,N-1$.

• 2.

Use fast Fourier transformation on $\{x(t_{\ell})\}$ to obtain $\{x(\Omega_{k})\}$ where $\{\Omega_{k}=k2\pi/(Nh)\}$ define a grid in the frequency space.

• 3.

Take the inverse Fourier transformation of $W(\Omega_{k})x(\Omega_{k})$ to obtain the discretized samples of Eqs. (86)–(90).

In the case of SLED, one needs to do this procedure only once, for $\xi_{\rm r}$, whereas for the SLN equation all five real-valued random variables $\xi_{\rm r},\xi_{\rm c}^{\rm R,I},$ and $\nu^{\rm R,I}$ are needed.

E.2 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 $\dot{\rho}_{\rm S}=\mathcal{L}\rho_{\rm S}$ where $\mathcal{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 1, 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 1 for the data relevant for Figs. 1 and 2.