Simulating an SNZ gate¶

This tutorial shows you step-by-step how to simulate an SNZ gate using the qruise-toolset.

The tutorial consists of the following:

1. What is an SNZ gate?
2. System Hamiltonian (Theory)
3. Loading parameters
4. Defining the Hamiltonian
5. Flux dependence of system energies
6. Defining the flux signal
7. Solving the Schrödinger equation
8. Plotting the qubit dynamics
9. Reducing the Hilbert space
10. Exploring flux and pulse duration dependence

1. What is an SNZ gate?¶

An SNZ (sudden net-zero) gate is a common two-qubit gate used in superconducting qubit systems with tunable frequencies, such as transmons. By applying a carefully designed flux pulse to one or both transmons, the system is temporarily brought into resonance, enabling controlled quantum transitions between specific states (e.g., $|11\rangle$ and $|02\rangle$). The net-zero flux pulse ensures that the qubit frequencies return to their original values at the end of the gate, thereby preserving coherence and preventing further unintended interactions.

SNZ retains the benefits of a conventional NZ (non-zero) gate, such as minimising low-frequency flux noise and leakage through destructive interference, and preventing long-term flux-control distortions for enhanced gate repeatability. In addition, the SNZ gate maximises intermediate leakage to the non-computational state for faster operation and simplifies tune-up thanks to the predictable relationship between conditional phase, leakage, and the corresponding pulse parameters.

To find out more about the SNZ gate, check out the original paper on arXiv.

2. System Hamiltonian (Theory)¶

The system we are interested in here consists of two coupled transmons. We can write a general Hamiltonian, $H$, given by:

$$ H =\sum_{k=1,2} H_{k}(t) + H_{12} \tag{1} $$

The first term includes the Hamiltonian for each of the individual qubits, $H_k$ ($k=1,2$):

1. Target qubit: $H_{1} = \omega_{1} a_1^{\dagger} a_1 +\frac{\alpha_{1}}{2} a_1^{\dagger} a_1^{\dagger} a_1 a_1$

2. Control qubit: $H_{2}(t) = \tilde{\omega}_{2}(\Phi_{2}(t)) a_2^{\dagger} a_2 +\frac{\alpha_{2}}{2} a_2^{\dagger} a_2^{\dagger} a_2 a_2$

Here, $a_k^\dagger$ and $a_k$ are the creation and annihilation operators for the $k^{\text{th}}$ qubit, respectively, and $\alpha_k$ is the anharmonicity. $\omega_1$ is the frequency of qubit 1, and $\tilde{\omega}_{2}(\Phi_{2e}(t))$ is the frequency of qubit 2 as a function of total flux across qubit 2, $\Phi_{2}(t) = \Phi_{2} + \Phi_{2e}(t)$, where $\Phi_{2}$ is the flux offset of qubit 2 and $\Phi_{2e}(t)$ is the external flux applied to it. $\tilde{\omega}_{2}(\Phi_{2}(t))$ is given by

$$ \tilde{\omega}_{2}(\Phi_{2}(t)) = \Big[ {\omega}_{2} -\alpha_2 \Big] f\left(\Phi_{2}(t)\right)^{1/2} + \alpha_2 , \tag{2} $$

where $\omega_2$ is the frequency of qubit 2 in the absence of an externally applied flux. $f(\Phi_{2}(t))$ is a time-dependent function quantifying the effect of the external flux on the qubit frequency given by

$$ f\left(\Phi_{2}(t)\right) = \sqrt{\cos ^2\left(\pi \frac{\Phi_{2}(t)}{\Phi_{0}}\right)+d_2^2 \sin ^2\left(\pi \frac{\Phi_{2}(t)}{\Phi_{0}}\right)}. \tag{3} $$

Here, $\Phi_{0}$ is the magnetic flux quantum, and $d_2$ is the junction asymmetry parameter of qubit 2.

The second term in the Hamiltonian, $H_{12}$, describes the coupling between the two qubits and is given by

$$ H_{12} = g_{12}(a_1^{\dagger} + a_1) (a_2^{\dagger}+ a_2), \tag{4} $$

where $g_{12}$ is the coupling strength between the two qubits.

---

Tip: Make sure that your environment is correctly set to handle float64 precision by setting JAX_ENABLE_X64=True or add the following codeblock to your scripts' preamble:

In [1]:

import jax

jax.config.update("jax_enable_x64", True)

import jax jax.config.update("jax_enable_x64", True)

---

3. Loading parameters¶

It's possible to define your system parameters directly in the code; however, it can often be more useful to load your parameters from somewhere external. For example, this could be the parameters obtained in the previous system calibration. This saves you having to update the code itself and helps to avoid errors that could be introduced by doing this.

