Superconducting quantum circuits (SCs) and circuit quantum electrodynamics (cQED) (Blais et al., 2004; Chiorescu et al., 2004; Wallraff et al., 2004; Devoret and Martinis, 2005; Wendin and Shumeiko, 2005; You and Nori, 2006; Clarke and Wilhelm, 2008; Schoelkopf and Girvin, 2008; Devoret and Schoelkopf, 2013; Kockum and Nori, 2019; Krantz et al., 2019; Kjaergaard et al., 2020; Martinis et al., 2020; Blais et al., 2021) have become promising quantum platforms for quantum information processing (Blais et al., 2007; Blais et al., 2020) to implement quantum algorithms (Montanaro, 2016). This method has also been suitable for quantum simulation (Buluta and Nori, 2009; Houck et al., 2012; Georgescu et al., 2014) and, lately, strong evidence has been provided of computational advantages over its classical counterpart (Arute et al., 2019; Wu et al., 2021).

The tailor-made feature offered by the quantum platform permits us to engineer integrated electronic devices that mimic the fundamental elements of cavity QED. For example, artificial atoms correspond to Josephson junction-based circuits behaving as non-linear oscillators whose energy spectrum exhibits noticeable large anharmonicity (Bouchiat et al., 1998; Mooij et al., 1999; Nakamura et al., 1999; Orlando et al., 1999; Martinis et al., 2002; Koch et al., 2007; Schreier et al., 2008; Manucharyan et al., 2009; Barends et al., 2013; Nguyen et al., 2019; Yan et al., 2020). At the same time, LC circuits and transmission line resonators correspond to the quantized field mode (Itoh, 1974; Göppl et al., 2008). The coupling between these subsystems with additional circuit elements as the superconducting quantum interference device (SQUID) allows access to tunable transition energy and switchable coupling strength (Srinivasan et al., 2011). Furthermore, the increasing understanding of the fabrication techniques leads to engineered quantum systems with longer coherence times that protect them from the action of the environment (Wendin and Shumeiko, 2005).

With these elements, we aimed to design a set of logic operations between these artificial atoms that correspond to one- and two-qubit quantum gates, respectively. In circuit quantum electrodynamics, a plethora of two-qubit quantum gates has been proposed, which includes the iSWAP (Bialczak et al., 2010; Dewes et al., 2012; Li et al., 2020), the controlled-not (CNOT) (Paraoanu, 2006; Rigetti and Devoret, 2010; Chow et al., 2011), the controlled-phase (CP) (Strauch et al., 2003; Yang and Han, 2005; DiCarlo et al., 2009; Yamamoto et al., 2010; Xiong et al., 2022), the cross-resonance gate (CR) (Paraoanu, 2006; Rigetti and Devoret, 2010; Chow et al., 2011), and the resonator-induced CP gate (Chow et al., 2013; Paik et al., 2016), among others. The physics underlying these implementations relies upon the tunability of control parameters such that the dynamics on the two-qubit subspace at some sweet spot represent the desired gate (Li et al., 2020). Likewise, it is also possible to engineer quantum gates regarding fixed-frequency qubits that are driven by external microwave pulses (Paraoanu, 2006; Rigetti and Devoret, 2010; Chow et al., 2011; Chow et al., 2013; Paik et al., 2016). The most common situation is where the effective two-qubit interaction appears in a dispersive approximation between the qubits and the mediator (resonator or coupler) (Hua et al., 2014; Srinivasa et al., 2016), leading to a slow gating time that is a detriment to gate performance.

