Efficient protocol for qubit initialization with a tunable environment

We propose an efficient qubit initialization protocol based on a dissipative environment that can be dynamically adjusted. Here, the qubit is coupled to a thermal bath through a tunable harmonic oscillator. On-demand initialization is achieved by sweeping the oscillator rapidly into resonance with the qubit. This resonant coupling with the engineered environment induces fast relaxation to the ground state of the system, and a consecutive rapid sweep back to off resonance guarantees weak excess dissipation during quantum computations. We solve the corresponding quantum dynamics using a Markovian master equation for the reduced density operator of the qubit-bath system. This allows us to optimize the parameters and the initialization protocol for the qubit. Our analytical calculations show that the ground-state occupation of our system is well protected during the fast sweeps of the environmental coupling and, consequently, we obtain an estimate for the duration of our protocol by solving the transition rates between the low-energy eigenstates with the Jacobian diagonalization method. Our results suggest that the current experimental state of the art for the initialization speed of superconducting qubits at a given fidelity can be considerably improved. The ability to accurately reset qubits is one of the fundamental criteria in implementing quantum computers. Jani Tuorila and co-workers from Aalto University (Finland) propose a qubit initialization protocol using a dynamically adjustable environment. The scheme differs from the conventional methods which reach for a given qubit state with microwave cooling or projective measurements of the qubit. The protocol relies on enhancing the dissipation to the qubit ground state by tuning an engineered environment on resonance with the qubit. The new theoretical study suggests that the experimental benchmark for the initialization speed at a given precision can be increased with this method almost by an order of magnitude. Such an improvement would present a significant step towards meeting the stringent operational conditions of a working largescale quantum computer.


INTRODUCTION
Preparation of a qubit into a well-defined initial state is one of the key requirements for any quantum computational algorithm. 1,2 The conventional passive initialization protocol relies on the relaxation of the qubit to a thermal state determined by the residual coupling to the environment. This protocol is inherently slow because the relaxation rate has to be minimized to decrease the probability of errors in a coherent quantum computation.
In the context of error-free quantum computing, a long initialization time would not present a problem since the quantum register has to be prepared only once in the beginning of the computation. However, realistic quantum computational devices also suffer from gate errors which have to be corrected with quantum-error-correction codes. 3 Such codes rely on logical qubits which consist of several physical qubits. The codes initiate from a predetermined state for the physical qubits and are being constantly executed during a computation. They also have strict requirements for the initialization and gate error thresholds for individual qubits, of the order of 10 −5 for the conventional concatenated codes. 4,5 In the more refined topological quantum error correction codes, 6 the logical error can be suppressed with stabilizing measurements, which increase the thresholds up to 10 −2 for the physical qubit operations and lead to improved protection of quantum information during the computation. Nevertheless, stabilizer codes such as surface 7,8 and color 9, 10 codes still require a continuous supply of measurement qubits in a known low-entropy state. Thus, initialization time is also an issue in large-scale quantum computing.
Fast and accurate qubit initialization remains a technological challenge in the superconducting qubit implementations that have shown great potential for scalability. 8 Major developments have been made with active protocols, such as initialization by successive projective measurements, [11][12][13] by Purcell-filtered cavity, 14 or by cooling with a coherent microwave drive. [15][16][17][18] The typical figure of merit of a protocol is the time τ 10 required for a qubit excitation to decay by an order of magnitude. In the reported experiments, 14 a thermal-equilibrium fidelity of 99.9% with τ 10 = 40 ns has been obtained with a protocol for a Purcellfiltered superconducting qubit with transition frequency ω q /(2π) = 5. 16 GHz and intrinsic relaxation time T 1 = 540 ns. However, the method uses a tuned qubit that is not desirable since it reserves a broad frequency band, and a high ground-state fidelity using a long-lived qubit remains to be demonstrated with this method. An initialization protocol based on coherent driving has reached a 99.5% ground-state fidelity with τ 10 = 1.4 μs, ω q /(2π) = 5 GHz, and T 1 = 37 μs, 17 without relying on qubit tuning. Even though this protocol meets the error thresholds of the topological codes, large-scale quantum computing calls for shorter initialization times. Furthermore, initialization fidelities greater than 99.99% are preferred in large-scale computations to reduce the number of physical qubits needed for a logical qubit.
Qubit initialization can also be achieved with feedback-based schemes, 11,[19][20][21][22] which rely on a single-shot measurement of the qubit state. In recent experiments, 19,20 the discrimination of the qubit state is predominantly limited by T 1 , with the current benchmark for the average assignment fidelity being 99.8% in 700-ns read-out time. 19 However, an increase of the fidelity by each order of magnitude requires yet another reduction of the read-out time or an increase of T 1 by an order of magnitude, which still remains an experimental challenge. Furthermore, the fidelity and speed of this method may not exceed that of the quantum measurement, and hence significantly contribute to the total error budget. Therefore, the search for a more efficient initialization scheme is of great importance.
In solid-state systems, an initialization protocol can potentially be realized by strongly coupling the qubit to a cold dissipative reservoir for fast initialization and by switching off the coupling for the actual computation. 23,24 If the bath has a lower effective temperature than the qubit, the ground-state fidelity of the qubit is increased. This type of set-up is a part of environmental quantum-state engineering by dissipation, where one aims to drive the system into a desired steady state by using a carefully tailored environment. [25][26][27][28][29][30] Such engineering has already been used in generation of coherent superposition states, [31][32][33] in creation of entanglement, [34][35][36] and in simulations of open quantum systems. 37,38 In this paper, we focus on a ground-state initialization proposal, 39,40 where a superconducting qubit and a lowtemperature resistive bath are coupled indirectly through two resonators as shown in Fig. 1. The effective inductance of the left resonator, which is capacitively coupled to the bath, can be dynamically adjusted, allowing control over its bare resonance frequency. If the left resonator is sufficiently detuned from the qubit, it shunts the noise of the resistive bath at the qubit frequency, and hence the decoherence of the qubit is dictated by its slow intrinsic relaxation rate. If the left resonator is in resonance with the qubit, the qubit couples strongly to the bath leading to an increased relaxation rate. We use two resonators because a single resonator cannot be easily detuned from the qubit by several linewidths due to the required strong coupling with the cold bath and because the second resonator, coupled directly to the qubit, can be used for qubit readout using the standard dispersive methods. 41 Previous calculations of static transition rates in this scenario indicate that the lifetime of the qubit can be controlled over several orders of magnitude. 39,40 However, an actual initialization protocol and its dynamics, speed, and fidelity have not been reported. Here, we develop a qubit initialization protocol and solve its corresponding quantum dynamics by using a Markovian master equation. We show analytically that the ground state is protected during the sweeps of the left resonator to and from the resonance, implying that the speed of the protocol is determined by the strong resonant coupling with the dissipative bath. At the resonance, we find an approximative analytic solution for the lowenergy eigenproblem using Jacobian diagonalization. 42,43 It yields a useful lower bound for the duration of the protocol at a given fidelity. Optimization of the physical circuit parameters suggests that, with present-day technology, the current benchmarks for fidelity and protocol speed 17 can be considerably improved with the help of our protocol.