We can load these parameters from a serialisable format like a JSON file. This is simple as it’s widely supported and easily parsed across many programming languages. Using serialisable files streamlines storing, sharing, and managing configurations, enabling quick updates without having to update the code.

• The sys_params.json file contains the intrinsic physical properties of the qubit system (e.g., frequency, anharmonicity, coupling strength, ...). These parameters are determined by hardware calibration (or learning process). It looks something like this:

In [2]:

{
    "q1": {
        "freq": 31730085801.256912,
        "anh": -2133496411.7573252,
        "phi0": 0.616981350058629,
        "fluxoffset": -1.732136553680111e-05,
        "d": -4.3105374201172966e-06,
    },
    "q2": {...},
    "q1-q2": {"coupling": 43503074.40583586},
}

{ "q1": { "freq": 31730085801.256912, "anh": -2133496411.7573252, "phi0": 0.616981350058629, "fluxoffset": -1.732136553680111e-05, "d": -4.3105374201172966e-06, }, "q2": {...}, "q1-q2": {"coupling": 43503074.40583586}, }

Out[2]:

{'q1': {'freq': 31730085801.256912,
  'anh': -2133496411.7573252,
  'phi0': 0.616981350058629,
  'fluxoffset': -1.732136553680111e-05,
  'd': -4.3105374201172966e-06},
 'q2': {Ellipsis},
 'q1-q2': {'coupling': 43503074.40583586}}

• The opt_params.json file contains the settings that optimise or simulate the system's behavior. These parameters do not define the system but are used during computational simulations (e.g., time intervals, volts-to-Hertz, ...). It looks like this:

In [3]:

{
    "res_flux": 0.11962114736749387,
    "v2hz": 1.089233172237356,
    "dt": 1e-11,
    "pulse_shape": "unipolar",
    "levels": 3,
}

{ "res_flux": 0.11962114736749387, "v2hz": 1.089233172237356, "dt": 1e-11, "pulse_shape": "unipolar", "levels": 3, }

Out[3]:

{'res_flux': 0.11962114736749387,
 'v2hz': 1.089233172237356,
 'dt': 1e-11,
 'pulse_shape': 'unipolar',
 'levels': 3}

In [4]:

import json

# read in JSON files and store the data in params dictionaries
with open("sys_params.json", "r") as file:
    sys_params = json.load(file)
with open("opt_params.json", "r") as file:
    opt_params = json.load(file)

# assign system parameters for each qubit
params_q1 = sys_params["q1"]
params_q2 = sys_params["q2"]

import json # read in JSON files and store the data in params dictionaries with open("sys_params.json", "r") as file: sys_params = json.load(file) with open("opt_params.json", "r") as file: opt_params = json.load(file) # assign system parameters for each qubit params_q1 = sys_params["q1"] params_q2 = sys_params["q2"]

4. Defining the Hamiltonian¶

Before we can code our Hamiltonian, we need to import and/or define the relevant operators. To simplify things, we'll use the following variables to represent key operator combinations:

• $\texttt{n\_k} = a^{\dagger}_{k} a_{k}$ (occupation number operator)
• $\texttt{anh\_k} = a^{\dagger}_{k} a^{\dagger}_{k} a_{k} a_{k}$ (anharmonicity operator)
• $\texttt{op\_kl} = (a^{\dagger}_{k} + a_{k} ) (a^{\dagger}_{l} + a_{l})$ (interaction operator)

where $k,l=1,2$.

In [5]:

import qutip as qt

num_subsystems = 2  # define number of subsystems, i.e., number of qubits

# define annihilation and identity operators for number of energy "levels" defined in JSON (for a transmon, "levels" = 3)
a = qt.destroy(opt_params["levels"])
Id = qt.identity(opt_params["levels"])

# define occupation number operator for each qubit
n_op = a.dag() @ a
n_1 = qt.tensor(n_op, Id)
n_2 = qt.tensor(Id, n_op)

# define anharmonicity operator (a^dagger*a^dagger*a*a) for each qubit
anh_op = a.dag() @ a.dag() @ a @ a
anh_1 = qt.tensor(anh_op, Id)
anh_2 = qt.tensor(Id, anh_op)

# define interaction operator, (a_1^dagger + a_1)(a_2^dagger + a_2)
op_a = a.dag() + a
op_12 = qt.tensor(op_a, op_a)

