High-fidelity dissipative engineering using parametric interactions

Established methods for dissipative state preparation typically rely on resolving resonances, limiting the target state fidelity due to competition between the stabilization mechanism and uncontrolled dissipation. We propose a protocol devoid of such constraints, using parametric couplings to engineer dissipation for preparation of any maximally entangled two-qubit state. Our scheme allows high-fidelity entanglement generation with short convergence time, continuous control of the target state in the stabilized manifold, and is realizable with state-of-the-art superconducting qubit technology.

Quantum state preparation and preservation are cornerstones of any quantum information platform. Standard methods of state preparation involve a set of unitary operations (or gates) on individual and multiple qubits to achieve a desired entangled state of the system. Such methods typically require multiple tightly synchronized pulse sequences and ancilla qubits, and complex algorithms to avoid unwanted interactions and contain the quantum information within the desired subspace. The prepared state also remains sensitive to environmental decoherence, which reduces the pure quantum state into a classical mixture and ultimately limits the power to harness quantum effects. Recent years have seen an emergence of an alternative approach that embraces the environment instead of hiding it [1,2]. The basic idea relies on engineering suitable interactions with the environment, that steer the reduced system towards a desired target state [3,4]. Besides engineering states inherently robust to dissipation, such dissipative state preparation precludes the need for any active control of the system and is relatively immune to initialization errors. Also, certain quantum operations become feasible only through dissipation; for instance, non-local interactions are essential for the stabilization of multi-partite entangled states such as the 3-qubit GHZ state [5,6]. This has led to interest in dissipative preparation techniques being adopted in a gamut of quantum information platforms such as cavity QED [7,8], trapped ions [9][10][11], superconducting qubits [12,13], Rydberg atoms [14,15], atomic ensembles [16], and NV centers [17].
In this Letter, we present a novel paradigm for engineering dissipation using parametric driving. We demonstrate that this platform enables high-fidelity entanglement generation and control in a circuit-QED setup [18]. In contrast to usual dissipation-engineering schemes, which rely on resonant driving, our scheme exhibits no tradeoff between fidelity and speed of stabilization protocol. It also leverages several technical and operational advantages offered by parametric driving, such as easy implementation relying on continuous-wave (CW) drives alone, strong coupling even in the non-resonant regimes unlike the strong-dispersive circuit-QED [19], and in-situ tunability of parametrically-mediated interactions. The framework of parametrically-engineered dissipation presented here constitutes a significant addition to the parametric quantum toolbox, which is being wielded in burgeoning applications, such as entangling gates [20][21][22], dynamical correction of qubit errors [23], nonreciprocal scattering [24] and synthetic magnetic fields [25], holonomic gates for continuous-variable quantum information [26] and even quantum annealing [27].
Parametrically-engineered dissipation: Our proposed scheme for quantum state engineering consists of a system of two qubits and a single mode cavity (resonator), described by a generic Hamiltonian of the form where X j = σ j + σ † j ; X r = (a + a † ). Each of the pairwise couplings, in addition to the usual static dispersive coupling, are parametrically modulated via a continuous wave pump asg jk (t) =ḡ jk + g + jk cos(ω . Different coupling terms can be activated in the circuit Hamiltonian by choosing the modulation frequencies as either the sum ω + jk = ω j + ω k or difference ω − jk = ω j − ω k of the resonance frequencies of the qubits and resonator. For instance, for preparing even-parity states, the qubit-qubit interactions are modulated at frequencies ω while qubit-resonator interactions are modulated at the four frequencies ω − 1r = ω 1 −ω r −χ 1 ±χ 2 and ω + 2r = ω 2 + ω r − χ 2 ± χ 1 . Here χ j =ḡ 2 jr /(ω r − ω j ) represent the dispersive shifts induced on the qubit states due to respective static qubit-resonator coupling strength g jr . In a rotating frame defined with respect to the dispersive qubit-resonator Hamiltonian, this leads to an effective Hamiltonian of the form [28], in the subspace truncated at the lowest two levels of the resonator. Similarly, for preparing odd-parity states,

leads to an interaction of the form
Henceforth, for simplicity of presentation, we shall assume g ± 12 = g and g ± jr = g r . It is worthwhile to highlight here the features of the Hamiltonians presented in Eqs. (1)-(2) that distinguish them from typical state preparation protocols. (i) As described above, the form of each of the coupling terms is uniquely set by the choice of frequency of the pump mediating the coupling. (ii) Additionally, the dispersive qubit-resonator couplings lead to a selective activation of the qubit-qubit couplings contingent on the photon occupation in the resonator. (iii) By virtue of the couplings being parametric, their strengths and phases are tunable through the amplitudes and phases of the pumps. This freedom of coupling parameters, combined with a decay on the resonator, provides a convenient method for realizing different parametrically-engineered collapse operators c p . This allows the stabilizion of distinct target states, i.e. c p |ξ = 0 [1], through the choice of ω jk and φ jk .
Stabilization Mechanism: We now describe how the above functionalities may be exploited to dissipatively stabilize any maximally-entangled two-qubit state. The generic coupling structure in the combined qubitresonator manifold, truncated at the first three resonator levels, is depicted in Figure 1. Note that while the qubit-qubit drives couple the states of same parity, the qubit-resonator drives mix the two parity manifolds. For the purposes of illustration, we now describe stabilization of an even-parity Bell state implemented using the effective interaction Hamiltonian of Eq. (1) in conjunction with the dissipation provided by the resonator. Assume that the two qubits have been stabilized into |Φ − = (|gg − |ee )/ √ 2. Under a bit flip, the parity of the state changes leading to a jump to the odd parity manifold spanned by {|ge , |eg } ⊗ |0 . Each of these odd-parity states is coherently coupled to the even-parity state |Φ + , 1 , either directly using the qubit-resonator drives (for |eg ) or through a combination of qubit-qubit and qubit-resonator drives (for |ge ), while exciting the resonator. The remaining qubit-qubit drive g − 12 pumps the population from |Φ + , 1 to |Φ − , 1 , which then decays down to the target state |Φ − , 0 as the resonator loses the photons at rate κ. Any population that enters the state |Φ − , 0 remains unaffected by the qubit-qubit drive (since it is off-resonant), and the qubit-resonator drive (since |Φ − is a dark state of the engineered collapse operator c p = σ 1 + σ † 2 ). A phase flip on a qubit moves the population to |Φ + ; this state strongly cou- Coupling diagram for two-qubit entanglement stabilization. Thick black arrows denote resonant interactions, while dotted gray arrows denote off-resonant couplings. Wavy arrows depict photon decay from the resonator. The target state (highlighted with the colored box) and the complementary state (phase-flipped but same parity) are denoted with Greek letters |ξ , |ξ respectively. The states in the opposite parity manifold are denoted with |C , |C respectively. ples to the odd-parity state |ge, 1 through the qubitresonator interaction. The excitation of the resonator induces a dispersive shift of the qubit frequency making the g + 12 drive off-resonant for this state; hence the dominant exit channel is the resonator decay which brings back the population to odd-parity manifold with no photons in the resonator. This is then repumped to |Φ − , 0 as described for the case of qubit decay. The generic coupling map in Fig. 1 describes the stabilization mechanism for all four Bell states, with an appropriate relabeling of states and optimal choice of pump phases as stated in Table I.
Steady state properties and optimization: The steadystate dynamics of the full qubits-resonator system can be described using a Liouvillianρ(t) = Lρ(t), where •} being the usual Lindblad superoperators describing Markovian decay dynamics for the two qubits and the resonator. Here κ denotes the bare resonance linewidth of the resonator, while (γ j 1 , γ j φ ) denote the qubit relaxation and dephasing rates respectively. The resonator acts as an engineered bath for the qubits when κ/γ j 1 1. Figure 2(a) shows the results of a full Liouvillian-based numerical optimization [28], presented as a function of g for four different values of g r and κ. A unique feature of our scheme is the simultaneous reduction in steady state error, ε ∞ , and preparation time, τ , scaling as 1/g when κ = 2g r = (3/2)g. To understand this further,  Fig. 1, calculated as a function of the parametric coupling strength g. The black lines show the results for optimal parameters, κ = 2gr = (3/2)g, with each of the colored dots indicating the value of g for which (gr, κ) are optimal. The solid lines are calculated for T1 = 100 µs, T2 = 200 µs ("best case"), while the dashed lines are calculated for T1 = 10 µs, T2 = 10 µs ("worst case"). The qubit decoherence rates affect the target fidelity, but not the convergence time to target state. (b) Analytical estimates of the total steady-state error, denoted by ε∞, and lifetime of the state |C, 0 decaying into |ξ, 0 , denoted by τ . Each plot also shows the result calculated in the limit κ g (κ-dominated) and in the limit κ g (g-dominated), for a fixed ratio gr/g = 3/4 [28].
we perform an analytical calculation of steady state error and stabilization time for |Φ − (starting from |eg ), using a Liouvillian truncated at lowest two resonator levels [28]. The solid-black curves in Fig. 2(b) present the results of this calculation, and verify that our scheme enables a simultaneous reduction of ε ∞ and τ such that ε ∞ /τ = constant, with the value of the constant set by the ratios g r /g and γ 1 /g. This is in contrast to typical stabilization protocols, which operate in either the κ-dominated regime [29] (decreasing ε ∞ , increasing τ ; dashed-black curves) or the g-dominated regime [30] (increasing ε ∞ , decreasing τ ; dotted-black curves), and are constrained to maintain the product ε ∞ τ = constant. Our proposed scheme realizes a parallel scaling between error and time by balancing the dissipation-induced errors and the coupling-induced errors for each value of κ/g. Furthermore, the minima of the solid-black curves in Fig. 2(b) correspond to the optimal ratios along the black line in Fig. 2(a), where the steady state error and the preparation time scale together as 1/g. Numerical studies also show that the metrics obtained above are robust to imperfections. High-fidelities in excess of 99% are achievable within few 100 ns, for as large as 50% deviations from optimal values of g r or κ, or for asymmetries in the qubit-qubit couplings g + 12 = g − 12 . The performance of the scheme is, however, more sensitive to asymmetries in the parametric qubit-resonator couplings g 1r = g 2r [28]. Other sources of error are the counterrotating terms, which lead to leakage out of |ξ, 0 and compromise the parallel scaling of error and stabilization time. Detailed numerical simulations including the dominant leakage channel in the zero-photon manifold show that the effect of such terms can be suppressed and fidelity can be recovered, either by employing sufficiently large dispersive shifts [28] or by the application of optimal control techniques [31,32]. Continuous-wave Coherent Control : In addition to their strength being tunable with pump amplitudes, the phases of the parametric couplings are determined by the phases of the pumps mediating the respective interactions. This provides a unique knob to implement rotations in the two-qubit stabilized subspace. For instance, Table I shows that the phase entering the collapse operator can be tuned continuously with the phases of the pumps mediating the qubit-resonator couplings. Figure  3 demonstrates how this can be exploited for continuous control of maximally-entangled even-parity states by tuning qubit-qubit coupling phases in tandem with qubitresonator coupling phases. Notably, this rotation of the target state within a given parity manifold maintains the total population and the purity in the manifold constant.
Circuit implementation: In Fig. 4(a), we present a circuit-QED implementation of the dissipative stabilization scheme presented here. In our proposed implementation, two superconducting qubits and a microwave LC resonator share a SQUID-based tunable inductance. The microwave resonator is coupled to an external superconducting transmission line via a coupling capacitor to provide the required linewidth κ. The tunable parametric coupling is realized through modulation of the SQUID inductance with an external flux line according to L sq = Φ 0 /(2πI c cos(πΦ/Φ 0 )), where I c is the critical current of the Josephson junctions, Φ is the flux threading the squid loop and Φ 0 is the magnetic flux quantum. This leads to a static coupling rate g jr between the qubit j and the resonator as where L j , ω j denote the inductance and plasma frequency of the j-th transmon, while L r , ω r denote the inductance and frequency of the resonator. We see that the coupling rate is proportional to the participation ratio between the SQUID and qubit/resonator inductance. For small flux modulation amplitudes, the parametric coupling rate resulting from a flux drive Φ(t) = δΦ jr cos(ω pj t + φ pj ) is given by g jr = (∂g jr /∂Φ)δΦ jr . Figure 4(b) shows a simulation of the qubit-qubit and qubit-resonator coupling rates as a function of flux, for a typical set of circuit-QED parameters. The static coupling rates can be easily tuned over more than 100 MHz, for a flux modulation within half-a-period of flux quantum. In the coupler design presented here, L sq L j to protect the qubits against flux noise from the SQUID.
Discussion: In summary, we have a proposed a scheme that implements fast and robust high-fidelity entanglement preparation using both unitary (qubit-qubit) and dissipative (qubit-resonator) parametric interactions. In particular, our scheme allows simultaneous optimization of fidelity and stabilization rate, since the target state is a dark state of the engineered collapse operator and the system Hamiltonian. Thus, there is no drive-dependent loss mechanism out of the target state that competes with the stabilization mechanism. This is in stark contrast to the usual dissipative state preparation protocols, where parasitic couplings out of the target state impose a hierarchy of coupling rates, and low error can be achieved only by driving the system slower than the effective linewidth of the dressed states to maintain resonant driving [33,34]. Such considerations become increasingly crucial as dissipative engineering protocols are extended to distributed quantum systems, where the speed of state transfer needs to beat the decoherence rate due to correlated noise [35]. Parametrically-engineered dissipation also opens possibilities for dissipation-mediated quantum control in the stabilized state space, by exploiting the phase-tunability of the target state.
In addition to the aforementioned unique features, our scheme offers distinct operational advantages over the previously considered schemes for entanglement stabilization in circuit-QED-like platforms. Since there is no direct driving of the bath resonator, the size of excitations in the resonator,n = a † a , remains small. This is a favorable situation from the point of view of avoiding measurement-induced dephasing of the qubits due to photon number fluctuations, which grows as Γ m ∼nκ [36]. Further, parametric qubit-resonator interactions have no number dependence; hence our scheme employs no photon-selectivity for shuttling excitations across the resonator ladder unlike the usual dispersively coupled schemes [29]. Driving number-selective transitions usually places stringent requirements on matching the dispersive shifts, |∆χ| χ 2 /κ √n , which ceases to be a constraint for parametrically-driven qubit-resonator transitions. Also, our scheme has no direct qubit driving and thus no power-dependent drive detunings due to Stark shifts. Finally, it enables convenient experimental design since, under optimal driving conditions, all the coherent rates (g, g r ) and the dominant dissipation rate (κ) scale linearly with each other.
The main text introduces the following Hamiltonian in the lab frame for parametrically-engineered stabilization with the time-dependent couplingsg The g jr terms correspond to static/dispersive couplings and the remaining terms correspond to time-dependent parametric drives.
To construct the effective Hamiltonians introduced in the main text as Eqs. (1) and (2), we first perform a Schrieffer-Wolff transformation on the dispersive part of the Hamiltonian to diagonalize the dispersive coupling between the qubits and resonator. Using the convention Z |g = + |g and Z |e = − |e , the Hamiltonian becomes with dispersive shifts χ j ≡ g 2 jr /∆ jr where ∆ jr ≡ ω r − ω j . We then move into the interaction frame defined by H 0 = (ω 1 /2) + χ 1 a † a Z 1 + (ω 2 /2) + χ 2 a † a Z 2 + ω r a † a via the unitary operator U = U f U 1 U 2 with the commuting factors We transform each of the three parametric coupling terms in Eq. (S3) individually, starting with the qubit-qubit coupling where we find To find U 1 σ 1 U † 1 , we write σ 1 in Dirac notation |g, q 2 , n e, q 2 , n| , where we have used the basis |q 1 , q 2 , n ≡ |q 1 ⊗ |q 2 ⊗ |n , with q 1,2 ∈ {g, e} and n ∈ N. This gives ,e} e 2iχ1nt |g, q 2 , n e, q 2 , n| .