An alternative way to improve gate performance would be to achieve larger coupling strength values or drive the resonator with a high power signal that yields gating time at a sub-nanosecond time scale (García-Ripoll et al., 2003; Campbell et al., 2010; Chen, 2011; Romero et al., 2012). Nevertheless, these extreme considerations affecting the coupling and drive intensity may provoke leakage outside the computational basis and introduce unwanted transitions that will reduce gate performance even more. A way to circumvent this problem relies upon the use of quantum control techniques that accelerate the gating time while maintaining bounded coupling strength between the subsystem and the intensity drive. Within the broad family of quantum control techniques, we assess the shortcuts to adiabaticity (STA) (Chen et al., 2010; Guéry-Odelin et al., 2019), which provide optimal and robust solutions for accelerating slow adiabatic processes while overcoming systematic errors induced by fabrication errors or the unavoidable interaction with the environment (Guéry-Odelin et al., 2019,and references therein). In this regard, STA protocols have been generalized to open quantum system dynamics (Dann et al., 2019; Alipour et al., 2020), and the dissipative qubit readouts have been improved by adding counter-diabatic terms (Yin et al., 2022) or through inverse engineering of the longitudinal coupling (Cárdenas-López and Chen, 2022) in cQED.

Motivated by these advancements, we propose a method to engineer the longitudinal coupling between two transmon qubits coupled to an LC resonator to accelerate the performance of a controlled-phase gate. By using the reverse-engineering approach, we eliminated the unwanted terms appearing on the time-evolution operator governing the system dynamics such that, at a gating time of approximately 2 nanoseconds, the operator represents the required gate. Moreover, we analyzed the performance of the quantum gate under the typical noise source appearing in superconducting circuits and observed no detrimental effects on the final state within the short gating time. Herein, we discuss an experimental implementation with consistent state-of-the-art cQED architecture.

Let us consider a quantized single field mode of frequency ω _r coupled longitudinally to a pair of two-level systems where each of them has transition frequency ω _q,ℓ , and we assume that they are coupled to the resonator through the time-dependent coupling strength λ _ℓ (t). This situation is described by the following Hamiltonian (ℏ ≡ 1): H t = ω r a ̂ † a ̂ + ∑ ℓ = 1,2 ω q , ℓ 2 σ ̂ ℓ z + λ ℓ t σ ̂ ℓ z a ̂ † + a ̂ , (1) where a ̂ ( a ̂ † ) is the annihilation (creation) operator for the bosonic mode and σ ̂ z corresponds to the z-component Pauli matrix describing the two-level system. The longitudinal coupling has been proposed and implemented in several quantum platforms, such as trapped ions (Milburn et al., 2000; Sørensen and Mølmer, 2000) and superconducting quantum circuits (Romero et al., 2012; Didier et al., 2015; Richer and DiVincenzo, 2016). The importance of this interaction stems from the field displacement conditional to the qubit state, with the applications ranging from qubit measurement (Cárdenas-López and Chen, 2022) to the generation of multi-partite quantum states. Here, by employing the technique of STA, we inversely engineered the coupling strength λ _ℓ (t) to reduce the gating time of a controlled-phase gate (Romero et al., 2012). To do so, we began with H in the interaction picture H I t = ∑ ℓ = 1,2 σ ̂ ℓ z Λ ℓ t a ̂ † + Λ ℓ ∗ t a ̂ , (2) where, for simplicity, we renamed the variables as follows: Λ ℓ ( t ) = λ ℓ ( t ) e i ω r t . We proceeded by writing the time evolution operator U I ( t ) = T ̂ exp ( − i ∫ H I ( s ) d s ) , where T ̂ is the time-ordered operator. Using the Baker–Campbell–Hausdorff formula, the time-evolution operator was divided into U I t = ∏ ℓ exp i A ℓ t σ ̂ ℓ z a ̂ ∏ ℓ exp i A ℓ ∗ t σ ̂ ℓ z a ̂ † ∏ ℓ , ℓ ′ exp i B ℓ , ℓ ′ t σ ̂ ℓ z σ ̂ ℓ ′ z , (3) where A ℓ ( t ) and B ℓ ( t ) are functions that can be obtained through the Schrödinger equation i ∂ t [ U I ( t ) ] = H I ( t ) U I ( t ) (see Supplementary Material for the detailed derivation). After some lengthy, but straightforward, calculations, we found that these functions needed to satisfy the following differential equations: A ̇ ℓ t = − λ ℓ t e i ω r t , B ̇ ℓ , ℓ ′ t = − i A ̇ ℓ ∗ t A ℓ ′ t . (4) For time-independent coefficients λ _ℓ , these equations reduced to A ℓ ( t ) = i λ ℓ ( e − i ω r t − 1 ) , and B ℓ , ℓ ′ ( t ) = λ ℓ λ ℓ ′ [ − i ( e i ω t − 1 ) − t ] / ω r , which vanished at t = 2π/ω _r . Nevertheless, in the context of the STA technique, it was possible to engineer λ _ℓ (t) such that the gate reduced to t _g = 1.89 (ns), while ensuring that the following final boundary conditions, A ℓ t g = 0 , B ℓ , ℓ ′ t g = − Θ ℓ , ℓ ′ , (5) were fulfilled at the gating time t _g . This gives the time-evolution operator as follows: U I t g = ∏ ℓ , ℓ ′ exp i Θ ℓ , ℓ ′ σ ̂ ℓ z σ ̂ ℓ ′ z . (6) Coming back to the Schrödinger picture, we obtained U t g = ∏ ℓ exp − i η ℓ σ ̂ ℓ z ∏ ℓ , ℓ ′ exp i Θ ℓ , ℓ ′ σ ̂ ℓ z σ ̂ ℓ ′ z , (7) with η _ℓ = ω _q,ℓ t _g /2. Thus, for a specific gating time t _g , and by engineering the coupling strength λ _ℓ (t), it was possible to express the time evolution operator U (t _g ) as a controlled phase gate.