import qutip as qt num_subsystems = 2 # define number of subsystems, i.e., number of qubits # define annihilation and identity operators for number of energy "levels" defined in JSON (for a transmon, "levels" = 3) a = qt.destroy(opt_params["levels"]) Id = qt.identity(opt_params["levels"]) # define occupation number operator for each qubit n_op = a.dag() @ a n_1 = qt.tensor(n_op, Id) n_2 = qt.tensor(Id, n_op) # define anharmonicity operator (a^dagger*a^dagger*a*a) for each qubit anh_op = a.dag() @ a.dag() @ a @ a anh_1 = qt.tensor(anh_op, Id) anh_2 = qt.tensor(Id, anh_op) # define interaction operator, (a_1^dagger + a_1)(a_2^dagger + a_2) op_a = a.dag() + a op_12 = qt.tensor(op_a, op_a)

We can now define our stationary ($H_0$) and time-dependent ($H_t$) Hamiltonian terms using these operators and the system parameters we loaded earlier. Taking only the time-independent terms from Equation 1, we end up with

$$ H_{0} = \underbrace{\omega_{1} a_{1}^{\dagger} a_{1} + \frac{\alpha_{1}}{2} a_{1}^{\dagger} a_{1}^{\dagger} a_{1} a_{1}}_{\text{H}_{1}} + \underbrace{\frac{\alpha_{2}}{2} a_{2}^{\dagger} a_{2}^{\dagger} a_{2} a_{2}}_{\tilde{\text{H}}_{2}} + \underbrace{g_{12} (a^{\dagger}_{1} + a_{1})(a^{\dagger}_{2} + a_{2})}_{\text{H}_{12}}. \tag{5} $$

Our time-dependent term is then given by

$$ H_{t} =a_{2}^{\dagger}a_{2}. \tag{6} $$

---

Note: $\tilde{\omega}_2(\Phi_{2e}(t))$ is omitted here because it represents the time-dependent frequency of qubit 2, which will be applied as a modulation factor to $ H_t = a_{2}^{\dagger} a_{2} $ later. This allows us to dynamically control the frequency based on the external flux.

---

In [6]:

import jax.numpy as jnp

# define individual stationary Hamiltonians
H_1 = params_q1["freq"] * n_1 + 0.5 * params_q1["anh"] * anh_1  # qubit 1
H_2 = 0.5 * params_q2["anh"] * anh_2  # qubit 2
H_12 = sys_params["q1-q2"]["coupling"] * op_12  # H_12 (coupling between qubit 1 and 2)

# combine stationary and time-dependent Hamiltonians
H_0 = jnp.asarray((H_1 + H_2 + H_12).full())
H_t = jnp.asarray(n_2.full())

import jax.numpy as jnp # define individual stationary Hamiltonians H_1 = params_q1["freq"] * n_1 + 0.5 * params_q1["anh"] * anh_1 # qubit 1 H_2 = 0.5 * params_q2["anh"] * anh_2 # qubit 2 H_12 = sys_params["q1-q2"]["coupling"] * op_12 # H_12 (coupling between qubit 1 and 2) # combine stationary and time-dependent Hamiltonians H_0 = jnp.asarray((H_1 + H_2 + H_12).full()) H_t = jnp.asarray(n_2.full())

Now we need to define and label our basis states. Since the 2nd excited level of the transmon is used, the basis state are a combination of two three-level systems $|00\rangle, |01\rangle, |02\rangle, |10\rangle, |11\rangle, |12\rangle, |20\rangle, |21\rangle, |22\rangle$.

In [7]:

from qruise.toolset.utils import computational_basis

# define basis states and label ("levels"=3, num_subsystems=2)
basis, labels = computational_basis([opt_params["levels"]] * num_subsystems)

from qruise.toolset.utils import computational_basis # define basis states and label ("levels"=3, num_subsystems=2) basis, labels = computational_basis([opt_params["levels"]] * num_subsystems)

5. Flux dependence of system energies¶

It's a good idea to verify that your previously defined resonant flux (loaded from opt_params.json and stored in opt_params["res_flux']) corresponds to a transition from the $|11\rangle$ to $|02\rangle$ state. We can do this by plotting the flux dependence of the computational basis energies and extracting the flux at which the $|11\rangle$ to $|02\rangle$ states intersect. We can also use this to determine a suitable flux range to ensure that no other states are involved in the gate.

We can calculate the computational basis energies by removing the coupling term ($H_{12}$) from the stationary Hamiltonian ($H_0$), and by introducing the $\Phi_{2}$ flux dependence as follows:

$$ H = H_{1} + \tilde{H}_{2} + \tilde{\omega}_{2}(\Phi_{2})H_{t}. \tag{7}$$

(Reminder: $\tilde{\omega}_{2}(\Phi_{2})$ is the flux-dependent frequency of qubit 2 as defined in Equation 2.)

In [8]:

# define flux-dependent frequency of qubit 2 as a function of flux
def tunable_frequency(freq, anh, fluxoffset, phi0, d, dcflux):
    flux = fluxoffset + dcflux
    f = jnp.sqrt(
        jnp.cos(jnp.pi * flux / phi0) ** 2 + d**2 * jnp.cos(jnp.pi * flux / phi0) ** 2
    )
    return (freq - anh) * f**0.5 + anh


# loop over flux values: control qubit 2
flux_values = jnp.linspace(-0.25, 0.25, 251)  # define range of flux values to explore
eigenval_flux = jax.vmap(
    lambda x: jnp.linalg.eigvals(jnp.asarray((H_1 + H_2).full()) + x * H_t)
)(
    tunable_frequency(
        **sys_params["q2"], dcflux=flux_values * opt_params["v2hz"]
    )  # calculate energies of qubit states in computational basis
)

# define flux-dependent frequency of qubit 2 as a function of flux def tunable_frequency(freq, anh, fluxoffset, phi0, d, dcflux): flux = fluxoffset + dcflux f = jnp.sqrt( jnp.cos(jnp.pi * flux / phi0) ** 2 + d**2 * jnp.cos(jnp.pi * flux / phi0) ** 2 ) return (freq - anh) * f**0.5 + anh # loop over flux values: control qubit 2 flux_values = jnp.linspace(-0.25, 0.25, 251) # define range of flux values to explore eigenval_flux = jax.vmap( lambda x: jnp.linalg.eigvals(jnp.asarray((H_1 + H_2).full()) + x * H_t) )( tunable_frequency( **sys_params["q2"], dcflux=flux_values * opt_params["v2hz"] ) # calculate energies of qubit states in computational basis )

We can then identify our flux range by plotting the state energies against the flux and selecting a range where only the $|11\rangle$ and $|02\rangle$ states intersect.

In [9]:

from qruise.toolset.plots import PlotUtil, qruise_colour_palette
from bokeh.palettes import Colorblind8

# plot frequency of qubit 2 as a function of flux
canvas = PlotUtil(
    x_axis_label="Control qubit flux [V]",
    y_axis_label="Energy [GHz/2π]",
    notebook=True,
    colour_palette=Colorblind8,
)
canvas.plot(
    flux_values,
    jnp.abs(eigenval_flux / 1e9 / (2.0 * jnp.pi)),
    labels=[str(label) for label in labels],
)
canvas.show_canvas()

from qruise.toolset.plots import PlotUtil, qruise_colour_palette from bokeh.palettes import Colorblind8 # plot frequency of qubit 2 as a function of flux canvas = PlotUtil( x_axis_label="Control qubit flux [V]", y_axis_label="Energy [GHz/2π]", notebook=True, colour_palette=Colorblind8, ) canvas.plot( flux_values, jnp.abs(eigenval_flux / 1e9 / (2.0 * jnp.pi)), labels=[str(label) for label in labels], ) canvas.show_canvas()

6. Defining the flux signal¶

For the flux signal ($\Phi_{2e}$), we will use a Heaviside step function, $\Theta$, to generate a square pulse:

$$ \Phi_{2e}(t) = \begin{cases} \Phi_{dc} \Theta(t) \Theta(t_{final} - t) & \text{for the unipolar case} \\ \\ \Phi_{dc} \Theta(t) \Theta(\frac{t_{final}}{2} - t) - \Phi_{dc} \Theta(t-\frac{t_{final}}{2}) \Theta(t_{final} - t) & \text{for the bipolar case} \end{cases} \tag{8} $$

In the unipolar case, the flux pulse has a constant magnitude of $\Phi_{dc}$ (the applied DC flux) between time $t=0$ and $t=t_{final}$, and is zero outside this range.

In the bipolar case, the flux pulse has a constant magnitude of $\Phi_{dc}$ between time $t=0$ and $t=\frac{t_{final}}{2}$, a constant magnitude of $-\Phi_{dc}$ between time $t=\frac{t_{final}}{2}$ and $t=t_{final}$, and is zero outside this range.

To define the pulse polarity, set opt_params["pulse shape"] = unipolar for unipolar or opt_params["pulse shape"] = bipolar for bipolar.

In [10]:

# set pulse shape to bipolar (net-zero)
opt_params["pulse_shape"] = "bipolar"

# set pulse shape to bipolar (net-zero) opt_params["pulse_shape"] = "bipolar"

To begin with, we need an estimate of the required pulse length to observe a half oscillation, i.e., a population transfer from the $|11\rangle$ to $|02\rangle$ state. We can do this by calculating the effective coupling strength between the qubits, $J_{12}$, which is given by

$$ J_{12} = \frac{2 \sqrt{2}g_{12}}{ 2 \pi} , \tag{9} $$