System
The above-discussed tunable-environment qubit can be conveniently studied using a lumped-element circuit model shown in Fig. 1a. Here, the superconducting qubit with the transition energy ħω q is coupled to a bosonic heat bath through two LC resonators. Both resonators are formed by a lumped capacitance C k and an inductance L k , where k = L,R refer to the left and right resonators, respectively. The bath arises from the resistance R at temperature T. The left resonator is coupled directly to the bath and to the right resonator through capacitances C E and C c , respectively. The inductance L L (t) of the left LC resonator is made tunable using a SQUID, the Josephson inductance of which is controlled by an external magnetic flux. As a consequence, the bare angular frequency of the left resonator ω L ðtÞ ¼ 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffi L L ðtÞC L p can also be adjusted with the external flux. The right resonator has a fixed bare angular frequency ω R ¼ 1= ffiffiffiffiffiffiffiffiffi L R C R p and is coupled to the qubit through the capacitance C q .
The Hamiltonian of the circuit can be written as 40 HðtÞ ¼Ĥ S ðtÞ þĤ E þĤ I ðtÞ; ð1Þ where the three terms describe the system, the resistive environment, and their interaction, respectively. The system Hamiltonian can be expressed aŝ H S ðtÞ ¼Ĥ 0 ðtÞ þĤ 1 ðtÞ; wherê frequency ω q and intrinsic relaxation rate κ I q is indirectly coupled to a thermal bath (brown), formed by a resistor R, through two LC resonators. By tuning the inductance L L (t) of the left resonator (orange), one can tune its bare resonance frequency and coupling strengths g LR (t) and κ L (t) with the right resonator (magenta) and the bath (temperature T), respectively. The right resonator has a bare angular frequency of ω R and is coupled to the qubit and an intrinsic bath with coupling strengths g Rq and κ I R . (b) Analogous cavity quantum electrodynamics (QED) set-up where a two-level atom is coupled to a thermal bath through two optical cavities. The coupling to the thermal bath is controlled by tuning the length ' L ðtÞ of the left cavity and Above,â L ,â R , andσ À are the annihilation operators for the left and right resonators and the qubit, respectively. In the following, we denote the eigenstates of the unperturbed Hamiltonian (3) with n L ; n R ; n q , where the occupation numbers of the left and right resonators can have values n L , n R = 0, 1, 2, … and that of the qubit assumes values n q = g,e, which stand for the ground and excited state, respectively. When the rotating wave approximation (RWA) applied in Eq. (5) is accurate, the total occupation number N = n L + n R + n q is a conserved quantity during a unitary time evolution. The HamiltonianĤ S ðtÞ describes a tripartite system formed by two harmonic resonators and a qubit. The right resonator is coupled bi-linearly to the left resonator and to the qubit with the respective coupling is the resonant coupling strength between the left and right resonators, C J is the capacitance of the Josephson junction, and E J and E C = e 2 /[2(C q + C J )] are the Josephson coupling energy and the charging energy per electron for the superconducting island, respectively. Coupling between the qubit and the left resonator is mediated by the right resonator and is, therefore, of second order in coupling frequencies g 0 LR and g Rq . Consequently, the right resonator acts as an additional filter for the thermal noise of the bath. In our analytic considerations, we apply the RWA for both of the coupling terms, cf. Eq. (5). The interaction with the bath is also bi-linear and described bŷ is the voltage over the left resonator, and δV res describes the voltage fluctuations over the resistor. The resistor HamiltonianĤ E is given by that of a bosonic thermal bath 44 . We do not express it explicitly since we are only interested in the transition rates, which are determined in thermal equilibrium by the Johnson-Nyquist spectrum of the voltage fluctuations: Furthermore, we neglect any filtering of this spectrum owing to the resistor itself by assuming that the frequencies relevant for the system obey ω ( 1/(RC R ) 40 .
In typical realizations of superconducting qubits, the temperature of the circuitry is below 100 mK and the resonator and qubit frequencies lie in the 5-10-GHz range. As a consequence, thermal occupations of the one and two excitation states for a detuned system are roughly below 10% and 1% relative to the groundstate occupation, respectively. Furthermore, during the protocol, non-adiabatic transitions between instantaneous eigenstates with different total occupation numbers N are suppressed owing to negligible matrix elements for such transitions. This allows the development of an approximative analytic model in the subspace spanned by the four lowest-energy instantaneous eigenstates of the Hamiltonian (2). We emphasize that the following discussion is not specific to the lumped-element model or superconducting qubits, but can be used rather generally for qubits with indirect adjustable coupling to a thermal bath, as shown in ref. 40 with the distributed circuit elements which are frequently used in contemporary circuit QED.
where ω mn (t) = ω n (t) − ω m (t), and ω k (t) are the eigenfrequencies corresponding to the eigenstates Ψ k ðtÞ j iof the HamiltonianĤ S ðtÞ, i.e.,Ĥ S ðtÞ Ψ k ðtÞ j i¼ hω k ðtÞ Ψ k ðtÞ j i, We thus observe that positive-frequency fluctuations in the environment cause emission in the small system. 45 Note, that in Eq. (8) we express the transition rates in units of Γ 0 , which equals the bare zero-temperature left-resonator transition rate for ω L (t) = ω R . Clearly, we can maximize the transition rates by maximizing C E and R with respect to C R and Z R , respectively. We further note that in principle C E should be added to C L to obtain the effective left-resonator capacitance, but if C E ( C L its effects on the eigenstates and eigenfrequencies of the system are negligible. 40 In Methods, we analytically solve the transition rates for the first three excited eigenstates employing the RWA. We first truncate the Hamiltonian (2) to the low-energy subspace spanned by the ground state and single-excitation states. Then, we diagonalize the static Jaynes-Cummings coupling in Eq. (5) between the right resonator and the qubit. As a result, the first three eigenstates of the Hamiltonian (2) are the eigenstates of the effective Hamiltonian matrix which is written in the basis formed by |1, 0, g〉 and the Jaynes-Cummings eigenstates | ± 〉 with the corresponding eigenfrequencies ω ± ¼ ω av Rq ± Ω Rq . We use the short-hand nota- Above, we have G Rq = g Rq and, thus, λ Rq is the usual perturbation parameter of the Jaynes-Cummings model that characterizes the strength of the coupling. Similarly, the coupling strengths between the left resonator and the states | ± 〉 are q . We treat the coupling to the left resonator perturbatively using Jacobian diagonalization. 42,43 Jacobian diagonalization of a hermitian matrix consists of a sequence of similarity transformations designed to annihilate one off-diagonal matrix element at a time. Optimal convergence is achieved by forming the subspace at each step of the basis states with the strongest coupling. Depending on the value of the control parameter ω L (t), we have either λ L+ (t) ≥ λL− (t) or λ L+ (t) < λ L− (t) and, thus, the Jacobian diagonalization is most efficiently done in a temporally piecewise manner. As a result, we obtain the transition rates where ω 1 ðtÞ ¼ ω av LÀ ðtÞ À Ω LÀ ðtÞ, ω 2 ðtÞ ¼ ω av L ± ðtÞ ∓ Ω L ± ðtÞ, and ω 3 ðtÞ ¼ ω av Lþ ðtÞ þ Ω Lþ ðtÞ. For transition rate (11), we use the top signs when λ L+ (t) ≥ λ L− (t), and bottom signs when λ L+ (t) < λ L− (t). Note that the transition rates between the first three excited states are zero in the RWA due to the selection rules for our environmental coupling term. However, we include all transition rates to the numerical solution, within the truncation. Furthermore, the principle of detailed balance Γ mn ¼ exp À hω mn = k B T ð Þ ½ Γ nm holds, which implies that the excitation rates are strongly suppressed in the low-temperature limit. In addition to the engineered environment described by the resistor R, the qubit and the right resonator typically dissipate energy to their own intrinsic environments with the transition rates κ I q and κ I R , respectively. We compare the analytic eigenfrequencies and transition rates with the corresponding numerical results in Fig. 2. For large detuning g 0 LR ; g Rq ( ω R À ω q À Á , we find that already the first diagonalization step in the Jacobian diagonalization procedure results in a very good agreement with the numerically obtained frequencies and rates. We also observe that when the left resonator is in resonance with either the right resonator or the qubit [δ L ± (t) ≈ 0], the transition rates from the resonant eigenstates to the ground state are equal. We take advantage of this fact in our protocol below. In the vicinity of the resonances, there exist regions where the transition rates from the nearly resonant states change by several orders of magnitude. The widths of these regions are directly proportional to the coupling term g 0 LR , and also to g Rq for the case of ω L (t) = ω − .

Protocol
Our proposed initialization protocol is depicted in Fig. 2c and proceeds as follows. In the beginning of the protocol (t = 0), the parameters of the set-up follow the hierarchy In particular, the qubit frequency is chosen to be the smallest in order to minimize the effects of the possible multi-photon processes in the qubit caused by the right-resonator at any stage of the protocol. For t = 0 and t = τ, where τ is the duration of the protocol, the coupling to the engineered environment should be minimal so that the intrinsic sources of dissipation are dominating the qubit decoherence. Therefore, we refer to the bare leftresonator frequency ω L (0) = ω L (τ) as the ''off'' state of our set-up, and choose its value such that g LR ð0Þ ( ω L ð0Þ À ω R . This guarantees that (see Figs. 1 and 2) since the widths of the regions for enhanced transition rates are proportional to g 0 LR . This way the added dissipative channel does not create excess decoherence for the qubit. We also choose g Rq ( ω R À ω q , which implies only weak hybridization between  Table 1. c Initialization protocol in terms of the control parameter ω L (t) (see text for details) Protocol for qubit initialization J Tuorila et al.
the qubit and the right-resonator. In our numerical simulations, we use the values shown in Table 1 for the relevant parameters.
In the first stage of the protocol, the left resonator is swept fast to resonance with the effective right resonator at Rq =δ Rq, which is only slightly hybridized with the qubit due to the dispersive coupling. As a consequence, the right resonator becomes strongly coupled to the cold bath resulting in an increase in its relaxation rate to the ground state by orders of magnitude (see Fig. 2). This operation point (denoted in Fig. 2 with ''op 1'') guarantees equal relaxation rates for both resonators, which is important as the relative occupations are typically not known in the beginning of the protocol and, also, because non-adiabatic transfer of occupation between the resonators can occur during the fast sweep. Any occupation in either resonator is then dissipated to the resistor during t : t 1 → t 2 . The wait time Δt 2 depends on the required protocol error α = 1 − P 0 (τ), where we use the notation Δt i = t i − t i−1 for the relevant protocol time intervals with t 0 = 0 and t 5 = τ, and define P 0 (τ) as the desired ground-state occupation and the end of the protocol.
After the first relaxation step (t ≥ t 2 ), the resonators are in the ground state within the protocol error α. Then, the left resonator is swept into resonance with the effective qubit at At this operation point (denoted in Fig. 2 with ''op 2''), the relaxation rate of the qubit is increased by several orders of magnitude and is equal to the left-resonator rate. The latter is important since at this operation point, the qubit is also hybridized with the left resonator and, consequently, any initial occupation in the qubit is partly transferred to the left resonator during the sweep. By waiting for t : t 3 → t 4 , the hybridized qubit and left resonator dissipate their energy to the resistor. The second wait time Δt 4 is also determined by the desired value of α. Finally, the left-resonator frequency is swept back to its high initial value ω L (τ) = ω L (0). In addition to the wait times Δt 2 and Δt 4 , the protocol duration τ ¼ P 5 i¼1 Δt i is determined by the three sweep times Δt 1 , Δt 3, and Δt 5 . For simplicity, we will assume in our discussions that the sweep times are equal, i.e., t s = Δt 1 = Δt 3 = Δt 5 .
The speed of the protocol for a given fidelity should be maximized for efficient use in quantum information processing. In general, a good initialization protocol should have κ I q τ ( 1. This way one can perform multiple initializations during a coherent quantum computation. In our protocol, this requires the minimization of the combined duration of the three sweeps of the control parameter ω L (t) and the relaxation intervals, during which the actual initialization takes place. The length of the relaxation intervals is set by the relaxation rates and the desired fidelity, implying that after they are optimized, the duration of the protocol can be shortened only by faster sweeping. However, it is well known that fast changes in parameters can induce nonadiabatic transitions between the instantaneous eigenstates of the system. [46][47][48][49] Our requirements for the optimal operation points (ω L ≈ ω ± ) guarantee that our protocol is robust with respect to changes between the relative occupations of the instantaneous eigenstates during the first two sweeps. After the second relaxation process, the system lies in the ground state within the desired error α. Thus, the essential sweep is the last one starting from the qubit resonance (ω L (t 4 ) = ω − ) and the ground state |Ψ 0 (t 4 )〉, and ending to the far off-resonant ground state |Ψ 0 (τ)〉. The relevant question is the following: how much of the ground state is excited during the final fast sweep?
We can calculate the transition probabilities P mn ≡ P m → n (τ ; t 4 ) = |〈Ψ n (τ)|Ψ m (t 4 )〉| 2 between the low-energy eigenstates in the sudden approximation using the RWA results obtained in the previous section. However, since we are interested in the groundstate sweep fidelity P S = P 00 and since in the RWA the ground state is unaffected by the change of parameters, we have to include the contributions arising from the counter-rotating terms. In Methods, we derive the worst-case estimate for the ground-state sweep fidelity in the sudden approximation. We obtain where we assume that the sweep duration Δt 5 → 0. Thus, the deviation from the perfect sweep fidelity is given by the difference between the perturbation parameters g LR (t)/[ω L (t) + ω R ] before and after the sweep. Order-of-magnitude estimates for typical superconducting circuit parameters give g LR /(2π)~100 MHz, ω R /(2π),ω q /(2π)~10 GHz and, accordingly, we have that the deviation from full fidelity is 1 − P S~[ g LR /(ω q + ω R )] 2~1 0 −4 . Typically, the deviation is a couple of orders of magnitude smaller since the difference between the perturbation parameters is very small. We can, therefore, assume in our analytic considerations that the ground state is well protected during the sweeps of parameters in the Hamiltonian, and that the protocol duration can be estimated solely in terms of the relaxation rates of the static stages of the protocol, i.e., τ ≈ Δt 2 + Δt 4 . In the experimental implementation of the protocol, any cross coupling between the qubit and the flux line, used to adjust the left-resonator frequency, should be considered together with excitations of the SQUID. However, these issues seem not to significantly affect the achievable speed of our initialization protocol. In Fig. 3, we compare the analytically obtained sweep fidelity from Eq. (15) with the sweep fidelity obtained by solving the instantaneous ground state |Ψ 0 (t)〉 numerically. We observe that if the RWA is valid, i.e., for g 0 LR ( ω R , the analytic result closely follows the numerical solution. For increasingly strong coupling, the second-order perturbation theory becomes insufficient which is visible as a large deviation between the analytic and numerical results. Even for g 0 LR ¼ 0:1ω R , however, the deviation from the perfect fidelity is of the order of 10 −7 , which indicates that the effects of the fast sweep on the ground-state fidelity of the protocol can be neglected. Markovian master equation for the time-dependent system In order to obtain more quantitative understanding of the protocol, we study its dynamics with the help of a Markovian master equation for the reduced system density operator ρ S ðtÞ ¼ Tr Eρ ðtÞ f g, where the total density operatorρðtÞ obeys GHz ω q /(2π) 9.5 GHz ω R /(2π) 10.0 GHz ω L (0)/(2π) 11.5 GHz Γ 0 3.1 × 10 −2 10 9 s −1 λ Rq 0.14 We also show the circuit parameters of the lumped-element model that produce the used angular frequencies Protocol for qubit initialization J Tuorila et al.
the von Neumann equation whereĤðtÞ is defined in Eq. (1). We first diagonalize the system HamiltonianĤ S ðtÞ, as defined in Eq. (2), in a time-independent basis {|n〉} with the time-dependent unitary transformation DðtÞ ¼ P n Ψ n ðtÞ j i n h j. After the transformation, the timeevolution ofρ 0 ðtÞ ¼D † ðtÞρðtÞDðtÞ is governed by the effective system Hamiltonian where the latter term causes non-adiabatic transitions. Such term always appears if one wishes to preserve the form of the von Neumann equation in a time-dependent unitary transformation. After the transformation, the derivation of the master equation proceeds in a conventional manner [50][51][52] : We assume that the initial state of the total system is uncorrelated, i.e.,ρð0Þ 1 ρ S ð0Þ ρ E ð0Þ, and that the bath is in a thermal state, described byρ E ð0Þ, throughout the temporal evolution. We consider only weak coupling to the environment and apply the standard Born and Markov approximations in the interaction picture, and subsequently trace over the environmental degrees of freedom.
In the Markov approximation, the correlation time of the environment is negligibly short and we can express the master equation forρ ′ S ðtÞ in the secular approximation as Above, we have defined the ladder operators between the static basis states asπ nm ¼ n j i m h j. The instantaneous transition rates Γ mn (t) are identical to those obtained with Fermi's golden rule in Eq. (8). 53 Similar to the case of a static Hamiltonian, 50, 51 the environmental decoherence is included in the Lindblad terms on the last three rows of the master equation. The first term models the transitions between the adiabatic states that dissipate energy to the environment, the second term represents absorption from the environment and the last term induces dephasing of the adiabatic states. We note that the secular approximation made above is justified if the relaxation rates are small compared to the minimum separation between the eigenfrequencies of the system, 50 i.e., we have Γ mn ðtÞ ( min i≠j ω i ðtÞ À ω j ðtÞ : We note that the adiabatic master equation above can be improved by making successive diagonalizing transformations for H eff ðtÞ. As a result, the master equation is represented in the basis of the so-called superadiabatic states, the time-dependence of which is typically suppressed after each diagonalizing transformation. Such adiabatic renormalization was first described for a general time-dependent quantum system by Berry, 54 and later applied to studies of dissipation in driven superconducting qubits. [55][56][57] The lowest-order superadiabatic correction was studied in refs. 58, 59. Protocol speed Let us study the decay of an excitation in our system by solving the master equation (18) for an initial state spanned by the lowenergy adiabatic states as Ψð0Þ j i¼ P 4 n¼0 a n ð0Þ Ψ n ð0Þ j i. We assume that the environment is so cold that thermal excitations of the system are negligible. This guarantees that the quantum stateρ S ðtÞ of the system remains in the low-energy subspace during the temporal evolution. In the beginning of a realistic initialization procedure, we may have no knowledge on the distribution of the occupations P n (0) = |a n (0)| 2 . Therefore, we choose the operation points of our protocol such that

Fig. 3
Ground-state sudden sweep fidelity 1 − P S as a function of the control parameter value ω L (τ) at the end of the protocol. We compare the analytically obtained sweep fidelity in Eq. (15) with that obtained by finding the ground state |Ψ 0 (t)〉 numerically. We show the data for a g 0 LR ¼ 0:001ω R ; b g 0 LR ¼ 0:0074ω R (value used in the simulations in Figs. 2, 5, and 6); and c g LR = 0.1ω R . The qubit frequency is ω q = 0.95ω R and t 2 t 3 ; t 4 ½ : which is the case for ω L ðt 2 Þ % ω þ % ω R þ g 2 Rq =δ Rq and ω L t 4 ð Þ % ω À % ω q À g 2 Rq =δ Rq , respectively. Before discussing the numerical results, we present a simple analytic estimate for the initialization fidelity. The general form for the excited-state occupations can be solved from the master equation (18) and written as By neglecting the small contributions of the fast sweeps, these can be readily written in terms of Eqs. (20) and (21). We find that the deviation from the perfect fidelity, αðτÞ ¼ 1 À P 0 ðτÞ ¼ P 3 i¼1 P i ðτÞ, can be written as αðτÞ αð0Þ where the relative initial occupations of the excited states have been defined as P i ð0Þ ¼ P i ð0Þ= P 3 j¼1 P j ð0Þ for i = 1, 2, 3. The tolerable deviation α(τ) from perfect fidelity at the end of the protocol is fixed in the beginning of the protocol. Since we do not know the relative occupations in the beginning of the initialization protocol, we have to wait the times Δt 2 ≈ ln[α(0)/α(τ)]/Γ 2 and Δt 4 ≈ ln[α(0)/α(τ)]/Γ 1 so that any excitation in the system is decayed down to the desired accuracy. As a consequence, we obtain an upper bound for the total wait time τ ¼ Δt 2 þ Δt 4 of the protocol as τ log 10 αð0Þ=αðτÞ ½ 2lnð10Þ Above, we denote α(τ) = 10 −β α(0) where β ≥ 0, and τ 10 = 2ln(10) [1 + (ω R /ω q ) 2 ]/Γ 0 . For ω R ≳ ω q , we have that τ 10 ≈ 4ln(10)/Γ 0, which sets the time scale for the decrease of α(τ) by an order of magnitude. We note that the secular approximation used in the derivation of the master equation requires through Eq. (19) that where the equality holds since in our case 0<g Rq ( δ Rq and ω q < ω R . Additionally, we required that in the off state of the protocol, the intrinsic rates dominate over those of the engineered environment, i.e., Γ 10 ð0Þ ¼ α q κ I q and Γ 20 ð0Þ ¼ α R κ I R , where α q ; α R ( 1. We show in Methods that these lead to the ( 1 and δ Lq (0) = ω L − ω q is the detuning between the left-resonator and the qubit in the off state of the protocol. This sets a lower bound for the decay time: where the numerical estimate is made for typical superconducting circuit parameters ω R /(2π) = 10 GHz, ω q /(2π) = 9.5 GHz, δ Lq (0)/(2π) = 2 GHz, and κ I q ¼ 10 4 s −1 . We have used β = 0.5 and confirmed with a classical calculation of the transition rate Γ 2 that the secular approximation holds within a relative error of 4%. We also set α q = 0.1. The lower bound for the decay time above would represent a significant improvement to that of the current experimental benchmark for qubit ground-state initialization protocol. 17 By choosing the off state infinitely far from the qubit (δ Lq → ∞), the decay time can be improved to τ 10 ≳ 90 ns. In the following, we set α(0) = 1 in order to obtain an estimate for the total wait time for initialization of a qubit excitation. For example, if one wishes to obtain a ground state fidelity of α(τ) = 10 −3 with our protocol, the wait time should be of the order of τ ¼ 3τ 10 .
Qubit temperature Above, we calculated the deviation α(τ) from the perfect fidelity by neglecting the intrinsic dissipation in the qubit. However, in addition to the temperature T of the resistor, the fidelity is reduced by the intrinsic temperature T q of the qubit. We estimate this effect by first defining the effective qubit temperature T eff in terms of the excited-state occupation P q ex of the qubit at the end of the protocol as By using the principle of detailed balance, we can solve P q ex in our four-state model from where P q ex ¼ 1 2 P 1 þ P 2 ð Þ and we have included the coupling κ I q between the qubit and its intrinsic environment. The occupation at the end of the protocol can be calculated at the second operation point, where we use Γ 10 = Γ 20 = (ω q /ω R ) 2 Γ 0 /2 = Γ 1 , Γ 30 = 0, and the detailed balance relation for the excitation and absorption rates. Above, we have neglected the intrinsic excitation of the right resonator. We note that with the above assumptions P 3 = 0 and, thus, P 0 ¼ 1 À 2P q ex . Thus, we can write the excited state population for the qubit as We show in Fig. 4 the above analytic excited qubit state occupation P q ex and the corresponding effective temperature T eff . For our example parameters (see Table 1) and for the intrinsic qubit temperature T q = 100 mK, the data show a saturation of the qubit excitation and the effective qubit temperature to values P q ex % 4 10 À6 and T eff ≈ 36 mK, respectively, for T ≲ 30 mK. In our simulations, we use the value T = 10 mK. In the remaining calculations, we have neglected the intrinsic qubit dissipation in order to keep our results independent on the temperature of the intrinsic environment of the qubit.

Numerical results
Let us compare the analytic model above with the numerical solution of the master equation (18). To this end, we solve the master equation by truncating to the subspace spanned by the five lowest adiabatic energy eigenstates which are obtained by diagonalizing the instantaneous HamiltonianĤ S ðtÞ, defined in Eq.
(2). In particular, the transition rates in Eq. (8) are calculated with the numerically obtained instantaneous eigenfrequencies and eigenstates. We have confirmed that the relative errors caused by the truncation in the ground-state occupation are of the order of 10 −7 or smaller. In Fig. 5, we present the dynamics of the occupations P i (t) for the four lowest-energy states. We study the decay of a single excitation with three different initial occupation probabilities. Our choice for the off-state (see Fig. 2) guarantees that the three degrees of freedom are initially weakly coupled and the states |Ψ 1 (0)〉, |Ψ 2 (0)〉, and |Ψ 3 (0)〉 can be well approximated by the first excited states of the qubit, the right resonator, and the left resonator, respectively. Protocol for qubit initialization J Tuorila et al.
In the numerical simulations, we choose the wait times at the operation points equal to those presented in the analytic model. If the excitation belongs initially to the qubit, the numerical results are in very good agreement with the analytic occupations obtained in Eq. (22). During the sweep to the qubit resonance, part of the initial qubit occupation is transferred non-adiabatically to the left resonator. However, the final numerical fidelity is close to the analytic estimate since at the second operation point the qubit and the left resonator become hybridized and the relaxation rates for the resulting two states are equal. This is clearly visible in Fig. 6 where we show the dependence of the protocol error α(τ) on the total duration τ. Similarly, if the right resonator has a nonvanishing initial occupation, part of it is transferred to the leftresonator during the first sweep. Again, this does not lead to major deviations from the analytic fidelity due to the equal relaxation rates at the first operation point.
We note that one has to be careful when tuning to the first operation point. If the operation point is above the resonance, i.e., ω L (t 1 ) > ω + , any occupation remaining in the state |Ψ 2 (t)〉 at t = t 2 will be transferred to the state |Ψ 3 (t)〉 in a Landau-Zener-type process when the system is subsequently swept across the avoided crossing. This leads to a decrease in the protocol fidelity, since at the second operation point the relaxation rate Γ 30 ≈ 10 −2 Γ 0 (see Fig. 2).
If the desired fidelity P 0 (τ) = 1 − α(τ) is close to unity, the wait times have to be long and, as a consequence, the choice for the initial location of the excitation does not have a large influence on the protocol time. The analytic estimate (24) serves as an upper bound for the total wait time τ, and we approach the upper bound in the case of the full qubit excitation. Fortunately, the numerically obtained fidelity is always higher than that given by the analytic upper bound for the wait time, since there always exist some residual relaxation to ground state from the nonresonant states. If the excitation is partly or completely in the right-resonator, the desired fidelity is reached faster than in the analytic prediction, as depicted in Fig. 6. This occurs because the part of the initial rightresonator occupation that has not decayed at the first operation point can continue the decay at the second operation point due to the hybridization with the left resonator, which causes occupation transfer during the first sweep. We also observe that for our choice of parameters the relatively strong coupling between the qubit and the right resonator, g Rq /δ Rq ≈ 0.14, causes oscillations in the protocol fidelity as a function of the protocol time τ.

DISCUSSION
In summary, we have proposed and modeled a qubit initialization protocol where the coupling between a superconducting qubit and an engineered environment can be externally controlled with a tunable resonator. Using experimentally realizable parameters, we have solved the time-dependent Markovian master equation Fig. 4 a Analytic excited state occupation P q ex ; and b the corresponding effective qubit temperature T eff after the initialization protocol as functions of the resistor temperature T, cf. Eq. (29). In b, we also plot the temperatures for the intrinsic qubit (red dashed line) the resistive bath (black dashed line). We use the intrinsic qubit temperature T q = 100 mK and the intrinsic relaxation rate κ I q ¼ 10 4 s −1 . Other parameters are shown in Table 1 Fig. 5 Dynamics of the occupation probabilities of the low-energy eigenstates. We show results for the initial probabilities P 0 (0) = P 3 (0) = 0, and a P 1 (0) = 0 and P 2 (0) = 1; b P 1 (0) = P 2 (0) = 0.5; c P 1 (0) = 1 and P 2 (0) = 0. Analytic occupation probabilities are shown with dashed lines, and obtained from Eq. (22) by assuming sudden sweeps. The used parameters are shown in Table 1. In the numerical solutions, the sweep time t s = 1 ns and α(τ) = 10 −5 for the protocol and shown that the tunable resonator can be used for fast and precise reset of a qubit excitation. We have also demonstrated that fast changes of the bare angular frequency of the tunable resonator do not reduce the final ground-state fidelity of the protocol. As a result, the duration of the protocol for a given fidelity can be estimated in terms of the decay rate of the tunable resonator, Γ 0 , with a simple analytic model. We also found that the present experimental state-of-the-art decay time 17 for groundstate qubit initialization may be decreased almost by an order of magnitude. Moreover, we observed that at dilution refridgerator temperatures (T≲30 mK) the effective qubit temperature can be reduced to one third of its intrinsic temperature of 100 mK, resulting in an excited qubit state occupation of roughly 10 −6 .
In the off-mode of the protocol, the coupling to the engineered environment has to be weak enough such that the decoherence during quantum computation is determined mainly by the intrinsic environment of the qubit. This sets a limitation for the protocol speed since the decay rate Γ 0 has to be lower than the intrinsic dissipation rates in the system. Furthermore, the assumption of weak coupling with the environment restricts the possible values of decay rates. In future work, the evolution of the reduced density operator should be solved in the regime of strong environmental coupling, which may lead to further improvements in the protocol speed. Also, the possibilities of combining our protocol with driven reset method 17 should be investigated for further improvements on the protocol duration and speed.

Jacobian diagonalization
We perturbatively solve the eigenproblem of the instantaneous Hamiltonian (2). In the RWA, the occupation numberN â † Lâ L þâ † Râ R þσ þσÀ is a conserved quantity becauseĤ S ;N Â Ã ¼ 0. This means thatĤ S ðtÞ andN have joint eigenstates. Eigenvalues of the occupation number are N = n L + n R + n q = 0, 1, 2, …, where n L and n R are the occupation numbers of the left and right resonators, respectively, and assume values n L , n R = 0, 1, 2, …. The number n q denotes the qubit occupation and assumes values 0 and 1. Each occupation number N belongs to (2N + 1) degenerate eigenstates of the form |n L , n R , n q 〉. These are also the eigenstates of the unperturbed HamiltonianĤ 0 ðtÞ. We note that if the counter-rotating terms neglected in the RWA are included, the occupation number is no longer conserved and the above arguments do not hold. In the numerical simulations we employ the non-RWA Hamiltonian.
We are interested in the transitions between the low-energy eigenstates. The N = 1 manifold consists of three states: {|0, 0, e〉, |0, 1, g〉, |1, 0, g〉} where, for clarity, we use the symbols g ↔ 0 and e ↔ 1 for the qubit occupation number n q in the ground and excited state, respectively. As a consequence, the eigenstates of the RWA Hamiltonian with occupation number N = 1 are linear combinations of these three basis states. Therefore, it is enough to find eigenstates for the truncated 3 × 3 Hamiltonian matrix We rely on approximative methods for an intuitive analytic solution for an arbitrary value of ω L (t). The challenge is to find the solution for nearly resonant states, for which non-degenerate perturbation theory fails. The coupling between the left-resonator and the qubit is of second order in the coupling frequencies g 0 LR and g Rq , which causes slow convergence in the conventional nearly degenerate and Brillouin-Wigner perturbation theories. Instead, we first diagonalize the Hamiltonian H S (t) in the N = 1 manifold of the static subspace formed by the right resonator and the qubit, resulting in the Hamiltonian in Eq. (9) which was written as where ω ± ¼ ω av Rq ± Ω Rq . We also recall the short-hand notations q . We then study the corrections caused by the coupling to the left resonator employing the Jacobian diagonalization 42,43 for the eigenproblem. We are especially interested in the transition rates when the left resonator is in resonance either with the right resonator or the qubit. In Eq. (31), HamiltonianĤ S is written in the basis spanned by {|−〉, |+〉, |1, 0, g〉}, where At this point, the left resonator is coupled to the hybridized qubit and right resonator states | ± 〉, which in the dispersive regime g Rq ( δ Rq À Á have the Jaynes-Cummings eigenfrequencies where we have anticipated our choice of parameters and assumed that δ Rq > 0. If the left resonator is far detuned, the perturbation parameters λ L ± ðtÞ ( 1 and hence the above basis states accurately approximate the true eigenstates. In the vicinity of resonances with λ Lþ ðtÞ ) 1 or λ LÀ ðtÞ ) 1, one has to take effects arising from the coupling terms into account by making subsequent diagonalizations in the resonant subspaces. However, if G L ± ( ω þ À ω À j jone can neglect the contributions caused by the off-resonant coupling term G ∓ . When λ L+ (t) ≥ λ L− (t), we make a diagonalizing Jacobian transformation for the Hamiltonian (31) in the subspace spanned by |+〉 and |1, 0, g〉 (see Fig. 7). We obtain where ω "=# ðtÞ ¼ ω av Lþ ðtÞ ± Ω Lþ ðtÞ and G "=#À ðtÞ ¼ G LÀ  Fig. 6 Qubit initialization error as a function of the protocol duration. We present a comparison between different initial states and sweep times. For P 2 ð0Þ ¼ 1, the initialization error is a strongly oscillating function of the protocol duration and its values are located within the shaded region. The physical parameters used are shown in Table 1. The decay time τ 10 = 300 ns are the states corresponding to the approximative eigenvalues hω "=# ðtÞ. Thus in Eq. (34), the HamiltonianĤ S ðtÞ is represented in the basis À j i; # ðtÞ j i; " ðtÞ j i f g . For all values of the control parameter ω L (t) of our protocol, the perturbation parameter λ "À ðtÞ ( 1 and, thus, further Jacobian transformations cause only minor corrections to the third excited eigenfrequency ω 3 ðtÞ % ω " ðtÞ and eigenstate Ψ 3 ðtÞ j i% " ðtÞ j i. In the vicinity of the qubit resonance, δ L− (t) ≈ 0, one has λ L+ (t) < λ L− (t). Consequently, we diagonalize the Hamiltonian (31) in the subspace spanned by the states |−〉 and |1, 0, g〉 and obtain where have neglected the small couplings between the new approximate eigenstates and defined the eigenfrequencies and the corresponding eigenstates Using Eq. (8), we obtain the transition rates (10)- (12). For the eigenfrequency (38) and eigenstate (41) of the second excited state, we use the top signs when λ L+ (t) ≥ λ L− (t) and bottom signs when λ L+ (t) < λ L− (t). Note also that the eigenstates (37) ( 1 for λ L+ (t) < λ L− (t). The Jacobian iteration could be continued further but by comparing these results with the numerical solution of the eigenproblem, we observe that the couplings between the states above are already so small that the Hamiltonian is accurately diagonalized.
Ground-state fidelity during sudden sweeps The ground-state sweep fidelity P S ¼ Ψ 0 ðτÞjΨ 0 ðt 4 Þ h i j j 2 can be calculated in the sudden approximation by finding the corrections to the RWA ground state Ψ RWA 0 ðtÞ ¼ 0; 0; g j i caused by the counter-rotating termsĤ 2 Since these change the occupation number by two, we should expand our low-energy basis to cover the eigenstates of the occupation numbers N = 0, 1, 2, 3, which gives altogether 16 basis states. However, the perturbation divides the eigenspace to two uncoupled sets, one of which is formed by the evenoccupation states and the other one by the odd-occupation states. In conclusion, we do not expect the ground-state sweep fidelity to depend on the matrix elements between the ground state and any eigenstate with odd parity (N = 1, 3, 5, …). Thus, the counter-rotating terms break the symmetry in the RWA Hamiltonian that leads to the conservation of the occupation number, and replace it with a weaker requirement of parity conservation.
We analytically calculate the correction arising from the counter-rotating termĤ 2 ðtÞ only for the ground state, which is created in our low-energy subspace by the off-resonant coupling between |0, 0, g〉 and the N = 2 eigenstates. Since we have g LR ðtÞ= ω L ðtÞ þ ω R ½ (1 and g Rq = ω q þ ω R À Á ( 1 during the whole protocol, we can treat the effects arising fromĤ 2 ðtÞ in the second-order non-degenerate perturbation theory. Since H 2 ðtÞ 0; 0; g j i¼ g LR ðtÞ 1; 1; g j iÀig Rq 0; 1; e j i , we obtain ω 0 ðtÞ ¼ À g 2 LR ðtÞ ω L ðtÞ þ ω R À g 2 Rq ω q þ ω R ; ð43Þ Ψ 0 ðtÞ j i¼ AðtÞ 0; 0; g j iÀ g LR ðtÞ ω L ðtÞ þ ω R 1; 1; g j iþ i g Rq ω q þ ω R 0; 1; e j i ; ð44Þ where the zero of energy is set by the RWA ground-state energy and in the ground state we have neglected the second-order terms, which are outside our low-energy subspace. The deviation from zero energy is caused by the Bloch-Siegert-type 60 non-resonant corrections to the RWA result. The normalization of the state is given by where the latter equality is written to second order in the small parameters g LR (t)/[ω L (t) + ω R ] and g Rq /(ω q + ω R ). Thus, the ground state can be approximately expressed as  Fig. 7 Energy levels in the Jacobian diagonalization scheme for the four lowest eigenstates. We show the Jaynes-Cummings diagonalization for the states |0, 1, g〉 and |0, 0, e〉, and the successive first diagonalizing transformation for the states |1, 0, g〉 and |+〉}, which is relevant for the case λ L+ ≥ λ L− (see text). The resulting approximative eigenfrequencies are ω 1 = ω − , ω 2 = ω ↓ , and ω 3 = ω ↑ . The angular frequencies ω ± , ω "=# , ω av ij , and Ω ij are defined in the text Protocol for qubit initialization J Tuorila et al.