In this section, we aimed to design a modulation for the coupling strength λ _ℓ (t) such that, at a desirable gating time t _f , the time-evolution operator in Eq. 7 represents a controlled-phase gate acting on the qubit subspace. In other words, we intend that U (t _g ) ≡ CZ ⊕I. To do so, we had to find the unitary transformation V ̂ ( t ) that expressed the Hamiltonian in its adiabatic frame (Guéry-Odelin et al., 2019). In this frame, the Hamiltonian is block diagonal. Conversely, dynamical evolution is always along the instantaneous eigenstates of the system without transitions. More precisely, V ̂ ( t ) must satisfy the following relation: H e ff t = V ̂ t H t − i ∂ t V ̂ † t = H q ⊕ H r , (8) where ∂ _t refers to the time-derivative acting over operators. In general, it is possible to express V ̂ ≡ ∏ ℓ ⁡ exp ( S ̂ ℓ ( t ) ) as a product of unitaries, where S ̂ ℓ is its ℓth generator. Though such generators can be calculated perturbatively through the Schrieffer–Wolff transformation (Schrieffer and Wolff, 1966), the generators S ̂ ℓ of the transformation for the longitudinal interaction (Chen et al., 2018) are exactly solvable (Čadež et al., 2014) and given by V ̂ t = ∏ ℓ e i θ ℓ t ⁡ exp − i g ̇ c , ℓ t ω r 2 σ ̂ ℓ z a ̂ † + a ̂ exp − g c , ℓ t ω r σ ̂ ℓ z a ̂ † − a ̂ , (9) where g _c,ℓ (t) is an auxiliary (classical) variable that eliminates the coupling strength in the Hamiltonian H ( t ) ; moreover, θ ℓ ( t ) = − ∫ 0 t L g , ℓ ( t ′ ) d t ′ is a phase relating the coupling strength λ _ℓ (t) with the classical variable g _c,ℓ (t) through the Lagrangian, L g , ℓ t = g ̇ c , ℓ 2 t ω r 3 − g c , ℓ 2 t ω r + 2 g c , ℓ t λ ℓ t ω r , (10) that has the following Euler–Lagrange equation: g ̈ c , ℓ t + ω r 2 g c , ℓ t − λ ℓ t = 0 , (11) from which the control problem is reduced to solving a set of N differential equations (N is the number of qubits). The initial conditions for these equation sets require that, at t = 0 and t = t _g , both adiabatic and lab frame must coincide (Chen et al., 2018; Guéry-Odelin et al., 2019; Cárdenas-López and Chen, 2022). Therefore, the initial conditions are given by g c , ℓ 0 = g c , ℓ t g = g ̇ c , ℓ 0 = g ̇ c , ℓ t g = 0 , ∀ ℓ ∈ 1,2 . (12) To assure smooth coupling strength modulation, we added the extra constraints g ̈ c , ℓ ( 0 ) = g ̈ c , ℓ ( t g ) = 0 . Furthermore, to guarantee that the time-evolution operator satisfied U _I (t _g ) = CZ, we demanded that the coupling strength λ _ℓ (t) met the conditions pointed out in Eq. 5 for any value of the angle Θ _ℓ,ℓ′. We solved this set of differential equations by proposing two different solutions, with the first of them corresponding to a trigonometric protocol pulse (Theis et al., 2018) of the form g c , ℓ Tri t = ∑ n = 0 M a n , ℓ ⁡ sin π n t t g k , (13) where M is the number coefficient required to fulfill all the initial conditions, and g _c,ℓ (t) and its (k − 1)th derivatives vanish at t = 0 and t = t _g for all values of n. Selecting k = 3 resulted in g _c,ℓ (t), g ̇ c , ℓ ( t ) , and g ̈ c , ℓ ( t ) vanishing at those times. This feature of the modulation permitted reduction of the number of coefficients to four, which corresponds to the real and imaginary part of A ℓ ( t ) and B ℓ , ℓ ′ ( t ) , respectively. Furthermore, we wanted to compare the trigonometric solution with a polynomial form of STA modulation for completeness (Chen et al., 2018): g c , ℓ STA t = ∑ n = 0 M ′ b n , ℓ t n + 3 . (14) The STA modulation satisfied construction of the conditions g c , ℓ ( 0 ) = g ̇ c , ℓ ( 0 ) = g ̈ c , ℓ ( 0 ) . In that case, we required seven coefficients to completely solve the set of differential equations. Engineering the auxiliary functions g c , ℓ ⋆ allowed us to find the coupling strength modulation through the relation λ ℓ ( t ) = g ̈ c , ℓ ⋆ ( t ) / ω r 2 + g c , ℓ ⋆ ( t ) , ⋆ = {Tri,STA}. Figure 1 shows the coupling strength modulation λ _ℓ (t), and the auxiliary functions A ℓ ( t ) and B _ℓ,ℓ′(t) as a function of the normalized time t/t _g for both anzats given in Eqs 13, 14 for different phases Θ _ℓ,ℓ′. From Figure 1, we can appreciate that neither modulation exhibited any kind of discontinuity, and the main difference between them is that the polynomial anzat takes negative values; in both cases, the variation of the final phase generated a family of shortcut coupling strength modulations {λ _ℓ (t)}. Next, we analyzed the system dynamics of the Hamiltonian H ( t ) using the coupling strength modulation λ _ℓ (t) obtained through the trigonometric solutions and STA.