where $g_{12}$ is the bare qubit-qubit coupling strength (see Equation 4). The required pulse length, $t_{final}$, to observe a half oscillation in a two-level system is then calculated as

$$ t_{final} = \frac{1}{2 J_{12}} . \tag{10} $$

As transmons are three-level systems, this value for $t_{final}$ will not be perfect, but it's a reasonable starting parameter.

In [11]:

# estimate required pulse length for a two-level system
J12 = 2 * jnp.sqrt(2) * sys_params["q1-q2"]["coupling"] / (2 * jnp.pi)
t_final = float(0.5 / J12)

dt = opt_params["dt"]  # set time interval for calculations

# estimate required pulse length for a two-level system J12 = 2 * jnp.sqrt(2) * sys_params["q1-q2"]["coupling"] / (2 * jnp.pi) t_final = float(0.5 / J12) dt = opt_params["dt"] # set time interval for calculations

Now we need to set the pulse parameters. For this, we use opt_params and params_q2 that we loaded in earlier.

In [12]:

# define pulse parameters using the qubit 2 values loaded in earlier
# False for a fixed value, True to allow for optimisation of the parameter later
parameters = {
    "dcflux": (opt_params["res_flux"], False),  # magnitude of DC flux (V)
    "tfinal": (t_final, False),  # pulse length (s)
    "freq": (params_q2["freq"], False),  # qubit frequency (Hz)
    "anh": (params_q2["anh"], False),  # anharmonicity (Hz)
    "fluxoffset": (params_q2["fluxoffset"], False),  # Phi_2
    "phi0": (
        params_q2["phi0"],
        False,
    ),  # magnetic flux quantum (includes Volts conversion factor)
    "d": (params_q2["d"], False),  # junction asymmetry parameter
    "v2hz": (
        opt_params["v2hz"],
        False,
    ),  # pulse amplitude conversion factor (Volts to Hz)
}

# define pulse parameters using the qubit 2 values loaded in earlier # False for a fixed value, True to allow for optimisation of the parameter later parameters = { "dcflux": (opt_params["res_flux"], False), # magnitude of DC flux (V) "tfinal": (t_final, False), # pulse length (s) "freq": (params_q2["freq"], False), # qubit frequency (Hz) "anh": (params_q2["anh"], False), # anharmonicity (Hz) "fluxoffset": (params_q2["fluxoffset"], False), # Phi_2 "phi0": ( params_q2["phi0"], False, ), # magnetic flux quantum (includes Volts conversion factor) "d": (params_q2["d"], False), # junction asymmetry parameter "v2hz": ( opt_params["v2hz"], False, ), # pulse amplitude conversion factor (Volts to Hz) }

We can then define the flux pulse, $\Phi_{2e}(t)$, using the Heaviside function as described in Equation 8 above.

In [13]:

from toolz import curry


# define flux Phi(t) using Heaviside step function for unipolar and bipolar cases (Equation 8)
def phi(t, dcflux, tfinal, fluxoffset, v2hz):
    if opt_params["pulse_shape"] == "unipolar":
        return (
            fluxoffset
            + dcflux * jnp.heaviside(t, 1) * jnp.heaviside(tfinal - t, 1) * v2hz
        )  # equation for unipolar pulse
    elif opt_params["pulse_shape"] == "bipolar":
        return (
            fluxoffset
            + dcflux * jnp.heaviside(t, 1) * jnp.heaviside(tfinal / 2 - t, 1) * v2hz
            - dcflux
            * jnp.heaviside(t - tfinal / 2, 1)
            * jnp.heaviside(tfinal - t, 1)
            * v2hz
        )  # equation for bipolar pulse
    else:
        raise ValueError("Invalid pulse shape, must be 'unipolar' or 'bipolar'")

from toolz import curry # define flux Phi(t) using Heaviside step function for unipolar and bipolar cases (Equation 8) def phi(t, dcflux, tfinal, fluxoffset, v2hz): if opt_params["pulse_shape"] == "unipolar": return ( fluxoffset + dcflux * jnp.heaviside(t, 1) * jnp.heaviside(tfinal - t, 1) * v2hz ) # equation for unipolar pulse elif opt_params["pulse_shape"] == "bipolar": return ( fluxoffset + dcflux * jnp.heaviside(t, 1) * jnp.heaviside(tfinal / 2 - t, 1) * v2hz - dcflux * jnp.heaviside(t - tfinal / 2, 1) * jnp.heaviside(tfinal - t, 1) * v2hz ) # equation for bipolar pulse else: raise ValueError("Invalid pulse shape, must be 'unipolar' or 'bipolar'")

It’s a good idea to check your pulse looks as expected before proceeding, so let’s quickly plot it.

In [14]:

# define time axis values
ts = jnp.arange(-0.05 * t_final, 1.05 * t_final, dt)

# plot flux pulse
canvas.init_canvas(
    x_axis_label="t [ns]",
    y_axis_label="Flux amplitude [V]",
    colour_palette=qruise_colour_palette,
)
canvas.plot(
    ts / 1e-9,
    phi(
        ts,
        parameters["dcflux"][0],
        parameters["tfinal"][0],
        parameters["fluxoffset"][0],
        parameters["v2hz"][0],
    ),
    labels=["Bipolar pulse"],
)
canvas.show_canvas()

# define time axis values ts = jnp.arange(-0.05 * t_final, 1.05 * t_final, dt) # plot flux pulse canvas.init_canvas( x_axis_label="t [ns]", y_axis_label="Flux amplitude [V]", colour_palette=qruise_colour_palette, ) canvas.plot( ts / 1e-9, phi( ts, parameters["dcflux"][0], parameters["tfinal"][0], parameters["fluxoffset"][0], parameters["v2hz"][0], ), labels=["Bipolar pulse"], ) canvas.show_canvas()

Great! Our bipolar pulse behaves as expected. We can now use it, together with Equations 2 and 3, to define the flux-dependent frequency of qubit 2.

In [15]:

# define flux-dependent frequency of qubit 2
@curry
def freq_flux(t, dcflux, tfinal, freq, anh, fluxoffset, phi0, d, v2hz):
    # define f(t) (Equation 3)
    def f(t):
        return jnp.sqrt(
            jnp.cos(jnp.pi * phi(t, dcflux, tfinal, fluxoffset, v2hz) / phi0) ** 2
            + d**2
            * jnp.sin(jnp.pi * phi(t, dcflux, tfinal, fluxoffset, v2hz) / phi0) ** 2
        )

    # return frequency of qubit 2 as a function of flux (Equation 2)
    return (freq - anh) * f(t) ** 0.5 + anh

# define flux-dependent frequency of qubit 2 @curry def freq_flux(t, dcflux, tfinal, freq, anh, fluxoffset, phi0, d, v2hz): # define f(t) (Equation 3) def f(t): return jnp.sqrt( jnp.cos(jnp.pi * phi(t, dcflux, tfinal, fluxoffset, v2hz) / phi0) ** 2 + d**2 * jnp.sin(jnp.pi * phi(t, dcflux, tfinal, fluxoffset, v2hz) / phi0) ** 2 ) # return frequency of qubit 2 as a function of flux (Equation 2) return (freq - anh) * f(t) ** 0.5 + anh

We now need to create a signal chain, add the flux component, and apply this to create the drive.

In [16]:

from qruise.toolset import Signal

sig = Signal()  # create instance of signal chain
sig.add("flux", freq_flux, parameters)  # add flux component to signal chain

# create drive after applying signal chain
drive, params = sig.function()
drive = jax.jit(drive)

from qruise.toolset import Signal sig = Signal() # create instance of signal chain sig.add("flux", freq_flux, parameters) # add flux component to signal chain # create drive after applying signal chain drive, params = sig.function() drive = jax.jit(drive)

We should also plot the drive (the flux-dependent frequency of qubit 2) to ensure it behaves as expected.

In [17]:

ts = jnp.arange(0.0, t_final, dt)

# plot drive
canvas.init_canvas(
    x_axis_label="t [ns]",
    y_axis_label="Flux-dependent frequency of qubit 2 [GHz/2π]",
    colour_palette=qruise_colour_palette,
)
canvas.plot(ts / 1e-9, drive(ts, params) / 1e9 / (2.0 * jnp.pi), labels=["Frequency"])
canvas.show_canvas()

ts = jnp.arange(0.0, t_final, dt) # plot drive canvas.init_canvas( x_axis_label="t [ns]", y_axis_label="Flux-dependent frequency of qubit 2 [GHz/2π]", colour_palette=qruise_colour_palette, ) canvas.plot(ts / 1e-9, drive(ts, params) / 1e9 / (2.0 * jnp.pi), labels=["Frequency"]) canvas.show_canvas()

Nice! The drive behaves as expected and we can incorporate it into our Hamiltonian (see Equation 7).

In [18]:

from qruise.toolset import Hamiltonian

H = Hamiltonian(H_0)  # create stationary Hamiltonian, H0
H.add_term(H_t, drive)  # add time-dependent part, Ht, with drive

from qruise.toolset import Hamiltonian H = Hamiltonian(H_0) # create stationary Hamiltonian, H0 H.add_term(H_t, drive) # add time-dependent part, Ht, with drive

7. Solving the Schrödinger equation¶

