The Optimization of Flux Trajectories for the Adiabatic Controlled-Z Gate on Split-Tunable Transmons (2024)

In this section, we briefly describe the model that the rest of this work is based on, which is used for both simulation and motivation of our optimization method. In addition, we review how a controlled phase gate naturally arises in the adiabatic control of the system when levels in the qutrit subspace of the transmon are considered. In our setup, two flux-tunable transmons are employed in the standard arrangement: the qubit frequency is tuned by applying magnetic flux through the symmetric split-junction (a dc-SQUID)[Fig. 1(a)]. We use simulation parameters similar to those in Ref. 10, in which Q1(Q2) has frequency $\omega/2\pi=5.889(5.031)$ GHz and anharmonicity $\alpha/2\pi=-324.3(-234.7)$ MHz, with coherence times $T1=25.5(48.8)$ $\mu$s and $T2=13.3(28.4)$ $\mu$s, with coupling strength between qubits $g/2\pi=24.7$ MHz. Though a common choice at present for a coupling setup is cQED⁶, for analytical simplicity we use direct capacitive coupling between the qubits, yet the strategies presented here are generalizeable to a cQED model (such as the generalized Tavis-Cummings Hamiltonian⁸). In this configuration, single-qubit microwave pulses are implemented through the capacitive coupling of each transmon to a drive line controlled by an arbitrary waveform generator[Fig. 1(b)]¹¹. The operational point for these microwave pulses occurs when there is no applied flux, i.e. $\varphi=0$; at this base point, the qubit frequencies of the transmons are furthest detuned from each other, thus minimizing the effect of the coupling as a source of decoherence. Fig. 1(c) shows the explicit effect of the flux detuning on the frequency for various initial states. In the tensor product eigenbases of the individual transmons (including coupling to their own drive lines), the Hamiltonian governing the capacitive interaction is⁶

where $g$ is the coupling strength between the transmons, $\omega^{j}_{1}$ and $\omega^{j}_{2}$ are the respective qubit frequencies of $j$th level of each transmon, $\mathcal{I}_{i}$ is the identity operator on the Hilbert space for the $i$th transmon, and we have used the usual oscillator annihilation and creation operators $a$ and $a^{\dagger}$. Since the coupling $g/2\pi$ is significantly smaller (roughly 200-fold at $\varphi=0$) than the qubit frequencies, the rotating wave approximation is justified, reducing the interaction term of Eq. 1 to