As mentioned previously, generating fast and high-fidelity two-qubit quantum gates is at the heart of quantum computations and their applications in simulations and information processing, among other disciplines. In particular, the controlled-phase gate, i.e., CZ, changes the phase of the quantum state conditionally. In the basis {|gg⟩, |ge⟩, |eg⟩, |ee⟩}, it has the following form: C Z Θ ℓ , ℓ ′ ≡ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 e i Θ ℓ , ℓ ′ . (15) The CZ gate acts as an entangling gate when the initial quantum state is of the form |Ψ⟩ = (α|g⟩ + β|e⟩)(a|g⟩ + b|e⟩). The change phase on the state |ee⟩ does not permit writing the state as a tensor product. Conversely, the entanglement of formation (EoF) in the state is different from zero. For instance, for a quantum state with balanced amplitudes α = β = a = b = 1 / 2 , the CZ gate with Θℓ,ℓ′ = π/2 generates a maximally entangled state. In particular, for a pure two-qubit state, the maximal entanglement is achieved when 2|αaβb (1 − exp iΘℓ,ℓ′)| = 1. Thus, having knowledge about the initial state of the system allowed us to use reverse-engineering to set an adequate angle Θℓ,ℓ′ that maximized the EoF of the system.

More specifically, with our reverse-engineering protocol, the system dynamics reduced to the time-evolution operator depicted in Eq. 7. During the time evolution, the states {|gg⟩, |ge⟩, |eg⟩, |ee⟩} accumulated different phases given by U t g = e i Θ ℓ , ℓ ′ diag e 2 i η , 1,1 , e − 2 i η × diag e 4 i Θ ℓ , ℓ ′ , 1,1 , e 4 i Θ ℓ , ℓ ′ , (16) where the equal qubit frequency, i.e., η _ℓ = η = ω _q t _g /2 was assumed. The first term on the time-evolution operator corresponds to the free terms of each qubit, whereas the second one is the effective ZZ interaction appearing due to factorization of the cavity degree of freedom. To obtain the typical controlled-phase gate CZ = diag [1, 1, 1–1], up to a global phase factor, we needed to solve 2 η + 4 Θ ℓ , ℓ ′ = 2 π n , − 2 η + 4 Θ ℓ , ℓ ′ = 2 n − 1 π , ∀ n ∈ Z . (17) For the shorter generation time, we chose n = 0 and obtained Θ _ℓ,ℓ′ = ±π/8, which resulted in U (t _g ) ≡ CZ = diag[1, 1, 1–1] up to a global time-dependent phase factor. Different values of Θ _ℓ,ℓ′ will accumulate relative phases on the other state component, thereby degrading the performance of the quantum gate.