We now want to simulate the population evolution for a given DC flux ($\Phi_{dc}$) and pulse duration ($t_{final}$), starting from the initial state $|11\rangle$.

To do this, we need to solve the Schrödinger equation, which requires us to specify the equations and parameters that govern the system. In the qruise-toolset, we combine these using a Problem. To instantiate the Problem, we need:

• the Hamiltonian (H)
• the initial qubit state (initial_state, here $|11\rangle$)
• the pulse parameters (params)
• the time interval of the simulation (0.0 to tfinal).

In [19]:

from qruise.toolset import Problem

# define initial qubit state (|11>)
initial_state = basis[labels.index((1, 1))]

# define Problem
prob = Problem(H, initial_state, params, (0.0, t_final))

from qruise.toolset import Problem # define initial qubit state (|11>) initial_state = basis[labels.index((1, 1))] # define Problem prob = Problem(H, initial_state, params, (0.0, t_final))

The user can then define the type of solver they want to use and the equation they want to solve (for example, the Schrödinger equation, the master equation, or the Lindblad, equation).

In this case, we’ll choose to use the piece-wise constant solver (PWCSolver) to solve the Schrödinger equation (sepwc). We can then calculate the wavefunction of the system at each timestamp.

In [20]:

from qruise.toolset.solvers.solvers import PWCSolver

# use PWCSolver to solve Schrödinger equation
solver = PWCSolver(dt=dt, store=True)
solver.set_system(prob.sepwc())
res = solver.evolve(*prob.problem())

from qruise.toolset.solvers.solvers import PWCSolver # use PWCSolver to solve Schrödinger equation solver = PWCSolver(dt=dt, store=True) solver.set_system(prob.sepwc()) res = solver.evolve(*prob.problem())

---

Note: Setting store=True in the PWCSolver ensures that the results are stored at each timestamp during the simulation.

---

8. Plotting the qubit dynamics¶

We can view the results of our simulation by calculating and plotting the expectation value of the transmon state populations against time.

In [21]:

from qruise.toolset.utils import get_population

# calculate expectation value of each qubit state
population = get_population(res)

# plot populations against time
canvas.init_canvas(x_axis_label="t [ns]", y_axis_label="Population")
canvas.plot(ts / 1e-9, population, labels=[str(label) for label in labels])
canvas.show_canvas()

from qruise.toolset.utils import get_population # calculate expectation value of each qubit state population = get_population(res) # plot populations against time canvas.init_canvas(x_axis_label="t [ns]", y_axis_label="Population") canvas.plot(ts / 1e-9, population, labels=[str(label) for label in labels]) canvas.show_canvas()

Congratulations! You’ve just simulated an SNZ gate!

You can see, however, that the $|02\rangle$ state population only reacehs a maximum of around 97%. We can improve this by optimising the flux and pulse duration.

9. Reduced Hilbert space¶

You can see in the figure above that only the $|11\rangle$ and $|02\rangle$ states are significantly involved in the gate dynamics, with minimal contributions from other states. We can therefore greatly simplify our simulation by reducing the Hilbert basis to only the $|11\rangle$ and $|02\rangle$ states. This reduces the computational task from a 9x9 matrix to a 2x2 one, making the calculation much faster.

In [22]:

from qruise.toolset.utils import project_op

# define states for reduced Hilbert space
states_labels = [(1, 1), (0, 2)]
projection_states = [basis[labels.index(label)] for label in states_labels]

# project Hamiltonian onto |11>,|02> subspace
H_0_red = project_op(H_0, projection_states)  # reduced stationary Hamiltonian
H_t_red = project_op(H_t, projection_states)  # reduced time-dependent Hamiltonian term
H_red = Hamiltonian(H_0_red, [(H_t_red, drive)])  # final Hamiltonian with drive

from qruise.toolset.utils import project_op # define states for reduced Hilbert space states_labels = [(1, 1), (0, 2)] projection_states = [basis[labels.index(label)] for label in states_labels] # project Hamiltonian onto |11>,|02> subspace H_0_red = project_op(H_0, projection_states) # reduced stationary Hamiltonian H_t_red = project_op(H_t, projection_states) # reduced time-dependent Hamiltonian term H_red = Hamiltonian(H_0_red, [(H_t_red, drive)]) # final Hamiltonian with drive

Now we can solve the Schrödinger equation in the reduced Hilbert space and plot the dynamics.

In [23]:

# define initial qubit state |11>
initial_state_red = jnp.array([1.0, 0.0], dtype=jnp.complex128)

# define Problem
prob_red = Problem(H_red, initial_state_red, params, (0.0, t_final))

# use PWCSolver to solve Schrödinger equation
solver = PWCSolver(dt=dt, store=True)
solver.set_system(prob_red.sepwc())
res_red = solver.evolve(*prob_red.problem())