S2. ANALYTICAL CALCULATION AND PARAMETER OPTIMIZATION
In this section we present an analytical estimation of the steady state error, extracted from the steady state density matrix ρ ss obtained from the Liouvillian as Lρ ss = 0. By truncating the resonator Hilbert space to include two lowest resonator levels, we obtain a closed set of equations, which can be further constrained by assuming that mixing of different states (off-diagonal elements of ρ ss ) due to qubit decays can be ignored in the regime of interest i.e. γ (g, g r , κ). Imposing a relation between qubit-qubit and qubit-resonator couplings, such as g r = (3/4)g, and performing an expansion to leading order in the relevant small parameter, we obtain the following expressions for steady state error in the κ-dominated (or strong dissipation) and g-dominated (or weak dissipation) regimes, respectively: where C = (4g 2 /κγ 1 ) is identified as the cooperativity of the system. Minimizing the total error as the sum of two contributions leads to an optimality condition for the ratio κ/g = (3/2) with the minimum error that scales inversely with parametric coupling strength g ε min ∞ ≈ 16.8 The optimal ratio between dissipative and coherent couplings is in excellent agreement with the value obtained using a full numerical optimization with the given choice of coupling ratio.
The calculation of convergence time τ is more complicated since it is derived from the relevant spectral gap which involves inverting a typically large matrix. For instance, even for a rather conservative system comprising two qubits and a resonator truncated at two levels, an exact calculation involves diagonalization of a 64 × 64 Liouvillian. Nonetheless, an analytical estimate and qualitative scaling with parameters can be obtained by identifying a minimal decoupled subspace including the target state, and self-consistently solving for the populations in different states spanning the subspace. The states {|C, 0 , |ξ, 1 , |ξ, 1 }, |ξ, 0 } form such a subspace (Fig. 1 of the main paper), in which we can write the rate of preparation of the target state |ξ, 0 aṡ P |ξ,0 (t) = Γ C,0 P |C,0 + Γξ ,1 P |ξ,1 + Γ ξ,1 P |ξ,1 .
where P |ψ is the population in state |ψ and Γ |ψ is the decay rate of |ψ into the target state |ξ, 0 . Following the same procedure as for the calculation of steady state error, for a fixed ratio of coupling parameters, the expressions for Γ can be expanded in a relevant small parameter to obtain the respective rates in a given regime. Here we report the value of Γ C,0 , calculated for g r = (3/4)g, (S17) The above expressions lead to a net preparation rate, as . (S18) Analytical estimate of preparation rate in Eq. (S18) serves as an upper bound for the exact result obtained from the Liouvillian gap, i.e. Γ eff ≥ Re[∆ L ]. This is because the total convergence rate is determined by a series of processes that shuffle excitations between multiple decoupled subspaces. It is also worth noting that, to leading order, qubit decays do not affect the convergence time of the scheme as also confirmed by numerical optimization.