Figure 2 shows the system dynamics of the Hamiltonian in Eq. 1 for two different trigonometric and STA modulations. We calculated the population of the pair of two-level systems, the resonator and the EoF, of the two-qubit subspace for different angles Θ _ℓ,ℓ′ when the system was initialized in the state |Ψ(0)⟩ = | + ⟩ ⊗| + ⟩ ⊗|0⟩, where | + 〉 = ( | g 〉 + | e 〉 ) / 2 is the eigenstate of x Pauli matrix σ ^x and |0⟩ is the vacuum state of the resonator, respectively. Figure 2 demonstrates no change in the population in the two-qubit subspace, implying no excitation exchange between the qubits and the resonator. On the other hand, in the resonator subspace during the system dynamics, the resonator populated higher excitation states due to the cavity displacement introduced by the interaction term of the Hamiltonian in Eq. 1. Nevertheless, at the end of the system dynamics, the resonator recovered the population in the ground state (Figures 2B–D). Therefore, it makes sense that STA was not adiabatic at all during the intermediate time, but the initial and final states coincided with those that were adiabatic due to the initial and final boundary conditions imposed previously.

Next, we turned to the correlations embedded on the two-qubit subspace through the EoF, which, for this system, have an exact form given by E o f | Ψ 〉 = − x ⁡ log 2 x − 1 − x log 2 1 − x , x = 1 + 1 − C 2 2 , (18) where the quantity C is termed concurrence defined as C = max(0, E ₄ − E ₃ − E ₂ − E ₁), where E _k are the eigenvalues (in increasing order) of the matrix R = ( σ y ⊗ σ y ) ( | Ψ 〉 〈 Ψ | ) * ( σ y ⊗ σ y ) | Ψ 〉 〈 Ψ | (Wootters, 1998). Figures 2C–F show the EoF as a function of the normalized gating time t/t _g for both trigonometric and STA modulations. By changing the angle Θ _ℓ,ℓ′ in both cases, we observed existence of the optimal value of π/8, where the EoF reached its maximal value for the state |Ψ(0)⟩ = | + ⟩ ⊗| + ⟩ ⊗|0⟩. A slight difference between the modulations occurred via the trigonometric-like pulses that exhibited a plateau in the EoF at the end of the dynamics, which did not occur with the STA modulation. The existence of this plateau indicates that only trigonometric-like modulation achieves optimal correlation values at times that are slightly shorter than t _g = 12π/ω _r ≡ 1.89 ns. However, we observed identical overall performance for both modulations. The next step in the characterization of our CZ gate relied on its performance under a loss mechanism, such as energy relaxation in the resonator and the two-level systems and depolarizing noise, respectively. In this scenario, the system dynamics were governed by the following master equation: ρ ̇ t = i H t , ρ t + ∑ ℓ = 1,2 γ ℓ D σ ℓ − ρ t + γ ϕ , ℓ D σ ℓ z ρ t + κ D a ̂ ρ t , (19) where ρ(t) is the density matrix describing the quantum state at time t and H ( t ) is the Hamiltonian in Eq. 1. Moreover, {γ _ℓ , γ _ϕ,ℓ } correspond to the relaxation and dephasing rate of the ℓth two-level system, respectively. Also, κ is the relaxation rate of the resonator. Finally, D [ L ] ρ ( t ) = L ρ ( t ) L † − { L † L , ρ ( t ) } / 2 is the Lindbladian superoperator describing the loss mechanism. We performed the numerical simulation by using experimental physical parameters for the decay rates, yielding γ _ℓ = 0.48 (MHz) (Saira et al., 2014), γ _ϕ,ℓ = 0.15 (MHz) (Saira et al., 2014), and κ = 21 (MHz) (Didier et al., 2015), respectively.