# calculate the expectation value of each qubit state
pop_red = get_population(res_red)

# plot dynamics
canvas.init_canvas(x_axis_label="t [ns]", y_axis_label="Population")
canvas.plot(ts / 1e-9, pop_red, labels=[str(x) for x in states_labels])
canvas.show_canvas()

# define initial qubit state |11> initial_state_red = jnp.array([1.0, 0.0], dtype=jnp.complex128) # define Problem prob_red = Problem(H_red, initial_state_red, params, (0.0, t_final)) # use PWCSolver to solve Schrödinger equation solver = PWCSolver(dt=dt, store=True) solver.set_system(prob_red.sepwc()) res_red = solver.evolve(*prob_red.problem()) # calculate the expectation value of each qubit state pop_red = get_population(res_red) # plot dynamics canvas.init_canvas(x_axis_label="t [ns]", y_axis_label="Population") canvas.plot(ts / 1e-9, pop_red, labels=[str(x) for x in states_labels]) canvas.show_canvas()

You can see from the figure that the dynamics are largely unchanged by moving into the reduced Hilbert space. This confirms the validity of this approach, which significantly enhances the computational efficiency.

10. Exploring flux and pulse duration dependence¶

We can now calculate the population of the $|02\rangle$ state as a function of applied DC flux, $\Phi_{dc}$, and pulse duration (t_final). This is achieved using an EnsembleProblem, which enables efficient computation by solving the Schrödinger equation across multiple values of flux and time simultaneously.

In [24]:

from copy import copy
from qruise.toolset import EnsembleProblem

# set simulation parameters
t_map = 120e-9  # total pulse duration
num_fluxes = 301  # number of flux values to be simulated
flux_amps = jnp.linspace(
    0.95 * opt_params["res_flux"], 1.05 * opt_params["res_flux"], num_fluxes
)  # create array of flux values

# set up ensemble problem
ensemb_params = copy(params)  # create a copy of params to avoid modifying the original
ensemb_params["sig1/flux/dcflux"] = flux_amps  # add flux array to params
ensemb_params["sig1/flux/tfinal"] = t_map  # set pulse duration for each simulation
ensemb_prob = EnsembleProblem(
    H_red, (initial_state_red, None), ensemb_params, (0.0, t_map)
)  # set up ensemble problem

result_flux_array = solver.ensemble_evolve(
    *ensemb_prob.problem()
)  # solve ensemble problem
pop = get_population(
    result_flux_array
)  # calculate expectation value for |02> state for each flux value and time

from copy import copy from qruise.toolset import EnsembleProblem # set simulation parameters t_map = 120e-9 # total pulse duration num_fluxes = 301 # number of flux values to be simulated flux_amps = jnp.linspace( 0.95 * opt_params["res_flux"], 1.05 * opt_params["res_flux"], num_fluxes ) # create array of flux values # set up ensemble problem ensemb_params = copy(params) # create a copy of params to avoid modifying the original ensemb_params["sig1/flux/dcflux"] = flux_amps # add flux array to params ensemb_params["sig1/flux/tfinal"] = t_map # set pulse duration for each simulation ensemb_prob = EnsembleProblem( H_red, (initial_state_red, None), ensemb_params, (0.0, t_map) ) # set up ensemble problem result_flux_array = solver.ensemble_evolve( *ensemb_prob.problem() ) # solve ensemble problem pop = get_population( result_flux_array ) # calculate expectation value for |02> state for each flux value and time

We can now plot our simulation as a 2D colour map.

In [25]:

from bokeh.plotting import figure, show

p = figure(
    x_axis_label="Pulse duration [ns]",
    y_axis_label="Pulse amplitude [V]",
    tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")],
)
p.x_range.range_padding = p.y_range.range_padding = 0

# must give a vector of image data for image parameter
p.image(
    image=[pop[:, :, 1]],
    x=t_map,
    y=float(flux_amps[0]),
    dw=t_map / 1e-9,
    dh=float(flux_amps[-1] - flux_amps[0]),
    palette="RdYlBu11",
    level="image",
)
p.grid.grid_line_width = 0.5

show(p)

from bokeh.plotting import figure, show p = figure( x_axis_label="Pulse duration [ns]", y_axis_label="Pulse amplitude [V]", tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")], ) p.x_range.range_padding = p.y_range.range_padding = 0 # must give a vector of image data for image parameter p.image( image=[pop[:, :, 1]], x=t_map, y=float(flux_amps[0]), dw=t_map / 1e-9, dh=float(flux_amps[-1] - flux_amps[0]), palette="RdYlBu11", level="image", ) p.grid.grid_line_width = 0.5 show(p)