S3. NUMERICAL OPTIMIZATION
At any instant in time during the operation of the stabilization protocol, the error from the target state is This error can be decomposed into a static and dynamical error as where the dynamical error decays exponentially to zero as t → ∞ with some characteristic lifetime τ leaving only the steady-state error ε ∞ . Both ε ∞ and τ may be accurately calculated by time-domain simulations of the master equation. Alternately, both ε ∞ and τ may be determined directly from the Liouvillian as ∞ = 1 − Tr(ρ ss I res ⊗ |ξ ξ|) and τ = ∆ −1 L = −Re[λ 1 ] −1 respectively, where Lρ ss = 0 and λ 1 is the lowest lying non-zero eigenvalue of the Liouvillian. We have confirmed that these two methods of calculating the performance metrics of our stabilization protocol are in excellent agreement, with relative errors on the order of 10 −6 for both ε ∞ and τ .

A. Robustness to imperfections
The scaling of performance metrics of our stabilization protocol as a function of g is demonstrated in Figure 2 in the main text. Specifically, the figure shows the effect of varying g while all other parameters are held constant. Increasing g improves the fidelity and the stabilization rate, both of which eventually saturate until the qubit-resonator coupling g r and resonator linewidth κ are also increased. For a specified g, both g r and κ have optimal values determined as κ opt = 2g opt r = (3/2)g. Figure S1 shows the robustness of the stabilization protocol to deviations of κ and g r from their respective optimal values. Large variations in either or both parameters have a limited effect on the achievable fidelities of the stabilization protocol. Figure S2 shows the effect of asymmetries between either the two qubit-qubit couplings g + 12 and g − 12 or between the two qubit-resonator couplings g 1r and g 2r . Of the two possible asymmetries, the case where g 1r = g 2r has a much larger negative impact on the achievable fidelities. This is expected as these couplings are used exclusively to engineer the collapse operator acting on the two-qubit subspace. It is worth noting that since parametric coupling strengths can be tuned with pump amplitudes, in practice the asymmetry between such couplings can be mitigated.    S2: (a) Steady-state error for asymmetric qubit-qubit couplings, calculated with g = (2/3)κ opt = (4/3)g opt r = 2π×50 MHz. Even large asymmetries have a limited effect on the achievable fidelities, and can in fact provide slightly better performance when g − 12 > g + 12 . (b) Steady-state error for asymmetric qubit-resonator couplings, for the same parameter values. In contrast to the qubit-qubit couplings, asymmetry in the qubit-resonator couplings has a detrimental effect on the achievable fidelity. This can be remedied by leveraging the tunable strength of the parametric couplings. The gray region corresponds to the parameter regime where counter-rotating terms induce a fidelity vs. speed tradeoff. The downturn of the curve in the region of large Ωχ/g is due to g becoming comparable to γ1. The dashed line is an approximation of the steady-state error derived from the Liouvillian [Eq. (S15)], in the limit γ1/g 1. (Inset) Time tε taken to reach a given fidelity threshold, for the same Ωχ = 2π × 20 MHz. When the threshold is near the maximum achievable fidelity for a given set of parameters, the nature of the scaling changes. When the threshold is sufficiently less than the achievable fidelity, then the threshold crossing will occur in the region where the dynamical error dominates and the error is decreasing exponentially. If however the threshold is close to the maximum fidelity, then the crossing will occur when the static error is comparable to or larger than the dynamic error and thus the absolute error is falling at a decreasing rate. The regions in which this behaviour is observed are the dotted portions of the curves, caused either by g being too large and thus the leakage term dominating, or by g being comparable to γ1.