Figure 3 shows the population of the two-qubit and resonator subspace, and the EoF as a function of the dimensionless gating time t/t _g for the density matrix ρ(t) solution of the master equation in Eq. 19. There was no appreciable detrimental effect on the two-qubit subspace, and the population changed slightly and only at the end of the system dynamics without significantly affecting the performance of the protocol. For the resonator population, slight distortions of the population evolution were observed, which were produced mainly by the dissipation on the resonator, since its presence altered the displacement trajectory of the phase space. For the evolution of the EoF, we also observed that, for the optimal angle Θ _ℓ,ℓ′ = π/8, there was no appreciable detrimental effect on the EoF at the end of the dynamics. However, for different angles, the generated correlation drastically decreased. In this situation, we concluded that the accumulated phase on each state component was affected by the depolarizing noise, since the relaxation losses did not abruptly change the population in the two-qubit state. In this sense, the controlled-phase gate was robust under relaxation noise, but fragile under depolarizing noise when the phase was not working at the optimal value Θ _ℓ,ℓ′ = π/8.

On the other hand, to demonstrate the robustness of the controlled-phase gate implementation over different initialization, we calculated the Schrödinger equation and the master equation [Eq. (19)] for the system prepared in the initial state |Ψ(0)⟩ = [cos(θ)|g⟩ + exp (iϕ) sin(θ)|e⟩] ⊗|e⟩ ⊗|0⟩ for different angles and phases in the range θ = (0, 2π) and ϕ = (0, 2π), and compared the result after applying the CZ gate. The comparison was carried out through the infidelity defined as 1 − F , where F ( ρ , σ ) = ( Tr [ ρ σ ρ ] ) 2 is the fidelity between the quantum states ρ and σ. Figure 4 shows the infidelity for those cases, demonstrating that the implementation of the controlled-phase gate using the reverse engineering method achieves good fidelity for a wide initial state preparation.