In what follows, by $\omega_{ij}$ we mean the the eigenvalues of $H$ corresponding to the $\ket{i,j}$ state, for which at $\varphi=0$, $\omega_{00}=0,\omega_{01}\approx\omega^{2}_{1},\omega_{10}\approx\omega^{1}_{1}$, and as stated, $\omega_{11}\approx\omega^{1}_{1}+\omega^{2}_{1}$. In addition, when discussing the matrix representation of an operator on a subspace of $\mathcal{H}$, we will work in a basis ordered by eigenvalues, which does not preserve standard qubit ordering when its domain is restricted to $\mathcal{H}_{qq}$. To naturally create a CZ operation, we shall work in a seven-dimensional subspace of $\mathcal{H}$, namely $\mathcal{H}_{Z}\subset\mathcal{H}$, spanned by the basis set $\mathcal{B}_{7}=\{\ket{i,j}|i,j\leq 2\}-\{\ket{0,0},\ket{2,2}\}$. We exclude $\ket{0,0}$ and $\ket{2,2}$ from this consideration since $H_{\mathrm{int}}$ as defined in Eq. 2 when restricted to the full qutrit subspace does not have any matrix elements that generate interaction in the first or ninth row or column (explicitly, $\bra{i,j}\!{H_{\mathrm{int}}}\!\ket{2,2}=\bra{2,2}\!{H_{\mathrm{int}}}\!\ket{i$ holds in the restriction to the qutrit subspace). For both notational simplicity and possible interpretational advantage, we write the Hamiltonian by means of block matrices, utilizing the well-known Pauli matrices and the lesser-known Gell-Mann matrices, whose use in a qutrit-like consideration is somewhat fundamental since, as the Pauli matrices span the Lie algebra $\mathfrak{su}(2)$, the eight Gell-Mann matrices span $\mathfrak{su}(3)$. For a complete list, description, and interpretation of these matrices, see Ref. 12; here we will only use them for convenience of notation (with the symbol $\lambda$). Eq. 2 under this restriction (and the earlier mentioned eigenvalue ordering) becomes

on $\mathcal{B}_{7}$ (for operator-matrix equivalence in some specified basis we use the symbol $\simeq$). This operator representation hints at contributions important to mitigate in the execution of CZ: for example, the $g\sigma_{1}$ term represents a swapping action between $\ket{1,0}$ and $\ket{0,1}$ excitations at the coupling strength. In fact, this term can be used to naturally induce an iSWAP gate¹³, however this should be avoided for a successful and coherent CZ. In addition, the $g\sqrt{2}\lambda_{1}$ and $g\sqrt{2}\lambda_{6}$ terms represent the main issue, i.e. swapping between $\ket{0,2}$ and $\ket{1,1}$, and $\ket{2,0}$ and $\ket{1,1}$ respectively; in fact only the latter will be an issue, since we shall drive the transmon close to a point at which swapping is highly likely. Lastly, we have another swapping interaction with strength $2g$ between $\ket{1,2}$ and $\ket{2,1}$, however ideally these should not matter much if leakage out of $\mathcal{H}_{qq}$ has been well-managed throughout the quantum gate sequence.

Fig. 2 shows the low, medium, and high frequency ranges of the spectral decomposition of $H_{\mathrm{int}}$, in which (a) and (b) show two avoided crossings useful for inducing various state changes. Most importantly, (c) gives an idea of how the CZ is implemented: the black dashed line, the sum of $\omega_{01}$ and $\omega_{10}$, deviates from $\omega_{11}$ as flux is increased, which occurs because of the downward push on $\ket{1,1}$ from $\ket{2,0}$. The optimal flux bias for CZ is $\mathcal{A}_{2}$, the avoided crossing between these two states, or just before it⁹, since being in the $\ket{1,1}$ state and moving towards $\mathcal{A}_{2}$ maximizes both accumulation of the difference in its frequency and that of the previously noted sum and the probability of remaining in $\ket{1,1}$ during the movement. As we have assumed coupling is minimal at zero flux bias, we can approximate the eigenstates of $H$ at $\varphi=0$ to be the individual transmon eigenstates. Parameterizing $H$ by $\varphi$, if $\partial_{\varphi}H\approx 0$, i.e. the trajectory $\varphi(t)$ is adiabatic, then the system stays in its instantaneous eigenstate¹⁴ at $\varphi(0)=0$. The rotating frame of the transmon at which the remainder of the quantum gate sequence is operated in is defined at the zero flux-bias point, so the Hamiltonian does have time-dependent components, but only along its diagonal. Since the Hamiltonian commutes with itself at different times, the Schrodinger equation and the adiabatic approximation dictates that the trajectory-dependent unitary $U_{A}:[0,\infty)^{T}\to\mathcal{B}(\mathcal{H}_{Z})$

(here $\mathcal{B}(\mathcal{H})$ denotes the Banach space of bounded linear operators on a Hilbert space $\mathcal{H}$) is applied to $\mathcal{H}_{Z}$ in the basis $\mathcal{B}_{7}$, where $T=\mathrm{dom}(\varphi)$ is the time interval over which $\varphi$ acts, and

where for some $t$, $\delta_{ij}(\varphi(t))=\omega_{ij}(0)-\omega_{ij}(\varphi(t))$ is the flux-dependent deviation of the eigenfrequency corresponding to $\ket{i,j}$. We conclude this section by noting that, if we choose some flux trajectory $\varphi_{\pi}$ such that $\Theta_{11}(\varphi_{\pi})-\Theta_{01}(\varphi_{\pi})-\Theta_{10}(\varphi_{\pi$, then the CZ operation is realized up to global phase as the product of the restriction of $U_{A}$ and single-qubit gates, namely

where $\mathcal{H}_{A}=\mathcal{H}_{Z}\,\cap\,\mathcal{H}_{qq}$, the equivalence relation $\sim_{\phi}$ is defined for operators $A,B\in\mathcal{B}(\mathcal{H})$ on any Hilbert space $\mathcal{H}$ by $A\sim_{\phi}B$ iff $A=e^{i\theta}B$ for some $\theta\in[0,2\pi)$, and we have used standard qubit rotation operators $R_{j}(\theta)=\exp(-i\pi\theta\sigma_{j}/2)$. Note that although Eq. 6 is written for $\mathcal{H}_{A}$, the three-dimensional subspace of $\mathcal{H}$ involving all two qubit states except for $\ket{0,0}$, it would also hold for the full $\mathcal{H}_{qq}\supset\mathcal{H}_{A}$ if the $\ket{0,0}$ state is included in $\mathcal{B}_{7}$, yet it is not relevant to $U_{A}$, which leaves $\ket{0,0}$ unchanged under this extension.

We next seek to mathematically quantify the leakage out of the computational subspace after the adiabatic transition; this manifests itself as error in the adiabatic approximation, which we derive first in the general case and then specify it to the leakages most important in this adiabatic passage. Given some time-dependent Hamiltonian $H(t)$ acting on a Hilbert space $\mathcal{H}$, we let $\{\ket{\phi_{n}(t)}\colon n\in\mathbbm{Z}^{+}\}$ be the set of satisfiers of the instantaneous eigenstate condition $H(t)\ket{\phi(t)}=E(t)\ket{\phi(t)}$. Since these form a complete orthonormal basis of $\mathcal{H}$, we can write a general solution to the Schrodinger equation as

for constants $c_{n}(t)\in\mathbbm{C}$. Applying the Schrodinger equation and taking the inner product with $\bra{\phi_{k}(t)}$, rearranging gives¹⁴

where $\mathcal{E}_{k}$ represents the error in the adiabatic approximation i.e. leakage for instantaneous eigenstate $\ket{\phi_{k}}$, defined by

This follows intuition as if $\mathcal{E}_{k}(t)\approx 0$, a system that starts in some $\ket{\phi_{k}(0)}=\ket{\psi(0)}$ will remain in the same eigenstate. Thus, to minimize leakage out of the computational subspace, we wish to minimize the inner product in Eq. 9. In the case we are interested in, after differentiating the instantaneous eigenstate condition and manipulating, the error can be rewritten as

In summary, minimizing the leakage involves, as expected, minimizing the time derivative of the flux trajectory, yet it is also important to keep energies of states at risk of coupling far from each other, as well as it is to minimize the corresponding matrix element of $\partial_{\varphi}H$.

The issue with the idealized CZ gate discussed in Sec. II above lies in the adiabatic approximation: perfect adiabaticity is not possible since $\partial_{t}H$ cannot non-trivially be the zero operator. In this section we shall discuss how to quantify adiabaticity for a trajectory $\varphi(t)$ using arguments related to the form of the Hamiltonian in Eq. 3. We shall work under the following assumptions and assertions for $\varphi\colon T\to\mathbbm{R}^{+}$:

1. 1.

   $\varphi\evaluated{}_{\partial T}=0$. Inaction at the endpoints ensures that the qubit frequency starts and ends up in the optimal single-qubit operational point where coupling between transmons is negligible.

2. 2.

   $\int_{T}{\!\zeta(\varphi(t))\,dt}\equiv\pi\pmod{2\pi}$, where $\zeta(\varphi)=\delta_{11}(\varphi)-\delta_{10}(\varphi)-\delta_{01}(\varphi)$. As mentioned in the previous section, this condition generates the relative phase applied to $\ket{1,1}$ during tuning.

3. 3.

   The set of implementable flux pulses is exactly $C_{0}(T)$. That is, any $\varphi(t)$ is physically continuous, and any $C_{0}$ function may be flux-implemented.

4. 4.

   $\varphi(t)$ does not induce a $\ket{1,1}\longleftrightarrow\ket{2,0}$ transition with high likelihood, i.e. $\bra{2,0(t)}\partial_{\varphi}H\ket{1,1(t)}\approx 0$.

5. 5.

   $\varphi(t)$ does not induce a $\ket{0,1}\longleftrightarrow\ket{1,0}$ transition with high likelihood, i.e. $\bra{1,0(t)}\partial_{\varphi}H\ket{0,1(t)}\approx 0$.

Incidentally, the last requirement is satisfied by keeping $\varphi(t)<\mathcal{A}_{1}\>\forall t\in T$ as shown in Fig. 2(a), and we can at least partially satisfy condition 4 by setting $\varphi(t)\leq\mathcal{A}_{2}$, thus we restrict the codomain of $\varphi$ to $[0,\mathcal{A}_{2}]$. Keeping these conditions in mind, we look to the $\mathcal{H}_{Z}$ Hamiltonian in Eq. 3, and transform into the rotating frame set by $H(\varphi=0)$, which has ‘central’ matrix element

where $\alpha_{1}$ and $\alpha_{2}$ are the anharmonicities of the transmons, $\delta$ is the deviation of the qubit frequency of transmon one, and $\Delta$ is the detuning between the qubit frequencies of the two transmons at $\varphi=0$. The eigenfrequency difference $\zeta(\varphi)$ is most efficiently raised when $\varphi$ is kept near the $\mathcal{A}_{2}$ avoided crossing as in Fig. 2(c), and most time should be spent there in the ideal case. This point is unfortunately the most dangerous with respect to condition 4, and we seek to develop a rough ansatz for how to construct $\varphi(t)$ such that this is carefully handled. We see that the unwanted transitions occur in the latter terms of Eq. 11, which, when restricted to the interaction subspace $\mathcal{H}_{I}$ spanned by $\{\ket{1,1},\ket{2,0}\}$, produces the propagator $U_{\mathrm{Ev}}(t)\in\mathcal{U}(\mathcal{H}_{I})$, where

where $\mathcal{T}$ is the Dyson time ordering operator¹⁵ (and we use the notation $\mathcal{U}(V)$ to signify the group of unitary operators on a vector space $V$). Ideally the time-evolution restricted to $\mathcal{H}_{I}$ is just the two-dimensional identity, which should not create any swapping interaction; in this case, clearly $\delta(\varphi)$ must be minimized to remove $\ket{2,0}$ amplification, and while the off-diagonal terms cannot be easily controlled, repeated integration by parts shows that their contribution is lowered with sufficiently distant initial qubit frequencies.

We address the issue of minimizing $\delta(\varphi)$ and $\partial_{t}H$ simultaneously by promoting the set of continuously differentiable choices for $\varphi$, that is $S=C_{1}\,\cap\,[0,\mathcal{A}_{2}]^{T}$, to a semi-normed associative algebra $\mathcal{N}=(S,\|\!\cdot\!\|_{\mathcal{N}}\colon\mathcal{N}\to\mathbbm{R},+,\cdot)$ under usual function addition and multiplication. We make a simple but non-trivial choice for $\|\cdot\|_{\mathcal{N}}$, where, for constants $\upsilon,\kappa$,

(note that the semi-norm may take on negative values) which defines a ‘large’ element as one which is generates a large change in frequency and its time-derivative, i.e. one that is diabatic. Of course, this is just one arbitrary way of defining diabaticity, but it provides a good starting point for finding sufficiently adiabatic trajectories. Note that the desired trajectory $\varphi$ must additionally satisfy condition 2, and if one imagines the norm on $\mathcal{N}$ as a single-valued uncountably-infinite dimensional curve, then we are looking for the minimum of the intersection of this curve and the restriction on $\zeta(\varphi)$. This is easily approached using the Calculus of Variations, in which constrained variation via Lagrange multipliers is performed (for a rigorous treatment of this procedure see Ref. 16). In the general case, for a Lagrange multiplier $\lambda$, we minimize the functional $F[\varphi]$ while satisfying constraint functional $G[\varphi]=c_{0}\in\mathbbm{R}$ by imposing $\delta(F+\lambda G)=0$ ($\delta$ being the first variation). In our case, $F[\varphi]=\|\varphi\|_{\mathcal{N}}$ and $G[\varphi]=\int_{T}{\!\zeta(\varphi(t))\,dt}\equiv\pi\pmod{2\pi}$. If the functional sum is treated as an effective classical action, then the effective Lagrangian is

where we have used the definition for $\delta(t)$ and removed the arbitrary multiplicative factor $\omega_{1}-\alpha_{1}$ by absorbing into the Lagrange multiplier. Applying the variation gives

a highly non-linear second order equation after absorbing more constants into $\lambda$. With this, we impose that $\varphi\evaluated{}_{\partial T}=\dot{\varphi}\evaluated{}_{\partial T}=0$ along with the constraint equation for the multiplier. This, in fact, is problematic enough to discourage physical realization of the optimal solution; the main issue is the trigonometric singularity on $\partial T$ induced by the boundary condition in the first term. We shall instead use Eq. 15 as another method of gauging how diabatic a trajectory $\varphi$ is; this equation, however, gives us a distance at each $t$ from an optimal solution curve, something more valuable than the norm presented above. As for the singularity, we note that only relative diabaticity is relevant to this discussion, and we shall thus consider trajectories that start at some small $\varepsilon\gtrsim 0$.

To proceed, we utilize Eq. 15 through a functional operator $\mathcal{D}\colon S^{\prime}\to C_{1}$ equal to the LHS, where $S^{\prime}=C_{1}\,\cap\,[0,\mathcal{A}_{2}]^{T^{\prime}}$, $T^{\prime}=[\varepsilon,\tau]$, and $\tau$ is the CZ gate time. As potential flux trajectories, we shall consider a few possible curves detailed below for demonstration (all temporal units are in n$s$):

1. 1.

   Standard Gaussian. For this we use the standard equation

   $$\mathscr{G}(t,\sigma,\tau/2)=\mathcal{P}_{\sigma,\tau}C_{\sigma,\tau}e^{-\frac$$

   (16)

   where $C_{\sigma,\tau}$ normalizes the usual part of the distribution to one, $\mathcal{P}_{\sigma,\mu}$ satisfies the $\pi\pmod{2\pi}$ requirement on $\zeta$, and $\mathcal{Z}_{\sigma,\tau}$ ensures that $\mathscr{G}(0,\sigma)=\mathscr{G}(\tau,\sigma)=0$. This will act as a control relative to the other pulses, as we expect to find that by increasing $\sigma$ (and subsequently $\tau$, otherwise there is a diabatic jump on the boundary) we can arbitrarily reduce the diabaticity, and that being the only parameter makes this a somewhat trivial case. From now on, we shall suppress the constants’ dependence on $\tau$.

2. 2.

   Isolated Mollifier. We introduce a mollifier curve, namely

   $$\mathscr{M}(t,\sigma,\tau/2)=\left\{\begin{matrix}\mathcal{P}_{\sigma}C_{$$

   (17)

   an often-used object when a point of nondifferentiability must be smoothed out (here $\mathcal{W}=[\tau/2-\sigma,\tau/2+\sigma]$). Like the Gaussian waveform, this has only one width-parameter, and therefore acts again as a trivial control with respect to optimization. This curve has the advantage that all of its derivatives are exactly zero on $\partial T$.

3. 3.

   Mollified Gaussian Convolving the two previous trajectories, we have a less trivial case, in which it is of interest where each curve is placed. This takes the following form:

   $$\mathscr{F}(t,\sigma,\tau/2,\mu)=\int_{T}{\!\mathscr{G}(t,\sigma,\tau/2)$$

   (18)

   Adjustment of the parameter $\mu$ changes where the convolution occurs; we have also suppressed any normalization constants which maintain the $\zeta$ requirement.

4. 4.

   Prepulsed Gaussian. We ‘prepulse’ the above Gaussian waveform, giving

   $$\mathscr{H}(t,\sigma,\tau/2,\mu)=\frac{1}{3}\mathscr{G}(t,\sigma,\mu)+\mathscr$$

   (19)

   by using another Gaussian at some point during the trajectory; this could relieve some of the diabaticity generated by the original curve since it concentrates some of the ‘mass’ to a different area. Though the arbitrary factor of 1/3 was chosen as a relative height between the prepulse and the main Gaussian, the important effects should at least somewhat be observed in this trial.

5. 5.

   Mollifier-Prepulsed Gaussian Similar to the previous curve, we again prepulse a Gaussian trajectory, however this time we do so using an isolated mollifier, which becomes

   $$\mathscr{J}(t,\sigma,\tau/2,\mu)=\frac{1}{6}\mathscr{M}(t,\sigma,\mu)+\mathscr$$

   (20)

These trajectories are shown in Fig. 3 with fixed example parameters.

While challenging to classify the results of a simulation under a specific trajectory, we can simplify the procedure by assuming perfect linearity on the restriction to $\mathcal{H}_{Z}$; we thus have a representation of the operator by temporally evolving the elements of the tensor product basis set. We can construct such a matrix representation using

where $\mathcal{B}_{8}=\mathcal{B}_{7}\,\cup\,\{\ket{0,0}\}$. We have included the basis vector $\ket{0,0}$ to draw an accurate comparison between the ideal CPHASE and $\mathcal{M}(\varphi)$ (note that $U_{A}$ as used here has been extended to act on a Hilbert space, $\mathcal{H}_{8}$, larger than $\mathcal{H}_{Z}$ which is spanned by $\mathcal{B}_{7}$). Defining the accuracy of $\mathcal{M}$ is equivalent to choosing a matrix norm, thus for simplicity we make use of the $\ell_{2}$ spectral norm¹⁷ to compute the distance from $\mathcal{M}$ to a perfect $CZ$. We define $\mathcal{F}:\mathcal{H}_{8}\to\mathbbm{R}^{+}$ such that