B. Effect of counter-rotating terms
The Hamiltonian considered in the main text ignores all non-resonant terms. Including the dominant counter-rotating correction to this Hamiltonian introduces a χ-dependent leakage out of the target state of the form H CR ≈ ge iφ ± 12 e i2Ωχt |ξ, 0 ξ, 0| + h.c.
where Ω χ = χ 1 ± χ 2 . The phase on the prefactor and the sign in the definition of Ω χ is determined by the pump frequencies; when preparing an even parity Bell state + is selected, and -for odd.
When considering the resonant Hamiltonian alone, increasing g does not reduce speed or fidelity as shown in Figure 2 of the main text. Instead, both the speed and achievable steady-state fidelity continue to improve, until they saturate due to processes mediated by g r or κ being the bottleneck. The scaling in the presence of χ-dependent leakage is different, since the effect of the leakage term does not saturate with increasing g. Therefore, in contrast to the resonant case where it is only necessary to optimize g r and κ, introduction of the leakage necessitates setting g to some optimal value g opt found numerically for a given Ω χ . Thus the parameter space can be split into two qualitatively different regimes of operation depending on the ratio g opt /g. When g opt /g 1 the leakage term has little effect and the scaling behavior of the resonant Hamiltonian is recovered, where ε ∞ and τ both scale as 1/g assuming g r and κ scale optimally with g. When g opt /g 1, the achievable steady state fidelity from the resonant Hamiltonian saturates and the scaling of ε ∞ is dominated by the increasingly strong leakage term. In this regime, the ratio ε/τ is no longer constant. These two regimes are demarcated in Figure S3.