In this section, we discuss a possible experimental implementation for a single-mode resonator coupled to a pair of transmon circuits, as denoted in Figure 5. The single-mode resonator was an LC circuit with capacitance C _r and inductance L _r . Moreover, we modeled each transmon circuit (Koch et al., 2007) as a capacitor C Σ ℓ parallel-connected to a Josephson junction with tunable energy E J ℓ ⁡ cos ( θ x , ℓ ) ; these subsystems were coupled through a SQUID. In this manner, we write the Lagrangian of the circuit in terms of the flux nodes of each device { ψ J 1 , ψ J 2 , ψ r } related with the voltage drop across their respective branch ψ ℓ = ∫ − ∞ t V ℓ ( x , t ′ ) d t ′ leading to L = ∑ ℓ = 1,2 C Σ ℓ 2 ψ ̇ J ℓ 2 + E J θ x , ℓ cos ψ J ℓ φ 0 + E J S ℓ φ x , ℓ cos ψ J ℓ − ψ r φ 0 + C r ψ ̇ r 2 2 − ψ r 2 2 L r , (20) where φ ₀ = ℏ/(2e) is the quantum magnetic flux and e is the electron charge. The canonical conjugate momenta are given by P = ∂ L c / ∂ [ ψ ̇ ] P J ℓ = C Σ ℓ ψ ̇ J ℓ , P r = C r ψ ̇ r . (21) We calculated the circuit Hamiltonian through the Legendre transformation H = ∑ k P k ψ k − L , arriving at H = ∑ ℓ = 1,2 P ℓ 2 2 C Σ ℓ − E J θ x , ℓ cos ψ J ℓ φ 0 − E J S ℓ φ x , ℓ cos ψ J ℓ − ψ r φ 0 + P r 2 2 C r + ψ r 2 2 L r . (22) The quantization of each degree of freedom was achieved by promoting the momentum (node charge) and the flux functions to quantum operators satisfying canonical commutation relations. For the transmon qubit, the charge is proportional to the number of Cooper-pairs on the superconducting island, i.e., P ℓ → − 2 e n ̂ ℓ and its conjugate variable corresponds to the superconducting phase drop θ ̂ J ℓ = ψ J ℓ / φ 0 satisfying commutation relation [ P ̂ ℓ , e ± i θ ̂ J ℓ ] = e ± i θ ̂ J ℓ . For the LC resonator, the quantized charge and phase operator satisfy [ θ ̂ r , n ̂ r ] = i . In this representation, we write the Hamiltonian as follows: H = ∑ ℓ = 1,2 4 E C J ℓ N ̂ J ℓ 2 − E J θ x , ℓ cos θ ̂ J ℓ − E J S ℓ φ x , ℓ cos θ ̂ J ℓ − θ ̂ r + 4 E C r N ̂ r + E L 2 θ ̂ r 2 . (23) Here, E C J ℓ = e 2 / 2 C Σ ℓ and E C r = e 2 / 2 C r are the charge energy of the transmon and the resonator, respectively. Moreover, E L = φ 0 2 / L r is the inductive energy of the oscillator. We then assumed that the SQUID works in the linear regime (Leib and Hartmann, 2014), meaning that most of the current flows through the transmon. Hence, the resonator phase is well located, allowing expansion of potential energy up to its leading order in θ ̂ r (Leib et al., 2012) H = ∑ ℓ = 1,2 4 E C J ℓ N ̂ J ℓ 2 − E J θ x , ℓ cos θ ̂ J ℓ − E J S ℓ φ x , ℓ cos θ ̂ J ℓ − E J S ℓ φ x , ℓ cos θ ̂ J ℓ θ ̂ r + 4 E C r N ̂ r + E L 2 θ ̂ r 2 . (24) We expressed the oscillator in terms of the harmonic oscillator basis described by the operators N ̂ r = E L 32 E C r 1 / 4 i a ̂ † − a ̂ ≡ N zpf i a ̂ † − a ̂ , (25) θ ̂ r = 2 E C r E L 1 / 4 a ̂ † + a ̂ ≡ θ zpf a ̂ † + a ̂ . which lead to the following Hamiltonian: H = ℏ ω r a ̂ † a ̂ + ∑ ℓ = 1,2 4 E C J ℓ N ̂ J ℓ 2 − E J θ x , ℓ , φ x , ℓ cos θ ̂ J ℓ − λ ̃ ℓ φ x , ℓ cos θ ̂ J ℓ a ̂ † + a ̂ , (26) where E J ( θ x , ℓ , φ x , ℓ ) = E J ( θ x , ℓ ) + E J S ℓ ( φ x , ℓ ) is the total Josephson energy of the transmon circuit, ω r = 8 E C r E L / ℏ is the resonator frequency, and λ ̃ ℓ ( φ x , ℓ ) = θ zpf E J S ℓ ( φ x , ℓ ) is the coupling strength between the ℓth transmon circuit and the resonator. For the next analysis, it was convenient to divide the circuit Hamiltonian into three parts: H T = ∑ ℓ = 1,2 4 E C J ℓ N ̂ J ℓ 2 − E J θ x , ℓ , φ x , ℓ cos θ ̂ J ℓ , (27) H r = ℏ ω r a ̂ † a ̂ , (28) H I = − ∑ ℓ = 1,2 λ ̃ ℓ φ x , ℓ cos θ ̂ J ℓ a ̂ † + a ̂ , (29) corresponding to the transmon, resonator, and interaction terms.

To illustrate this in the two-level approximation of the transmon circuit, the two-level approximation H I turns on a longitudinal oscillator–qubit interaction of the form σ ^z . In the transmon regime ( E C J ℓ ≪ E J ( θ x , ℓ , φ x , ℓ ) ) , it is possible to neglect the charge term on Hamiltonian H T . In this approximation, we have that H T ∝ cos θ ̂ J ℓ . For a single transmon, in the charge basis, we chose the states | g 〉 ≡ | 0 J ℓ 〉 and | e 〉 ≡ | 1 J ℓ 〉 as the computational basis of the transmon, which leads to: H T = − ∑ ℓ = 1,2 E J θ x , ℓ , φ x , ℓ 2 σ ℓ x , (30) H I = ∑ ℓ = 1,2 λ ℓ φ x , ℓ σ ℓ x a ̂ † + a ̂ , (31) where σ ^x = |e⟩⟨g| + |g⟩⟨e| corresponds to the x-component Pauli matrix, and λ ℓ ( φ x , ℓ ) = λ ̃ ℓ ( φ x , ℓ ) ⟨ g | cos ( θ ̂ J ℓ ) | e ⟩ is the dressed coupling strength between the ℓth transmon and the resonator. In the diagonal basis of the transmon, we obtained H T = − ∑ ℓ = 1,2 ℏ ω q , ℓ 2 σ ℓ z , (32) H I = ∑ ℓ = 1,2 λ ℓ φ x , ℓ σ ℓ z a ̂ † + a ̂ , (33) where the qubit frequency is defined as ω q , ℓ = E J ( θ x , ℓ , φ x , ℓ ) / ℏ . Figure 6 shows the coefficients α k = Tr [ σ k ⁡ cos ( θ ̂ J ) ] as a function of both external magnetic fluxes θ _x,ℓ and φ _x,ℓ , respectively. We observed that the only relevant contribution corresponds to the σ ^z operator. Thus, our final Hamiltonian reads: H t = ω r a ̂ † a ̂ + ∑ ℓ = 1,2 ω q , ℓ 2 σ ̂ ℓ z + λ ℓ φ x , ℓ σ ̂ ℓ z a ̂ † + a ̂ . (34)

It is worth mentioning that it is possible to perform single qubit gates on this architecture just by choosing a value of the external magnetic flux φ _x,ℓ such that we switch off the interaction term between the resonator and the transmon circuits. On the other hand, the anharmonicity of the transmon circuit is proportional to its charge energy E _C = 238 MHz, which, for a typical transmon, could be small enough to consider the additional level. In our proposal, however, as the artificial atom is coupled longitudinally to the field mode, we may not expect either population leakage or transitions induced by them. We demonstrated calculation of the protocol to generate the CZ gate by including the third energy level, as depicted in Figure 7. As expected, we did not observe substantial modification to the system dynamics.

In summary, we have proposed a reverse-engineering method to modulate the longitudinal interaction between a pair of two-level systems with a cavity field to implement a controlled-phase gate using a pair of superconducting artificial atoms. We have constrained the modulation of the coupling strength to delete the unwanted terms appearing on the time-evolution operator, resulting in the two-qubit quantum gate that operates at 2 nanoseconds, which is within the current state-of-the-art in superconducting quantum circuits. The shorter generation time allows the gate to be insensitive to the unavoidable effects of the environment. The method can be extended to the fast control of multi-particle and other physical setups, such as trapped ions and quantum dots. The use of STA paves the way to implementation of faster and more resilient quantum gates by providing a key factor for quantum computation, quantum information processing, and quantum simulation.